Quadrature representation doesn't produce correct code for oriented interior facet integrals
Issue #25
resolved
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)
Comments (5)
-
-
reporter 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"
-
- changed milestone to 1.5
-
- changed status to resolved
-
assigned issue to
This issue is fixed. Just the test that needs to be added.
-
- removed milestone
Removing milestone: 1.5 (automated comment)
- Log in to comment
Thanks for reporting Lawrence. Could you give us a small test too so that this regression does not happen again?