Check geometry lowering of facet normal on quads/hexes

[Jan Blechta]

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/.

[Lawrence Mitchell]

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:

from firedrake import *

for quad in [False, True]:
    mesh = UnitSquareMesh(16, 16, quadrilateral=quad)
    V = FunctionSpace(mesh, "CG", 1)
    x, y = SpatialCoordinate(mesh)
    u = interpolate(x**2 + y**2, V)
    n = FacetNormal(mesh)
    # 4 is the plane y == 1
    j = assemble(dot(grad(u), n)*ds(4))
    k = assemble(dot(grad(u), as_vector([0, 1]))*ds(4))
    print("%s, grad(u).n: %s, grad(u).(0, 1): %s" % ({False: "triangle",
                                                      True: "quadrilateral"}[quad],
                                                     j, k))

=>

triangle, grad(u).n: 1.9375, grad(u).(0, 1): 1.9375
quadrilateral, grad(u).n: 1.9375, grad(u).(0, 1): 1.9375

[Michal Habera]

Another suggestion, the integral of dot(n, n) is (wrongly) independent on the measure (size of ds(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

[Jan Blechta]

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.

[Michal Habera]

@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).

[Ivan Yashchuk]

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.

[Michal Habera]

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.

[Jan Blechta]

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.]])

[Jan Blechta]

https://bitbucket.org/fenics-project/dolfin/issues/991

[Jan Blechta]

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.

[Jan Blechta]

Seems like a bug in ufc_geometry.h. Fix in https://bitbucket.org/fenics-project/ffc/commits/00276c2c03e0987b65356357d51014ae9fac36b3

[Jan Blechta]

Seems like it fixes DG on quads too.

[Jan Blechta]

Seems resolved.