make purpose of SeismicMesh clearer

[nschloe]

This is about the README. The premise of SeismicMesh is unclear:

2D/3D triangular meshing for a slab of Earth based on modifications to the DistMesh algorithm. SeismicMesh is distributed under GPL3.

First off, I think it's tetrahedral meshing in 3D, right? What does is "meshing for a slab of Earth" mean, how is it different from meshing for the air around a wing? What exactly does SeismicMesh do differently than gmsh, or what's it's specialty? I'm guessing it has to do with the segy files. How is this useful for seismic computations?

[krober10nd]

It's both 2D and 3D quality triangular meshing in serial and distributed memory parallel. I can make the purpose a bit clearer of course on the README but I thought this was the purpose of the paper.

[nschloe]

The paper should definitely be linked in the readme, but the first paragraph in the readme should definitely make clear what SeismicMesh is and what it can do. How is it better than the existing mesh generators out there, why did you write it?

[nschloe]

I personally am curious what makes SeismicMesh particular to seismic problems.

[krober10nd]

Sure, okay I agree in can be made clearer. I'll work on that.

Here are answers in regard to existing mesh generations like gmsh vs. `SeismicMesh:

• The origin behind SeismicMesh arose from a project-wide need at my university to automatically and approximately reproducibly script the generation of 2D/3D quality meshes intended for acoustic and elastic FEM wave propagators directly from seismic velocity models (e.g., scalar fields defined on relatively large structured grids). We needed a standard and open source tool everyone could use that perhaps just knew a little Python and weren't very familiar with mesh generation tools. These considerations largely deterred us from gmsh or cgal.

• Following that, we wanted a tool that doesn't require you to learn a whole new meta-programming language (as is the case with gmsh, cgal). Users run a simple Python script, change a handful of plain-text options, and it directly outputs 2D/3D meshes that are adapted to variations in the seismic velocity model and can be reasonably ensured numerically stable through some provided mesh sizing function capabilities. This greatly reduces potential simulation problems using the generated mesh. As an aside, if you want mesh refinement in the interior of the domain with gmsh or cgal or tetgen, you must define your own sizing function. However, each user may do that differently so that didn't bode well for our users too well.

• Compared to gmsh and tetgen, anecdotally DistMesh tends to often produce far higher quality meshes given enough meshing iterations especially when combined with some mesh optimization strategy like Laplacian smoothing or sliver removal (in fact you even allude to this point here: https://github.com/nschloe/optimesh). To this end, the parallelization here accelerates those meshing iterations leading to perhaps higher quality meshes faster. Perhaps the same parallelization could be used with optmimesh, I don't know. Further, the sliver removal algorithm I programmed in here bears strong resemblance to the one used in cgal but isn't buried in thousands of lines of cpp. In the future, I actually plan to enable users to just pass a 3D mesh and it will remove those slivers if they exist. To my knowledge, there are currently no simple 3D sliver removal tools out there in Python.

• Besides this, DistMesh remains a popular algorithm to build meshes for a variety of application and has merit on its own as a mesh generator. The big downside of DistMesh has always been it is serial, quite slow because of the re-triangulation call each iteration, was filled with memory bottlenecks , and in 3D slivers appeared ruining numerical simulations. I have parallelized this algorithm and avoided those problems in this implementation, which I think is a contribution. In fact, I specifically designed the meshgen class to operate standalone-- just like the original DistMesh only requiring a signed distance function and sizing function (e.g., https://github.com/krober10nd/SeismicMesh/blob/par3d/examples/example_3D_SDF2.py). Furthermore, in an effort to future proof it and to enable other geophysical applications like OceanMesh2Dhttps://github.com/CHLNDDEV/OceanMesh2D_ which I'm currently porting to Python pyOceanMesh https://github.com/CHLNDDEV/pyOceanMesh to use it in a plug-and-go type fashion.

Here are my comments on your question about your query why SeismicMesh is particular to seismic problems:

• The specific connection to seismology is the provided mesh sizing function module, which maps variations in seismic velocity to triangular mesh sizes and does it in a way that the mesh size transitions are graded and the Courant number can be bounded above by a constant. In seismic applications, the interior of the domain must be adapted so that the wavelengths of a source frequency (a controlled explosion) are well-resolved. But the problem is the transitions of the material properties are very rapid so the problem requires mesh gradation and some checks to make sure the CFL limit isn't violated (since we almost exclusively rely on explicit schemes). Further, users often want domain extensions and have to apply a special treatment of the mesh sizes and material properties in this extension subregion.

• For seismic domains with irregular surface boundaries (in 2D or 3D) like surveys over land, the geometry information describing the free surface boundary doesn't always exist and thus has to be created. This can be a big pain and bottleneck for automation. But the good news is these boundaries are represented by thresholding material properties (like density or seismic velocity). SeismicMesh can threshold these fields and build signed distance functions that can then be used to build meshes domains with.

Bottom line: the are many general purpose use cases of SeismicMesh besides in seismology but currently that's what it's used for.

[krober10nd]

I took a stab at improving the purpose and clarity in #45, which can be seen here https://github.com/krober10nd/SeismicMesh/tree/purpose

[krober10nd]

I did a couple more revisions of the README and merged it in. I personally think it's clearer now.

[krober10nd]

Revised it a little more today.

[nschloe]

These aforementioned mesh generation software generally require a sizing function to be defined on a mesh, which creates a circular problem of building a mesh to build a mesh.

That's not true. CGAL does mesh refinement based on general functions (see, e.g., https://github.com/nschloe/pygalmesh#local-refinement) and so does gmsh.

So if I get it right, SeismicMesh is basically a big distmesh implementation that can also ready seismic velocity model files and create a mesh sizing field from it. Is that correct?

[krober10nd]

That's not true. CGAL does mesh refinement based on general functions (see, e.g., https://github.com/nschloe/pygalmesh#local-refinement) and so does gmsh.

That's fair and I'll edit the language to reflect your point. However, I meant more along the lines of how does a user build a sizing function from a given seismic velocity model.

For analytical sizing functions, it's straightforward. For user-defined fields, the gmsh 4.6.0 user guide syas:

The PostView field specifies an explicit background mesh in the form of a scalar post-processing view (see Post-processing commands, and File formats) in which the nodal values are the target element sizes. This method is very general but it requires a first (usually rough) mesh and a way to compute the target sizes on this mesh (usually through an error estimation procedure, in an iterative process of mesh adaptation). Warning: only parsed (.pos) files can currently be used as background meshes (.msh files cannot be used, since the mesh used to define the field will be destroyed during the meshing process). (Note that you can also load a background mesh directly from the command line using the -bgm option (see Command-line options), or in the GUI by selecting ‘Apply as background mesh’ in the post-processing view option menu.)

Hence where my comment a mesh to define a mesh.

So if I get it right, SeismicMesh is basically a big distmesh implementation that can also ready seismic velocity model files and create a mesh sizing field from it. Is that correct?

That's correct. A parallel distmesh with mesh improvement and mesh sizing field creation from seismic velocity model files.

[nschloe]

That's correct.

Good. Can SeismicMesh generate meshes for non-rectangular domains, too?

[krober10nd]

Good. Can SeismicMesh generate meshes for non-rectangular domains, too?

Yep.

There's two examples in the folder examples that show how you can mesh non-rectangular domains (2D and 3D).

I use iso contours of seismic velocity to define a 0-level set, which defines the meshing domain.