Ferrite says that detJ is negative, but Gmsh says that all element have a positive detJ

[Drachta]

Hello,

I'm new in FEM and in Ferrite. So far I have been able to use this library and Gmsh with different geometry and it worked fine, but this one gives me trouble:

Here is the error message:

Info    : Reading 'folder\file.msh'...
Info    : 18 entities
Info    : 507 nodes
Info    : 1021 elements
Info    : Done reading 'folder\file.msh'
ERROR: LoadError: ArgumentError: det(J) is not positive: det(J) = -1.1052236851259822
Stacktrace:
 [1] throw_detJ_not_pos(detJ::Float64)
   @ Ferrite C:\some_folder\.julia\packages\Ferrite\qYVTc\src\FEValues\common_values.jl:17
 [2] reinit!(cv::CellScalarValues{2, Float64, RefTetrahedron}, x::Vector{Vec{2, Float64}})
   @ Ferrite C:\some_folder\.julia\packages\Ferrite\qYVTc\src\FEValues\cell_values.jl:163
 [3] reinit!(cv::CellScalarValues{2, Float64, RefTetrahedron}, ci::CellIterator{2, Triangle, Float64, DofHandler{2, Triangle, Float64}})
   @ Ferrite C:\some_folder\.julia\packages\Ferrite\qYVTc\src\iterators.jl:130
 [4] top-level scope
   @ d:\some\folder\temp.jl:17
in expression starting at d:\some\folder\temp.jl:16

Here is a minimal geo file that reproduce the problem:

// Gmsh project created on Wed Mar 02 09:48:56 2022
//+
Point(1) = {0, 0, 0, 1.0};
//+
Point(2) = {0, 10, 0, 1.0};
//+
Point(3) = {10, 10, 0, 1.0};
//+
Point(4) = {10, 20, 0, 1.0};
//+
Point(5) = {20, 20, 0, 1.0};
//+
Point(6) = {20, 5, 0, 1.0};
//+
Point(7) = {30, 5, 0, 1.0};
//+
Point(8) = {30, 0, 0, 1.0};
//+
Line(1) = {1, 2};
//+
Line(2) = {2, 3};
//+
Line(3) = {3, 4};
//+
Line(4) = {4, 5};
//+
Line(5) = {5, 6};
//+
Line(6) = {7, 7};
//+
Line(7) = {6, 7};
//+
Line(8) = {7, 8};
//+
Line(9) = {8, 1};
//+
Curve Loop(1) = {2, 3, 4, 5, 7, 8, 9, 1};
//+
Plane Surface(1) = {1};

And here is a minimal Julia code to reproduce the error:

using Ferrite
using FerriteGmsh

folder = "some folder here"
grid = saved_file_to_grid(folder*"file.msh")

dim = 2
ip = Lagrange{dim, RefTetrahedron, 1}() # or RefCube depending on mesh type
qr = QuadratureRule{dim, RefTetrahedron}(2) # or RefCube depending on mesh type
cellvalues = CellScalarValues(qr, ip);

dh = DofHandler(grid)
push!(dh, :u, 1)
close!(dh);

@inbounds for cell in CellIterator(dh)
    reinit!(cellvalues, cell)
end

I have tried with different mesh type and size. I also tried feeding the geo file to Ferrite instead of a msh file, but the problem persist. Quad gives a higher detJ, and reducing the mesh size also gives higher detJ, but detJ stays negative. I can't really decrease the mesh size further.

Is there anything else I can try?

[fredrikekre]

Looks like the noder order is messed up:

julia> grid.cells[1]
Triangle((301, 405, 180))

julia> grid.nodes[[301, 405, 180]]
3-element Vector{Node{2, Float64}}:
 Node{2, Float64}([11.817618796000435, 18.15055360025356])
 Node{2, Float64}([11.491358029469588, 19.16755446688042])
 Node{2, Float64}([12.413274345295044, 19.1037038750435])

which gives roughly

2 ----- 3
 \     /
  \   /
   \1/

(clockwise), but Ferrite wants anti-clockwise:

3 ----- 2
 \     /
  \   /
   \1/

I guess something goes wrong in the mesh transfer @koehlerson ?

[koehlerson]

I'm not sure what we can do about it, since the gmsh docs say that the elements are provided anti-clockwise as we expect them: https://gmsh.info/doc/texinfo/gmsh.html#Low-order-elements

That's the reason why all other tested triangle meshes worked ootb

[koehlerson]

https://gitlab.onelab.info/gmsh/gmsh/-/issues/1169

Could you try to reverse the mesh:

ReverseMesh Surface{1};

This works for me

// Gmsh project created on Wed Mar 02 09:48:56 2022
//+
Point(1) = {0, 0, 0, 1.0};
//+
Point(2) = {0, 10, 0, 1.0};
//+
Point(3) = {10, 10, 0, 1.0};
//+
Point(4) = {10, 20, 0, 1.0};
//+
Point(5) = {20, 20, 0, 1.0};
//+
Point(6) = {20, 5, 0, 1.0};
//+
Point(7) = {30, 5, 0, 1.0};
//+
Point(8) = {30, 0, 0, 1.0};
//+
Line(1) = {1, 2};
//+
Line(2) = {2, 3};
//+
Line(3) = {3, 4};
//+
Line(4) = {4, 5};
//+
Line(5) = {5, 6};
//+
Line(6) = {7, 7};
//+
Line(7) = {6, 7};
//+
Line(8) = {7, 8};
//+
Line(9) = {8, 1};
//+
Curve Loop(1) = {2, 3, 4, 5, 7, 8, 9, 1};
//+
Plane Surface(1) = {1};
//+
ReverseMesh Surface{1};

So following the logic of gmsh the following goes wrong: gmsh provides clockwise or anti-clockwise counting based on your surface definition, which is clockwise and thus the elements are defined in this way as well.

Unfortunately, I don't know how we can check this beforehand and call the appropriate API function to reverse the mesh

[fredrikekre]

Ah, that makes sense I guess. Perhaps there is a way to check the normal of the surface or something?

[koehlerson]

yes, there is

help?> gmsh.model.get_normal
  gmsh.model.getNormal(tag, parametricCoord)

  Get the normal to the surface with tag tag at the parametric coordinates parametricCoord. parametricCoord are given
  by pairs of u and v coordinates, concatenated: [p1u, p1v, p2u, ...]. normals are returned as triplets of x, y, z
  components, concatenated: [n1x, n1y, n1z, n2x, ...].

  Return normals.

[fredrikekre]

Then we can probably check that and revert the order ourselves if needed?

[koehlerson]

There is also

help?> gmsh.model.mesh.reverse
  gmsh.model.mesh.reverse(dimTags = Tuple{Cint,Cint}[])

  Reverse the orientation of the elements in the entities dimTags. If dimTags is empty, reverse the orientation of the
  elements in the whole mesh.

But I think there is no easy accessible global criterion to check the definition. So probably best way is, as you say, to check it in the translate_elements function and revert it manually

[Drachta]

Could you try to reverse the mesh:

ReverseMesh Surface{1};

Thank you very much, this solved my problem!

Should I close the issue, or maybe you want to leave it open to automate the solution ?

[koehlerson]

leave it open, I will close it as soon as we have an automated check!

[DRollin]

Hi, I had a similar issue once regarding the ordering of vertices of a polyhedron and solved it by calculating the volume. The formula I used gives a negative volume in case of a clock-wise ordering. So, at least if all vertices are ordered in the same way, this might help. (I attached a PDF with the formula) Centroid_And_Volume_Of_A_Polyhedron.pdf

[koehlerson]

help?> gmsh.model.mesh.set_outward_orientation
  gmsh.model.mesh.setOutwardOrientation(tag)

  Set meshing constraints on the bounding surfaces of the volume of tag tag so that all surfaces are oriented with outward pointing normals;
  and if a mesh already exists, reorient it. Currently only available with the OpenCASCADE kernel, as it relies on the STL triangulation.

Calling this would be a no brainer, but I don't know how to obtain the tag of the bounding surface or more specifically how to distinguish any surface from the bounding surface

[DRollin]

I just re-read your message and noticed: gmsh.model.mesh.set_outward_orientation(tag) requires the tag of a volume, not its surface. So, calling the method for all volumes should be relatively easy, right? Or did I missunderstand something?