Is mesh distortion okay?

[gatling-nrl]

I'm working on a 2D problem in the unit square, where I've defined several disk subdomains. I'm running this against thousands of arrangements of the disks, each time generating the mesh with gmsh/occ and it is sometimes complaining, as in

Warning : Surface mesh: worst distortion = -0.395997 (avg = 0.998517, 1 elements with jac. < 0); worst gamma = 0.0309604
Warning : Surface mesh: worst distortion = -1.31986 (avg = 0.999582, 1 elements with jac. < 0); worst gamma = 0.00291112
Warning : Surface mesh: worst distortion = -0.349159 (avg = 0.986488, 2 elements with jac. < 0); worst gamma = 0.0273358
Warning : Surface mesh: worst distortion = -1.20468 (avg = 0.994867, 2 elements with jac. < 0); worst gamma = 0.0273358
Warning : Surface mesh: worst distortion = -1.5496 (avg = 0.998435, 2 elements with jac. < 0); worst gamma = 0.000741227
Warning : Surface mesh: worst distortion = -24.3248 (avg = 0.99927, 2 elements with jac. < 0); worst gamma = 3.39695e-05
Warning : Surface mesh: worst distortion = -0.471796 (avg = 0.987441, 1 elements with jac. < 0); worst gamma = 0.0324295

However skfem seems "okay" with these meshes. I'm getting solutions that "look" right. I don't know that much about meshing, so I'm not sure how to interpret these warnings in the context of FEM and google isn't much help at all here.

I realize this question is 1/2 gmsh and 1/2 skfem, but I figured I would start here... do these warnings matter? Any advice how to get rid of them?

[gatling-nrl]

One simple way to stop the warnings has been to add gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 25).

From the gmsh docs:

Second, if Mesh.MeshSizeFromCurvature is set to a positive value (it is set to 0 by default), the mesh will be adapted with respect to the curvature of the model entities, the value giving the target number of elements per 2 Pi radians.

But I would still like to know how the warning impacts skfem.

[gdmcbain]

If it doesn't take too much longer, i would be issuing this option for sure. I expect it will pay off in accuracy and stability.

[kinnala]

Possible to draw one of the meshes? Triangular P1 element? P1 FEM is pretty robust against bad meshes and it is actually really hard to find mesh families where FEM is not converging.

[kinnala]

Cannot see anything suspicious although it's hard to see the small elements. If you want to investigate further, you can export basis.mapping.detA to VTK file and zoom in using Paraview. Command is something like m.save('file.vtk', cell_data={'detA': basis.mapping.detA). Make sure that whatever you put to cell_data is 1d numpy array with one value per element.

[kinnala]

Btw, you have multiple triangle meshes inside the meshio object. How does that look after skfem.io.from_meshio? Does it successfully read all triangles? I have not tried such mesh I think.

[gdmcbain]

meshio>=4 stores cells as a list of homogeneous blocks of cells. We adapted to this promptly.

The cell_data needs to have the same shape; I.e. the values of the dict should be lists of values.

When skfem exports to meshio, it just puts all the domain cells in one block; Gmsh (MSH 4) has some internal rationale for grouping the cells into blocks.

[gatling-nrl]

@kinnala it definitely reads all the triangles.

[gdmcbain]

Yep, this is because of

scikit-fem/skfem/io/meshio.py

Line 50 in bdaf712

cells = m.cells_dict

The meshio.Mesh.cell_dict method regroups all the cells in the list of cell locks by type. It returns a dictionary which will only have one entry for three-noded triangles.

[gatling-nrl]

To turn it high enough to get rid of all of the warnings made the mesh very dense and slow. It turned out to be better to specify a background field according to what I new about the geometry. See gmsh tutorial #10. The gmsh docs are hard to read, but the section on the different field types is important, and combining several together with a min field eliminated all of the warnings.

Note that if you are going after a disk (as I was) adding a temporary point at the center and using the distance field against that point is more effective than trying to use the circle perimiter of the disk or the area of the disk itself. See t10 for example. Make sure not to leave the temporary point in the mesh as an orphan!

[gdmcbain]

Yes. Gmsh's background fields are powerful, offering a lot of control.

Another approach to getting a good mesh is adaptive mesh refinement. This can be useful particularly if it's not obvious beforehand where the solution will need a finer mesh for resolution. Also it's more automatic and optimal. See ex22: Adaptive Poisson equation.

[gatling-nrl]

Does this work starting with an arbitrary meshio mesh? I.e., I create a coarse mesh useful to describe the geometry of the domain (regions of constant conductivity in my current application) and then apply ex22 to that mesh to adaptively refine it? For example, starting mesh might look like one of those posted earlier in this discussion.

[gdmcbain]

I have done that kind of thing in the past, yes, begin with a coarse mesh and have the algorithm tailor a mesh to the problem. This was particularly effective in a problem of forced convection for which the variation in temperature was largely confinely to thin wakes.

A difficulty with adaptive mesh refinement in the original example here is that the boundaries aren't polygonal and if they're initially defined as polygons with few edges, extra work will need to be done to get the refined mesh to respect the underlying curvilinear geometry. There are ways to do this, as in MeshTri1.init_circle:

scikit-fem/skfem/mesh/mesh_tri_1.py

Line 196 in bdaf712

tmp[:, D] = tmp[:, D] / np.linalg.norm(tmp[:, D], axis=0)