- changed title to Quad/hex meshes need ordering check
- edited description
Quad/hex meshes need ordering check
Assembly on quads and hexes works only by a fluke - generated meshes are already ordered but imported or MeshEditor
ed 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))
Comments (7)
-
reporter -
reporter - edited description
-
Note that generic hex meshes are not consistently orientable. See https://arxiv.org/abs/1512.02137
-
reporter Fix in branch
jan/fix-issue-997
. -
reporter Fix in pull request #476.
-
reporter - changed status to resolved
-
reporter 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.
- Log in to comment