Meshes

Generic classes

class stagpyviz.UnstructuredSphere(*args, **kwargs)

Unstructured grid class for spherical meshes. Provides methods for Cartesian - Spherical coordinates transformations, as well as vector and gradient transformations. Inherits from pyvista.UnstructuredGrid, so all methods and properties of the parent class are available.

Attributes:

points_spherical

The points of the mesh in spherical coordinates. Shape is (number_of_points, 3) with columns corresponding to (radius, colatitude, longitude).

Type:: numpy.ndarray

centroids

The centroids of the cells in Cartesian coordinates. Shape is (number_of_cells, 3) with columns corresponding to (x, y, z).

Type:: numpy.ndarray

centroids_spherical

The centroids of the cells in spherical coordinates. Shape is (number_of_cells, 3) with columns corresponding to (radius, colatitude, longitude).

Type:: numpy.ndarray

Methods:

Transform Cartesian coordinates to Spherical coordinates such that:

\[\begin{split}\begin{split} r &= \sqrt{x^2 + y^2 + z^2} \\ \theta &= \arctan\left(\frac{\sqrt{x^2 + y^2}}{z}\right) \\ \phi &= \arctan\left(\frac{y}{x}\right) \end{split}\end{split}\]

Parameters:

x (numpy.ndarray|float) – The x coordinate(s).
y (numpy.ndarray|float) – The y coordinate(s).
z (numpy.ndarray|float) – The z coordinate(s).

Returns:

A tuple containing the spherical coordinates (R, theta, phi). R is the radial distance, theta is the colatitude (angle from the \(z\)-axis), and phi is the longitude (angle from the \(x\)-axis in the \(xy\)-plane).

Return type:

tuple[numpy.ndarray|float, numpy.ndarray|float, numpy.ndarray|float]

cell_data_to_point_data(pass_cell_data: bool = False) → None

Convert all cell data to point data on the mesh.

Parameters:: pass_cell_data (bool) – If True, the cell data will be kept in the mesh after conversion. If False, the cell data will be removed from the mesh after conversion.

create_cell_field(bs: int = 1, dtype: dtype = <class 'numpy.float64'>) → ndarray

Create a cell field initialized with zeros.

Parameters:

bs (int) – The number of components of the field. Default is 1 (scalar field).
dtype (numpy.dtype) – The data type of the field. Default is np.float64.

Returns:

A cell field initialized with zeros. Scalar fields have shape (number_of_cells,), vector fields have shape (number_of_cells, bs).

Return type:

numpy.ndarray

create_point_field(bs: int = 1, dtype: dtype = <class 'numpy.float64'>) → ndarray

Create a point field initialized with zeros.

Parameters:

bs (int) – The number of components of the field. Default is 1 (scalar field).
dtype (numpy.dtype) – The data type of the field. Default is np.float64.

Returns:

A point field initialized with zeros. Scalar fields have shape (number_of_points,), vector fields have shape (number_of_points, bs).

Return type:

numpy.ndarray

is_cell_field(field: ndarray) → bool

Check if the input field is a cell field.

Parameters:: field (numpy.ndarray) – The field to check.
Returns:: True if the field is a cell field, False otherwise.
Return type:: bool

is_point_field(field: ndarray) → bool

Check if the input field is a point field.

Parameters:: field (numpy.ndarray) – The field to check.
Returns:: True if the field is a point field, False otherwise.
Return type:: bool

point_data_to_cell_data(pass_point_data: bool = False) → None

Convert all point data to cell data on the mesh.

Parameters:: pass_point_data (bool) – If True, the point data will be kept in the mesh after conversion. If False, the point data will be removed from the mesh after conversion.

point_field_to_cell_field(field: ndarray) → ndarray

Transform a point field to a cell field by interpolating the nodal values at elements centroid.

Parameters:: field (numpy.ndarray) – The field to transform. Shape should be (number_of_points,). If a cell field is provided, it will be returned as is.
Returns:: The transformed cell field. Shape is (number_of_cells,).
Return type:: numpy.ndarray

rotation_matrix(theta: ndarray, phi: ndarray) → ndarray

Computes the rotation matrix for Cartesian - Spherical vector transformations at given colatitude (\(\theta\)) and longitude (\(\phi\)) angles such that:

\[\begin{split}\boldsymbol{R} = \begin{bmatrix} \sin(\theta) \cos(\phi) & \sin(\theta) \sin(\phi) & \cos(\theta) \\ \cos(\theta) \cos(\phi) & \cos(\theta) \sin(\phi) & -\sin(\theta) \\ -\sin(\phi) & \cos(\phi) & 0 \end{bmatrix}\end{split}\]

Parameters:

theta (numpy.ndarray) – The colatitude angles (angle from the \(z\)-axis).
phi (numpy.ndarray) – The longitude angles (angle from the \(x\)-axis in the \(xy\)-plane).

Returns:

The rotation matrix for Cartesian - Spherical vector transformations at the given angles. Shape is (N, 3, 3) where N is the length of the input angle arrays.

Return type:

numpy.ndarray

rotation_matrix_centroids() → ndarray: Computes the rotation matrix for Cartesian - Spherical vector transformations at the centroids of the cells of the mesh. Calls rotation_matrix using the attribute centroids_spherical.

rotation_matrix_vertices() → ndarray: Computes the rotation matrix for Cartesian - Spherical vector transformations at the vertices of the mesh. Calls rotation_matrix using the attribute points_spherical.

Transform Spherical coordinates to Cartesian coordinates such that:

\[\begin{split}\begin{split} x &= r \sin(\theta) \cos(\phi) \\ y &= r \sin(\theta) \sin(\phi) \\ z &= r \cos(\theta) \end{split}\end{split}\]

Parameters:

R (numpy.ndarray|float) – The radial distance(s).
lat (numpy.ndarray|float) – The colatitude(s) (angle(s) from the \(z\)-axis).
lon (numpy.ndarray|float) – The longitude(s) (angle(s) from the \(x\)-axis in the \(xy\)-plane).

Returns:

A tuple containing the Cartesian coordinates (x, y, z).

Return type:

tuple[numpy.ndarray|float, numpy.ndarray|float, numpy.ndarray|float]

vector_cartesian_to_spherical(v_cartesian: ndarray) → ndarray

Transform a vector field from Cartesian to Spherical coordinates such that:

\[\mathbf{v}_{r} = \boldsymbol{R} \mathbf{v}_{x}\]

where \(\mathbf{v}_{r}\) is the vector field in spherical coordinates, \(\mathbf{v}_{x}\) is the vector field in Cartesian coordinates, and \(\boldsymbol{R}\) is the rotation matrix of the transformomation between Cartesian and Spherical coordinates (see rotation_matrix).

Parameters:: v_cartesian (numpy.ndarray) – The vector field in Cartesian coordinates. Shape should be (N, 3) where N can be either the number of points or the number of cells of the mesh.
Returns:: The vector field in Spherical coordinates. Shape is the same as the input field.
Return type:: numpy.ndarray

vector_spherical_to_cartesian(v_spherical: ndarray) → ndarray

Transform a vector field from Spherical to Cartesian coordinates such that:

\[\mathbf{v}_{x} = \boldsymbol{R}^T \mathbf{v}_{r}\]

where \(\mathbf{v}_{x}\) is the vector field in Cartesian coordinates, \(\mathbf{v}_{r}\) is the vector field in spherical coordinates, and \(\boldsymbol{R}^T\) is the transpose of the rotation matrix of the transformomation between Cartesian and Spherical coordinates (see rotation_matrix).

Parameters:: v_spherical (numpy.ndarray) – The vector field in Spherical coordinates. Shape should be (N, 3) where N can be either the number of points or the number of cells of the mesh.
Returns:: The vector field in Cartesian coordinates. Shape is the same as the input field.
Return type:: numpy.ndarray

Specific mesh types

class stagpyviz.ShellMesh(*args, **kwargs)

Class representing the surface mesh of a spherical shell, which can be loaded from a VTU file or created from a point cloud. Inherits from UnstructuredSphere and pyvista.UnstructuredGrid.

The mesh is represented as a collection of triangular facets defining \(\mathcal P_1\) elements in \(\mathbb R^3\) using the class P1_2D_R3.

The constructor can be called in three ways:

With a single argument that is a path to a VTU file containing an unstructured grid. The mesh will be loaded from the file, and point and cell data will be copied if available.
With a single argument that is a pyvista.UnstructuredGrid object. The mesh will be created from the provided unstructured grid, and point and cell data will be copied if available.
With a single argument that is a numpy.ndarray of shape (number_of_points, 3) containing the coordinates of points. The mesh will be created as the convex hull of the provided points, and a "neighbors" array representing the connectivity of the facets will be added to the cell data.

Attributes:

elements

The isoparametric element class representing triangular facets in 3D space.

Type:: P1_2D_R3

points_normal

An array of shape (number_of_points, 3) containing the normal vectors at each point. The normals are evaluated as the normalized position vectors of the points.

Type:: numpy.ndarray

cells_normal

An array of shape (number_of_cells, 3) containing the normal vectors at each cell (facet). The normals are evaluated as the normalized position vectors of the cell centroids.

Type:: numpy.ndarray

neighbors

An array of shape (number_of_cells, 3) containing the indices of neighboring cells across each facet. This attribute is only available if the mesh was created from a points array or if the VTU file from which the mesh was loaded contained this information, and is stored in the cell data under the key "neighbors".

Type:: numpy.ndarray

cells_area

An array of shape (number_of_cells,) containing the area of each triangular facet. The area is computed using the determinant of the Jacobian of the transformation from the reference element to the physical element such that: \(A_e = \frac{1}{2} |\det(J_e)|\).

Type:: numpy.ndarray

Methods:

connectivity(*args, **kwargs) → tuple[ndarray, ndarray]

Overload of pyvista connectivity filter to return an array of region IDs for each point in the mesh instead of a new mesh containing only the connected region. Non connected points will be assigned a region ID of -1. See pyvista connectivity documentation for details on the supported arguments and options.

Returns:: A tuple of arrays containing the region IDs for each point and cell in the mesh. Points and cells that are not part of any connected region are assigned a value of -1.
Return type:: tuple[numpy.ndarray, numpy.ndarray]

integrate_over_cell(field: ndarray, rule: str = '1pt') → ndarray

Compute the numerical integral of a field over each cell of the mesh using a 1 point quadrature rule. The field can be either a point field (defined at mesh points) or a cell field (defined at mesh cells). If the field is a point field, it will be interpolated to cell centroids using the shape functions of the elements evaluated at the cell centroids before integration. The integral over each cell is computed such that:

\[I_e = \int_{\Omega_e} \phi_e \, dS \approx \phi_e A_e\]

where \(\phi_e\) is the value of the field at the cell centroid and \(A_e\) is the area of the cell.

Parameters:: field (numpy.ndarray) – A 1D array containing the values of the field to integrate, either at points or at cells.
Returns:: A 1D array of shape (number_of_cells,) containing the integral of the field over each cell.
Return type:: numpy.ndarray

class stagpyviz.YinYangMesh(*args, **kwargs)

A class to represent an unstructured mesh of a sphere reconstructed from a StagYY Yin-Yang grid. The mesh is discretized as a collection of wedge (prismatic) elements \(\mathcal W_1\).

To reconstruct the mesh from a raw binary file containing the Yin-Yang grid we process as follows:

Extract the yin and yang grids
Determine the indices of the points that are not in the overlapping region between the two grids
Exploiting the fact that the grid is structured in the radial direction even after removing the overlapping points and that the number of points per radial layer is always the same, we reshape the coordinates of the points in the Yin and Yang grids into 2D arrays of shape (points_per_layer, n_radial_layers)
Generate a surface mesh of the outer shell of the sphere using the points of the last radial layer of both grids using the class ShellMesh
Generate the volume mesh by extruding the surface mesh radially and connecting the nodes of the surface mesh to the nodes of the inner layers by constructing wedge elements.

The mesh can be generated directly from the raw binary file or from a VTU file containing the mesh and the fields.

Parameters:

rawbin (pathlib.Path|str) – The path to the raw binary file containing the Yin-Yang grid (output of StagYY)
scaling (Scaling) – Scaling object containing the scaling factors and units to apply to the coordinates of the points in the mesh.

Attributes:

yin

Dictionary containing the coordinates of the points in the Yin grid. The keys are "x", "y", "z" for Cartesian coordinates and "R", "theta", "phi" for spherical coordinates. Contain the entire grid including the overlapping region. Arrays of shape (nx, ny, nz).

Type:: dict
Canonical:: stagpyviz.YinYangMesh.yin

yang

Dictionary containing the coordinates of the points in the Yang grid. The keys are "x", "y", "z" for Cartesian coordinates and "R", "theta", "phi" for spherical coordinates. Contain the entire grid including the overlapping region. Arrays of shape (nx, ny, nz).

Type:: dict
Canonical:: stagpyviz.YinYangMesh.yang

good_indices

Boolean array of shape (nx*ny*nz,) to flag the points that are not in the overlapping region between the Yin and Yang grids.

Type:: numpy.ndarray
Canonical:: stagpyviz.YinYangMesh.good_indices

points_per_layer

Number of points per radial layer in the Yin (or Yang) grid after removing the overlapping points.

Type:: int
Canonical:: stagpyviz.YinYangMesh.points_per_layer

surface_mesh

A surface mesh of the outer shell of the sphere generated from the points of the last radial layer of both grids.

Type:: ShellMesh
Canonical:: stagpyviz.YinYangMesh.surface_mesh

elements

An instance of the class Wedge3D to represent the wedge elements of the mesh and perform operations on them.

Type:: Wedge3D
Canonical:: stagpyviz.YinYangMesh.elements

grid_dimensions

Tuple of the dimensions of the Yin (or Yang) grid in the x, y and z directions.

Type:: tuple[int,int,int]
Canonical:: stagpyviz.YinYangMesh.grid_dimensions

grid_npoints

Total number of points in the Yin (or Yang) grid.

Type:: int
Canonical:: stagpyviz.YinYangMesh.grid_npoints

yin_radial_idx

2D array of shape (points_per_layer, n_radial_layers) containing the indices of the points in the Yin grid reshaped radially.

Type:: numpy.ndarray
Canonical:: stagpyviz.YinYangMesh.yin_radial_idx

yang_radial_idx

2D array of shape (points_per_layer, n_radial_layers) containing the indices of the points in the Yang grid reshaped radially.

Type:: numpy.ndarray
Canonical:: stagpyviz.YinYangMesh.yang_radial_idx

surface_idx

1D array of shape (2*points_per_layer,) containing the indices of the points in the surface mesh (last radial layer of both grids).

Type:: numpy.ndarray
Canonical:: stagpyviz.YinYangMesh.surface_idx

surface_cells

1D array of shape (number_of_surface_cells,) containing the indices of the cells in the surface mesh (last radial layer of both grids).

Type:: numpy.ndarray
Canonical:: stagpyviz.YinYangMesh.surface_cells

cells_Jacobian

3D array of shape (number_of_cells, 3, 3) containing the Jacobian matrices of the transformation from the reference element to the physical element for each cell in the mesh.

Type:: numpy.ndarray
Canonical:: stagpyviz.YinYangMesh.cells_Jacobian

Methods:

add_field(name: str, values: ndarray) → None

Add a field defined on the Yin and Yang grids to the mesh. The field can be either a scalar field or a vector field.

Parameters:

name (str) – The name of the field to be added to the mesh. This field will then be accessible in the mesh as self[name].
values (numpy.ndarray) – The field values defined on the Yin and Yang grids. For a scalar field, the array should have shape (nx, ny, nz, 2), where the first three dimensions correspond to the grid dimensions and the last dimension corresponds to the Yin and Yang grids (0 for Yin and 1 for Yang). For a vector field, the array should have shape (3, nx, ny, nz, 2), where the first dimension corresponds to the vector components, the next three dimensions correspond to the grid dimensions and the last dimension corresponds to the Yin and Yang grids (0 for Yin and 1 for Yang).

add_fields(fields: dict) → None

Add multiple fields defined on the Yin and Yang grids to the mesh.

Parameters:: fields (dict) – A dictionary containing the fields to be added to the mesh. The keys of the dictionary are the names of the fields and the values are the field values defined on the Yin and Yang grids. The field values should be in the same format as described in the add_field method.

compute_gradient(field: ndarray) → ndarray

Compute the Cartesian gradient of a scalar field defined on the mesh using the shape functions of the wedge elements.

Parameters:: field (numpy.ndarray) – The scalar field defined on the mesh for which the gradient is to be computed. The array should have shape (number_of_points,).
Returns:: The gradient of the field at the centroids of the elements. The array has shape (number_of_cells, 3), where the last dimension corresponds to the x, y and z components of the gradient.
Return type:: numpy.ndarray

get_radial_cells(k: int) → ndarray

Get the cells indices of the cells in the k-th radial layer of the reconstructed unstructured mesh. Useful to extract cell values at specific radial layers.

Parameters:: k (int) – The index of the radial layer for which to get the cells indices. The index starts from 1 for the innermost layer and goes up to n-1 for the outermost layer. Should never be outside of the range [1, n-1] where n is the number of radial layers in the grid (grid_dimensions[2]).
Returns:: A 1D array of shape (number_of_cells_in_layer,) containing the indices of the cells in the k-th radial layer of both grids.
Return type:: numpy.ndarray

get_radial_indices(k: int) → ndarray

Get the nodes indices of the points in the k-th radial layer of the reconstructed unstructured mesh. Useful to extract nodal values at specific radial layers.

Parameters:: k (int) – The index of the radial layer for which to get the points indices. The index starts from 0 for the innermost layer and goes up to n-1 for the outermost layer. Should never be outside of the range [0, n-1] where n is the number of radial layers in the grid (grid_dimensions[2]).
Returns:: A 1D array of shape (2*:py:attr:points_per_layer,) containing the indices of the points in the k-th radial layer of both grids.
Return type:: numpy.ndarray

integrate_over_cell(field: ndarray, rule: str = '1pt') → ndarray

Integrate a field over cells for each cell. The integration can be performed using different quadrature rules specified by the rule parameter, by default the 1 point quadrature rule is used (see Wedge3D.quadrature_1pt). Integration is perfomed using integrate_field.

Parameters:

field (np.ndarray) – The field to be integrated. The array should have shape (number_of_points, ...) for a point field or (number_of_cells, ...) for a cell field.
rule (str) –
Quadrature rule to be used. Available:
- "1pt": 1 point quadrature rule using Wedge3D.quadrature_1pt. Cell fields automatically use this method.
- "3x2pt": 3x2 points rule using method quadrature_6pt.

Returns:

The integrated field over the cells of shape (number_of_cells, ...).

Return type:

numpy.ndarray

reconstruct_velocity(velocity_raw: ndarray) → None

Reconstruct the velocity field from the raw velocity components defined on the Yin and Yang grids. First the raw coordinates of the points are transformed such that:

\[\begin{split}\begin{split} r &= z_{\text{raw}} + r_{\text{cmb}} \\ \theta &= \frac{\pi}{4} - x_{\text{raw}} \\ \phi &= y_{\text{raw}} - \frac{3\pi}{4} \end{split}\end{split}\]

Then the velocity components are transformed from the Yin-Yang grid coordinates to Cartesian coordinates using the following relations:

\[\begin{split}\begin{split} \mathbf{v}_{\text{yin}} &= \begin{bmatrix} v_0 \sin(\theta)\cos(\phi) - v_1\sin(\phi) + v_2\cos(\theta)\cos(\phi) \\ v_0 \sin(\theta)\sin(\phi) + v_1\cos(\phi) + v_2\cos(\theta)\sin(\phi) \\ -v_0 \cos(\theta) + v_2\sin(\theta) \end{bmatrix} \\ \mathbf{v}_{\text{yang}} &= \begin{bmatrix} -v_{\text{yin}_x} \\ v_{\text{yin}_z} \\ v_{\text{yin}_y} \end{bmatrix} \end{split}\end{split}\]

where \(v_0\), \(v_1\) and \(v_2\) are the raw velocity components defined on the Yin and Yang grids, \(\mathbf{v}_{\text{yin}}\) and \(\mathbf{v}_{\text{yang}}\) are the velocity vectors in Cartesian coordinates on the Yin and Yang grids respectively.

Parameters:: velocity_raw (numpy.ndarray) – The raw velocity components defined on the Yin and Yang grids. The array should have shape (3, nx, ny, nz, 2), where the first dimension corresponds to the velocity components, the next three dimensions correspond to the grid dimensions and the last dimension corresponds to the Yin and Yang grids (0 for Yin and 1 for Yang).

Warning

This is an in-place operation that modifies the input array velocity_raw.

reshape_radially(field: ndarray) → ndarray

Reshape a field defined on the Yin or Yang grid into a 2D array of shape (points_per_layer, n_radial_layers).

Parameters:: field (numpy.ndarray) – 1D array of shape (npoints,) containing the values of the field defined on the Yin or Yang grid.
Returns:: 2D array of shape (points_per_layer, n_radial_layers) containing the values of the field reshaped radially.
Return type:: numpy.ndarray