Boundary Quadrature element

FIAT/polynomial_set.py

@@ -69,7 +69,7 @@ def tabulate_new(self, pts):
     def tabulate(self, pts, jet_order=0):
         """Returns the values of the polynomial set."""
         base_vals = self.expansion_set._tabulate(self.embedded_degree, pts, order=jet_order)
-        result = {alpha: numpy.dot(self.coeffs, base_vals[alpha]) for alpha in base_vals}
+        result = {alpha: numpy.tensordot(self.coeffs, base_vals[alpha], (-1, 0)) for alpha in base_vals}
         return result
 
     def get_expansion_set(self):

finat/element_factory.py

@@ -149,13 +149,14 @@ def convert(element, **kwargs):
 @convert.register(finat.ufl.FiniteElement)
 def convert_finiteelement(element, **kwargs):
     cell = as_fiat_cell(element.cell)
-    if element.family() == "Quadrature":
+    if element.family() in {"Quadrature", "Boundary Quadrature"}:
         degree = element.degree()
-        scheme = element.quadrature_scheme()
+        scheme = element.quadrature_scheme() or "default"
         if degree is None or scheme is None:
             raise ValueError("Quadrature scheme and degree must be specified!")
 
-        return finat.make_quadrature_element(cell, degree, scheme), set()
+        codim = 1 if element.family() == "Boundary Quadrature" else 0
+        return finat.make_quadrature_element(cell, degree, scheme, codim), set()
     lmbda = supported_elements[element.family()]
     if element.family() == "Real" and element.cell.cellname() in {"quadrilateral", "hexahedron"}:
         lmbda = None

finat/fiat_elements.py

@@ -116,9 +116,8 @@ def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
                 continue
 
             derivative = sum(alpha)
-            table_roll = fiat_table.reshape(
-                space_dimension, value_size, len(ps.points)
-            ).transpose(1, 2, 0)
+            shp = (space_dimension, value_size, *ps.points.shape[:-1])
+            table_roll = np.moveaxis(fiat_table.reshape(shp), 0, -1)
 
             exprs = []
             for table in table_roll:

finat/point_set.py

@@ -33,8 +33,7 @@ def points(self):
     @property
     def dimension(self):
         """Point dimension."""
-        _, dim = self.points.shape
-        return dim
+        return self.points.shape[-1]
 
     @property
     @abc.abstractmethod
@@ -129,8 +128,7 @@ def points(self):
 
     @cached_property
     def indices(self):
-        N, _ = self._points_expr.shape
-        return (gem.Index(extent=N),)
+        return tuple(gem.Index(extent=N) for N in self._points_expr.shape[:-1])
 
     @cached_property
     def expression(self):
@@ -159,7 +157,7 @@ def points(self):
 
     @cached_property
     def indices(self):
-        return (gem.Index(extent=len(self.points)),)
+        return tuple(gem.Index(extent=N) for N in self.points.shape[:-1])
 
     @cached_property
     def expression(self):
@@ -224,3 +222,40 @@ def almost_equal(self, other, tolerance=1e-12):
             len(self.factors) == len(other.factors) and \
             all(s.almost_equal(o, tolerance=tolerance)
                 for s, o in zip(self.factors, other.factors))
+
+
+class MappedPointSet(AbstractPointSet):
+
+    def __init__(self, cell, ps):
+        self.cell = cell
+        self.ps = ps
+
+    @cached_property
+    def entities(self):
+        to_int = lambda x: sum(x) if isinstance(x, tuple) else x
+        top = self.cell.topology
+        return [(dim, entity)
+                for dim in sorted(top)
+                for entity in sorted(top[dim])
+                if to_int(dim) == self.ps.dimension]
+
+    @cached_property
+    def points(self):
+        ref_pts = self.ps.points
+        pts = [self.cell.get_entity_transform(dim, entity)(ref_pts)
+               for dim, entity in self.entities]
+        return numpy.concatenate(pts)
+
+    @cached_property
+    def indices(self):
+        return (gem.Index(extent=len(self.entities)), *self.ps.indices)
+
+    @cached_property
+    def expression(self):
+        raise NotImplementedError("Should not use MappedPointSet like this")
+
+    def almost_equal(self, other, tolerance=1e-12):
+        """Approximate numerical equality of point sets"""
+        return type(self) is type(other) and \
+            self.cell == other.cell and \
+            self.ps.almost_equal(other.ps, tolerance=tolerance)

finat/quadrature.py

@@ -1,6 +1,6 @@
 import hashlib
 from abc import ABCMeta, abstractmethod
-from functools import cached_property, reduce
+from functools import cached_property
 
 import gem
 import numpy
@@ -163,4 +163,4 @@ def point_set(self):
 
     @cached_property
     def weight_expression(self):
-        return reduce(gem.Product, (q.weight_expression for q in self.factors))
+        return gem.Product(*(q.weight_expression for q in self.factors))

finat/quadrature_element.py

@@ -1,5 +1,4 @@
-from finat.point_set import UnknownPointSet
-from functools import reduce
+from finat.point_set import UnknownPointSet, MappedPointSet
 
 import numpy
 
@@ -13,7 +12,7 @@
 from finat.quadrature import make_quadrature, AbstractQuadratureRule
 
 
-def make_quadrature_element(fiat_ref_cell, degree, scheme="default"):
+def make_quadrature_element(fiat_ref_cell, degree, scheme="default", codim=0):
     """Construct a :class:`QuadratureElement` from a given a reference
     element, degree and scheme.
 
@@ -23,9 +22,16 @@ def make_quadrature_element(fiat_ref_cell, degree, scheme="default"):
         integrate exactly.
     :param scheme: The quadrature scheme to use - e.g. "default",
         "canonical" or "KMV".
+    :param codim: The codimension of the quadrature scheme.
     :returns: The appropriate :class:`QuadratureElement`
     """
-    rule = make_quadrature(fiat_ref_cell, degree, scheme=scheme)
+    if codim:
+        sd = fiat_ref_cell.get_spatial_dimension()
+        rule_ref_cell = fiat_ref_cell.construct_subelement(sd - codim)
+    else:
+        rule_ref_cell = fiat_ref_cell
+
+    rule = make_quadrature(rule_ref_cell, degree, scheme=scheme)
     return QuadratureElement(fiat_ref_cell, rule)
 
 
@@ -42,8 +48,6 @@ def __init__(self, fiat_ref_cell, rule):
         self.cell = fiat_ref_cell
         if not isinstance(rule, AbstractQuadratureRule):
             raise TypeError("rule is not an AbstractQuadratureRule")
-        if fiat_ref_cell.get_spatial_dimension() != rule.point_set.dimension:
-            raise ValueError("Cell dimension does not match rule's point set dimension")
         self._rule = rule
 
     @cached_property
@@ -64,10 +68,18 @@ def formdegree(self):
 
     @cached_property
     def _entity_dofs(self):
-        # Inspired by ffc/quadratureelement.py
+        top = self.cell.get_topology()
         entity_dofs = {dim: {entity: [] for entity in entities}
-                       for dim, entities in self.cell.get_topology().items()}
-        entity_dofs[self.cell.get_dimension()] = {0: list(range(self.space_dimension()))}
+                       for dim, entities in top.items()}
+        ps = self._rule.point_set
+        num_pts = len(ps.points)
+        to_int = lambda x: sum(x) if isinstance(x, tuple) else x
+        cur = 0
+        for dim in sorted(top):
+            if to_int(dim) == ps.dimension:
+                for entity in sorted(top[dim]):
+                    entity_dofs[dim][entity].extend(range(cur, cur + num_pts))
+                    cur += num_pts
         return entity_dofs
 
     def entity_dofs(self):
@@ -76,9 +88,15 @@ def entity_dofs(self):
     def space_dimension(self):
         return numpy.prod(self.index_shape, dtype=int)
 
+    @cached_property
+    def _point_set(self):
+        ps = self._rule.point_set
+        sd = self.cell.get_spatial_dimension()
+        return ps if ps.dimension == sd else MappedPointSet(self.cell, ps)
+
     @property
     def index_shape(self):
-        ps = self._rule.point_set
+        ps = self._point_set
         return tuple(index.extent for index in ps.indices)
 
     @property
@@ -87,7 +105,7 @@ def value_shape(self):
 
     @cached_property
     def fiat_equivalent(self):
-        ps = self._rule.point_set
+        ps = self._point_set
         if isinstance(ps, UnknownPointSet):
             raise ValueError("A quadrature element with rule with runtime points has no fiat equivalent!")
         weights = getattr(self._rule, 'weights', None)
@@ -107,8 +125,13 @@ def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
         :param ps: the point set object.
         :param entity: the cell entity on which to tabulate.
         '''
-        if entity is not None and entity != (self.cell.get_dimension(), 0):
-            raise ValueError('QuadratureElement does not "tabulate" on subentities.')
+        rule_dim = self._rule.point_set.dimension
+        if entity is None:
+            entity = (rule_dim, 0)
+        entity_dim, entity_id = entity
+        if entity_dim != rule_dim:
+            raise ValueError(f"Cannot tabulate QuadratureElement of dimension {rule_dim}"
+                             f" on subentities of dimension {entity_dim}.")
 
         if order:
             raise ValueError("Derivatives are not defined on a QuadratureElement.")
@@ -118,23 +141,24 @@ def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
 
         # Return an outer product of identity matrices
         multiindex = self.get_indices()
-        product = reduce(gem.Product, [gem.Delta(q, r)
-                                       for q, r in zip(ps.indices, multiindex)])
+        fid = ps.indices
+        if len(multiindex) > len(fid):
+            fid = (entity_id, *fid)
+        product = gem.Delta(fid, multiindex)
 
-        dim = self.cell.get_spatial_dimension()
-        return {(0,) * dim: gem.ComponentTensor(product, multiindex)}
+        sd = self.cell.get_spatial_dimension()
+        return {(0,) * sd: gem.ComponentTensor(product, multiindex)}
 
     def point_evaluation(self, order, refcoords, entity=None):
         raise NotImplementedError("QuadratureElement cannot do point evaluation!")
 
     @property
     def dual_basis(self):
-        ps = self._rule.point_set
+        ps = self._point_set
         multiindex = self.get_indices()
         # Evaluation matrix is just an outer product of identity
         # matrices, evaluation points are just the quadrature points.
-        Q = reduce(gem.Product, (gem.Delta(q, r)
-                                 for q, r in zip(ps.indices, multiindex)))
+        Q = gem.Delta(ps.indices, multiindex)
         Q = gem.ComponentTensor(Q, multiindex)
         return Q, ps
 

finat/sympy2gem.py

@@ -27,13 +27,13 @@ def sympy2gem_expr(node, self):
 @sympy2gem.register(sympy.Add)
 @sympy2gem.register(symengine.Add)
 def sympy2gem_add(node, self):
-    return reduce(gem.Sum, map(self, node.args))
+    return gem.Sum(*map(self, node.args))
 
 
 @sympy2gem.register(sympy.Mul)
 @sympy2gem.register(symengine.Mul)
 def sympy2gem_mul(node, self):
-    return reduce(gem.Product, map(self, node.args))
+    return gem.Product(*map(self, node.args))
 
 
 @sympy2gem.register(sympy.Pow)

finat/tensor_product.py

@@ -1,4 +1,3 @@
-from functools import reduce
 from itertools import chain, product
 from operator import methodcaller
 
@@ -32,11 +31,11 @@ def __init__(self, factors):
 
     @cached_property
     def cell(self):
-        return TensorProductCell(*[fe.cell for fe in self.factors])
+        return TensorProductCell(*(fe.cell for fe in self.factors))
 
     @cached_property
     def complex(self):
-        return TensorProductCell(*[fe.complex for fe in self.factors])
+        return TensorProductCell(*(fe.complex for fe in self.factors))
 
     @property
     def degree(self):
@@ -69,7 +68,7 @@ def space_dimension(self):
 
     @property
     def index_shape(self):
-        return tuple(chain(*[fe.index_shape for fe in self.factors]))
+        return tuple(chain.from_iterable(fe.index_shape for fe in self.factors))
 
     @property
     def value_shape(self):
@@ -117,24 +116,20 @@ def _merge_evaluations(self, factor_results):
         # multiindex describing the value shape of the subelement.
         zetas = [fe.get_value_indices() for fe in self.factors]
 
+        multiindex = tuple(chain(*alphas, *zetas))
         result = {}
         for derivative in range(order + 1):
             for Delta in mis(dimension, derivative):
                 # Split the multiindex for the subelements
                 deltas = [Delta[s] for s in dim_slices]
-                # GEM scalars (can have free indices) for collecting
-                # the contributions from the subelements.
-                scalars = []
-                for fr, delta, alpha, zeta in zip(factor_results, deltas, alphas, zetas):
-                    # Turn basis shape to free indices, select the
-                    # right derivative entry, and collect the result.
-                    scalars.append(gem.Indexed(fr[delta], alpha + zeta))
-                # Multiply the values from the subelements and wrap up
-                # non-point indices into shape.
-                result[Delta] = gem.ComponentTensor(
-                    reduce(gem.Product, scalars),
-                    tuple(chain(*(alphas + zetas)))
-                )
+                # Multiply the values from the subelements
+                # Turn basis shape to free indices, select the
+                # right derivative entry.
+                scalar = gem.Product(*(gem.Indexed(fr[delta], alpha + zeta)
+                                       for fr, delta, alpha, zeta
+                                       in zip(factor_results, deltas, alphas, zetas)))
+                # Wrap up non-point indices into shape.
+                result[Delta] = gem.ComponentTensor(scalar, multiindex)
         return result
 
     def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
@@ -179,12 +174,11 @@ def dual_basis(self):
         # Naming as _merge_evaluations above
         alphas = [factor.get_indices() for factor in self.factors]
         zetas = [factor.get_value_indices() for factor in self.factors]
-        # Index the factors by so that we can reshape into index-shape
-        # followed by value-shape
-        qis = [q[alpha + zeta] for q, alpha, zeta in zip(qs, alphas, zetas)]
         Q = gem.ComponentTensor(
-            reduce(gem.Product, qis),
-            tuple(chain(*(alphas + zetas)))
+            # Index the factors by so that we can reshape into index-shape
+            # followed by value-shape
+            gem.Product(*(q[alpha + zeta] for q, alpha, zeta in zip(qs, alphas, zetas))),
+            tuple(chain(*alphas, *zetas))
         )
         return Q, ps
 

finat/tensorfiniteelement.py

@@ -1,4 +1,3 @@
-from functools import reduce
 from itertools import chain
 
 import numpy
@@ -134,8 +133,7 @@ def _tensorise(self, scalar_evaluation):
         tensor_vi = tuple(gem.Index(extent=d) for d in self._shape)
 
         # Couple new basis function and value indices
-        deltas = reduce(gem.Product, (gem.Delta(j, k)
-                                      for j, k in zip(tensor_i, tensor_vi)))
+        deltas = gem.Delta(tensor_i, tensor_vi)
 
         if self._transpose:
             index_ordering = tensor_i + scalar_i + tensor_vi + scalar_vi
@@ -163,8 +161,7 @@ def dual_basis(self):
         tensor_i = tuple(gem.Index(extent=d) for d in self._shape)
         tensor_vi = tuple(gem.Index(extent=d) for d in self._shape)
         # Couple new basis function and value indices
-        deltas = reduce(gem.Product, (gem.Delta(j, k)
-                                      for j, k in zip(tensor_i, tensor_vi)))
+        deltas = gem.Delta(tensor_i, tensor_vi)
         if self._transpose:
             index_ordering = tensor_i + scalar_i + tensor_vi + scalar_vi
         else: