Version:
This tutorial is presented as Python code running inside a Jupyter Notebook, the recommended way to use Salvus. To run it yourself you can copy/type each individual cell or directly download the full notebook, including all required files.

Basic Meshing API

The most basic mesh to run with Salvus

This exercise shows the most basic mesh with all necessary commands to put material properties and boundaries ready to run with Salvus.
Copy
# initialize notebook
%matplotlib inline
%config Completer.use_jedi = False
import matplotlib.pyplot as plt
import numpy as np

plt.rcParams["figure.figsize"] = (10, 8)

2D

Set basic input parameters for a rectangular domain with homogeneous material properties.
vp = 1700  # p-wave velocity in m/s
vs = 1000  # s-wave velocity in m/s
rho = 2700  # density in kg / m^3

max_x = 10000  # domain size in x direction
max_y = 5000  # domain size in y direction

fmax = 1.0  # maximum frequency in Hz
elements_per_wavelength = 2.0  # resolution criterion
Compute element size and number of elements needed according to the resolution criterion.
hmax = vs / fmax / elements_per_wavelength
nelem_x = int(max_x / hmax)
nelem_y = int(max_y / hmax)
In the mesher, we build unstructured meshes by combining multiple structured grids. For the simplest case here, we build a single structured grid and plot it.
from salvus.mesh.structured_grid_2D import StructuredGrid2D

sg = StructuredGrid2D.rectangle(
    nelem_x=nelem_x, nelem_y=nelem_y, max_x=max_x, max_y=max_y
)
sg.plot(show=True)
Then, we convert the structured grid to an unstructured mesh (trivial in this special case, but uses the same general routines).
m = sg.get_unstructured_mesh()
Attach the material properties, constant here, but using the general array storage on element nodes.
m.attach_field("VP", vp * np.ones([m.nelem, m.nodes_per_element]))
m.attach_field("VS", vs * np.ones([m.nelem, m.nodes_per_element]))
m.attach_field("RHO", rho * np.ones([m.nelem, m.nodes_per_element]))

# this is necessary to tell salvus that this is elastic material everywhere.
# Setting it to 1 would make it fluid.
m.attach_field("fluid", np.zeros(m.nelem))
Find outer Boundaries of the domain assuming the domain is rectangular.
m.find_side_sets(mode="cartesian")
Open a preview inline, showing all parameters stored on the mesh.
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7f21bb732f90>
The interfaces are very much the same as in 2D, so less comments.
vp = 1700  # p-wave velocity in m/s
vs = 1000  # s-wave velocity in m/s
rho = 2700  # density in kg / m^3

max_x = 10000  # domain size in x direction
max_y = 5000  # domain size in y direction
max_z = 3000  # domain size in y direction

fmax = 1.0  # maximum frequency in Hz
elements_per_wavelength = 2.0  # resolution criterion
hmax = vs / fmax / elements_per_wavelength
nelem_x = int(max_x / hmax)
nelem_y = int(max_y / hmax)
nelem_z = int(max_z / hmax)
Build a structured Grid. Don't plot it, because in 3D matplotlib is to slow.
from salvus.mesh.structured_grid_3D import StructuredGrid3D

sg = StructuredGrid3D.cube(
    nelem_x=nelem_x,
    nelem_y=nelem_y,
    nelem_z=nelem_z,
    max_x=max_x,
    max_y=max_y,
    max_z=max_z,
)
m = sg.get_unstructured_mesh()
m.attach_field("VP", vp * np.ones([m.nelem, m.nodes_per_element]))
m.attach_field("VS", vs * np.ones([m.nelem, m.nodes_per_element]))
m.attach_field("RHO", rho * np.ones([m.nelem, m.nodes_per_element]))
m.attach_field("fluid", np.zeros(m.nelem))
m.find_side_sets(mode="cartesian")
m
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7f20912577d0>
The general strategy in the mesher is to compose unstructured meshes from multiple structured grids. For example, a mesh with a doubling can be thought of as two structured rectilinear grids connected by a doubling grid. The doubling is also considered as a structured grid, with some elements masked out and nodes moved. Nodes that appear multiple times are reduced to single node by lexicographic sorting when converting to an unstructured mesh.
from salvus.mesh.skeleton import Skeleton
Start by creating three structured grids
max_x = 1000.0
y = np.array([0.0, 500.0, 750.0, 1000.0])
nelem_x = 5
nelem_y = np.array([2, 2])

sg1 = StructuredGrid2D.rectangle(
    nelem_x, nelem_y[0], min_x=0, max_x=max_x, min_y=y[0], max_y=y[1]
)

sg2 = StructuredGrid2D.cartesian_doubling_layer(
    nelem_x, min_x=0, max_x=max_x, min_y=y[1], max_y=y[2]
)

sg3 = StructuredGrid2D.rectangle(
    nelem_x * 2, nelem_y[1], min_x=0, max_x=max_x, min_y=y