SeismicMesh 1.2.0

2D/3D serial and parallel triangular mesh generation and mesh improvement for seismology

SeismicMesh: Mesh generation for Seismology in Python

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

Quick Installation

For installation, SeismicMesh needs [CGAL](https://www.cgal.org/):

sudo apt install libcgal-dev

After that, SeismicMesh can be [installed from the Python Package Index](https://pypi.org/project/SeismicMesh/), so with:

pip install -U SeismicMesh

For more detailed information about installation and requirements see:

Install - How to install SeismicMesh.

Usage

The user can quickly build quality 3D meshes from seismic velocity models like so:

import numpy as np
import zipfile

from mpi4py import MPI
import meshio
import SeismicMesh

comm = MPI.COMM_WORLD
size = comm.Get_size()
rank = comm.Get_rank()

# The velocity model was downloaded from here:
# https://s3.amazonaws.com/open.source.geoscience/open_data/seg_eage_models_cd/Salt_Model_3D.tar.gz

# Dimensions of model
nx, ny, nz = 676, 676, 210

path = "velocity_models/Salt_Model_3D/3-D_Salt_Model/VEL_GRIDS/"
# Extract Saltf@@ from SALTF.ZIP
zipfile.ZipFile(path + "SALTF.ZIP", "r").extract("Saltf@@", path=path)

# Load data into a numpy array
with open(path + "Saltf@@", "r") as file:
    vp = np.fromfile(file, dtype=np.dtype("float32").newbyteorder(">"))
    vp = vp.reshape(nx, ny, nz, order="F")
    vp = np.flipud(vp.transpose((2, 0, 1)))  # z, x and then y

# Bounding box describing domain extents (corner coordinates)
bbox = (-4200, 0, 0, 13520, 0, 13520)

# Construct mesh sizing object from velocity model
ef = SeismicMesh.MeshSizeFunction(
    bbox=bbox,
    velocity_grid=vp,
    dt=0.001,
    freq=2,
    wl=5,
    grade=0.25,
    hmin=150,
    hmax=5e3,
    domain_ext=250,
    padstyle="linear_ramp",
)

# Build mesh size function (in parallel)
ef = ef.build()

# Construct a mesh generator object
mshgen = SeismicMesh.MeshGenerator(ef)

# Build the mesh
points, cells = mshgen.build(max_iter=75, axis=1)

# Do mesh improvement in serial to bound lower dihedral angle to >= 5 degrees
points, cells = mshgen.build(
    points=points, mesh_improvement=True, max_iter=50, min_dh_bound=5,
)

if rank == 0:
    # Write to disk (see meshio package for more details)
    meshio.write_points_cells(
        "EAGE_Salt.vtk", points / 1000.0, [("tetra", cells)],
    )

Getting Started

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.

Gallery: