I have the following script to create mesh from boolean operations on cylinders and bricks.
import sys
import numpy as np
sys.path.append("/opt/Coreform-Cubit-2025.3/bin")
import cubit
last_id = 0
LX = 20
LY = 20
lxs = np.arange(-0.5*LX+2.5, 0.5*LX, 5)
lys = np.arange(-0.5*LY+2.5, 0.5*LY, 5)
radius = 1.5
height = 50
for x in lxs:
for y in lys:
cubit.cmd(f"create cylinder radius {radius} height {height}")
last_id += 1
cubit.cmd(f"volume {last_id} move x {x} y {y} z 50")
cubit.cmd("create brick x 20 y 20 z 5")
last_id += 1
cubit.cmd(f"volume {last_id} move z 77.5")
cubit.cmd("unite all")
cubit.cmd("create brick x 20 y 20 z 80")
last_id = cubit.parse_cubit_list('volume', 'all')[-1]
cubit.cmd(f"volume {last_id} move z 40")
volume_ids = cubit.parse_cubit_list('volume', 'all')
cubit.cmd(f"remove overlap volume {volume_ids[0]} {volume_ids[1]} modify larger")
curves = cubit.parse_cubit_list('curve', 'all')
circle_arcs = []
for curve in curves:
# print(cubit.get_curve_center(curve))
if np.isclose(cubit.get_arc_length(curve), 2 * np.pi * radius):
circle_arcs.append(curve)
cubit.cmd(f"modify curve {' '.join(map(str, circle_arcs))} blend radius 0.025")
volumes = cubit.parse_cubit_list('volume', 'all')
cubit.cmd(f"volume {volumes[0]} name 'vol1'")
cubit.cmd(f"volume {volumes[1]} name 'vol2'")
surfaces = cubit.parse_cubit_list('surface', 'all')
for surf in surfaces:
centroid = cubit.get_surface_centroid(surf)
area = cubit.get_surface_area(surf)
if np.all(np.isclose(centroid[:2], [0, 0])):
normal = cubit.get_surface_normal(surf)
print(surf, "centroid:", centroid, "area:", area, "normal:", normal)
cubit.cmd("vol all scheme tetmesh")
cubit.cmd("mesh volume all")
filename = "mesh.bdf"
cubit.cmd(f'export nastran {filename} overwrite')
Can someone explain why at the intersection two surface ids are provided, and each have slightly different centroids?
`
55 centroid: (-3.120176742220682e-15, 4.019035726780615e-15, 75.00000000000122) area: 283.10133735992366 normal: (-0.0, -0.0, -1.0)
100 centroid: (-5.55674100982099e-16, 1.4934047829731165e-15, 75.00000000000074) area: 283.10133735992366 normal: (0.0, 0.0, 1.0)
`
I understand that the centroids are almost equal given float precision >>> np.finfo(float).eps np.float64(2.220446049250313e-16)
, but I did not expect that two distinct surfaces would be reported. A similar operation in GMSH produces only one distinct surface hence my confusion.