Check geometry lowering of facet normal on quads/hexes
Here might be a wrong assumption taken.
There's an evidence that it does not work: https://www.allanswered.com/post/qvqrl/how-do-i-integrate-including-the-facet-normal-with-quadrilateral-elements/.
Comments (13)
-
-
Another suggestion, the integral of
dot(n, n)
is (wrongly) independent on the measure (size ofds(1)
).from dolfin import * for ctype in [CellType.Type.triangle, CellType.Type.quadrilateral]: mesh = RectangleMesh.create(MPI.comm_world, [Point(0, 0), Point(3.5, 1)], [16, 16], ctype) # Mark the top edge ff = MeshFunction("size_t", mesh, 1, 0) CompiledSubDomain("near(x[1], 1.0)").mark(ff, 1) ds = Measure("ds", subdomain_data=ff) # Compute functionals n = FacetNormal(mesh) j = assemble(dot(n, n)*ds(1)) print("CellType.%s, n.n: %.3e" % (ctype, j))
=>
CellType.Type.triangle, n.n: 3.500e+00 CellType.Type.quadrilateral, n.n: 1.000e+00
-
reporter Michal's example fails in the same way with
parameters['form_compiler']['representation'] = 'tsfc'
. That suggests that the problem is in DOLFIN. Perhaps entity ordering. -
@ivanyashchuk Do you know about this issue? It could well be the reason why DG Poisson demo on quad mesh produced wrong results (facet normal on quads seems to be crapped).
-
I don't know about the issue. The dx measure is incorrect as well.
from dolfin import * for ctype in [CellType.Type.triangle, CellType.Type.quadrilateral]: mesh = RectangleMesh.create(MPI.comm_world, [Point(0, 0), Point(3.5, 2)], [16, 16], ctype) # Compute functionals j = assemble(1*dx(mesh)) print("CellType.%s, 1.dx: %.3e" % (ctype, j))
CellType.Type.triangle, 1.dx: 7.000e+00 CellType.Type.quadrilateral, 1.dx: 1.000e+00
I don't remember noticing similar behavior in September. I think dx measure worked for quads correctly and I think I've tried meshes of different size. Now I am not so sure if I actually did it.
-
Good point, the problem is in Mesh generation via RectangleMesh, the Points are ignored here https://bitbucket.org/fenics-project/dolfin/src/773843a63c4f7d6e293cbdbbcc946c09b773d7b0/dolfin/generation/RectangleMesh.cpp?at=master&fileviewer=file-view-default#RectangleMesh.cpp-259 still assumes (1.0, 1.0) mesh.
-
reporter You have hit here an orthogonal bug
mesh = RectangleMesh.create(MPI.comm_world, [Point(0, 0), Point(6666666, 77777777)], [1, 1], CellType.Type.quadrilateral) print(mesh.coordinates())
array([[ 0., 0.], [ 1., 0.], [ 0., 1.], [ 1., 1.]])
-
reporter -
reporter The original Nate's problem is fixed by adding the line
parameters['form_compiler']['representation'] = 'tsfc'
. I guess that UFLACS does not handle correctly code generation for non-constant Jacobian, which is used in geometry lowering for facet normal. -
reporter Seems like a bug in
ufc_geometry.h
. Fix in https://bitbucket.org/fenics-project/ffc/commits/00276c2c03e0987b65356357d51014ae9fac36b3 -
reporter Seems like it fixes DG on quads too.
-
reporter - removed responsible
- removed milestone
-
reporter - changed status to resolved
Seems resolved.
- Log in to comment
FWIW, here's evidence to suggest it is perhaps not a geometry lowering issue. In that the posted code works fine using a tsfc->firedrake pipeline:
=>