Add tests for 3D zany elements (#4459)

* TSFC: physical_normals

* Test 3D zany elements

Author: pbrubeck

tests/firedrake/regression/test_projection_zany.py

@@ -24,11 +24,17 @@ def jrc(a, b, n):
 
 
 @pytest.fixture
-def hierarchy(request):
+def hierarchy_2d(request):
     msh = UnitSquareMesh(2, 2)
     return MeshHierarchy(msh, 4)
 
 
+@pytest.fixture
+def hierarchy_3d(request):
+    n = 2
+    return [UnitCubeMesh(n*2**r, n*2**r, n*2**r) for r in range(3)]
+
+
 def do_projection(mesh, el_type, degree):
     V = FunctionSpace(mesh, el_type, degree)
 
@@ -57,6 +63,12 @@ def do_projection(mesh, el_type, degree):
     return sqrt(assemble(inner(fh - f, fh - f) * dx))
 
 
+def run_convergence_test(mh, el, degree, convrate):
+    err = [do_projection(m, el, degree) for m in mh]
+    conv = np.diff(-np.log2(err))
+    assert (conv > convrate).all()
+
+
 @pytest.mark.parametrize(('el', 'deg', 'convrate'),
                          [('Johnson-Mercier', 1, 1.8),
                           ('Morley', 2, 2.4),
@@ -69,10 +81,15 @@ def do_projection(mesh, el_type, degree):
                           ('Bell', 5, 4.7),
                           ('Argyris', 5, 5.8),
                           ('Argyris', 6, 6.7)])
-def test_projection_zany_convergence_2d(hierarchy, el, deg, convrate):
-    diff = np.array([do_projection(m, el, deg) for m in hierarchy[2:]])
-    conv = np.log2(diff[:-1] / diff[1:])
-    assert (np.array(conv) > convrate).all()
+def test_projection_zany_convergence_2d(hierarchy_2d, el, deg, convrate):
+    run_convergence_test(hierarchy_2d[2:], el, deg, convrate)
+
+
+@pytest.mark.parametrize(('el', 'deg', 'convrate'),
+                         [('Johnson-Mercier', 1, 1.8),
+                          ('Morley', 2, 2.4)])
+def test_projection_zany_convergence_3d(hierarchy_3d, el, deg, convrate):
+    run_convergence_test(hierarchy_3d, el, deg, convrate)
 
 
 @pytest.mark.parametrize(('element', 'degree'),
@@ -81,9 +98,9 @@ def test_projection_zany_convergence_2d(hierarchy, el, deg, convrate):
                           ('HCT', 4),
                           ('Argyris', 5),
                           ('Argyris', 6)])
-def test_mass_conditioning(element, degree, hierarchy):
+def test_mass_conditioning(element, degree, hierarchy_2d):
     mass_cond = []
-    for msh in hierarchy[1:4]:
+    for msh in hierarchy_2d[1:4]:
         V = FunctionSpace(msh, element, degree)
         u = TrialFunction(V)
         v = TestFunction(V)

tsfc/fem.py

@@ -155,36 +155,33 @@ def config(self):
         config["interface"] = self.interface
         return config
 
-    def cell_size(self):
-        return self.interface.cell_size(self.mt.restriction)
-
-    def jacobian_at(self, point):
-        ps = PointSingleton(point)
-        expr = Jacobian(extract_unique_domain(self.mt.terminal))
-        assert ps.expression.shape == (extract_unique_domain(expr).topological_dimension(), )
+    def translate_point_expression(self, expr, point=None):
         if self.mt.restriction == '+':
             expr = PositiveRestricted(expr)
         elif self.mt.restriction == '-':
             expr = NegativeRestricted(expr)
+
+        cell = self.interface.fiat_cell
+        sd = cell.get_spatial_dimension()
+        if point is None:
+            point, = cell.make_points(sd, 0, sd+1)
         config = {"point_set": PointSingleton(point)}
         config.update(self.config)
         config.update(use_canonical_quadrature_point_ordering=False)  # quad point ordering not relevant.
         context = PointSetContext(**config)
         expr = self.preprocess(expr, context)
         return map_expr_dag(context.translator, expr)
 
+    def cell_size(self):
+        return self.interface.cell_size(self.mt.restriction)
+
+    def jacobian_at(self, point):
+        expr = Jacobian(extract_unique_domain(self.mt.terminal))
+        return self.translate_point_expression(expr, point=point)
+
     def detJ_at(self, point):
         expr = JacobianDeterminant(extract_unique_domain(self.mt.terminal))
-        if self.mt.restriction == '+':
-            expr = PositiveRestricted(expr)
-        elif self.mt.restriction == '-':
-            expr = NegativeRestricted(expr)
-        config = {"point_set": PointSingleton(point)}
-        config.update(self.config)
-        config.update(use_canonical_quadrature_point_ordering=False)  # quad point ordering not relevant.
-        context = PointSetContext(**config)
-        expr = self.preprocess(expr, context)
-        return map_expr_dag(context.translator, expr)
+        return self.translate_point_expression(expr, point=point)
 
     def reference_normals(self):
         cell = self.interface.fiat_cell
@@ -210,30 +207,25 @@ def physical_tangents(self):
 
     def physical_normals(self):
         cell = self.interface.fiat_cell
-        if not (isinstance(cell, UFCSimplex) and cell.get_dimension() == 2):
+        sd = cell.get_spatial_dimension()
+        num_faces = len(cell.get_topology()[sd-1])
+        if isinstance(cell, UFCSimplex) and sd == 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)]
+            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)
+        else:
             raise NotImplementedError("Can't do physical normals on that cell yet")
 
-        num_edges = len(cell.get_topology()[1])
-        pts = self.physical_tangents()
-        return gem.ListTensor([[pts[i, 1], -1*pts[i, 0]] for i in range(num_edges)])
-
     def physical_edge_lengths(self):
         expr = ufl.classes.CellEdgeVectors(extract_unique_domain(self.mt.terminal))
-        if self.mt.restriction == '+':
-            expr = PositiveRestricted(expr)
-        elif self.mt.restriction == '-':
-            expr = NegativeRestricted(expr)
-
-        cell = self.interface.fiat_cell
-        sd = cell.get_spatial_dimension()
-        num_edges = len(cell.get_topology()[1])
-        expr = ufl.as_vector([ufl.sqrt(ufl.dot(expr[i, :], expr[i, :])) for i in range(num_edges)])
-        config = {"point_set": PointSingleton(cell.make_points(sd, 0, sd+1)[0])}
-        config.update(self.config)
-        config.update(use_canonical_quadrature_point_ordering=False)  # quad point ordering not relevant.
-        context = PointSetContext(**config)
-        expr = self.preprocess(expr, context)
-        return map_expr_dag(context.translator, expr)
+        expr = ufl.as_vector([ufl.sqrt(ufl.dot(expr[i, :], expr[i, :])) for i in range(expr.ufl_shape[0])])
+        return self.translate_point_expression(expr)
 
     def physical_points(self, point_set, entity=None):
         """Converts point_set from reference to physical space"""