SeismicMesh
2D/3D serial and parallel triangular mesh generation tool for finite element methods.
Install / Use
/learn @krober10nd/SeismicMeshREADME
SeismicMesh: Triangular Mesh generation in Python
SeismicMesh is a Python package for simplex mesh generation in two or three dimensions. As an implementation of DistMesh, it produces high-geometric quality meshes at the expense of speed. For increased efficiency, the core package is written in C++, works in parallel, and uses the Computational Geometry Algorithms Library. SeismicMesh can also produce mesh-density functions from seismological data to be used in the mesh generator.
SeismicMesh is distributed under the GPL3 license and more details can be found in our short paper.
Table of contents
<!--ts-->- Installation
- Contributing
- Codebase
- Citing
- Getting help
- Examples
- Parallelism
- Performance comparison
- Changelog
Installation
For installation, SeismicMesh needs CGAL:
sudo apt install libcgal-dev
After that, SeismicMesh can be installed from the Python Package Index (pypi), so with:
pip install -U SeismicMesh
If you'd like to read and write velocity models from segy/h5 format, you can install like:
pip install -U SeismicMesh[io]
For more detailed information about installation and requirements see:
Install - How to install SeismicMesh.
Contributing
All contributions are welcome!
To contribute to the software:
- Fork the repository.
- Clone the forked repository, add your contributions and push the changes to your fork.
- 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.
Codebase
Here is a visual overview of the repository. An interactive version of this image can be found here: https://octo-repo-visualization.vercel.app/?repo=krober10nd%2FSeismicMesh
Citing
You may use the following BibTeX entry:
@article{Roberts2021,
doi = {10.21105/joss.02687},
url = {https://doi.org/10.21105/joss.02687},
year = {2021},
publisher = {The Open Journal},
volume = {6},
number = {57},
pages = {2687},
author = {Keith J. Roberts and Rafael dos Santos Gioria and William J. Pringle},
title = {SeismicMesh: Triangular meshing for seismology},
journal = {Journal of Open Source Software}
}
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.
BP2004
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 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",
)
Note SeismicMesh can also be used to write velocity models to disk in a hdf5 format using the function write_velocity_model. Following the previous example above with the BP2004 velocity model, we create an hdf5 file with a domain pad of 1000 m.
from SeismicMesh import write_velocity_model
# 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)
write_velocity_model(
fname,
ofname="bp2004_velocity_model", # how the file will be called (with a .hdf5 extension)
bbox=bbox,
domain_pad=500, # the width of the domain pad in meters
pad_style="edge", # how the velocity data will be extended into the layer
units="m-s", # the units that the velocity model is in.
)
EAGE
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/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,
