I am trying to import a 3D mesh consisting of hexahedrons, quadrangles and lines created in Gmsh (Gmsh script provided after the MWE). The following MWE tries to read the mesh and calculate areas of surfaces, lengths of sides, and the volume. I am using dolfin 2019.2.0.dev0
MWE:
from dolfin import *
import meshio
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True}
import numpy
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
return out_mesh
############
msh = meshio.read("mesh.msh")
############
hexa_mesh = create_mesh(msh, "hexahedron", True)
quad_mesh = create_mesh(msh, "quad", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", hexa_mesh)
meshio.write("surface.xdmf", quad_mesh)
meshio.write("mf.xdmf", line_mesh)
mesh = Mesh()
xdmf = XDMFFile(mesh.mpi_comm(),"mesh.xdmf")
xdmf.read(mesh) # Generates error
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
xdmf.close()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("surface.xdmf") as infile:
infile.read(mvc, "name_to_read")
sf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
xdmf.close()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-2)
with XDMFFile("mf.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
def assemble_edge(tag):
mv = MeshView.create(mf, tag)
return assemble(Constant(1.0)*dx(domain=mv))
dx_custom = Measure("ds", domain=mesh, subdomain_data=sf)
dv_custom = Measure("dx", domain=mesh, subdomain_data=cf)
print(f"Volume = {assemble(Constant(1.0)*dv_custom(19))}")
print(f"frontbottomline = {assemble_edge(4)}")
print(f"fronttopline = {assemble_edge(3)}")
print(f"frontleftline = {assemble_edge(1)}")
print(f"frontrightline = {assemble_edge(2)}")
print(f"lefttopline = {assemble_edge(5)}")
print(f"leftbottomline = {assemble_edge(6)}")
print(f"righttopline = {assemble_edge(7)}")
print(f"rightbottomline = {assemble_edge(8)}")
print(f"backtopline = {assemble_edge(9)}")
print(f"backbottomline ={assemble_edge(10)}")
print(f"backleftline = {assemble_edge(11)}")
print(f"backrightline = {assemble_edge(12)}")
print(f"frontsurface = {assemble(Constant(1.0)*dx_custom(13))}")
print(f"leftsurface = {assemble(Constant(1.0)*dx_custom(14))}")
print(f"rightsurface = {assemble(Constant(1.0)*dx_custom(15))}")
print(f"topsurface = {assemble(Constant(1.0)*dx_custom(16))}")
print(f"bottomsurface = {assemble(Constant(1.0)*dx_custom(17))}")
print(f"backsurface = {assemble(Constant(1.0)*dx_custom(18))}")
This generates the following error:
-------------------------------------------------------------------------
*** Error: Unable to order hexahedron cell.
*** Reason: Cell is not orderable.
*** Where: This error was encountered inside HexahedronCell.cpp.
Here is the Gmsh script:
SetFactory("OpenCASCADE");
//+
Point(1) = {0.0, 0.25, 0.0, 1.0};
//+
Point(2) = {9.0, 0.25, 0.0, 1.0};
//+
Point(3) = {11.0, 0.25, 0.0, 1.0};
//+
Point(4) = {20.0, 0.25, 0.0, 1.0};
//+
Point(5) = {0.0, -0.25, 0.0, 1.0};
//+
Point(6) = {9.0, -0.25, 0.0, 1.0};
//+
Point(7) = {11.0, -0.25, 0.0, 1.0};
//+
Point(8) = {20.0, -0.25, 0.0, 1.0};
//+
Point(9) = {10.0, 0.25-0.017455, 0.0, 1.0};
//+
Point(10) = {10.0, -0.25-0.017455, 0.0, 1.0};
//+
Line(1) = {1, 2};
//+
Line(2) = {3, 4};
//+
Line(3) = {4, 8};
//+
Line(4) = {8, 7};
//+
Line(5) = {1, 5};
//+
Line(6) = {5, 6};
//+
Point(11) = {0.0, 0.25, 0.5, 1.0};
//+
Point(12) = {9.0, 0.25, 0.5, 1.0};
//+
Point(13) = {11.0, 0.25, 0.5, 1.0};
//+
Point(14) = {20.0, 0.25, 0.5, 1.0};
//+
Point(15) = {0.0, -0.25, 0.5, 1.0};
//+
Point(16) = {9.0, -0.25, 0.5, 1.0};
//+
Point(17) = {11.0, -0.25, 0.5, 1.0};
//+
Point(18) = {20.0, -0.25, 0.5, 1.0};
//+
Point(19) = {10.0, 0.25-0.017455, 0.5, 1.0};
//+
Point(20) = {10.0, -0.25-0.017455, 0.5, 1.0};
//+
Line(9) = {11, 12};
//+
Line(10) = {13, 14};
//+
Line(11) = {14, 18};
//+
Line(12) = {18, 17};
//+
Line(13) = {11, 15};
//+
Line(14) = {15, 16};
//+
Line(17) = {14, 4};
//+
Line(18) = {18, 8};
//+
Line(19) = {15, 5};
//+
Line(20) = {11, 1};
//+
Line(21) = {12, 19};
//+
Line(22) = {19, 13};
//+
Line(23) = {2, 9};
//+
Line(24) = {9, 3};
//+
Line(25) = {16, 20};
//+
Line(26) = {20, 17};
//+
Line(27) = {6, 10};
//+
Line(28) = {10, 7};
//+
Line(29) = {9, 19};
//+
Line(30) = {10, 20};
//+
Physical Curve("fronttop", 3) = {9, 21, 22, 10};
//+
Physical Curve("frontbottom", 4) = {14, 25, 26, 12};
//+
Physical Curve("backbottom", 10) = {6, 27, 28, 4};
//+
Physical Curve("backtop", 9) = {1, 23, 24, 2};
//+
Physical Curve("frontleft", 1) = {13};
//+
Physical Curve("leftbottom", 6) = {19};
//+
Physical Curve("lefttop", 5) = {20};
//+
Physical Curve("backleft", 11) = {5};
//+
Line(31) = {2, 12};
//+
Line(32) = {3, 13};
//+
Line(33) = {6, 16};
//+
Line(34) = {7, 17};
//+
Physical Curve("righttop", 7) = {17};
//+
Physical Curve("rightbottom", 8) = {18};
//+
Physical Curve("backright", 12) = {3};
//+
Physical Curve("frontright", 2) = {11};//+
Curve Loop(1) = {2, -17, -10, -32};
//+
Plane Surface(1) = {1};
//+
Curve Loop(2) = {24, 32, -22, -29};
//+
Plane Surface(2) = {2};
//+
Curve Loop(3) = {23, 29, -21, -31};
//+
Plane Surface(3) = {3};
//+
Curve Loop(4) = {1, 31, -9, 20};
//+
Plane Surface(4) = {4};
//+
Physical Surface("top", 16) = {4, 3, 2, 1};
//+
Curve Loop(5) = {14, -33, -6, -19};
//+
Plane Surface(5) = {5};
//+
Curve Loop(6) = {33, 25, -30, -27};
//+
Plane Surface(6) = {6};
//+
Curve Loop(7) = {28, 34, -26, -30};
//+
Plane Surface(7) = {7};
//+
Curve Loop(8) = {4, 34, -12, 18};
//+
Plane Surface(8) = {8};
//+
Physical Surface("bottom", 17) = {5, 6, 7, 8};
//+
Curve Loop(9) = {9, 21, 22, 10, 11, 12, -26, -25, -14, -13};
//+
Plane Surface(9) = {9};
//+
Physical Surface("front", 13) = {9};
//+
Curve Loop(10) = {2, 3, 4, -28, -27, -6, -5, 1, 23, 24};
//+
Plane Surface(10) = {10};
//+
Physical Surface("back", 18) = {10};
//+
Curve Loop(11) = {20, 5, -19, -13};
//+
Plane Surface(11) = {11};
//+
Physical Surface("left", 14) = {11};
//+
Curve Loop(12) = {18, -3, -17, 11};
//+
Plane Surface(12) = {12};
//+
Physical Surface("right", 15) = {12};
//+
Surface Loop(1) = {4, 10, 1, 12, 8, 7, 9, 3, 2, 6, 5, 11};
//+
Volume(1) = {1};
//+
Physical Volume("volume", 19) = {1};
For meshing, I selected the 2D algorithm to be Frontal-Delaunay for quads (experimental), the 3D algorithm is Delaunay, 2D recombination algorithm is Blossom. I checked the box for Recombine all triangular meshes and selected the subdivision algorithm to be All Hexas.
I am sorry if this is a trivial question. I haven’t worked with hexahedral meshes before. I searched the group but I found nothing that solves this issue. I am not sure if I am creating the mesh incorrectly in Gmsh or if I am doing something wrong while importing it to FEniCS.
Please advise. Thanks in advance!