Combine two extruded volumes where one has a hole

Question

I have tried to combine two volumes where one has a cylindrical hole. Both volumes were extruded from a surface. I'm using pygmsh and the code looks like this:

import pygmsh

L = 3
W = 3
c1 = [W/2, L/2, 0]
r = 0.075

resolution1 = 0.11 
resolution2 = 0.3 

geom = pygmsh.occ.Geometry()

model = geom.__enter__()

circle1 = model.add_disk(c1, r, mesh_size=resolution1)

points = [model.add_point((0, 0, 0), mesh_size=resolution2),
          model.add_point((W, 0, 0), mesh_size=resolution2),
          model.add_point((W, L, 0), mesh_size=resolution2),
          model.add_point((0, L, 0), mesh_size=resolution2)]

points2 = [model.add_point((0, 0, 2), mesh_size=resolution2),
          model.add_point((W, 0, 2), mesh_size=resolution2),
          model.add_point((W, L, 2), mesh_size=resolution2),
          model.add_point((0, L, 2), mesh_size=resolution2)]

channel_lines = [model.add_line(points[i], points[i+1])
                 for i in range(-1, len(points)-1)]

channel_lines2 = [model.add_line(points2[i], points2[i+1])
                 for i in range(-1, len(points2)-1)]

channel_loop = model.add_curve_loop(channel_lines)

channel_loop2 = model.add_curve_loop(channel_lines2)

plane_surface = model.add_plane_surface(channel_loop) 

plane_surface1 = model.boolean_difference(plane_surface, circle1)
plane_surface2 = model.add_plane_surface(channel_loop2)

s = model.extrude(plane_surface1, [0.0, 0.0, 2.0], num_layers=2)
t = model.extrude(plane_surface2, [0.0, 0.0, 2.0], num_layers=2)

model.boolean_fragments(s[1], t[1])

mesh = model.generate_mesh()

However I get a message that some node wasn't found:

Exception                                 Traceback (most recent call last)
Cell In[11], line 1
----> 1 mesh = model.generate_mesh()

File ~\env\lib\site-packages\pygmsh\common\geometry.py:374, in CommonGeometry.generate_mesh(self, dim, order, algorithm, verbose)
    371 if algorithm:
    372     gmsh.option.setNumber("Mesh.Algorithm", algorithm)
--> 374 gmsh.model.mesh.generate(dim)
    376 # setOrder() after generate(), see
    377 # <https://github.com/nschloe/pygmsh/issues/515#issuecomment-1020106499>
    378 if order is not None:

File ~\env\lib\site-packages\gmsh.py:2006, in model.mesh.generate(dim)
   2002 lib.gmshModelMeshGenerate(
   2003     c_int(dim),
   2004     byref(ierr))
   2005 if ierr.value != 0:
-> 2006     raise Exception(logger.getLastError())

Exception: Could not find extruded node (1.483310929953277, 1.573119593413637, 4) in surface 13

Can someone help me understand why it does not work? I do not get any problems if I seperate the two volumes, for example by putting

points2 = [model.add_point((0, 0, 3), mesh_size=resolution2),
          model.add_point((W, 0, 3), mesh_size=resolution2),
          model.add_point((W, L, 3), mesh_size=resolution2),
          model.add_point((0, L, 3), mesh_size=resolution2)]

However I need one mesh and not two seperate ones.

Answer 1

Adding a volume (box) in between the two extruded meshes and then combining them with fragments did the job for me.

A working code would be the following:

import pygmsh
import numpy as np

L = 3
W = 3
H = 3
H2 = 2
H3 = 1
c1 = np.array([W/2, L/2, H2+H3])

r = 0.075

resolution1 = 0.1 
resolution2 = 0.5  

geom = pygmsh.occ.Geometry()

model = geom.__enter__()  

points = [model.add_point((0, 0, H2+H3), mesh_size=resolution2),
          model.add_point((W, 0, H2+H3), mesh_size=resolution2),
          model.add_point((W, L, H2+H3), mesh_size=resolution2),
          model.add_point((0, L, H2+H3), mesh_size=resolution2)]

channel_lines = [model.add_line(points[i], points[i+1])
                 for i in range(-1, len(points)-1)]

channel_loop = model.add_curve_loop(channel_lines)

point1 = model.add_point(c1+np.array([r,0,0]), mesh_size=resolution1)
point3 = model.add_point(c1+np.array([-r,0,0]), mesh_size=resolution1)
point5 = model.add_point(c1, mesh_size=resolution1)

arc1 = model.add_circle_arc(point1,point5,point3)
arc2 = model.add_circle_arc(point3,point5,point1)

circ_loop = model.add_curve_loop([arc1,arc2])

rectangle = model.add_plane_surface(channel_loop, holes = [circ_loop]) 

points2 = [model.add_point((0, 0, 0), mesh_size=resolution2),
          model.add_point((W, 0, 0), mesh_size=resolution2),
          model.add_point((W, L, 0), mesh_size=resolution2),
          model.add_point((0, L, 0), mesh_size=resolution2)]

channel_lines2 = [model.add_line(points2[i], points2[i+1])
                 for i in range(-1, len(points2)-1)]

channel_loop2 = model.add_curve_loop(channel_lines2)

rectangle2 = model.add_plane_surface(channel_loop2)

cuboid = model.extrude(rectangle, [0, 0, H], num_layers = H)
cuboid2 = model.extrude(rectangle2, [0, 0, H2], num_layers = H2)
cuboid3 = model.add_box([0, 0, H2], [W, L, H3])

vol = model.boolean_fragments(cuboid[1],cuboid3)
vol2 = model.boolean_fragments(cuboid2[1], vol)

mesh = geom.generate_mesh()