#!/usr/bin/env python3
gmsh_surface_to_3d_periodic_copy.py (edge-seeding + fiber-profile capture)
STEP → (fragment + heal + OPTIONAL equal partitions) → surface mesh
→ declare periodic min->max faces (±X/±Y/±Z) on dims 2/1/0 (identical nodes)
→ generate 3D tets (boundary respected) → delete 2D/1D elements
→ export orphan INP with only NODE + C3D elements
(optional) export partitioned STEP for downstream tools
import sys, os, glob, math, gmsh
STEP_PATH_DEFAULT = r"C:\temp1"
OUT_DIR_DEFAULT = r"C:\temp1"
OUT_BASENAME = “Composite Part”
Default number of equal parts along X, Y, Z (same on all axes).
EQUAL_PARTS_DEFAULT = 5
-------------------- small utils --------------------
def _set_periodic_faces(master_face, slave_face, dxyz):
dx, dy, dz = dxyz
T = [ 1,0,0,dx,
0,1,0,dy,
0,0,1,dz,
0,0,0, 1 ]
gmsh.model.mesh.setPeriodic(2, [slave_face], [master_face], T)
def _safe_set(name, value):
try:
gmsh.option.setNumber(name, value)
except Exception:
pass
def _split_paths(spec):
for sep in (';', ',', '|'):
if sep in spec:
return [s.strip().strip('"') for s in spec.split(sep) if s.strip()]
return [spec]
def _gather_step_inputs(step_path, steps_flag):
c=[]
if steps_flag:
for p in _split_paths(steps_flag):
if os.path.isdir(p):
c += sorted(glob.glob(os.path.join(p, "*.stp"))+
glob.glob(os.path.join(p, "*.step")))
else:
c.append(p)
else:
if os.path.isdir(step_path):
c = sorted(glob.glob(os.path.join(step_path, "*.stp"))+
glob.glob(os.path.join(step_path, "*.step")))
else:
c = [step_path]
seen, files=set(),[]
for p in c:
q=os.path.abspath(p)
if os.path.isfile(q) and q not in seen:
seen.add(q); files.append(q)
return files
def _parse_cli(argv):
step_path = STEP_PATH_DEFAULT
target_size = 0.3
out_dir = OUT_DIR_DEFAULT
out_base = OUT_BASENAME
periodic = "xyz" # e.g. 'x', 'yz', 'xyz'
steps_flag = None
keep3dopt = 0 # 0 strict boundary; 1 allow interior optimize
export_step = 0 # export partitioned STEP
parts_equal = EQUAL_PARTS_DEFAULT
pos, flags = [], {}
for a in argv[1:]:
if a.startswith("--"):
k,_,v=a.partition("="); flags[k]=v
else:
pos.append(a)
if len(pos)>=1: step_path = pos[0]
if len(pos)>=2: target_size = float(pos[1])
if "--steps" in flags and flags["--steps"]: steps_flag = flags["--steps"]
if "--outdir" in flags and flags["--outdir"]: out_dir = flags["--outdir"]
if "--basename" in flags and flags["--basename"]: out_base = flags["--basename"]
if "--periodic" in flags and flags["--periodic"]: periodic = flags["--periodic"].lower()
if "--keep3dopt" in flags and flags["--keep3dopt"]: keep3dopt = int(flags["--keep3dopt"])
if "--exportSTEP" in flags and flags["--exportSTEP"]: export_step = int(flags["--exportSTEP"])
if "--parts" in flags and flags["--parts"]:
try:
parts_equal = max(1, int(flags["--parts"]))
except Exception:
parts_equal = EQUAL_PARTS_DEFAULT
# NEW flags
return (step_path, target_size, out_dir, out_base, periodic, steps_flag,
keep3dopt, export_step, parts_equal,
)
def _clean_inp_to_orphan_mesh(inp_path,
keep_element_prefixes=("C3D",),
drop_keywords=("NSET","ELSET","SURFACE")):
try:
lines=open(inp_path,"r").readlines()
except IOError:
return
def _is_hdr(s): return s.lstrip().startswith("*")
def _HDR(s): return s.lstrip().upper()
drop_hdrs = tuple("*"+k.upper() for k in drop_keywords)
keep_pref = tuple(p.upper() for p in keep_element_prefixes) if keep_element_prefixes else None
out=[]; i=0
while i<len(lines):
ln = lines[i]
if _is_hdr(ln):
H=_HDR(ln)
if H.startswith(drop_hdrs):
i+=1
while i<len(lines) and not _is_hdr(lines[i]): i+=1
continue
if H.startswith("*ELEMENT") and keep_pref is not None:
header = H.replace(" ",""); et=None
for tok in header.split(","):
if tok.startswith("TYPE="): et = tok[5:]; break
if et and any(et.startswith(p) for p in keep_pref):
out.append(lines[i]); i+=1
while i<len(lines) and not _is_hdr(lines[i]):
out.append(lines[i]); i+=1
else:
i+=1
while i<len(lines) and not _is_hdr(lines[i]): i+=1
continue
out.append(ln); i+=1
if H.startswith("*NODE"):
while i<len(lines) and not _is_hdr(lines[i]):
out.append(lines[i]); i+=1
continue
out.append(ln); i+=1
try:
open(inp_path,"w").writelines(out)
except IOError:
pass
---------- geometry helpers ----------
def _bb(dim, tag):
x0,y0,z0,x1,y1,z1 = gmsh.model.getBoundingBox(dim, tag)
return (0.5*(x0+x1), 0.5*(y0+y1), 0.5*(z0+z1)), (x1-x0, y1-y0, z1-z0)
def _edges_of_face(face_tag):
return [t for (d,t) in gmsh.model.getBoundary([(2,face_tag)], oriented=False, recursive=False) if d==1]
def _face_pairs_by_axis(axis, bbox, h):
x0,y0,z0,x1,y1,z1 = bbox
gmin = (x0,y0,z0)[axis]; gmax = (x1,y1,z1)[axis]
span_tol = max(0.25*h, 1e-9)
pos_tol = max(0.10*h, 1e-9)
faces = gmsh.model.getEntities(2)
mins, maxs = [], []
for _, s in faces:
(cx,cy,cz), (sx,sy,sz) = _bb(2, s)
if (sx,sy,sz)[axis] > span_tol:
continue
if abs((cx,cy,cz)[axis] - gmin) <= pos_tol: mins.append(s)
elif abs((cx,cy,cz)[axis] - gmax) <= pos_tol: maxs.append(s)
if axis==0: key=lambda s: (_bb(2,s)[0][1], _bb(2,s)[0][2])
elif axis==1: key=lambda s: (_bb(2,s)[0][0], _bb(2,s)[0][2])
else: key=lambda s: (_bb(2,s)[0][0], _bb(2,s)[0][1])
mins.sort(key=key); maxs.sort(key=key)
return mins, maxs
---------- periodic helpers ----------
def _sorted_edges_of_face(face_tag, axis):
edges = _edges_of_face(face_tag)
def _key(e):
(cx,cy,cz), _ = _bb(1, e)
if axis == 0: return (cy, cz)
if axis == 1: return (cx, cz)
return (cx, cy)
return sorted(edges, key=_key)
def _points_of_edge(edge_tag):
bnd = gmsh.model.getBoundary([(1, edge_tag)], oriented=False, recursive=False)
return [t for (d,t) in bnd if d == 0]
def _set_periodic_all(master_face, slave_face, dxyz, axis):
_set_periodic_faces(master_face, slave_face, dxyz)
dx,dy,dz = dxyz
T = [1,0,0,dx, 0,1,0,dy, 0,0,1,dz, 0,0,0,1]
m_edges = _sorted_edges_of_face(master_face, axis)
s_edges = _sorted_edges_of_face(slave_face, axis)
if len(m_edges) != len(s_edges):
raise RuntimeError("Edge count mismatch on faces %d vs %d" % (master_face, slave_face))
for me, se in zip(m_edges, s_edges):
gmsh.model.mesh.setPeriodic(1, [se], [me], T)
m_pts = sorted(_points_of_edge(me))
s_pts = sorted(_points_of_edge(se))
n = min(len(m_pts), len(s_pts))
for mp, sp in zip(m_pts[:n], s_pts[:n]):
gmsh.model.mesh.setPeriodic(0, [sp], [mp], T)
def _copy_all_opposite_faces(axes=‘xyz’):
x0,y0,z0,x1,y1,z1 = gmsh.model.getBoundingBox(-1,-1)
bbox = (x0,y0,z0,x1,y1,z1)
Lx,Ly,Lz = x1-x0, y1-y0, z1-z0
base = min(max(Lx,1e-12), max(Ly,1e-12), max(Lz,1e-12))
h = max(base/25.0, 1e-4)
gmsh.model.occ.removeAllDuplicates()
gmsh.model.occ.synchronize()
axes = axes.lower().replace('all','xyz')
axes = ''.join(ch for ch in axes if ch in 'xyz')
total = 0
for ax, dxyz in zip((0,1,2), [(Lx,0,0),(0,Ly,0),(0,0,Lz)]):
if 'xyz'[ax] not in axes:
continue
mins, maxs = _face_pairs_by_axis(ax, bbox, h)
if len(mins) != len(maxs):
raise RuntimeError("Opposite face count mismatch on axis %s: %d vs %d"
% ('xyz'[ax], len(mins), len(maxs)))
for mFace, sFace in zip(mins, maxs):
_set_periodic_all(mFace, sFace, dxyz, axis=ax)
total += 1
return total
def _prune_tiny_volumes(min_span=None):
x0,y0,z0,x1,y1,z1 = gmsh.model.getBoundingBox(-1,-1)
Lx,Ly,Lz = max(x1-x0,1e-12), max(y1-y0,1e-12), max(z1-z0,1e-12)
if min_span is None:
min_span = 1e-4 * min(Lx, Ly, Lz)
tiny = []
for _, v in gmsh.model.getEntities(3):
a,b,c,d,e,f = gmsh.model.getBoundingBox(3, v)
if min(d-a, e-b, f-c) < min_span:
tiny.append((3, v))
if tiny:
gmsh.model.occ.remove(tiny)
gmsh.model.occ.synchronize()
print("[PRUNE] Removed %d tiny volumes (min_span=%.3g)" % (len(tiny), min_span))
def _remove_surface_and_edge_elements():
for dim in (2,1):
for _, tag in gmsh.model.getEntities(dim):
eTypes, eTags, _ = gmsh.model.mesh.getElements(dim, tag)
if eTypes:
allET = [t for arr in eTags for t in arr]
if allET:
gmsh.model.mesh.removeElements(dim, tag, allET)
---------------- Equal partitions + fragment ----------------
def _equal_internal_planes(x0, x1, parts):
if parts is None or parts <= 1:
return []
L = max(x1 - x0, 1e-12)
return [x0 + (i * L) / float(parts) for i in range(1, parts)]
def _add_equal_partitions_and_fragment(bbox, parts_equal):
if parts_equal is None or parts_equal <= 1:
return 0
x0,y0,z0,x1,y1,z1 = bbox
Lx, Ly, Lz = x1-x0, y1-y0, z1-z0
pad = 1e-6 + 0.02*max(Lx,Ly,Lz)
xs = _equal_internal_planes(x0, x1, parts_equal)
ys = _equal_internal_planes(y0, y1, parts_equal)
zs = _equal_internal_planes(z0, z1, parts_equal)
epsX = 1e-6 * max(Lx, 1.0)
epsY = 1e-6 * max(Ly, 1.0)
epsZ = 1e-6 * max(Lz, 1.0)
xs = [u for u in xs if (x0 + epsX) < u < (x1 - epsX)]
ys = [u for u in ys if (y0 + epsY) < u < (y1 - epsY)]
zs = [u for u in zs if (z0 + epsZ) < u < (z1 - epsZ)]
if not (xs or ys or zs):
return 0
def _rectYZ_at_x(xc):
p1 = gmsh.model.occ.addPoint(xc, y0-pad, z0-pad)
p2 = gmsh.model.occ.addPoint(xc, y1+pad, z0-pad)
p3 = gmsh.model.occ.addPoint(xc, y1+pad, z1+pad)
p4 = gmsh.model.occ.addPoint(xc, y0-pad, z1+pad)
l1 = gmsh.model.occ.addLine(p1,p2); l2 = gmsh.model.occ.addLine(p2,p3)
l3 = gmsh.model.occ.addLine(p3,p4); l4 = gmsh.model.occ.addLine(p4,p1)
cl = gmsh.model.occ.addCurveLoop([l1,l2,l3,l4])
return gmsh.model.occ.addPlaneSurface([cl])
def _rectXZ_at_y(yc):
p1 = gmsh.model.occ.addPoint(x0-pad, yc, z0-pad)
p2 = gmsh.model.occ.addPoint(x1+pad, yc, z0-pad)
p3 = gmsh.model.occ.addPoint(x1+pad, yc, z1+pad)
p4 = gmsh.model.occ.addPoint(x0-pad, yc, z1+pad)
l1 = gmsh.model.occ.addLine(p1,p2); l2 = gmsh.model.occ.addLine(p2,p3)
l3 = gmsh.model.occ.addLine(p3,p4); l4 = gmsh.model.occ.addLine(p4,p1)
cl = gmsh.model.occ.addCurveLoop([l1,l2,l3,l4])
return gmsh.model.occ.addPlaneSurface([cl])
def _rectXY_at_z(zc):
p1 = gmsh.model.occ.addPoint(x0-pad, y0-pad, zc)
p2 = gmsh.model.occ.addPoint(x1+pad, y0-pad, zc)
p3 = gmsh.model.occ.addPoint(x1+pad, y1+pad, zc)
p4 = gmsh.model.occ.addPoint(x0-pad, y1+pad, zc)
l1 = gmsh.model.occ.addLine(p1,p2); l2 = gmsh.model.occ.addLine(p2,p3)
l3 = gmsh.model.occ.addLine(p3,p4); l4 = gmsh.model.occ.addLine(p4,p1)
cl = gmsh.model.occ.addCurveLoop([l1,l2,l3,l4])
return gmsh.model.occ.addPlaneSurface([cl])
cut_surfs=[]
for xc in xs: cut_surfs.append((2,_rectYZ_at_x(xc)))
for yc in ys: cut_surfs.append((2,_rectXZ_at_y(yc)))
for zc in zs: cut_surfs.append((2,_rectXY_at_z(zc)))
gmsh.model.occ.synchronize()
vols = [(3,v) for (_,v) in gmsh.model.getEntities(3)]
try:
out,_ = gmsh.model.occ.fragment(vols, cut_surfs)
cnt = len([1 for d,_ in out if d==3])
finally:
gmsh.model.occ.synchronize()
try:
gmsh.model.occ.removeAllDuplicates()
finally:
gmsh.model.occ.synchronize()
_prune_tiny_volumes()
return cnt
------------------------------ main ------------------------------
def main():
(step_path, target_size, out_dir, out_base, periodic, steps_flag,
keep3dopt, export_step, parts_equal,
) = _parse_cli(sys.argv)
step_files = _gather_step_inputs(step_path, steps_flag)
if not step_files: raise RuntimeError("No STEP files found.")
if not os.path.isdir(out_dir): os.makedirs(out_dir)
OUT_MSH = os.path.join(out_dir, out_base + ".msh")
OUT_INP = os.path.join(out_dir, out_base + ".inp")
OUT_STP = os.path.join(out_dir, out_base + "_partitioned.step")
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 1)
gmsh.model.add("STEPModel")
gmsh.logger.start()
print("[IMPORT] STEP files:")
for p in step_files:
print(" -", p)
gmsh.model.occ.importShapes(p)
gmsh.model.occ.synchronize()
vols = gmsh.model.getEntities(3)
if not vols: raise RuntimeError("No 3D volumes in geometry.")
# imprint touching faces + heal
try:
gmsh.model.occ.fragment(vols, [])
finally:
try:
gmsh.model.occ.removeAllDuplicates()
finally:
gmsh.model.occ.synchronize()
# bbox + default h
xmin,ymin,zmin,xmax,ymax,zmax = gmsh.model.getBoundingBox(-1,-1)
Lx,Ly,Lz = xmax-xmin, ymax-ymin, zmax-zmin
if target_size is None:
base = min(max(Lx,1e-12), max(Ly,1e-12), max(Lz,1e-12))
target_size = max(base/25.0, 1e-4)
# OPTIONAL: add equal partitions (for conformal sub-regions)
nvol_after = _add_equal_partitions_and_fragment((xmin,ymin,zmin,xmax,ymax,zmax), parts_equal)
if nvol_after:
print("[PARTITION] Equal parts/axis:", parts_equal, "| total 3D entities now:", nvol_after)
# ---------- 2D boundary settings (deterministic) ----------
pts = gmsh.model.getEntities(0)
if pts: gmsh.model.mesh.setSize(pts, target_size)
# Turn ON curvature-based sizing (for round fibers) and enforce elems per 2*pi
_safe_set("Mesh.CharacteristicLengthFromCurvature", 1)
_safe_set("Mesh.MinimumElementsPerTwoPi", 16) # default 24; bump with --m2pi
_safe_set("Mesh.Algorithm", 5) # 2D Delaunay
_safe_set("Mesh.Optimize", 0)
_safe_set("Mesh.Smoothing", 0)
_safe_set("Mesh.RecombineAll", 0)
# ---------- Declare periodicity (faces+edges+points) then generate 2D ----------
if periodic:
try:
n_pairs = _copy_all_opposite_faces(periodic)
if n_pairs == 0:
print("[WARN] No opposite planar face pairs found to declare periodicity.")
except RuntimeError as e:
print("[PERIODIC ERROR]", str(e))
# fall back to non-periodic if pairing fails
# Optional background field to make elements even smaller near *all* curves
gmsh.model.mesh.generate(2)
# ---------- 3D fill (Tet4) ----------
_safe_set("Mesh.ElementOrder", 1)
_safe_set("Mesh.SecondOrderLinear", 0)
_safe_set("Mesh.Algorithm3D", 4)
if keep3dopt:
_safe_set("Mesh.Optimize", 0)
_safe_set("Mesh.OptimizeNetgen", 0)
else:
_safe_set("Mesh.Optimize", 0)
_safe_set("Mesh.OptimizeNetgen", 0)
gmsh.model.mesh.generate(3)
# remove 2D/1D so only 3D solids remain
_remove_surface_and_edge_elements()
# Export meshes
_safe_set("Mesh.SaveAll", 0)
gmsh.option.setNumber("Mesh.SaveGroupsOfElements", 0)
gmsh.option.setNumber("Mesh.SaveGroupsOfNodes", 0)
gmsh.write(OUT_MSH)
gmsh.write(OUT_INP)
# keep only *NODE + *ELEMENT (C3D*)
_clean_inp_to_orphan_mesh(OUT_INP, keep_element_prefixes=("C3D",),
drop_keywords=("NSET","ELSET","SURFACE"))
# (optional) export the partitioned CAD
if export_step:
try:
gmsh.write(OUT_STP)
print("[STEP] Exported partitioned STEP:", OUT_STP)
except Exception:
print("[STEP] Export failed (writer not available for this geometry).")
print("\n==== Summary ====")
print("STEPs:", len(step_files))
print("h : {:.6g}".format(target_size))
print("BBox : Lx={:.6g}, Ly={:.6g}, Lz={:.6g}".format(Lx, Ly, Lz))
print("Equal parts/axis:", parts_equal)
print("Copied axes (periodic):", periodic if periodic else "none")
print("3D optimize:", "ON" if keep3dopt else "OFF (strict)")
print("Wrote:\n - {}\n - {}".format(OUT_MSH, OUT_INP))
for msg in gmsh.logger.get():
if ("volumes left to mesh" in msg) or ("Error" in msg) or ("Warning" in msg):
print("[GMSH]", msg)
gmsh.logger.stop()
gmsh.finalize()
if name == “main”:
main()
