Use elementwise size in node vertex consistency check (#316)

* use elementwise size in node vertex consistency check

* add comment explaining maximum(..., 1)

* use test based on maximum absolute vertex coordinate instead of element size, and add tests

* assume absolute and relative tolerances are the same for now

* add test for false negatives

* zero-size -> zero-D

* improve vertex consistency check error message

* use a specific exception type for consistency check

* improve tests

* just use relative tolerance and skip dim==0 case

* only do consistency check in debug run

* Reenable 0D tests, remove 1e-5 crit tol scaling, take into account abs tol in node-vert-consistency test tests

Co-authored-by: Andreas Kloeckner <[email protected]>

Author: majosm

meshmode/__init__.py

@@ -25,6 +25,7 @@
 .. exception:: Error
 .. exception:: DataUnavailable
 .. exception:: FileExistsError
+.. exception:: InconsistentVerticesError
 """
 
 from builtins import FileExistsError  # noqa: F401
@@ -39,6 +40,14 @@ class DataUnavailable(Error):
     pass
 
 
+class InconsistentVerticesError(Error):
+    """
+    Raised when an element's local-to-global mapping does not map the unit vertices
+    to the corresponding values in the mesh's *vertices* array.
+    """
+    pass
+
+
 def _acf():
     """A tiny undocumented function to pass to tests that take an ``actx_factory``
     argument when running them from the command line.

meshmode/mesh/__init__.py

@@ -1218,28 +1218,37 @@ def _mesh_group_node_vertex_error(mesh, mgrp):
     return map_vertices - grp_vertices
 
 
-def _test_node_vertex_consistency_resampling(mesh, mgrp, tol):
+def _test_node_vertex_consistency_resampling(mesh, igrp, tol):
     if mesh.vertices is None:
         return True
 
+    mgrp = mesh.groups[igrp]
+
     if mgrp.nelements == 0:
         return True
 
     per_vertex_errors = _mesh_group_node_vertex_error(mesh, mgrp)
-    per_element_vertex_errors = np.sqrt(
-        np.sum(np.sum((per_vertex_errors)**2, axis=0), axis=-1)
-    )
+    per_element_vertex_errors = np.max(
+        np.max(np.abs(per_vertex_errors), axis=-1), axis=0)
 
     if tol is None:
         tol = 1e3 * np.finfo(per_element_vertex_errors.dtype).eps
 
-    from meshmode.mesh.processing import find_bounding_box
+    grp_vertices = mesh.vertices[:, mgrp.vertex_indices]
+
+    coord_scales = np.max(np.max(np.abs(grp_vertices), axis=-1), axis=0)
 
-    bbox_min, bbox_max = find_bounding_box(mesh)
-    size = la.norm(bbox_max-bbox_min)
+    per_element_tols = tol + tol * coord_scales
 
-    max_el_vertex_error = np.max(per_element_vertex_errors)
-    assert max_el_vertex_error < tol*size, max_el_vertex_error
+    elements_above_tol, = np.where(per_element_vertex_errors >= per_element_tols)
+    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"{per_element_vertex_errors[i_grp_elem]} >= "
+            f"{per_element_tols[i_grp_elem]}")
 
     return True
 
@@ -1248,10 +1257,12 @@ def _test_node_vertex_consistency(mesh, tol):
     """Ensure that order of by-index vertices matches that of mapped
     unit vertices.
     """
+    if not __debug__:
+        return True
 
-    for mgrp in mesh.groups:
+    for igrp, mgrp in enumerate(mesh.groups):
         if isinstance(mgrp, _ModepyElementGroup):
-            assert _test_node_vertex_consistency_resampling(mesh, mgrp, tol)
+            assert _test_node_vertex_consistency_resampling(mesh, igrp, tol)
         else:
             from warnings import warn
             warn("not implemented: node-vertex consistency check for '%s'"

test/test_mesh.py

@@ -24,6 +24,7 @@
 """
 
 from functools import partial
+from dataclasses import replace
 import numpy as np
 import numpy.linalg as la
 import pytest
@@ -1130,6 +1131,156 @@ def test_mesh_element_group_constructor(shape_cls, group_cls):
 # }}}
 
 
+# {{{ test_node_vertex_consistency_check
+
+def test_node_vertex_consistency_check(actx_factory):
+    actx = actx_factory()
+
+    from meshmode import InconsistentVerticesError
+
+    dtype = np.float64
+    tol = 1e3 * np.finfo(dtype).eps
+
+    # Mesh bounds invariance
+
+    def gen_rect_mesh_with_perturbed_vertices(
+            a, b, nelements_per_axis, perturb_amount):
+        mesh_unperturbed = mgen.generate_regular_rect_mesh(
+            a=a, b=b, nelements_per_axis=nelements_per_axis)
+        return mesh_unperturbed.copy(  # noqa: F841
+            vertices=(
+                mesh_unperturbed.vertices
+                + perturb_amount*np.ones(mesh_unperturbed.vertices.shape)),
+            node_vertex_consistency_tolerance=tol)
+
+    # Find critical perturbation amount "w" such that vertices shifted by w works
+    # but 10*w doesn't
+    h = 1
+    nelems = 8
+    size = h*nelems
+    crit_perturb = None
+    for p in range(32):
+        w = h*10**(p-16)
+        try:
+            gen_rect_mesh_with_perturbed_vertices(
+                a=(-size/2,), b=(size/2,), nelements_per_axis=(nelems,),
+                perturb_amount=w)
+        except InconsistentVerticesError:
+            if p > 0:
+                crit_perturb = w/10
+            break
+    if crit_perturb is None:
+        raise RuntimeError("failed to find critical vertex perturbation amount")
+
+    # Scale invariance
+    nelems = 8
+    for p in range(32):
+        h = 10**(p-16)
+        size = h*nelems
+        perturb = crit_perturb*h
+        gen_rect_mesh_with_perturbed_vertices(
+            a=(-size/2,)*3, b=(size/2,)*3, nelements_per_axis=(nelems,)*3,
+            perturb_amount=perturb)
+        if 10*perturb > tol:
+            with pytest.raises(InconsistentVerticesError):
+                gen_rect_mesh_with_perturbed_vertices(
+                    a=(-size/2,)*3, b=(size/2,)*3, nelements_per_axis=(nelems,)*3,
+                    perturb_amount=10*perturb)
+
+    # Translation invariance
+    h = 1
+    nelems = 8
+    size = h*nelems
+    for p in range(10):
+        shift = 10**p-1
+        perturb = crit_perturb*(size/2 + shift)/(size/2)
+        gen_rect_mesh_with_perturbed_vertices(
+            a=(shift-size/2,)*3, b=(shift+size/2,)*3,
+            nelements_per_axis=(nelems,)*3,
+            perturb_amount=perturb)
+        with pytest.raises(InconsistentVerticesError):
+            gen_rect_mesh_with_perturbed_vertices(
+                a=(shift-size/2,)*3, b=(shift+size/2,)*3,
+                nelements_per_axis=(nelems,)*3,
+                perturb_amount=10*perturb)
+
+    # Aspect ratio invariance
+    h = 1
+    nelems = 8
+    size = h*nelems
+    for p in range(10):
+        h_x = 10**p * h
+        size_x = h_x * nelems
+        perturb = crit_perturb*h_x
+        gen_rect_mesh_with_perturbed_vertices(
+            a=(-size_x/2,) + (-size/2,)*2,
+            b=(size_x/2,) + (size/2,)*2,
+            nelements_per_axis=(nelems,)*3,
+            perturb_amount=perturb)
+        with pytest.raises(InconsistentVerticesError):
+            gen_rect_mesh_with_perturbed_vertices(
+                a=(-size_x/2,) + (-size/2,)*2,
+                b=(size_x/2,) + (size/2,)*2,
+                nelements_per_axis=(nelems,)*3,
+                perturb_amount=10*perturb)
+
+    # Mesh size relative to element size invariance
+    h = 1
+    for p in range(5):
+        nelems = 2**(5-p)
+        size = h*nelems
+        perturb = crit_perturb
+        gen_rect_mesh_with_perturbed_vertices(
+            a=(-size/2,)*3, b=(size/2,)*3, nelements_per_axis=(nelems,)*3,
+            perturb_amount=perturb)
+        with pytest.raises(InconsistentVerticesError):
+            gen_rect_mesh_with_perturbed_vertices(
+                a=(-size/2,)*3, b=(size/2,)*3, nelements_per_axis=(nelems,)*3,
+                perturb_amount=10*perturb)
+
+    # Zero-D elements
+    h = 1
+    nelems = 7
+    size = h*nelems
+    vol_mesh = mgen.generate_regular_rect_mesh(
+        a=(-size/2,), b=(size/2,),
+        nelements_per_axis=(nelems,))
+    from meshmode.discretization import Discretization
+    group_factory = default_simplex_group_factory(1, 1)
+    vol_discr = Discretization(actx, vol_mesh, group_factory)
+    from meshmode.discretization.connection import (
+        FACE_RESTR_ALL,
+        make_face_restriction)
+    make_face_restriction(
+        actx, vol_discr, group_factory, FACE_RESTR_ALL, per_face_groups=False)
+
+    # Zero-D elements at the origin
+    h = 1
+    nelems = 8
+    size = h*nelems
+    vol_mesh = mgen.generate_regular_rect_mesh(
+        a=(-size/2,), b=(size/2,),
+        nelements_per_axis=(nelems,))
+    group_factory = default_simplex_group_factory(1, 1)
+    vol_discr = Discretization(actx, vol_mesh, group_factory)
+    make_face_restriction(
+        actx, vol_discr, group_factory, FACE_RESTR_ALL, per_face_groups=False)
+
+    # Element vertex indices rotated
+    with pytest.raises(InconsistentVerticesError):
+        vol_mesh_unrotated = mgen.generate_regular_rect_mesh(
+            a=(-1,)*2, b=(1,)*2,
+            nelements_per_axis=(8,)*2)
+        vol_mesh = vol_mesh_unrotated.copy(  # noqa: F841
+            groups=[
+                replace(
+                    grp,
+                    vertex_indices=np.roll(grp.vertex_indices, 1, axis=1))
+                for grp in vol_mesh_unrotated.groups])
+
+# }}}
+
+
 # {{{ mesh boundary gluing
 
 @pytest.mark.parametrize("use_tree", [False, True])