Hi All,
I am trying to load the results of a previous simulation using hdf5. I am able to load the solution and the mesh but the boundary labels are not working. There is a cell mismatch error (see below). I have used this method successfully before but this was 2D and possibly a few releases ago. The current simulation is 3D and uses hexahedral elements.
I am confused as to what is causing the dimension mismatch. When I check they are the same.
Kind Regards,
Nick
*** -------------------------------------------------------------------------
*** Error: Unable to read meshfunction topology.
*** Reason: Cell dimension mismatch.
*** Where: This error was encountered inside HDF5File.cpp.
*** Process: 0
*** DOLFIN version: 2018.1.0
*** Git changeset: 948dc42cc4e06ed9227d0201ad50f94ac94cbf9f
*** -------------------------------------------------------------------------
#Test code copied below.
from dolfin import *
from mshr import *
class Top(SubDomain):
def inside(self,x,on_boundary):
return near(x[2],1)
class Bottom(SubDomain):
def inside(self,x,on_boundary):
return near(x[2],0)
class LHS(SubDomain):
def inside(self,x,on_boundary):
return near(x[0],0.0) and between(x[1],(0.0,1)) and between(x[2],(0.0,1))
class RHS(SubDomain):
def inside(self,x,on_boundary):
return near(x[0],1) and between(x[1],(0.0,1)) and between(x[2],(0.0,1))
class FRONT(SubDomain):
def inside(self,x,on_boundary):
return between(x[0],(0.0,1)) and near(x[1],0.0) and between(x[2],(0.0,1))
class BACK(SubDomain):
def inside(self,x,on_boundary):
return between(x[0],(0.0,1)) and near(x[1],1) and between(x[2],(0.0,1))
#Output problem setup
mesh = UnitCubeMesh.create(10,5,5,dolfin.cpp.mesh.CellType.Type(5))
boundaries = MeshFunction(“size_t”,mesh,mesh.topology().dim()-1)
top = Top()
bottom = Bottom()
lhs = LHS()
rhs = RHS()
front = FRONT()
back = BACK()
boundaries.set_all(0)
top.mark(boundaries,1)
bottom.mark(boundaries,2)
lhs.mark(boundaries,3)
rhs.mark(boundaries,3)
front.mark(boundaries,4)
back.mark(boundaries,4)
cell = mesh.ufl_cell()
VE_u = VectorElement(“CG”,cell,2)
VE_g = VectorElement(“CG”,cell,1)
FE_p = FiniteElement(“CG”,cell,1)
element = MixedElement([VE_u,VE_g,FE_p])
V = FunctionSpace(mesh, element)
ugp1 = Function(V)
output_file = HDF5File(mesh.mpi_comm(),“Test.h5”, “w”)
output_file.write(ugp1, “solution”)
output_file.write(mesh, “mesh”)
output_file.write(boundaries, “boundaries”)
output_file.close()
##########################
#Input problem setup
mesh2 = Mesh()
InputFileLoc = “Test”
input_file = HDF5File(MPI.comm_world,InputFileLoc+".h5", “r”)
input_file.read(mesh2, “mesh”,False)
cell2 = mesh2.ufl_cell()
VE_u2 = VectorElement(“CG”,cell2,2)
VE_g2 = VectorElement(“CG”,cell2,1)
FE_p2 = FiniteElement(“CG”,cell2,1)
element2 = MixedElement([VE_u2,VE_g2,FE_p2])
V = FunctionSpace(mesh2, element)
input_file.read(ugp1, “solution”)
ugp0 = Function(V)
ugp0.assign(ugp1)
#Error occurs here when the new boundary labels are updated
boundaries2 = MeshFunction(“size_t”,mesh2,mesh.topology().dim()-1)
input_file.read(boundaries2,“boundaries”)