Test that connections that should be are actually permutations

meshmode/discretization/connection/direct.py

@@ -250,28 +250,36 @@ def _resample_point_pick_indices(self, actx: ArrayContext,
         :class:`pyopencl.array.Array` containing the index subset.
         """
 
-        mat = actx.to_numpy(actx.thaw(
-                self._resample_matrix(actx, to_group_index, ibatch_index)))
-
-        nrows, ncols = mat.shape
-        result = np.zeros(nrows, dtype=self.to_discr.mesh.element_id_dtype)
+        ibatch = self.groups[to_group_index].batches[ibatch_index]
+        from_grp = self.from_discr.groups[ibatch.from_group_index]
 
         if tol_multiplier is None:
             tol_multiplier = 50
 
-        tol = np.finfo(mat.dtype).eps * tol_multiplier
+        tol = np.finfo(ibatch.result_unit_nodes.dtype).eps * tol_multiplier
 
-        for irow in range(nrows):
-            one_indices, = np.where(np.abs(mat[irow] - 1) < tol)
-            zero_indices, = np.where(np.abs(mat[irow]) < tol)
+        dim, ntgt_nodes = ibatch.result_unit_nodes.shape
+        if dim == 0:
+            assert ntgt_nodes == 1
+            return actx.freeze(actx.from_numpy(np.array([0], dtype=np.int32)))
 
-            if len(one_indices) != 1:
-                return None
-            if len(zero_indices) != ncols - 1:
+        dist_vecs = (ibatch.result_unit_nodes.reshape(dim, -1, 1)
+                - from_grp.unit_nodes.reshape(dim, 1, -1))
+        dists = np.sqrt(np.sum(dist_vecs**2, axis=0))
+
+        result = np.zeros(ntgt_nodes, dtype=self.to_discr.mesh.element_id_dtype)
+
+        if tol_multiplier is None:
+            tol_multiplier = 500
+
+        for irow in range(ntgt_nodes):
+            close_indices, = np.where(dists[irow] < tol)
+
+            if len(close_indices) != 1:
                 return None
 
-            one_index, = one_indices
-            result[irow] = one_index
+            close_index, = close_indices
+            result[irow] = close_index
 
         return actx.freeze(actx.from_numpy(result))
 

meshmode/discretization/connection/opposite_face.py

@@ -193,7 +193,7 @@ def get_map_jacobian(unit_nodes):
             df_inv_resid[1] = 1/det * (-ata[0, 1] * atb[0] + ata[0, 0]*atb[1])
 
         else:
-            # The boundary of a 3D mesh is 2D, so that's the
+            # The boundary of a 3D mesh is 2D, so that (the case right above) is the
             # highest-dimensional case we genuinely care about.
             #
             # This stinks, performance-wise, because it's not vectorized.
@@ -222,6 +222,47 @@ def get_map_jacobian(unit_nodes):
         max_resid = np.max(np.abs(resid))
 
         if max_resid < tol:
+            if 0:
+                near_zero = np.abs(src_unit_nodes) < 1e-5
+                near_one = np.abs(src_unit_nodes-1) < 1e-5
+                near_mone = np.abs(src_unit_nodes+1) < 1e-5
+                reasonable = near_zero | near_one | near_mone
+                if np.any(~reasonable):
+                    unreasonable_els, unreasonable_nodes = np.where(
+                            np.any(~reasonable, axis=0))
+                    mapped = apply_map(src_unit_nodes)
+                    # print(src_unit_nodes[~reasonable])
+                    urel = unreasonable_els[0]
+                    tgt_urel = tgt_bdry_nodes[:, urel]
+                    src_urel = src_bdry_nodes[:, urel]
+
+                    dist_vecs = tgt_urel.reshape(3, -1, 1) - src_urel.reshape(3, 1, -1)
+                    dists = np.sqrt(np.sum(dist_vecs**2, axis=0))
+
+                    mat = (dists < 1e-12)
+                    from_indices = np.array([np.where(row)[0][0] for row in mat])
+
+                    sunodes = src_bdry_discr.groups[i_src_grp].unit_nodes
+
+                    print("GOOD UNIT NODES")
+                    print(sunodes[:, from_indices])
+                    print("BAD UNIT NODES")
+                    print(src_unit_nodes[:, urel])
+                    print()
+
+                    fixed_src_unit_nodes = src_unit_nodes.copy()
+                    fixed_src_unit_nodes[:, urel] = sunodes[:, from_indices]
+                    mapped_fixed = apply_map(fixed_src_unit_nodes)
+
+                    print("FIXED RESIDUAL")
+                    print(np.max(np.abs(mapped_fixed[:, urel]
+                        - tgt_bdry_nodes[:, urel])))
+                    print("UNFIXED RESIDUAL")
+                    print(np.max(np.abs(mapped[:, urel] - tgt_bdry_nodes[:, urel])))
+
+
+                    pu.db
+
             logger.debug("_make_cross_face_batches: gauss-newton: done, "
                     "final residual: %g", max_resid)
             break

meshmode/mesh/generation.py

@@ -1193,13 +1193,17 @@ def generate_warped_rect_mesh(dim, order, *, nelements_side=None,
 
     def m(x):
         result = np.empty_like(x)
-        result[0] = (
-                1.5*x[0] + np.cos(x[0])
-                + 0.1*np.sin(10*x[1]))
-        result[1] = (
-                0.05*np.cos(10*x[0])
-                + 1.3*x[1] + np.sin(x[1]))
-        if len(x) == 3:
+        if len(x) >= 2:
+            result[0] = (
+                    1.5*x[0] + np.cos(x[0])
+                    + 0.1*np.sin(10*x[1]))
+            result[1] = (
+                    0.05*np.cos(10*x[0])
+                    + 1.3*x[1] + np.sin(x[1]))
+        else:
+            result[0] = 1.5*x[0] + np.cos(x[0])
+
+        if len(x) >= 3:
             result[2] = x[2] + np.sin(x[0] / 2) / 2
         return result
 

test/test_connection.py

@@ -0,0 +1,109 @@
+__copyright__ = "Copyright (C) 2020 Andreas Kloeckner"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+from arraycontext import _acf  # noqa: F401
+from functools import partial
+import numpy as np  # noqa: F401
+import numpy.linalg as la  # noqa: F401
+
+import meshmode         # noqa: F401
+from arraycontext import (  # noqa
+        pytest_generate_tests_for_pyopencl_array_context
+        as pytest_generate_tests)
+
+from meshmode.mesh import SimplexElementGroup, TensorProductElementGroup
+from meshmode.discretization.poly_element import (
+        PolynomialWarpAndBlendGroupFactory,
+        PolynomialEquidistantSimplexGroupFactory,
+        LegendreGaussLobattoTensorProductGroupFactory,
+        PolynomialRecursiveNodesGroupFactory,
+        )
+from meshmode.discretization import Discretization
+from meshmode.discretization.connection import FACE_RESTR_ALL
+import meshmode.mesh.generation as mgen
+
+import pytest
+
+import logging
+logger = logging.getLogger(__name__)
+
+
+def connection_is_permutation(actx, conn):
+    for i_tgrp, cgrp in enumerate(conn.groups):
+        for i_batch, batch in enumerate(cgrp.batches):
+            if not len(batch.from_element_indices):
+                continue
+
+            point_pick_indices = conn._resample_point_pick_indices(
+                    actx, i_tgrp, i_batch)
+
+            if point_pick_indices is None:
+                return False
+
+    return True
+
+
+@pytest.mark.parametrize("group_factory", [
+        PolynomialWarpAndBlendGroupFactory,
+        PolynomialEquidistantSimplexGroupFactory,
+        LegendreGaussLobattoTensorProductGroupFactory,
+        partial(PolynomialRecursiveNodesGroupFactory, family="lgl"),
+        ])
+@pytest.mark.parametrize("dim", [2, 3])
+@pytest.mark.parametrize("order", [1, 2, 3, 4, 5])
+def test_bdry_restriction_is_permutation(actx_factory, group_factory, dim, order):
+    """Check that restriction to the boundary and opposite-face swap
+    for the element groups, orders and dimensions above is actually just
+    indirect access.
+    """
+    actx = actx_factory()
+
+    if group_factory is LegendreGaussLobattoTensorProductGroupFactory:
+        group_cls = TensorProductElementGroup
+    else:
+        group_cls = SimplexElementGroup
+
+    mesh = mgen.generate_warped_rect_mesh(dim, order=order, n=5,
+            group_cls=group_cls)
+
+    vol_discr = Discretization(actx, mesh, group_factory(order))
+    from meshmode.discretization.connection import (
+            make_face_restriction, make_opposite_face_connection)
+    bdry_connection = make_face_restriction(
+            actx, vol_discr, group_factory(order),
+            FACE_RESTR_ALL)
+
+    assert connection_is_permutation(actx, bdry_connection)
+
+    opp_face = make_opposite_face_connection(actx, bdry_connection)
+    assert connection_is_permutation(actx, opp_face)
+
+
+if __name__ == "__main__":
+    import sys
+    if len(sys.argv) > 1:
+        exec(sys.argv[1])
+    else:
+        from pytest import main
+        main([__file__])
+
+# vim: fdm=marker