2D/3D serial and parallel triangular mesh generation for seismology
Project description
Create highquality 2D/3D meshes from seismic velocity models.
SeismicMesh: Mesh generation for Seismology in Python
SeismicMesh is a tool to generate 2D/3D triangular meshes that are used for acoustic and elastic wave propagators discretized with the finite element method. With short scripts, variable resolution meshes are built that have size transitions which reflect variations in Pwave or Swave velocities. Seismic Pwave/Swave data are typically provided on regular Cartesian grids for global and regional domains.
Generating variable resolution unstructured meshes of complex and largescale 2D/3D geophysical domains with popular mesh generation tools such as gmsh or cgal requires deciding how to size elements in the domain, which can become a complex task. Typically users must either create their own mesh sizing function that reflect the geophysical data in the domain or rely on analytical mesh sizing functions. However, for seismological problems with mesh size variations mostly in the interior of the domain, mesh sizing heuristics like boundary layer/attractor adaptation or characteristic size calculations may not be directly relevant. Instead, in a typical seismologial domain, variations in mesh size generally reflect internal material properties such as Pwave or Swave velocity, which are used to costeffectively model waves by ensuring there are sufficient number of gridpoints per wavelength.
Thus in this work we provide an approach to build large scale 2D/3D mesh sizing functions with the mesh sizing function module. This tool maps variations in seismic velocities from a seismic velocity model to triangular mesh sizes. Importantly, the sizing module can ensure that mesh size transitions vary smoothly (e.g., are graded) and an estimate of the Courant number can be bounded above by a constantamongst other capabilitieswhich are important considerations for accurate and successful simulations with finite elements. Our sizing functions are defined on the same regular Cartesian grids as the original seismic velocity model bypassing the need for the user to create their own sizing function on a triangular mesh. Their structured nature also enables efficient performance.
SeismicMesh supports both 2D and 3D triangular meshing in either serial or using distributed memory parallelism relying on the Message Passing Interface. In 3D mesh generation, a mesh improvement method (sliver removal) is used to ensure a minimum quality bound (e.g. minimum dihedral angle) can be enforced and will lead to numerically stable simulations.
Our mesh generation approach provided in this package can be operated standalone (e.g., without the sizing function module). It is based off modifications to the DistMesh algorithm. Thus in its most basic operation, SeismicMesh can mesh any domain that can be defined by a signed distance function with mesh sizes that follow variations described by a userdefined mesh sizing function
SeismicMesh is distributed under the GPL3 license and more details can be found in our short paper.
Installation
For installation, SeismicMesh needs CGAL and pybind11:
sudo apt install libcgaldev python3pybind11
After that, SeismicMesh can be installed from the Python Package Index (pypi), so with:
pip install U SeismicMesh
For more detailed information about installation and requirements see:
Install  How to install SeismicMesh.
Examples
The user can quickly build quality 2D/3D meshes from seismic velocity models in serial/parallel.
WARNING: To run the code snippet below you must download the 2D BP2004 seismic velocity model and then you must uncompress it (e.g., gunzip). This file can be downloaded from here
from mpi4py import MPI import meshio from SeismicMesh import get_sizing_function_from_segy, generate_mesh, Rectangle 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 SEGY 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) 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", )
WARNING: To run the code snippet below you must download (and uncompress) the 3D EAGE seismic velocity model from (WARNING: File is ~500 MB) here
WARNING: Computationaly demanding! Running this example takes around 3 minutes in serial and requires around 2 GB of RAM due to the 3D nature of the problem and the domain size.
from mpi4py import MPI import zipfile import meshio from SeismicMesh import ( get_sizing_function_from_segy, generate_mesh, sliver_removal, Cube, ) comm = MPI.COMM_WORLD # Bounding box describing domain extents (corner coordinates) bbox = (4200.0, 0.0, 0.0, 13520.0, 0.0, 13520.0) # Desired minimum mesh size in domain. hmin = 150.0 # This file is in a big Endian binary format, so we must tell the program the shape of the velocity model. path = "Salt_Model_3D/3D_Salt_Model/VEL_GRIDS/" if comm.rank == 0: # Extract binary file Saltf@@ from SALTF.ZIP zipfile.ZipFile(path + "SALTF.ZIP", "r").extract("Saltf@@", path=path) fname = path + "Saltf@@" # Dimensions of model (number of grid points in z, x, and y) nx, ny, nz = 676, 676, 210 cube = Cube(bbox) # A graded sizing function is created from the velocity model along with a signed distance function by passing # the velocity grid that we created above. # More details can be found here: https://seismicmesh.readthedocs.io/en/master/api.html ef = get_sizing_function_from_segy( fname, bbox, hmin=hmin, dt=0.001, freq=2, wl=5, grade=0.15, hmax=5e3, domain_pad=250, pad_style="linear_ramp", nz=nz, nx=nx, ny=ny, byte_order="big", axes_order=(2, 0, 1), # order for EAGE (x, y, z) to default order (z,x,y) axes_order_sort="F", # binary is packed in a FORTRANstyle ) points, cells = generate_mesh(domain=cube, h0=hmin, edge_length=ef, max_iter=75) # For 3D mesh generation, we provide an implementation to bound the minimum dihedral angle:: points, cells = sliver_removal( points=points, bbox=bbox, h0=hmin, domain=cube, edge_length=ef ) # Meshes can be written quickly to disk using meshio and visualized with ParaView:: if comm.rank == 0: # NOTE: SeismicMesh outputs assumes the domain is (z,x,y) so for visualization # in ParaView, we swap the axes so it appears as in the (x,y,z) plane. meshio.write_points_cells( "EAGE_Salt.vtk", points[:, [1, 2, 0]] / 1000.0, [("tetra", cells)], )
The user can still specify their own signed distance functions and sizing functions to generate_mesh
(in serial or parallel) just like the original DistMesh algorithm. Try the codes below!
# Mesh a unit cylinder from mpi4py import MPI from numpy import array, maximum, sqrt, zeros_like import meshio from SeismicMesh import generate_mesh, sliver_removal comm = MPI.COMM_WORLD hmin = 0.10 bbox = (1.0, 1.0, 1.0, 1.0, 1.0, 1.0) def cylinder(p): r, z = sqrt(p[:, 0] ** 2 + p[:, 1] ** 2), p[:, 2] d1, d2, d3 = r  1.0, z  1.0, z  1.0 d4, d5 = sqrt(d1 ** 2 + d2 ** 2), sqrt(d1 ** 2 + d3 ** 2) d = maximum.reduce([d1, d2, d3]) ix = (d1 > 0) * (d2 > 0) d[ix] = d4[ix] ix = (d1 > 0) * (d3 > 0) d[ix] = d5[ix] return d def fh(p): # Note: for parallel execution this logic is required # since the decomposition of the sizing function passes a tuple to fh if type(p) == tuple: h = zeros_like(p[0]) + hmin else: h = array([hmin] * len(p)) return h points, cells = generate_mesh( bbox=bbox, domain=cylinder, h0=hmin, edge_length=fh, max_iter=100, ) points, cells = sliver_removal( points=points, domain=cylinder, edge_length=fh, h0=hmin, min_dh_angle_bound=5.0, bbox=bbox, ) if comm.rank == 0: meshio.write_points_cells( "Cylinder.vtk", points, [("tetra", cells)], file_format="vtk", )
import meshio import SeismicMesh disk = SeismicMesh.Disk([0.0, 0.0], 1.0) points, cells = SeismicMesh.generate_mesh(domain=disk, edge_length=0.1) meshio.write_points_cells( "disk.vtk", points, [("triangle", cells)], file_format="vtk", )
import meshio import SeismicMesh bbox = (0.0, 1.0, 0.0, 1.0, 0.0, 1.0) cube = SeismicMesh.Cube(bbox) points, cells = SeismicMesh.generate_mesh(domain=cube, edge_length=0.05) points, cells = SeismicMesh.sliver_removal(points=points, domain=cube, edge_length=0.05) meshio.write_points_cells( "cube.vtk", points, [("tetra", cells)], file_format="vtk", )
import meshio import SeismicMesh bbox = (0.0, 1.0, 0.0, 1.0, 0.0, 1.0) cube = SeismicMesh.Cube(bbox) points, cells = SeismicMesh.generate_mesh(domain=cube, edge_length=0.05) meshio.write_points_cells( "cube.vtk", points, [("tetra", cells)], file_format="vtk", )
# mesh a torus import numpy as np import meshio import SeismicMesh bbox = (1.0, 1.0, 1.0, 1.0, 10.0, 1.0) hmin = 0.05 def length(x): return np.sum(np.abs(x) ** 2, axis=1) ** (1.0 / 2) def Torus(p, t=(0.5, 0.2)): xz = np.column_stack((p[:, 0], p[:, 2])) q = np.column_stack((length(xz)  t[0], p[:, 1])) return length(q)  t[1] points, cells = SeismicMesh.generate_mesh( domain=Torus, edge_length=hmin, bbox=bbox, verbose=2, max_iter=100 ) points, cells = SeismicMesh.sliver_removal( points=points, domain=Torus, edge_length=hmin, bbox=bbox, verbose=2 ) meshio.write_points_cells( "torus.vtk", points, [("tetra", cells)], )
# mesh a prism import numpy as np import meshio import SeismicMesh bbox = (1.0, 1.0, 1.0, 1.0, 1.0, 1.0) hmin = 0.05 def length(x): return np.sum(np.abs(x) ** 2, axis=1) ** (1.0 / 2) def sdTriPrism(p, h=[0.5, 0.5]): q = np.abs(p) return np.maximum( q[:, 2]  h[1], np.maximum(q[:, 0] * 0.866025 + p[:, 1] * 0.5, p[:, 1])  h[0] * 0.5, ) points, cells = SeismicMesh.generate_mesh( bbox=bbox, domain=sdTriPrism, edge_length=hmin, verbose=2, max_iter=100, ) points, cells = SeismicMesh.sliver_removal( bbox=bbox, points=points, domain=sdTriPrism, edge_length=hmin ) meshio.write_points_cells( "prism.vtk", points, [("tetra", cells)], file_format="vtk", )
# Compute the union of several SDFs to create more complex geometries import meshio import SeismicMesh h = 0.10 rect0 = SeismicMesh.Rectangle((0.0, 1.0, 0.0, 0.5)) rect1 = SeismicMesh.Rectangle((0.0, 0.5, 0.0, 1.0)) disk0 = SeismicMesh.Disk([0.5, 0.5], 0.5) union = SeismicMesh.Union([rect0, rect1, disk0]) points, cells = SeismicMesh.generate_mesh(domain=union, edge_length=h) meshio.write_points_cells( "Lshape_wDisk.vtk", points, [("triangle", cells)], file_format="vtk", )
# Compute the intersection of several SDFs to create more complex geometries import meshio import SeismicMesh h = 0.05 rect0 = SeismicMesh.Rectangle((0.0, 1.0, 0.0, 1.0)) disk0 = SeismicMesh.Disk([0.25, 0.25], 0.5) disk1 = SeismicMesh.Disk([0.75, 0.75], 0.5) intersection = SeismicMesh.Intersection([rect0, disk0, disk1]) points, cells = SeismicMesh.generate_mesh(domain=intersection, edge_length=h) meshio.write_points_cells( "Leaf.vtk", points, [("triangle", cells)], file_format="vtk", )
# Compute the difference of two SDFs to create more complex geometries. import meshio import SeismicMesh h = 0.05 rect0 = SeismicMesh.Rectangle((0.0, 1.0, 0.0, 1.0)) disk0 = SeismicMesh.Disk([0.5, 0.5], 0.1) disk1 = SeismicMesh.Disk([0.75, 0.75], 0.20) difference = SeismicMesh.Difference([rect0, disk0, disk1]) points, cells = SeismicMesh.generate_mesh(domain=difference, edge_length=h) meshio.write_points_cells( "Hole.vtk", points, [("triangle", cells)], file_format="vtk", )
import meshio import SeismicMesh h = 0.10 cube0 = SeismicMesh.Cube((0.0, 1.0, 0.0, 1.0, 0.0, 1.0)) ball1 = SeismicMesh.Ball([0.5, 0.0, 0.5], 0.30) ball2 = SeismicMesh.Ball([0.5, 0.5, 0.0], 0.30) ball3 = SeismicMesh.Ball([0.0, 0.5, 0.5], 0.30) ball4 = SeismicMesh.Ball([0.5, 0.5, 0.5], 0.45) difference = SeismicMesh.Difference([cube0, ball1, ball2, ball3, ball4]) points, cells = SeismicMesh.generate_mesh(domain=difference, edge_length=h, verbose=1) points, cells = SeismicMesh.sliver_removal( points=points, domain=difference, edge_length=h, verbose=1 ) meshio.write_points_cells( "Ball_wHoles.vtk", points, [("tetra", cells)], file_format="vtk", )
How does performance and cell quality compare to gmsh
and cgal
mesh generators?
Here we use SeismicMesh 3.0.4, pygalmesh 0.8.2, and pygmsh 7.0.0 (more details in the benchmarks folder).
Some key findings:
 Mesh generation in 2D and 3D using analytical sizing functions is quickest when using
gmsh
followed bycgal
and thenSeismicMesh
.  However, using mesh sizing functions defined on gridded interpolants significantly slow down both
gmsh
andcgal
. In these cases,SeismicMesh
andgmsh
perform similarly both outperformingcgal
's 3D mesh generator in terms of mesh generation time. SeismicMesh
produces often comparable or higher mean cell qualities in 2D/3D than eithergmsh
orcgal
and this may have implications on mesh improvement strategies as higher minimum quality may be realizable with some common mesh improvement strategies (e.g., NetGen) All methods produce 3D triangulations that have a minimum dihedral angle > 10 degrees enabling stable numerical simulation.
 Head over to the
benchmarks
folder for more detailed information on these experiments.

In the figure for the panels that show cell quality, solid lines indicate the mean and dashed lines indicate the minimum cell quality in the mesh.

Note: it's important to point out here that a significant speedup can be achieved for moderate to large problems using the parallel capabilities provided in
SeismicMesh
.
Changelog
The format is based on Keep a Changelog, and this project (tries to) adhere to Semantic Versioning.
Unreleased
[3.0.6]  20201021
Fixed
 Silence messages about pfix when verbose=0
Added
 Added more examples on README
 New unions/intersections/differences with several SDF primivitives
 Automatic corner constraints in serial
[3.0.5]  20201018
Fixed
 Preserving fixed points in serial.
 Units in kms detection warning bug.
 Docstring fixes to
generate_mesh
 Improved mesh quality in 3D
Added
 Automatic corner point extraction for cubes and rectangles.
 More support for reading binary files packed in a binary format.
 Check to make sure bbox is composed of all floats.
[3.0.4]  20201012
Added
 Improve conformity of levelset in final mesh through additional set of Newton boundary projection iterations.
More information
All other information is available at: https://seismicmesh.readthedocs.io
Getting started  Learn the basics about the program and the application domain.
Tutorials  Tutorials that will guide you through the main features.
Project details
Release history Release notifications  RSS feed
Download files
Download the file for your platform. If you're not sure which to choose, learn more about installing packages.
Filename, size  File type  Python version  Upload date  Hashes 

Filename, size SeismicMesh3.0.6cp37cp37mmacosx_10_15_x86_64.whl (2.3 MB)  File type Wheel  Python version cp37  Upload date  Hashes View 
Filename, size SeismicMesh3.0.6.tar.gz (6.0 MB)  File type Source  Python version None  Upload date  Hashes View 
Hashes for SeismicMesh3.0.6cp37cp37mmacosx_10_15_x86_64.whl
Algorithm  Hash digest  

SHA256  87c242300c8719faa07328b95966421601e65fd13a7b7ea0755a3ddc721444ed 

MD5  5e68e98cbc2664c8e84f063d7777a0f1 

BLAKE2256  e087e3f85922575630a5ee789b7a258688f82cf3e3ab323f5d5292dbf622218a 