Implement HDivTrace as a macroelement

Author: pbrubeck

FIAT/bubble.py

@@ -7,7 +7,6 @@
 # SPDX-License-Identifier:    LGPL-3.0-or-later
 
 from FIAT.lagrange import Lagrange
-from FIAT.hierarchical import IntegratedLegendre
 from FIAT.restricted import RestrictedElement
 from itertools import chain
 
@@ -16,11 +15,7 @@ class CodimBubble(RestrictedElement):
     """Bubbles of a certain codimension."""
 
     def __init__(self, ref_el, degree, codim, variant=None):
-        if variant and variant.startswith("integral"):
-            CG = IntegratedLegendre
-        else:
-            CG = Lagrange
-        element = CG(ref_el, degree, variant=variant)
+        element = Lagrange(ref_el, degree, variant=variant)
 
         cell_dim = ref_el.get_dimension()
         assert cell_dim == max(element.entity_dofs().keys())

FIAT/discontinuous_lagrange.py

@@ -14,6 +14,7 @@
 from FIAT.barycentric_interpolation import LagrangePolynomialSet, get_lagrange_points
 from FIAT.polynomial_set import mis
 from FIAT.check_format_variant import parse_lagrange_variant
+from FIAT.hierarchical import Legendre
 
 
 def make_entity_permutations(dim, npoints):
@@ -215,6 +216,8 @@ class DiscontinuousLagrange(finite_element.CiarletElement):
                               macroelement for Scott-Vogelius.
     """
     def __new__(cls, ref_el, degree, variant="equispaced"):
+        if variant and variant.startswith("integral"):
+            return Legendre(ref_el, degree, variant=variant)
         if degree == 0:
             splitting, _ = parse_lagrange_variant(variant, discontinuous=True)
             if splitting is None and not ref_el.is_macrocell():
@@ -230,7 +233,7 @@ def __init__(self, ref_el, degree, variant="equispaced"):
             dual = BrokenLagrangeDualSet(ref_el, degree, point_variant=point_variant)
         else:
             dual = DiscontinuousLagrangeDualSet(ref_el, degree, point_variant=point_variant)
-        if ref_el.shape == LINE:
+        if ref_el.shape == LINE and not ref_el.is_trace():
             # In 1D we can use the primal basis as the expansion set,
             # avoiding any round-off coming from a basis transformation
             points = get_lagrange_points(dual)

FIAT/dual_set.py

@@ -277,15 +277,15 @@ def merge_entities(nodes, ref_el, entity_ids, entity_permutations):
     parent_cell = ref_el.get_parent()
     if parent_cell is None:
         return nodes, ref_el, entity_ids, entity_permutations
-    parent_ids = {}
+    parent_top = parent_cell.get_topology()
+    parent_ids = {dim: {entity: [] for entity in parent_top[dim]} for dim in parent_top}
     parent_permutations = None
     parent_to_children = ref_el.get_parent_to_children()
 
-    if all(isinstance(node, functional.PointEvaluation) for node in nodes):
+    if all(isinstance(node, functional.PointEvaluation) for node in nodes) and not ref_el.is_trace():
         # Merge Lagrange dual with lexicographical reordering
         parent_nodes = []
         for dim in sorted(parent_to_children):
-            parent_ids[dim] = {}
             for entity in sorted(parent_to_children[dim]):
                 cur = len(parent_nodes)
                 for child_dim, child_entity in parent_to_children[dim][entity]:
@@ -296,9 +296,7 @@ def merge_entities(nodes, ref_el, entity_ids, entity_permutations):
         # Merge everything else with the same node ordering
         parent_nodes = nodes
         for dim in sorted(parent_to_children):
-            parent_ids[dim] = {}
             for entity in sorted(parent_to_children[dim]):
-                parent_ids[dim][entity] = []
                 for child_dim, child_entity in parent_to_children[dim][entity]:
                     parent_ids[dim][entity].extend(entity_ids[child_dim][child_entity])
 

FIAT/expansions.py

@@ -268,7 +268,7 @@ def __init__(self, ref_el, scale=None, variant=None):
         base_ref_el = reference_element.default_simplex(sd)
         base_verts = base_ref_el.get_vertices()
 
-        self.affine_mappings = [reference_element.make_affine_mapping(
+        self.affine_mappings = [get_affine_mapping(
                                 ref_el.get_vertices_of_subcomplex(top[sd][cell]),
                                 base_verts) for cell in top[sd]]
         if scale is None:
@@ -345,6 +345,8 @@ def distance(alpha, beta):
 
     def _tabulate(self, n, pts, order=0):
         """A version of tabulate() that also works for a single point."""
+        from FIAT.hdiv_trace import TraceError
+        from FIAT.polynomial_set import mis
         pts = numpy.asarray(pts)
         unique = self.continuity is not None and order == 0
         cell_point_map = compute_cell_point_map(self.ref_el, pts, unique=unique)
@@ -354,6 +356,28 @@ def _tabulate(self, n, pts, order=0):
         if not self.ref_el.is_macrocell():
             return phis[0]
 
+        if self.ref_el.is_trace():
+            parent = self.ref_el.get_parent()
+            tdim = self.ref_el.get_spatial_dimension()
+            gdim = parent.get_spatial_dimension()
+            for cell in phis:
+                # Promote facet keys to cell keys
+                phi = phis[cell][(0,) * tdim]
+                # Raise TraceError on gradient tabulations
+                msg = "Gradients on trace elements are not well-defined."
+                phis[cell] = {alpha: phi if sum(alpha) == 0 else TraceError(msg)
+                              for i in range(order+1)
+                              for alpha in mis(gdim, i)}
+            if tdim == 0 and len(phis) == 0:
+                # Hack for TensorProduct HDivTrace: do not raise TraceError on the interval
+                for cell in parent.topology[gdim]:
+                    phis[cell] = {(0,)*gdim: numpy.zeros((1, 1))}
+            elif sum(len(cell_point_map[cell]) for cell in cell_point_map) < len(pts):
+                # Raise TraceError when interior points fail to be binned on facets
+                for cell in parent.topology[gdim]:
+                    msg = "The HDivTrace element can only be tabulated on facets."
+                    phis[cell] = {(0,)*gdim: TraceError(msg)}
+
         if pts.dtype == object:
             # If binning is undefined, scale by the characteristic function of each subcell
             Xi = compute_partition_of_unity(self.ref_el, pts, unique=unique)
@@ -362,13 +386,14 @@ def _tabulate(self, n, pts, order=0):
                     phi[alpha] *= Xi[cell]
         elif not unique:
             # If binning is not unique, divide by the multiplicity of each point
-            mult = numpy.zeros(pts.shape[:-1])
+            mult = numpy.zeros(pts.shape[:-1], dtype=int)
             for cell, ipts in cell_point_map.items():
                 mult[ipts] += 1
-            for cell, ipts in cell_point_map.items():
-                phi = phis[cell]
-                for alpha in phi:
-                    phi[alpha] /= mult[None, ipts]
+            if (mult != 1).any():
+                for cell, ipts in cell_point_map.items():
+                    phi = phis[cell]
+                    for alpha in phi:
+                        phi[alpha] /= mult[None, ipts]
 
         # Insert subcell tabulations into the corresponding submatrices
         idx = lambda *args: args if args[1] is Ellipsis else numpy.ix_(*args)
@@ -377,6 +402,9 @@ def _tabulate(self, n, pts, order=0):
         result = {}
         base_phi = tuple(phis.values())[0]
         for alpha in base_phi:
+            if isinstance(base_phi[alpha], TraceError):
+                result[alpha] = base_phi[alpha]
+                continue
             dtype = base_phi[alpha].dtype
             result[alpha] = numpy.zeros((num_phis, *pts.shape[:-1]), dtype=dtype)
             for cell in cell_point_map:
@@ -497,8 +525,7 @@ def get_dmats(self, degree, cell=0):
     def tabulate(self, n, pts):
         if len(pts) == 0:
             return numpy.array([])
-        sd = self.ref_el.get_spatial_dimension()
-        return self._tabulate(n, pts)[(0,) * sd]
+        return tuple(self._tabulate(n, pts).values())[0]
 
     def tabulate_derivatives(self, n, pts):
         from FIAT.polynomial_set import mis
@@ -591,13 +618,30 @@ def __init__(self, ref_el, **kwargs):
         super().__init__(ref_el, **kwargs)
 
 
+def get_affine_mapping(xs, ys):
+    """Constructs (A,b) such that x --> A * x + b is the affine
+    mapping from the simplex defined by xs to the simplex defined by ys.
+    Uses least-squares when the simplex dimension does not match the spatial
+    dimension.
+    """
+    X = numpy.asarray(xs)
+    Y = numpy.asarray(ys)
+    A = X[1:] - X[:1]
+    B = Y[1:] - Y[:1]
+    if A.shape[0] == A.shape[1]:
+        C = numpy.linalg.solve(A, B)
+    else:
+        C, *_ = numpy.linalg.lstsq(A, B)
+    b = Y[0] - numpy.dot(X[0], C)
+    return C.T, b
+
+
 def polynomial_dimension(ref_el, n, continuity=None):
     """Returns the dimension of the space of polynomials of degree no
     greater than n on the reference complex."""
     if ref_el.get_shape() == reference_element.POINT:
         if n > 0:
             raise ValueError("Only degree zero polynomials supported on point elements.")
-        return 1
     top = ref_el.get_topology()
     if continuity == "C0":
         space_dimension = sum(math.comb(n - 1, dim) * len(top[dim]) for dim in top)