Support physical_normals on manifolds

.github/workflows/core.yml

@@ -218,6 +218,8 @@ jobs:
             --extra-index-url https://download.pytorch.org/whl/cpu \
             "$(echo ./firedrake-repo/dist/firedrake-*.tar.gz)[ci]"
 
+          : # UNDO ME
+          pip install -v --no-deps --ignore-installed git+https://github.com/firedrakeproject/fiat.git@pbrubeck/zany-manifold
           firedrake-clean
           pip list
 

tsfc/fem.py

@@ -268,19 +268,28 @@ def physical_tangents(self):
         return rts @ jac.T
 
     def physical_normals(self):
+        domain = extract_unique_domain(self.mt.terminal)
+        gdim = domain.geometric_dimension
+        tdim = domain.topological_dimension
         cell = self.interface.fiat_cell
-        sd = cell.get_spatial_dimension()
-        num_faces = len(cell.get_topology()[sd-1])
-        if isinstance(cell, UFCSimplex) and sd == 2:
+        num_faces = len(cell.get_topology()[tdim-1])
+        if isinstance(cell, UFCSimplex) and tdim == 1:
+            pts = self.physical_tangents()
+            return gem.ListTensor([[pts[0, j] for j in range(gdim)] for i in range(num_faces)])
+        elif isinstance(cell, UFCSimplex) and tdim == 2 and gdim == 2:
             pts = self.physical_tangents()
-            return gem.ListTensor([[pts[i, 1], -1*pts[i, 0]] for i in range(num_faces)])
-        elif isinstance(cell, UFCSimplex) and sd == 3:
-            t = ufl.classes.CellEdgeVectors(extract_unique_domain(self.mt.terminal))
-            edges = cell.get_connectivity()[(sd-1, 1)]
+            return gem.ListTensor([[pts[i, 1], -pts[i, 0]] for i in range(num_faces)])
+        elif isinstance(cell, UFCSimplex) and gdim == 3:
+            t = ufl.classes.CellEdgeVectors(domain)
             normalize = lambda x: x / ufl.sqrt(ufl.dot(x, x))
-            expr = ufl.as_tensor([-2.0*normalize(ufl.cross(t[edges[i][0], :], t[edges[i][1], :]))
-                                  for i in range(num_faces)])
-            return self.translate_point_expression(expr)
+            if tdim == 2:
+                R = ufl.cross(t[0, :], t[1, :])
+                exprs = [normalize(ufl.cross(t[i, :], R)) for i in range(num_faces)]
+            else:
+                edges = cell.get_connectivity()[(tdim-1, 1)]
+                exprs = [-2.0*normalize(ufl.cross(t[edges[i][0], :], t[edges[i][1], :]))
+                         for i in range(num_faces)]
+            return self.translate_point_expression(ufl.as_tensor(exprs))
         else:
             raise NotImplementedError("Can't do physical normals on that cell yet")
 

tsfc/kernel_interface/common.py

@@ -346,6 +346,10 @@ def set_quad_rule(params, cell, integral_type, functions):
         finat_elements = set(create_element(e) for e in elements if e.family() != "Real")
         fiat_cells = [fiat_cell] + [finat_el.complex for finat_el in finat_elements]
         if any(c.is_macrocell() for c in fiat_cells):
+            if len(set(c.get_spatial_dimension() for c in fiat_cells)) > 1:
+                dimension = fiat_cell.get_dimension()
+                fiat_cells = [c if c.get_dimension() == dimension else
+                              c.construct_subcomplex(dimension) for c in fiat_cells]
             fiat_cell = max_complex(fiat_cells)
         integration_dim, _ = lower_integral_type(fiat_cell, integral_type)
         quad_rule = fem.get_quadrature_rule(fiat_cell, integration_dim, quadrature_degree, scheme)