Modules

Here we document the public API

SeimsicMesh.geometry

Routines to perform geometrical/topological operations and calculate things on meshes.

SeismicMesh.geometry.calc_re_ratios(vertices, entities, dim=2)

Calculate radius edge ratios–mesh quality metric

Parameters:
  • vertices (numpy.ndarray[float x dim]) – point coordinates of the mesh vertices’
  • entities (numpy.ndarray[int x (dim + 1)]) – mesh connectivity table
Returns:

ce_ratios: array of radius-to-edge ratios
Return type:

ce_ratios: numpy.ndarray[float x 1]
SeismicMesh.geometry.do_any_overlap(vertices, entities, dim=2)

Check if any entities connected to boundary of the mesh overlap ignoring self-intersections. This routine checks only the 1-ring around each entity for potential intersections using barycentric coordinates.

Parameters:
  • vertex (numpy.ndarray[float x dim]) – vertex coordinates of mesh
  • entities (numpy.ndarray[int x (dim+1)]) – mesh connectivity
  • dim (int, optional) – dimension of mesh
Returns:

intersections: a list of 2-tuple of entity indices that intersect
Return type:

List[tuple(num_intersections x 2)]
SeismicMesh.geometry.linter(vertices, entities, dim=2, min_qual=0.1)

Remove and check mesh for geometric and toplogical defects.

Parameters:
  • vertex (numpy.ndarray[float x dim]) – vertex coordinates of mesh
  • entities (numpy.ndarray[int x (dim+1)]) – mesh connectivity
  • dim (int, optional) – dimension of mesh
  • min_qual (float, optional) – minimum geometric quality to consider “poor” quality
Return vertices:
 

updated mesh vertices
Return type:

numpy.ndarray[float x dim]
Returns:

entities: updated mesh connectivity
Return type:

numpy.ndarray[int x (dim+1)]
SeismicMesh.geometry.laplacian2(vertices, entities, max_iter=20, tol=0.01)

Move vertices to the average position of their connected neighbors with the goal to hopefully improve geometric entity quality.

Parameters:
  • vertices (numpy.ndarray[float x dim]) – vertex coordinates of mesh
  • entities (numpy.ndarray[int x (dim+1)]) – the mesh connectivity
  • max_iter (int, optional) – maximum number of iterations to perform
  • tol (float, optional) – iterations will cease when movement < tol
Return vertices:
 

updated vertices of mesh
Return type:

numpy.ndarray[float x dim]
Returns:

entities: updated mesh connectivity
Return type:

numpy.ndarray[int x (dim+1)]
SeismicMesh.geometry.vertex_to_entities(vertices, entities, dim=2)

Determine which elements are connected to which vertices.

Parameters:
  • vertices (numpy.ndarray[float x dim]) – point coordinates of mesh vertices
  • entities (numpy.ndarray[int x (dim + 1)]) – mesh connectivity
  • dim (int, optional) – dimension of mesh
Returns:

vtoe: indices of entities connected to each vertex
Return type:

numpy.ndarray[int x 1]
Returns:

vtoe_pointer: indices into vtoe such that vertex v is connected to vtoe[vtoe_pointer[v]:vtoe_pointer[v+1]] entities
Return type:

numpy.ndarray[int x 1]
SeismicMesh.geometry.remove_external_entities(vertices, entities, extent, dim=2)

Remove entities with all dim+1 vertices outside block.

Parameters:
  • vertices (numpy.ndarray[float x dim]) – point coordinates of mesh
  • entities (numpy.ndarray[int x (dim + 1)]) – mesh connectivity
  • extent (numpy.ndarray[tuple(float x (2*dim))]) – coords. of the local block extents
  • dim (int, optional) – dimension of mesh
Returns:

vertices_new: point coordinates of mesh w/ removed entities
Return type:

numpy.ndarray[float x dim]
Returns:

entities_new: mesh connectivity w/ removed entities
Return type:

numpy.ndarray[int x (dim +1)]
Returns:

jx: mapping from old point indexing to new point indexing
Return type:

numpy.ndarray[int x 1]
SeismicMesh.geometry.unique_rows(A, return_index=False, return_inverse=False)

Similar to MATLAB’s unique(A, ‘rows’), this returns B, I, J where B is the unique rows of A and I and J satisfy A = B[J,:] and B = A[I,:]

Parameters:
  • A (numpy.ndarray[int/float x N]) – array of data
  • return_index (boolean, optional) – whether to return the indices of unique data
  • return_inverse (boolean, optional) – whether to return the inverse mapping back to A from B.
Returns:

B: array of data with duplicates removed
Return type:

numpy.ndarray[int/float x N]
Returns:

I: array of indices to unique data B.
Return type:

numpy.ndarray[int x 1]
Returns:

J: array of indices to A from B.
Return type:

numpy.ndarray[int x 1]
SeismicMesh.geometry.fix_mesh(p, t, ptol=2e-13, dim=2, delete_unused=False)
Remove duplicated/unused vertices and entities and
ensure orientation of entities is CCW.
Parameters:
  • p (numpy.ndarray[float x dim]) – point coordinates of mesh
  • t (numpy.ndarray[int x (dim + 1)]) – mesh connectivity
  • ptol (float, optional) – point tolerance to detect duplicates
  • dim (int, optional) – dimension of mesh
  • delete_unused (boolean, optional) – flag to delete disjoint vertices.
Returns:

p: updated point coordinates of mesh
Return type:

numpy.ndarray[float x dim]
Returns:

t: updated mesh connectivity
Return type:

numpy.ndarray[int x (dim+1)]
SeismicMesh.geometry.simp_vol(p, t)

Signed volumes of the simplex elements in the mesh.

Parameters:
  • p (numpy.ndarray[float x dim]) – point coordinates of mesh
  • t (numpy.ndarray[int x (dim + 1)]) – mesh connectivity
Returns:

volume: signed volume of entity/simplex.
Return type:

numpy.ndarray[float x 1]
SeismicMesh.geometry.simp_qual(p, t)

Simplex quality radius-to-edge ratio

Parameters:
  • p (numpy.ndarray[float x dim]) – vertex coordinates of mesh
  • t (numpy.ndarray[int x (dim + 1)]) – mesh connectivity
Returns:

signed mesh quality: signed mesh quality (1.0 is perfect)
Return type:

numpy.ndarray[float x 1]
SeismicMesh.geometry.get_centroids(vertices, entities, dim=2)

Calculate the centroids of all the entities.

Parameters:
  • vertex (numpy.ndarray[float x dim]) – vertex coordinates of mesh
  • entities (numpy.ndarray[int x (dim+1)]) – mesh connectivity
  • dim (int, optional) – dimension of mesh
Returns:

centroids of entities
Return type:

numpy.ndarray[float x dim]
SeismicMesh.geometry.get_edges(entities, dim=2)

Get the undirected edges of mesh in no order (NB: are repeated)

Parameters:
  • entities (numpy.ndarray[int x (dim+1)]) – the mesh connectivity
  • dim (int, optional) – dimension of the mesh
Returns:

edges: the edges that make up the mesh
Return type:

numpy.ndarray[`int`x 2]
SeismicMesh.geometry.get_facets(entities)

Gets the four facets of each tetrahedral.

Parameters:entities (numpy.ndarray[int x (dim+1)]) – the mesh connectivity
Returns:facets: facets of a tetrahedral entity.
Return type:numpy.ndarray[int x 4]
SeismicMesh.geometry.get_boundary_vertices(entities, dim=2)

Get the indices of the mesh representing boundary vertices.

Parameters:
  • entities (numpy.ndarray[`int `x (dim+1)]) – the mesh connectivity
  • dim (int, optional) – dimension of the mesh
Returns:

indices: indices into the vertex array that are on the boundary.
Return type:

numpy.ndarray[float x dim]
SeismicMesh.geometry.get_boundary_entities(vertices, entities, dim=2)

Determine the entities that lie on the boundary of the mesh.

Parameters:
  • vertices (numpy.ndarray[float x dim]) – vertex coordinates of mesh
  • entities (numpy.ndarray[int x (dim+1)]) – the mesh connectivity
  • dim (int, optional) – dimension of the mesh
Returns:

bele: indices of entities on the boundary of the mesh.
Return type:

numpy.ndarray[int x 1]
SeismicMesh.geometry.delete_boundary_entities(vertices, entities, dim=2, min_qual=0.1)

Delete boundary entities with poor geometric quality (i.e., < min. quality)

Parameters:
  • vertices (numpy.ndarray[float x dim]) – vertex coordinates of mesh
  • entities (numpy.ndarray[int x (dim+1)]) – the mesh connectivity
  • dim (int, optional) – dimension of the mesh
  • min_qual (float, optional) – minimum geometric quality to consider “poor” quality
Returns:

vertices: updated vertex array of mesh
Return type:

numpy.ndarray[int x dim]
Returns:

entities: update mesh connectivity
Return type:

numpy.ndarray[int x (dim+1)]
SeismicMesh.geometry.get_boundary_edges(entities, dim=2)

Get the boundary edges of the mesh. Boundary edges only appear (dim-1) times

Parameters:
  • entities (numpy.ndarray[int x (dim+1)]) – the mesh connectivity
  • dim (int, optional) – dimension of the mesh
Returns:

boundary_edges: the edges that make up the boundary of the mesh
Return type:

numpy.ndarray[int x 2]
SeismicMesh.geometry.get_boundary_facets(entities)

Get the facets that represent the boundary of a 3D mesh.

Parameters:entities (numpy.ndarray[int x (dim+1)]) – the mesh connectivity
Returns:boundary_facets: facets on the boundary of a 3D mesh.
Return type:numpy.ndarray[int x 4]
SeismicMesh.geometry.get_winded_boundary_edges(entities)

Order the boundary edges of the mesh in a winding fashiono

Parameters:entities (numpy.ndarray[int x (dim+1)]) – the mesh connectivity
Returns:boundary_edges: the edges that make up the boundary of the mesh in a winding order
Return type:numpy.ndarray[int x 2]
SeismicMesh.geometry.vertex_in_entity3(vertex, entity)

Does the 3D vertex lie in the entity defined by vertices (x1,y1,z1,x2,y2,z2,x3,y3,z3)?

Parameters:
  • vertex (numpy.ndarray[float x dim]) – vertex coordinates of mesh
  • entity (numpy.ndarray[int x (dim+1)]) – connectivity of an entity
Returns:

vertex_in_entity3: logical flag indicating if it is or isn’t.
Return type:

boolean

SeimsicMesh.sizing

Function to build a \(f(h)\) mesh sizing function from a seismic velocity model. Assumes the domain can be represented by a rectangle (2D) or cube (3D) and thus builds a \(f(d)\) accordingly.

SeismicMesh.sizing.get_sizing_function_from_segy(filename, bbox, comm=None, **kwargs)

Build a mesh size function from a seismic velocity model.

Parameters:
  • filename (string) – The name of a SEG-y or binary file containing a seismic velocity model
  • bbox (tuple with size (2*dim). For example, in 2D (zmin, zmax, xmin, xmax)) – Bounding box containing domain extents of the velocity model contained in filename.
  • **kwargs – See below
Keyword Arguments:
 
  • hmin (float) –
    Minimum edge length in the domain (default==150 m)
  • hmax (float) –
    Maximum edge length in the domain (default==10,000 m)
  • wl (int) –
    Number of vertices per wavelength for a given 𝑓𝑚𝑎𝑥 (default==0 vertices)
  • freq (float) –
    𝑓𝑚𝑎𝑥 in hertz for which to estimate wl (default==2 Hertz)
  • grad (float) –
    Resolution in meters nearby sharp gradients in velociy (default==0 m)
  • grade (float) –
    Maximum allowable variation in mesh size in decimal percent (default==0.0)
  • space_order (int) –
    Simulation will be attempted with a mesh using the polynomial order space_order of the basis functions (default==1)
  • dt (float) –
    Theoretical maximum stable timestep in seconds given Courant number Cr (default==0.0 s)
  • cr_max (float) –
    The mesh simulated with this dt has this maximum Courant number (default==1.0)
  • pad_style (string) –
    The method (edge, linear_ramp, constant) to pad velocity in the domain pad region (default==None)
  • domain_pad (float) –
    The width of the domain pad in -z, +x, -x, +y, -y directions (default==0.0 m).
  • units (string) –
    The units of the seismic velocity model (default=’m-s’)
  • nz (int) –
    REQUIRED FOR BINARY VELOCITY MODEL. The number of grid points in the z-direction in the velocity model.
  • ny (int) –
    REQUIRED FOR BINARY VELOCITY MODEL. The number of grid points in the y-direction in the velocity model.
  • nx (int) –
    REQUIRED FOR BINARY VELOCITY MODEL. The number of grid points in the x-direction in the velocity model.
  • byte_order (string) –
    REQUIRED FOR BINARY VELOCITY MODEL. The order of bytes in a 3D sesimic velocity model (little or big).
Returns:

a SizeFunction object with a obj.bbox field and an obj.eval method.
Return type:

a SizeFunction object
SeismicMesh.sizing.write_velocity_model(filename, ofname=None, comm=None, **kwargs)

Reads and then writes a velocity model as a hdf5 file

Parameters:
  • filename (string) – filename of a velocity model (either .segy or .bin)
  • ofname (string, optional) – filename of output hdf5 file
  • comm (MPI4py communicator, optional) – MPI communicator
  • **kwargs – See below
Keyword Arguments:
 
  • nz (int) –
    REQUIRED FOR BINARY VELOCITY MODEL. The number of grid points in the z-direction in the velocity model.
  • ny (int) –
    REQUIRED FOR BINARY VELOCITY MODEL. The number of grid points in the y-direction in the velocity model.
  • nx (int) –
    REQUIRED FOR BINARY VELOCITY MODEL. The number of grid points in the x-direction in the velocity model.
  • byte_order (string) –
    REQUIRED FOR BINARY VELOCITY MODEL. The order of bytes in a 3D sesimic velocity model (little or big).
  • bbox (tuple) –
    Coordinates of the velocity model’s domain extents. Only required if padding the domain.
  • domain_pad (float) –
    Width of the domain pad in meters.
  • pad_style (string) –
    Type of padding.
SeismicMesh.sizing.plot_sizing_function(cell_size, stride=1, comm=None)

Plot the mesh size function in 2D

Parameters:
  • cell_size (a callable function object) – a callable function that takes a point and gives a size
  • stride (int, optional) – skip stride points to save on memory when plotting
  • comm (MPI4py communicator, optional) – MPI communicator

SeimsicMesh.generation

Functions to build and improve a simplical mesh that conforms to the signed distance function \(f(d)\) and \(f(h)\).

SeismicMesh.generation.generate_mesh(domain, edge_length, comm=None, **kwargs)

Generate a 2D/3D mesh using callbacks to a sizing function edge_length and signed distance function domain

Parameters:
  • domain (A geometry.Rectangle/Cube/Disk object or a function object.) – A function that takes a point and returns the signed nearest distance to the domain boundary Ω.
  • edge_length (A SizeFunction object, a function object, or a scalar value.) – Edge lengths in the domain (e.g., triangular edge lengths assuming all triangles are equilateral).
  • comm (MPI4py communicator object, optional) – MPI4py communicator
  • **kwargs – See below
Keyword Arguments:
 
  • h0 (float) –
    The minimum edge length in the domain. REQUIRED IF USING A VARIABLE RESOLUTION EDGE LENGTH
  • bbox (tuple) –
    Bounding box containing domain extents. REQUIRED IF NOT USING edge_length
  • verbose (int) –
    Output to the screen verbose (default==1). If verbose`==1 only start and end messages are written, `verbose`==0, no messages are written besides errors, `verbose > 1 all messages are written.
  • max_iter (float) –
    Maximum number of meshing iterations. (default==50)
  • seed (float or int) –
    Psuedo-random seed to initialize meshing points. (default==0)
  • perform_checks (boolean) –
    Whether or not to perform mesh linting/mesh cleanup. (default==False)
  • pfix (array-like) –
    An array of points to constrain in the mesh. (default==None)
  • axis (int) –
    The axis to decompose the mesh (1,2, or 3). (default==1)
  • delta_t (float) –
    Psuedo-timestep to control movement of points (default=0.10)
Returns:

points: vertex coordinates of mesh
Return type:

points: (numpy.ndarray[float x dim])
Returns:

t: mesh connectivity table.
Return type:

t: numpy.ndarray[int x (dim + 1)]
SeismicMesh.generation.sliver_removal(points, domain, edge_length, comm=None, **kwargs)

Improve an existing 3D mesh by removing degenerate cells. commonly referred to as slivers.

Parameters:
  • points (np.ndarray) – An array of points that describe the vertices of an existing (higher-quality) mesh.
  • domain (A geometry.Rectangle/Cube/Disk object or a function object.) – A function that takes a point and returns the signed nearest distance to the domain boundary Ω
  • edge_length (A SizeFunction object or a function object.) – A function that can evalulate a point and return an edge length (e.g. length of the triangular edge)
  • comm (MPI4py communicator object, optional) – MPI4py communicator
  • **kwargs – See below
Keyword Arguments:
 
  • h0 (float) –
    The minimum edge length in the domain. REQUIRED IF USING A VARIABLE RESOLUTION EDGE LENGTH.
  • verbose (int) –
    Output to the screen verbose (default==1). If verbose`==1 only start and end messages are written, `verbose`==0, no messages are written besides errors, `verbose > 1 all messages are written.
  • max_iter (float) –
    Maximum number of meshing iterations. (default==50)
  • perform_checks (boolean) –
    Whether or not to perform mesh linting/mesh cleanup. (default==False)
  • pfix (array-like) –
    An array of points to constrain in the mesh. (default==None)
  • axis (int) –
    The axis to decompose the mesh (1,2, or 3). (default==1)
  • min_dh_angle_bound (float) –
    The minimum allowable dihedral angle bound. (default==10 degrees)
  • max_dh_angle_bound (float) –
    The maximum allowable dihedral angle bound. (default==170 degrees)