from __future__ import annotations
import importlib.util
import inspect
import operator
import tempfile
import warnings
from contextlib import contextmanager
from pathlib import Path
from typing import Iterable, Literal
import cadquery as cq
import gmsh
import numpy as np
from cadquery import importers
from cadquery.occ_impl.geom import Vector
from cadquery.occ_impl.importers.assembly import importStep as importStepAssembly
from cadquery.occ_impl.shapes import setThreads
from OCP.OSD import OSD_Parallel, OSD_ThreadPool
from OCP.TopoDS import TopoDS_Iterator
from cad_to_dagmc import __version__
# cad-to-dagmc-mesher surface tolerances when the caller sets none: the chordal
# deflection in cm (after scale_factor) and the angular deflection in
# radians. On large imported B-spline models they kept the worst solid's volume
# within 0.2% and lost no particles, while an angular tolerance of 0.5 let a
# thin wall's two faces cross.
DEFAULT_TOLERANCE = 0.5
DEFAULT_ANGULAR_TOLERANCE = 0.2
[docs]
class PyMoabNotFoundError(ImportError):
"""Raised when pymoab is not installed but the pymoab backend is requested."""
def __init__(self, message=None):
if message is None:
message = (
"pymoab is not installed. pymoab/MOAB is not available on PyPI so it "
"cannot be included as a dependency of cad-to-dagmc.\n\n"
"You can install pymoab via one of these methods:\n"
" 1. From conda-forge: conda install -c conda-forge moab\n"
" 2. From extra index: pip install --extra-index-url https://shimwell.github.io/wheels moab\n"
" 3. From source: https://bitbucket.org/fathomteam/moab\n\n"
"Alternatively, use the h5py backend (the default) which does not require pymoab:\n"
" export_dagmc_h5m_file(..., h5m_backend='h5py')"
)
super().__init__(message)
[docs]
class CadToDagmcMesherNotFoundError(ImportError):
"""Raised when cad-to-dagmc-mesher is not installed but its backend is requested."""
def __init__(self, message=None):
if message is None:
message = (
"cad-to-dagmc-mesher is not installed. It is not available on "
"conda-forge so it cannot be included as a dependency of the "
"cad-to-dagmc conda package.\n\n"
"Install it with pip, which works alongside a conda installation:\n"
" pip install cad-to-dagmc-mesher\n\n"
"Alternatively, use the gmsh meshing backend, which is always "
"available:\n"
" export_dagmc_h5m_file(..., meshing_backend='gmsh')"
)
super().__init__(message)
def _cad_to_dagmc_mesher_is_available() -> bool:
"""Return True when the cad-to-dagmc-mesher package can be imported."""
return importlib.util.find_spec("cad_to_dagmc_mesher") is not None
def write_vtk(filename, vertices, tetrahedra):
"""Write a tetrahedral mesh to an ASCII VTK legacy file.
The output is a pure tetrahedron UNSTRUCTURED_GRID in the same legacy
format that gmsh writes today, so it can be read back with
openmc.UnstructuredMesh(filename, library="moab"). The MOAB reader does
not require the GLOBAL_ID POINT_DATA/CELL_DATA blocks that MOAB itself
writes when it exports a mesh, so they are intentionally omitted. This
was confirmed with a round trip transport test (see
tests/test_write_vtk.py::test_write_vtk_openmc_moab_round_trip): a mesh
written without GLOBAL_ID loads in MOAB and tallies identically to one
written with it.
Args:
filename: Output file path.
vertices: Sequence of [x, y, z] coordinates (list or numpy array).
tetrahedra: Sequence of [v0, v1, v2, v3] zero-based vertex indices
(list or numpy array).
"""
n_tets = len(tetrahedra)
# Stream the point/cell blocks with writelines() over generators. This
# keeps memory bounded (nothing bigger than one line is materialised at a
# time, matching the old per-line writes) while letting the C-level
# writelines do the looping, which matters for the large meshes the mesher
# can produce.
with open(filename, "w") as f:
f.write("# vtk DataFile Version 2.0\n")
f.write("Unstructured mesh\n")
f.write("ASCII\n")
f.write("DATASET UNSTRUCTURED_GRID\n")
f.write(f"POINTS {len(vertices)} double\n")
f.writelines(f"{v[0]} {v[1]} {v[2]}\n" for v in vertices)
f.write(f"CELLS {n_tets} {n_tets * 5}\n")
f.writelines(f"4 {t[0]} {t[1]} {t[2]} {t[3]}\n" for t in tetrahedra)
f.write(f"CELL_TYPES {n_tets}\n")
f.writelines("10\n" for _ in range(n_tets))
def combine_tet_meshes(tet_data):
"""Combine per-solid tetrahedral meshes into a single mesh.
``cad_to_dagmc_mesher.cad.mesh_assembly`` returns a ``tet_data`` dict
mapping ``solid_id`` to ``{"vertices": (n, 3) array, "tetrahedra":
(m, 4) array, ...}`` where each solid's tetrahedra index into that
solid's own vertex list. To write a single unstructured grid the vertex
arrays are concatenated and each solid's tetrahedra are offset by the
running vertex count so they index into the combined vertex array.
Solids that were imprinted against each other are meshed from the same
surface mesh on the face they share, so each of them carries its own copy
of the vertices on that face. Concatenating the arrays leaves those copies
as separate points and the result is a mesh that looks like one body but is
really several that only touch, which matters to anything solving on it. A
finite element solve on such a mesh treats the shared face as an exterior
boundary on both sides and nothing crosses it.
The duplicate coordinates are bit for bit identical, so they are merged
here on an exact match. Nothing is merged on proximity, which means solids
that do not touch keep their own vertices and stay as separate bodies, as
they should.
Args:
tet_data: Mapping of solid_id -> dict with "vertices" and
"tetrahedra" entries, as returned by mesh_assembly.
Returns:
(vertices, tetrahedra): a single (N, 3) float array of vertex
coordinates with no exact duplicates, and a single (M, 4) int array of
zero-based tetrahedron vertex indices.
"""
all_vertices = []
all_tetrahedra = []
offset = 0
for solid_id in tet_data:
verts = np.asarray(tet_data[solid_id]["vertices"], dtype=float).reshape(-1, 3)
tets = np.asarray(tet_data[solid_id]["tetrahedra"], dtype=np.int64).reshape(-1, 4)
all_vertices.append(verts)
all_tetrahedra.append(tets + offset)
offset += len(verts)
if not all_vertices:
return np.empty((0, 3), dtype=float), np.empty((0, 4), dtype=np.int64)
vertices = np.vstack(all_vertices)
tetrahedra = np.vstack(all_tetrahedra)
# first_index maps each kept vertex back to where it first appeared, so the
# vertices stay in the order the solids were meshed in rather than being
# sorted, and remap turns the old indices into the new ones
_, first_index, inverse = np.unique(
vertices, axis=0, return_index=True, return_inverse=True
)
keep = np.sort(first_index)
remap = np.empty(len(keep), dtype=np.int64)
remap[np.argsort(first_index)] = np.arange(len(keep))
return vertices[keep], remap[inverse.reshape(-1)][tetrahedra]
def resolve_imprint(
imprint: Literal["partial", "full", False] | None = "partial",
imprint_threads: int | Literal["all"] = "all",
) -> int:
"""Validate the mode and thread limit, including when imprinting is disabled.
Return the thread count. ``"all"`` actively requests the
logical processor count, not the caller's current OCC pool size. The
caller applies the limit only while imprinting and restores it afterwards.
"""
if imprint is not False and imprint not in ("partial", "full", None):
raise ValueError(
"imprint must be 'partial', 'full', None or False, "
f"got {imprint!r}"
)
return _resolve_threads(imprint_threads, "imprint_threads")
def resolve_mesh_threads(mesh_threads: int | Literal["all"] = "all") -> None:
"""Validate the cad-to-dagmc-mesher meshing thread count.
It is passed to the mesher as given, so this only raises for a value that
is neither ``"all"`` nor a positive int, before any work is done.
"""
if mesh_threads != "all":
_resolve_threads(mesh_threads, "mesh_threads")
def _resolve_threads(value, name: str) -> int:
"""A thread count from ``"all"`` (every OCC logical processor, at least
one) or a positive int; *name* is the argument named in the error."""
if value == "all":
return max(1, OSD_Parallel.NbLogicalProcessors_s())
try:
# operator.index accepts int-likes such as numpy integers; bool is
# an int subclass but never a meaningful thread count.
threads = None if isinstance(value, bool) else operator.index(value)
except TypeError:
threads = None
if threads is None or threads < 1:
raise ValueError(f"{name} must be 'all' or a positive int, got {value!r}")
return threads
def _resolve_glue(glue: str | None) -> str | None:
"""Validate the glue mode, matching cad-to-dagmc-mesher."""
if glue not in ("partial", "full", None):
raise ValueError(
f"imprint must be 'partial', 'full' or None, got {glue!r}"
)
return glue
@contextmanager
def thread_limit(threads: int | None):
"""Limit the threads OpenCASCADE uses, restoring the limit afterwards.
cadquery's setThreads sets the size of the OpenCASCADE thread pool that
the boolean operations behind imprinting run on. The pool is process wide,
so the previous size is put back on the way out and the cadquery
operations that follow are left running on as many threads as before.
Args:
threads: the number of threads to allow, or None to leave the pool
alone.
"""
if threads is None:
yield
return
previous = OSD_ThreadPool.DefaultPool_s().NbThreads()
setThreads(threads)
try:
yield
finally:
setThreads(previous)
def imprint_assembly(assembly, threads: int | None = None, glue: str | None = "partial"):
"""Imprint a CadQuery assembly into a connected compound.
Uses cadquery's BOPAlgo_Builder based imprint (CadQuery/cadquery#2069),
passing glue through to it. CadQuery 2.8.0's legacy MakeConnected imprint
always uses partial glue (BOPAlgo_GlueShift), so "partial" works without
a warning. "full" and None raise RuntimeError with upgrade advice.
Glue modes are the same as the export methods' imprint configuration.
Args:
assembly: the cadquery assembly to imprint.
threads: the number of threads to imprint with. Fewer threads lowers
the peak RAM of the imprint at the cost of speed. Defaults to None
which leaves the thread count as it is.
glue: the glue option for the imprint, "partial", "full" or None.
Defaults to "partial". Explicit None disables glue, not imprinting;
it does not select the default.
Returns:
(imprinted_shape, imprinted_solids_with_original_ids)
"""
glue = _resolve_glue(glue)
# Imprinting needs at least two solids to do anything. Skipping it for a
# single solid is not just an optimization: the BOPAlgo_Builder based
# imprint returns a Null shape when given fewer than two arguments.
id_map = {}
for obj, name, loc, _ in assembly:
for solid in obj.moved(loc).Solids():
id_map[solid] = name
if len(id_map) < 2:
solids = list(id_map)
compound = cq.occ_impl.shapes.Compound.makeCompound(solids)
return compound, {s: (id_map[s],) for s in solids}
imprint = cq.occ_impl.assembly.imprint
if "glue" not in inspect.signature(imprint).parameters:
if glue != "partial":
raise RuntimeError(
f"CadQuery {cq.__version__} cannot honour imprint={glue!r}: "
"its legacy imprint always uses partial glue (BOPAlgo_GlueShift). "
"Use imprint='partial', "
"or upgrade CadQuery to a version "
"including CadQuery/cadquery#2069."
)
with thread_limit(threads):
return imprint(assembly)
with thread_limit(threads):
return imprint(assembly, glue=glue)
def define_moab_core_and_tags():
"""Creates a MOAB Core instance which can be built up by adding sets of
triangles to the instance
Returns:
(pymoab Core): A pymoab.core.Core() instance
(pymoab tag_handle): A pymoab.core.tag_get_handle() instance
"""
try:
from pymoab import core, types
except ImportError as e:
raise PyMoabNotFoundError() from e
# create pymoab instance
moab_core = core.Core()
tags = {}
sense_tag_name = "GEOM_SENSE_2"
sense_tag_size = 2
tags["surf_sense"] = moab_core.tag_get_handle(
sense_tag_name,
sense_tag_size,
types.MB_TYPE_HANDLE,
types.MB_TAG_SPARSE,
create_if_missing=True,
)
tags["category"] = moab_core.tag_get_handle(
types.CATEGORY_TAG_NAME,
types.CATEGORY_TAG_SIZE,
types.MB_TYPE_OPAQUE,
types.MB_TAG_SPARSE,
create_if_missing=True,
)
tags["name"] = moab_core.tag_get_handle(
types.NAME_TAG_NAME,
types.NAME_TAG_SIZE,
types.MB_TYPE_OPAQUE,
types.MB_TAG_SPARSE,
create_if_missing=True,
)
tags["geom_dimension"] = moab_core.tag_get_handle(
types.GEOM_DIMENSION_TAG_NAME,
1,
types.MB_TYPE_INTEGER,
types.MB_TAG_DENSE,
create_if_missing=True,
)
# Global ID is a default tag, just need the name to retrieve
tags["global_id"] = moab_core.tag_get_handle(types.GLOBAL_ID_TAG_NAME)
return moab_core, tags
[docs]
def vertices_to_h5m(
vertices: list[tuple[float, float, float]] | list[Vector],
triangles_by_solid_by_face: dict[int, dict[int, list[list[int]]]],
material_tags: list[str],
h5m_filename: str = "dagmc.h5m",
implicit_complement_material_tag: str | None = None,
method: str = "h5py",
):
"""Converts vertices and triangle sets into a tagged h5m file compatible
with DAGMC enabled neutronics simulations
Args:
vertices: List of vertex coordinates as (x, y, z) tuples or CadQuery vectors
triangles_by_solid_by_face: Dict mapping solid_id -> face_id -> list of triangles
material_tags: List of material tag names, one per solid
h5m_filename: Output filename for the h5m file
implicit_complement_material_tag: Optional material tag for implicit complement
method: Backend to use for writing h5m file ('pymoab' or 'h5py')
"""
if method == "pymoab":
return _vertices_to_h5m_pymoab(
vertices=vertices,
triangles_by_solid_by_face=triangles_by_solid_by_face,
material_tags=material_tags,
h5m_filename=h5m_filename,
implicit_complement_material_tag=implicit_complement_material_tag,
)
elif method == "h5py":
return _vertices_to_h5m_h5py(
vertices=vertices,
triangles_by_solid_by_face=triangles_by_solid_by_face,
material_tags=material_tags,
h5m_filename=h5m_filename,
implicit_complement_material_tag=implicit_complement_material_tag,
)
else:
raise ValueError(f"method must be 'pymoab' or 'h5py', not '{method}'")
def _vertices_to_h5m_pymoab(
vertices: list[tuple[float, float, float]] | list[Vector],
triangles_by_solid_by_face: dict[int, dict[int, list[list[int]]]],
material_tags: list[str],
h5m_filename: str = "dagmc.h5m",
implicit_complement_material_tag: str | None = None,
):
"""PyMOAB backend for vertices_to_h5m."""
try:
from pymoab import types
except ImportError as e:
raise PyMoabNotFoundError() from e
if len(material_tags) != len(triangles_by_solid_by_face):
msg = f"The number of material_tags provided is {len(material_tags)} and the number of sets of triangles is {len(triangles_by_solid_by_face)}. You must provide one material_tag for every triangle set"
raise ValueError(msg)
# limited attribute checking to see if user passed in a list of CadQuery vectors
if (
hasattr(vertices[0], "x")
and hasattr(vertices[0], "y")
and hasattr(vertices[0], "z")
):
vertices_floats = []
for vert in vertices:
vertices_floats.append((vert.x, vert.y, vert.z))
else:
vertices_floats = vertices
face_ids_with_solid_ids = {}
for solid_id, triangles_on_each_face in triangles_by_solid_by_face.items():
for face_id, triangles_on_face in triangles_on_each_face.items():
if face_id in face_ids_with_solid_ids:
face_ids_with_solid_ids[face_id].append(solid_id)
else:
face_ids_with_solid_ids[face_id] = [solid_id]
moab_core, tags = define_moab_core_and_tags()
# Add the vertices once at the start
all_moab_verts = moab_core.create_vertices(vertices)
volume_sets_by_solid_id = {}
for material_tag, (solid_id, triangles_on_each_face) in zip(
material_tags, triangles_by_solid_by_face.items()
):
volume_set = moab_core.create_meshset()
volume_sets_by_solid_id[solid_id] = volume_set
added_surfaces_ids = {}
for material_tag, (solid_id, triangles_on_each_face) in zip(
material_tags, triangles_by_solid_by_face.items()
):
volume_set = volume_sets_by_solid_id[solid_id]
moab_core.tag_set_data(tags["global_id"], volume_set, solid_id)
moab_core.tag_set_data(tags["geom_dimension"], volume_set, 3)
moab_core.tag_set_data(tags["category"], volume_set, "Volume")
group_set = moab_core.create_meshset()
moab_core.tag_set_data(tags["category"], group_set, "Group")
moab_core.tag_set_data(tags["name"], group_set, f"mat:{material_tag}")
moab_core.tag_set_data(tags["global_id"], group_set, solid_id)
# moab_core.tag_set_data(tags["geom_dimension"], group_set, 4)
for face_id, triangles_on_face in triangles_on_each_face.items():
if face_id not in added_surfaces_ids:
face_set = moab_core.create_meshset()
moab_core.tag_set_data(tags["global_id"], face_set, face_id)
moab_core.tag_set_data(tags["geom_dimension"], face_set, 2)
moab_core.tag_set_data(tags["category"], face_set, "Surface")
if len(face_ids_with_solid_ids[face_id]) == 2:
other_solid_id = face_ids_with_solid_ids[face_id][1]
other_volume_set = volume_sets_by_solid_id[other_solid_id]
sense_data = np.array(
[other_volume_set, volume_set], dtype="uint64"
)
else:
sense_data = np.array([volume_set, 0], dtype="uint64")
moab_core.tag_set_data(tags["surf_sense"], face_set, sense_data)
# Collect only the vertices that lie on triangles on this face
face_vertices_set = set()
for triangle in triangles_on_face:
face_vertices_set.update(triangle)
face_vertices_list = sorted(face_vertices_set)
# Only add these to the MOAB face
moab_verts = [all_moab_verts[ii] for ii in face_vertices_list]
moab_core.add_entity(face_set, moab_verts)
for triangle in triangles_on_face:
tri = (
all_moab_verts[int(triangle[0])],
all_moab_verts[int(triangle[1])],
all_moab_verts[int(triangle[2])],
)
moab_triangle = moab_core.create_element(types.MBTRI, tri)
moab_core.add_entity(face_set, moab_triangle)
added_surfaces_ids[face_id] = face_set
else:
face_set = added_surfaces_ids[face_id]
other_solid_id = face_ids_with_solid_ids[face_id][0]
other_volume_set = volume_sets_by_solid_id[other_solid_id]
sense_data = np.array([other_volume_set, volume_set], dtype="uint64")
moab_core.tag_set_data(tags["surf_sense"], face_set, sense_data)
moab_core.add_parent_child(volume_set, face_set)
moab_core.add_entity(group_set, volume_set)
if implicit_complement_material_tag:
group_set = moab_core.create_meshset()
moab_core.tag_set_data(tags["category"], group_set, "Group")
moab_core.tag_set_data(
tags["name"], group_set, f"mat:{implicit_complement_material_tag}_comp"
)
moab_core.tag_set_data(tags["geom_dimension"], group_set, 4)
moab_core.add_entity(
group_set, volume_set
) # volume is arbitrary but should exist in moab core
all_sets = moab_core.get_entities_by_handle(0)
file_set = moab_core.create_meshset()
moab_core.add_entities(file_set, all_sets)
# makes the folder if it does not exist
if Path(h5m_filename).parent:
Path(h5m_filename).parent.mkdir(parents=True, exist_ok=True)
# moab_core.write_file only accepts strings
if isinstance(h5m_filename, Path):
moab_core.write_file(str(h5m_filename))
else:
moab_core.write_file(h5m_filename)
print(f"written DAGMC file {h5m_filename}")
return h5m_filename
def _vertices_to_h5m_h5py(
vertices: list[tuple[float, float, float]] | list[Vector],
triangles_by_solid_by_face: dict[int, dict[int, list[list[int]]]],
material_tags: list[str],
h5m_filename: str = "dagmc.h5m",
implicit_complement_material_tag: str | None = None,
):
"""H5PY backend for vertices_to_h5m.
Creates an h5m file compatible with DAGMC using h5py directly,
without requiring pymoab.
"""
from datetime import datetime
import h5py
if len(material_tags) != len(triangles_by_solid_by_face):
msg = f"The number of material_tags provided is {len(material_tags)} and the number of sets of triangles is {len(triangles_by_solid_by_face)}. You must provide one material_tag for every triangle set"
raise ValueError(msg)
# Convert CadQuery vectors to floats if needed
if (
hasattr(vertices[0], "x")
and hasattr(vertices[0], "y")
and hasattr(vertices[0], "z")
):
vertices_floats = [(vert.x, vert.y, vert.z) for vert in vertices]
else:
vertices_floats = vertices
# Build face_ids_with_solid_ids to track shared faces
face_ids_with_solid_ids = {}
for solid_id, triangles_on_each_face in triangles_by_solid_by_face.items():
for face_id in triangles_on_each_face:
if face_id in face_ids_with_solid_ids:
face_ids_with_solid_ids[face_id].append(solid_id)
else:
face_ids_with_solid_ids[face_id] = [solid_id]
# Collect all unique faces and their triangles
all_faces = {} # face_id -> list of triangles
for solid_id, triangles_on_each_face in triangles_by_solid_by_face.items():
for face_id, triangles_on_face in triangles_on_each_face.items():
if face_id not in all_faces:
all_faces[face_id] = triangles_on_face
# Convert vertices to numpy array
vertices_arr = np.asarray(vertices_floats, dtype=np.float64)
num_vertices = len(vertices_arr)
# Collect all triangles
all_triangles = []
for face_id in sorted(all_faces.keys()):
all_triangles.extend(all_faces[face_id])
all_triangles = np.asarray(all_triangles, dtype=np.int64)
num_triangles = len(all_triangles)
# Create the h5m file
# makes the folder if it does not exist
if Path(h5m_filename).parent:
Path(h5m_filename).parent.mkdir(parents=True, exist_ok=True)
with h5py.File(h5m_filename, "w") as f:
tstt = f.create_group("tstt")
# Global ID counter - starts at 1
global_id = 1
# === NODES ===
nodes_group = tstt.create_group("nodes")
coords = nodes_group.create_dataset("coordinates", data=vertices_arr)
coords.attrs.create("start_id", global_id)
global_id += num_vertices
# Node tags
node_tags = nodes_group.create_group("tags")
node_tags.create_dataset("GLOBAL_ID", data=np.full(num_vertices, -1, dtype=np.int32))
# === ELEMENTS ===
elements = tstt.create_group("elements")
# Element type enum
elems = {
"Edge": 1, "Tri": 2, "Quad": 3, "Polygon": 4, "Tet": 5,
"Pyramid": 6, "Prism": 7, "Knife": 8, "Hex": 9, "Polyhedron": 10,
}
tstt["elemtypes"] = h5py.enum_dtype(elems)
# History
now = datetime.now().astimezone()
tstt.create_dataset(
"history",
data=[
"cad_to_dagmc".encode("ascii"),
__version__.encode("ascii"),
now.strftime("%m/%d/%y").encode("ascii"),
now.strftime("%H:%M:%S").encode("ascii"),
],
)
# Triangles
tri3_group = elements.create_group("Tri3")
tri3_group.attrs.create("element_type", elems["Tri"], dtype=tstt["elemtypes"])
# Node indices are 1-based in h5m
connectivity = tri3_group.create_dataset(
"connectivity",
data=all_triangles + 1,
dtype=np.uint64,
)
triangle_start_id = global_id
connectivity.attrs.create("start_id", triangle_start_id)
global_id += num_triangles
# Triangle tags
tags_tri3 = tri3_group.create_group("tags")
tags_tri3.create_dataset("GLOBAL_ID", data=np.full(num_triangles, -1, dtype=np.int32))
# === SETS ===
# Plan out the entity set structure:
# For each solid: 1 volume set, N surface sets (one per face), 1 group set (material)
# Plus: 1 file set at the end, optionally 1 implicit complement group
solid_ids = list(triangles_by_solid_by_face)
# Assign set IDs
sets_start_id = global_id
# Map solid_id -> volume_set_id
volume_set_ids = {}
# Map face_id -> surface_set_id
surface_set_ids = {}
# Map solid_id -> group_set_id
group_set_ids = {}
current_set_id = sets_start_id
# First, assign IDs to all surfaces (one per unique face)
for face_id in sorted(all_faces.keys()):
surface_set_ids[face_id] = current_set_id
current_set_id += 1
# Then assign IDs to volumes
for solid_id in solid_ids:
volume_set_ids[solid_id] = current_set_id
current_set_id += 1
# Then assign IDs to groups (materials)
for solid_id in solid_ids:
group_set_ids[solid_id] = current_set_id
current_set_id += 1
# Implicit complement group (if requested)
implicit_complement_set_id = None
if implicit_complement_material_tag:
implicit_complement_set_id = current_set_id
current_set_id += 1
# File set
file_set_id = current_set_id
current_set_id += 1
global_id = current_set_id
# === TAGS ===
tstt_tags = tstt.create_group("tags")
# Collect tagged set IDs for CATEGORY (all entities)
# and GEOM_DIMENSION (only surfaces and volumes - not groups, to match pymoab)
category_set_ids = []
categories = []
geom_dim_set_ids = []
geom_dimensions = []
# Volumes first (to match pymoab ordering)
for solid_id in solid_ids:
category_set_ids.append(volume_set_ids[solid_id])
categories.append("Volume")
geom_dim_set_ids.append(volume_set_ids[solid_id])
geom_dimensions.append(3)
# Groups (CATEGORY only - pymoab doesn't set geom_dimension on groups)
# Note: Groups COULD have geom_dimension=4 set, but pymoab doesn't do this
for solid_id in solid_ids:
category_set_ids.append(group_set_ids[solid_id])
categories.append("Group")
# Surfaces
for face_id in sorted(all_faces.keys()):
category_set_ids.append(surface_set_ids[face_id])
categories.append("Surface")
geom_dim_set_ids.append(surface_set_ids[face_id])
geom_dimensions.append(2)
# Implicit complement (CATEGORY only)
if implicit_complement_material_tag:
category_set_ids.append(implicit_complement_set_id)
categories.append("Group")
# CATEGORY tag
# Note: We use opaque dtype (|V32) to match pymoab output exactly.
# A string dtype (|S32) would also work and be more readable in h5dump,
# but we match pymoab for maximum compatibility.
cat_group = tstt_tags.create_group("CATEGORY")
cat_group.attrs.create("class", 1, dtype=np.int32)
cat_group.create_dataset("id_list", data=np.array(category_set_ids, dtype=np.uint64))
# Create opaque 32-byte type to match pymoab's H5T_OPAQUE
opaque_dt = h5py.opaque_dtype(np.dtype("V32"))
cat_group["type"] = opaque_dt
# Encode category strings as 32-byte null-padded values
cat_values = np.array([s.encode("ascii").ljust(32, b"\x00") for s in categories], dtype="V32")
cat_group.create_dataset("values", data=cat_values)
# GEOM_DIMENSION tag
# Note: We only tag surfaces (dim=2) and volumes (dim=3), not groups.
# Groups COULD have geom_dimension=4, but pymoab doesn't set this.
geom_group = tstt_tags.create_group("GEOM_DIMENSION")
geom_group["type"] = np.dtype("i4")
geom_group.attrs.create("class", 1, dtype=np.int32)
geom_group.attrs.create("default", -1, dtype=geom_group["type"])
geom_group.attrs.create("global", -1, dtype=geom_group["type"])
geom_group.create_dataset("id_list", data=np.array(geom_dim_set_ids, dtype=np.uint64))
geom_group.create_dataset("values", data=np.array(geom_dimensions, dtype=np.int32))
# GEOM_SENSE_2 tag (only for surfaces)
surface_ids_list = [surface_set_ids[fid] for fid in sorted(all_faces.keys())]
gs2_group = tstt_tags.create_group("GEOM_SENSE_2")
gs2_dtype = np.dtype("(2,)u8")
gs2_group["type"] = gs2_dtype
gs2_group.attrs.create("class", 1, dtype=np.int32)
gs2_group.attrs.create("is_handle", 1, dtype=np.int32)
gs2_group.create_dataset("id_list", data=np.array(surface_ids_list, dtype=np.uint64))
# Build sense data for each surface
sense_values = []
for face_id in sorted(all_faces.keys()):
solids_for_face = face_ids_with_solid_ids[face_id]
if len(solids_for_face) == 2:
# Shared face - both volumes
vol1 = volume_set_ids[solids_for_face[0]]
vol2 = volume_set_ids[solids_for_face[1]]
sense_values.append([vol1, vol2])
else:
# Single volume
vol = volume_set_ids[solids_for_face[0]]
sense_values.append([vol, 0])
if sense_values:
gs2_values = np.zeros((len(sense_values),), dtype=[("f0", "<u8", (2,))])
gs2_values["f0"] = np.array(sense_values, dtype=np.uint64)
gs2_space = h5py.h5s.create_simple((len(sense_values),))
gs2_arr_type = h5py.h5t.array_create(h5py.h5t.NATIVE_UINT64, (2,))
gs2_dset = h5py.h5d.create(gs2_group.id, b"values", gs2_arr_type, gs2_space)
gs2_dset.write(h5py.h5s.ALL, h5py.h5s.ALL, gs2_values, mtype=gs2_arr_type)
gs2_dset.close()
# GLOBAL_ID tag - store as sparse tag with id_list and values
# This stores the user-facing IDs for surfaces and volumes
gid_ids = []
gid_values = []
# Surfaces get their face_id as global_id
for face_id in sorted(all_faces.keys()):
gid_ids.append(surface_set_ids[face_id])
gid_values.append(face_id)
# Volumes get their solid_id as global_id
for solid_id in solid_ids:
gid_ids.append(volume_set_ids[solid_id])
gid_values.append(solid_id)
# Groups also get the solid_id
for solid_id in solid_ids:
gid_ids.append(group_set_ids[solid_id])
gid_values.append(solid_id)
gid_group = tstt_tags.create_group("GLOBAL_ID")
gid_group["type"] = np.dtype("i4")
gid_group.attrs.create("class", 2, dtype=np.int32)
gid_group.attrs.create("default", -1, dtype=gid_group["type"])
gid_group.attrs.create("global", -1, dtype=gid_group["type"])
gid_group.create_dataset("id_list", data=np.array(gid_ids, dtype=np.uint64))
gid_group.create_dataset("values", data=np.array(gid_values, dtype=np.int32))
# NAME tag (for groups - material names)
name_ids = []
name_values = []
for solid_id, mat_tag in zip(solid_ids, material_tags):
name_ids.append(group_set_ids[solid_id])
name_values.append(f"mat:{mat_tag}")
if implicit_complement_material_tag:
name_ids.append(implicit_complement_set_id)
name_values.append(f"mat:{implicit_complement_material_tag}_comp")
name_group = tstt_tags.create_group("NAME")
name_group.attrs.create("class", 1, dtype=np.int32)
name_group.create_dataset("id_list", data=np.array(name_ids, dtype=np.uint64))
name_group["type"] = h5py.opaque_dtype(np.dtype("S32"))
name_group.create_dataset("values", data=name_values, dtype=name_group["type"])
# Other standard tags (empty but needed)
for tag_name in ["DIRICHLET_SET", "MATERIAL_SET", "NEUMANN_SET"]:
tag_grp = tstt_tags.create_group(tag_name)
tag_grp["type"] = np.dtype("i4")
tag_grp.attrs.create("class", 1, dtype=np.int32)
tag_grp.attrs.create("default", -1, dtype=tag_grp["type"])
tag_grp.attrs.create("global", -1, dtype=tag_grp["type"])
# FACETING_TOLERANCE tag — stored on the root meshset via the
# "global" dataset so DAGMC's GeomQueryTool reads a valid value.
# Without this, DAGMC reads uninitialised memory and particle
# tracking fails with lost particles at curved surface boundaries.
ft_grp = tstt_tags.create_group("FACETING_TOLERANCE")
ft_type = np.dtype("f8")
ft_grp["type"] = ft_type
ft_grp.attrs.create("class", 2, dtype=np.int32)
# Compute a representative faceting tolerance from the mesh extent.
_diag = np.linalg.norm(vertices_arr.max(axis=0) - vertices_arr.min(axis=0))
_facet_tol = max(_diag * 1e-3, 1e-3)
# MOAB's mhdf reader expects "default" and "global" as HDF5
# datasets (not attributes). Store them both ways for compat.
ft_grp.create_dataset("default", data=np.array([_facet_tol], dtype=ft_type))
ft_grp.create_dataset("global", data=np.array([_facet_tol], dtype=ft_type))
# Also store as sparse tag data on root meshset (handle 0).
ft_grp.create_dataset("id_list", data=np.array([0], dtype=np.uint64))
ft_grp.create_dataset("values", data=np.array([_facet_tol], dtype=ft_type))
# === SETS structure ===
sets_group = tstt.create_group("sets")
# Build contents, parents, children, and list arrays
contents = []
list_rows = []
parents_list = []
children_list = []
# Track triangle ranges per face
tri_offset = 0
face_triangle_ranges = {}
for face_id in sorted(all_faces.keys()):
tris = all_faces[face_id]
face_triangle_ranges[face_id] = (tri_offset, len(tris))
tri_offset += len(tris)
# Track vertices per face
face_vertex_sets = {}
for face_id, tris in all_faces.items():
verts = set()
for tri in tris:
verts.update(tri)
face_vertex_sets[face_id] = sorted(verts)
contents_end = -1
children_end = -1
parents_end = -1
# Surface sets
for face_id in sorted(all_faces.keys()):
# Content: vertices + triangles for this face
verts = face_vertex_sets[face_id]
tri_start, tri_count = face_triangle_ranges[face_id]
# Add individual vertex handles (1-based IDs)
# Don't assume vertices are contiguous - store each one
for v in verts:
contents.append(v + 1) # 1-based vertex ID
# Add individual triangle handles
for i in range(tri_count):
contents.append(triangle_start_id + tri_start + i)
contents_end = len(contents) - 1
# Parent-child: surface is child of volume(s)
solids_for_face = face_ids_with_solid_ids[face_id]
for solid_id in solids_for_face:
parents_list.append(volume_set_ids[solid_id])
parents_end = len(parents_list) - 1
# flags: 2 = MESHSET_SET (handles, not ranges)
list_rows.append([contents_end, children_end, parents_end, 2])
# Volume sets (empty contents, but have surface children)
for solid_id in solid_ids:
# Volumes have no direct content
# Children are the surfaces
faces_in_solid = list(triangles_by_solid_by_face[solid_id].keys())
for face_id in faces_in_solid:
children_list.append(surface_set_ids[face_id])
children_end = len(children_list) - 1
# flags: 2 = handle-based (0b0010)
list_rows.append([contents_end, children_end, parents_end, 2])
# Group sets (contain volume handles)
for solid_id in solid_ids:
contents.append(volume_set_ids[solid_id])
contents_end = len(contents) - 1
list_rows.append([contents_end, children_end, parents_end, 2])
# Implicit complement group
if implicit_complement_material_tag:
# Add the last volume to the implicit complement group
contents.append(volume_set_ids[solid_ids[-1]])
contents_end = len(contents) - 1
list_rows.append([contents_end, children_end, parents_end, 2])
# File set (contains everything)
contents.extend([1, file_set_id - 1]) # range of all entities
contents_end = len(contents) - 1
list_rows.append([contents_end, children_end, parents_end, 10])
# Write sets datasets
sets_group.create_dataset("contents", data=np.array(contents, dtype=np.uint64))
if children_list:
sets_group.create_dataset("children", data=np.array(children_list, dtype=np.uint64))
else:
sets_group.create_dataset("children", data=np.array([], dtype=np.uint64))
if parents_list:
sets_group.create_dataset("parents", data=np.array(parents_list, dtype=np.uint64))
else:
sets_group.create_dataset("parents", data=np.array([], dtype=np.uint64))
lst = sets_group.create_dataset("list", data=np.array(list_rows, dtype=np.int64))
lst.attrs.create("start_id", sets_start_id)
# Set tags (GLOBAL_ID for each set)
sets_tags = sets_group.create_group("tags")
set_global_ids = []
# Surface global IDs
set_global_ids.extend(sorted(all_faces))
# Volume global IDs
set_global_ids.extend(solid_ids)
# Group global IDs
set_global_ids.extend(solid_ids)
# Implicit complement
if implicit_complement_material_tag:
set_global_ids.append(-1)
# File set
set_global_ids.append(-1)
sets_tags.create_dataset("GLOBAL_ID", data=np.array(set_global_ids, dtype=np.int32))
# Max ID attribute
tstt.attrs.create("max_id", np.uint64(global_id - 1))
print(f"written DAGMC file {h5m_filename}")
return h5m_filename
[docs]
def get_volumes(gmsh, assembly, method="file", scale_factor=1.0):
if method == "in memory":
volumes = gmsh.model.occ.importShapesNativePointer(assembly.wrapped._address())
elif method == "file":
with tempfile.NamedTemporaryFile(suffix=".brep") as temp_file:
if isinstance(assembly, cq.Assembly):
assembly.toCompound().exportBrep(temp_file.name)
else:
assembly.exportBrep(temp_file.name)
volumes = gmsh.model.occ.importShapes(temp_file.name)
# updating the model to ensure the entities in the geometry are found
gmsh.model.occ.synchronize()
if scale_factor != 1.0:
dim_tags = gmsh.model.getEntities(3)
gmsh.model.occ.dilate(
dim_tags, 0.0, 0.0, 0.0, scale_factor, scale_factor, scale_factor
)
# update the model to ensure the scaling factor has been applied
gmsh.model.occ.synchronize()
return gmsh, volumes
[docs]
def init_gmsh():
# gmsh is a global singleton. If a previous session was left initialized
# (for example by an export that errored part way through, or an earlier
# call that did not finalize) then adding a new model here would leave the
# stale models from that session alive, leaking memory and growing the
# session on every call (see issue #187). Finalize any pre-existing
# session first so we always start from a clean, single-model state.
if gmsh.isInitialized():
gmsh.finalize()
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 1)
gmsh.model.add(f"made_with_cad_to_dagmc_package_{__version__}")
return gmsh
[docs]
def set_sizes_for_mesh(
gmsh,
min_mesh_size: float | None = None,
max_mesh_size: float | None = None,
mesh_algorithm: int = 1,
set_size: dict[int | str, float] | None = None,
original_set_size: dict[int | str, float] | None = None,
threads: int = 0,
*,
mesh_algorithm_3d: int = 1,
):
"""Sets up the mesh sizes for each volume in the mesh.
Args:
occ_shape: the occ_shape of the Brep file to convert
min_mesh_size: the minimum mesh element size to use in Gmsh. Passed
into gmsh.option.setNumber("Mesh.MeshSizeMin", min_mesh_size)
max_mesh_size: the maximum mesh element size to use in Gmsh. Passed
into gmsh.option.setNumber("Mesh.MeshSizeMax", max_mesh_size)
mesh_algorithm: The Gmsh 2D surface mesh algorithm number. Passed into
gmsh.option.setNumber("Mesh.Algorithm", mesh_algorithm)
mesh_algorithm_3d: The Gmsh 3D volume mesh algorithm number. Passed into
gmsh.option.setNumber("Mesh.Algorithm3D", mesh_algorithm_3d).
Defaults to 1 (Delaunay); 10 selects HXT.
set_size: a dictionary of volume ids (int) and target mesh sizes
(floats) to set for each volume, passed to gmsh.model.mesh.setSize.
threads: the number of threads for Gmsh to use. Passed into
gmsh.option.setNumber("General.NumThreads", threads). 0 uses
all available cores (default), 1 uses a single thread.
Returns:
The resulting gmsh object and volumes
"""
if min_mesh_size and max_mesh_size and min_mesh_size > max_mesh_size:
raise ValueError(
f"min_mesh_size must be less than or equal to max_mesh_size. Currently min_mesh_size is set to {min_mesh_size} and max_mesh_size is set to {max_mesh_size}"
)
if min_mesh_size:
gmsh.option.setNumber("Mesh.MeshSizeMin", min_mesh_size)
if max_mesh_size:
gmsh.option.setNumber("Mesh.MeshSizeMax", max_mesh_size)
gmsh.option.setNumber("Mesh.Algorithm", mesh_algorithm)
gmsh.option.setNumber("Mesh.Algorithm3D", mesh_algorithm_3d)
gmsh.option.setNumber("General.NumThreads", threads)
if set_size:
volumes = gmsh.model.getEntities(3)
available_volumes = [volume[1] for volume in volumes]
print("volumes", volumes)
# Ensure all volume IDs in set_size exist in the available volumes
for volume_id in set_size:
if volume_id not in available_volumes:
raise ValueError(
f"volume ID of {volume_id} set in set_sizes but not found in available volumes {volumes}"
)
# Warn if any set_size values fall outside the global min/max range
# Use original_set_size keys (which may be material tag strings) for
# user-friendly warnings, falling back to resolved volume IDs
warn_items = original_set_size.items() if original_set_size else set_size.items()
for key, size in warn_items:
if min_mesh_size is not None and size < min_mesh_size:
warnings.warn(
f"set_size for {key} is {size} which is below "
f"min_mesh_size of {min_mesh_size}. The mesh size will be "
f"clamped to {min_mesh_size}. Try reducing min_mesh_size to "
f"encompass the set_size value."
)
if max_mesh_size is not None and size > max_mesh_size:
warnings.warn(
f"set_size for {key} is {size} which is above "
f"max_mesh_size of {max_mesh_size}. The mesh size will be "
f"clamped to {max_mesh_size}. Try enlarging max_mesh_size to "
f"encompass the set_size value."
)
# Step 1: Preprocess boundaries to find the smallest size for shared surfaces
boundary_sizes = {} # Dictionary to store the minimum mesh size for each boundary
for volume_id, size in set_size.items():
boundaries = gmsh.model.getBoundary(
[(3, volume_id)], recursive=True
) # dim must be set to 3
print(f"Boundaries for volume {volume_id}: {boundaries}")
for boundary in boundaries:
boundary_key = (boundary[0], boundary[1]) # (dimension, tag)
if boundary_key in boundary_sizes:
# If the boundary is already processed, keep the smaller size
boundary_sizes[boundary_key] = min(boundary_sizes[boundary_key], size)
else:
boundary_sizes[boundary_key] = size
# Step 2: Apply mesh sizes to all boundaries
for boundary, size in boundary_sizes.items():
gmsh.model.mesh.setSize([boundary], size)
print(f"Set mesh size {size} for boundary {boundary}")
return gmsh
def _generate_volume_mesh(
gmsh,
volume_mesh_options: dict[str, float | str] | None = None,
*,
surfaces_meshed: bool = False,
):
"""Apply volume options only after the surface mesh has been generated.
Keep the single generate(3) call when no overrides are provided. The DAGMC
export already meshes surfaces for the h5m file, so it can reuse that mesh.
Gmsh validates option names and values; errors propagate to the exporter,
which finalizes the session on failure as well as on success.
"""
if volume_mesh_options:
if not surfaces_meshed:
gmsh.model.mesh.generate(2)
for name, value in volume_mesh_options.items():
if isinstance(value, str):
gmsh.option.setString(name, value)
else:
gmsh.option.setNumber(name, value)
gmsh.model.mesh.generate(3)
[docs]
def mesh_to_vertices_and_triangles(
dims_and_vol_ids,
):
"""Converts gmsh volumes into vertices and triangles for each face.
Args:
volumes: the volumes in the gmsh file, found with gmsh.model.occ.importShapes
Returns:
vertices and triangles (grouped by solid then by face)
"""
n = 3 # number of verts in a triangles
triangles_by_solid_by_face = {}
for dim_and_vol in dims_and_vol_ids:
# removes all groups so that the following getEntitiesForPhysicalGroup
# command only finds surfaces for the volume
face_groups = gmsh.model.getPhysicalGroups(2)
if face_groups: # Only remove if 2D groups exist
gmsh.model.removePhysicalGroups(face_groups)
vol_id = dim_and_vol[1]
entities_in_volume = gmsh.model.getAdjacencies(3, vol_id)
surfaces_in_volume = entities_in_volume[1]
ps = gmsh.model.addPhysicalGroup(2, surfaces_in_volume)
gmsh.model.setPhysicalName(2, ps, f"surfaces_on_volume_{vol_id}")
groups = gmsh.model.getPhysicalGroups()
group = groups[0]
# for group in groups:
dim = group[0]
tag = group[1]
surfaces = gmsh.model.getEntitiesForPhysicalGroup(dim, tag)
# nodes_in_all_surfaces = []
nodes_in_each_surface = {}
for surface in surfaces:
_, _, nodeTags = gmsh.model.mesh.getElements(2, surface)
nodeTags = nodeTags[0].tolist()
shifted_node_tags = []
for nodeTag in nodeTags:
shifted_node_tags.append(nodeTag - 1)
grouped_node_tags = [
shifted_node_tags[i : i + n]
for i in range(0, len(shifted_node_tags), n)
]
nodes_in_each_surface[surface] = grouped_node_tags
triangles_by_solid_by_face[vol_id] = nodes_in_each_surface
_, all_coords, _ = gmsh.model.mesh.getNodes()
vertices = [all_coords[i : i + n].tolist() for i in range(0, len(all_coords), n)]
return vertices, triangles_by_solid_by_face
def get_ids_from_assembly(assembly: cq.assembly.Assembly):
ids = []
for obj, name, loc, _ in assembly:
ids.append(name)
return ids
def get_ids_from_imprinted_assembly(solid_id_dict):
ids = []
for id in list(solid_id_dict.values()):
ids.append(id[0])
return ids
def check_material_tags(material_tags, iterable_solids):
if material_tags:
if len(material_tags) != len(iterable_solids):
msg = (
"When setting material_tags the number of material_tags \n"
"should be equal to the number of volumes in the CAD \n"
f"geometry {len(iterable_solids)} volumes found in model \n"
f"and {len(material_tags)} material_tags found"
)
raise ValueError(msg)
for material_tag in material_tags:
if not isinstance(material_tag, str):
msg = "material_tags should be an iterable of strings."
raise TypeError(msg)
if len(material_tag) > 28:
msg = (
f"Material tag {material_tag} is too long. DAGMC will truncate this material tag "
f"to 28 characters. The resulting tag in the h5m file will be {material_tag[:28]}"
)
warnings.warn(msg)
def order_material_ids_by_brep_order(original_ids, scrambled_id, material_tags):
material_tags_in_brep_order = []
for brep_id in scrambled_id:
id_of_solid_in_org = original_ids.index(brep_id)
material_tags_in_brep_order.append(material_tags[id_of_solid_in_org])
return material_tags_in_brep_order
def resolve_unstructured_volumes(
unstructured_volumes: Iterable[int | str],
volumes: list[tuple[int, int]],
material_tags: list[str],
) -> list[int]:
"""Resolve a mixed list of volume IDs (int) and material tags (str) to volume IDs.
Args:
unstructured_volumes: An iterable containing volume IDs (int) or material tag
names (str). Material tags are resolved to all volume IDs that have that tag.
volumes: List of (dim, volume_id) tuples from GMSH, where the order corresponds
to the order of material_tags.
material_tags: List of material tags in the same order as volumes.
Returns:
A list of unique volume IDs (int) corresponding to the input.
Raises:
ValueError: If a material tag string is not found in material_tags.
TypeError: If an element is neither int nor str.
"""
resolved_ids = []
# Build a mapping from material tag to volume IDs
# volumes is a list of (dim, volume_id), and material_tags has the same order
material_to_volume_ids: dict[str, list[int]] = {}
for (_, volume_id), material_tag in zip(volumes, material_tags):
if material_tag not in material_to_volume_ids:
material_to_volume_ids[material_tag] = []
material_to_volume_ids[material_tag].append(volume_id)
for item in unstructured_volumes:
if isinstance(item, int):
resolved_ids.append(item)
elif isinstance(item, str):
if item not in material_to_volume_ids:
available_tags = sorted(set(material_tags))
raise ValueError(
f"Material tag '{item}' not found. "
f"Available material tags are: {available_tags}"
)
resolved_ids.extend(material_to_volume_ids[item])
else:
raise TypeError(
f"unstructured_volumes must contain int (volume ID) or str (material tag), "
f"got {type(item).__name__}"
)
# Remove duplicates while preserving order
seen = set()
unique_ids = []
for vol_id in resolved_ids:
if vol_id not in seen:
seen.add(vol_id)
unique_ids.append(vol_id)
return unique_ids
def resolve_set_size(
set_size: dict[int | str, float],
volumes: list[tuple[int, int]],
material_tags: list[str],
) -> dict[int, float]:
"""Resolve a set_size dict with int or str keys to int keys only.
Args:
set_size: A dictionary mapping volume IDs (int) or material tag names (str)
to mesh sizes (float). Material tags are resolved to all volume IDs
that have that tag.
volumes: List of (dim, volume_id) tuples from GMSH, where the order corresponds
to the order of material_tags.
material_tags: List of material tags in the same order as volumes.
Returns:
A dictionary mapping volume IDs (int) to mesh sizes (float).
Raises:
ValueError: If a material tag string is not found in material_tags,
or if a volume ID is specified multiple times with different sizes.
TypeError: If a key is neither int nor str.
"""
resolved: dict[int, float] = {}
# Build a mapping from material tag to volume IDs
material_to_volume_ids: dict[str, list[int]] = {}
for (_, volume_id), material_tag in zip(volumes, material_tags):
if material_tag not in material_to_volume_ids:
material_to_volume_ids[material_tag] = []
material_to_volume_ids[material_tag].append(volume_id)
for key, size in set_size.items():
if isinstance(key, int):
volume_ids = [key]
elif isinstance(key, str):
if key not in material_to_volume_ids:
available_tags = sorted(set(material_tags))
raise ValueError(
f"Material tag '{key}' not found in set_size. "
f"Available material tags are: {available_tags}"
)
volume_ids = material_to_volume_ids[key]
else:
raise TypeError(
f"set_size keys must be int (volume ID) or str (material tag), "
f"got {type(key).__name__}"
)
for vol_id in volume_ids:
if vol_id in resolved:
if resolved[vol_id] != size:
raise ValueError(
f"Volume ID {vol_id} specified multiple times with different sizes: "
f"{resolved[vol_id]} and {size}. "
f"Each volume can only have one mesh size."
)
else:
resolved[vol_id] = size
return resolved
[docs]
def export_gmsh_object_to_dagmc_h5m_file(
material_tags: list[str] | None = None,
implicit_complement_material_tag: str | None = None,
filename: str = "dagmc.h5m",
h5m_backend: str = "h5py",
) -> str:
"""
Exports a GMSH object to a DAGMC-compatible h5m file. Note gmsh should
be initialized by the user prior and the gmsh model should be meshed before
calling this. Also users should ensure that the gmsh model is finalized.
Args:
material_tags: A list of material tags corresponding to the volumes in the GMSH object.
implicit_complement_material_tag: The material tag for the implicit complement (void space).
filename: The name of the output h5m file. Defaults to "dagmc.h5m".
h5m_backend: Backend for writing h5m file, 'pymoab' or 'h5py'. Defaults to 'h5py'.
Returns:
str: The filename of the generated DAGMC h5m file.
Raises:
ValueError: If the number of material tags does not match the number of volumes in the GMSH object.
"""
if material_tags is None:
material_tags = _get_material_tags_from_gmsh()
dims_and_vol_ids = gmsh.model.getEntities(3)
if len(dims_and_vol_ids) != len(material_tags):
msg = f"Number of volumes {len(dims_and_vol_ids)} is not equal to number of material tags {len(material_tags)}"
raise ValueError(msg)
vertices, triangles_by_solid_by_face = mesh_to_vertices_and_triangles(
dims_and_vol_ids=dims_and_vol_ids
)
h5m_filename = vertices_to_h5m(
vertices=vertices,
triangles_by_solid_by_face=triangles_by_solid_by_face,
material_tags=material_tags,
h5m_filename=filename,
implicit_complement_material_tag=implicit_complement_material_tag,
method=h5m_backend,
)
return h5m_filename
def _get_material_tags_from_gmsh() -> list[str]:
"""Gets the Physical groups of 3D groups from the GMSH object and returns
their names."""
# Get all 3D physical groups (volumes)
volume_groups = gmsh.model.getPhysicalGroups(3)
material_tags = []
# Get the name for each physical group
for dim, tag in volume_groups:
name = gmsh.model.getPhysicalName(dim, tag)
material_tags.append(name)
print(f"Material tag: {name}")
print(f"Material tags: {material_tags}")
return material_tags
[docs]
def export_gmsh_file_to_dagmc_h5m_file(
gmsh_filename: str,
material_tags: list[str] | None = None,
implicit_complement_material_tag: str | None = None,
dagmc_filename: str = "dagmc.h5m",
h5m_backend: str = "h5py",
) -> str:
"""Saves a DAGMC h5m file of the geometry GMsh file. This function
initializes and finalizes Gmsh.
Args:
gmsh_filename (str): the filename of the GMSH mesh file.
material_tags (list[str]): the names of the DAGMC
material tags to assign. These will need to be in the same
order as the volumes in the GMESH mesh and match the
material tags used in the neutronics code (e.g. OpenMC).
implicit_complement_material_tag (str | None, optional):
the name of the material tag to use for the implicit
complement (void space). Defaults to None which is a vacuum.
dagmc_filename (str, optional): Output filename. Defaults to "dagmc.h5m".
h5m_backend (str, optional): Backend for writing h5m file, 'pymoab' or 'h5py'.
Defaults to 'h5py'.
Returns:
str: The filename of the generated DAGMC h5m file.
Raises:
ValueError: If the number of material tags does not match the number of volumes in the GMSH object.
"""
gmsh.initialize()
gmsh.open(gmsh_filename)
if material_tags is None:
material_tags = _get_material_tags_from_gmsh()
dims_and_vol_ids = gmsh.model.getEntities(3)
if len(dims_and_vol_ids) != len(material_tags):
msg = f"Number of volumes {len(dims_and_vol_ids)} is not equal to number of material tags {len(material_tags)}"
raise ValueError(msg)
vertices, triangles_by_solid_by_face = mesh_to_vertices_and_triangles(
dims_and_vol_ids=dims_and_vol_ids
)
gmsh.finalize()
h5m_filename = vertices_to_h5m(
vertices=vertices,
triangles_by_solid_by_face=triangles_by_solid_by_face,
material_tags=material_tags,
h5m_filename=dagmc_filename,
implicit_complement_material_tag=implicit_complement_material_tag,
method=h5m_backend,
)
return h5m_filename
[docs]
class CadToDagmc:
"""Converts Step files and CadQuery parts to a DAGMC h5m file"""
def __init__(self):
self.parts = []
self.material_tags = []
# The cad-to-dagmc-mesher refinement report from the last export, or
# None when no refinement ran.
self.refinement = None
[docs]
def add_stp_file(
self,
filename: str,
scale_factor: float = 1.0,
material_tags: list[str] | str | None = None,
) -> int:
"""Loads the parts from stp file into the model.
Args:
filename: the filename used to save the html graph.
material_tags: the names of the DAGMC material tags to assign.
Can be a list of strings (one per volume), or one of the
special strings "assembly_names" or "assembly_materials" to
automatically extract tags from the STEP file's assembly
structure (if the STEP file contains named parts or materials).
When using a list, tags must be in the same order as the
volumes in the geometry.
scale_factor: a scaling factor to apply to the geometry that can be
used to increase the size or decrease the size of the geometry.
Useful when converting the geometry to cm for use in neutronics
simulations.
Returns:
int: number of volumes in the stp file.
"""
# If using assembly_names or assembly_materials, try to load as assembly
if material_tags in ("assembly_names", "assembly_materials"):
assembly = cq.Assembly()
importStepAssembly(assembly, str(filename))
# add_cadquery_object scales the solids after placing them at their
# assembly locations, so nested parts and part locations scale too
return self.add_cadquery_object(
cadquery_object=assembly,
material_tags=material_tags,
scale_factor=scale_factor,
)
# Default behavior: load as compound/solid
part = importers.importStep(str(filename)).val()
if scale_factor == 1.0:
scaled_part = part
else:
scaled_part = part.scale(scale_factor)
return self.add_cadquery_object(
cadquery_object=scaled_part, material_tags=material_tags
)
[docs]
def add_cadquery_object(
self,
cadquery_object: (
cq.assembly.Assembly
| cq.occ_impl.shapes.Compound
| cq.occ_impl.shapes.Solid
),
material_tags: list[str] | str,
scale_factor: float = 1.0,
) -> int:
"""Loads the parts from CadQuery object into the model.
Args:
cadquery_object: the cadquery object to convert, can be a CadQuery assembly
cadquery workplane or a cadquery solid
material_tags (Optional list[str]): the names of the
DAGMC material tags to assign. These will need to be in the
same order as the volumes in the geometry added (STP file and
CadQuery objects) and match the material tags used in the
neutronics code (e.g. OpenMC).
scale_factor: a scaling factor to apply to the geometry that can be
used to increase the size or decrease the size of the geometry.
Useful when converting the geometry to cm for use in neutronics
simulations.
Returns:
int: number of volumes in the stp file.
"""
if isinstance(material_tags, str) and material_tags not in [
"assembly_materials",
"assembly_names",
]:
raise ValueError(
f"If material_tags is a string it must be 'assembly_materials' or 'assembly_names' but got {material_tags}"
)
if isinstance(cadquery_object, cq.assembly.Assembly):
# look for materials in each part of the assembly
if material_tags == "assembly_materials":
material_tags = []
for child in _get_all_leaf_children(cadquery_object):
if child.material is not None and child.material.name is not None:
# count solids in this child to repeat the tag appropriately
child_shape = child.toCompound() if hasattr(child, 'toCompound') else child.obj
if child_shape is not None:
child_solids = child_shape.Solids() if hasattr(child_shape, 'Solids') else []
else:
child_solids = []
for _ in child_solids:
material_tags.append(str(child.material.name))
else:
raise ValueError(
f"Not all parts in the assembly have materials assigned.\n"
f"When adding to an assembly include material=cadquery.Material('material_name')\n"
f"Missing material tag for child: {child}.\n"
"Please assign material tags to all parts or provide material_tags argument when adding the assembly.\n"
)
print("material_tags found from assembly materials:", material_tags)
elif material_tags == "assembly_names":
material_tags = []
for child in _get_all_leaf_children(cadquery_object):
# count solids in this child to repeat the tag appropriately
child_shape = child.toCompound() if hasattr(child, 'toCompound') else child.obj
if child_shape is not None:
child_solids = child_shape.Solids() if hasattr(child_shape, 'Solids') else []
else:
child_solids = []
# parts always have a name as cq will auto assign one
for _ in child_solids:
material_tags.append(child.name)
print("material_tags found from assembly names:", material_tags)
cadquery_compound = cadquery_object.toCompound()
else:
cadquery_compound = cadquery_object
if isinstance(
cadquery_compound, (cq.occ_impl.shapes.Compound, cq.occ_impl.shapes.Solid)
):
iterable_solids = cadquery_compound.Solids()
else:
iterable_solids = cadquery_compound.val().Solids()
scaled_iterable_solids = _scale_together(iterable_solids, scale_factor)
check_material_tags(material_tags, scaled_iterable_solids)
if material_tags:
self.material_tags = self.material_tags + material_tags
self.parts = self.parts + scaled_iterable_solids
return len(scaled_iterable_solids)
[docs]
def export_unstructured_mesh_file(
self,
filename: str = "umesh.vtk",
min_mesh_size: float = 1,
max_mesh_size: float = 5,
mesh_algorithm: int = 1,
method: str = "file",
scale_factor: float = 1.0,
imprint: Literal["partial", "full", False] | None = "partial",
set_size: dict[int | str, float] | None = None,
volumes: Iterable[int] | None = None,
threads: int = 0,
meshing_backend: str | None = None,
target_edge_length: float | None = None,
tet_volumes: Iterable[str] | None = None,
tolerance: float = DEFAULT_TOLERANCE,
angular_tolerance: float = DEFAULT_ANGULAR_TOLERANCE,
*,
imprint_threads: int | Literal["all"] = "all",
mesh_threads: int | Literal["all"] = "all",
mesh_algorithm_3d: int = 1,
volume_mesh_options: dict[str, float | str] | None = None,
):
"""
Exports an unstructured mesh file in VTK format for use with
openmc.UnstructuredMesh. Compatible with the MOAB unstructured mesh
library. Example useage openmc.UnstructuredMesh(filename="umesh.vtk",
library="moab").
The mesh can be produced either with gmsh or with the
cad-to-dagmc-mesher backend. The gmsh backend uses the min/max mesh
size and set_size arguments, while the cad-to-dagmc-mesher backend
uses target_edge_length (and optionally tet_volumes) to control the
tetrahedra. gmsh is used unless meshing_backend or one of the
cad-to-dagmc-mesher specific arguments is provided.
Parameters:
-----------
filename : str, optional
The name of the output file. Default is "umesh.vtk".
min_mesh_size: the minimum mesh element size to use in Gmsh. Passed
into gmsh.option.setNumber("Mesh.MeshSizeMin", min_mesh_size)
max_mesh_size: the maximum mesh element size to use in Gmsh. Passed
into gmsh.option.setNumber("Mesh.MeshSizeMax", max_mesh_size)
mesh_algorithm: The Gmsh 2D surface mesh algorithm number. Passed into
gmsh.option.setNumber("Mesh.Algorithm", mesh_algorithm)
mesh_algorithm_3d: Gmsh 3D volume mesh algorithm, passed to
Mesh.Algorithm3D. Defaults to 1 (Delaunay); 10 selects HXT.
Only used by the gmsh backend.
volume_mesh_options: Gmsh option names mapped to numeric or string
values, applied after surface meshing and before volume meshing.
Overrides the corresponding sizing/algorithm options for the
volume step only. None or an empty dict keeps the single
generate(3) call. Only used by the gmsh backend.
method: the method to use to import the geometry into gmsh. Options
are 'file' or 'in memory'. 'file' is the default and will write
the geometry to a temporary file before importing it into gmsh.
'in memory' will import the geometry directly into gmsh but
requires the version of OpenCASCADE used to build gmsh to be
the same as the version used by CadQuery. This is possible to
ensure when installing the package with Conda but harder when
installing from PyPI.
scale_factor: a scaling factor to apply to the geometry that can be
used to enlarge or shrink the geometry. Useful when converting
the geometry to cm for use in neutronics.
imprint: "partial" (default) for non-overlapping contact, "full" for
exactly coincident whole faces, None to imprint without glue, or
False to skip imprinting entirely. None does not mean the default.
CadQuery 2.8.0 supports only partial glue; other modes raise
RuntimeError before imprinting multiple solids. A single solid
skips imprinting, but arguments are always validated. True,
integers, and dictionaries are not accepted.
imprint_threads: keyword-only thread limit for imprinting: "all"
(default, all OCC logical processors) or a positive int.
Fewer threads reduce peak RAM at the cost of speed. The previous
OCC pool is restored before meshing, including after errors.
Validated even when imprint=False; independent of Gmsh threads.
mesh_threads: keyword-only thread count for the cad-to-dagmc-mesher
meshing that follows the imprint: "all" (default) or a positive
int. Fewer threads reduce peak RAM at the cost of speed, and
give the same mesh. Validated with any backend; the gmsh
backend warns that it is unused (see threads).
set_size: a dictionary mapping volume IDs (int) or material tag names
(str) to target mesh sizes (floats). Material tags are resolved to
all volume IDs that have that tag. Only used by the gmsh backend.
volumes: a list of volume ids (int) to include in the mesh. If left
as default (None) then all volumes will be included. Only used by
the gmsh backend.
threads: the number of threads for Gmsh to use. 0 uses all
available cores (default), 1 uses a single thread.
meshing_backend: the backend used to generate the tetrahedra, either
"gmsh" or "cad-to-dagmc-mesher". If not set, the backend is
auto-selected: "cad-to-dagmc-mesher" when target_edge_length or
tet_volumes is provided, otherwise "gmsh".
target_edge_length: the target tetrahedron edge length used by the
cad-to-dagmc-mesher backend. Required when meshing_backend is
"cad-to-dagmc-mesher".
tet_volumes: an iterable of material tag names identifying which
volumes to fill with tetrahedra when using the
cad-to-dagmc-mesher backend. Defaults to all volumes.
tolerance: linear deflection tolerance for the surface mesh, used by
the cad-to-dagmc-mesher backend. This is in the units of the
SCALED geometry, since scale_factor is applied before meshing,
so scale it alongside scale_factor. With scale_factor=100 the
0.01 default is a 0.1 mm deflection, which on a large model can
produce a very fine mesh and exhaust memory. The same applies to
min_mesh_size/max_mesh_size/set_size for the gmsh backend.
angular_tolerance: angular deflection tolerance for the surface mesh,
used by the cad-to-dagmc-mesher backend. An angle, so unaffected
by scale_factor.
Returns:
--------
filename : str
The filename of the written unstructured mesh file.
"""
# gmesh writes out a vtk file that is accepted by openmc.UnstructuredMesh
# The library argument must be set to "moab"
if Path(filename).suffix != ".vtk":
raise ValueError("Unstructured mesh filename must have a .vtk extension")
resolved_imprint_threads = resolve_imprint(imprint, imprint_threads)
resolve_mesh_threads(mesh_threads)
if meshing_backend is None:
# Auto-select the backend: the tet arguments are specific to
# cad-to-dagmc-mesher, everything else defaults to gmsh.
if target_edge_length is not None or tet_volumes is not None:
meshing_backend = "cad-to-dagmc-mesher"
else:
meshing_backend = "gmsh"
print(f"Using meshing backend: {meshing_backend}")
if meshing_backend not in ("gmsh", "cad-to-dagmc-mesher"):
raise ValueError(
f'meshing_backend "{meshing_backend}" not supported. '
'Available options are "gmsh" or "cad-to-dagmc-mesher"'
)
if meshing_backend == "cad-to-dagmc-mesher":
unused_params = []
if mesh_algorithm_3d != 1:
unused_params.append("mesh_algorithm_3d")
if volume_mesh_options is not None:
unused_params.append("volume_mesh_options")
if unused_params:
warnings.warn(
"The following parameters are ignored when using "
f"cad-to-dagmc-mesher backend: {', '.join(unused_params)}"
)
return self._export_unstructured_mesh_file_with_mesher(
filename=filename,
target_edge_length=target_edge_length,
tet_volumes=tet_volumes,
tolerance=tolerance,
angular_tolerance=angular_tolerance,
imprint=imprint,
imprint_threads=imprint_threads,
mesh_threads=mesh_threads,
scale_factor=scale_factor,
)
if mesh_threads != "all":
warnings.warn(
"The following parameters are ignored when using GMSH backend: "
"mesh_threads (use threads)"
)
assembly = cq.Assembly()
for part in self.parts:
assembly.add(part)
if imprint is not False:
print("Imprinting assembly for unstructured mesh generation")
imprinted_assembly, _ = imprint_assembly(
assembly, threads=resolved_imprint_threads, glue=imprint
)
else:
imprinted_assembly = assembly
# gmsh is a global singleton; finalize the session on every exit path
# (including a mid-mesh exception) so repeated calls don't accumulate
# models. gmsh_session_started is only set once init_gmsh() has bound
# the local gmsh name, keeping the finally safe if init_gmsh() itself
# raises. See issue #187.
gmsh_session_started = False
try:
gmsh = init_gmsh()
gmsh_session_started = True
gmsh, volumes_in_model = get_volumes(
gmsh, imprinted_assembly, method=method, scale_factor=scale_factor
)
# Resolve any material tag strings in set_size to volume IDs
resolved_set_size = None
if set_size:
resolved_set_size = resolve_set_size(
set_size, volumes_in_model, self.material_tags
)
gmsh = set_sizes_for_mesh(
gmsh=gmsh,
min_mesh_size=min_mesh_size,
max_mesh_size=max_mesh_size,
mesh_algorithm=mesh_algorithm,
mesh_algorithm_3d=mesh_algorithm_3d,
set_size=resolved_set_size,
original_set_size=set_size,
threads=threads,
)
if volumes:
for volume_id in volumes_in_model:
if volume_id[1] not in volumes:
gmsh.model.occ.remove([volume_id], recursive=True)
gmsh.option.setNumber("Mesh.SaveAll", 1)
gmsh.model.occ.synchronize()
# Clear the mesh
gmsh.model.mesh.clear()
gmsh.option.setNumber(
"Mesh.SaveElementTagType", 3
) # Save only volume elements
_generate_volume_mesh(gmsh, volume_mesh_options)
# makes the folder if it does not exist
if Path(filename).parent:
Path(filename).parent.mkdir(parents=True, exist_ok=True)
# gmsh.write only accepts strings
if isinstance(filename, Path):
gmsh.write(str(filename))
else:
gmsh.write(filename)
return filename
finally:
if gmsh_session_started and gmsh.isInitialized():
gmsh.finalize()
def _export_unstructured_mesh_file_with_mesher(
self,
filename: str,
target_edge_length: float | None,
tet_volumes: Iterable[str] | None,
tolerance: float,
angular_tolerance: float,
imprint: Literal["partial", "full", False] | None,
imprint_threads: int | Literal["all"],
mesh_threads: int | Literal["all"] = "all",
scale_factor: float = 1.0,
) -> str:
"""Write an unstructured .vtk volume mesh using cad-to-dagmc-mesher.
Meshes the assembly with cad-to-dagmc-mesher, combines the per-solid
tetrahedra into a single mesh, and writes it as a legacy VTK file
readable by openmc.UnstructuredMesh(filename, library="moab").
"""
if target_edge_length is None:
raise ValueError(
"target_edge_length is required when meshing_backend is "
'"cad-to-dagmc-mesher"'
)
assembly = _build_assembly(
self.parts, scale_factor, names=_solid_names(self.material_tags)
)
# Default to tetrahedralising every volume. tet_volumes is matched
# against material tags by the mesher, so pass the material tags.
if tet_volumes is None:
tet_volumes = list(self.material_tags)
else:
tet_volumes = list(tet_volumes)
_, _, _, tet_data = _mesh_with_cad_to_dagmc_mesher(
assembly=assembly,
material_tags=self.material_tags,
tolerance=tolerance,
angular_tolerance=angular_tolerance,
tet_volumes=tet_volumes,
target_edge_length=target_edge_length,
imprint=imprint,
imprint_threads=imprint_threads,
mesh_threads=mesh_threads,
)
if not tet_data:
raise ValueError(
"cad-to-dagmc-mesher produced no tetrahedra. Check that "
"tet_volumes contains valid material tags and that "
"target_edge_length is set."
)
tet_vertices, tetrahedra = combine_tet_meshes(tet_data)
if Path(filename).parent:
Path(filename).parent.mkdir(parents=True, exist_ok=True)
write_vtk(filename, tet_vertices, tetrahedra)
print(f"written unstructured mesh file {filename}")
return filename
[docs]
def export_gmsh_mesh_file(
self,
filename: str = "mesh.msh",
min_mesh_size: float | None = None,
max_mesh_size: float | None = None,
mesh_algorithm: int = 1,
dimensions: int = 2,
method: str = "file",
scale_factor: float = 1.0,
imprint: Literal["partial", "full", False] | None = "partial",
set_size: dict[int | str, float] | None = None,
threads: int = 0,
*,
imprint_threads: int | Literal["all"] = "all",
mesh_algorithm_3d: int = 1,
volume_mesh_options: dict[str, float | str] | None = None,
):
"""Saves a GMesh msh file of the geometry in either 2D surface mesh or
3D volume mesh.
Args:
filename
min_mesh_size: the minimum size of mesh elements to use.
max_mesh_size: the maximum size of mesh elements to use.
mesh_algorithm: the Gmsh 2D surface mesh algorithm to use.
mesh_algorithm_3d: Gmsh 3D volume mesh algorithm, passed to
Mesh.Algorithm3D. Defaults to 1 (Delaunay); 10 selects HXT.
volume_mesh_options: Gmsh option names mapped to numeric or string
values, applied after surface meshing and before volume meshing.
Overrides the corresponding sizing/algorithm options for the
volume step only. Ignored unless dimensions=3. None or an empty
dict keeps the single generate(dimensions) call.
dimensions: The number of dimensions, 2 for a surface mesh 3 for a
volume mesh. Passed to gmsh.model.mesh.generate()
method: the method to use to import the geometry into gmsh. Options
are 'file' or 'in memory'. 'file' is the default and will write
the geometry to a temporary file before importing it into gmsh.
'in memory' will import the geometry directly into gmsh but
requires the version of OpenCASCADE used to build gmsh to be
the same as the version used by CadQuery. This is possible to
ensure when installing the package with Conda but harder when
installing from PyPI.
scale_factor: a scaling factor to apply to the geometry that can be
used to enlarge or shrink the geometry. Useful when converting
the geometry to cm for use in neutronics.
imprint: "partial" (default) for non-overlapping contact, "full" for
exactly coincident whole faces, None to imprint without glue, or
False to skip imprinting entirely. None does not mean the default.
CadQuery 2.8.0 supports only partial glue; other modes raise
RuntimeError before imprinting multiple solids. A single solid
skips imprinting, but arguments are always validated. True,
integers, and dictionaries are not accepted.
imprint_threads: keyword-only thread limit for imprinting: "all"
(default, all OCC logical processors) or a positive int.
Fewer threads reduce peak RAM at the cost of speed. The previous
OCC pool is restored before meshing, including after errors.
Validated even when imprint=False; independent of Gmsh threads.
set_size: a dictionary mapping volume IDs (int) or material tag names
(str) to target mesh sizes (floats). Material tags are resolved to
all volume IDs that have that tag.
threads: the number of threads for Gmsh to use. 0 uses all
available cores (default), 1 uses a single thread.
"""
resolved_imprint_threads = resolve_imprint(imprint, imprint_threads)
assembly = cq.Assembly()
for part in self.parts:
assembly.add(part)
if imprint is not False:
print("Imprinting assembly for mesh generation")
imprinted_assembly, _ = imprint_assembly(
assembly, threads=resolved_imprint_threads, glue=imprint
)
else:
imprinted_assembly = assembly
# gmsh is a global singleton; finalize the session on every exit path
# (including a mid-mesh exception) so repeated calls don't accumulate
# models. gmsh_session_started is only set once init_gmsh() has bound
# the local gmsh name, keeping the finally safe if init_gmsh() itself
# raises. See issue #187.
gmsh_session_started = False
try:
gmsh = init_gmsh()
gmsh_session_started = True
gmsh, volumes = get_volumes(
gmsh, imprinted_assembly, method=method, scale_factor=scale_factor
)
# Resolve any material tag strings in set_size to volume IDs
resolved_set_size = None
if set_size:
resolved_set_size = resolve_set_size(
set_size, volumes, self.material_tags
)
gmsh = set_sizes_for_mesh(
gmsh=gmsh,
min_mesh_size=min_mesh_size,
max_mesh_size=max_mesh_size,
mesh_algorithm=mesh_algorithm,
mesh_algorithm_3d=mesh_algorithm_3d,
set_size=resolved_set_size,
original_set_size=set_size,
threads=threads,
)
if dimensions == 3:
_generate_volume_mesh(gmsh, volume_mesh_options)
else:
gmsh.model.mesh.generate(dimensions)
# makes the folder if it does not exist
if Path(filename).parent:
Path(filename).parent.mkdir(parents=True, exist_ok=True)
# gmsh.write only accepts strings
if isinstance(filename, Path):
gmsh.write(str(filename))
else:
gmsh.write(filename)
print(f"written GMSH mesh file {filename}")
finally:
if gmsh_session_started and gmsh.isInitialized():
gmsh.finalize()
[docs]
def export_dagmc_h5m_file(
self,
filename: str = "dagmc.h5m",
implicit_complement_material_tag: str | None = None,
scale_factor: float = 1.0,
imprint: Literal["partial", "full", False] | None = "partial",
*,
imprint_threads: int | Literal["all"] = "all",
mesh_threads: int | Literal["all"] = "all",
**kwargs,
) -> str:
"""Saves a DAGMC h5m file of the geometry
Args:
filename: the filename to use for the saved DAGMC file.
implicit_complement_material_tag: the name of the material tag to use
for the implicit complement (void space).
scale_factor: a scaling factor to apply to the geometry. All the
linear mesh sizing arguments (min_mesh_size, max_mesh_size and
set_size for gmsh, tolerance and target_edge_length for
cad-to-dagmc-mesher) are in the units of
the SCALED geometry, so the same number means the same thing on
the output mesh whichever backend is used. For example with
scale_factor=100 (m to cm) a tolerance of 0.5 is a 5 mm
deflection. Note this means the defaults get finer as
scale_factor grows: the cad-to-dagmc-mesher tolerance default of
0.01 is a 0.1 mm deflection at scale_factor=100, which on a large
model can produce a very fine mesh and exhaust memory, so scale
the tolerance along with the geometry. angular_tolerance is an
angle and so is unaffected by scaling.
imprint: "partial" (default) for non-overlapping contact, "full" for
exactly coincident whole faces, None to imprint without glue, or
False to skip imprinting entirely. None does not mean the default.
CadQuery 2.8.0 supports only partial glue; other modes raise
RuntimeError before imprinting multiple solids. A single solid
skips imprinting, but arguments are always validated. True,
integers, and dictionaries are not accepted.
imprint_threads: keyword-only thread limit for imprinting: "all"
(default, all OCC logical processors) or a positive int.
Fewer threads reduce peak RAM at the cost of speed. The previous
OCC pool is restored before meshing, including after errors.
Validated even when imprint=False; independent of Gmsh threads.
mesh_threads: keyword-only thread count for the cad-to-dagmc-mesher
meshing that follows the imprint: "all" (default) or a positive
int. Fewer threads reduce peak RAM at the cost of speed, and
give the same mesh. Validated with any backend; the gmsh
backend warns that it is unused (it takes threads instead).
**kwargs: Backend-specific parameters:
Backend selection:
- meshing_backend (str, optional): explicitly specify 'gmsh' or
'cad-to-dagmc-mesher'. If not provided, backend is
auto-selected based on other arguments: cad-to-dagmc-mesher
arguments select 'cad-to-dagmc-mesher', gmsh-specific
arguments select 'gmsh'. Defaults to 'cad-to-dagmc-mesher' if
no backend-specific arguments are given, falling back to
'gmsh' when cad-to-dagmc-mesher is not installed.
- h5m_backend (str, optional): 'pymoab' or 'h5py' for writing h5m files.
Defaults to 'h5py'.
For GMSH backend:
- min_mesh_size (float): minimum mesh element size
- max_mesh_size (float): maximum mesh element size
- mesh_algorithm (int): GMSH 2D surface mesh algorithm (default: 1)
- mesh_algorithm_3d (int): GMSH 3D volume mesh algorithm, passed
to Mesh.Algorithm3D (default: 1, Delaunay; 10 selects HXT).
- volume_mesh_options (dict[str, float | str]): Gmsh options
applied after surface meshing and before volume meshing.
Overrides the corresponding sizing/algorithm options for the
volume step only. Only used when unstructured_volumes is set.
Defaults to None (no overrides).
- method (str): import method 'file' or 'in memory' (default: 'file')
- set_size (dict[int | str, float]): volume IDs (int) or material tag
names (str) mapped to target mesh sizes. Material tags are resolved
to all volume IDs that have that tag.
- unstructured_volumes (Iterable[int | str]): volume IDs (int) or material
tag names (str) for unstructured mesh. Material tags are resolved to
all volume IDs that have that tag. Can mix ints and strings.
- umesh_filename (str): filename for unstructured mesh (default: 'umesh.vtk')
- threads (int): number of threads for Gmsh to use. 0 uses all
available cores (default), 1 uses a single thread.
For cad-to-dagmc-mesher backend:
- tolerance (float): surface meshing tolerance in cm, after
scale_factor (default: 0.5, a 5 mm deflection).
- angular_tolerance (float): surface angular tolerance (default: 0.2)
- tet_volumes (Iterable[str]): material tag names of the volumes to
fill with tetrahedra for an unstructured volume mesh.
- target_edge_length (float): target tetrahedron edge length. Both
tet_volumes and target_edge_length must be given together to write
a volume mesh; when they are, the return value is a
(dagmc_filename, umesh_filename) tuple.
- umesh_filename (str): filename for the unstructured volume mesh
(default: 'umesh.vtk').
After meshing, ``self.refinement`` holds the mesher's
auto-refine report (rounds, convergence, remaining folds,
final per-solid tolerances), or None if no refinement ran.
Returns:
str: the filename(s) for the files created.
Raises:
ValueError: If invalid parameter combinations are used.
"""
resolved_imprint_threads = resolve_imprint(imprint, imprint_threads)
resolve_mesh_threads(mesh_threads)
# Define all acceptable kwargs
gmsh_keys = {
"min_mesh_size",
"max_mesh_size",
"mesh_algorithm",
"mesh_algorithm_3d",
"volume_mesh_options",
"set_size",
"umesh_filename",
"method",
"unstructured_volumes",
"threads",
}
cad_to_dagmc_mesher_keys = {"tolerance", "angular_tolerance", "tet_volumes", "target_edge_length"}
all_acceptable_keys = gmsh_keys | cad_to_dagmc_mesher_keys | {"meshing_backend", "h5m_backend"}
# A report from an earlier export must not outlive it.
self.refinement = None
# Check for invalid kwargs
invalid_keys = set(kwargs.keys()) - all_acceptable_keys
if invalid_keys:
raise ValueError(
f"Invalid keyword arguments: {sorted(invalid_keys)}\n"
f"Acceptable arguments are: {sorted(all_acceptable_keys)}"
)
# Handle meshing_backend - either from kwargs or auto-detect
meshing_backend = kwargs.pop("meshing_backend", None)
# Handle h5m_backend - pymoab or h5py
h5m_backend = kwargs.pop("h5m_backend", "h5py")
if meshing_backend is None:
# Auto-select meshing_backend based on kwargs.
# umesh_filename is accepted by both the gmsh and the
# cad-to-dagmc-mesher backends, so it is not gmsh specific, but it
# is still part of gmsh_keys because it selects gmsh when nothing
# else narrows the choice.
gmsh_only_keys = gmsh_keys - {"umesh_filename"}
has_gmsh = any(key in kwargs for key in gmsh_keys)
has_mesher = any(key in kwargs for key in cad_to_dagmc_mesher_keys)
if has_mesher:
provided_gmsh = [key for key in sorted(gmsh_only_keys) if key in kwargs]
if provided_gmsh:
provided_mesher = [
key for key in sorted(cad_to_dagmc_mesher_keys) if key in kwargs
]
raise ValueError(
"Ambiguous backend: both cad-to-dagmc-mesher and GMSH-specific arguments provided.\n"
f"cad-to-dagmc-mesher-specific arguments: {sorted(cad_to_dagmc_mesher_keys)}\n"
f"GMSH-specific arguments: {sorted(gmsh_only_keys)}\n"
f"Provided cad-to-dagmc-mesher arguments: {provided_mesher}\n"
f"Provided GMSH arguments: {provided_gmsh}\n"
"Please provide only one backend's arguments."
)
if not _cad_to_dagmc_mesher_is_available():
provided_mesher = [
key for key in sorted(cad_to_dagmc_mesher_keys) if key in kwargs
]
raise CadToDagmcMesherNotFoundError(
f"The arguments {provided_mesher} select the "
"cad-to-dagmc-mesher meshing backend, but "
"cad-to-dagmc-mesher is not installed. It is not available "
"on conda-forge so it has to be installed separately.\n\n"
"Either install it:\n"
" pip install cad-to-dagmc-mesher\n\n"
"or use the gmsh backend, which sizes the mesh with its own "
"arguments such as min_mesh_size and max_mesh_size:\n"
" export_dagmc_h5m_file(..., meshing_backend='gmsh', "
"max_mesh_size=...)"
)
meshing_backend = "cad-to-dagmc-mesher"
elif has_gmsh:
meshing_backend = "gmsh"
elif _cad_to_dagmc_mesher_is_available():
meshing_backend = "cad-to-dagmc-mesher" # default
else:
# cad-to-dagmc-mesher is a dependency of the pip package but is
# not on conda-forge, so a conda installation can be without it.
# A call that names no backend at all has expressed no
# preference, so fall back to gmsh, which is always present,
# rather than failing. Warn so the substitution is visible and
# so the remedy is to hand, since pip installing the mesher
# works alongside a conda installation.
warnings.warn(
"No meshing backend was given so the cad-to-dagmc-mesher "
"backend would be used, but cad-to-dagmc-mesher is not "
"installed. Falling back to the gmsh backend. "
"cad-to-dagmc-mesher is not available on conda-forge, "
"install it with pip, which works alongside a conda "
"installation:\n"
" pip install cad-to-dagmc-mesher\n\n"
"Pass meshing_backend='gmsh' to select the gmsh "
"backend explicitly and silence this warning."
)
meshing_backend = "gmsh"
# Validate meshing backend
if meshing_backend not in ["gmsh", "cad-to-dagmc-mesher"]:
raise ValueError(
f'meshing_backend "{meshing_backend}" not supported. '
'Available options are "gmsh" or "cad-to-dagmc-mesher"'
)
print(f"Using meshing backend: {meshing_backend}")
# Initialize variables to avoid unbound errors
tolerance = 0.1
angular_tolerance = 0.1
min_mesh_size = None
max_mesh_size = None
mesh_algorithm = 1
mesh_algorithm_3d = 1
volume_mesh_options = None
method = "file"
set_size = None
unstructured_volumes = None
umesh_filename = "umesh.vtk"
threads = 0
tet_data = None
# Extract backend-specific parameters with defaults
if meshing_backend == "gmsh":
# GMSH parameters
min_mesh_size = kwargs.get("min_mesh_size")
max_mesh_size = kwargs.get("max_mesh_size")
mesh_algorithm = kwargs.get("mesh_algorithm", 1)
mesh_algorithm_3d = kwargs.get("mesh_algorithm_3d", 1)
volume_mesh_options = kwargs.get("volume_mesh_options")
method = kwargs.get("method", "file")
set_size = kwargs.get("set_size")
unstructured_volumes = kwargs.get("unstructured_volumes")
umesh_filename = kwargs.get("umesh_filename", "umesh.vtk")
threads = kwargs.get("threads", 0)
# Warn about unused cad-to-dagmc-mesher parameters
non_gmsh_params = [
"tolerance",
"angular_tolerance",
"tet_volumes",
"target_edge_length",
]
unused_params = [param for param in non_gmsh_params if param in kwargs]
if mesh_threads != "all":
unused_params.append("mesh_threads")
if unused_params:
warnings.warn(
f"The following parameters are ignored when using GMSH backend: "
f"{', '.join(unused_params)}"
)
elif meshing_backend == "cad-to-dagmc-mesher":
tolerance = kwargs.get("tolerance", DEFAULT_TOLERANCE)
angular_tolerance = kwargs.get("angular_tolerance", DEFAULT_ANGULAR_TOLERANCE)
unused_params = [
param
for param in ("mesh_algorithm_3d", "volume_mesh_options")
if param in kwargs
]
if unused_params:
warnings.warn(
"The following parameters are ignored when using "
f"cad-to-dagmc-mesher backend: {', '.join(unused_params)}"
)
assembly = cq.Assembly()
for part in self.parts:
assembly.add(part)
original_ids = get_ids_from_assembly(assembly)
# both id lists should be the same length as each other and the same
# length as the self.material_tags
if len(original_ids) != len(self.material_tags):
msg = f"Number of volumes {len(original_ids)} is not equal to number of material tags {len(self.material_tags)}"
raise ValueError(msg)
# The gmsh backend opens a gmsh session (gmsh is a global singleton).
# Wrap the whole meshing and export in try/finally so the session is
# always finalized - on every return path and even if meshing raises
# part way through. Without this, repeated calls accumulate gmsh models
# in the session (see issue #187). gmsh_session_started is only set once
# init_gmsh() has run, so the finally never touches the (function-local)
# gmsh name before it is bound and never finalizes a session the caller
# may own when using a non-gmsh backend.
gmsh_session_started = False
try:
# Use gmsh
if meshing_backend == "gmsh":
# If assembly is not to be imprinted, pass through the assembly as-is
if imprint is not False:
print("Imprinting assembly for mesh generation")
imprinted_assembly, imprinted_solids_with_org_id = (
imprint_assembly(
assembly, threads=resolved_imprint_threads, glue=imprint
)
)
scrambled_ids = get_ids_from_imprinted_assembly(
imprinted_solids_with_org_id
)
material_tags_in_brep_order = order_material_ids_by_brep_order(
original_ids, scrambled_ids, self.material_tags
)
else:
material_tags_in_brep_order = self.material_tags
imprinted_assembly = assembly
check_material_tags(material_tags_in_brep_order, self.parts)
# Start generating the mesh
gmsh = init_gmsh()
gmsh_session_started = True
gmsh, volumes = get_volumes(
gmsh, imprinted_assembly, method=method, scale_factor=scale_factor
)
# Resolve any material tag strings in set_size to volume IDs
resolved_set_size = None
if set_size:
resolved_set_size = resolve_set_size(
set_size, volumes, material_tags_in_brep_order
)
gmsh = set_sizes_for_mesh(
gmsh=gmsh,
min_mesh_size=min_mesh_size,
max_mesh_size=max_mesh_size,
mesh_algorithm=mesh_algorithm,
mesh_algorithm_3d=mesh_algorithm_3d,
set_size=resolved_set_size,
original_set_size=set_size,
threads=threads,
)
gmsh.model.mesh.generate(2)
vertices, triangles_by_solid_by_face = mesh_to_vertices_and_triangles(
dims_and_vol_ids=volumes
)
elif meshing_backend == "cad-to-dagmc-mesher":
tet_volumes_arg = kwargs.get("tet_volumes", kwargs.get("unstructured_volumes"))
target_edge_length = kwargs.get("target_edge_length")
umesh_filename = kwargs.get("umesh_filename", umesh_filename)
# A volume (tet) mesh needs BOTH tet_volumes and
# target_edge_length. Passing only one (or asking for a
# umesh_filename without them) is a user error: fail fast with a
# clear message rather than silently writing no .vtk and
# returning a bare string instead of the (h5m, vtk) tuple.
wants_umesh = (
bool(tet_volumes_arg)
or target_edge_length is not None
or "umesh_filename" in kwargs
)
if wants_umesh and not (tet_volumes_arg and target_edge_length):
raise ValueError(
"Writing an unstructured volume mesh with the "
"cad-to-dagmc-mesher backend requires BOTH tet_volumes "
"(material tag names) and target_edge_length. Got "
f"tet_volumes={tet_volumes_arg!r}, "
f"target_edge_length={target_edge_length!r}."
)
# scale_factor is applied to the geometry before meshing so the
# h5m and .vtk match the gmsh backend (which scales).
mesher_assembly = _build_assembly(
self.parts, scale_factor,
names=_solid_names(self.material_tags),
)
vertices, triangles_by_solid_by_face, material_tags_in_brep_order, tet_data = (
_mesh_with_cad_to_dagmc_mesher(
assembly=mesher_assembly,
material_tags=self.material_tags,
tolerance=tolerance,
angular_tolerance=angular_tolerance,
tet_volumes=tet_volumes_arg,
target_edge_length=target_edge_length,
imprint=imprint,
imprint_threads=imprint_threads,
mesh_threads=mesh_threads,
report=self,
)
)
else:
raise ValueError(
f'meshing_backend {meshing_backend} not supported. '
'Available options are "gmsh" or "cad-to-dagmc-mesher"'
)
dagmc_filename = vertices_to_h5m(
vertices=vertices,
triangles_by_solid_by_face=triangles_by_solid_by_face,
material_tags=material_tags_in_brep_order,
h5m_filename=filename,
implicit_complement_material_tag=implicit_complement_material_tag,
method=h5m_backend,
)
if meshing_backend == "gmsh" and unstructured_volumes:
# Resolve any material tag strings to volume IDs
unstructured_volumes = resolve_unstructured_volumes(
unstructured_volumes, volumes, material_tags_in_brep_order
)
# remove all the unused occ volumes, this prevents them being meshed
for volume_dim, volume_id in volumes:
if volume_id not in unstructured_volumes:
gmsh.model.occ.remove(
[(volume_dim, volume_id)], recursive=True
)
gmsh.option.setNumber("Mesh.SaveAll", 1)
gmsh.model.occ.synchronize()
# removes all the 2D groups so that 2D faces are not included in the vtk file
all_2d_groups = gmsh.model.getPhysicalGroups(2)
for entry in all_2d_groups:
gmsh.model.removePhysicalGroups([entry])
_generate_volume_mesh(gmsh, volume_mesh_options, surfaces_meshed=True)
gmsh.option.setNumber(
"Mesh.SaveElementTagType", 3
) # Save only volume elements
gmsh.write(umesh_filename)
return dagmc_filename, umesh_filename
# The cad-to-dagmc-mesher backend produces the tetrahedra itself
# (when tet_volumes + target_edge_length are given). Combine the
# per-solid tet meshes and write a .vtk unstructured volume mesh
# without going through gmsh. Keying on the user's request (both
# tet args, guaranteed present together by the check above) rather
# than on tet_data means a mesher that unexpectedly yields no tets
# raises here instead of silently returning a bare string.
if meshing_backend == "cad-to-dagmc-mesher" and tet_volumes_arg and target_edge_length:
if not tet_data:
raise ValueError(
"cad-to-dagmc-mesher produced no tetrahedra despite "
f"tet_volumes={tet_volumes_arg!r} and "
f"target_edge_length={target_edge_length!r}. Check that "
"tet_volumes contains valid material tags."
)
tet_vertices, tetrahedra = combine_tet_meshes(tet_data)
if Path(umesh_filename).parent:
Path(umesh_filename).parent.mkdir(parents=True, exist_ok=True)
write_vtk(umesh_filename, tet_vertices, tetrahedra)
print(f"written unstructured mesh file {umesh_filename}")
return dagmc_filename, umesh_filename
return dagmc_filename
finally:
if gmsh_session_started and gmsh.isInitialized():
gmsh.finalize()
def _scale_together(shapes, scale_factor: float = 1.0) -> list:
"""Scale shapes as one, keeping any faces they share shared.
Solids that were imprinted against each other share the face between them,
and the mesher relies on that sharing to give the face a sense for both
volumes. Scaling each solid on its own transforms that face once per solid,
so every solid ends up with its own copy and the interface becomes two
coincident one sided DAGMC surfaces. Scaling them together as one compound
transforms each shared face once.
Shape.scale returns a new shape (it does not mutate in place), so the
shapes passed in are left untouched.
Args:
shapes: the shapes to scale.
scale_factor: the scaling factor to apply.
Returns:
The scaled shapes, in the same order as the shapes passed in.
"""
shapes = list(shapes)
if scale_factor == 1.0:
return shapes
scaled = cq.Compound.makeCompound(shapes).scale(scale_factor)
# The direct children of the compound, rather than Solids(), so that every
# shape passed in comes back as exactly one shape in the same position.
children = []
iterator = TopoDS_Iterator(scaled.wrapped)
while iterator.More():
children.append(cq.Shape.cast(iterator.Value()))
iterator.Next()
return children
def _build_assembly(parts, scale_factor: float = 1.0, names=None):
"""Build a CadQuery assembly from parts, optionally scaling them.
The parts are scaled together with _scale_together so faces they share
stay shared, and the original parts in self.parts are left untouched so
repeated exports stay consistent.
names, when given, labels each child. cad-to-dagmc-mesher's SolidConfig
addresses solids by assembly child name, so naming them is what lets us ask
for per-solid meshing rather than the positional material_tags API.
"""
assembly = cq.Assembly()
for index, scaled in enumerate(_scale_together(parts, scale_factor)):
if names is None:
assembly.add(scaled)
else:
assembly.add(scaled, name=names[index])
return assembly
def _solid_names(material_tags) -> list[str]:
"""Unique per-solid names for the mesher's SolidConfig.
Material tags cannot be used directly: several solids commonly share one tag
and SolidConfig needs a distinct name per solid. Suffixing the index keeps
them unique while staying readable, which matters because the mesher prints
these names when it reports what it refined.
The name also carries the input index, so the real tags can be reattached
afterwards whatever order the mesher returns solids in.
"""
return [f"{tag}#{i}" for i, tag in enumerate(material_tags)]
def _mesh_with_cad_to_dagmc_mesher(
assembly, material_tags, tolerance, angular_tolerance,
tet_volumes, target_edge_length, imprint, imprint_threads="all",
report=None, mesh_threads="all",
):
"""Mesh using cad-to-dagmc-mesher and return vertices_to_h5m-compatible output.
Returns ``(vertices, triangles_by_solid_by_face, material_tags,
tet_data)`` where ``tet_data`` is the per-solid tetrahedral mesh dict
(``{solid_id: {"vertices": ..., "tetrahedra": ..., ...}}``) or ``None``
when no solids were volume-meshed. Volume meshing only happens when both
``tet_volumes`` and ``target_edge_length`` are supplied.
When ``report`` is given, its ``refinement`` attribute is set to the
mesher's auto-refine report.
"""
try:
from cad_to_dagmc_mesher.cad import (
OverlappingSolidsError,
SolidConfig,
mesh_assembly,
)
except ImportError as e:
raise CadToDagmcMesherNotFoundError() from e
# Address solids individually rather than through the positional
# material_tags API. Two reasons:
#
# - the mesher only refines per solid, and only guards against refinement
# making things worse, on its solid_config path. The material_tags path
# refines the whole assembly and returns whatever the last round produced,
# so on geometry whose fold count grows with refinement it can hand back a
# mesh worse than the one it started from.
# - tags are reattached by name here instead of by position, so the mapping
# survives the mesher returning solids in a different order.
names = _solid_names(material_tags)
tag_by_name = dict(zip(names, material_tags))
tet_tags = set(tet_volumes or ())
configs = [
SolidConfig(
name=name,
tolerance=tolerance,
angular_tolerance=angular_tolerance,
target_edge_length=(
target_edge_length
if target_edge_length is not None and tag_by_name[name] in tet_tags
else None
),
)
for name in names
]
# The mesher owns its imprint and meshing thread scopes; do not limit
# anything here. A mesher older than mesh_threads meshes on all threads,
# so it is only passed when set, and a version that cannot honour it says
# so rather than ignoring it.
threads = {}
if mesh_threads != "all":
if "mesh_threads" not in inspect.signature(mesh_assembly).parameters:
raise RuntimeError(
"The installed cad-to-dagmc-mesher cannot limit its meshing "
f"threads, so mesh_threads={mesh_threads!r} cannot be honoured. "
"Leave mesh_threads='all', or upgrade cad-to-dagmc-mesher to a "
"version that takes mesh_threads."
)
threads["mesh_threads"] = mesh_threads
try:
result = mesh_assembly(
assembly,
solid_config=configs,
imprint=imprint,
imprint_threads=imprint_threads,
**threads,
)
refinement = result.get("refinement")
if report is not None:
report.refinement = refinement
return (
result["vertices"],
result["triangles_by_solid_by_face"],
[tag_by_name[name] for name in result["material_tags"]],
result.get("tet_data"),
)
except OverlappingSolidsError as e:
# Imprinting fused two or more solids together, which means they overlap.
# That geometry cannot be written as valid DAGMC: a region inside two
# volumes has no single material, and a DAGMC surface separates at most
# two volumes. It is better to say so than to mesh it, because the result
# is silently wrong rather than obviously broken -- it is watertight, it
# transports, and the material in the shared region is whichever volume
# DAGMC happens to resolve first.
#
# The mesher names the solids it could not find, but those are the
# synthetic per-solid names from _solid_names, so they are re-expressed
# here as the material tags the caller actually supplied.
fused = sorted({
tag_by_name[name]
for name in names
if name in str(e) and name in tag_by_name
})
detail = (
f" The solids tagged {fused} could not be told apart after imprinting."
if fused else ""
)
raise ValueError(
"Overlapping solids in the CAD." + detail + " Overlapping geometry is "
"invalid for DAGMC: a region inside two volumes has no single "
"material, so transport through it is not well defined. Resolve the "
"overlaps before meshing, for example by boolean-subtracting the "
"inner solid from the outer one "
"(outer = outer.cut(inner)) so the two share a surface instead of a "
"volume."
) from e
def _get_all_leaf_children(assembly):
"""Recursively yield all leaf children (parts, not assemblies) from a CadQuery assembly."""
for child in assembly.children:
# If the child is itself an assembly, recurse
if hasattr(child, "children") and len(child.children) > 0:
yield from _get_all_leaf_children(child)
else:
yield child
