# Copyright 2026, UChicago Argonne, LLC
# All Rights Reserved
# Software Name: graintrace
# By: Argonne National Laboratory
# OPEN SOURCE LICENSE (MIT)
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
"""Reconstruct FF-HEDM microstructures via NEPER Voronoi tessellation (VoronoiMeshBuilder)."""
from __future__ import annotations
import os
import shutil
import subprocess
from pathlib import Path
import numpy as np
import pandas as pd
[docs]
class VoronoiMeshBuilder:
"""Reconstruct a 3D microstructure from far-field (FF) HEDM grain centroids.
Wraps NEPER (Voronoi/CVT tessellation) and optionally GMSH to turn a per-grain FF
CSV (positions, Euler orientations, elastic strains) into a ``.tess``/``.msh`` mesh,
per-grain orientation and initial elastic-strain files, and a voxelized
reconstruction: the inputs for a CPFE run or a grain graph. See
``examples/demonstrate_farfield.py`` and the ``/ff-reconstruction`` skill.
"""
def __init__(
self,
input_csv,
output_dir,
bounding_box,
dim=3,
weighted=False,
auto_fix_bbox=False,
bbox_fix_mode=None,
bbox_tolerance=0.0, # percentage tolerance (e.g., 5.0 = 5%)
auto_rotate=False,
rotate_angles=(0, 0, 0), # rotation angles (rotate around X -> Y -> Z)
rotate_convention="xyz",
angle_identifier=None,
orientation_descriptor="euler-bunge",
orientation_active_convention=False,
unit="deg", # rotation unit ('deg' or 'rad')
elastic_strain_identifier=None,
ori_rotmat_identifier=None,
strain_unit="microstrain",
env=None,
neper_path=None,
auto_install=False,
):
"""Configure an FF-HEDM Voronoi reconstruction from a per-grain CSV.
Args:
input_csv: Path to the FF grain CSV (columns X, Y[, Z], optional GrainRadius, plus the angle and elastic-strain columns named below).
output_dir: Directory where NEPER inputs, the tessellation, orientation/strain files, and the voxel reconstruction are written.
bounding_box: Reconstruction domain as [xlo, xhi, ylo, yhi] (2D) or [xlo, xhi, ylo, yhi, zlo, zhi] (3D) in micrometers.
dim: Spatial dimension, 2 or 3.
weighted: If True, use a Laguerre (weighted) tessellation with per-grain weights from GrainRadius (requires a GrainRadius column).
auto_fix_bbox: If True, reconcile the data with bounding_box using bbox_fix_mode before tessellating.
bbox_fix_mode: How to reconcile when auto_fix_bbox is True: "remove_points" (drop out-of-box points, production) or "extend_bounding_box" (grow the box, debug).
bbox_tolerance: Percentage tolerance used by the bbox fix (e.g. 5.0 = 5%).
auto_rotate: If True, align the point cloud to its PCA principal axes instead of applying rotate_angles.
rotate_angles: Explicit sample-tilt correction as xyz Euler angles applied in order X -> Y -> Z; units set by unit.
rotate_convention: Euler rotation convention for rotate_angles; only "xyz" is supported.
angle_identifier: The 3 column names holding the per-grain Euler angles, or None to skip orientation parsing.
orientation_descriptor: NEPER orientation descriptor for the input angles (e.g. "euler-bunge").
orientation_active_convention: If True, treat orientations as the active convention when applying the rotation matrix.
unit: Angle/rotation unit, "deg" or "rad"; must match the actual data units.
elastic_strain_identifier: The 9 column names (row-major) of the per-grain elastic strain tensor, or None to skip strain.
ori_rotmat_identifier: The 9 column names for the rotated orientation matrix written into self.data; defaults to O11..O33.
strain_unit: Unit of the elastic strain columns, "microstrain" or "strain".
env: Pre-built subprocess environment for NEPER; if None it is resolved automatically via check_dependencies.
neper_path: Explicit path to the NEPER binary; if None it is resolved from the NEPER env var, tools.json, or PATH.
auto_install: If True, attempt an opt-in local ~/.local NEPER source build when NEPER cannot be found.
"""
self.input_csv = input_csv
self.output_dir = output_dir
self.bounding_box = bounding_box
self.dim = dim
self.weighted = weighted
self.auto_fix_bbox = auto_fix_bbox
self.bbox_fix_mode = bbox_fix_mode
self.bbox_tolerance = bbox_tolerance
self.rotate_convention = rotate_convention
self.ori_descriptor = orientation_descriptor
self.angle_id = angle_identifier
self.rotate_matrix = np.eye(3)
self.auto_rotate = auto_rotate
self.orientation_active_convention = orientation_active_convention
self.data = None
self.elastic_strain_id = elastic_strain_identifier
self.strain_unit = strain_unit
self.ori_rotmat_id = ori_rotmat_identifier
if env is None:
self.neper_bin, self.env = self.check_dependencies(
neper_path=neper_path, auto_install=auto_install
)
else:
self.env = env
self.neper_bin = neper_path or "neper"
if unit not in ("deg", "rad"):
raise ValueError("Rotation unit must be 'deg' or 'rad'.")
self.rotate_angles = rotate_angles
self.unit = unit
if self.rotate_convention != "xyz":
raise ValueError("currently only rotate_convention = 'xyz' is supported.")
valid_modes = ["remove_points", "extend_bounding_box"]
if self.auto_fix_bbox:
if self.bbox_fix_mode not in valid_modes:
raise ValueError(
f"Invalid bbox_fix_mode='{self.bbox_fix_mode}'. "
f"When auto_fix_bbox=True, bbox_fix_mode must be one of {valid_modes}."
)
else:
if self.bbox_fix_mode is not None:
print(
f"Warning: bbox_fix_mode='{self.bbox_fix_mode}' "
f"is ignored since auto_fix_bbox=False."
)
if self.elastic_strain_id is not None:
if len(self.elastic_strain_id) != 9:
raise ValueError(
"elastic_strain_identifier must contain exactly 9 components."
)
if self.strain_unit not in ("microstrain", "strain"):
raise ValueError(
f"Invalid strain_unit '{self.strain_unit}'. Must be 'microstrain' or 'strain'."
)
if self.ori_rotmat_id is not None:
if len(self.ori_rotmat_id) != 9:
raise ValueError(
"ori_rotmat_identifier must contain exactly 9 components."
)
else:
self.ori_rotmat_id = [
"O11",
"O12",
"O13",
"O21",
"O22",
"O23",
"O31",
"O32",
"O33",
]
if self.angle_id is not None:
if not isinstance(self.angle_id, (list, tuple)) or len(self.angle_id) != 3:
raise ValueError(
"angle_identifier must contain 3 elements if provided."
)
if self.dim not in [2, 3]:
raise ValueError("Dimension 'dim' must be either 2 or 3.")
[docs]
def check_dependencies(self, neper_path=None, auto_install=False):
"""Locate NEPER (user-installed) and build its subprocess environment.
gmsh ships with graintrace via pip, so no gmsh handling is needed here.
Returns ``(neper_bin, env)``. Raises ``RuntimeError`` with install
guidance when NEPER cannot be found and ``auto_install`` is False.
"""
# pylint: disable=import-outside-toplevel # avoid an import cycle
from graintrace.neper_env import resolve_neper_env
neper_bin, env = resolve_neper_env(
neper_path=neper_path, auto_install=auto_install
)
self.env = env
return neper_bin, env
[docs]
def rotate(self):
"""Align data to PCA principal axes (auto_rotate) or apply explicit xyz Euler rotation."""
if self.data is None:
raise RuntimeError("No data loaded. Call read_input() first.")
df = self.data.copy()
# Automatic alignment using PCA
if self.auto_rotate:
print("Auto-rotation enabled: aligning data to principal axes.")
eigenvectors, _ = self.compute_principal_axes()
coords = (
df[["X", "Y", "Z"]].to_numpy()
if self.dim == 3
else df[["X", "Y"]].to_numpy()
)
mean = np.mean(coords, axis=0)
centered = coords - mean
rotated = np.dot(centered, eigenvectors)
# Ensure right-handed orientation
if np.linalg.det(eigenvectors) < 0:
eigenvectors[:, 2] *= -1
rotated = np.dot(centered, eigenvectors)
# Align PCA axes' signs to global XYZ
ref_axis = np.eye(3)
signs = np.sign(np.diag(np.dot(eigenvectors.T, ref_axis)))
eigenvectors *= signs
rotated = np.dot(centered, eigenvectors)
if self.dim == 3:
df[["X", "Y", "Z"]] = rotated
else:
df[["X", "Y"]] = rotated[:, :2]
self.data = df
self.rotate_matrix = eigenvectors
print("Data rotated to align with principal axes.\n")
return
# Manual rotation using scipy
print("Manual rotation:")
# pylint: disable=import-outside-toplevel # scipy imported lazily
from scipy.spatial.transform import Rotation as R
ang = np.array(self.rotate_angles, dtype=float)
if self.unit == "deg":
ang = np.radians(ang)
if self.rotate_convention == "xyz":
rot = R.from_euler("xyz", ang, degrees=False)
self.rotate_matrix = rot.as_matrix()
if self.dim == 3:
df[["X", "Y", "Z"]] = np.dot(
df[["X", "Y", "Z"]].to_numpy(), self.rotate_matrix.T
)
else:
df[["X", "Y"]] = np.dot(
df[["X", "Y"]].to_numpy(), self.rotate_matrix[:2, :2].T
)
self.data = df
print(f"Rotated dataset by rotate_angles {self.rotate_angles} ({self.unit}).\n")
[docs]
def compute_principal_axes(self):
"""Compute PCA principal axes (2D/3D, right-handed). Returns (eigenvectors, eigenvalues)."""
# pylint: disable=import-outside-toplevel # sklearn is a heavy optional dep
from sklearn.decomposition import PCA
if self.data is None or self.data.empty:
raise RuntimeError("Data not loaded. Call read_input() first.")
df = self.data.copy()
coords = (
df[["X", "Y"]].to_numpy()
if self.dim == 2
else df[["X", "Y", "Z"]].to_numpy()
)
mean = np.mean(coords, axis=0)
centered = coords - mean
pca = PCA(n_components=self.dim)
pca.fit(centered)
eigenvectors = pca.components_.T # columns = principal axes
eigenvalues = pca.explained_variance_
ratios = pca.explained_variance_ratio_
print(ratios, " -> ", np.round(ratios * 100, 1), "%")
# Ensure right-handed orientation
if self.dim == 3 and np.linalg.det(eigenvectors) < 0:
eigenvectors[:, 2] *= -1
# Renormalize
u, _, vT = np.linalg.svd(eigenvectors)
eigenvectors = np.dot(u, vT)
print(f"Principal axes (columns):\n{eigenvectors}")
return eigenvectors, eigenvalues
[docs]
def plot_centroids(self, plot_box=True, save_path="centroid_plot.png"):
"""Plot centroid data in three orthogonal views (XY, YZ, XZ)."""
# pylint: disable=import-outside-toplevel # matplotlib is a heavy optional dep
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
if self.data is None or self.data.empty:
raise RuntimeError("No centroid data loaded. Call read_input() first.")
df = self.data.copy()
bb = self.bounding_box
if self.dim == 2:
xmin, xmax, ymin, ymax = map(float, bb)
fig, ax = plt.subplots(figsize=(5, 5))
ax.scatter(df["X"], df["Y"], s=10, color="black", alpha=0.6)
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_aspect("equal")
if plot_box:
ax.add_patch(
Rectangle(
(bb[0], bb[2]),
bb[1] - bb[0],
bb[3] - bb[2],
fill=False,
lw=1,
ls="-",
ec="r",
)
)
fig.tight_layout()
fig.savefig(save_path, dpi=300)
plt.close(fig)
print(f"2D centroid plot saved: {save_path}")
return
# 3D case: three projections
xmin, xmax, ymin, ymax, zmin, zmax = map(float, bb)
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
# Z view (X-Y)
if plot_box:
axes[0].add_patch(
Rectangle(
(xmin, ymin),
xmax - xmin,
ymax - ymin,
fill=False,
lw=1,
ls="-",
ec="r",
)
)
axes[0].scatter(df["X"], df["Y"], s=10, color="black", alpha=0.6)
axes[0].set_xlabel("X")
axes[0].set_ylabel("Y")
axes[0].set_title("Z-view (X-Y)")
# X view (Y-Z)
if plot_box:
axes[1].add_patch(
Rectangle(
(ymin, zmin),
ymax - ymin,
zmax - zmin,
fill=False,
lw=1,
ls="-",
ec="r",
)
)
axes[1].scatter(df["Y"], df["Z"], s=10, color="black", alpha=0.6)
axes[1].set_xlabel("Y")
axes[1].set_ylabel("Z")
axes[1].set_title("X-view (Y-Z)")
# Y view (X-Z)
if plot_box:
axes[2].add_patch(
Rectangle(
(xmin, zmin),
xmax - xmin,
zmax - zmin,
fill=False,
lw=1,
ls="-",
ec="r",
)
)
axes[2].scatter(df["X"], df["Z"], s=10, color="black", alpha=0.6)
axes[2].set_xlabel("X")
axes[2].set_ylabel("Z")
axes[2].set_title("Y-view (X-Z)")
for ax in axes:
ax.set_aspect("equal")
fig.tight_layout()
fig.savefig(save_path, dpi=300)
plt.close(fig)
[docs]
def remove_outside_points(self, bbox=None, tolerance=0.0):
"""Remove out-of-bounds points and their associated per-grain data from self.data."""
if self.data is None or self.data.empty:
raise RuntimeError("No data loaded. Call read_input() first.")
bbox = np.array(bbox if bbox is not None else self.bounding_box, dtype=float)
tol = tolerance / 100.0
# shrink inward
for i in range(0, len(bbox), 2):
low, high = bbox[i], bbox[i + 1]
size = high - low
bbox[i] = low + tol * size
bbox[i + 1] = high - tol * size
df = self.data
if self.dim == 3:
print(
f"\nBounding box: "
f"X=({bbox[0]:.4f},{bbox[1]:.4f}), "
f"Y=({bbox[2]:.4f},{bbox[3]:.4f}), "
f"Z=({bbox[4]:.4f},{bbox[5]:.4f})"
)
print(
f" Data extents before: "
f"X=({df['X'].min():.4f},{df['X'].max():.4f}), "
f"Y=({df['Y'].min():.4f},{df['Y'].max():.4f}), "
f"Z=({df['Z'].min():.4f},{df['Z'].max():.4f})"
)
else:
print(
f"\nBounding box: "
f"X=({bbox[0]:.4f},{bbox[1]:.4f}), "
f"Y=({bbox[2]:.4f},{bbox[3]:.4f})"
)
print(
f"Data extents before: "
f"X=({df['X'].min():.4f},{df['X'].max():.4f}), "
f"Y=({df['Y'].min():.4f},{df['Y'].max():.4f})"
)
# mask for in-box points
if self.dim == 2:
mask = (
(df["X"] >= bbox[0])
& (df["X"] <= bbox[1])
& (df["Y"] >= bbox[2])
& (df["Y"] <= bbox[3])
)
else:
mask = (
(df["X"] >= bbox[0])
& (df["X"] <= bbox[1])
& (df["Y"] >= bbox[2])
& (df["Y"] <= bbox[3])
& (df["Z"] >= bbox[4])
& (df["Z"] <= bbox[5])
)
before = len(df)
self.data = df.loc[mask].reset_index(drop=True)
removed = before - len(self.data)
print(
f"\n\nRemoved {removed} out-of-bounds points; "
f"{len(self.data)} remain.\n\n"
)
dfa = self.data
if self.dim == 3:
print(
f"Data extents after: "
f"X=({dfa['X'].min():.4f},{dfa['X'].max():.4f}), "
f"Y=({dfa['Y'].min():.4f},{dfa['Y'].max():.4f}), "
f"Z=({dfa['Z'].min():.4f},{dfa['Z'].max():.4f})"
)
else:
print(
f"Data extents after: "
f"X=({dfa['X'].min():.4f},{dfa['X'].max():.4f}), "
f"Y=({dfa['Y'].min():.4f},{dfa['Y'].max():.4f})"
)
[docs]
def validate_bounding_box(self):
"""Validate or correct the bounding box (remove_points or extend_bounding_box)."""
print("Validating bounding box...")
if self.data is None:
raise RuntimeError(
"validate_bounding_box() called before reading input data."
)
if not self.auto_fix_bbox:
print("Auto-fix disabled; using provided bounding box as-is.")
else:
df = self.data.copy()
if self.bbox_fix_mode == "remove_points":
print(
f"Auto-fix enabled: removing points outside bounding box "
f"with {self.bbox_tolerance}% tolerance."
)
self.remove_outside_points(self.bounding_box, self.bbox_tolerance)
elif self.bbox_fix_mode == "extend_bounding_box":
print(
f"Auto-fix enabled: extending bounding box to enclose all points "
f"with +{self.bbox_tolerance}% margin."
)
bbox = np.array(self.bounding_box, dtype=float)
tol = self.bbox_tolerance / 100.0
if self.dim == 2:
data_extents = [
df["X"].min(),
df["X"].max(),
df["Y"].min(),
df["Y"].max(),
]
else:
data_extents = [
df["X"].min(),
df["X"].max(),
df["Y"].min(),
df["Y"].max(),
df["Z"].min(),
df["Z"].max(),
]
# Expand bounding box to include data + margin
new_bbox = bbox.copy()
for i in range(0, len(bbox), 2):
low, high = bbox[i], bbox[i + 1]
data_min, data_max = data_extents[i], data_extents[i + 1]
low_new = min(low, data_min)
high_new = max(high, data_max)
size = high_new - low_new
new_bbox[i] = low_new - tol * size
new_bbox[i + 1] = high_new + tol * size
print(f"-> Original bounding box: {bbox.tolist()}")
print(f"-> Data extents: {data_extents}")
print(f"-> Extended bounding box: {new_bbox.tolist()}")
self.bounding_box = new_bbox.tolist()
else:
raise ValueError(
f"Unexpected bbox_fix_mode='{self.bbox_fix_mode}'. "
"Valid options are: ['remove_points', 'extend_bounding_box']."
)
df = self.data
bx = self.bounding_box
if self.dim == 2:
min_x, max_x = df["X"].min(), df["X"].max()
min_y, max_y = df["Y"].min(), df["Y"].max()
if bx[0] > min_x or bx[1] < max_x or bx[2] > min_y or bx[3] < max_y:
raise ValueError(
f"Bounding box invalid for 2D data.\n"
f"Data range X=({min_x}, {max_x}), Y=({min_y}, {max_y})\n"
f"Bounding box X=({bx[0]}, {bx[1]}), Y=({bx[2]}, {bx[3]})"
)
else:
min_x, max_x = df["X"].min(), df["X"].max()
min_y, max_y = df["Y"].min(), df["Y"].max()
min_z, max_z = df["Z"].min(), df["Z"].max()
if (
bx[0] > min_x
or bx[1] < max_x
or bx[2] > min_y
or bx[3] < max_y
or bx[4] > min_z
or bx[5] < max_z
):
raise ValueError(
f"Bounding box invalid for 3D data.\n"
f"Data range X=({min_x}, {max_x}), Y=({min_y}, {max_y}), Z=({min_z}, {max_z})\n"
f"Bounding box X=({bx[0]}, {bx[1]}), Y=({bx[2]}, {bx[3]}), Z=({bx[4]}, {bx[5]})"
)
[docs]
def build_voronoi(
self,
option: str = "voronoi",
generate_mesh: bool = False,
relative_el_size: float = None,
morphoalgo: str = "praxis",
mesh_quality_min: float = 0.9,
tesr_size: list = None,
CVT_iter: int = 1000,
):
"""Build a Voronoi (or Laguerre) tessellation using Neper.
Reads the input CSV (if not already loaded), validates/fixes the bounding box, runs NEPER to produce the tessellation, applies the rotation to the cell
orientations and strains, evaluates the output, optionally meshes it, and writes the reformatted voxel reconstruction CSV/VTK plus the CPFE strain files.
Args:
option: Tessellation morphology. "voronoi" seeds a Voronoi tessellation at the input coordinates (and applies the GrainRadius weights when ``weighted``); "centroidal" seeds a CVT relaxation at them.
"centroid" instead treats the coordinates as an optimization *target* -- NEPER moves seeds until the cell centroids match them,
starting from its own packing rather than from the data -- and "centroidsize" targets centroid and size together (requires GrainRadius).
Note "centroid" and "centroidal" are distinct NEPER morphologies, not aliases.
generate_mesh: If True, also generate a NEPER/GMSH tet (or tri) .msh mesh; this is a fallback only, the default CPFE mesh is SCULPT hex.
relative_el_size: Mesh element size relative to grain size, passed to generate_mesh; defaults to 1.0 when None.
morphoalgo: NEPER morphology optimization algorithm, e.g. "praxis", "subplex", or "lloyd".
mesh_quality_min: Minimum acceptable mesh element quality, passed to generate_mesh.
tesr_size: Voxel raster resolution [nx, ny, nz] for the .tesr output; defaults to [20, 20, 20] when None.
CVT_iter: Number of CVT/morphology optimization iterations.
"""
if tesr_size is None:
tesr_size = [20, 20, 20]
valid_options = ["voronoi", "centroid", "centroidal", "centroidsize"]
if option not in valid_options:
raise ValueError(
f"Invalid option='{option}'. Must be one of {valid_options}."
)
if self.data is None or self.data.empty:
self.read_input()
self.validate_bounding_box()
df = self.data.copy()
os.makedirs(self.output_dir, exist_ok=True)
input_path = os.path.join(self.output_dir, "points.dat")
weight_path = os.path.join(self.output_dir, "weights.dat")
orientation_path = os.path.join(self.output_dir, "orientations.dat")
if self.dim == 2:
coord_cols = ["X", "Y"]
else:
coord_cols = ["X", "Y", "Z"]
# File creation depending on option
if option == "centroidsize":
if "GrainRadius" not in df.columns:
raise ValueError("centroidsize mode requires 'GrainRadius' column.")
# combine coordinates + weight into one file
df_out = df[coord_cols + ["GrainRadius"]]
df_out.to_csv(input_path, sep=" ", index=False, header=False)
print(f"Neper centroidsize input file created: {input_path}")
else:
df[coord_cols].to_csv(input_path, sep=" ", index=False, header=False)
print(f"Neper coordinate file created: {input_path}")
if self.weighted:
if "GrainRadius" not in df.columns:
raise ValueError("Weighted mode requires 'GrainRadius' column.")
if self.dim == 3:
df["Weight"] = (4.0 / 3.0) * np.pi * (df["GrainRadius"] ** 3)
weight_type = "volume"
else:
df["Weight"] = np.pi * (df["GrainRadius"] ** 2)
weight_type = "area"
df["Weight"] /= df["Weight"].sum()
df["Weight"].to_csv(weight_path, sep=" ", index=False, header=False)
print(f"Neper weight file created ({weight_type}-based): {weight_path}")
# Orientation (3-column)
if all(c in df.columns for c in self.angle_id):
df[self.angle_id].to_csv(
orientation_path, sep=" ", index=False, header=False
)
print(f"Neper orientation file created: {orientation_path}")
ori_args = ["-ori", f"file({orientation_path},des={self.ori_descriptor})"]
else:
ori_args = []
# Compute domain
if self.dim == 2:
xmin, xmax, ymin, ymax = self.bounding_box
sx, sy = xmax - xmin, ymax - ymin
tx, ty = xmin, ymin
domain_arg = f"square({sx},{sy}):translate({tx},{ty})"
else:
xmin, xmax, ymin, ymax, zmin, zmax = self.bounding_box
sx, sy, sz = xmax - xmin, ymax - ymin, zmax - zmin
tx, ty, tz = xmin, ymin, zmin
domain_arg = f"cube({sx},{sy},{sz}):translate({tx},{ty},{tz})"
# Morphological arguments
if option == "centroidsize":
morpho_args = [
"-morpho",
f"centroidsize:file({input_path})",
"-morphooptistop",
f"iter={CVT_iter}",
"-morphooptialgo",
morphoalgo,
]
print(
f"\nCentroid-size tessellation ({self.dim}D): combined coordinate+size input."
)
elif option == "centroid":
morpho_args = [
"-morpho",
f"centroid:file({input_path})",
"-morphooptiini",
"coo:LLLFP2011",
"-morphooptistop",
f"iter={CVT_iter}",
"-morphooptialgo",
morphoalgo,
]
print(
f"\nCentroid tessellation ({self.dim}D): only coordinate input. Weighted option will be ignored."
)
elif self.weighted:
morphoalgo = "lloyd"
morpho_args = [
"-morpho",
f"{option}",
"-morphooptiini",
f"coo:file({input_path}),weight:file({weight_path})",
"-morphooptistop",
f"iter={CVT_iter}",
"-morphooptialgo",
morphoalgo,
]
print(
f"\nWeighted tessellation ({self.dim}D): using weight file {weight_path}."
)
else:
morpho_args = [
"-morpho",
f"{option}",
"-morphooptiini",
f"coo:file({input_path})",
"-morphooptistop",
f"iter={CVT_iter}",
"-morphooptialgo",
morphoalgo,
]
print(f"\nUnweighted tessellation ({self.dim}D): standard Poisson-Voronoi.")
tess_name = os.path.join(self.output_dir, "reconstruction")
neper_cmd = (
[
self.neper_bin,
"-T",
"-n",
str(len(df)),
"-reg",
str(1),
"-dim",
str(self.dim),
"-domain",
domain_arg,
"-oridescriptor",
"rotmat",
"-o",
tess_name,
"-format",
"tess,geo,tesr",
"-tesrsize",
f"{tesr_size[0]},{tesr_size[1]},{tesr_size[2]}",
"-tesrformat",
"ascii",
]
+ morpho_args
+ ori_args
)
print("\n=== Running Neper Tessellation ===\n")
print("> " + " ".join(neper_cmd))
log_path = os.path.join(self.output_dir, "neper_voronoi_builder.log")
with open(log_path, "w", encoding="utf-8") as logf:
print("\n=== Running Neper Tessellation ===", file=logf)
print("> " + " ".join(neper_cmd), file=logf)
logf.flush()
subprocess.run(
neper_cmd,
check=True,
env=self.env,
stderr=subprocess.STDOUT,
)
print(f"Voronoi tessellation completed: {tess_name}.tess\n")
print("Evaluating tessellation correctness, in figures folder\n")
print("\n=== Updating cell properties with rotation ===\n")
self.apply_rotation_to_properties(
tess_file=tess_name + ".tess",
)
# apply_rotation_to_properties rotates coords + the .tess/.ori matrices, but
# orientations.dat (written above) and the reformatted Euler are still in the
# unrotated frame. When a sample rotation is applied, re-derive them so the
# orientations.dat consumers (CPFE via euler_to_mrp, VoxelMeshBuilder) share
# the rotated mesh frame. No-op when no rotation is applied.
if self.angle_id is not None and not np.allclose(self.rotate_matrix, np.eye(3)):
self._rewrite_orientations_dat_rotated(orientation_path)
self.evaluate_output_voronoi(tess_name + ".tess")
if generate_mesh:
self.generate_mesh(
tess_file=tess_name + ".tess",
output_name=tess_name,
format_type=["msh", "vtk"],
relative_cl=relative_el_size if relative_el_size else 1.0,
mesh_quality_min=mesh_quality_min,
)
print("\n=== Reformatting .tesr file ===\n")
self.reformat_tesr_file(
tesr_file=tess_name + ".tesr", orientation_file=orientation_path
)
[docs]
def build_graph(
self,
device: str = "cpu",
option: str = "centroid",
CVT_iter: int = 100,
morphoalgo: str = "praxis",
visualize2D: bool = False,
visualize3D: bool = False,
):
"""Build a Voronoi tessellation and convert it to a graph representation.
Runs build_voronoi (no mesh) then parses the resulting .tess into a per-cell graph via NeperTessToGraphNN, registering the loaded dataframe features on it.
Args:
device: Torch device for the graph tensors, e.g. "cpu" or "cuda".
option: Tessellation seeding mode passed to build_voronoi: "voronoi", "centroid", or "centroidsize".
CVT_iter: Number of CVT/morphology optimization iterations passed to build_voronoi.
morphoalgo: NEPER morphology optimization algorithm passed to build_voronoi.
visualize2D: If True, save a 2D rendering of the graph under output_dir/figures/gnn.
visualize3D: If True, save a 3D rendering of the graph under output_dir/figures/gnn.
Returns:
The per-cell graph object produced by NeperTessToGraphNN.build_cell_graph.
"""
self.build_voronoi(
generate_mesh=False,
option=option,
relative_el_size=None,
morphoalgo=morphoalgo,
CVT_iter=CVT_iter,
)
# pylint: disable=import-outside-toplevel # torch/tess_to_gnn are heavy optional deps
from .tess_to_gnn import NeperTessToGraphNN
import torch
print("\n=== Building Graph Neural Network representation ===\n")
parser = NeperTessToGraphNN(
tess_path=os.path.join(self.output_dir, "reconstruction.tess"),
device=device,
dtype=torch.float64,
)
parser.register_dataframe_features(data=self.data, verbose=False)
print(self.data.columns)
graph = parser.build_cell_graph()
if visualize3D:
os.makedirs(self.output_dir + "/figures/gnn", exist_ok=True)
parser.visualize_graph_3D(
graph, outpath=self.output_dir + "/figures/gnn/graph_3D.png"
)
if visualize2D:
os.makedirs(self.output_dir + "/figures/gnn", exist_ok=True)
parser.visualize_graph_2D(
graph, outpath=self.output_dir + "/figures/gnn/graph_2D.png"
)
return graph
def _rewrite_orientations_dat_rotated(self, orientation_path: str):
"""Rewrite orientations.dat (Euler, degrees) in the rotated sample frame.
Mirrors the rotation ``apply_rotation_to_properties`` applies to the matrices,
but stays entirely in the ``orientation_helper`` convention that the
orientations.dat consumers use, so no NEPER matrix-convention dependency
enters. ``self.data`` Euler are degrees after ``read_input``.
"""
# pylint: disable=import-outside-toplevel # heavy/optional deps
import torch
from .orientation_helper import euler_to_matrix, matrix_to_euler
convention = self.ori_descriptor.replace("euler-", "")
euler = self.data[self.angle_id].to_numpy(dtype=float)
g0 = euler_to_matrix(torch.tensor(euler), convention, "degrees")
Rr = torch.tensor(self.rotate_matrix, dtype=g0.dtype)
if self.orientation_active_convention:
g_rot = torch.einsum("ij,njk->nik", Rr, g0) # Rr @ g
else:
g_rot = torch.einsum("nij,kj->nik", g0, Rr) # g @ Rr^T
new_euler = matrix_to_euler(g_rot, convention, "degrees").numpy()
self.data[self.angle_id] = new_euler
np.savetxt(orientation_path, new_euler, fmt="%.10g")
print(f"Rewrote {orientation_path} in the rotated sample frame.")
[docs]
def apply_rotation_to_properties(self, tess_file: str):
"""Apply the rotation matrix to cell orientations (and strains) in the .tess file."""
if not hasattr(self, "data") or self.data is None or self.data.empty:
raise RuntimeError(
"No data loaded. Run read_input() first before applying rotation to properties."
)
if not os.path.exists(tess_file):
raise FileNotFoundError(f"Tessellation file not found: {tess_file}")
transpose = self.orientation_active_convention
# Backup file
backup_path = tess_file + ".bak"
if not os.path.exists(backup_path):
shutil.copy(tess_file, backup_path)
# Read file and locate *ori section
with open(tess_file, "r", encoding="utf-8") as f:
lines = f.readlines()
ori_start, ori_end = None, None
for i, line in enumerate(lines):
if line.strip().lower().startswith("*ori"):
ori_start = i
continue
if (
ori_start is not None
and line.strip().startswith("*")
and i > ori_start + 1
):
ori_end = i
break
if ori_start is None:
raise RuntimeError(f"No *ori section found in {tess_file}.")
if ori_end is None:
ori_end = len(lines)
# skip descriptor line
ori_lines = lines[ori_start + 2 : ori_end]
num_cells = len(ori_lines)
if num_cells == 0:
raise RuntimeError("No orientation data lines found after *ori header.")
# No orientation provided -> identity orientation for all cells
if self.angle_id is None:
print("No orientation data provided; assigning identity matrices.")
identity_rot = np.eye(3)
ori_rot = np.tile(identity_rot, (num_cells, 1, 1))
else:
ori_data = np.array(
[list(map(float, ln.split())) for ln in ori_lines], dtype=float
)
if ori_data.shape[1] != 9:
raise ValueError(
f"Each *ori line must have 9 numbers; found {ori_data.shape[1]}."
"this indicates corrupted .tess file."
)
ori_mats = ori_data.reshape((-1, 3, 3))
if transpose:
ori_rot = np.array([np.dot(self.rotate_matrix, R) for R in ori_mats])
else:
ori_rot = np.array(
[np.dot(self.rotate_matrix, R.T).T for R in ori_mats]
)
self.data[self.ori_rotmat_id] = ori_rot.reshape(-1, 9)
# Write back new orientations
new_ori_lines = [
" ".join(f"{val:.8f}" for val in R.flatten()) + "\n" for R in ori_rot
]
lines[ori_start + 2 : ori_end] = new_ori_lines
with open(tess_file, "w", encoding="utf-8") as f:
f.writelines(lines)
print(f"\nUpdated *ori section in {tess_file}")
print(f"by applying the rotation matrix:\n{self.rotate_matrix}")
# Apply rotation for elastic strain if available
if self.elastic_strain_id is not None:
print("\nApply rotation to elastic strain tensors\n")
strain_tensors = (
self.data[self.elastic_strain_id]
.to_numpy(dtype=float)
.reshape(-1, 3, 3)
)
# Rotate: e' = R e R^T
R = self.rotate_matrix
strain_rot = np.array([np.dot(R, np.dot(E, R.T)) for E in strain_tensors])
self.data[self.elastic_strain_id] = strain_rot.reshape(-1, 9)
else:
print("No elastic strain identifier provided; skipping strain rotation.")
[docs]
def generate_mesh(
self,
tess_file: str,
output_name: str = "voronoi",
format_type=("msh4",),
relative_cl: float = 1.0,
mesh_quality_min: float = 0.9,
interface_type: str = "continuous",
partition: int = 16,
):
"""Generate an FE mesh from a Neper .tess (tri for 2D, tet for 3D, 2nd order)."""
print("\n=== Generating Mesh from Tessellation ===")
# Normalize formats
if isinstance(format_type, str):
fmt_arg = format_type
elif isinstance(format_type, (list, tuple)):
fmt_arg = ",".join(str(f).strip() for f in format_type)
else:
raise TypeError(
f"format_type must be str/list/tuple, got {type(format_type)}"
)
if not os.path.exists(tess_file):
raise FileNotFoundError(f"Tessellation file not found: {tess_file}")
element_type = "tri" if self.dim == 2 else "tet"
neper_cmd = [
self.neper_bin,
"-M",
tess_file,
"-elttype",
element_type,
"-rcl",
str(relative_cl),
"-order",
"2", # always 2nd order
"-meshqualmin",
str(mesh_quality_min),
"-interface",
interface_type,
"-format",
fmt_arg,
"-o",
output_name,
"-part",
str(partition),
]
print(
"\n=== Running Neper Meshing ===\n"
"(this could take a while for large system)...\n"
)
print("> " + " ".join(neper_cmd))
print("\n")
log_path = os.path.join(self.output_dir, "neper_voronoi_mesh.log")
with open(log_path, "w", encoding="utf-8") as logf:
print("> " + " ".join(neper_cmd), file=logf)
logf.flush()
subprocess.run(
neper_cmd,
check=True,
env=self.env,
stdout=logf,
stderr=subprocess.STDOUT,
)
for fmt in fmt_arg.split(","):
print(f"Generated mesh file: {output_name}.{fmt}")
print("\n")
[docs]
def export_cell_properties(
self,
tess_file: str,
output_cell="id,vol,w,x,y,z,radeq",
output_seed="id,w,x,y,z",
):
"""Compute Neper cell/seed stats, merge them, and export CPFE strain files."""
print("\n=== Computing Voronoi Cell and Seed Properties ===")
if not os.path.exists(tess_file):
raise FileNotFoundError(f"Tessellation file not found: {tess_file}")
# Run Neper to generate .statcell and .statseed
neper_cmd = [
self.neper_bin,
"-T",
"-loadtess",
tess_file,
"-statcell",
output_cell,
"-statseed",
output_seed,
"-format",
"ori",
"-oridescriptor",
"rotmat",
]
print("> " + " ".join(neper_cmd))
subprocess.run(neper_cmd, check=True, env=self.env)
base = Path(tess_file).with_suffix("")
stat_cell_file = base.with_suffix(".stcell")
stat_seed_file = base.with_suffix(".stseed")
out_file = Path(str(base) + "_out.csv")
if not stat_cell_file.exists():
raise FileNotFoundError(f"Missing .stcell file: {stat_cell_file}")
if not stat_seed_file.exists():
raise FileNotFoundError(f"Missing .stseed file: {stat_seed_file}")
cols_cell = output_cell.split(",")
cols_seed = output_seed.split(",")
df_cell = pd.read_csv(stat_cell_file, sep=r"\s+", header=None, names=cols_cell)
df_seed = pd.read_csv(stat_seed_file, sep=r"\s+", header=None, names=cols_seed)
df_cell = df_cell.add_suffix("_cell")
df_seed = df_seed.add_suffix("_seed")
df_merged = pd.merge(
df_cell, df_seed, left_on="id_cell", right_on="id_seed", how="outer"
)
df_merged.to_csv(out_file, sep=",", index=False)
print(f"Exported merged cell/seed data to {out_file}\n")
# Export elastic strain (.ee)
num_cells = len(df_cell)
ee_file = Path(self.output_dir) / (base.name + ".ee")
if (
self.elastic_strain_id is not None
and hasattr(self, "data")
and not self.data.empty
):
strain_data = self.data[self.elastic_strain_id].to_numpy(dtype=float)
else:
strain_data = np.zeros((num_cells, 9), dtype=float)
# Write .ee file
with open(ee_file, "w", encoding="utf-8") as f:
for row in strain_data:
f.write(" ".join(f"{v:.8e}" for v in row) + "\n")
print(f"Exported elastic strain tensors to {ee_file}\n")
# Write .csv file for MOOSE
df_strain = pd.read_csv(ee_file, sep=r"\s+", header=None)
moose_file = Path(self.output_dir) / (base.name + "_cpfe_ee.csv")
moose_dat = df_seed[["x_seed", "y_seed", "z_seed"]].copy()
moose_dat = pd.merge(moose_dat, df_strain, left_index=True, right_index=True)
moose_dat.to_csv(moose_file, sep=",", index=False, header=False)
[docs]
def evaluate_output_voronoi(self, tess_file: str, length_norm: bool = False):
"""Evaluate Neper Voronoi output: compare cell vs seed props, save figures."""
# pylint: disable=import-outside-toplevel # matplotlib is a heavy optional dep
import matplotlib.pyplot as plt
if not os.path.exists(tess_file):
raise FileNotFoundError(f"Tessellation file not found: {tess_file}")
self.export_cell_properties(tess_file)
outputdir = self.output_dir
os.makedirs(os.path.join(outputdir, "figures"), exist_ok=True)
df_path = Path(str(Path(tess_file).with_suffix("")) + "_out.csv")
df = pd.read_csv(df_path)
self.plot_centroids(
save_path=os.path.join(outputdir, "figures", "centroids.png")
)
# Consistency check between IDs
id_cell = np.array(df["id_cell"])
id_seed = np.array(df["id_seed"])
if not np.array_equal(id_cell, id_seed):
print("Warning: id_cell and id_seed are not identical!")
dfdat = self.data.copy()
if self.dim == 2:
coord_cols = ["X", "Y"]
elif self.dim == 3:
coord_cols = ["X", "Y", "Z"]
else:
raise ValueError("Unsupported dimension (must be 2 or 3).")
_ = dfdat[coord_cols]
seed_w = np.array(df["w_seed"]) # input weight
cell_xyz = np.array(df[["x_cell", "y_cell", "z_cell"]])
input_xyz = np.array(dfdat[coord_cols])
cell_radeq = np.array(df["radeq_cell"])
seed_radeq = self.data["GrainRadius"].to_numpy()
bb = self.bounding_box
if self.dim == 3:
xmin, xmax, ymin, ymax, zmin, zmax = bb
Lbox = np.array([xmax - xmin, ymax - ymin, zmax - zmin], dtype=float)
else:
xmin, xmax, ymin, ymax = bb
Lbox = np.array([xmax - xmin, ymax - ymin], dtype=float)
eps = 1e-12
Lbox = 1.0
if length_norm:
Lbox = np.maximum(Lbox, eps)
# Position error (normalized by box lengths per-axis)
norm_cell = cell_xyz / Lbox
norm_seed = input_xyz / Lbox
delta = norm_cell - norm_seed
error_xyz = np.linalg.norm(delta, axis=1)
# Equivalent radius error (%)
error_radeq = (cell_radeq - seed_radeq) / seed_radeq * 100.0
# Error histograms
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.hist(error_xyz, bins=20)
if length_norm:
plt.xlabel("|(Xi_cell - Xi_input)/box_Xi| (norm-2)")
else:
plt.xlabel("|Xi_cell - Xi_input| (norm-2) (input length units)")
plt.ylabel("Frequency")
plt.subplot(1, 2, 2)
plt.hist(error_radeq, bins=20)
plt.xlabel("Relative Equivalent-Radius Error (%)")
plt.ylabel("Frequency")
plt.suptitle("Voronoi Cell vs Seed Point Errors")
plt.tight_layout()
plt.savefig(os.path.join(outputdir, "figures", "voronoi_error.png"), dpi=300)
plt.close()
# Distribution comparison plots
plt.figure(figsize=(12, 5))
plt.subplot(1, 3, 1)
plt.hist(seed_w, bins=20)
plt.xlabel("Seed Weight (from Grain Equivalent Volume)")
plt.ylabel("Frequency")
plt.subplot(1, 3, 2)
plt.hist(seed_radeq, bins=20)
plt.xlabel("Seed Equivalent Radius")
plt.ylabel("Frequency")
plt.subplot(1, 3, 3)
plt.hist(cell_radeq, bins=20)
plt.xlabel("Cell Equivalent Radius")
plt.ylabel("Frequency")
plt.tight_layout()
plt.savefig(
os.path.join(outputdir, "figures", "distribution_comparison.png"), dpi=300
)
plt.close()
print(f"Evaluation figures saved in {os.path.join(outputdir, 'figures')}")