Tensor Product Additions

meshmode/discretization/poly_element.py

@@ -41,19 +41,25 @@
 .. autoclass:: PolynomialRecursiveNodesElementGroup
 .. autoclass:: PolynomialEquidistantSimplexElementGroup
 .. autoclass:: PolynomialGivenNodesElementGroup
+
+.. autoclass:: InterpolatoryQuadratureTensorProductElementGroup
 .. autoclass:: LegendreGaussLobattoTensorProductElementGroup
 .. autoclass:: EquidistantTensorProductElementGroup
 
 Group factories
 ^^^^^^^^^^^^^^^
 
 .. autoclass:: ElementGroupFactory
+.. autoclass:: OrderAndTypeBasedGroupFactory
+
 .. autoclass:: InterpolatoryQuadratureSimplexGroupFactory
 .. autoclass:: QuadratureSimplexGroupFactory
 .. autoclass:: PolynomialWarpAndBlendGroupFactory
 .. autoclass:: PolynomialRecursiveNodesGroupFactory
 .. autoclass:: PolynomialEquidistantSimplexGroupFactory
 .. autoclass:: PolynomialGivenNodesGroupFactory
+
+.. autoclass:: InterpolatoryQuadratureTensorProductGroupFactory
 .. autoclass:: LegendreGaussLobattoTensorProductGroupFactory
 """
 
@@ -346,46 +352,103 @@ def unit_nodes(self):
 
 # {{{ concrete element groups for tensor product elements
 
-class LegendreGaussLobattoTensorProductElementGroup(PolynomialElementGroupBase):
+class _TensorProductElementGroupBase(PolynomialElementGroupBase):
+    # {{{
+
+    def basis_1d(self):
+        from modepy.modes import jacobi
+        from functools import partial
+        return tuple(partial(jacobi, 0, 0, i) for i in range(self.order + 1))
+
+    def grad_basis_1d(self):
+        from modepy.modes import grad_jacobi
+        from functools import partial
+        return tuple(partial(grad_jacobi, 0, 0, i) for i in range(self.order + 1))
+
+    def unit_nodes_1d(self):
+        raise NotImplementedError
+
+    # }}}
+
     def is_orthogonal_basis(self):
         return True
 
     def basis(self):
-        from modepy.modes import tensor_product_basis, jacobi
-        from functools import partial
-        return tensor_product_basis(
-                self.dim, tuple(
-                    partial(jacobi, 0, 0, i)
-                    for i in range(self.order+1)))
+        from modepy.modes import tensor_product_basis
+        return tensor_product_basis(self.dim, self.basis_1d())
 
     def grad_basis(self):
-        raise NotImplementedError()
+        from modepy.modes import grad_tensor_product_basis
+        return grad_tensor_product_basis(self.dim,
+                self.basis_1d(), self.grad_basis_1d())
 
     @property
     @memoize_method
     def unit_nodes(self):
         from modepy.nodes import tensor_product_nodes
-        from modepy.quadrature.jacobi_gauss import legendre_gauss_lobatto_nodes
-        return tensor_product_nodes(
-                self.dim, legendre_gauss_lobatto_nodes(self.order))
+        return tensor_product_nodes(self.dim, self.unit_nodes_1d())
 
     @memoize_method
     def from_mesh_interp_matrix(self):
+        from modepy.modes import legendre_tensor_product_basis
         meg = self.mesh_el_group
         return mp.resampling_matrix(
-                self.basis(),
+                legendre_tensor_product_basis(self.dim, self.order),
                 self.unit_nodes,
                 meg.unit_nodes)
 
 
-class EquidistantTensorProductElementGroup(
-        LegendreGaussLobattoTensorProductElementGroup):
+class InterpolatoryQuadratureTensorProductElementGroup(
+        _TensorProductElementGroupBase):
+    """Elemental discretization supplying a high-order quadrature rule
+    with a number of nodes matching the number of polynomials in the tensor
+    product basis, hence usable for differentiation and interpolation.
+
+    No interpolation nodes are present on the boundary of the hypercube.
+    """
+
+    @memoize_method
+    def _quadrature_rule_1d(self):
+        from modepy.quadrature import LegendreGaussQuadrature
+        return LegendreGaussQuadrature(self.order)
+
+    def unit_nodes_1d(self):
+        return self._quadrature_rule().nodes
+
     @property
     @memoize_method
-    def unit_nodes(self):
-        from modepy.nodes import tensor_product_nodes, equidistant_nodes
-        return tensor_product_nodes(
-                self.dim, equidistant_nodes(1, self.order)[0])
+    def weights(self):
+        from modepy.nodes import tensor_product_nodes
+        return tensor_product_nodes(self.dim, self._quadrature_rule().weights)
+
+
+class LegendreGaussLobattoTensorProductElementGroup(
+        _TensorProductElementGroupBase):
+    """Elemental discretization supplying a high-order quadrature rule
+    with a number of nodes matching the number of polynomials in the tensor
+    product basis, hence usable for differentiation and interpolation.
+    Nodes are present on the boundary of the hypercube.
+
+    Uses :func:`~modepy.quadrature.jacobi_gauss.legendre_gauss_lobatto_nodes`.
+    """
+
+    def unit_nodes_1d(self):
+        from modepy.quadrature.jacobi_gauss import legendre_gauss_lobatto_nodes
+        return legendre_gauss_lobatto_nodes(self.order)
+
+
+class EquidistantTensorProductElementGroup(_TensorProductElementGroupBase):
+    """Elemental discretization supplying a high-order quadrature rule
+    with a number of nodes matching the number of polynomials in the tensor
+    product basis, hence usable for differentiation and interpolation.
+    Nodes are present on the boundary of the hypercube.
+
+    Uses :func:`~modepy.equidistant_nodes`.
+    """
+
+    def unit_nodes_1d(self):
+        from modepy.nodes import equidistant_nodes
+        return equidistant_nodes(1, self.order)[0]
 
 # }}}
 
@@ -483,15 +546,26 @@ def __call__(self, mesh_el_group, index):
 
         return PolynomialGivenNodesElementGroup(
                 mesh_el_group, self.order, self.unit_nodes, index)
+
 # }}}
 
 
+# {{{ group factories for tensor products
+
+class InterpolatoryQuadratureTensorProductGroupFactory(
+        HomogeneousOrderBasedGroupFactory):
+    mesh_group_class = _MeshTensorProductElementGroup
+    group_class = InterpolatoryQuadratureTensorProductElementGroup
+
+
 class LegendreGaussLobattoTensorProductGroupFactory(
         HomogeneousOrderBasedGroupFactory):
     mesh_group_class = _MeshTensorProductElementGroup
     group_class = LegendreGaussLobattoTensorProductElementGroup
 
 # }}}
 
+# }}}
+
 
 # vim: fdm=marker

meshmode/mesh/__init__.py

@@ -366,7 +366,7 @@ def face_vertex_indices(self):
                 (1, 2, 3)
                 )
         else:
-            raise NotImplementedError("dim=%d" % self.dim)
+            raise NotImplementedError(f"dimension {self.dim}")
 
     def vertex_unit_coordinates(self):
         from modepy.tools import unit_vertices
@@ -412,10 +412,33 @@ def __init__(self, order, vertex_indices, nodes,
                 element_nr_base, node_nr_base, unit_nodes)
 
     def face_vertex_indices(self):
-        raise NotImplementedError()
+        if self.dim == 1:
+            return (
+                (0,),
+                (1,),
+                )
+        elif self.dim == 2:
+            return (
+                    (0, 1),
+                    (3, 2),
+                    (2, 0),
+                    (1, 3),
+                    )
+        elif self.dim == 3:
+            return (
+                    (0, 2, 1, 3),
+                    (4, 6, 5, 7),
+                    (0, 4, 5, 1),
+                    (2, 6, 7, 3),
+                    (0, 4, 6, 2),
+                    (1, 5, 7, 3),
+                    )
+        else:
+            raise NotImplementedError(f"dimension {self.dim}")
 
     def vertex_unit_coordinates(self):
-        raise NotImplementedError()
+        from pytools import generate_nonnegative_integer_tuples_below as gnitb
+        return np.array(list(gnitb(2, self.dim)), dtype=np.float64)*2.0 - 1.0
 
 # }}}
 
@@ -971,15 +994,22 @@ def __ne__(self, other):
 
 # {{{ node-vertex consistency test
 
-def _test_node_vertex_consistency_simplex(mesh, mgrp, tol):
+def _test_node_vertex_consistency_resampling(mesh, mgrp, tol):
     if mesh.vertices is None:
         return True
 
     if mgrp.nelements == 0:
         return True
 
+    if isinstance(mgrp, SimplexElementGroup):
+        basis = mp.simplex_best_available_basis(mgrp.dim, mgrp.order)
+    elif isinstance(mgrp, TensorProductElementGroup):
+        basis = mp.legendre_tensor_product_basis(mgrp.dim, mgrp.order)
+    else:
+        raise TypeError(f"unsupported group type: {type(mgrp).__name__}")
+
     resampling_mat = mp.resampling_matrix(
-            mp.simplex_best_available_basis(mgrp.dim, mgrp.order),
+            basis,
             mgrp.vertex_unit_coordinates().T,
             mgrp.unit_nodes)
 
@@ -1013,8 +1043,8 @@ def _test_node_vertex_consistency(mesh, tol):
     """
 
     for mgrp in mesh.groups:
-        if isinstance(mgrp, SimplexElementGroup):
-            assert _test_node_vertex_consistency_simplex(mesh, mgrp, tol)
+        if isinstance(mgrp, (SimplexElementGroup, TensorProductElementGroup)):
+            assert _test_node_vertex_consistency_resampling(mesh, mgrp, tol)
         else:
             from warnings import warn
             warn("not implemented: node-vertex consistency check for '%s'"

meshmode/mesh/generation.py

@@ -767,7 +767,7 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
 
     if dim == 1:
         if mesh_type is not None:
-            raise ValueError("unsupported mesh_type")
+            raise ValueError(f"unsupported mesh type: '{mesh_type}'")
 
         for i in range(shape[0]-1):
             # a--b
@@ -801,7 +801,7 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
             pass
 
         else:
-            raise ValueError("unsupported mesh_type")
+            raise ValueError(f"unsupported mesh type: '{mesh_type}'")
 
         for i in range(shape[0]-1):
             for j in range(shape[1]-1):
@@ -957,7 +957,8 @@ def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), n=(5, 5), order=1,
       on [a,b].
     :param boundary_tag_to_face: an optional dictionary for tagging boundaries.
         See :func:`generate_box_mesh`.
-    :param mesh_type: See :func:`generate_box_mesh`.
+    :param group_cls: see :func:`generate_box_mesh`.
+    :param mesh_type: see :func:`generate_box_mesh`.
     """
     if min(n) < 2:
         raise ValueError("need at least two points in each direction")
@@ -975,15 +976,15 @@ def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), n=(5, 5), order=1,
 
 # {{{ generate_warped_rect_mesh
 
-def generate_warped_rect_mesh(dim, order, n):
+def generate_warped_rect_mesh(dim, order, n, *, group_cls=None):
     """Generate a mesh of a warped line/square/cube. Mainly useful for testing
     functionality with curvilinear meshes.
     """
 
     assert dim in [1, 2, 3]
     mesh = generate_regular_rect_mesh(
             a=(-0.5,)*dim, b=(0.5,)*dim,
-            n=(n,)*dim, order=order)
+            n=(n,)*dim, order=order, group_cls=group_cls)
 
     def m(x):
         result = np.empty_like(x)