Negative element volume and Jacobian in UnitCubeMesh

I was trying to export a UnitCubeMesh for use in a different solver but it appears that some of the tetrahedral elements are inverted. When I checked the mesh in ParaView using the MeshQuality and CellSize filters it seems that half the elements have negative volume (see below) and the total volume is zero. Does anyone know why this might be happening and is there a way to fix it?

When I try to open the exported mesh in Abaqus it looks like this:

Here is the code that I used to generate the mesh:

from dolfin import *
import numpy as np

parameters["reorder_dofs_serial"]=False

mesh=UnitCubeMesh(2,2,2)
# mesh=BoxMesh(Point(0,0,0),Point(1,1,1),2,2,2)
nodes=mesh.coordinates()
element=mesh.cells()

np.savetxt("nodes_N2.dat",nodes,delimiter=",")
np.savetxt("element_N2.dat",element,delimiter=",")

file=File("mesh_N2.pvd");
file<<mesh;
1 Like

I tried using the latest dolfinx and it appears to be fixed there.

Here is the code to generate the mesh:

import dolfinx
import dolfinx.io
from mpi4py import MPI

mesh = dolfinx.UnitCubeMesh(MPI.COMM_WORLD, 5, 5, 5)

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "w") as f:
    f.write_mesh(mesh)

f = dolfinx.io.VTKFile("mesh.pvd")
f.write(mesh)

Comparing dolfin and dolfinx the element connectivity is slightly different (see below). It seems to me that this is a bug in dolfin that has been fixed in dolfinx. Is there a fix that can be applied to the stable dolfin version?

dolfin input:

import dolfin

mesh = dolfin.UnitCubeMesh(1, 1, 1)

with dolfin.XDMFFile("dolfin.xdmf") as f:
    f.write(mesh)

f = dolfin.File("dolfin.pvd")
f.write(mesh)

mesh.cells()

dolfin output:

array([[0, 1, 3, 7],
       [0, 1, 5, 7],
       [0, 4, 5, 7],
       [0, 2, 3, 7],
       [0, 4, 6, 7],
       [0, 2, 6, 7]], dtype=uint32)

dolfinx input:

import dolfinx
import dolfinx.io
from mpi4py import MPI

mesh = dolfinx.UnitCubeMesh(MPI.COMM_WORLD, 1, 1, 1)

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "dolfinx.xdmf", "w") as f:
    f.write_mesh(mesh)

f = dolfinx.io.VTKFile("dolfinx.pvd")
f.write(mesh)

mesh.topology.create_connectivity_all()
mesh.topology.connectivity(3,0)

dolfinx output:

<AdjacencyList> with 6 nodes
  0: [0 1 3 7 ]
  1: [0 1 7 5 ]
  2: [0 5 7 4 ]
  3: [0 3 2 7 ]
  4: [0 6 4 7 ]
  5: [0 2 6 7 ]