Hello,
I’m using OCC to construct 3D geometries, and I’d like to apply prismatic layers to selected boundaries. While testing this, I’ve noticed an issue when combining boundary renaming with boolean operations. Specifically, if I rename a boundary and use an operator such as + (union), - (difference), or * (intersection), the meshing step fails.
Either operation on its own works fine. I can rename boundaries without boolean operations, or I can use boolean operations without renaming, and the mesh is generated correctly. However, doing both together consistently causes a meshing error.
As a temporary workaround, I’m currently meshing the geometry without renaming the boundaries and then modifying the boundary names after the mesh has been created. I was wondering whether there is a more robust solution to this within the NGSolve toolchain.
I’ve attached a small script below that reproduces the issue.
Thank you in advance for any advice.
from ngsolve import *
from netgen.occ import *
from netgen.meshing import BoundaryLayerParameters
from pathlib import Path
# Configuration flags
replace_bc = True
quarter_geom = False
using_substraction = False
using_intersection = True
L = 1
# Define geometries
box = Box(Pnt(-L, -L, -L), Pnt(L, L, L))
inner_box = Box(Pnt(-0.5*L, -0.5*L, -0.5*L), Pnt(0.5*L, 0.5*L, 0.5*L))
quarter = Box(Pnt(0, 0, -2*L), Pnt(2*L, 2*L, 2*L))
three_quarter = Box(Pnt(-2*L, -2*L, -2*L), Pnt(2*L, 2*L, 2*L)) - quarter
# Set materials and boundary conditions
quarter.mat("quarter")
quarter.bc("cut")
box.mat("outer_box")
box.bc("outer")
inner_box.mat("inner_box")
inner_box.bc("inner_bc1")
if replace_bc:
inner_box.faces.Min(Z).name = "inner_bc2"
inner_box.faces.Max(Z).name = "inner_bc2"
box = box - inner_box
if quarter_geom:
if using_substraction:
box_quarter = box - three_quarter
inner_box_quarter = inner_box - three_quarter
elif using_intersection:
box_quarter = box * quarter
inner_box_quarter = inner_box * quarter
object3D = Glue([box_quarter, inner_box_quarter])
else:
object3D = Glue([box, inner_box])
geometry = OCCGeometry(object3D)
# Boundary layer parameters
B = BoundaryLayerParameters(boundary="inner_bc1", thickness=[0.1], domain="inner_box")
# Generate mesh
ngmesh = geometry.GenerateMesh(boundary_layers=[B])
mesh = Mesh(ngmesh)
print("The mesh has", mesh.ne, "elements")
# Save mesh
file_path = Path().absolute()
mesh_save_path = Path(file_path) / "meshes"
mesh_save_path.mkdir(parents=True, exist_ok=True)
mesh.ngmesh.Save(str(Path(mesh_save_path) / "meshvol.vol"))
print("Mesh exported")
print("Replaced boundary condition names?:", replace_bc)
print("Quarter geometry used?:", quarter_geom)
