Minimal mesh quality of DistMesh algorithm

[jorgensd]

As I've mentioned in previous comments, it seems like the DistMesh algorithm produces cells where the average mesh quality is very high. However, it seems like the minimal quality of cells in a mesh is quite alot worse than for GMSH. In the setting of 3D you have shown that the sliver removal algorithm can increase the mesh quality of these cells quite efficiently.

However, in the 2D setting, how can one control the minimum mesh quality? Take for instance the BP2004 benchmark. Running the benchmark produces the following minimal, avg and max mesh quality:

1.466166516027528066e-01
9.593859657712772160e-01
9.999999866475057786e-01

while gmsh produces:

5.564115344181080891e-01
9.413313063736762354e-01
9.999983294113985455e-01

I've attached the figures of the worst cell in gmsh and seismic mesh below, where the color indicate the radius ratio (Seismic Mesh 6.82 on the scale [1,infty), Gmsh 1.79 on the same scale). Allthough the cell is far from collapsed, I think it is very interesting in the perspective of FEM and multigrid solvers.

GMSH:

SeismicMesh:

[jorgensd]

As a side note, I ran then BP2004 benchmark on my Desktop computer with 64 GB RAM (Intel® Core™ i9-9900K CPU @ 3.60GHz × 16 ) and obtained the following timings: Seismic mesh number of cells and timings where:

[85557, 163342, 272707, 451197, 758263]
[2.9803144931793213, 5.7315192222595215, 9.913894414901733, 16.317556619644165, 27.86741805076599]

while GMSH

[113184, 216634, 360605, 596192, 1003911]
[8.692016124725342, 15.265772342681885, 24.70887804031372, 40.25092577934265, 68.42505860328674]

confirming a similar trend to what you have shown in your studies (just with quite a reduced runtime).

[krober10nd]

Hey Jorgen,

So in regard to minimal cell qualities in 2D, they almost exclusively occur on the boundary of the domain and there's typically less than 10 or so interspersed throughout the domain.

I either delete them https://github.com/krober10nd/SeismicMesh/blob/b45ac8dd9ee3368bf4a701073e1aa3876ba7e7d8/SeismicMesh/geometry/utils.py#L434

or/and I run Laplacian smoothing https://github.com/krober10nd/SeismicMesh/blob/b45ac8dd9ee3368bf4a701073e1aa3876ba7e7d8/SeismicMesh/geometry/utils.py#L487

running the example with Laplacian smoothing at the end produces a minimal cell quality of 0.53/1 for me

 from mpi4py import MPI
 import numpy as np
 import meshio

 from SeismicMesh import (
     get_sizing_function_from_segy,
     generate_mesh,
     Rectangle,
     geometry,
 )

 comm = MPI.COMM_WORLD

 """
 Build a mesh of the BP2004 benchmark velocity model in serial or parallel
 Takes roughly 1 minute with 2 processors and less than 1 GB of RAM.
 """

 # Name of SEG-Y file containg velocity model.
 fname = "vel_z6.25m_x12.5m_exact.segy"

 # Bounding box describing domain extents (corner coordinates)
 bbox = (-12000.0, 0.0, 0.0, 67000.0)

 # Desired minimum mesh size in domain
 hmin = 75.0

 rectangle = Rectangle(bbox)

 # Construct mesh sizing object from velocity model
 ef = get_sizing_function_from_segy(
     fname,
     bbox,
     hmin=hmin,
     wl=10,
     freq=2,
     dt=0.001,
     grade=0.15,
     domain_pad=1e3,
     pad_style="edge",
 )

 points, cells = generate_mesh(domain=rectangle, edge_length=ef, h0=hmin)

 points, cells = geometry.laplacian2(points, cells)
 print(np.amin(geometry.simp_qual(points, cells)))

 if comm.rank == 0:
     # Write the mesh in a vtk format for visualization in ParaView
     # NOTE: SeismicMesh outputs assumes the domain is (z,x) so for visualization
     # in ParaView, we swap the axes so it appears as in the (x,z) plane.
     meshio.write_points_cells(
         "BP2004.vtk",
         points[:, [1, 0]] / 1000,
         [("triangle", cells)],
         file_format="vtk",
     )

[krober10nd]

Good to hear about the timings! I'm running on a 2015 MacBook Pro that's been through quite a bit of debugging and mesh generation in its lifetime so I'd expect the results to be better on newer machines.

Parallel speed-ups in 2D are better as well (regardless of the machine).

[jorgensd]

I would probably remove CGAL from the EAGE plot, as it hdes any difference between GMSH and SeismicMesh, and isnt really competitive for this example. What do you think at @nschloe? I think you could just mention it in the text, and keep the bench as a standalone function people could try if they want to

[krober10nd]

Yes, it does make the difference between gmsh and SM difficult to see.

[jorgensd]

I added the laplace2 (in addition to the delete_boundary_entities) to the benchbp2004 bench, and Im attaching the results: I would probably include the laplace smoothing, as it puts SeismicMesh closer to gmsh wrt. minimal mesh quality (and even better when it comes to max quality).

[krober10nd]

Ah neat. Yea I can make it a user option inside the 'mesh_generate' by default for 2d. It's fairly quick.

Wonder what happened in the drop out case.

[jorgensd]

Eage benchmark without CGAL Here we observe a quite different trend than in the 2D case. Any suggestions as to why?

[krober10nd]

I'll inspect more. But it's likely the unique edges calculation is more expensive in 3d. The other factor is that 3d sliver removal takes a couple seconds.

Perhaps Gmsh is more optimized in 3d as well.

[krober10nd]

I added the laplace2 (in addition to the delete_boundary_entities) to the benchbp2004 bench, and Im attaching the results:

Just curious, where did you put the delete_boundary_entities step that produced this figure?

[jorgensd]

I added the laplace2 (in addition to the delete_boundary_entities) to the benchbp2004 bench, and Im attaching the results:

Just curious, where did you put the delete_boundary_entities step that produced this figure?

It was placed before the laplace smoothing

[krober10nd]

In regard to the EAGE benchmark (sans CGAL) this is what I get with the laplacian branch.

The primary slow down in 3d for SM is sliver removal.

what's your gmsh and pygmsh versions? I'm using pygmsh 7.0.0 using Gmsh 4.7.0.

[krober10nd]

I remade the benchmark figure #145 with 3.1.4.

I'll update the paper with the new conclusions from the figure in the ongoing #140

[jorgensd]

Very tiny detail (feel free to deny changing it. In the BP benchmark, I would suggest keeping the y-axis from 0 to 1, to illustrate that both GMSH and SeismicMesh produce very high quality cells (and it would align nicely with the other figures).

[krober10nd]

Okay this is addressed in #147 I redid the axis for the disk as well as BP2004.