Tensor product 1d nodes & basis functions

.ci/install-for-firedrake.sh

@@ -5,25 +5,21 @@ sudo apt update
 sudo apt upgrade -y
 sudo apt install time
 
-# Otherwise we get missing CL symbols
-export PYOPENCL_CL_PRETEND_VERSION=2.1
-sudo apt install -y pocl-opencl-icd ocl-icd-opencl-dev
-
 . /home/firedrake/firedrake/bin/activate
 grep -v loopy requirements.txt > /tmp/myreq.txt
 
 # no need for these in the Firedrake tests
 sed -i "/boxtree/ d" /tmp/myreq.txt
 sed -i "/sumpy/ d" /tmp/myreq.txt
 sed -i "/pytential/ d" /tmp/myreq.txt
+sed -i "/pyopencl/ d" /tmp/myreq.txt
 
 # This shouldn't be necessary, but...
 # https://github.com/inducer/meshmode/pull/48#issuecomment-687519451
 pip install pybind11
-
 pip install pytest
-
 pip install -r /tmp/myreq.txt
+pip install pyopencl[pocl]
 
 # Context: https://github.com/OP2/PyOP2/pull/605
 python -m pip install --force-reinstall git+https://github.com/inducer/pytools.git

.github/workflows/ci.yml

@@ -150,14 +150,13 @@ jobs:
             run: |
                 # texlive is needed for sphinxcontrib-tikz
                 sudo apt-get update
-                sudo apt-get install -y --no-install-recommends texlive-pictures
+                sudo apt-get install -y --no-install-recommends texlive-pictures graphviz pdf2svg
 
                 curl -L -O https://gitlab.tiker.net/inducer/ci-support/raw/main/ci-support.sh
                 . ci-support.sh
 
                 export EXTRA_INSTALL="sphinxcontrib-tikz"
                 build_py_project_in_conda_env
-                conda install graphviz pdf2svg
 
                 # Work around
                 # intersphinx inventory 'https://firedrakeproject.org/objects.inv' not fetchable

meshmode/discretization/connection/direct.py

@@ -700,7 +700,7 @@ def group_pick_knl(is_surjective: bool):
                                 from_element_indices[iel],
                                 dof_pick_lists[dof_pick_list_indices[iel], idof]
                             ]
-                        { if_present })
+                        {if_present})
                 """,
                 [
                     lp.GlobalArg("ary", None,

meshmode/discretization/poly_element.py

@@ -41,6 +41,7 @@
         InterpolatoryElementGroupBase)
 
 import modepy as mp
+from modepy import Basis
 
 __doc__ = """
 Group types
@@ -518,7 +519,9 @@ def from_mesh_interp_matrix(self):
 
 class TensorProductElementGroupBase(PolynomialElementGroupBase,
         HypercubeElementGroupBase):
-    def __init__(self, mesh_el_group, order, index=None, *, basis, unit_nodes):
+    def __init__(self, mesh_el_group: _MeshTensorProductElementGroup,
+                 order: int, index=None, *, basis: Basis,
+                 unit_nodes: np.ndarray) -> None:
         """
         :arg basis: a :class:`modepy.TensorProductBasis`.
         :arg unit_nodes: unit nodes for the tensor product, obtained by
@@ -530,16 +533,43 @@ def __init__(self, mesh_el_group, order, index=None, *, basis, unit_nodes):
             raise ValueError("basis dimension does not match element group: "
                     f"expected {mesh_el_group.dim}, got {basis._dim}.")
 
+        if isinstance(basis, mp.TensorProductBasis):
+            for b in basis.bases:
+                if b._dim != 1:
+                    raise NotImplementedError(
+                        "All bases used to construct the tensor "
+                        "product must be of dimension 1. Support "
+                        "for higher-dimensional component bases "
+                        "does not yet exist.")
+
         if unit_nodes.shape[0] != mesh_el_group.dim:
             raise ValueError("unit node dimension does not match element group: "
-                    f"expected {mesh_el_group.dim}, got {unit_nodes.shape[0]}.")
+                    f"expected {self.mesh_el_group.dim}, "
+                    f"got {unit_nodes.shape[0]}.")
+
+        # NOTE there are cases where basis is a 1D `_SimplexONB` object. We wrap
+        # in a TensorProductBasis object if this is the case
+        if not isinstance(basis, mp.TensorProductBasis):
+            if basis._dim == 1 and unit_nodes.shape[0] == 1:
+                basis = mp.TensorProductBasis([basis])
+            else:
+                raise ValueError("`basis` is not a TensorProductBasis object, "
+                                 "and `basis` and `unit_nodes` are not both of "
+                                 "dimension 1. Found `basis` dim = {basis._dim}, "
+                                 "`unit_nodes` dim = {unit_nodes.shape[0]}.")
 
         self._basis = basis
+        self._bases_1d = basis.bases[0]
         self._nodes = unit_nodes
 
     def basis_obj(self):
         return self._basis
 
+    def bases_1d(self):
+        """Return 1D component bases used to construct the tensor product basis.
+        """
+        return self._bases_1d
+
     @memoize_method
     def quadrature_rule(self):
         basis_fcts = self._basis.functions
@@ -549,6 +579,11 @@ def quadrature_rule(self):
                          np.ones(len(basis_fcts)))
         return mp.Quadrature(nodes, weights, exact_to=self.order)
 
+    @property
+    @memoize_method
+    def unit_nodes_1d(self):
+        return self._nodes[0][:self.order + 1].reshape(1, self.order + 1)
+
     def discretization_key(self):
         # FIXME?
         # The unit_nodes numpy array isn't hashable, and comparisons would
@@ -604,9 +639,10 @@ class LegendreGaussLobattoTensorProductElementGroup(
     def __init__(self, mesh_el_group, order, index=None):
         from modepy.quadrature.jacobi_gauss import legendre_gauss_lobatto_nodes
         unit_nodes_1d = legendre_gauss_lobatto_nodes(order)
-        unit_nodes = mp.tensor_product_nodes([unit_nodes_1d] * mesh_el_group.dim)
+        unit_nodes = mp.tensor_product_nodes([unit_nodes_1d]*mesh_el_group.dim)
 
-        super().__init__(mesh_el_group, order, index=index, unit_nodes=unit_nodes)
+        super().__init__(mesh_el_group, order, index=index,
+                         unit_nodes=unit_nodes)
 
     def discretization_key(self):
         return (type(self), self.dim, self.order)
@@ -624,9 +660,10 @@ class EquidistantTensorProductElementGroup(LegendreTensorProductElementGroup):
     def __init__(self, mesh_el_group, order, index=None):
         from modepy.nodes import equidistant_nodes
         unit_nodes_1d = equidistant_nodes(1, order)[0]
-        unit_nodes = mp.tensor_product_nodes([unit_nodes_1d] * mesh_el_group.dim)
+        unit_nodes = mp.tensor_product_nodes([unit_nodes_1d]*mesh_el_group.dim)
 
-        super().__init__(mesh_el_group, order, index=index, unit_nodes=unit_nodes)
+        super().__init__(mesh_el_group, order, index=index,
+                         unit_nodes=unit_nodes)
 
     def discretization_key(self):
         return (type(self), self.dim, self.order)

meshmode/mesh/__init__.py

@@ -513,6 +513,11 @@ class TensorProductElementGroup(_ModepyElementGroup):
     r"""Inherits from :class:`MeshElementGroup`."""
     _modepy_shape_cls: ClassVar[Type[mp.Shape]] = mp.Hypercube
 
+    def is_affine(self):
+        # Tensor product mappings are generically bilinear.
+        # FIXME: Are affinely mapped ones a 'juicy' enough special case?
+        return False
+
 # }}}