salvus.mesh.algorithms.unstructured_mesh.utils

def add_element_nodal_fields(
    meshes: typing.List[UnstructuredMesh],
) -> 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: np.ndarray, points: np.ndarray
) -> np.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 create_rotation_matrices(
    vectors_from: npt.NDArray, vectors_to: npt.NDArray
) -> npt.NDArray: ...

Create the rotation matrices for rotating vectors_from onto vectors_to. Will always create 3D rotation matrices. If either vectors_from or vectors_to is passed in shape (…, 2), zeros are appended to create 3D vectors in the last index. Follows general broadcasting rules, output shape will be np.broadcast(vectors_A, vectors_B).shape[:-1] + (3, 3). If drop_dimension is True, and one of the original input shapes was 2 dimensional, the returned rotation matrices will be 2D. This might be incorrect, if rotation in 3D is required to fully make the rotation.

def disconnect_along_side_set(
    mesh: UnstructuredMesh,
    side_set: str,
    disconnection_decider: typing.Callable[
        [UnstructuredMesh, npt.NDArray, npt.NDArray, npt.NDArray, npt.NDArray],
        typing.Iterable[bool],
    ],
) -> UnstructuredMesh: ...

Disconnect mesh along internal side set, turning it into two free surfaces.

The disconnection is made by using a callback able to tell for any element if it is on side A or side B of the disconnection, in relation to a point.

def edge_lengths(mesh: UnstructuredMesh) -> np.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: UnstructuredMesh) -> np.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: UnstructuredMesh) -> np.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: UnstructuredMesh,
    max_radius: float = 6371000.0,
    exclude_filter: typing.Optional[typing.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: UnstructuredMesh,
    ds: xr.Dataset,
    pars: typing.Union[str, typing.List[str]],
    max_tree_doublings: int = 4,
    verbose: bool = False,
) -> xr.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: _ElementCollectionProtocol,
    points: np.ndarray,
    max_tree_doublings: typing.Union[str, int] = "auto",
    allow_points_outside_mesh: bool = True,
    verbose: bool = False,
    element_restrict: typing.Optional[np.ndarray] = None,
    point_restrict: typing.Optional[np.ndarray] = None,
) -> typing.Tuple[np.ndarray, np.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: UnstructuredMesh,
    m_fine: UnstructuredMesh,
    number_of_neighbours: int = 8,
    verify: bool = True,
) -> np.ndarray: ...

Find the hierarchical map between two meshes.

def get_internal_side_set_facets(
    mesh: UnstructuredMesh, side_set: str
) -> typing.Tuple[npt.NDArray, npt.NDArray]: ...

Get the internal facets of a sideset, i.e. those that have two surrounding elements, and the respective elements.

def get_interpolation_coefficients(
    mesh: _ElementCollectionProtocol, points: np.ndarray
) -> npt.NDArray: ...

Get the interpolation coefficients for a series of reference coordinates.

def get_side_set_elevations(
    mesh: UnstructuredMesh, side_set: str, points: npt.NDArray
) -> npt.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.

For Cartesian and 2-D spherical domains 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. For 3-D spherical domains the result is an approximation built by projecting the points to a set of triangles that span each element.

def get_surface_unit_vectors(
    points: npt.NDArray, mesh: UnstructuredMesh, side_set: str
) -> npt.NDArray: ...

Get the surface unit vectors of a mesh relative to a side set, for specific points.

def get_topological_facets(mesh: UnstructuredMesh) -> npt.NDArray: ...

Get an array of node ids corresponding to each facet in the mesh.

The facets returned here are “topological” in that they only consider the mesh’s first-order connectivity. Returns an array of shape (n_elem, n_facet_per_elem, n_nodes_per_facet). Facet ordering is as followings:

2D Quad: bottom, top, right, left. 3D Hex: Left, right, bottom, top, back, front.

def interpolate_from_element_nodes(
    ref: npt.NDArray, values: npt.NDArray, n_dim: int, num_threads: int = 16
) -> typing.Tuple[npt.NDArray, npt.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: npt.NDArray,
    points: npt.NDArray,
    max_iter: int = 100,
    rtol: float = 1e-08,
    atol: float = 1e-11,
    atol_ref: float = 0.001,
    num_threads: int = 16,
) -> typing.Tuple[npt.NDArray[np.float_], npt.NDArray[np.int_]]: ...

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

def name_free_side_set(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: UnstructuredMesh,
    block: typing.Union[int, typing.List[int]],
    top_side_set: str,
    bot_side_set: str,
) -> typing.Tuple[npt.NDArray, npt.NDArray, npt.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: typing.Union[str, pathlib.Path],
    fields: typing.Optional[typing.Sequence[str]] = None,
) -> typing.Dict[str, np.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.

def remove_from_side_set(
    mesh: UnstructuredMesh,
    side_set: str,
    indices_to_remove: typing.Optional[npt.NDArray] = None,
    indices_to_retain: typing.Optional[npt.NDArray] = None,
) -> None: ...

Remove or retain facets from side set by index, operating on the passed mesh. Facets are indexed according to their order of appearance in mesh.side_sets[side_set].

def retain_unique_facets(
    mesh: UnstructuredMesh, side_set: str, verbosity: int = 0
) -> None: ...

Retain only facets in side set that don’t appear in other sidesets, deleting the side set if it becomes empty.

def uniquefy_side_sets(
    mesh: UnstructuredMesh,
    order: typing.Optional[typing.Iterable[str]] = None,
    verbosity: int = 0,
) -> None: ...

Retain only unique facets in side sets, deleting empty side sets.

The order influences what elements are retained. The algorithm will remove facets from sidesets in this order, i.e. the first side set has all facets subtracted that occur in other side sets. Side sets appearing earlier will have more facets removed. Empty side sets will not be retained.