Skip to content

2D/3D serial and parallel triangular mesh generation tool for Finite Element Methods that solve acoustic/elastic wave equations.

License

Notifications You must be signed in to change notification settings

Shwe-creator/SeismicMesh

 
 

Repository files navigation

SeismicMesh

Create high-quality 2D/3D meshes from seismic velocity models.

CircleCI CodeCov Code style: black PyPI pyversions PyPi downloads ReadTheDocs Zenodo PyPi GPL status

SeismicMesh: Triangular Mesh generation in Python

SeismicMesh is a tool to generate unstructured and high geometric quality 2D/3D triangular meshes that can be used for acoustic and elastic wave propagators discretized with the finite element method. SeismicMesh however can also be used as a general purpose simplicial mesh generation tool in 2D/3D. With short scripts, variable resolution meshes are built that have size transitions which reflect variations in P-wave or S-wave velocities from seismic velocity models. Seismic P-wave/S-wave data are typically provided on regular Cartesian grids for global and regional domains.

Generating variable resolution unstructured meshes of complex and large-scale 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 P-wave or S-wave velocity, which are used to cost-effectively model waves by ensuring there are sufficient number of grid-points 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 constant--amongst other capabilities--which 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 user-defined 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 libcgal-dev python3-pybind11

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.

Contributing

All contributions are welcome!

To contribute to the software:

  1. Fork the repository.
  2. Clone the forked repository, add your contributions and push the changes to your fork.
  3. Create a Pull request

Before creating the pull request, make sure that the tests pass by running

tox

Some things that will increase the chance that your pull request is accepted:

  • Write tests.
  • Add Python docstrings that follow the Sphinx.
  • Write good commit and pull request messages.

Problems?

If something isn't working as it should or you'd like to recommend a new addition/feature to the software, please let me know by starting an issue through the issues tab. I'll try to get to it as soon as possible.

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

Above shows the mesh in ParaView that results from running the code below

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 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)

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.

Above shows the mesh in ParaView that results from running the code below.

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/3-D_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 FORTRAN-style
)

points, cells = generate_mesh(domain=cube, 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, 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 but now with quality bounds in 3D. Try the codes below!

Cylinder

# Mesh a cylinder
from mpi4py import MPI
import meshio

import SeismicMesh

comm = MPI.COMM_WORLD

hmin = 0.10

cylinder = SeismicMesh.Cylinder(h=1.0, r=0.5)

points, cells = SeismicMesh.generate_mesh(
    domain=cylinder,
    edge_length=hmin,
)

points, cells = SeismicMesh.sliver_removal(
    points=points,
    domain=cylinder,
    edge_length=hmin,
)

if comm.rank == 0:
    meshio.write_points_cells(
        "Cylinder.vtk",
        points,
        [("tetra", cells)],
        file_format="vtk",
    )

Disk

# mesh a disk
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",
)

Square

# mesh a square/rectangle
import meshio
import SeismicMesh

bbox = (0.0, 1.0, 0.0, 1.0)
square = SeismicMesh.Rectangle(bbox)
points, cells = SeismicMesh.generate_mesh(domain=square, edge_length=0.05)
meshio.write_points_cells(
    "square.vtk",
    points,
    [("triangle", cells)],
    file_format="vtk",
)

Cube

# mesh a cuboid/cube
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",
)

Torus

# mesh a torus
import meshio
import SeismicMesh

hmin = 0.10

torus = SeismicMesh.Torus(r1=1.0, r2=0.5)
points, cells = SeismicMesh.generate_mesh(
    domain=torus,
    edge_length=hmin,
)
points, cells = SeismicMesh.sliver_removal(
    points=points, domain=torus, edge_length=hmin
)
meshio.write_points_cells(
    "torus.vtk",
    points,
    [("tetra", cells)],
)

Torus

# mesh a prism
import meshio

import SeismicMesh

hmin = 0.05
prism = SeismicMesh.Prism(b=0.5, h=0.5)

points, cells = SeismicMesh.generate_mesh(
    domain=prism,
    edge_length=hmin,
)
points, cells = SeismicMesh.sliver_removal(
    points=points, domain=prism, edge_length=hmin
)
meshio.write_points_cells(
    "prism.vtk",
    points,
    [("tetra", cells)],
    file_format="vtk",
)

Union

# 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",
)

Leaf

# 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",
)

Hole

# 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",
)

Cube wHoles

# Compute the difference of several SDFs in 3D
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(
    "Cube_wHoles.vtk",
    points,
    [("tetra", cells)],
    file_format="vtk",
)

Immersed disk

# Immerse a subdomain so that it's boundary is conforming in the mesh.
import numpy as np

import meshio

import SeismicMesh

box0 = SeismicMesh.Rectangle((-1.25, 0.0, -0.250, 1.250))
disk0 = SeismicMesh.Disk([-0.5, 0.5], 0.25)

hmin = 0.10


fh = lambda p: 0.05 * np.abs(disk0.eval(p)) + hmin

points, cells = SeismicMesh.generate_mesh(
    domain=box0,
    edge_length=fh,
    h0=hmin,
    subdomains=[disk0],
    max_iter=100,
)
meshio.write_points_cells(
    "Square_wsubdomain.vtk",
    points,
    [("triangle", cells)],
    file_format="vtk",
)

How does performance and cell quality compare to Gmsh and CGAL mesh generators?

Here we use SeismicMesh 3.1.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 but a closer competition for CGAL and SeismicMesh.
  • However, using mesh sizing functions defined on gridded interpolants significantly slow down both Gmsh and CGAL. In these cases, SeismicMesh and Gmsh perform similarly both outperforming CGAL's 3D mesh generator in terms of mesh generation time.
  • All methods produce 3D triangulations that have a minimum dihedral angle > 10 degrees enabling stable numerical simulation (not shown)
  • Head over to the benchmarks folder for more detailed information on these experiments.

Summary of the benchmarks

  • 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 speed-up 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 adheres to Semantic Versioning.

Unreleased

  • Support for constraining/immersing subdomains represented as signed distance functions.

[3.1.4] - 2020-11-15

  • Laplacian smoothing at termination for 2D meshing...significantly improves minimum cell quality.
  • Made hmin a field of the SizeFunction class, which implies the user no longer needs to pass h0 to generate_mesh or sliver_removal.

[3.1.3] - 2020-11-06

Fixed

  • Cylinder radius and height are now correct.
  • Torus, Prism, and Cylinder now have dim tag.

Improved

  • More control over the grad option in the mesh sizing function.

[3.1.2] - 2020-11-04

Improved

  • Faster calculation of boundary vertices.
  • More robust sliver removal in 3D.

Fixed

  • Corners are only constrained for constant resolution meshes

[3.1.0] - 2020-10-28

Added

  • New geometric primitives--torus, wedge/prism, and cylinder.
  • Updated images on README.

Fixed

  • Only constrain corners near 0-level set.
  • Bug fix to 3D binary velocity reading.

[3.0.6] - 2020-10-21

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] - 2020-10-18

Fixed

  • Preserving fixed points in serial.
  • Units in km-s 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] - 2020-10-12

Added

  • Improve conformity of level-set 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.

About

2D/3D serial and parallel triangular mesh generation tool for Finite Element Methods that solve acoustic/elastic wave equations.

Resources

License

Code of conduct

Stars

Watchers

Forks

Packages

No packages published

Languages

  • Python 67.7%
  • C++ 27.4%
  • TeX 3.7%
  • Other 1.2%