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:
  1. slowness constant within the voxels of the mesh (the default)

  2. 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:

Mesh2d

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 of np.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:

Mesh3d

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 of csr_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 of np.ndarray) – Coordinates of segments forming raypaths (if return_rays is True)

  • L (list of csr_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