Wrong mesh for quadratic xdmf triangles

I regularly import gmsh meshes via meshio to use them in dolfin. This works fine for linear triangles. But I realized wrong/strange results in my calculations for quadratic ones.

See the the two attached xdmf-xml files (same issue for xdmf-hd5) for a minimal example. Each file contains a quadratic unit triangle (x y = 0 0 -- 1 0 -- 0 1) with the same geometry but different numbering. When imported to paraview, both look fine.

2          1     
|\         |\
| \        | \
4  5  vs.  5  2
|   \      |   \
|    \     |    \
1--3--0    3--4--0

[012345]   [031452]

The [012345] numbering works fine, I think. But I realized an issue with the [031452] numbering:

wrong area calculation (1/6 instead of 1/2)
overlapping mesh when refined

from dolfin import *
from matplotlib.pyplot import show

def from_xdmf(xdmf_file):
    mesh = Mesh()
    with XDMFFile(xdmf_file) as f:
        f.read(mesh)
    return mesh

def area(mesh):
    return assemble(1 * Measure('dx', mesh))

print('Expected  :  ', 0.5)
print('ordering 1:  ', area(from_xdmf('ordering012345.xdmf'))) # ~0.5
print('ordering 2:  ', area(from_xdmf('ordering031452.xdmf'))) # ~0.166

plot(refine(from_xdmf('ordering031452.xdmf')), marker='o')
show()

Comments (1)