diff --git a/packages/loop_common/src/loop_common/geometry/_bounding_box.py b/packages/loop_common/src/loop_common/geometry/_bounding_box.py index 65d584fb..b1a284a6 100644 --- a/packages/loop_common/src/loop_common/geometry/_bounding_box.py +++ b/packages/loop_common/src/loop_common/geometry/_bounding_box.py @@ -8,7 +8,6 @@ # from LoopStructural.utils.exceptions import LoopValueError from loop_common.math import rng -from loop_common.supports import StructuredGrid logger = getLogger(__name__) @@ -701,6 +700,9 @@ def structured_grid( name: str = "bounding_box", local_coordinates: bool = False, ): + # imported here because loop_common.supports builds on loop_common.geometry + from loop_common.supports import StructuredGrid + # python is passing a reference to the cell_data, vertex_data dicts so we need to # copy them to make sure that different instances of StructuredGrid are not sharing the same # underlying objects diff --git a/packages/loop_common/src/loop_common/geometry/_structured_grid_3d.py b/packages/loop_common/src/loop_common/geometry/_structured_grid_3d.py index a151dfb5..4dd39258 100644 --- a/packages/loop_common/src/loop_common/geometry/_structured_grid_3d.py +++ b/packages/loop_common/src/loop_common/geometry/_structured_grid_3d.py @@ -39,6 +39,9 @@ def __init__(self, origin=None, nsteps=None, step_vector=None, rotation_xy=None) self._rotation_xy[2, 2] = 1 self.rotation_xy = rotation_xy + def _on_geometry_change(self): + """Hook called after origin, step_vector, nsteps or maximum change; subclasses override.""" + @property def volume(self): return np.prod(self.maximum - self.origin) @@ -69,6 +72,7 @@ def nsteps(self, nsteps): change_factor = nsteps / self.nsteps self._step_vector /= change_factor self._nsteps = nsteps + self._on_geometry_change() @property def nsteps_cells(self): @@ -99,6 +103,7 @@ def step_vector(self, step_vector): newsteps = self._nsteps / change_factor self._nsteps = np.ceil(newsteps).astype(int) self._step_vector = step_vector + self._on_geometry_change() @property def origin(self): @@ -115,6 +120,7 @@ def origin(self, origin): logger.error(f"Cannot resize the grid. The proposed number of steps is {self._nsteps}, these must be all > 0") raise ValueError("Cannot resize the grid.") self._origin = origin + self._on_geometry_change() @property def maximum(self): @@ -130,6 +136,7 @@ def maximum(self, maximum): if np.any(~(self._nsteps > 0)): logger.error(f"Cannot resize the grid. The proposed number of steps is {self._nsteps}, these must be all > 0") raise ValueError("Cannot resize the grid.") + self._on_geometry_change() @property def n_nodes(self): @@ -211,27 +218,6 @@ def check_position(self, pos: np.ndarray) -> np.ndarray: raise ValueError("Position array needs to be a list of points or a point") return pos - def neighbour_global_indexes(self, mask=None, **kwargs): - indexes = None - if "indexes" in kwargs: - indexes = kwargs["indexes"] - if "indexes" not in kwargs: - gi = np.arange(self.n_nodes) - indexes = self.global_index_to_node_index(gi) - edge_mask = ( - (indexes[:, 0] > 0) & (indexes[:, 0] < self.nsteps[0] - 1) - & (indexes[:, 1] > 0) & (indexes[:, 1] < self.nsteps[1] - 1) - & (indexes[:, 2] > 0) & (indexes[:, 2] < self.nsteps[2] - 1) - ) - indexes = indexes[edge_mask, :].T - if indexes.ndim != 2: - logger.error("indexes.ndim = %s, expected 2", indexes.ndim) - return - if mask is None: - mask = np.array([[-1, 0, 1, -1, 0, 1, -1, 0, 1], [1, 1, 1, 0, 0, 0, -1, -1, -1]]) - neighbours = indexes[:, None, :] + mask[:, :, None] - return (neighbours[0, :, :] + self.nsteps[0, None, None] * neighbours[1, :, :]).astype(np.int64) - def cell_corner_indexes(self, cell_indexes: np.ndarray) -> np.ndarray: corner_indexes = np.zeros((cell_indexes.shape[0], 8, 3), dtype=np.int64) xcorner = np.array([0, 1, 0, 1, 0, 1, 0, 1]) diff --git a/packages/loop_common/src/loop_common/supports/_3d_base_structured.py b/packages/loop_common/src/loop_common/supports/_3d_base_structured.py index 56f653f9..55f606c1 100644 --- a/packages/loop_common/src/loop_common/supports/_3d_base_structured.py +++ b/packages/loop_common/src/loop_common/supports/_3d_base_structured.py @@ -2,6 +2,7 @@ import numpy as np +from loop_common.geometry import StructuredGrid3DGeometry from loop_common.logging import get_logger as getLogger from . import SupportType @@ -15,9 +16,14 @@ class LoopException(Exception): """Custom exception for LoopStructural errors.""" +class BaseStructuredSupport(StructuredGrid3DGeometry, BaseSupport): + """A structured grid used as an interpolation support. -class BaseStructuredSupport(BaseSupport): - """ """ + The grid geometry (origin, step vector, node counts, rotation, node and cell indexing) + comes from :class:`loop_common.geometry.StructuredGrid3DGeometry`; this class adds the + interpolation machinery. Note the constructor takes the number of **cells** + (``nsteps``) whereas the geometry class takes the number of nodes. + """ dimension = 3 @@ -33,71 +39,29 @@ def __init__( Parameters ---------- origin - 3d list or numpy array - nsteps - 3d list or numpy array of ints + nsteps - 3d list or numpy array of ints, the number of cells in each direction step_vector - 3d list or numpy array of int """ - # the geometry in the mesh can be calculated from the - # nsteps, step vector and origin - # we use property decorators to update these when different parts of - # the geometry need to change - # inisialise the private attributes - # cast to numpy array, to allow list like input - if origin is None: - origin = np.zeros(3) if nsteps is None: nsteps = np.array([10, 10, 10]) - if step_vector is None: - step_vector = np.ones(3) - origin = np.array(origin) nsteps = np.array(nsteps) - step_vector = np.array(step_vector) - - self.type = SupportType.BaseStructured - if np.any(step_vector == 0): - logger.warning(f"Step vector {step_vector} has zero values") if np.any(nsteps == 0): raise LoopException("nsteps cannot be zero") if np.any(nsteps < 0): raise LoopException("nsteps cannot be negative") - # if np.any(nsteps < 3): - # raise LoopException( - # "step vector cannot be less than 3. Try increasing the resolution of the interpolator" - # ) - self._nsteps = np.array(nsteps, dtype=int) + 1 - self._step_vector = np.array(step_vector) - self._origin = np.array(origin) + StructuredGrid3DGeometry.__init__( + self, + origin=origin, + nsteps=nsteps.astype(int) + 1, + step_vector=step_vector, + rotation_xy=rotation_xy, + ) + self.type = SupportType.BaseStructured self.supporttype = "Base" - self._rotation_xy = np.zeros((3, 3)) - self._rotation_xy[0, 0] = 1 - self._rotation_xy[1, 1] = 1 - self._rotation_xy[2, 2] = 1 - self.rotation_xy = rotation_xy self.interpolator = None - @property - def volume(self): - return np.prod(self.maximum - self.origin) - - def set_nelements(self, nelements) -> int: - box_vol = self.volume - ele_vol = box_vol / nelements - # calculate the step vector of a regular cube - step_vector = np.zeros(3) - - step_vector[:] = ele_vol ** (1.0 / 3.0) - - # number of steps is the length of the box / step vector - nsteps = np.ceil((self.maximum - self.origin) / step_vector).astype(int) - self.nsteps = nsteps - return self.n_elements - - def to_dict(self): - return { - "origin": self.origin, - "nsteps": self.nsteps, - "step_vector": self.step_vector, - "rotation_xy": self.rotation_xy, - } + def _on_geometry_change(self): + self.onGeometryChange() @abstractmethod def onGeometryChange(self): @@ -106,120 +70,6 @@ def onGeometryChange(self): def associateInterpolator(self, interpolator): self.interpolator = interpolator - @property - def nsteps(self): - return self._nsteps - - @nsteps.setter - def nsteps(self, nsteps): - # if nsteps changes we need to change the step vector - change_factor = nsteps / self.nsteps - self._step_vector /= change_factor - self._nsteps = nsteps - self.onGeometryChange() - - @property - def nsteps_cells(self): - return self.nsteps - 1 - - @property - def rotation_xy(self): - return self._rotation_xy - - @rotation_xy.setter - def rotation_xy(self, rotation_xy): - if rotation_xy is None: - return - if isinstance(rotation_xy, (float, int)): - rotation_xy = np.array( - [ - [ - np.cos(np.deg2rad(rotation_xy)), - -np.sin(np.deg2rad(rotation_xy)), - 0, - ], - [ - np.sin(np.deg2rad(rotation_xy)), - np.cos(np.deg2rad(rotation_xy)), - 0, - ], - [0, 0, 1], - ] - ) - rotation_xy = np.array(rotation_xy) - if rotation_xy.shape != (3, 3): - raise ValueError(f"Rotation matrix should be 3x3, not {rotation_xy.shape}") - self._rotation_xy = rotation_xy - - @property - def step_vector(self): - return self._step_vector - - @step_vector.setter - def step_vector(self, step_vector): - change_factor = step_vector / self._step_vector - newsteps = self._nsteps / change_factor - self._nsteps = np.ceil(newsteps).astype(int) - self._step_vector = step_vector - self.onGeometryChange() - - @property - def origin(self): - return self._origin - - @origin.setter - def origin(self, origin): - origin = np.array(origin) - length = self.maximum - origin - length /= self.step_vector - self._nsteps = np.ceil(length).astype(np.int64) - self._nsteps[self._nsteps == 0] = ( - 3 # need to have a minimum of 3 elements to apply the finite difference mask - ) - if np.any(~(self._nsteps > 0)): - logger.error( - f"Cannot resize the interpolation support. The proposed number of steps is {self._nsteps}, these must be all > 0" - ) - raise ValueError("Cannot resize the interpolation support.") - self._origin = origin - self.onGeometryChange() - - @property - def maximum(self): - return self.origin + self.nsteps_cells * self.step_vector - - @maximum.setter - def maximum(self, maximum): - """ - update the number of steps to fit new boundary - """ - maximum = np.array(maximum, dtype=float) - length = maximum - self.origin - length /= self.step_vector - self._nsteps = np.ceil(length).astype(np.int64) - self._nsteps[self._nsteps == 0] = 3 - if np.any(~(self._nsteps > 0)): - logger.error( - f"Cannot resize the interpolation support. The proposed number of steps is {self._nsteps}, these must be all > 0" - ) - raise ValueError("Cannot resize the interpolation support.") - self.onGeometryChange() - - @property - def n_nodes(self): - return np.prod(self.nsteps) - - @property - def n_elements(self): - return np.prod(self.nsteps_cells) - - @property - def elements(self): - global_index = np.arange(self.n_elements) - cell_indexes = self.global_index_to_cell_index(global_index) - - return self.global_node_indices(self.cell_corner_indexes(cell_indexes)) - def __str__(self): return ( f"LoopStructural interpolation support: {self.supporttype} \n" @@ -230,68 +80,9 @@ def __str__(self): f"Degrees of freedon {self.n_nodes}" ) - @property - def nodes(self): - max = self.origin + self.nsteps_cells * self.step_vector - if np.any(np.isnan(self.nsteps)): - raise ValueError("Cannot resize mesh nsteps is NaN") - if np.any(np.isnan(self.origin)): - raise ValueError("Cannot resize mesh origin is NaN") - - x = np.linspace(self.origin[0], max[0], self.nsteps[0]) - y = np.linspace(self.origin[1], max[1], self.nsteps[1]) - z = np.linspace(self.origin[2], max[2], self.nsteps[2]) - xx, yy, zz = np.meshgrid(x, y, z, indexing="ij") - return np.array([xx.flatten(order="F"), yy.flatten(order="F"), zz.flatten(order="F")]).T - - def rotate(self, pos): - """ """ - return np.einsum("ijk,ik->ij", self.rotation_xy[None, :, :], pos) - - def position_to_cell_index(self, pos: np.ndarray) -> tuple[np.ndarray, np.ndarray]: - """Get the indexes (i,j,k) of a cell - that a point is inside - - - Parameters - ---------- - pos : np.array - Nx3 array of xyz locations - - Returns - ------- - np.ndarray - N,3 i,j,k indexes of the cell that the point is in - """ - pos = self.check_position(pos) - inside = self.inside(pos) - cell_indexes = np.zeros((pos.shape[0], 3), dtype=int) - - x = pos[:, 0] - self.origin[None, 0] - y = pos[:, 1] - self.origin[None, 1] - z = pos[:, 2] - self.origin[None, 2] - cell_indexes[inside, 0] = np.floor(x[inside] / self.step_vector[None, 0]).astype(int) - cell_indexes[inside, 1] = np.floor(y[inside] / self.step_vector[None, 1]).astype(int) - cell_indexes[inside, 2] = np.floor(z[inside] / self.step_vector[None, 2]).astype(int) - - # Boundary points at maximum extents belong to the last cell. - cell_indexes[inside, 0] = np.clip(cell_indexes[inside, 0], 0, self.nsteps_cells[0] - 1) - cell_indexes[inside, 1] = np.clip(cell_indexes[inside, 1], 0, self.nsteps_cells[1] - 1) - cell_indexes[inside, 2] = np.clip(cell_indexes[inside, 2], 0, self.nsteps_cells[2] - 1) - - return cell_indexes, inside - def position_to_cell_global_index(self, pos): _ix, _iy, _iz = self.position_to_cell_index(pos) - def inside(self, pos): - # check whether point is inside box - pos = self.check_position(pos) - extent = np.maximum(np.abs(self.maximum - self.origin), 1.0) - tol = np.finfo(float).eps * 32.0 * extent - inside = np.all((pos >= (self.origin - tol)) & (pos <= (self.maximum + tol)), axis=1) - return inside - def check_position(self, pos: np.ndarray) -> np.ndarray: """[summary] @@ -325,64 +116,6 @@ def check_position(self, pos: np.ndarray) -> np.ndarray: raise ValueError("Position array needs to be a list of points or a point") return pos - def _global_indicies(self, indexes: np.ndarray, nsteps: np.ndarray) -> np.ndarray: - """ - Convert from cell indexes to global cell index - - Parameters - ---------- - indexes - - Returns - ------- - - """ - if len(indexes.shape) == 1: - raise ValueError("Cell indexes needs to be Nx3") - if indexes.shape[-1] != 3: - raise ValueError("Last dimensions should be ijk indexing") - original_shape = indexes.shape - indexes = indexes.reshape(-1, 3) - gi = ( - indexes[:, 0] - + nsteps[None, 0] * indexes[:, 1] - + nsteps[None, 0] * nsteps[None, 1] * indexes[:, 2] - ) - return gi.reshape(original_shape[:-1]) - - def cell_corner_indexes(self, cell_indexes: np.ndarray) -> np.ndarray: - """ - Returns the indexes of the corners of a cell given its location xi, - yi, zi - - Parameters - ---------- - x_cell_index - y_cell_index - z_cell_index - - Returns - ------- - - """ - - corner_indexes = np.zeros((cell_indexes.shape[0], 8, 3), dtype=int) - - xcorner = np.array([0, 1, 0, 1, 0, 1, 0, 1]) - ycorner = np.array([0, 0, 1, 1, 0, 0, 1, 1]) - zcorner = np.array([0, 0, 0, 0, 1, 1, 1, 1]) - corner_indexes[:, :, 0] = ( - cell_indexes[:, None, 0] + corner_indexes[:, :, 0] + xcorner[None, :] - ) - corner_indexes[:, :, 1] = ( - cell_indexes[:, None, 1] + corner_indexes[:, :, 1] + ycorner[None, :] - ) - corner_indexes[:, :, 2] = ( - cell_indexes[:, None, 2] + corner_indexes[:, :, 2] + zcorner[None, :] - ) - - return corner_indexes - def position_to_cell_corners(self, pos): """Get the global indices of the vertices (corners) of the cell containing each point. @@ -438,80 +171,6 @@ def node_indexes_to_position(self, node_indexes: np.ndarray) -> np.ndarray: xyz = xyz.reshape(original_shape) return xyz - def global_index_to_cell_index(self, global_index): - """ - Convert from global indexes to xi,yi,zi - - Parameters - ---------- - global_index - - Returns - ------- - - """ - # determine the ijk indices for the global index. - # remainder when dividing by nx = i - # remained when dividing modulus of nx by ny is j - cell_indexes = np.zeros((global_index.shape[0], 3), dtype=int) - cell_indexes[:, 0] = global_index % self.nsteps_cells[0, None] - cell_indexes[:, 1] = global_index // self.nsteps_cells[0, None] % self.nsteps_cells[1, None] - cell_indexes[:, 2] = ( - global_index // self.nsteps_cells[0, None] // self.nsteps_cells[1, None] - ) - return cell_indexes - - def global_index_to_node_index(self, global_index): - """ - Convert from global indexes to xi,yi,zi - - Parameters - ---------- - global_index - - Returns - ------- - - """ - # determine the ijk indices for the global index. - # remainder when dividing by nx = i - # remained when dividing modulus of nx by ny is j - original_shape = global_index.shape - global_index = global_index.reshape(-1) - local_indexes = np.zeros((global_index.shape[0], 3), dtype=int) - local_indexes[:, 0] = global_index % self.nsteps[0, None] - local_indexes[:, 1] = global_index // self.nsteps[0, None] % self.nsteps[1, None] - local_indexes[:, 2] = global_index // self.nsteps[0, None] // self.nsteps[1, None] - return local_indexes.reshape(*original_shape, 3) - - def global_node_indices(self, indexes) -> np.ndarray: - """ - Convert from node indexes to global node index - - Parameters - ---------- - indexes - - Returns - ------- - - """ - return self._global_indicies(indexes, self.nsteps) - - def global_cell_indices(self, indexes) -> np.ndarray: - """ - Convert from cell indexes to global cell index - - Parameters - ---------- - indexes - - Returns - ------- - - """ - return self._global_indicies(indexes, self.nsteps_cells) - @property def element_size(self): return np.prod(self.step_vector) diff --git a/packages/loop_common/tests/data/structured_support_reference.npz b/packages/loop_common/tests/data/structured_support_reference.npz new file mode 100644 index 00000000..43a4b079 Binary files /dev/null and b/packages/loop_common/tests/data/structured_support_reference.npz differ diff --git a/packages/loop_common/tests/test_structured_support_reference.py b/packages/loop_common/tests/test_structured_support_reference.py new file mode 100644 index 00000000..459b3c4d --- /dev/null +++ b/packages/loop_common/tests/test_structured_support_reference.py @@ -0,0 +1,238 @@ +"""Pin the behaviour of the structured interpolation supports. + +``StructuredGrid``, ``TetMesh`` and ``P2TetMesh`` share ``BaseStructuredSupport``. This test +records a wide set of their outputs (geometry, node/element arrays, point location, evaluation, +regularisation operators, setters, ``onGeometryChange`` calls, constructor validation) on a few +grids and compares them with ``data/structured_support_reference.npz``. + +The reference was generated from the code *before* ``BaseStructuredSupport`` was changed to extend +``StructuredGrid3DGeometry``, so it checks that refactor (and later ones) did not change behaviour. + +Regenerate after an intended change, and review the diff: + + python packages/loop_common/tests/test_structured_support_reference.py --update +""" +import sys +import warnings +from pathlib import Path + +import numpy as np +import pytest +from loop_common.supports import P2TetMesh, StructuredGrid, TetMesh + +REFERENCE = Path(__file__).parent / "data" / "structured_support_reference.npz" + +# (origin, cells per axis, step vector). Step vectors are floats on purpose. +GRIDS = { + "unit": ([0.0, 0.0, 0.0], [4, 4, 4], [1.0, 1.0, 1.0]), + "aniso": ([10.0, -3.0, 2.0], [6, 4, 3], [0.5, 2.0, 1.5]), + "small": ([0.0, 0.0, 0.0], [2, 2, 2], [1.0, 1.0, 1.0]), +} +CLASSES = {"StructuredGrid": StructuredGrid, "TetMesh": TetMesh, "P2TetMesh": P2TetMesh} +GEOMETRY_ATTRIBUTES = ( + "nsteps", + "nsteps_cells", + "maximum", + "origin", + "step_vector", + "n_nodes", + "n_elements", + "volume", + "nodes", + "elements", + "rotation_xy", + "element_size", + "element_scale", +) +REGULARISATION = {"dxx": 1, "dyy": 1, "dzz": 1, "dxy": 1, "dxz": 1, "dyz": 1} + + +def _build(cls, origin, cells, step): + if cls is P2TetMesh: + return cls(origin=origin, nsteps=cells, step_vector=step) + return cls(origin=origin, nsteps_cells=cells, step_vector=step) + + +def _flatten(value, prefix, out): + """Store ``value`` under ``prefix`` (tuples and dicts become ``prefix.``).""" + if isinstance(value, dict): + for key in sorted(value): + _flatten(value[key], f"{prefix}.{key}", out) + elif isinstance(value, (tuple, list)) and not all(np.isscalar(v) for v in value): + for i, item in enumerate(value): + _flatten(item, f"{prefix}.{i}", out) + else: + out[prefix] = np.asarray(value) + + +def _record(out, key, call): + try: + _flatten(call(), key, out) + except Exception as exc: # noqa: BLE001 - the exception type is part of the behaviour + out[key] = np.asarray(f"EXC:{type(exc).__name__}") + + +def collect() -> dict[str, np.ndarray]: + """Run every probe and return ``{key: array}``.""" + out: dict[str, np.ndarray] = {} + rng = np.random.default_rng(3) + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + for cls_name, cls in CLASSES.items(): + for grid_name, (origin, cells, step) in GRIDS.items(): + prefix = f"{cls_name}/{grid_name}" + grid = _build(cls, origin, cells, step) + calls = [] + original = type(grid).onGeometryChange + type(grid).onGeometryChange = lambda self, _o=original, _c=calls: ( + _c.append(1), + _o(self), + )[1] + try: + _probe(out, prefix, grid, origin, cells, step, rng, calls) + finally: + type(grid).onGeometryChange = original + for cells in ([0, 3, 3], [-1, 3, 3]): + _record( + out, + f"validation/{cells}", + lambda cells=cells: StructuredGrid( + origin=[0.0, 0.0, 0.0], nsteps_cells=cells, step_vector=[1.0, 1.0, 1.0] + ), + ) + return out + + +def _probe(out, prefix, grid, origin, cells, step, rng, calls): + for name in GEOMETRY_ATTRIBUTES: + _record(out, f"{prefix}/{name}", lambda name=name: getattr(grid, name)) + _record(out, f"{prefix}/to_dict", grid.to_dict) + + low = np.array(origin) - 1.0 + high = np.array(origin) + np.array(cells) * np.array(step) + 1.0 + anywhere = rng.uniform(low, high, (120, 3)) + # points exactly on the boundary: every corner and the centre of every face + corners = np.array(np.meshgrid(*zip(origin, grid.maximum), indexing="ij")).reshape(3, -1).T + centre = 0.5 * (np.array(origin) + grid.maximum) + faces = [] + for axis in range(3): + for bound in (origin[axis], grid.maximum[axis]): + point = centre.copy() + point[axis] = bound + faces.append(point) + anywhere = np.vstack([anywhere, corners, np.array(faces)]) + interior = rng.uniform(np.array(origin) + 1e-3, grid.maximum - 1e-3, (80, 3)) + _record(out, f"{prefix}/inside", lambda: grid.inside(anywhere)) + _record(out, f"{prefix}/position_to_cell_index", lambda: grid.position_to_cell_index(anywhere)) + _record(out, f"{prefix}/position_to_cell_corners", lambda: grid.position_to_cell_corners(interior)) + _record(out, f"{prefix}/position_to_cell_vertices", lambda: grid.position_to_cell_vertices(interior)) + _record(out, f"{prefix}/check_position_1d", lambda: grid.check_position(np.array([1.0, 2.0, 3.0]))) + _record(out, f"{prefix}/check_position_list", lambda: grid.check_position([[1.0, 2.0, 3.0]])) + _record(out, f"{prefix}/rotate", lambda: grid.rotate(anywhere)) + + cell_index = np.arange(grid.n_elements) + node_index = np.arange(grid.n_nodes) + _record(out, f"{prefix}/global_index_to_cell_index", lambda: grid.global_index_to_cell_index(cell_index)) + node_ijk = grid.global_index_to_node_index(node_index) + _record(out, f"{prefix}/global_index_to_node_index", lambda: node_ijk) + _record(out, f"{prefix}/node_indexes_to_position", lambda: grid.node_indexes_to_position(node_ijk[:20])) + if cls_has_ijk(node_ijk): + _record(out, f"{prefix}/global_node_indices", lambda: grid.global_node_indices(node_ijk)) + _record( + out, + f"{prefix}/global_cell_indices", + lambda: grid.global_cell_indices(grid.global_index_to_cell_index(cell_index)), + ) + + values = rng.normal(size=grid.n_nodes) + _record(out, f"{prefix}/evaluate_value", lambda: grid.evaluate_value(interior, values)) + _record(out, f"{prefix}/evaluate_gradient", lambda: grid.evaluate_gradient(interior, values)) + _record(out, f"{prefix}/get_element_for_location", lambda: grid.get_element_for_location(interior)) + _record( + out, + f"{prefix}/get_element_gradient_for_location", + lambda: grid.get_element_gradient_for_location(interior), + ) + _record(out, f"{prefix}/neighbour_global_indexes", lambda: grid.neighbour_global_indexes()) + if hasattr(grid, "get_operators"): + _record( + out, + f"{prefix}/get_operators", + lambda: {k: v[:2] for k, v in grid.get_operators(REGULARISATION).items()}, + ) + _record(out, f"{prefix}/vtk_n_cells", lambda: np.array([grid.vtk().n_cells])) + + # setters, applied in sequence; each records the state it leaves behind + before = len(calls) + + def state(): + return {"nsteps": grid.nsteps, "step_vector": grid.step_vector, "origin": grid.origin} + + grid.set_nelements(200) + _record(out, f"{prefix}/after_set_nelements", lambda: {**state(), "n_elements": grid.n_elements}) + grid.nsteps = np.array(grid.nsteps) + 1 + _record(out, f"{prefix}/after_nsteps", state) + grid.step_vector = np.array(grid.step_vector) * 2 + _record(out, f"{prefix}/after_step_vector", state) + grid.origin = np.array(grid.origin) - 1 + _record(out, f"{prefix}/after_origin", state) + grid.maximum = np.array(grid.maximum) + 3 + _record(out, f"{prefix}/after_maximum", lambda: {**state(), "maximum": grid.maximum}) + out[f"{prefix}/onGeometryChange_calls"] = np.array([len(calls) - before]) + grid.rotation_xy = 30 + _record(out, f"{prefix}/after_rotation", lambda: grid.rotation_xy) + + +def cls_has_ijk(node_ijk) -> bool: + """The tet meshes return something other than an (N, 3) array here; the ijk helpers do not apply.""" + return isinstance(node_ijk, np.ndarray) and node_ijk.ndim == 2 and node_ijk.shape[-1] == 3 + + +def _assert_same(key, got, want): + if want.dtype.kind in "US": + assert got.dtype.kind in "US" and str(got) == str(want), f"{key}: {got} != {want}" + return + assert got.dtype.kind not in "US", f"{key}: now raises {got}" + assert got.shape == want.shape, f"{key}: shape {got.shape} != {want.shape}" + assert got.dtype.kind == want.dtype.kind, f"{key}: dtype {got.dtype} != {want.dtype}" + if want.dtype.kind == "f": + np.testing.assert_allclose(got, want, rtol=1e-9, atol=1e-12, err_msg=key) + else: + np.testing.assert_array_equal(got, want, err_msg=key) + + +@pytest.fixture(scope="module") +def current(): + return collect() + + +@pytest.fixture(scope="module") +def reference(): + with np.load(REFERENCE) as data: + return {key: data[key] for key in data.files} + + +def test_same_probes_are_recorded(current, reference): + assert sorted(current) == sorted(reference) + + +@pytest.mark.parametrize("cls_name", list(CLASSES)) +def test_structured_support_behaviour_is_unchanged(cls_name, current, reference): + keys = [k for k in reference if k.startswith(f"{cls_name}/")] + assert keys + for key in keys: + _assert_same(key, current[key], reference[key]) + + +def test_constructor_validation_is_unchanged(current, reference): + for key in [k for k in reference if k.startswith("validation/")]: + _assert_same(key, current[key], reference[key]) + + +if __name__ == "__main__": + if sys.argv[1:] != ["--update"]: + raise SystemExit(__doc__) + REFERENCE.parent.mkdir(exist_ok=True) + result = collect() + np.savez_compressed(REFERENCE, **result) + print(f"wrote {len(result)} arrays to {REFERENCE}")