Implement the face mass matrices for tensor products

modepy/matrices.py

@@ -242,7 +242,7 @@ def mass_matrix(basis, nodes):
     return la.inv(inverse_mass_matrix(basis, nodes))
 
 
-class _FaceMap:
+class _SimplexFaceMap:
     def __init__(self, face_vertices):
         """
         :arg face_vertices: an array of shape ``[dim, npts]``, where *npts*
@@ -252,33 +252,68 @@ def __init__(self, face_vertices):
         if npts != vol_dim:
             raise ValueError("face_vertices has wrong shape")
 
-        self.origin = face_vertices[:, 0]
-        self.span = face_vertices[:, 1:] - self.origin[:, np.newaxis]
+        self.origin = face_vertices[:, 0].reshape(-1, 1)
+        self.span = face_vertices[:, 1:] - self.origin
 
         self.face_dim = vol_dim - 1
 
     def __call__(self, points):
-        return (self.origin[:, np.newaxis]
-                + np.einsum("ad,dn->an", self.span, points*0.5 + 0.5))
+        return self.origin + np.einsum("ad,dn->an", self.span, points*0.5 + 0.5)
 
 
-def modal_face_mass_matrix(trial_basis, order, face_vertices, test_basis=None):
+class _HypercubeFaceMap:
+    def __init__(self, face_vertices):
+        """
+        :arg face_vertices: an array of shape ``[dim, npts]``, where *npts*
+            should equal `2**(dim - 1)`.
+        """
+        vol_dim, npts = face_vertices.shape
+        if npts != 2**(vol_dim-1):
+            raise ValueError("face_vertices has wrong shape")
+
+        self.origin = face_vertices[:, 0].reshape(-1, 1)
+        self.span = face_vertices[:, -2:0:-1] - self.origin
+
+        self.face_dim = vol_dim - 1
+
+    def __call__(self, points):
+        return self.origin + np.einsum("ad,dn->an", self.span, points*0.5 + 0.5)
+
+
+def modal_face_mass_matrix(trial_basis, order, face_vertices,
+        test_basis=None, domain=None):
     """
     :arg face_vertices: an array of shape ``[dim, npts]``, where *npts*
         should equal *dim*.
+    :arg domain: identifier for the reference element, can be one of
+        `"simplex"` or `"hypercube"`.
 
     .. versionadded :: 2016.1
+
+    .. versionchanged:: 2020.5
+
+        Added *domain* parameter and support for :math:`[-1, 1]^d` domains.
     """
 
     if test_basis is None:
         test_basis = trial_basis
 
-    fmap = _FaceMap(face_vertices)
+    if domain is None:
+        domain = "simplex"
+
+    if domain == "simplex":
+        from modepy.quadrature.grundmann_moeller import \
+                GrundmannMoellerSimplexQuadrature
+        fmap = _SimplexFaceMap(face_vertices)
+        quad = GrundmannMoellerSimplexQuadrature(order, fmap.face_dim)
+    elif domain == "hypercube":
+        from modepy.quadrature import LegendreGaussTensorProductQuadrature
+        fmap = _HypercubeFaceMap(face_vertices)
+        quad = LegendreGaussTensorProductQuadrature(fmap.face_dim, order)
+    else:
+        raise ValueError(f"unknown domain: '{domain}'")
 
-    from modepy.quadrature.grundmann_moeller import GrundmannMoellerSimplexQuadrature
-    quad = GrundmannMoellerSimplexQuadrature(order, fmap.face_dim)
     assert quad.exact_to > order*2
-
     mapped_nodes = fmap(quad.nodes)
 
     nrows = len(test_basis)
@@ -296,7 +331,7 @@ def modal_face_mass_matrix(trial_basis, order, face_vertices, test_basis=None):
 
 
 def nodal_face_mass_matrix(trial_basis, volume_nodes, face_nodes, order,
-        face_vertices, test_basis=None):
+        face_vertices, test_basis=None, domain=None):
     """
     :arg face_vertices: an array of shape ``[dim, npts]``, where *npts*
         should equal *dim*.
@@ -307,13 +342,22 @@ def nodal_face_mass_matrix(trial_basis, volume_nodes, face_nodes, order,
     if test_basis is None:
         test_basis = trial_basis
 
-    fmap = _FaceMap(face_vertices)
+    if domain is None:
+        domain = "simplex"
+
+    if domain == "simplex":
+        fmap = _SimplexFaceMap(face_vertices)
+    elif domain == "hypercube":
+        fmap = _HypercubeFaceMap(face_vertices)
+    else:
+        raise ValueError(f"unknown domain: '{domain}'")
 
     face_vdm = vandermonde(trial_basis, fmap(face_nodes))  # /!\ non-square
     vol_vdm = vandermonde(test_basis, volume_nodes)
 
     modal_fmm = modal_face_mass_matrix(
-            trial_basis, order, face_vertices, test_basis=test_basis)
+            trial_basis, order, face_vertices,
+            test_basis=test_basis, domain=domain)
     return la.inv(vol_vdm.T).dot(modal_fmm).dot(la.pinv(face_vdm))
 
 

modepy/modes.py

@@ -658,6 +658,10 @@ def tensor_product_basis(dims, basis_1d):
 
     .. versionadded:: 2017.1
     """
+    if dims == 0:
+        # NOTE: using to maintain consistency in the 0d case
+        return simplex_onb(dims, len(basis_1d))
+
     from pytools import generate_nonnegative_integer_tuples_below as gnitb
     return tuple(
             _TensorProductBasisFunction(order, [basis_1d[i] for i in order])
@@ -673,6 +677,10 @@ def grad_tensor_product_basis(dims, basis_1d, grad_basis_1d):
 
     .. versionadded:: 2020.2
     """
+    if dims == 0:
+        # NOTE: using to maintain consistency in the 0d case
+        return grad_simplex_onb(dims, len(basis_1d))
+
     from pytools import (
             wandering_element,
             generate_nonnegative_integer_tuples_below as gnitb)

modepy/nodes.py

@@ -302,6 +302,10 @@ def tensor_product_nodes(dims, nodes_1d):
 
     .. versionadded:: 2017.1
     """
+    if dims == 0:
+        # NOTE: using this to maintain consistency in the 0d case
+        return warp_and_blend_nodes(dims, 1)
+
     nnodes_1d = len(nodes_1d)
     result = np.empty((dims,) + (nnodes_1d,) * dims)
     for d in range(dims):

modepy/quadrature/__init__.py

@@ -82,3 +82,29 @@ def __init__(self, quad, left, right):
         Quadrature.__init__(self,
                 left + (quad.nodes+1) / 2 * length,
                 quad.weights * half_length)
+
+
+class TensorProductQuadrature(Quadrature):
+    def __init__(self, dims, quad):
+        """
+        :arg quad: a :class:`Quadrature` class for one-dimensional intervals.
+        """
+
+        from modepy.nodes import tensor_product_nodes
+        x = tensor_product_nodes(dims, quad.nodes)
+        from itertools import product
+        w = np.fromiter(
+                (np.prod(w) for w in product(quad.weights, repeat=dims)),
+                dtype=np.float,
+                count=quad.weights.size**dims)
+        assert w.size == x.shape[1]
+
+        super().__init__(x, w)
+        self.exact_to = quad.exact_to
+
+
+class LegendreGaussTensorProductQuadrature(TensorProductQuadrature):
+    def __init__(self, dims, N, backend=None):      # noqa: N803
+        from modepy.quadrature.jacobi_gauss import LegendreGaussQuadrature
+        super().__init__(
+                dims, LegendreGaussQuadrature(N, backend=backend))

modepy/tools.py

@@ -146,6 +146,10 @@ def inverse(self):
         """The inverse :class:`AffineMap` of *self*."""
         return AffineMap(la.inv(self.a), -la.solve(self.a, self.b))
 
+# }}}
+
+
+# {{{ simplex coordinate mapping
 
 EQUILATERAL_TO_UNIT_MAP = {
         1: AffineMap([[1]], [0]),
@@ -229,6 +233,20 @@ def barycentric_to_equilateral(bary):
 # }}}
 
 
+# {{{ hypercube coordinate mapping
+
+def hypercube_unit_vertices(dims):
+    from pytools import flatten, generate_nonnegative_integer_tuples_below as gnitb
+    vertices_01 = np.fromiter(
+            flatten(gnitb(2, dims)),
+            dtype=np.float64,
+            count=dims * 2**dims)
+
+    return -1.0 + 2.0 * vertices_01.reshape(-1, dims)
+
+# }}}
+
+
 def pick_random_simplex_unit_coordinate(rng, dims):
     offset = 0.05
     base = -1 + offset