salvus.mesh.unstructured_mesh_utils

def add_element_nodal_fields(
    meshes: List[salvus.mesh.unstructured_mesh.UnstructuredMesh],
) -> salvus.mesh.unstructured_mesh.UnstructuredMesh:
    ...

Sum the element nodal fields with the same name across several meshes.

Meshes must have identical points, element nodal fields names, and connectivity. Only fields with the same names across all mesh instances will be summed.

def compute_topological_facet_normals(
    topological_facets: numpy.ndarray, points: numpy.ndarray
) -> numpy.ndarray:
    ...

Compute unit normal vectors for each topological facet on each element.

As topological facets are required for this function, it is not suitable for general normal vector computation (see the salvus.fem.jacobian module for this use case). This function is useful, however, if one wants to cheaply compute the general orientation of a facet.

def edge_lengths(
    mesh: salvus.mesh.unstructured_mesh.UnstructuredMesh,
) -> numpy.ndarray:
    ...

Compute the edge lengths along each edge of each element.

The ordering, in the reference element, is consistent with the DMPLEX ordering

2-D: bottom, right, top, left 3-D: bottom left, bottom back, bottom right, bottom front, top front, top right, top back, top left, front right, front left, back left, back right

def edge_lengths_projected(
    mesh: salvus.mesh.unstructured_mesh.UnstructuredMesh,
) -> numpy.ndarray:
    ...

Compute the edge lengths for each coordinate direction.

The ordering, in the reference element, is consistent with the DMPLEX ordering

2-D: bottom, right, top, left 3-D: bottom left, bottom back, bottom right, bottom front, top front, top right, top back, top left, front right, front left, back left, back right

def edges(
    mesh: salvus.mesh.unstructured_mesh.UnstructuredMesh,
) -> numpy.ndarray:
    ...

Get an ordered list of the first-order edges for each element.

The ordering, in the reference element, is consistent with the DMPLEX ordering

2-D: bottom, right, top, left 3-D: bottom left, bottom back, bottom right, bottom front, top front, top right, top back, top left, front right, front left, back left, back right

def extract_conservative_bm_file(
    mesh: salvus.mesh.unstructured_mesh.UnstructuredMesh,
    max_radius: float = 6371000.0,
    exclude_filter: Optional[Tuple[str, int]] = None,
) -> str:
    ...

Get a BM file containing minimum parameter values and discontinuities.

Useful to extract a bm file the preserves discontinuities, and contains the min parameter value in each layer. Can be used, along with a mesh-to-mesh interpolation routine, as a background model when a) re-meshing for a higher frequency, or b) re-meshing because material velocities have reduced below some threshold.

The successful use of this routines requires that the parameter "z_node_1D" be present as an element nodal field in the mesh, and that it represents the normalized z-coordinate. The normalized value of this parameter will be clipped to a maximum value of 1.0, in case intermediate mesh manipulations resulted in the stretching of the 1-D radial values.

def extract_model_to_regular_grid(
    mesh: salvus.mesh.unstructured_mesh.UnstructuredMesh,
    ds: xarray.core.dataset.Dataset,
    pars: Union[str, List[str]],
    max_tree_doublings: int = 4,
    verbose: bool = False,
) -> xarray.core.dataset.Dataset:
    ...

Return interpolated model values at locations defined by an xarray dataset.

It is often useful to visualize slices of a 3-D model, or to just in general have a regularly-gridded representation of a model for analysis. This function allows one to generate such a representation. As input, it takes a mesh, n xarray dataset, and a list of parameters to extract. The xarray dataset is likely the only parameter which is not self-explanatory. Here one must pass a dataset with one of the following sets of coordinate dimensions:

{"x", "y"},
{"x", "y", "z"},
{"latitude", "longitude", "radius"}, or
{"latitude", "longitude", "depth"}.

An error will be thrown if the dataset’s coordinate dimensions do not match exactly one of the above. Additionally, if the “depth” variant is chosen, a “radius_in_meters” global dataset attribute must also be present. Parameters also must, of course, exist in the mesh.

This runtime of this function will scale with both a) the number of elements in the mesh, and b) the number of point locations in the xarray dataset. Good performance on most machines is expected (for example, extracting 10,000 points from a 100k element mesh should take less than a second on most machines), but in principle the cost of this function can be as large as you like. Passing verbose = True here could give you insight on how long a longer-running call might take, as it will provide a progress bar of the algorithm’s current state.

For some troublesome points, it may be a nontrivial task to find an enclosing element in the mesh. In this case, max_tree_doublings controls the maximum amount of times that a search will be retried for a delinquent point; in each pass the number of elements searched will be doubled. If an element truly is outside of the mesh (as it may well be when interpolating from a spherical domain), extracted values at those points after max_tree_doublings tries will be marked with np.nan — no extrapolation is performed. This convention was chosen to match the standard xarray and CF convention for missing data.

The grids extracted can be visualized in a number of ways, including with xarray directly and with PyGMT. A nice example of how to visualize something with xarray using a map can be found here: http://xarray.pydata.org/en/stable/plotting.html#maps.

def get_enclosing_elements(
    mesh: salvus.mesh.unstructured_mesh_utils._ElementCollectionProtocol,
    points: numpy.ndarray,
    max_tree_doublings: Union[str, int] = "auto",
    allow_points_outside_mesh: bool = True,
    verbose: bool = False,
    element_restrict: Optional[numpy.ndarray] = None,
    point_restrict: Optional[numpy.ndarray] = None,
) -> Tuple[numpy.ndarray, numpy.ndarray]:
    ...

Get the indices of the enclosing element given an array of points.

It is often helpful to know within which element a given spatial point lies. Given an array of points and a mesh, this function will return indices that correspond to the elements containing each point in the input array. The second element of the returned tuple contains the point’s position in reference coordinates with respect to the corresponding element ID. If an enclosing element is not found and no interrupt is passed, the returned element index will be -(closest_centroid_element_id + 1), and the reference coordinates will be outside of the canonical interval (outside [-1, +1]). This allows for extrapolation if required.

The internal algorithm will begin by trying to quickly extract the closest element to a given point. The points that remain unfound after one pass of this are re-injected into the algorithm with the search size being doubled each time. This continues recursively until max_tree_doublings is reached. For well behaved meshes, the enclosing element should just be found in at most a few passes (likely just one). If you are having issues with unclaimed points, you can try increasing the iteration count. Note: be careful to ensure that most of your points are actually in the domain! The algorithm can get slow at higher iteration counts, so it pays to match up your input point extents with the mesh extents as much as possible.

The behavior of the algorithm in the case where a point is outside all elements can be controlled with the allow_points_outside_mesh flag. If this is set to false, then the function will throw if a single point remains unclaimed after max_tree_doublings. This may suggest that either a) the point is indeed outside the mesh (in which case no extrapolation is performed), or b) more iterations are needed.

def get_hierarchical_map(
    m_coarse: salvus.mesh.unstructured_mesh.UnstructuredMesh,
    m_fine: salvus.mesh.unstructured_mesh.UnstructuredMesh,
    nneighbour: int = 8,
    check_m_fine_is_refinement: bool = True,
) -> numpy.ndarray:
    ...

create a hierarchical map from another unstructured mesh to this one, i.e. other is a refinement of self.

def get_interpolation_coefficients(
    mesh: salvus.mesh.unstructured_mesh_utils._ElementCollectionProtocol,
    points: numpy.ndarray,
) -> numpy.ndarray:
    ...

Get the interpolation coefficients for a series of reference coordinates.

def get_side_set_elevations(
    mesh: salvus.mesh.unstructured_mesh.UnstructuredMesh,
    side_set: str,
    points: numpy.ndarray,
) -> numpy.ndarray:
    ...

Get the elevation values on a side set.

This routine:

- Finds the elements that enclose a set of points
- Interpolates the vertical coordinate value from the nodes of element
  to those points.

This is done using the Lagrange basis of the element itself, so the elevation values returned are an accurate representation of the side set’s elevation as discretized by a given mesh.

def interpolate_from_element_nodes(
    ref: numpy.ndarray,
    values: numpy.ndarray,
    n_dim: int,
    num_threads: int = 32,
) -> Tuple[numpy.ndarray, numpy.ndarray]:
    ...

Interpolate values from element nodes to reference coordiantes.

This function also returns polynomials evaluated at each point in reference coordinates as these are useful for further calculations.

def inverse_coordinate_transform(
    control_nodes: numpy.ndarray,
    points: numpy.ndarray,
    max_iter: int = 100,
    rtol: float = 1e-08,
    atol: float = 1e-11,
    atol_ref: float = 0.001,
    num_threads: int = 32,
) -> Tuple[numpy.ndarray, numpy.ndarray]:
    ...

Find the reference coordinates of a point w.r.t. a given element geometry.

def name_free_side_set(
    mesh: salvus.mesh.unstructured_mesh.UnstructuredMesh, name: str
) -> None:
    ...

Assign a side set to all facets on a mesh’s surface that aren’t in one.

def normalize_block_coordinates(
    mesh: salvus.mesh.unstructured_mesh.UnstructuredMesh,
    block: Union[int, List[int]],
    top_side_set: str,
    bot_side_set: str,
) -> numpy.ndarray:
    ...

Normalize vertical coordiantes between two side sets.

It is sometimes desirable to rescale the vertical coordinates of a collection of elements to the range [0, 1] to, for instance, interpolate material parameters defined in terms of thickness or relative coordinates. This function performs that rescaling.

def read_model_from_h5(
    filename: Union[str, pathlib.Path], fields: Optional[Sequence[str]] = None
) -> Dict[str, numpy.ndarray]:
    ...

Read a model from a mesh file and return it as a dictionary. The function supports reading only a subset of fields using the optional argument fields.