GenerateVolumeMesh() error after merging two meshes.

from netgen.occ import * from netgen import meshing from netgen.meshing import Mesh, FaceDescriptor, Element3D,Element2D, Element1D geo1 = OCCGeometry("merge_cube1.step") m1 = geo1.GenerateMesh (maxh=1, minh=0.5, perfstepsend= meshing.MeshingStep.MESHSURFACE) geo2 = OCCGeometry("merge_cube2.step") m2 = geo2.GenerateMesh (maxh=1, minh=0.5, perfstepsend= meshing.MeshingStep.MESHSURFACE) print ("***************************") print ("** merging suface meshes **") print ("***************************") # create an empty mesh mesh = Mesh() # a face-descriptor stores properties associated with a set of surface elements # bc .. boundary condition marker, # domin/domout .. domain-number in front/back of surface elements (0 = void), # surfnr .. number of the surface described by the face-descriptor fd_outside = mesh.Add (FaceDescriptor(bc=1,domin=1,surfnr=1)) fd_inside = mesh.Add (FaceDescriptor(bc=2,domin=2,domout=1,surfnr=2)) # copy all boundary points from first mesh to new mesh. # pmap1 maps point-numbers from old to new mesh pmap1 = { } for e in m1.Elements2D(): for v in e.vertices: if (v not in pmap1): pmap1[v] = mesh.Add (m1[v]) # copy surface elements from first mesh to new mesh # we have to map point-numbers: for e in m1.Elements2D(): mesh.Add(Element2D (fd_outside, [pmap1[v] for v in e.vertices])) # same for the second mesh: pmap2 = { } for e in m2.Elements2D(): for v in e.vertices: if (v not in pmap2): pmap2[v] = mesh.Add (m2[v]) for e in m2.Elements2D(): mesh.Add (Element2D (fd_inside, [pmap2[v] for v in e.vertices])) print ("******************") print ("** merging done **") print ("******************") mesh.GenerateVolumeMesh()