Test that connections that should be permutations actually *are* permutations

examples/moving-geometry.py

@@ -140,7 +140,7 @@ def run(actx, *,
     # a bit of work when reconstructing after a time step
 
     if group_factory_name == "warp_and_blend":
-        group_factory_cls = poly.PolynomialWarpAndBlendGroupFactory
+        group_factory_cls = poly.PolynomialWarpAndBlend2DRestrictingGroupFactory
 
         unit_nodes = mp.warp_and_blend_nodes(ambient_dim - 1, mesh_order)
     elif group_factory_name == "quadrature":

examples/plot-connectivity.py

@@ -20,10 +20,10 @@ def main():
 
     from meshmode.discretization import Discretization
     from meshmode.discretization.poly_element import \
-            PolynomialWarpAndBlendGroupFactory
+            PolynomialWarpAndBlend3DRestrictingGroupFactory
 
     discr = Discretization(
-            actx, mesh, PolynomialWarpAndBlendGroupFactory(order))
+            actx, mesh, PolynomialWarpAndBlend3DRestrictingGroupFactory(order))
 
     from meshmode.discretization.visualization import make_visualizer
     vis = make_visualizer(actx, discr, order)

examples/simple-dg.py

@@ -74,8 +74,9 @@ def __init__(self, actx, mesh, order):
 
         from meshmode.discretization import Discretization
         from meshmode.discretization.poly_element import \
-                PolynomialWarpAndBlendGroupFactory
-        self.group_factory = PolynomialWarpAndBlendGroupFactory(order=order)
+                PolynomialWarpAndBlend2DRestrictingGroupFactory
+        self.group_factory = PolynomialWarpAndBlend2DRestrictingGroupFactory(
+                order=order)
         self.volume_discr = Discretization(actx, mesh, self.group_factory)
 
         assert self.volume_discr.dim == 2

experiments/refinement-playground.py

@@ -186,9 +186,10 @@ def refine_and_generate_chart_function(mesh, filename, function):
 
     from meshmode.discretization import Discretization
     from meshmode.discretization.poly_element import \
-            PolynomialWarpAndBlendGroupFactory
+            default_simplex_group_factory
     discr = Discretization(
-            cl_ctx, mesh, PolynomialWarpAndBlendGroupFactory(order))
+            cl_ctx, mesh, default_simplex_group_factory(
+                base_dim=mesh.dim, order=order))
     from meshmode.discretization.visualization import make_visualizer
     vis = make_visualizer(queue, discr, order)
     remove_if_exists("connectivity2.vtu")

meshmode/discretization/connection/direct.py

@@ -22,6 +22,7 @@
 
 
 import numpy as np
+import numpy.linalg as la
 
 import loopy as lp
 from pytools import memoize_in, keyed_memoize_method
@@ -250,28 +251,33 @@ 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 = la.norm(dist_vecs, axis=0, ord=2)
+
+        result = np.zeros(ntgt_nodes, dtype=self.to_discr.mesh.element_id_dtype)
+
+        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

@@ -70,7 +70,6 @@ def _make_cross_face_batches(actx,
 
     dim = src_grp.dim
     _, nelements, ntgt_unit_nodes = tgt_bdry_nodes.shape
-    assert tgt_bdry_nodes.shape == src_bdry_nodes.shape
 
     # {{{ invert face map (using Gauss-Newton)
 
@@ -101,7 +100,9 @@ def apply_map(unit_nodes):
         one_deviation = np.abs(np.sum(intp_coeffs, axis=0) - 1)
         assert (one_deviation < tol).all(), np.max(one_deviation)
 
-        return np.einsum("fet,aef->aet", intp_coeffs, src_bdry_nodes)
+        mapped = np.einsum("fet,aef->aet", intp_coeffs, src_bdry_nodes)
+        assert tgt_bdry_nodes.shape == mapped.shape
+        return mapped
 
     def get_map_jacobian(unit_nodes):
         # unit_nodes: (dim, nelements, ntgt_unit_nodes)

meshmode/discretization/poly_element.py

@@ -53,7 +53,8 @@
 
 .. autoclass:: InterpolatoryQuadratureSimplexElementGroup
 .. autoclass:: QuadratureSimplexElementGroup
-.. autoclass:: PolynomialWarpAndBlendElementGroup
+.. autoclass:: PolynomialWarpAndBlend2DRestrictingElementGroup
+.. autoclass:: PolynomialWarpAndBlend3DRestrictingElementGroup
 .. autoclass:: PolynomialRecursiveNodesElementGroup
 .. autoclass:: PolynomialEquidistantSimplexElementGroup
 .. autoclass:: PolynomialGivenNodesElementGroup
@@ -80,7 +81,8 @@
 
 .. autoclass:: InterpolatoryQuadratureSimplexGroupFactory
 .. autoclass:: QuadratureSimplexGroupFactory
-.. autoclass:: PolynomialWarpAndBlendGroupFactory
+.. autoclass:: PolynomialWarpAndBlend2DRestrictingGroupFactory
+.. autoclass:: PolynomialWarpAndBlend3DRestrictingGroupFactory
 .. autoclass:: PolynomialRecursiveNodesGroupFactory
 .. autoclass:: PolynomialEquidistantSimplexGroupFactory
 .. autoclass:: PolynomialGivenNodesGroupFactory
@@ -289,6 +291,16 @@ class PolynomialWarpAndBlendElementGroup(_MassMatrixQuadratureElementGroup):
 
     Uses :func:`modepy.warp_and_blend_nodes`.
     """
+    def __init__(self, mesh_el_group, order, index):
+        from warnings import warn
+        warn("PolynomialWarpAndBlendElementGroup is deprecated, since "
+                "the facial restrictions of the 3D nodes are not the 2D nodes. "
+                "It will go away in 2022. "
+                "Use PolynomialWarpAndBlend2DRestrictingElementGroup or "
+                "PolynomialWarpAndBlend3DRestrictingElementGroup instead.",
+                DeprecationWarning, stacklevel=2)
+        super().__init__(mesh_el_group, order, index)
+
     @property
     @memoize_method
     def _interp_nodes(self):
@@ -304,6 +316,66 @@ def _interp_nodes(self):
         return result
 
 
+class PolynomialWarpAndBlend2DRestrictingElementGroup(
+        _MassMatrixQuadratureElementGroup):
+    """Elemental discretization with a number of nodes matching the number of
+    polynomials in :math:`P^k`, hence usable for differentiation and
+    interpolation. Interpolation nodes edge-clustered for avoidance of Runge
+    phenomena. Nodes are present on the boundary of the simplex.
+    Provides nodes in two and fewer dimensions, based on the 2D
+    warp-and-blend nodes and their facial restrictions.
+
+    Uses :func:`modepy.warp_and_blend_nodes`.
+    """
+    @property
+    @memoize_method
+    def _interp_nodes(self):
+        dim = self.mesh_el_group.dim
+        if self.order == 0:
+            result = mp.warp_and_blend_nodes(dim, 1)
+            result = np.mean(result, axis=1).reshape(-1, 1)
+        elif dim >= 3:
+            raise ValueError(
+                    "PolynomialWarpAndBlend2DRestrictingElementGroup does not "
+                    f"provide nodes in {dim}D")
+        else:
+            result = mp.warp_and_blend_nodes(dim, self.order)
+
+        dim2, _ = result.shape
+        assert dim2 == dim
+        return result
+
+
+class PolynomialWarpAndBlend3DRestrictingElementGroup(
+        _MassMatrixQuadratureElementGroup):
+    """Elemental discretization with a number of nodes matching the number of
+    polynomials in :math:`P^k`, hence usable for differentiation and
+    interpolation. Interpolation nodes edge-clustered for avoidance of Runge
+    phenomena. Nodes are present on the boundary of the simplex.
+    Provides nodes in two and fewer dimensions, based on the 3D
+    warp-and-blend nodes and their facial restrictions.
+
+    Uses :func:`modepy.warp_and_blend_nodes`.
+    """
+    @property
+    @memoize_method
+    def _interp_nodes(self):
+        dim = self.mesh_el_group.dim
+        if self.order == 0:
+            result = mp.warp_and_blend_nodes(dim, 1)
+            result = np.mean(result, axis=1).reshape(-1, 1)
+        else:
+            result = mp.warp_and_blend_nodes(3, self.order)
+            tol = np.finfo(result.dtype).eps * 50
+            for d in range(dim, 3):
+                wanted = np.abs(result[d] - (-1)) < tol
+                result = result[:, wanted]
+
+            result = result[:dim]
+
+        return result
+
+
 class PolynomialRecursiveNodesElementGroup(_MassMatrixQuadratureElementGroup):
     """Elemental discretization with a number of nodes matching the number of
     polynomials in :math:`P^k`, hence usable for differentiation and
@@ -609,10 +681,33 @@ class QuadratureSimplexGroupFactory(HomogeneousOrderBasedGroupFactory):
 
 
 class PolynomialWarpAndBlendGroupFactory(HomogeneousOrderBasedGroupFactory):
+    def __init__(self, order):
+        from warnings import warn
+        warn("PolynomialWarpAndBlendGroupFactory is deprecated, since "
+                "the facial restrictions of the 3D nodes are not the 2D nodes. "
+                "It will go away in 2022. "
+                "Use PolynomialWarpAndBlend2DRestrictingGroupFactory or "
+                "PolynomialWarpAndBlend3DRestrictingGroupFactory instead.",
+                DeprecationWarning, stacklevel=2)
+
+        super().__init__(order)
+
     mesh_group_class = _MeshSimplexElementGroup
     group_class = PolynomialWarpAndBlendElementGroup
 
 
+class PolynomialWarpAndBlend2DRestrictingGroupFactory(
+        HomogeneousOrderBasedGroupFactory):
+    mesh_group_class = _MeshSimplexElementGroup
+    group_class = PolynomialWarpAndBlend2DRestrictingElementGroup
+
+
+class PolynomialWarpAndBlend3DRestrictingGroupFactory(
+        HomogeneousOrderBasedGroupFactory):
+    mesh_group_class = _MeshSimplexElementGroup
+    group_class = PolynomialWarpAndBlend3DRestrictingElementGroup
+
+
 class PolynomialRecursiveNodesGroupFactory(HomogeneousOrderBasedGroupFactory):
     def __init__(self, order, family):
         super().__init__(order)
@@ -674,4 +769,28 @@ class LegendreGaussLobattoTensorProductGroupFactory(
 # }}}
 
 
+# undocumented for now, mainly for internal use
+def default_simplex_group_factory(base_dim, order):
+    """
+    :arg base_dim: The dimension of the 'base' discretization to be used.
+        The returned group factory will also support creating lower-dimensional
+        discretizations.
+    """
+
+    try:
+        # recursivenodes is only importable in Python 3.8 since
+        # it uses :func:`math.comb`, so need to check if it can
+        # be imported.
+        import recursivenodes  # noqa: F401
+    except ImportError:
+        # If it cannot be imported, use warp-and-blend nodes.
+        if base_dim <= 2:
+            return PolynomialWarpAndBlend2DRestrictingGroupFactory(order)
+        elif base_dim == 3:
+            return PolynomialWarpAndBlend3DRestrictingGroupFactory(order)
+        else:
+            raise ValueError(f"no usable set of nodes found for {base_dim}D")
+
+    return PolynomialRecursiveNodesGroupFactory(order, family="lgl")
+
 # vim: fdm=marker

meshmode/dof_array.py

@@ -667,7 +667,7 @@ def flat_norm(ary, ord=None) -> float:
         actx = ary.array_context
         return la.norm(
                 [
-                    actx.np.linalg.norm(actx.np.ravel(ary, order="A"), ord=ord)
+                    actx.np.linalg.norm(actx.np.ravel(subary, order="A"), ord=ord)
                     for _, subary in serialize_container(ary)],
                 ord=ord)
 

meshmode/interop/firedrake/connection.py

@@ -39,8 +39,7 @@
 from meshmode.mesh.processing import get_simplex_element_flip_matrix
 
 from meshmode.discretization.poly_element import (
-    PolynomialWarpAndBlendGroupFactory,
-    PolynomialRecursiveNodesGroupFactory,
+    default_simplex_group_factory,
     ElementGroupFactory)
 from meshmode.discretization import (
     Discretization, InterpolatoryElementGroupBase)
@@ -623,14 +622,7 @@ def build_connection_from_firedrake(actx, fdrake_fspace, grp_factory=None,
         a :class:`~meshmode.discretization.poly_element.ElementGroupFactory`
         whose group class is a subclass of
         :class:`~meshmode.discretization.InterpolatoryElementGroupBase`.
-        If *None*, and :mod:`recursivenodes` can be imported,
-        a :class:`~meshmode.discretization.poly_element.\
-PolynomialRecursiveNodesGroupFactory` with ``"lgl"`` nodes is used.
-        Note that :mod:`recursivenodes` may not be importable
-        as it uses :func:`math.comb`, which is new in Python 3.8.
-        In the case that :mod:`recursivenodes` cannot be successfully
-        imported, a :class:`~meshmode.discretization.poly_element.\
-PolynomialWarpAndBlendGroupFactory` is used.
+        If *None*, and a default factory is automatically selected.
     :arg restrict_to_boundary: (optional)
         If not *None*, then must be one of the following:
 
@@ -659,6 +651,16 @@ def build_connection_from_firedrake(actx, fdrake_fspace, grp_factory=None,
                          "must be be "
                          "'Discontinuous Lagrange', not '%s'."
                          % ufl_elt.family())
+
+    # If only converting a portion of the mesh near the boundary, get
+    # *cells_to_use* as described in
+    # :func:`meshmode.interop.firedrake.mesh.import_firedrake_mesh`
+    cells_to_use = _get_cells_to_use(fdrake_fspace.mesh(),
+                                     restrict_to_boundary)
+
+    mm_mesh, orient = import_firedrake_mesh(fdrake_fspace.mesh(),
+                                            cells_to_use=cells_to_use)
+
     # Make sure grp_factory is the right type if provided, and
     # uses an interpolatory class.
     if grp_factory is not None:
@@ -676,17 +678,9 @@ def build_connection_from_firedrake(actx, fdrake_fspace, grp_factory=None,
                             % type(grp_factory.group_class))
     # If not provided, make one
     else:
-        degree = ufl_elt.degree()
-        try:
-            # recursivenodes is only importable in Python 3.8 since
-            # it uses :func:`math.comb`, so need to check if it can
-            # be imported
-            import recursivenodes  # noqa : F401
-            family = "lgl"  # L-G-Legendre
-            grp_factory = PolynomialRecursiveNodesGroupFactory(degree, family)
-        except ImportError:
-            # If cannot be imported, uses warp-and-blend nodes
-            grp_factory = PolynomialWarpAndBlendGroupFactory(degree)
+        grp_factory = default_simplex_group_factory(
+                mm_mesh.dim, ufl_elt.degree())
+
     # validate restrict_to_boundary, if present
     if restrict_to_boundary is not None:
         firedrake_bdy_ids = fdrake_fspace.mesh().exterior_facets.unique_markers
@@ -718,15 +712,7 @@ def build_connection_from_firedrake(actx, fdrake_fspace, grp_factory=None,
                             "restrict_to_boundary. Must be an int, a tuple "
                             "of ints, or the string 'on_boundary'.")
 
-    # If only converting a portion of the mesh near the boundary, get
-    # *cells_to_use* as described in
-    # :func:`meshmode.interop.firedrake.mesh.import_firedrake_mesh`
-    cells_to_use = _get_cells_to_use(fdrake_fspace.mesh(),
-                                     restrict_to_boundary)
-
     # Create to_discr
-    mm_mesh, orient = import_firedrake_mesh(fdrake_fspace.mesh(),
-                                            cells_to_use=cells_to_use)
     to_discr = Discretization(actx, mm_mesh, grp_factory)
 
     # get firedrake unit nodes and map onto meshmode reference element