Integrate hypercube changes to fuse

FIAT/discontinuous_pc.py

@@ -7,19 +7,13 @@
 # Modified by David A. Ham (david.ham@imperial.ac.uk), 2018
 
 from FIAT import finite_element, polynomial_set, dual_set, functional
-from FIAT.reference_element import (Point,
-                                    DefaultLine,
-                                    UFCInterval,
-                                    UFCQuadrilateral,
-                                    UFCHexahedron,
-                                    UFCTriangle,
-                                    UFCTetrahedron,
-                                    make_affine_mapping,
+from FIAT.reference_element import (make_affine_mapping,
                                     flatten_reference_cube,
                                     cell_to_simplex)
 from FIAT.P0 import P0Dual
 import numpy as np
 
+
 class DPC0(finite_element.CiarletElement):
     def __init__(self, ref_el):
         flat_el = flatten_reference_cube(ref_el)

FIAT/dual_set.py

@@ -35,6 +35,9 @@ def __init__(self, nodes, ref_el, entity_ids, entity_permutations=None):
     def __iter__(self):
         return iter(self.nodes)
 
+    def __len__(self):
+        return len(self.nodes)
+
     def get_nodes(self):
         return self.nodes
 

FIAT/expansions.py

@@ -436,26 +436,38 @@ def tabulate_jumps(self, n, points, order=0):
         """
 
         from FIAT.polynomial_set import mis
+        sd = self.ref_el.get_spatial_dimension()
         num_members = self.get_num_members(n)
         cell_node_map = self.get_cell_node_map(n)
+        cell_point_map = compute_cell_point_map(self.ref_el, points, unique=False)
+
+        num_jumps = 0
+        facet_point_map = {}
+        for facet in self.ref_el.get_interior_facets(sd-1):
+            try:
+                cells = self.ref_el.connectivity[(sd-1, sd)][facet]
+                ipts = list(set.intersection(*(set(cell_point_map[c]) for c in cells)))
+                if ipts != ():
+                    facet_point_map[facet] = ipts
+                    num_jumps += len(ipts)
+            except KeyError:
+                pass
 
-        top = self.ref_el.get_topology()
-        sd = self.ref_el.get_spatial_dimension()
-        interior_facets = self.ref_el.get_interior_facets(sd-1)
         derivs = {cell: self._tabulate_on_cell(n, points, order=order, cell=cell)
-                  for cell in top[sd]}
+                  for cell in cell_point_map}
+
         jumps = {}
         for r in range(order+1):
             cur = 0
             alphas = mis(sd, r)
-            jumps[r] = numpy.zeros((num_members, len(alphas)*len(interior_facets), len(points)))
-            for facet in interior_facets:
-                facet_verts = set(top[sd-1][facet])
-                c0, c1 = tuple(k for k in top[sd] if facet_verts < set(top[sd][k]))
+            jumps[r] = numpy.zeros((num_members, len(alphas) * num_jumps))
+            for facet, ipts in facet_point_map.items():
+                c0, c1 = self.ref_el.connectivity[(sd-1, sd)][facet]
                 for alpha in alphas:
-                    jumps[r][cell_node_map[c1], cur] += derivs[c1][alpha]
-                    jumps[r][cell_node_map[c0], cur] -= derivs[c0][alpha]
-                    cur = cur + 1
+                    ijump = range(cur, cur + len(ipts))
+                    jumps[r][numpy.ix_(cell_node_map[c1], ijump)] += derivs[c1][alpha][:, ipts]
+                    jumps[r][numpy.ix_(cell_node_map[c0], ijump)] -= derivs[c0][alpha][:, ipts]
+                    cur += len(ipts)
         return jumps
 
     def get_dmats(self, degree, cell=0):

FIAT/finite_element.py

@@ -134,6 +134,9 @@ def __init__(self, poly_set, dual, order, formdegree=None, mapping="affine", ref
         ref_complex = ref_complex or poly_set.get_reference_element()
         super().__init__(ref_el, dual, order, formdegree, mapping, ref_complex)
 
+        if len(poly_set) != len(dual):
+            raise ValueError(f"Dimension of function space is {len(poly_set)}, but got {len(dual)} nodes.")
+
         # build generalized Vandermonde matrix
         old_coeffs = poly_set.get_coeffs()
         dualmat = dual.to_riesz(poly_set)

FIAT/macro.py

@@ -412,21 +412,35 @@ class CkPolynomialSet(polynomial_set.PolynomialSet):
 
     :arg ref_el: The simplicial complex.
     :arg degree: The polynomial degree.
-    :kwarg order: The order of continuity across subcells.
+    :kwarg order: The order of continuity across facets or a dict of dicts
+        mapping dimension and entity_id to the order of continuity at each entity.
     :kwarg vorder: The order of super-smoothness at interior vertices.
     :kwarg shape: The value shape.
     :kwarg variant: The variant for the underlying ExpansionSet.
     :kwarg scale: The scale for the underlying ExpansionSet.
     """
-    def __init__(self, ref_el, degree, order=1, vorder=0, shape=(), **kwargs):
+    def __init__(self, ref_el, degree, order=1, vorder=None, shape=(), **kwargs):
         from FIAT.quadrature_schemes import create_quadrature
+
+        if not isinstance(order, (int, dict)):
+            raise TypeError(f"'order' must be either an int or dict, not {type(order).__name__}")
+
+        sd = ref_el.get_spatial_dimension()
+        if isinstance(order, int):
+            order = {sd-1: dict.fromkeys(ref_el.get_interior_facets(sd-1), order)}
+        if vorder is not None:
+            order[0] = dict.fromkeys(ref_el.get_interior_facets(0), vorder)
+        elif 0 not in order:
+            order[0] = {}
+
+        if not all(k in {0, sd-1} for k in order):
+            raise NotImplementedError("Only face or vertex constraints have been implemented.")
+
         expansion_set = expansions.ExpansionSet(ref_el, **kwargs)
         k = 1 if expansion_set.continuity == "C0" else 0
-        num_members = expansion_set.get_num_members(degree)
 
-        sd = ref_el.get_spatial_dimension()
+        # Impose C^forder continuity across interior facets
         facet_el = ref_el.construct_subelement(sd-1)
-
         phi_deg = 0 if sd == 1 else degree - k
         phi = polynomial_set.ONPolynomialSet(facet_el, phi_deg)
         Q = create_quadrature(facet_el, 2 * phi_deg)
@@ -435,30 +449,33 @@ def __init__(self, ref_el, degree, order=1, vorder=0, shape=(), **kwargs):
         weights = numpy.multiply(phi_at_qpts, qwts)
 
         rows = []
-        for facet in ref_el.get_interior_facets(sd-1):
-            jumps = expansion_set.tabulate_normal_jumps(degree, qpts, facet, order=order)
-            for r in range(k, order+1):
+        for facet in order[sd-1]:
+            forder = order[sd-1][facet]
+            jumps = expansion_set.tabulate_normal_jumps(degree, qpts, facet, order=forder)
+            for r in range(k, forder+1):
                 num_wt = 1 if sd == 1 else expansions.polynomial_dimension(facet_el, degree-r)
-                rows.append(numpy.tensordot(weights[:num_wt], jumps[r], axes=(-1, -1)).reshape(-1, num_members))
+                rows.append(numpy.tensordot(weights[:num_wt], jumps[r], axes=(-1, -1)))
+
+        # Impose C^vorder super-smoothness at interior vertices
+        # C^forder automatically gives C^{forder+dim-1} at the interior vertex
+        verts = numpy.asarray(ref_el.get_vertices())
+        for vorder in set(order[0].values()):
+            vids = [i for i in order[0] if order[0][i] == vorder]
+            facets = chain.from_iterable(ref_el.connectivity[(0, sd-1)][v] for v in vids)
+            forder = min(order[sd-1][f] for f in facets)
+            sorder = forder + sd - 1
+            if vorder > sorder:
+                jumps = expansion_set.tabulate_jumps(degree, verts[vids], order=vorder)
+                rows.extend(numpy.vstack(jumps[r].T) for r in range(sorder+1, vorder+1))
 
         if len(rows) > 0:
+            for row in rows:
+                row *= 1 / max(numpy.max(abs(row)), 1)
             dual_mat = numpy.vstack(rows)
             coeffs = polynomial_set.spanning_basis(dual_mat, nullspace=True)
         else:
             coeffs = numpy.eye(expansion_set.get_num_members(degree))
 
-        # Impose C^vorder super-smoothness at interior vertices
-        # C^order automatically gives C^{order+dim-1} at the interior vertex
-        sorder = order + sd - 1
-        if vorder > sorder:
-            verts = ref_el.get_vertices()
-            points = [verts[i] for i in ref_el.get_interior_facets(0)]
-            jumps = expansion_set.tabulate_jumps(degree, points, order=vorder)
-            for r in range(sorder+1, vorder+1):
-                dual_mat = numpy.dot(numpy.vstack(jumps[r].T), coeffs.T)
-                nsp = polynomial_set.spanning_basis(dual_mat, nullspace=True)
-                coeffs = numpy.dot(nsp, coeffs)
-
         if shape != ():
             m, n = coeffs.shape
             ncomp = numpy.prod(shape)

FIAT/polynomial_set.py

@@ -104,6 +104,9 @@ def take(self, items):
         return PolynomialSet(self.ref_el, self.degree, self.embedded_degree,
                              self.expansion_set, new_coeffs)
 
+    def __len__(self):
+        return self.num_members
+
 
 class ONPolynomialSet(PolynomialSet):
     """Constructs an orthonormal basis out of expansion set by having an