Implement HDivTrace as a macroelement

FIAT/barycentric_interpolation.py

@@ -30,7 +30,7 @@ def barycentric_interpolation(nodes, wts, dmat, pts, order=0):
         sp_simplify = numpy.vectorize(simplify)
     else:
         sp_simplify = lambda x: x
-    phi = numpy.add.outer(-nodes, pts.flatten())
+    phi = numpy.add.outer(-nodes.flatten(), pts.flatten())
     with numpy.errstate(divide='ignore', invalid='ignore'):
         numpy.reciprocal(phi, out=phi)
         numpy.multiply(phi, wts[:, None], out=phi)
@@ -49,7 +49,7 @@ def barycentric_interpolation(nodes, wts, dmat, pts, order=0):
 def make_dmat(x):
     """Returns Lagrange differentiation matrix and barycentric weights
     associated with x[j]."""
-    dmat = numpy.add.outer(-x, x)
+    dmat = numpy.add.outer(-x.flatten(), x.flatten())
     numpy.fill_diagonal(dmat, 1.0)
     wts = numpy.prod(dmat, axis=0)
     numpy.reciprocal(wts, out=wts)
@@ -59,22 +59,38 @@ def make_dmat(x):
 
 
 class LagrangeLineExpansionSet(expansions.LineExpansionSet):
-    """Lagrange polynomial expansion set for given points the line."""
+    """Lagrange polynomial expansion set for given points on the line."""
     def __init__(self, ref_el, pts):
+        if ref_el.get_spatial_dimension() != 1:
+            raise Exception("Must have a line")
+        pts = numpy.asarray(pts)
         self.points = pts
-        self.x = numpy.array(pts, dtype="d").flatten()
+
         self.cell_node_map = expansions.compute_cell_point_map(ref_el, pts, unique=False)
         self.dmats = [None for _ in self.cell_node_map]
         self.weights = [None for _ in self.cell_node_map]
         self.nodes = [None for _ in self.cell_node_map]
+        self.affine_mappings = {}
         for cell, ibfs in self.cell_node_map.items():
-            self.nodes[cell] = self.x[ibfs]
-            self.dmats[cell], self.weights[cell] = make_dmat(self.nodes[cell])
+            x = pts[ibfs]
+            if ref_el.is_trace():
+                verts = ref_el.get_vertices_of_subcomplex(ref_el.topology[1][cell])
+                A, = numpy.diff(verts, axis=0)
+                A /= numpy.linalg.norm(A)
+                b = -numpy.dot(numpy.sum(verts, axis=0)/2, A.T)
+                self.affine_mappings[cell] = (A, b)
+                x = numpy.add(numpy.dot(x, A.T), b)
+            self.nodes[cell] = x
+            self.dmats[cell], self.weights[cell] = make_dmat(x)
 
         self.degree = max(len(wts) for wts in self.weights)-1
         self.recurrence_order = self.degree + 1
-        super().__init__(ref_el)
-        self.continuity = None if len(self.x) == sum(len(xk) for xk in self.nodes) else "C0"
+        self.ref_el = ref_el
+        self.variant = None
+        self.scale = 1
+        self.continuity = None if len(pts) == sum(len(xk) for xk in self.nodes) else "C0"
+        self._dmats_cache = {}
+        self._cell_node_map_cache = {}
 
     def get_num_members(self, n):
         return len(self.points)
@@ -89,6 +105,12 @@ def get_dmats(self, degree, cell=0):
         return [self.dmats[cell].T]
 
     def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None):
+        try:
+            A, b = self.affine_mappings[cell]
+            ref_pts = numpy.add(numpy.dot(pts.reshape(-1, A.shape[-1]), A.T), b)
+            pts = ref_pts.reshape(*pts.shape[:-1], -1)
+        except KeyError:
+            pass
         return barycentric_interpolation(self.nodes[cell], self.weights[cell], self.dmats[cell], pts, order=order)
 
 

FIAT/bubble.py

@@ -14,8 +14,8 @@
 class CodimBubble(RestrictedElement):
     """Bubbles of a certain codimension."""
 
-    def __init__(self, ref_el, degree, codim):
-        element = Lagrange(ref_el, degree)
+    def __init__(self, ref_el, degree, codim, variant=None):
+        element = Lagrange(ref_el, degree, variant=variant)
 
         cell_dim = ref_el.get_dimension()
         assert cell_dim == max(element.entity_dofs().keys())
@@ -29,12 +29,12 @@ def __init__(self, ref_el, degree, codim):
 class Bubble(CodimBubble):
     """The bubble finite element: the dofs of the Lagrange FE in the interior of the cell"""
 
-    def __init__(self, ref_el, degree):
-        super().__init__(ref_el, degree, codim=0)
+    def __init__(self, ref_el, degree, variant=None):
+        super().__init__(ref_el, degree, codim=0, variant=variant)
 
 
 class FacetBubble(CodimBubble):
     """The facet bubble finite element: the dofs of the Lagrange FE in the interior of the facets"""
 
-    def __init__(self, ref_el, degree):
-        super().__init__(ref_el, degree, codim=1)
+    def __init__(self, ref_el, degree, variant=None):
+        super().__init__(ref_el, degree, codim=1, variant=variant)

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():

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:
@@ -354,6 +354,9 @@ def _tabulate(self, n, pts, order=0):
         if not self.ref_el.is_macrocell():
             return phis[0]
 
+        if self.ref_el.is_trace():
+            phis = trace_tabulation(self.ref_el, cell_point_map, order, pts, phis)
+
         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 +365,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 +381,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], Exception):
+                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 +504,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 +597,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)
+    DX = X[1:] - X[:1]
+    DY = Y[1:] - Y[:1]
+    if DX.shape[0] == DX.shape[1]:
+        AT = numpy.linalg.solve(DX, DY)
+    else:
+        AT, *_ = numpy.linalg.lstsq(DX, DY)
+    b = Y[0] - numpy.dot(X[0], AT)
+    return AT.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)
@@ -679,7 +702,10 @@ def compute_cell_point_map(ref_el, pts, unique=True, tol=1E-12):
         return {cell: Ellipsis for cell in sorted(top[sd])}
 
     # The distance to the nearest cell is equal to the distance to the parent cell
-    best = ref_el.get_parent().distance_to_point_l1(pts, rescale=True)
+    parent = ref_el
+    while parent.get_parent() is not None:
+        parent = parent.get_parent()
+    best = parent.distance_to_point_l1(pts, rescale=True)
     tol = best + tol
 
     cell_point_map = {}
@@ -735,3 +761,29 @@ def compute_partition_of_unity(ref_el, pt, unique=True, tol=1E-12):
         mult = sum(masks)
         masks = [m / mult for m in masks]
     return masks
+
+
+def trace_tabulation(ref_el, cell_point_map, order, pts, phis):
+    """Lift trace tabulations into the cells and raise TraceError on invalid tabulations."""
+    from FIAT.polynomial_set import mis
+    from FIAT.hdiv_trace import TraceError
+    parent = ref_el.get_parent()
+    tdim = ref_el.get_spatial_dimension()
+    gdim = parent.get_spatial_dimension()
+    for cell in phis:
+        # Lift 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 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] = {alpha: TraceError(msg)
+                          for i in range(order+1)
+                          for alpha in mis(gdim, i)}
+    return phis