Quadrature representation doesn't produce correct code for oriented interior facet integrals

[Lawrence Mitchell]

Marie fixed this issue in f63d6985d5 (merged as 80efbc2f34). However, with current master (362a2a), The interior facet orientation code is no longer restricted. I don't know quite what's gone on, but this fix is easy.

diff --git a/ffc/quadrature/quadraturegenerator.py b/ffc/quadrature/quadraturegenerator.py
index f0d97a3..85a9f68 100644
--- a/ffc/quadrature/quadraturegenerator.py
+++ b/ffc/quadrature/quadraturegenerator.py
@@ -187,7 +187,7 @@ def _tabulate_tensor(ir, prefix, parameters):
             jacobi_code += "\n"
             jacobi_code += format["compute_jacobian_inverse"](tdim, gdim, r=_r)
             if oriented:
-                jacobi_code += format["orientation"](tdim, gdim)
+                jacobi_code += format["orientation"](tdim, gdim, r=_r)
             jacobi_code += "\n"
         jacobi_code += "\n\n" + format["facet determinant"](tdim, gdim, r="+")
         jacobi_code += "\n\n" + format["generate normal"](tdim, gdim, integral_type)

[Marie Elisabeth Rognes]

Thanks for reporting Lawrence. Could you give us a small test too so that this regression does not happen again?

[Lawrence Mitchell]

So this is, I think, easier to test in dolfin than FFC (I'm not really sure how to do the latter).

Here's a patch against dolfin that adds another test to the manifolds test suite:

diff --git a/test/unit/fem/python/manifolds.py b/test/unit/fem/python/manifolds.py
index 0cb5314..6395055 100644
--- a/test/unit/fem/python/manifolds.py
+++ b/test/unit/fem/python/manifolds.py
@@ -273,6 +273,82 @@ class ManifoldBasisEvaluation(unittest.TestCase):
                     self.assertAlmostEqual(abs(derivs_rot2-derivs_cmp).max(), 0.0, 10)
                     self.assertAlmostEqual(abs(values_cmp-values_rot).max(), 0.0, 10)

+
+class ManifoldInteriorFacet(unittest.TestCase):
+
+    def test_contravariant_piola_facet_integral(self):
+        import math
+        phi = (1 + math.sqrt(5)) / 2
+
+        # Surface of icosahedron
+        vertices = np.array([[-1, phi, 0],
+                             [1, phi, 0],
+                             [-1, -phi, 0],
+                             [1, -phi, 0],
+                             [0, -1, phi],
+                             [0, 1, phi],
+                             [0, -1, -phi],
+                             [0, 1, -phi],
+                             [phi, 0, -1],
+                             [phi, 0, 1],
+                             [-phi, 0, -1],
+                             [-phi, 0, 1]])
+
+        faces = np.array([[0, 11, 5],
+                          [0, 5, 1],
+                          [0, 1, 7],
+                          [0, 7, 10],
+                          [0, 10, 11],
+                          [1, 5, 9],
+                          [5, 11, 4],
+                          [11, 10, 2],
+                          [10, 7, 6],
+                          [7, 1, 8],
+                          [3, 9, 4],
+                          [3, 4, 2],
+                          [3, 2, 6],
+                          [3, 6, 8],
+                          [3, 8, 9],
+                          [4, 9, 5],
+                          [2, 4, 11],
+                          [6, 2, 10],
+                          [8, 6, 7],
+                          [9, 8, 1]], dtype=np.uintp)
+
+        me = MeshEditor()
+
+        mesh = Mesh()
+
+        me.open(mesh, 2, 3)
+
+        me.init_vertices(vertices.shape[0])
+
+        me.init_cells(faces.shape[0])
+
+        for i, points in enumerate(vertices):
+            me.add_vertex(i, points)
+
+        for i, cell in enumerate(faces):
+            me.add_cell(i, cell)
+
+        me.close()
+
+        mesh.order()
+
+        mesh.init_cell_orientations(Expression(('x[0]', 'x[1]', 'x[2]')))
+
+        V = FunctionSpace(mesh, 'RT', 1)
+
+        u = project(Expression(('x[0]', '-x[1]', '0')), V)
+        n = FacetNormal(mesh)
+
+        pos = inner(u('+'), n('+'))*dS
+        neg = inner(u('-'), n('-'))*dS
+
+        parameters['form_compiler']['representation'] = 'quadrature'
+        self.assertAlmostEqual(assemble(pos + neg), 0.0, 10)
+
+
 if __name__ == "__main__":
     print ""
     print "Testing solving and evaluate basis over manifolds"

[Marie Elisabeth Rognes]

This issue is fixed. Just the test that needs to be added.

[Martin Sandve Alnæs]

Removing milestone: 1.5 (automated comment)