Module tmesh
Raytracing on unstructured triangular and tetrahedral meshes
This module contains two classes to perform traveltime computation and raytracing on unstructured meshes:
Mesh2d for 2D media
Mesh3d for 3D media
- Three algorithms are implemented
the Shortest-Path Method
the Fast-Sweeping Method
the Dynamic Shortest-Path Method
- Slowness model can be defined in two ways:
slowness constant within the voxels of the mesh (the default)
slowness defined at nodes of the mesh
This code is part of ttcr ( https://github.com/groupeLIAMG/ttcr )
- class ttcrpy.tmesh.Mesh2d
class to perform raytracing with triangular meshes
- nparams
total number of parameters for grid
- Type:
int
- n_threads
number of threads for raytracing
- Type:
int
- Constructor
- Mesh2d(nodes, triangles, n_threads=1, cell_slowness=1, method='FSM', aniso='iso', eps=1e-15, maxit=20, process_obtuse=1, n_secondary=5, n_tertiary=2, radius_factor_tertiary=2, tt_from_rp=0) Mesh2d
- Parameters:
nodes (np.ndarray, shape (nnodes, 2)) – node coordinates
triangles (np.ndarray of int, shape (ntriangles, 3)) – indices of nodes forming the triangles
n_threads (int) – number of threads for raytracing (default is 1)
cell_slowness (bool) – slowness defined for cells (True) or nodes (False) (default is 1)
method (string) –
- raytracing method (default is FSM)
’FSM’ : fast marching method
’SPM’ : shortest path method
’DSPM’ : dynamic shortest path
aniso (string) –
- type of anisotropy (implemented only for the SPM method)
’iso’ : isotropic medium
’elliptical’ : elliptical anisotropy
’tilted_elliptical’ : tilted elliptical anisotropy
’weakly_anelliptical’ : Weakly-Anelliptical formulation of B. Rommel
eps (double) – convergence criterion (FSM) (default is 1e-15)
maxit (int) – max number of sweeping iterations (FSM) (default is 20)
process_obtuse (bool) – use method of Qian et al (2007) to improve accuracy for triangles with obtuse angle (default is True)
n_secondary (int) – number of secondary nodes (SPM) (default is 5)
n_tertiary (int) – number of tertiary nodes (DSPM) (default is 2)
radius_factor_tertiary (double) – multiplication factor used to compute radius of sphere around source that includes tertiary nodes (DSPM). The radius is the average edge length multiplied by this factor (default is 2)
tt_from_rp (bool) – compute traveltimes using raypaths (default is False)
Notes
For raytracing in anisotropic media, the convention for inputting slowness depends on the model. For elliptical anisotropy, the method set_slowness is used to input horizontal slowness, while for weakly anelliptical anisotropy, the method is used to input vertical slowness.
- static builder(filename, n_threads, cell_slowness, method, eps, maxit, process_obtuse, n_secondary, n_tertiary, radius_factor_tertiary, tt_from_rp)
Build instance of Mesh2d from VTK file
- Parameters:
filename (str) – Name of file holding a vtkUnstructuredGrid. The grid must have point or cell attribute named either ‘Slowness’, ‘slowness’, ‘Velocity’, ‘velocity’, or ‘P-wave velocity’. All cells must be of type vtkTriangle
Constructor (Other parameters are defined in) –
- Returns:
mesh – mesh instance
- Return type:
- get_grid_traveltimes(thread_no=0)
Obtain traveltimes computed at primary grid nodes
- Parameters:
thread_no (int) – thread used to computed traveltimes (default is 0)
- Returns:
tt – traveltimes
- Return type:
np ndarray, shape (nnodes,)
- get_number_of_cells()
- Returns:
number of cells in grid
- Return type:
int
- get_number_of_nodes()
- Returns:
number of nodes in grid
- Return type:
int
- n_threads
number of threads for raytracing
- Type:
int
- nparams
total number of parameters for mesh
- Type:
int
- raytrace(source, rcv, slowness=None, thread_no=None, aggregate_src=False, compute_L=False, return_rays=False) tt, rays
Perform raytracing
- Parameters:
source (2D np.ndarray with 2 or 3 columns) – see notes below
rcv (2D np.ndarray with 2 columns) – Columns correspond to x and z coordinates
slowness (np ndarray, (None by default)) – slowness at grid nodes or cells (depending on cell_slowness) if None, slowness must have been assigned previously
thread_no (int (None by default)) – Perform calculations in thread number “thread_no” if None, attempt to run in parallel if warranted by number of sources and value of n_threads in constructor
aggregate_src (bool (False by default)) – if True, all source coordinates belong to a single event
compute_L (bool (False by default)) – Compute matrices of partial derivative of travel time w/r to slowness
return_rays (bool (False by default)) – Return raypaths
- Returns:
tt (np.ndarray) – travel times for the appropriate source-rcv (see Notes below)
rays (
list
ofnp.ndarray
) – Coordinates of segments forming raypaths (if return_rays is True)
Notes
- If source has 2 columns:
Columns correspond to x and z coordinates
Origin time (t0) is 0 for all points
- If source has 3 columns:
1st column corresponds to origin times
2nd & 3rd columns correspond to x and z coordinates
source and rcv can contain the same number of rows, each row corresponding to a source-receiver pair, or the number of rows may differ if aggregate_src is True or if all rows in source are identical.
- set_s2(s2)
Assign energy-velocity parameter s2 to grid
- Parameters:
s2 (np ndarray, shape (nparams, )) –
- set_s4(s4)
Assign energy-velocity parameter s4 to grid
- Parameters:
s4 (np ndarray, shape (nparams, )) –
- set_slowness(slowness)
Assign slowness to grid
- Parameters:
slowness (np ndarray, shape (nparams, )) –
- set_tilt_angle(theta)
Assign tilted elliptical anisotropy angle to grid
- Parameters:
theta (np ndarray, shape (nparams, )) –
- set_traveltime_from_raypath(ttrp)
Set option to compute traveltime using raypath
- Parameters:
ttrp (bool) – option value
- set_use_thread_pool(use_thread_pool)
Set option to use thread pool instead of parallel loop
- Parameters:
use_thread_pool (bool) – option value
- set_velocity(velocity)
Assign velocity to grid
- Parameters:
velocity (np ndarray, shape (nparams, )) –
- set_xi(xi)
Assign elliptical anisotropy ratio to grid
- Parameters:
xi (np ndarray, shape (nparams, )) –
- to_vtk(fields, filename)
Save mesh variables and/or raypaths to VTK format
- Parameters:
fields (dict) – dict of variables to save to file. Variables should be np.ndarray of size equal to either the number of nodes of the number of cells of the mesh, or a list of raypath coordinates.
filename (str) – Name of file without extension for saving (extension vtu will be added). Raypaths are saved in separate files, and filename will be appended by the dict key and have a vtp extension.
Notes
VTK files can be visualized with Paraview (https://www.paraview.org)
- class ttcrpy.tmesh.Mesh3d
class to perform raytracing with tetrahedral meshes
- nparams
total number of parameters for grid
- Type:
int
- n_threads
number of threads for raytracing
- Type:
int
- Constructor
- Mesh3d(nodes, tetra, n_threads, cell_slowness, method, gradient_method, tt_from_rp, process_vel, eps, maxit, min_dist, n_secondary, n_tertiary, radius_factor_tertiary, translate_grid=False) Mesh3d
- Parameters:
nodes (np.ndarray, shape (nnodes, 3)) – node coordinates
tetra (np.ndarray of int, shape (ntetra, 4)) – indices of nodes forming the tetrahedra
n_threads (int) – number of threads for raytracing (default is 1)
cell_slowness (bool) – slowness defined for cells (True) or nodes (False) (default is 1)
method (string) –
- raytracing method (default is FSM)
’FSM’ : fast marching method
’SPM’ : shortest path method
’DSPM’ : dynamic shortest path
gradient_method (int) –
- method to compute traveltime gradient (default is 1)
0 : least-squares first-order
1 : least-squares second-order
2 : Averaging-Based method
tt_from_rp (bool) – compute traveltimes from raypaths (FSM or DSPM only) (default is 1)
process_vel (bool) – process velocity instead of slowness at nodes when interpolating and computing matrix of partial derivative of traveltime w/r to model parameters (interpolation: for cell_slowness == False or FSM) (defauls is False)
eps (double) – convergence criterion (FSM) (default is 1e-15)
maxit (int) – max number of sweeping iterations (FSM) (default is 20)
min_dist (double) – tolerance for backward raytracing (default is 1e-5)
n_secondary (int) – number of secondary nodes (SPM & DSPM) (default is 2)
n_tertiary (int) – number of tertiary nodes (DSPM) (default is 2)
radius_factor_tertiary (double) – multiplication factor used to compute radius of sphere around source that includes tertiary nodes (DSPM). The radius is the average edge length multiplied by this factor (default is 3)
translate_grid (bool) – Translate the grid such that origin is (0, 0, 0) to perform computations, which may increase accuracy when large values, e.g. UTM coordinates, are used. When raytracing, src and rcv should be given in the original system, and output raypath coordinates are also given in the original system (default if False)
- static builder(filename, n_threads, cell_slowness, method, gradient_method, tt_from_rp, process_vel, eps, maxit, min_dist, n_secondary, n_tertiary, radius_factor_tertiary, translate_grid=0)
Build instance of Mesh3d from VTK file
- Parameters:
filename (str) – Name of file holding a vtkUnstructuredGrid. The grid must have point or cell attribute named either ‘Slowness’, ‘slowness’, ‘Velocity’, ‘velocity’, or ‘P-wave velocity’. All cells must be of type vtkTetra
Constructor (Other parameters are defined in) –
- Returns:
mesh – mesh instance
- Return type:
- compute_D(coord)
Return matrix of interpolation weights for velocity data points constraint
- Parameters:
coord (np.ndarray, shape (npts, 3)) – coordinates of data points
- Returns:
D – Matrix of interpolation weights
- Return type:
scipy csr_matrix, shape (npts, nparams)
- compute_K(order=2, taylor_order=2, weighting=True, squared=True, s0inside=False, additional_points=0)
Compute smoothing matrices (spatial derivative)
- Parameters:
order (int) – order of derivative (1 or 2, 2 by default)
taylor_order (int) – order of taylors series expansion (1 or 2, 2 by default)
weighting (bool) – apply inverse distance weighting (True by default)
squared (bool) – Second derivative evaluated by taking the square of first derivative. Applied only if order == 2 (True by default)
s0inside (bool) – (experimental) ignore slowness value at local node (value is a filtered estimate) (False by default)
additional_points (int) – use additional points to compute derivatives (minimum sometimes yield noisy results when rays are close to domain limits) (0 by default)
- Returns:
Kx, Ky, Kz – matrices for derivatives along x, y, & z
- Return type:
tuple
ofcsr_matrix
- data_kernel_straight_rays(Tx, Rx) L
Raytracing with straight rays in 3D
- Parameters:
Tx (np.ndarray) –
- source coordinates, nTx by 3
1st column contains X coordinates,
2nd contains Y coordinates
3rd contains Z coordinates
Rx (np.ndarray) –
- receiver coordinates, nTx by 3
1st column contains X coordinates,
2nd contains Y coordinates
3rd contains Z coordinates
- Returns:
L – data kernel matrix (tt = L*slowness)
- Return type:
scipy csr_matrix
Note
Tx and Rx should contain the same number of rows, each row corresponding to a source-receiver pair
- get_grid_traveltimes(thread_no=0)
Obtain traveltimes computed at primary grid nodes
- Parameters:
thread_no (int) – thread used to computed traveltimes (default is 0)
- Returns:
tt – traveltimes
- Return type:
np ndarray, shape (nnodes,)
- get_number_of_cells()
- Returns:
number of cells in grid
- Return type:
int
- get_number_of_nodes()
- Returns:
number of nodes in grid
- Return type:
int
- get_s0(hypo, slowness=None)
Return slowness at source points
- Parameters:
hypo (np.ndarray with 5 columns) –
- hypo holds source information, i.e.
1st column is event ID number
2nd column is origin time
3rd column is source easting
4th column is source northing
5th column is source elevation
slowness (np ndarray, shape (nparams, ) (optional)) – slowness at grid nodes or cells (depending on cell_slowness)
- Returns:
s0 – slowness at source points
- Return type:
np.ndarray
- is_outside(pts)
Check if points are outside grid
- Parameters:
pts (np ndarray, shape(npts, 3)) – coordinates of points to check
- Returns:
True if at least one point outside grid
- Return type:
bool
- n_threads
number of threads for raytracing
- Type:
int
- nparams
total number of parameters for mesh
- Type:
int
- raytrace(source, rcv, slowness=None, thread_no=None, aggregate_src=False, compute_L=False, return_rays=False) tt, rays, L
Perform raytracing
- Parameters:
source (2D np.ndarray with 3, 4 or 5 columns) – see notes below
rcv (2D np.ndarray with 3 columns) – Columns correspond to x, y and z coordinates
slowness (np ndarray, (None by default)) – slowness at grid nodes or cells (depending on cell_slowness) if None, slowness must have been assigned previously
thread_no (int (None by default)) – Perform calculations in thread number “thread_no” if None, attempt to run in parallel if warranted by number of sources and value of n_threads in constructor
aggregate_src (bool (False by default)) – if True, all source coordinates belong to a single event
compute_L (bool (False by default)) – Compute matrices of partial derivative of travel time w/r to slowness (or velocity if process_vel == True in constructor)
return_rays (bool (False by default)) – Return raypaths
- Returns:
tt (np.ndarray) – travel times for the appropriate source-rcv (see Notes below)
rays (
list
ofnp.ndarray
) – Coordinates of segments forming raypaths (if return_rays is True)L (
list
ofcsr_matrix
or scipy csr_matrix) – Matrix of partial derivative of travel time w/r to slowness. if input argument source has 5 columns or if slowness is defined at nodes, L is a list of matrices and the number of matrices is equal to the number of sources otherwise, L is a single csr_matrix
Notes
- If source has 3 columns:
Columns correspond to x, y and z coordinates
Origin time (t0) is 0 for all points
- If source has 4 columns:
1st column corresponds to origin times
2nd, 3rd & 4th columns correspond to x, y and z coordinates
- If source has 5 columns:
1st column corresponds to event ID
2nd column corresponds to origin times
3rd, 4th & 5th columns correspond to x, y and z coordinates
For the latter case (5 columns), source and rcv should contain the same number of rows, each row corresponding to a source-receiver pair. For the 2 other cases, source and rcv can contain the same number of rows, each row corresponding to a source-receiver pair, or the number of rows may differ if aggregate_src is True or if all rows in source are identical.
- set_slowness(slowness)
Assign slowness to grid
- Parameters:
slowness (np ndarray, shape (nparams, )) –
- set_traveltime_from_raypath(ttrp)
Set option to compute traveltime using raypath
- Parameters:
ttrp (bool) – option value
- set_use_thread_pool(use_thread_pool)
Set option to use thread pool instead of parallel loop
- Parameters:
use_thread_pool (bool) – option value
- set_velocity(velocity)
Assign velocity to grid
- Parameters:
velocity (np ndarray, shape (nparams, )) –
- to_vtk(fields, filename)
Save mesh variables and/or raypaths to VTK format
- Parameters:
fields (dict) – dict of variables to save to file. Variables should be np.ndarray of size equal to either the number of nodes of the number of cells of the mesh, or a list of raypath coordinates.
filename (str) – Name of file without extension for saving (extension vtu will be added). Raypaths are saved in separate files, and filename will be appended by the dict key and have a vtp extension.
Notes
VTK files can be visualized with Paraview (https://www.paraview.org)
- ttcrpy.tmesh.set_verbose(v)
Set verbosity level for C++ code
- Parameters:
v (int) – verbosity level