MeshQuad._splitquads sometimes misorients second triangle

[gdmcbain]

Consider a mesh with a single quadrilateral:

0---1
|   |
2---3

i.e.

m = MeshQuad(np.array([[0, 1, 0, 1], [1, 1, 0, 0]]), np.array([[2, 3, 1, 0]]).T)

This is just the reference Q1 with differently numbered points (but ordered correctly in m.t).

Then

m._splitquads().t

returns

array([[0, 0],
       [2, 1],
       [3, 3]])

which has the second triangle clockwise.

This doesn't happen for the default quad:

3---2
|   |
0---1

MeshQuad()._splitquads().t returns

 array([[0, 1],
       [1, 2],
       [3, 3]])

both anticlockwise.

[gdmcbain]

It looks like this happens in the constructor for the new MeshTri, in https://github.com/kinnala/scikit-fem/blob/740021f05497a05be6a677866328d6679c66e635/skfem/mesh/mesh2d/mesh_tri.py#L130

where in all calls from _splitquads, sort_t==True.

[gdmcbain]

The sorting of the indices of the vertices of each triangle looks suspicious:

https://github.com/kinnala/scikit-fem/blob/740021f05497a05be6a677866328d6679c66e635/skfem/mesh/mesh2d/mesh_tri.py#L228

[gdmcbain]

I'm not sure what the intention of the sorting is here, but perhaps what's in order is to cycle each column so that it begins with its least entry? i.e. entry self.t.argmin(axis=0).

[gdmcbain]

numpy.roll looks like it might do that.

[gdmcbain]

As in:

    self.t = np.roll(self.t, self.t.argmin(axis=0), 0)

[kinnala]

Fixing this by removing the sort would require quite a lot of work, see https://github.com/kinnala/scikit-fem/pull/73#issuecomment-430503834

Do you find it problematic that the input orientation is not preserved?

[gdmcbain]

Ah, I'm not sure if it's problematic but it is a surprise. I had been working from the beginning on the false assumption that all two-dimensional finite element meshes had the nodal-element connectivity described so that the mapping from each element to the reference element had positive Jacobian. That implies kind of the opposite rule: that each internal edge has opposite senses in the two elements that it belongs to. I'm not sure where I picked that up; I suppose I associated it with the idea that in the application of Green's (and thus also Gauss's or Stokes's) theorem, the contributions from all internal edges would exactly cancel, leaving only the boundary edges to balance the summed area integrals over the faces.

I fear that in following this false fundamental assumption I may have committed many other errors.

[gdmcbain]

Hmm, no, that's not it: the two triangles produced by MeshQuad()._splitquads() do have the same orientation... I'm not understanding...

[gdmcbain]

Do you find it problematic that the input orientation is not preserved?

So where it could be problematic is in interoperation with other finite element tools (meshers, postprocessors, ). How I hit this was in constructing a MeshQuad by hand (all cells described with nodes anticlockwise), then _splitquads and .save to MEDIT .mesh and loadmesh into FreeFem++ to compare with a solution there. This failed as FreeFem++ aborted on finding triangles with negative Jacobian.

[kinnala]

Many FE codes indeed use such rules, i.e. det(J) should be positive. Here instead we use abs(det(J)) if the sign does not matter. Moreover if it does, such as when transforming Raviart-Thomas basis between reference and global elements using Piola transform, we can leave out abs and the orientation of the global basis is done automatically. I'll try to find a reference.

[kinnala]

https://arxiv.org/abs/1409.4618 Section 2.3.

[kinnala]

I believe in the case of ex02 the failure is due to this part unless you sort:

https://github.com/kinnala/scikit-fem/blob/bd544491e8670ed28ddb8f104c2ee25c2488c895/skfem/element.py#L259

Tangents are rotated 90 degrees counterclockwise. Unless neighboring elements agree on the direction of tangent we get non-continuous normal gradients in Morley element.

[gdmcbain]

Thanks for the reference (Anjam & Valdman 2015) and clear explanations; I'll close this, concluding that the bug was in my faulty assumption.