SubMesh on 3D hexahedral mesh

mesh
1

Hi,

I am trying use the SubMesh function on a 3D hexahedral mesh, but I get the following error: Unable to order hexahedron cell. *** Reason: Cell is not orderable.

Does anyone know a workaround for this error? I know SubMesh works for tetrahedral meshes but I have some other functions in my code which relies on using a hexmesh which is why I am trying to debug this.

MWE:

from dolfin import *

#Define geometry
L = 10; W = L; H = 5
Nx = 30; Ny = Nx; Nz = 15

mesh = BoxMesh.create(MPI.comm_world,[Point(0, 0, 0), Point(L, H, W)],[Nx, Ny, Nz],CellType.Type.hexahedron)

#Select subdomain
class SUBDOMAIN(SubDomain):
    def inside(self, x, on_boundary):
        return (((1-0.1)*H) - x[1] >= -1e-16 and ((0.1)*H) - x[1] <= -1e-16)


OptDom = OptimizationDomain()

domains = MeshFunction('size_t', mesh,3)
domains.set_all(0)
OptDom.mark(domains, 1)
#Fails here
OptMesh = SubMesh(mesh, domains, 1)

I have tried to create the mesh in GMSH and split the domain, however I can’t seem to import hexahedral meshes into FEniCS.

I am using DOLFIN version 2019.1.0 on MacOS 10.15

Thanks in advance,
Kiri