Quad/hex meshes need ordering check

[Jan Blechta]

Assembly on quads and hexes works only by a fluke - generated meshes are already ordered but imported or MeshEditored quad/hex meshes are not guaranteed to be reorderable, moreover it is not checked that the crucial ordering property is fulfilled. Once the property is violated, it leads to a catastrophe, see MWE.

Solution is to implement QuadrilateralCell::ordered() and HexahedronCell::ordered() and let (Quadrilateral|Hexahedron)Cell::order() just call ordered() (rather than to do any ordering).

from dolfin import *

mesh = Mesh()
e = MeshEditor()
e.open(mesh, "quadrilateral", 2, 2, degree=1)
e.init_vertices(6)
e.init_cells(2)
e.add_vertex(0, Point(0, 0))
e.add_vertex(1, Point(1, 0))
e.add_vertex(2, Point(2, 0))
e.add_vertex(3, Point(0, 1))
e.add_vertex(4, Point(1, 1))
e.add_vertex(5, Point(2, 1))
e.add_cell(0, [0, 1, 3, 4])
e.add_cell(1, [4, 5, 1, 2])  #bad
#e.add_cell(1, [1, 2, 4, 5])  # good
e.close()

Q1 = FunctionSpace(mesh, "Q", 1)
u = Function(Q1)
zero = (u('+')-u('-'))**2*dS  # should be zero by definition
for i in range(Q1.dim()):
    u.vector().zero()
    u.vector()[i] = 1
    print("zero = ", assemble(zero))

[Lawrence Mitchell]

Note that generic hex meshes are not consistently orientable. See https://arxiv.org/abs/1512.02137

[Jan Blechta]

Fix in branch jan/fix-issue-997.

[Jan Blechta]

Fix in pull request #476.

[Jan Blechta]

Fixed at f01c707d62681e2b13e84aca1d18fd090a794300

[Jan Blechta]

As a follow-up, the following test shows that mesh with bad ordering obtained through IO also raises.

Run this on 2017.2.0/SWIG to produce a bad mesh

from dolfin import *
import numpy as np

mesh = Mesh()
e = MeshEditor()
e.open(mesh, "quadrilateral", 2, 2, 1)
e.init_vertices(6)
e.init_cells(2)
e.add_vertex(0, Point(0, 0))
e.add_vertex(1, Point(1, 0))
e.add_vertex(2, Point(2, 0))
e.add_vertex(3, Point(0, 1))
e.add_vertex(4, Point(1, 1))
e.add_vertex(5, Point(2, 1))
e.add_cell(0, np.array([0, 1, 3, 4], dtype=np.uintp))
e.add_cell(1, np.array([4, 5, 1, 2], dtype=np.uintp))
e.close()

with XDMFFile('mesh.xdmf') as f:
    f.write(mesh)

and then run on fixed version (2018.1.0 and later)

from dolfin import *
mesh = Mesh()
with XDMFFile(mesh.mpi_comm(), 'mesh.xdmf') as f:
    f.read(mesh)

and expect that error is raised:

*** Error:   Unable to order quadrilateral cell.
*** Reason:  Cell is not orderable.
*** Where:   This error was encountered inside QuadrilateralCell.cpp.