Make HDiv matrices less magical

.github/workflows/test.yml

@@ -36,6 +36,7 @@ jobs:
       PYTEST_MPI_MAX_NPROCS: 8
       # TODO: reduce this 
       FIREDRAKE_RUN_SPLIT_TESTS_TIMEOUT: 3600
+      LOOPY_NO_CACHE: 1
     steps:
       - name: Fix HOME
         # For unknown reasons GitHub actions overwrite HOME to /github/home

fuse/cells.py

@@ -1401,7 +1401,7 @@ def __init__(self, cell, name=None, renumber=False):
         #     breakpoint()
 
     def cellname(self):
-        return "FUSE_" + self.name
+        return "FUSE_" + self.name + repr(self.cell)
 
     def construct_subelement(self, dimension, e_id=0, o=None):
         """Constructs the reference element of a cell
@@ -1441,7 +1441,7 @@ def __init__(self, cell, name=None):
         super(CellComplexToFiatTensorProduct, self).__init__(*fiat_factors)
 
     def cellname(self):
-        return "FUSE_" + self.name
+        return "FUSE_" + self.name + repr(self.cell)
 
     def construct_subelement(self, dimension):
         """Constructs the reference element of a cell
@@ -1473,7 +1473,7 @@ def __init__(self, cell, product):
         super(CellComplexToFiatHypercube, self).__init__(product.get_spatial_dimension(), product)
 
     def cellname(self):
-        return "FUSE_" + self.name
+        return "FUSE_" + self.name + repr(self.cell)
 
     def construct_subelement(self, dimension):
         """Constructs the reference element of a cell
@@ -1540,7 +1540,7 @@ def to_fiat(self):
         return self.cell_complex.to_fiat(name=self.cellname)
 
     def __repr__(self):
-        return "FUSE_" + super(CellComplexToUFL, self).__repr__()
+        return "FUSE_" + super(CellComplexToUFL, self).__repr__() + repr(self.cell_complex)
 
     def reconstruct(self, **kwargs):
         """Reconstruct this cell, overwriting properties by those in kwargs."""

fuse/element_construction.py

@@ -453,7 +453,7 @@ def construct_interval_cgN(deg, cell=None):
     vert = cell.vertices()[0]
 
     xs = [DOF(DeltaPairing(), PointKernel(()))]
-    dg0 = ElementTriple(vert, (P0, CellL2, C0), DOFGenerator(xs, S1, S1))
+    dg0 = ElementTriple(vert, (P0, C0, Fid), DOFGenerator(xs, S1, S1))
     v_xs = [immerse(cell, dg0, TrH1)]
     v_dofs = [DOFGenerator(v_xs, get_cyc_group(len(cell.vertices())), S1)]
 
@@ -467,7 +467,7 @@ def construct_interval_cgN(deg, cell=None):
     else:
         centre_dof = []
 
-    cg = ElementTriple(cell, (Pk, CellH1, C0), v_dofs + sym_dofs + centre_dof)
+    cg = ElementTriple(cell, (Pk, C0, Fid), v_dofs + sym_dofs + centre_dof)
     assert len(cg.generate()) == deg + 1
     return cg
 
@@ -477,7 +477,7 @@ def construct_interval_dgN_integral(deg, cell=None):
         cell = line()
     Pk = PolynomialSpace(deg)
     dofs = lagrange_facet_fns(cell, deg, interior=True, vector=False)
-    dg = ElementTriple(cell, (Pk, CellL2, C0), dofs)
+    dg = ElementTriple(cell, (Pk, C0, Fid), dofs)
     assert len(dg.generate()) == deg + 1
     return dg
 

fuse/enriched.py

@@ -19,18 +19,20 @@ def __init__(self, A, B, flat=False, symmetric=None, matrices=True):
             raise ValueError("EnrichedElement should only be used for Tensor product elements. Use + between triples for enrichment.")
         self.A = A
         self.B = B
-        self.spaces = (A.spaces[0] + B.spaces[0], A.spaces[1], max([A.spaces[2], B.spaces[2]]))
+        if A.spaces[2] != B.spaces[2]:
+            raise NotImplementedError("Two FUSE triples with different pullbacks cannot be added.")
+        self.spaces = (A.spaces[0] + B.spaces[0], max(A.spaces[1], B.spaces[1]), A.spaces[2])
 
         self.DOFGenerator = [A.DOFGenerator, B.DOFGenerator]
         if A.cell.flat != B.cell.flat:
             raise ValueError("Tensor products must both be flat or both not flat for enrichment.")
         self.cell = A.cell
-        # Derived, not trusted: an enriched element whose cell is flat but
-        # which reports flat=False silently skips the axis-swap fill and the
-        # DOF regrouping.
         self.flat = flat or self.cell.flat
         if hasattr(A, "unflat_cell"):
             self.unflat_cell = A.unflat_cell
+        if getattr(A, "trace", None) is not getattr(B, "trace", None):
+            raise ValueError("Cannot enrich elements with different traces.")
+        self.trace = getattr(A, "trace", None)
         # See TensorProductTriple.__init__ for the meaning of ``symmetric``.
         self.requested_symmetric = symmetric
         self.symmetric = True if symmetric is None else symmetric
@@ -65,6 +67,8 @@ def setup_matrices(self):
         self.B.to_ufl()
         dofs = self.generate()
         dof_keys, key_to_index = self._axis_key_maps(dofs)
+        # Reset closure failures
+        self._closure_failures = set()
         oriented_mats_by_entity, flat_by_entity = self._initialise_entity_dicts(dofs, tensor=(not self.flat))
         if self.flat:
             cell = self.A.unflat_cell
@@ -95,7 +99,8 @@ def setup_matrices(self):
                                                     [np.zeros((b_sub_mat.shape[0], a_sub_mat.shape[1])), b_sub_mat]])
                         sub_mat[o][np.ix_(ent_dofs, ent_dofs)] = np.matmul(sub_mat[o][np.ix_(ent_dofs, ent_dofs)], combined_sub_mat)
                     if self.flat:
-                        self._fill_face_axis_swaps(dim, ent_dofs, sub_mat, dof_keys, key_to_index)
+                        entity = self.cell.d_entities(total_dim)[e]
+                        self._fill_axis_permutations(entity, dim, ent_dofs, sub_mat, dof_keys, key_to_index)
 
         self.matrices = oriented_mats_by_entity
         self.reversed_matrices = self.reverse_dof_perms(self.matrices)