From 4ecbecaaff90ad262de48d6c45249a57a81540f5 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sat, 27 Apr 2024 16:08:38 +0200 Subject: [PATCH] mesh: make check_mesh_consistency raise instead of assert --- meshmode/mesh/__init__.py | 208 ++++++++++++++++++++++++++++---------- meshmode/mesh/io.py | 4 +- 2 files changed, 157 insertions(+), 55 deletions(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 39e4021c..23dd9442 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -44,6 +44,7 @@ .. autoclass:: Mesh .. autofunction:: make_mesh .. autofunction:: check_mesh_consistency +.. autofunction:: is_mesh_consistent .. autoclass:: NodalAdjacency .. autoclass:: FacialAdjacencyGroup @@ -761,42 +762,149 @@ def check_mesh_consistency( * The facial adjacency shapes and dtypes. * The mesh orientation using :func:`~meshmode.mesh.processing.find_volume_mesh_element_orientations`. + + :arg node_vertex_consistency_tolerance: If *False*, do not check for + consistency between vertex and nodal data. If *None*, a default tolerance + based on the :class:`~numpy.dtype` of the *vertices* array will be used. + Otherwise, the given value is used as the tolerance. + :arg skip_element_orientation_test: If *False*, check that element + orientation is positive in volume meshes (i.e. ones where ambient and + topological dimension match). + + :raises InconsistentMeshError: when the mesh is found to be inconsistent in + some fashion. """ + from meshmode import ( + InconsistentAdjacencyError, InconsistentArrayDTypeError, + InconsistentMeshError) + if node_vertex_consistency_tolerance is not False: - assert _test_node_vertex_consistency( - mesh, node_vertex_consistency_tolerance) + _test_node_vertex_consistency(mesh, tol=node_vertex_consistency_tolerance) + + for i, g in enumerate(mesh.groups): + if g.vertex_indices is None: + continue - for g in mesh.groups: - if g.vertex_indices is not None: - assert g.vertex_indices.dtype == mesh.vertex_id_dtype + if g.vertex_indices.dtype != mesh.vertex_id_dtype: + raise InconsistentArrayDTypeError( + f"Group '{i}' attribute 'vertex_indices' has incorrect dtype: " + f"{g.vertex_indices.dtype!r} (expected mesh 'vertex_id_dtype' = " + f"{mesh.vertex_id_dtype!r})") nodal_adjacency = mesh._nodal_adjacency if nodal_adjacency: - assert nodal_adjacency.neighbors_starts.shape == (mesh.nelements + 1,) - assert len(nodal_adjacency.neighbors.shape) == 1 - assert nodal_adjacency.neighbors_starts.dtype == mesh.element_id_dtype - assert nodal_adjacency.neighbors.dtype == mesh.element_id_dtype + if nodal_adjacency.neighbors_starts.shape != (mesh.nelements + 1,): + raise InconsistentAdjacencyError( + "Nodal adjacency 'neighbors_starts' has incorrect shape: " + f"'{nodal_adjacency.neighbors_starts.shape}' (expected " + f"nelements + 1 = {mesh.nelements + 1})") + + if len(nodal_adjacency.neighbors.shape) != 1: + raise InconsistentAdjacencyError( + "Nodal adjacency 'neighbors' have incorrect dim: " + f"{nodal_adjacency.neighbors.shape} (expected ndim = 1)") + + if nodal_adjacency.neighbors_starts.dtype != mesh.element_id_dtype: + raise InconsistentArrayDTypeError( + "Nodal adjacency 'neighbors_starts' has incorrect dtype: " + f"{nodal_adjacency.neighbors_starts.dtype!r} (expected mesh " + f"'element_id_dtype' = {mesh.element_id_dtype!r})") + + if nodal_adjacency.neighbors.dtype != mesh.element_id_dtype: + raise InconsistentArrayDTypeError( + "Nodal adjacency 'neighbors' has incorrect dtype: " + f"{nodal_adjacency.neighbors.dtype!r} (expected mesh " + f"'element_id_dtype' = {mesh.element_id_dtype!r})") facial_adjacency_groups = mesh._facial_adjacency_groups if facial_adjacency_groups: - assert len(facial_adjacency_groups) == len(mesh.groups) + if len(facial_adjacency_groups) != len(mesh.groups): + raise InconsistentAdjacencyError( + "Facial adjacency groups do not match mesh groups: " + f"{len(facial_adjacency_groups)} (expected {len(mesh.groups)})") + + for igrp, fagrp_list in enumerate(facial_adjacency_groups): + for ifagrp, fagrp in enumerate(fagrp_list): + if len(fagrp.elements.shape) != 1: + raise InconsistentAdjacencyError( + f"Facial adjacency {ifagrp} for group {igrp} has incorrect " + f"'elements' shape: {fagrp.elements.shape} " + "(expected ndim = 1)") - for fagrp_list in facial_adjacency_groups: - for fagrp in fagrp_list: nfagrp_elements, = fagrp.elements.shape - assert fagrp.element_faces.dtype == mesh.face_id_dtype - assert fagrp.element_faces.shape == (nfagrp_elements,) + if fagrp.element_faces.shape != (nfagrp_elements,): + raise InconsistentAdjacencyError( + f"Facial adjacency {ifagrp} for group {igrp} has incorrect " + f"'element_faces' shape: {fagrp.element_faces.shape} " + f"(expected 'elements.shape' = {fagrp.elements.shape})") + + if fagrp.element_faces.dtype != mesh.face_id_dtype: + raise InconsistentArrayDTypeError( + f"Facial adjacency {ifagrp} for group {igrp} has " + "incorrect 'element_faces' dtype: " + f"{fagrp.element_faces.dtype!r} (expected mesh " + f"'face_id_dtype' = {mesh.face_id_dtype!r})") if isinstance(fagrp, InteriorAdjacencyGroup): - assert fagrp.neighbors.dtype == mesh.element_id_dtype - assert fagrp.neighbors.shape == (nfagrp_elements,) - assert fagrp.neighbor_faces.dtype == mesh.face_id_dtype - assert fagrp.neighbor_faces.shape == (nfagrp_elements,) + if fagrp.neighbors.dtype != mesh.element_id_dtype: + raise InconsistentArrayDTypeError( + f"Facial adjacency {ifagrp} for group {igrp} has " + "incorrect 'neighbors' dtype: " + f"{fagrp.neighbors.dtype!r} (expected mesh " + f"'element_id_dtype' = {mesh.element_id_dtype!r})") + + if fagrp.neighbor_faces.dtype != mesh.face_id_dtype: + raise InconsistentArrayDTypeError( + f"Facial adjacency {ifagrp} for group {igrp} has " + "incorrect 'neighbor_faces' dtype: " + f"{fagrp.neighbor_faces.dtype!r} (expected mesh " + f"'face_id_dtype' = {mesh.face_id_dtype!r})") + + if fagrp.neighbors.shape != (nfagrp_elements,): + raise InconsistentAdjacencyError( + f"Facial adjacency {ifagrp} for group {igrp} has " + "incorrect 'neighbors' shape: " + f"{fagrp.neighbors.shape} (expected " + f"'elements.shape' = {fagrp.elements.shape})") + + if fagrp.neighbor_faces.shape != (nfagrp_elements,): + raise InconsistentAdjacencyError( + f"Facial adjacency {ifagrp} for group {igrp} has " + "incorrect 'neighbor_faces' shape: " + f"{fagrp.neighbor_faces.shape} (expected " + f"'elements.shape' = {fagrp.elements.shape})") from meshmode.mesh.processing import test_volume_mesh_element_orientations - if mesh.dim == mesh.ambient_dim and not skip_element_orientation_test: - assert test_volume_mesh_element_orientations(mesh) + if not skip_element_orientation_test: + if mesh.dim == mesh.ambient_dim: + if not test_volume_mesh_element_orientations(mesh): + raise InconsistentMeshError("Mesh has inconsistent orientations") + else: + warn("Cannot check element orientation for a mesh with " + "mesh.dim != mesh.ambient_dim", stacklevel=2) + + +def is_mesh_consistent( + mesh: "Mesh", + *, + node_vertex_consistency_tolerance: Optional[ + Union[Literal[False], float]] = None, + skip_element_orientation_test: bool = False, + ) -> bool: + """A boolean version of :func:`check_mesh_consistency`.""" + + from meshmode import InconsistentMeshError + + try: + check_mesh_consistency( + mesh, + node_vertex_consistency_tolerance=node_vertex_consistency_tolerance, + skip_element_orientation_test=skip_element_orientation_test) + except InconsistentMeshError: + return False + else: + return True def make_mesh( @@ -857,12 +965,8 @@ def make_mesh( :arg skip_tests: a flag used to skip any mesh consistency checks. This can be set to *True* in special situation, e.g. when loading a broken mesh that will be fixed later. - :arg node_vertex_consistency_tolerance: If *False*, do not check for - consistency between vertex and nodal data. If *None*, a default tolerance - based on the :class:`~numpy.dtype` of the *vertices* array will be used. - :arg skip_element_orientation_test: If *False*, check that element - orientation is positive in volume meshes (i.e. ones where ambient and - topological dimension match). + :arg node_vertex_consistency_tolerance: see :func:`check_mesh_consistency`. + :arg skip_element_orientation_test: see :func:`check_mesh_consistency`. """ vertex_id_dtype = np.dtype(vertex_id_dtype) if vertex_id_dtype.kind not in {"i", "u"}: @@ -1145,8 +1249,8 @@ def dim(self) -> int: def nvertices(self) -> int: """Number of vertices in the mesh, if available.""" if self.vertices is None: - from meshmode import DataUnavailable - raise DataUnavailable("vertices") + from meshmode import DataUnavailableError + raise DataUnavailableError("vertices") return self.vertices.shape[-1] @@ -1175,8 +1279,8 @@ def base_node_nrs(self): def vertex_dtype(self): """The :class:`~numpy.dtype` of the :attr:`~Mesh.vertices` array, if any.""" if self.vertices is None: - from meshmode import DataUnavailable - raise DataUnavailable("vertices") + from meshmode import DataUnavailableError + raise DataUnavailableError("vertices") return self.vertices.dtype @@ -1187,17 +1291,17 @@ def nodal_adjacency(self) -> NodalAdjacency: This property gets the :attr:`Mesh._nodal_adjacency` of the mesh. If the attribute value is *None*, the adjacency is computed and cached. - :raises DataUnavailable: if the nodal adjacency cannot be obtained. + :raises DataUnavailableError: if the nodal adjacency cannot be obtained. """ - from meshmode import DataUnavailable + from meshmode import DataUnavailableError nodal_adjacency = self._nodal_adjacency if nodal_adjacency is False: - raise DataUnavailable("Nodal adjacency is not available") + raise DataUnavailableError("Nodal adjacency is not available") if nodal_adjacency is None: if not self.is_conforming: - raise DataUnavailable( + raise DataUnavailableError( "Nodal adjacency can only be computed for conforming meshes" ) @@ -1254,17 +1358,17 @@ def facial_adjacency_groups( Note that element groups are not necessarily geometrically contiguous like the figure may suggest. - :raises DataUnavailable: if the facial adjacency cannot be obtained. + :raises DataUnavailableError: if the facial adjacency cannot be obtained. """ - from meshmode import DataUnavailable + from meshmode import DataUnavailableError facial_adjacency_groups = self._facial_adjacency_groups if facial_adjacency_groups is False: - raise DataUnavailable("Facial adjacency is not available") + raise DataUnavailableError("Facial adjacency is not available") if facial_adjacency_groups is None: if not self.is_conforming: - raise DataUnavailable( + raise DataUnavailableError( "Facial adjacency can only be computed for conforming meshes" ) @@ -1333,14 +1437,17 @@ def _mesh_group_node_vertex_error(mesh, mgrp): return map_vertices - grp_vertices -def _test_node_vertex_consistency_resampling(mesh, igrp, tol): +def _test_group_node_vertex_consistency_resampling( + mesh: Mesh, igrp: int, *, tol: Optional[float] = None) -> None: if mesh.vertices is None: - return True + return mgrp = mesh.groups[igrp] if mgrp.nelements == 0: - return True + return + + from meshmode import InconsistentVerticesError per_vertex_errors = _mesh_group_node_vertex_error(mesh, mgrp) per_element_vertex_errors = np.max( @@ -1359,32 +1466,27 @@ def _test_node_vertex_consistency_resampling(mesh, igrp, tol): if len(elements_above_tol) > 0: i_grp_elem = elements_above_tol[0] ielem = i_grp_elem + mesh.base_element_nrs[igrp] - from meshmode import InconsistentVerticesError + raise InconsistentVerticesError( - f"vertex consistency check failed for element {ielem}; " + f"Vertex consistency check failed for element {ielem}; " f"{per_element_vertex_errors[i_grp_elem]} >= " f"{per_element_tols[i_grp_elem]}") - return True +def _test_node_vertex_consistency( + mesh: Mesh, *, tol: Optional[float] = None) -> None: + """Ensure that order of by-index vertices matches that of mapped unit vertices. -def _test_node_vertex_consistency(mesh, tol): - """Ensure that order of by-index vertices matches that of mapped - unit vertices. + :raises InconsistentVerticesError: if the vertices are not consistent. """ - if not __debug__: - return True - for igrp, mgrp in enumerate(mesh.groups): if isinstance(mgrp, _ModepyElementGroup): - assert _test_node_vertex_consistency_resampling(mesh, igrp, tol) + _test_group_node_vertex_consistency_resampling(mesh, igrp, tol=tol) else: warn("Not implemented: node-vertex consistency check for " f"groups of type '{type(mgrp).__name__}'.", stacklevel=3) - return True - # }}} diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index 8c923bba..ce36a5d4 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -441,12 +441,12 @@ def group_to_json(group): "dim": group.dim, } - from meshmode import DataUnavailable + from meshmode import DataUnavailableError def nodal_adjacency_to_json(mesh): try: na = mesh.nodal_adjacency - except DataUnavailable: + except DataUnavailableError: return None return {