import os
import re
import logging
import numpy as np
from collections import defaultdict
from .units import length_scale_from_um
from .entities import Ellipsoid, Cuboid
from .initializations import NodeSets, validate_matrix_mapping
from typing import Dict, Any, Optional, Mapping, Union
[docs]
def write_dump(Ellipsoids: list[Any], sim_box: Any) -> None:
"""
Write .dump files into a sub-directory "dump_files" for visualization or reuse in Kanapy
Parameters
----------
Ellipsoids : list
Contains information of ellipsoids such as positions, axes lengths, and tilt angles
sim_box : Cuboid
Contains information of the dimensions of the simulation box
Returns
-------
None
Notes
-----
The function creates a 'dump_files' directory in the current working directory if it does not exist
and writes a file named 'particle.<timestep>.dump', where <timestep> is given by sim_box.sim_ts.
The generated .dump files can be read by visualization software like OVITO or imported back into Kanapy
to avoid repeating the packing simulation. Each file contains simulation domain bounds and ellipsoid
attributes including position, axes lengths, and quaternion orientation.
"""
num_particles = len(Ellipsoids)
cwd = os.getcwd()
output_dir = os.path.join(cwd, 'dump_files')
dump_outfile = os.path.join(output_dir, 'particle') # output dump file
if not os.path.exists(output_dir):
os.makedirs(output_dir)
with open(dump_outfile + ".{0}.dump".format(sim_box.sim_ts), 'w') as f:
f.write('ITEM: TIMESTEP\n')
f.write('{0}\n'.format(sim_box.sim_ts))
f.write('ITEM: NUMBER OF ATOMS\n')
f.write('{0}\n'.format(num_particles))
f.write('ITEM: BOX BOUNDS ff ff ff\n')
f.write('{0} {1}\n'.format(sim_box.left, sim_box.right))
f.write('{0} {1}\n'.format(sim_box.bottom, sim_box.top))
f.write('{0} {1}\n'.format(sim_box.back, sim_box.front))
f.write(
'ITEM: ATOMS id x y z AsphericalShapeX AsphericalShapeY AsphericalShapeZ OrientationX OrientationY ' +
'OrientationZ OrientationW\n')
for ell in Ellipsoids:
qw, qx, qy, qz = ell.quat
f.write('{0} {1} {2} {3} {4} {5} {6} {7} {8} {9} {10} {11}\n'.format(
ell.id, ell.x, ell.y, ell.z, ell.a, ell.b, ell.c, qx, qy, qz, qw, ell.phasenum))
[docs]
def read_dump(file: Union[str, os.PathLike[str]]) -> tuple[Any, list[Any]]:
"""
Read a .dump file to extract simulation box and ellipsoid information for voxelization
Parameters
----------
file : str
Path to the .dump file containing ellipsoid data generated in the packing routine
Returns
-------
sim_box : Cuboid
Cuboid object representing the RVE (simulation box)
Ellipsoids : list of Ellipsoid
List of ellipsoid objects representing the grains
Raises
------
FileNotFoundError
If the specified .dump file does not exist
Notes
-----
The function reads the simulation box dimensions and ellipsoid attributes including
position, axes lengths, orientation, and phase number. Duplicate ellipsoids are identified
based on naming convention in the dump file
"""
print(' Reading the .dump file for particle information')
try:
# Read the Simulation box dimensions
with open(file, 'r+') as fd:
lookup = "ITEM: NUMBER OF ATOMS"
lookup2 = "ITEM: BOX BOUNDS ff ff ff"
for num, lines in enumerate(fd, 1):
if lookup in lines:
number_particles = int(next(fd))
par_line_num = num + 7
if lookup2 in lines:
valuesX = re.findall(r'\S+', next(fd))
RVE_minX, RVE_maxX = list(map(float, valuesX))
valuesY = re.findall(r'\S+', next(fd))
RVE_minY, RVE_maxY = list(map(float, valuesY))
valuesZ = re.findall(r'\S+', next(fd))
RVE_minZ, RVE_maxZ = list(map(float, valuesZ))
except FileNotFoundError:
raise FileNotFoundError('.dump file not found, execute packing with option save_files=True')
# Create an instance of simulation box
sim_box = Cuboid(RVE_minX, RVE_maxY, RVE_maxX, RVE_minY, RVE_maxZ, RVE_minZ)
# Read the particle shape & position information
# Create instances for ellipsoids & assign values from dump files
Ellipsoids = []
with open(file, "r") as f:
count = 0
for num, lines in enumerate(f, 1):
if num >= par_line_num:
count += 1
values: list = re.findall(r'\S+', lines)
int_values = list(map(float, values[1:]))
values = [values[0]] + int_values
iden = count # ellipsoid 'id'
a, b, c = values[4], values[5], values[6] # Semi-major length, Semi-minor length 1 & 2
x, y, z = values[1], values[2], values[3]
qx, qy, qz, qw = values[7], values[8], values[9], values[10]
ip = int(values[11])
quat = np.array([qw, qx, qy, qz])
ellipsoid = Ellipsoid(iden, x, y, z, a, b, c, quat, phasenum=ip) # instance of Ellipsoid class
# Find the original particle if the current is duplicate
for c in values[0]:
if c == '_':
split_str = values[0].split("_")
original_id = int(split_str[0])
ellipsoid.duplicate = original_id
break
else:
continue
Ellipsoids.append(ellipsoid)
return sim_box, Ellipsoids
def _abaqus_phase_options(props_file, crystal_plasticity, nphases):
"""Validate explicit phase options; None retains legacy material selection."""
if props_file is None and not isinstance(crystal_plasticity, list):
return None, None
if nphases > 1:
if props_file is not None and not isinstance(props_file, list):
raise ValueError('props_file must be a list with one entry per phase.')
if not isinstance(crystal_plasticity, list):
raise ValueError('crystal_plasticity must be a list with one entry per phase.')
files = props_file if isinstance(props_file, list) else [props_file] * nphases
cp = crystal_plasticity if isinstance(crystal_plasticity, list) else [
props_file is not None if crystal_plasticity is None else crystal_plasticity] * nphases
if len(files) != nphases or len(cp) != nphases:
raise ValueError('props_file and crystal_plasticity must contain one entry per phase.')
if any(not isinstance(value, (bool, np.bool_)) for value in cp):
raise ValueError('crystal_plasticity entries must be booleans.')
if props_file is not None and any(flag and inc is None for flag, inc in zip(cp, files)):
raise ValueError('Each crystal_plasticity phase requires a props_file; None requires False.')
return files, cp
[docs]
def export2abaqus(
nodes: Any,
file: Union[str, os.PathLike[str]],
grain_dict: Mapping[Any, Any],
voxel_dict: Mapping[Any, Any],
units: str = 'mm',
gb_area=None, dual_phase=False, thermal=False,
ialloy=None, grain_phase_dict=None,
crystal_plasticity=False, phase_props=None,
boundary_conditions: Optional[Dict[str, Any]] = None, props_file=None):
"""
Create an ABAQUS input file with microstructure morphology information including nodes, elements, and element sets
Parameters
----------
nodes : list
List of node objects or coordinates
file : str
Output file path for the ABAQUS input file
grain_dict : dict
Dictionary containing grain information
voxel_dict : dict
Dictionary containing voxel-to-grain mapping
units : str, optional
Output length units: 'µm', 'um', 'mm', or 'm' (default is 'mm').
Input node coordinates are always in micrometers.
gb_area : optional
Grain boundary area information
dual_phase : bool, optional
If True, element sets with phase numbers are defined and assigned to materials "PHASE_{phase_id}MAT"
If False, each grain refers to a material "GRAIN_{grain_id}MAT"
thermal : bool, optional
If True, thermal properties are included
ialloy : optional
Alloy information for grains, used if dual_phase is False
grain_phase_dict : dict, optional
Mapping from grain ID to phase ID
crystal_plasticity : bool or list of bool, optional
Per-phase CP selection. With multiphase props_file, a list of the same
length is required. CP uses grain materials and requires dual_phase=False.
props_file : str, path-like, list, or None, optional
One include per phase: CP constants belong in _mat.inp; standard J2
property includes are written under PHASE{id}_MAT in this geometry file.
None entries require crystal_plasticity=False and produce only section
references, leaving material definitions to the user in CAE or _geom.inp.
Grain 0 always references PHASE0_MAT, which receives an empty material
definition when its include is None so it remains visible in CAE.
phase_props : optional
Phase properties for crystal plasticity
boundary_conditions : dict, optional
Dictionary defining boundary conditions, keys may include type, loading_direction, value, etc.
Notes
-----
If dual_phase is False, a "_mat.inp" file with the same name is expected to define alloy number
and Euler angles for each grain. Element sets are written for grains or phases depending on dual_phase.
The function handles nodes, elements, sets, and optionally boundary conditions for ABAQUS simulations.
"""
from kanapy import __version__
scale_fact = length_scale_from_um(units)
phase_map = dict(grain_phase_dict or {})
validate_matrix_mapping(phase_map)
if not dual_phase and 0 in grain_dict:
if any(gid != 0 and gid not in phase_map for gid in grain_dict):
raise ValueError('Matrix exports require an explicit phase ID for every grain.')
phase_map[0] = 0
validate_matrix_mapping(phase_map)
nphases = max(phase_map.values(), default=0) + 1
if dual_phase:
nphases = max(grain_dict, default=0) + 1
if isinstance(props_file, list):
nphases = max(nphases, len(props_file))
if isinstance(crystal_plasticity, list):
nphases = max(nphases, len(crystal_plasticity))
files, cp = _abaqus_phase_options(props_file, crystal_plasticity, nphases)
if dual_phase and cp is not None and any(cp):
raise ValueError('Crystal plasticity requires grain-wise sets (dual_phase=False).')
if not dual_phase and 0 in grain_dict and cp is not None and cp[0]:
raise ValueError('Grain 0 requires standard plasticity in phase 0.')
grain_phase_dict = phase_map or None
def write_node_set(name, nset):
"""
Write a node set to a file with formatting compatible with ABAQUS input
Parameters
----------
name : str
Name of the node set to be written
nset : list of int
List of node indices belonging to the node set
Notes
-----
Nodes are written in rows of 16 entries. Node indices are incremented by 1
to match ABAQUS indexing convention (1-based)
"""
f.write(name)
for i, val in enumerate(nset[:-1], start=1):
if i % 16 == 0:
f.write(f'{val + 1}\n')
else:
f.write(f'{val + 1}, ')
f.write(f'{nset[-1] + 1}\n')
def write_grain_sets():
"""
Write element sets for grains and assign materials in an ABAQUS input file
Notes
-----
For each grain in grain_dict, an element set named "GRAIN{grain_id}_SET" is written.
Elements are written in rows of 16 entries. Each set is assigned to a material:
- If grain_phase_dict is None or the phase ID is less than nall, the grain is assigned
"GRAIN{grain_id}_MAT"
- Otherwise, the grain is assigned "PHASE{phase_id}_MAT"
The function updates ph_set with phase IDs used for dual-phase assignments
"""
for k, v in grain_dict.items():
f.write('*ELSET, ELSET=GRAIN{0}_SET\n'.format(k))
for enum, el in enumerate(v[:-1], start=1):
if enum % 16 != 0:
f.write('%d, ' % el)
else:
f.write('%d\n' % el)
f.write('%d\n' % v[-1])
for k in grain_dict.keys():
pid = phase_map.get(k, 0)
use_cp = cp[pid] if cp is not None else (grain_phase_dict is None or pid < nall)
if k != 0 and use_cp:
f.write(
'*Solid Section, elset=GRAIN{0}_SET, material=GRAIN{1}_MAT\n'
.format(k, k))
else:
f.write(
'*Solid Section, elset=GRAIN{0}_SET, material=PHASE{1}_MAT\n'
.format(k, pid))
ph_set.add(pid)
return
def write_phase_sets():
"""
Write element sets for phases and assign phase materials in an ABAQUS input file
Notes
-----
For each phase in grain_dict, an element set named "PHASE{phase_id}_SET" is written.
Elements are written in rows of 16 entries. Each set is assigned to a material
"PHASE{phase_id}_MAT". The function updates ph_set with all phase IDs used.
"""
for k, v in grain_dict.items():
if len(v) == 0:
continue
f.write('*ELSET, ELSET=PHASE{0}_SET\n'.format(k))
for enum, el in enumerate(v, start=1):
if enum % 16 != 0:
if enum == len(v):
f.write('%d\n' % el)
else:
f.write('%d, ' % el)
else:
f.write('%d\n' % el)
for k in grain_dict.keys():
if len(grain_dict[k]) == 0:
continue
f.write(
'*Solid Section, elset=PHASE{0}_SET, material=PHASE{1}_MAT\n'
.format(k, k))
ph_set.add(k)
return
def write_surface_sets():
"""
Write surface element sets for the simulation domain in an ABAQUS input file
Notes
-----
The function identifies the maximum faces along X, Y, and Z directions using the node grid.
For each face, an ELSET and corresponding SURFACE block is written. Element IDs are sorted
and written in rows of 16 entries per line. The helper function _emit handles formatting
and writing of each face. Surface sets are named "_Surf-{idx}_{label}" and the corresponding
SURFACE block is named "Surf-{idx}" with type ELEMENT.
"""
# 1) Find grid dimensions
xs = np.unique(nodes[:, 0])
ys = np.unique(nodes[:, 1])
zs = np.unique(nodes[:, 2])
nx, ny, nz = len(xs) - 1, len(ys) - 1, len(zs) - 1
# 2) Precompute face‐ID sets in C‐order flattening (Z fastest)
B = ny * nz
# Xmax: last block of size B
x_max = set(range((nx - 1) * B + 1, nx * B + 1))
# Ymax: in each x‐block, offset by (ny-1)*nz, vary z=0…nz-1
y_offset = (ny - 1) * nz
y_max = {
x * B + y_offset + (z + 1)
for x in range(nx)
for z in range(nz)
}
# Zmax: in each xy‐slice, ID = x*B + y*nz + nz
z_max = {
x * B + y * nz + nz
for x in range(nx)
for y in range(ny)
}
# 3) Helper to write one ELSET/SURFACE block
def _emit(face_ids, idx, label):
"""
Write a single surface element set and corresponding SURFACE block to a file
Parameters
----------
face_ids : set of int
Set of element IDs belonging to the surface
idx : int
Index of the surface for labeling
label : str
Label indicating the surface direction or type
Notes
-----
Element IDs are sorted and written in rows of 16 entries per line.
The ELSET is named "_Surf-{idx}_{label}" and the corresponding SURFACE block
is named "Surf-{idx}" with type ELEMENT.
"""
f.write(f'*ELSET, ELSET=_Surf-{idx}_{label}, internal, instance=PART-1-1\n')
# chunk 16 IDs per line
ids = sorted(face_ids)
for i in range(0, len(ids), 16):
chunk = ids[i:i + 16]
f.write(', '.join(str(e) for e in chunk) + ',\n')
f.write(f'*Surface, type=ELEMENT, name=Surf-{idx}\n')
f.write(f'_Surf-{idx}_{label}, {label}\n\n')
f.write('\n')
# 4) Emit all three faces
_emit(x_max, 1, 'S3') # X = max
_emit(y_max, 2, 'S2') # Y = max
_emit(z_max, 3, 'S6') # Z = max
def write_boundary_conditions():
"""
Write boundary conditions for the ABAQUS input file based on the loading direction
Notes
-----
Depending on the value of loading_direction ('x', 'y', 'z', 'xy', 'xz', 'yz'),
corresponding nodes are fixed in specific displacement/rotation directions.
Each boundary condition block is preceded by a comment indicating its name and type.
"""
f.write('** BOUNDARY CONDITIONS\n')
f.write('**\n')
if loading_direction == 'x':
f.write('** Name: F0yzFix Type: Displacement/Rotation\n')
f.write('*Boundary\nF0yz, 1, 1\n')
f.write('** Name: E0y0Fix Type: Displacement/Rotation\n')
f.write('*Boundary\nE0y0, 3, 3\n')
f.write('** Name: E00zFix Type: Displacement/Rotation\n')
f.write('*Boundary\nE00z, 2, 2\n')
elif loading_direction == 'y':
f.write('** Name: Fx0zFix Type: Displacement/Rotation\n')
f.write('*Boundary\nFx0z, 2, 2\n')
f.write('** Name: Ex00Fix Type: Displacement/Rotation\n')
f.write('*Boundary\nEx00, 3, 3\n')
f.write('** Name: E00zFix Type: Displacement/Rotation\n')
f.write('*Boundary\nE00z, 1, 1\n')
elif loading_direction == 'z':
f.write('** Name: Fxy0Fix Type: Displacement/Rotation\n')
f.write('*Boundary\nFxy0, 3, 3\n')
f.write('** Name: Ex00Fix Type: Displacement/Rotation\n')
f.write('*Boundary\nEx00, 2, 2\n')
f.write('** Name: E0y0Fix Type: Displacement/Rotation\n')
f.write('*Boundary\nE0y0, 1, 1\n')
elif loading_direction == 'xy':
f.write('** Name: Fx0zFix Type: Displacement/Rotation\n')
f.write('*Boundary\nFx0z, 1, 1\n')
f.write('*Boundary\nFx0z, 2, 2\n')
f.write('** Name: Fx1zFix Type: Displacement/Rotation\n')
f.write('*Boundary\nFx1z, 2, 2\n')
f.write('** Name: Fxy0Fix Type: Displacement/Rotation\n')
f.write('*Boundary\nFxy0, 3, 3\n')
elif loading_direction == 'xz':
f.write('** Name: F0yzFix Type: Displacement/Rotation\n')
f.write('*Boundary\nF0yz, 1, 1\n')
f.write('*Boundary\nF0yz, 3, 3\n')
f.write('** Name: Fx1zFix Type: Displacement/Rotation\n')
f.write('*Boundary\nFx1z, 1, 1\n')
f.write('** Name: Fx0zFix Type: Displacement/Rotation\n')
f.write('*Boundary\nFx0z, 2, 2\n')
elif loading_direction == 'yz':
f.write('** Name: Fxy0Fix Type: Displacement/Rotation\n')
f.write('*Boundary\nFxy0, 2, 2\n')
f.write('*Boundary\nFxy0, 3, 3\n')
f.write('** Name: Fxy1Fix Type: Displacement/Rotation\n')
f.write('*Boundary\nFxy1, 3, 3\n')
f.write('** Name: F0yzFix Type: Displacement/Rotation\n')
f.write('*Boundary\nF0yz, 1, 1\n')
def write_load():
"""
Write loading conditions to the ABAQUS input file based on load type and direction
Notes
-----
Supports 'strain' or 'stress' loading types. For strain loading, the displacement
is calculated using logarithmic strain based on the edge length and stretch ratio.
For stress loading, a pressure load is applied to the corresponding surface set.
The function uses mappings from loading_direction to node or surface sets and
writes appropriate *Boundary or *Dload blocks in the ABAQUS input file.
"""
if load_type == 'strain':
displacement_bc_map = {
'x' : ('F1yz', 1, 'disX' ),
'y' : ('Fx1z', 2, 'disY' ),
'z': ('Fxy1', 3, 'disZ' ),
'xy': ('Fx1z', 1, 'shearXY'),
'xz': ('F1yz', 3, 'shearXZ'),
'yz': ('Fxy1', 2, 'shearYZ'),
}
direction = loading_direction.lower()
if direction in displacement_bc_map.keys():
set_name, dof, bc_name = displacement_bc_map[direction]
# dof=1→X, 2→Y, 3→Z, 4→XY, 5→XZ, 6→YZ,
vstrain = float(mag)
true_strain = vstrain # ε = 0.2 for provided value =0.2
displacement = edge_lengths[direction] * (np.exp(true_strain) - 1.0) # Logarithmic strain
print(
f"Applied component: {direction}, "
f"Stretch ratio λ: {vstrain:.6f}, "
f"Log strain ε (per rule): {true_strain:.6f}, "
f"Edge length: {edge_lengths[direction]:.6f} mm, "
f"Displacement: {displacement:.6f} mm"
)
f.write(f'** Name: {bc_name} Type: Displacement/Rotation\n')
f.write('*Boundary\n')
f.write(f'{set_name}, {dof}, {dof}, {displacement:.6f}\n')
elif load_type == 'stress':
load_bc_map = {
'x' : ('SURF-1', 1, 'loadX' ),
'y' : ('SURF-2', 2, 'loadY' ),
'z' : ('SURF-3', 3, 'loadZ' ),
'xy': ('SURF-2', 1, 'loadXY'),
'xz': ('SURF-1', 3, 'loadXZ'),
'yz': ('SURF-3', 2, 'loadYZ'),
}
direction = loading_direction.lower()
if direction in load_bc_map:
set_name, dof, bc_name = load_bc_map[direction]
# dof=1→X, 2→Y, 3→Z, 4→XY, 5→XZ, 6→YZ,
vstress = mag
f.write('** LOADS\n')
f.write('**\n')
f.write(f'** Name: {bc_name} Type: Pressure\n')
f.write('*Dload\n')
f.write(f'{set_name}, P, {-vstress:.6f}\n')
def write_periodic_load():
"""
Write periodic loading conditions for stress or strain to the ABAQUS input file
Raises
------
ValueError
If loading_direction is unsupported for the selected load_type
or if load_type is not 'stress' or 'strain'
Notes
-----
For stress loading, the force is calculated as stress multiplied by the face area and
applied using *Cload. For strain loading, the displacement is calculated using
logarithmic strain based on edge length and applied using *Boundary.
Maps from loading_direction to node or surface IDs are used for correct application.
"""
if load_type == 'stress':
f.write('** BOUNDARY CONDITIONS\n')
f.write('** \n')
f.write('** LOADS \n')
f.write('** \n')
f.write('** Name: Load Type: Stress BC\n')
f.write('*Cload \n')
# map each index to its Abaqus CLOAD pattern
stress_map = {
'x': "V101,1", # σ_xx → force along x
'y': "V011,2", # σ_yy → force along y
'z': "V000,3", # σ_zz → force along z
'xy': "V011,1", # τ_xy → force along x (plane normal y)
'xz': "V101,3", # τ_xz → force along z (plane normal x)
'yz': "V000,2", # τ_yz → force along y (plane normal z)
}
d = loading_direction.lower()
if d not in stress_map or d not in face_areas:
raise ValueError(f"Unsupported loading_direction for stress: {loading_direction!r}")
area = float(face_areas[d])
stress_value = float(mag) # user-provided stress (units e.g., MPa)
force_value = stress_value * area # Force = Stress × Area
# Write the force to CLOAD
f.write(f'{stress_map[d]},{force_value}\n')
f.write('** \n')
elif load_type == 'strain':
f.write('** BOUNDARY CONDITIONS\n')
f.write('** \n')
f.write('** Displacements \n')
f.write('** \n')
f.write('** Name: Dis Type: Displacement BC\n')
f.write('*Boundary \n')
# map each index to (node, direction) for the *Boundary card
strain_map = {
'x': ("V101", 1),
'y': ("V011", 2),
'z': ("V000", 3),
'xy': ("V011", 1),
'xz': ("V101", 3),
'yz': ("V000", 2),
}
d = loading_direction.lower()
if d not in strain_map:
raise ValueError(f"Unsupported loading_direction for strain: {loading_direction!r}")
node, dof = strain_map[d]
#v_lambda = float(mag) # stretch ratio λ (e.g., 1.2 for +20%)
#eps_true = v_lambda - 1.0 # your convention
eps_true = float(mag) # true strain (e.g., 0.2 for +20%)
disp = edge_lengths[d] * (np.exp(eps_true) - 1.0)
f.write(f'{node}, {dof}, {dof}, {disp:.6f}\n')
else:
raise ValueError(f"Unknown load_type: {load_type!r}")
def _parse_boundary_conditions():
"""
Parse and validate the boundary_conditions dictionary
Parameters
----------
None
Uses the `boundary_conditions` variable from the outer scope
Returns
-------
apply_bc : bool
Whether boundary conditions should be applied
periodic_bc : bool
Whether periodic boundary conditions are requested
load_type : str
'strain' or 'stress' depending on type_bc or components_bc
load_case : str
'uni-axial' for normal loads or 'shear' for shear loads
direction : str
Loading direction: 'x', 'y', 'z', 'xy', 'xz', or 'yz'
magnitude : float
Magnitude of the load or displacement
Raises
------
ValueError
If `boundary_conditions` is not a dict, has invalid format, or contains invalid entries
Notes
-----
The function validates that only one degree of freedom has a non-zero magnitude (for stress)
or a numeric value (for strain), and maps the index to the corresponding direction and load type.
"""
if boundary_conditions is None:
return False, False, 'strain', 'uni-axial', 'x', 0.0
if not isinstance(boundary_conditions, dict):
raise ValueError("boundary_conditions must be a dict (or None).")
apply_bc = bool(boundary_conditions.get("apply_bc", False))
periodic_bc = bool(boundary_conditions.get("periodic_bc", False))
type_bc = boundary_conditions.get("type_bc", None)
components_bc = boundary_conditions.get("components_bc", None)
if components_bc is None or not (isinstance(components_bc, (list, tuple)) and len(components_bc) == 6):
raise ValueError("`components_bc` must be a list/tuple of length 6.")
if type_bc is not None:
t = str(type_bc).strip().lower()
if t in ("stress", "force"):
load_type = "stress"
if t == 'force':
logging.warning('Argument "force" will be interpreted as stress.')
elif t in ("strain", "displacement"):
load_type = "strain"
if t == 'displacement':
logging.warning('Argument "displacement" will be interpreted as strain.')
else:
raise ValueError("type_bc must be one of: 'stress'|'strain'.")
else:
load_type = "strain" if any(v == '*' for v in components_bc) else "stress"
logging.warning(f'load_type not specified, assuming {load_type}.')
# Validate + extract idx & magnitude
if load_type == "strain":
numeric_positions = []
for i, v in enumerate(components_bc):
if v == '*':
continue
if isinstance(v, (int, float)):
numeric_positions.append(i)
else:
raise ValueError(
f"Strain BC: entries must be either a number (the controlled DOF) or '*' (free). Got {v!r} at index {i}.")
if len(numeric_positions) != 1 or any(v != '*' for j, v in enumerate(components_bc) if j != (numeric_positions[0])):
raise ValueError("Strain BC: provide exactly one numeric magnitude and five '*' (free DOFs).")
idx = numeric_positions[0]
magnitude = float(components_bc[idx])
else: # stress/force
if any(v == '*' for v in components_bc):
raise ValueError("Stress BC: '*' is not allowed. Use 0 to indicate a free DOF.")
if not all(isinstance(v, (int, float)) for v in components_bc):
bad = [(i, v) for i, v in enumerate(components_bc) if not isinstance(v, (int, float))]
raise ValueError(f"Stress BC: all entries must be numeric. Bad entries: {bad}")
nz_indices = [i for i, v in enumerate(components_bc) if float(v) != 0.0]
if len(nz_indices) != 1 or any(float(v) != 0.0 for j, v in enumerate(components_bc) if j != nz_indices[0]):
raise ValueError(
"Stress BC: provide exactly one non-zero numeric magnitude and five zeros (free DOFs).")
idx = nz_indices[0]
magnitude = float(components_bc[idx])
# 3) Map index → direction & bc_type
dir_map = {0: 'x', 1: 'y', 2: 'z', 3: 'xy', 4: 'xz', 5: 'yz'}
direction = dir_map[idx]
load_case = 'uni-axial' if idx < 3 else 'shear'
return apply_bc, periodic_bc, load_type, load_case, direction, magnitude
# Parse boundary_conditions into control mode, loading mode, direction, and magnitude
if not boundary_conditions or not boundary_conditions.get("apply_bc", False):
apply_bc = False
periodic_bc = bool(boundary_conditions.get("periodic_bc", False)) if isinstance(boundary_conditions,
dict) else False
load_type = bc_type = loading_direction = None
mag = None
else:
apply_bc, periodic_bc, load_type, bc_type, loading_direction, mag = _parse_boundary_conditions()
print(
f"Control mode: {load_type} | "
f"Boundary condition type: {periodic_bc} | "
f"Loading mode: {bc_type} | "
f"Applied component: {loading_direction} | "
f"Magnitude value: {mag}"
)
print('')
print(f'Writing RVE as ABAQUS file "{file}"')
if gb_area is None:
if thermal:
print('Using brick element type C3D8T for coupled structural-thermal analysis.')
eltype = 'C3D8T'
else:
print('Using brick element type C3D8.')
eltype = 'C3D8'
else:
print('Using tet element type SFM3D4.')
eltype = 'SFM3D4'
if type(ialloy) is not list:
ialloy = [ialloy]
nall = len(ialloy)
ph_set = set()
nsets = NodeSets(nodes)
# Calculate RVE edge lengths and face areas in the selected output units
edge_lengths = {
'x': (max(nodes[:, 0]) - min(nodes[:, 0])) * scale_fact,
'y': (max(nodes[:, 1]) - min(nodes[:, 1])) * scale_fact,
'z': (max(nodes[:, 2]) - min(nodes[:, 2])) * scale_fact,
'xy': (max(nodes[:, 1]) - min(nodes[:, 1])) * scale_fact,
'xz': (max(nodes[:, 0]) - min(nodes[:, 0])) * scale_fact,
'yz': (max(nodes[:, 2]) - min(nodes[:, 2])) * scale_fact,
}
face_areas = {
# normal stresses
'x': edge_lengths['y'] * edge_lengths['z'], # yz face
'y': edge_lengths['x'] * edge_lengths['z'], # xz face
'z': edge_lengths['x'] * edge_lengths['y'], # xy face
# shear stresses (choose the plane normal consistent with mapping above)
'xy': edge_lengths['x'] * edge_lengths['z'], # plane normal y → xz face
'xz': edge_lengths['y'] * edge_lengths['z'], # plane normal x → yz face
'yz': edge_lengths['x'] * edge_lengths['y'], # plane normal z → xy face
}
#####################################
#### Start writing the .inp file ####
#####################################
with open(file, 'w') as f:
f.write(f'** Input file generated by kanapy version {__version__}\n')
f.write('** Nodal coordinates scale in mm\n')
f.write('*HEADING\n')
f.write('*PREPRINT,ECHO=NO,HISTORY=NO,MODEL=NO,CONTACT=NO\n')
f.write('**\n')
f.write('** PARTS\n')
f.write('**\n')
f.write('*Part, name=PART-1\n')
f.write('*Node\n')
for k, v in enumerate(nodes):
f.write('{0}, {1}, {2}, {3}\n'.format(k + 1, v[0] * scale_fact, v[1] * scale_fact, v[2] * scale_fact))
f.write('*ELEMENT, TYPE={0}\n'.format(eltype))
if gb_area is None:
for k, v in voxel_dict.items():
f.write('{0}, {1}, {2}, {3}, {4}, {5}, {6}, {7}, {8}\n'.format(
k, v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7]))
if dual_phase: # Why don't we write both phase sets and grain sets anyway?
write_phase_sets()
else:
write_grain_sets()
else:
fcList = {}
fcNum = 0
gr_fcs = defaultdict(list)
for gid, ginfo in enumerate(gb_area):
for fc, conn in ginfo.items():
if fc not in fcList.keys():
fcNum += 1
fcList[fc] = fcNum
f.write('%d,%d,%d,%d,%d\n' % (fcNum, conn[0], conn[1], conn[2], conn[3]))
gr_fcs[gid].append(fcNum)
elif fc in fcList.keys():
f.write('%d,%d,%d,%d,%d\n' % (fcList[fc], conn[0], conn[1], conn[2], conn[3]))
gr_fcs[gid].append(fcNum)
for gid, fcs in gr_fcs.items():
f.write('*ELSET, ELSET=GRAIN{}_SET\n'.format(gid))
for enum, el in enumerate(fcs, start=1):
if enum % 16 != 0:
if enum == len(fcs):
f.write('%d\n' % el)
else:
f.write('%d, ' % el)
else:
if enum == len(fcs):
f.write('%d\n' % el)
else:
f.write('%d\n' % el)
for gid, fcs in gr_fcs.items():
f.write('*SURFACE SECTION, ELSET=GRAIN{}_SET\n'.format(gid))
if periodic_bc and apply_bc:
f.write('**** ======================================================== \n')
f.write('**** Left to Right \n')
# LeftToRight
f.write('**** \n')
f.write('**** X-DIR \n')
for i in range(len(nsets.F0yzP)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.F1yzP[i] + 1) + ',1, 1 \n')
f.write(str(nsets.F0yzP[i] + 1) + ',1,-1 \n')
f.write(str(nsets.V101 + 1) + ',1,-1 \n')
f.write(str(nsets.V001 + 1) + ',1, 1 \n')
f.write('**** \n')
f.write('**** Y-DIR \n')
for i in range(len(nsets.F0yzP)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.F1yzP[i] + 1) + ',2, 1 \n')
f.write(str(nsets.F0yzP[i] + 1) + ',2,-1 \n')
f.write(str(nsets.V101 + 1) + ',2,-1 \n')
f.write(str(nsets.V001 + 1) + ',2, 1 \n')
f.write('**** \n')
f.write('**** Z-DIR \n')
for i in range(len(nsets.F0yzP)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.F1yzP[i] + 1) + ',3, 1 \n')
f.write(str(nsets.F0yzP[i] + 1) + ',3,-1 \n')
f.write(str(nsets.V101 + 1) + ',3,-1 \n')
f.write(str(nsets.V001 + 1) + ',3, 1 \n')
f.write('**** Bottom to Top \n')
# BottomToTop
f.write('**** X-DIR \n')
for i in range(len(nsets.Fx0zP)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.Fx0zP[i] + 1) + ',1, 1 \n')
f.write(str(nsets.Fx1zP[i] + 1) + ',1,-1 \n')
f.write(str(nsets.V001 + 1) + ',1,-1 \n')
f.write(str(nsets.V011 + 1) + ',1, 1 \n')
f.write('**** \n')
f.write('**** Y-DIR \n')
for i in range(len(nsets.Fx0zP)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.Fx0zP[i] + 1) + ',2, 1 \n')
f.write(str(nsets.Fx1zP[i] + 1) + ',2,-1 \n')
f.write(str(nsets.V001 + 1) + ',2,-1 \n')
f.write(str(nsets.V011 + 1) + ',2, 1 \n')
f.write('**** \n')
f.write('**** Z-DIR \n')
for i in range(len(nsets.Fx0zP)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.Fx0zP[i] + 1) + ',3, 1 \n')
f.write(str(nsets.Fx1zP[i] + 1) + ',3,-1 \n')
f.write(str(nsets.V001 + 1) + ',3,-1 \n')
f.write(str(nsets.V011 + 1) + ',3, 1 \n')
f.write('**** Front to Rear \n')
# FrontToRear
f.write('**** \n')
f.write('**** X-DIR \n')
for i in range(len(nsets.Fxy0P)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.Fxy0P[i] + 1) + ',1,-1 \n')
f.write(str(nsets.Fxy1P[i] + 1) + ',1, 1 \n')
f.write(str(nsets.V001 + 1) + ',1,-1 \n')
f.write(str(nsets.V000 + 1) + ',1, 1 \n')
f.write('**** \n')
f.write('**** Y-DIR \n')
for i in range(len(nsets.Fxy0P)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.Fxy0P[i] + 1) + ',2,-1 \n')
f.write(str(nsets.Fxy1P[i] + 1) + ',2, 1 \n')
f.write(str(nsets.V001 + 1) + ',2,-1 \n')
f.write(str(nsets.V000 + 1) + ',2, 1 \n')
f.write('**** \n')
f.write('**** Z-DIR \n')
for i in range(len(nsets.Fxy0P)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.Fxy0P[i] + 1) + ',3,-1 \n')
f.write(str(nsets.Fxy1P[i] + 1) + ',3, 1 \n')
f.write(str(nsets.V001 + 1) + ',3,-1 \n')
f.write(str(nsets.V000 + 1) + ',3, 1 \n')
f.write('**** ======================================================== \n')
f.write('**** Edges\n')
f.write('**** \n')
# Edges in x-y Plane
# Right top edge to left top edge
f.write('**** X-DIR \n')
for i in range(len(nsets.E11zP)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.E11zP[i] + 1) + ',1, 1 \n')
f.write(str(nsets.E01zP[i] + 1) + ',1,-1 \n')
f.write(str(nsets.V101 + 1) + ',1,-1 \n')
f.write(str(nsets.V001 + 1) + ',1, 1 \n')
f.write('**** Y-DIR \n')
for i in range(len(nsets.E11zP)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.E11zP[i] + 1) + ',2, 1 \n')
f.write(str(nsets.E01zP[i] + 1) + ',2,-1 \n')
f.write(str(nsets.V101 + 1) + ',2,-1 \n')
f.write(str(nsets.V001 + 1) + ',2, 1 \n')
f.write('**** Z-DIR \n')
for i in range(len(nsets.E11zP)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.E11zP[i] + 1) + ',3, 1 \n')
f.write(str(nsets.E01zP[i] + 1) + ',3,-1 \n')
f.write(str(nsets.V101 + 1) + ',3,-1 \n')
f.write(str(nsets.V001 + 1) + ',3, 1 \n')
# Right bottom edge to left bottom edge
f.write('**** X-DIR \n')
for i in range(len(nsets.E10zP)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.E10zP[i] + 1) + ',1, 1 \n')
f.write(str(nsets.E00zP[i] + 1) + ',1,-1 \n')
f.write(str(nsets.V101 + 1) + ',1,-1 \n')
f.write(str(nsets.V001 + 1) + ',1, 1 \n')
f.write('**** Y-DIR \n')
for i in range(len(nsets.E10zP)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.E10zP[i] + 1) + ',2, 1 \n')
f.write(str(nsets.E00zP[i] + 1) + ',2,-1 \n')
f.write(str(nsets.V101 + 1) + ',2,-1 \n')
f.write(str(nsets.V001 + 1) + ',2, 1 \n')
f.write('**** Z-DIR \n')
for i in range(len(nsets.E10zP)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.E10zP[i] + 1) + ',3, 1 \n')
f.write(str(nsets.E00zP[i] + 1) + ',3,-1 \n')
f.write(str(nsets.V101 + 1) + ',3,-1 \n')
f.write(str(nsets.V001 + 1) + ',3, 1 \n')
# Left top edge to left bottom edge
f.write('**** X-DIR \n')
for i in range(len(nsets.E01zP)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.E01zP[i] + 1) + ',1, 1 \n')
f.write(str(nsets.E00zP[i] + 1) + ',1,-1 \n')
f.write(str(nsets.V011 + 1) + ',1,-1 \n')
f.write(str(nsets.V001 + 1) + ',1, 1 \n')
f.write('**** Y-DIR \n')
for i in range(len(nsets.E01zP)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.E01zP[i] + 1) + ',2, 1 \n')
f.write(str(nsets.E00zP[i] + 1) + ',2,-1 \n')
f.write(str(nsets.V011 + 1) + ',2,-1 \n')
f.write(str(nsets.V001 + 1) + ',2, 1 \n')
f.write('**** Z-DIR \n')
for i in range(len(nsets.E01zP)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.E01zP[i] + 1) + ',3, 1 \n')
f.write(str(nsets.E00zP[i] + 1) + ',3,-1 \n')
f.write(str(nsets.V011 + 1) + ',3,-1 \n')
f.write(str(nsets.V001 + 1) + ',3, 1 \n')
# Edges in y-z Plane
# Top back edge to top front edge
f.write('**** X-DIR \n')
for i in range(len(nsets.Ex10P)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.Ex10P[i] + 1) + ',1, 1 \n')
f.write(str(nsets.Ex11P[i] + 1) + ',1,-1 \n')
f.write(str(nsets.V000 + 1) + ',1,-1 \n')
f.write(str(nsets.V001 + 1) + ',1, 1 \n')
f.write('**** Y-DIR \n')
for i in range(len(nsets.Ex10P)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.Ex10P[i] + 1) + ',2, 1 \n')
f.write(str(nsets.Ex11P[i] + 1) + ',2,-1 \n')
f.write(str(nsets.V000 + 1) + ',2,-1 \n')
f.write(str(nsets.V001 + 1) + ',2, 1 \n')
f.write('**** Z-DIR \n')
for i in range(len(nsets.Ex10P)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.Ex10P[i] + 1) + ',3, 1 \n')
f.write(str(nsets.Ex11P[i] + 1) + ',3,-1 \n')
f.write(str(nsets.V000 + 1) + ',3,-1 \n')
f.write(str(nsets.V001 + 1) + ',3, 1 \n')
# Bottom back edge to bottom front edge
f.write('**** X-DIR \n')
for i in range(len(nsets.Ex00P)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.Ex00P[i] + 1) + ',1, 1 \n')
f.write(str(nsets.Ex01P[i] + 1) + ',1,-1 \n')
f.write(str(nsets.V000 + 1) + ',1,-1 \n')
f.write(str(nsets.V001 + 1) + ',1, 1 \n')
f.write('**** Y-DIR \n')
for i in range(len(nsets.Ex00P)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.Ex00P[i] + 1) + ',2, 1 \n')
f.write(str(nsets.Ex01P[i] + 1) + ',2,-1 \n')
f.write(str(nsets.V000 + 1) + ',2,-1 \n')
f.write(str(nsets.V001 + 1) + ',2, 1 \n')
f.write('**** Z-DIR \n')
for i in range(len(nsets.Ex00P)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.Ex00P[i] + 1) + ',3, 1 \n')
f.write(str(nsets.Ex01P[i] + 1) + ',3,-1 \n')
f.write(str(nsets.V000 + 1) + ',3,-1 \n')
f.write(str(nsets.V001 + 1) + ',3, 1 \n')
# top front edge to bottom front edge
f.write('**** X-DIR \n')
for i in range(len(nsets.Ex11P)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.Ex11P[i] + 1) + ',1, 1 \n')
f.write(str(nsets.Ex01P[i] + 1) + ',1,-1 \n')
f.write(str(nsets.V011 + 1) + ',1,-1 \n')
f.write(str(nsets.V001 + 1) + ',1, 1 \n')
f.write('**** Y-DIR \n')
for i in range(len(nsets.Ex11P)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.Ex11P[i] + 1) + ',2, 1 \n')
f.write(str(nsets.Ex01P[i] + 1) + ',2,-1 \n')
f.write(str(nsets.V011 + 1) + ',2,-1 \n')
f.write(str(nsets.V001 + 1) + ',2, 1 \n')
f.write('**** Z-DIR \n')
for i in range(len(nsets.Ex11P)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.Ex11P[i] + 1) + ',3, 1 \n')
f.write(str(nsets.Ex01P[i] + 1) + ',3,-1 \n')
f.write(str(nsets.V011 + 1) + ',3,-1 \n')
f.write(str(nsets.V001 + 1) + ',3, 1 \n')
# Edges in x-z Plane
# Rear right edge to rear left edge
f.write('**** X-DIR \n')
for i in range(len(nsets.E1y0P)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.E1y0P[i] + 1) + ',1, 1 \n')
f.write(str(nsets.E0y0P[i] + 1) + ',1,-1 \n')
f.write(str(nsets.V101 + 1) + ',1,-1 \n')
f.write(str(nsets.V001 + 1) + ',1, 1 \n')
f.write('**** Y-DIR \n')
for i in range(len(nsets.E1y0P)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.E1y0P[i] + 1) + ',2, 1 \n')
f.write(str(nsets.E0y0P[i] + 1) + ',2,-1 \n')
f.write(str(nsets.V101 + 1) + ',2,-1 \n')
f.write(str(nsets.V001 + 1) + ',2, 1 \n')
f.write('**** Z-DIR \n')
for i in range(len(nsets.E1y0P)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.E1y0P[i] + 1) + ',3, 1 \n')
f.write(str(nsets.E0y0P[i] + 1) + ',3,-1 \n')
f.write(str(nsets.V101 + 1) + ',3,-1 \n')
f.write(str(nsets.V001 + 1) + ',3, 1 \n')
# Front right edge to front left edge
f.write('**** X-DIR \n')
for i in range(len(nsets.E1y1P)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.E1y1P[i] + 1) + ',1, 1 \n')
f.write(str(nsets.E0y1P[i] + 1) + ',1,-1 \n')
f.write(str(nsets.V101 + 1) + ',1,-1 \n')
f.write(str(nsets.V001 + 1) + ',1, 1 \n')
f.write('**** Y-DIR \n')
for i in range(len(nsets.E1y1P)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.E1y1P[i] + 1) + ',2, 1 \n')
f.write(str(nsets.E0y1P[i] + 1) + ',2,-1 \n')
f.write(str(nsets.V101 + 1) + ',2,-1 \n')
f.write(str(nsets.V001 + 1) + ',2, 1 \n')
f.write('**** Z-DIR \n')
for i in range(len(nsets.E1y1P)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.E1y1P[i] + 1) + ',3, 1 \n')
f.write(str(nsets.E0y1P[i] + 1) + ',3,-1 \n')
f.write(str(nsets.V101 + 1) + ',3,-1 \n')
f.write(str(nsets.V001 + 1) + ',3, 1 \n')
# Top front edge to bottom front edge
f.write('**** X-DIR \n')
for i in range(len(nsets.E0y0P)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.E0y0P[i] + 1) + ',1, 1 \n')
f.write(str(nsets.E0y1P[i] + 1) + ',1,-1 \n')
f.write(str(nsets.V000 + 1) + ',1,-1 \n')
f.write(str(nsets.V001 + 1) + ',1, 1 \n')
f.write('**** Y-DIR \n')
for i in range(len(nsets.E0y0P)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.E0y0P[i] + 1) + ',2, 1 \n')
f.write(str(nsets.E0y1P[i] + 1) + ',2,-1 \n')
f.write(str(nsets.V000 + 1) + ',2,-1 \n')
f.write(str(nsets.V001 + 1) + ',2, 1 \n')
f.write('**** Z-DIR \n')
for i in range(len(nsets.E0y0P)):
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.E0y0P[i] + 1) + ',3, 1 \n')
f.write(str(nsets.E0y1P[i] + 1) + ',3,-1 \n')
f.write(str(nsets.V000 + 1) + ',3,-1 \n')
f.write(str(nsets.V001 + 1) + ',3, 1 \n')
f.write('**** ======================================================== \n')
f.write('**** Corners \n')
# V3 (V111) to V4 (V011)
f.write('**** X-DIR \n')
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.V111 + 1) + ',1, 1 \n')
f.write(str(nsets.V011 + 1) + ',1,-1 \n')
f.write(str(nsets.V101 + 1) + ',1,-1 \n')
f.write(str(nsets.V001 + 1) + ',1, 1 \n')
f.write('**** y-DIR \n')
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.V111 + 1) + ',2, 1 \n')
f.write(str(nsets.V011 + 1) + ',2,-1 \n')
f.write(str(nsets.V101 + 1) + ',2,-1 \n')
f.write(str(nsets.V001 + 1) + ',2, 1 \n')
f.write('**** z-DIR \n')
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.V111 + 1) + ',3, 1 \n')
f.write(str(nsets.V011 + 1) + ',3,-1 \n')
f.write(str(nsets.V101 + 1) + ',3,-1 \n')
f.write(str(nsets.V001 + 1) + ',3, 1 \n')
# H4 (V010) to V4 (V011)
f.write('**** X-DIR \n')
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.V010 + 1) + ',1, 1 \n')
f.write(str(nsets.V011 + 1) + ',1,-1 \n')
f.write(str(nsets.V000 + 1) + ',1,-1 \n')
f.write(str(nsets.V001 + 1) + ',1, 1 \n')
f.write('**** y-DIR \n')
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.V010 + 1) + ',2, 1 \n')
f.write(str(nsets.V011 + 1) + ',2,-1 \n')
f.write(str(nsets.V000 + 1) + ',2,-1 \n')
f.write(str(nsets.V001 + 1) + ',2, 1 \n')
f.write('**** z-DIR \n')
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.V010 + 1) + ',3, 1 \n')
f.write(str(nsets.V011 + 1) + ',3,-1 \n')
f.write(str(nsets.V000 + 1) + ',3,-1 \n')
f.write(str(nsets.V001 + 1) + ',3, 1 \n')
# H3 (V110) to V3 (V111)
f.write('**** X-DIR \n')
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.V110 + 1) + ',1, 1 \n')
f.write(str(nsets.V111 + 1) + ',1,-1 \n')
f.write(str(nsets.V000 + 1) + ',1,-1 \n')
f.write(str(nsets.V001 + 1) + ',1, 1 \n')
f.write('**** y-DIR \n')
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.V110 + 1) + ',2, 1 \n')
f.write(str(nsets.V111 + 1) + ',2,-1 \n')
f.write(str(nsets.V000 + 1) + ',2,-1 \n')
f.write(str(nsets.V001 + 1) + ',2, 1 \n')
f.write('**** z-DIR \n')
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.V110 + 1) + ',3, 1 \n')
f.write(str(nsets.V111 + 1) + ',3,-1 \n')
f.write(str(nsets.V000 + 1) + ',3,-1 \n')
f.write(str(nsets.V001 + 1) + ',3, 1 \n')
# H2 (V100) to V2 (V101)
f.write('**** X-DIR \n')
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.V100 + 1) + ',1, 1 \n')
f.write(str(nsets.V101 + 1) + ',1,-1 \n')
f.write(str(nsets.V000 + 1) + ',1,-1 \n')
f.write(str(nsets.V001 + 1) + ',1, 1 \n')
f.write('**** y-DIR \n')
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.V100 + 1) + ',2, 1 \n')
f.write(str(nsets.V101 + 1) + ',2,-1 \n')
f.write(str(nsets.V000 + 1) + ',2,-1 \n')
f.write(str(nsets.V001 + 1) + ',2, 1 \n')
f.write('**** z-DIR \n')
f.write('*Equation \n')
f.write('4 \n')
f.write(str(nsets.V100 + 1) + ',3, 1 \n')
f.write(str(nsets.V101 + 1) + ',3,-1 \n')
f.write(str(nsets.V000 + 1) + ',3,-1 \n')
f.write(str(nsets.V001 + 1) + ',3, 1 \n')
f.write('*End Part\n')
f.write('**\n')
f.write('** ASSEMBLY\n')
f.write('**\n')
f.write('*Assembly, name=Assembly\n')
f.write('**\n')
f.write('*Instance, name=PART-1-1, part=PART-1\n')
f.write('*End Instance\n')
f.write('**\n')
f.write('** DEFINE NODE SETS\n')
f.write('** 1. VERTICES\n')
f.write(f'*Nset, nset=V000, instance=PART-1-1\n')
f.write(f'{nsets.V000 + 1}\n')
f.write(f'*Nset, nset=V001, instance=PART-1-1\n')
f.write(f'{nsets.V001 + 1}\n')
f.write(f'*Nset, nset=V010, instance=PART-1-1\n')
f.write(f'{nsets.V010 + 1}\n')
f.write(f'*Nset, nset=V100, instance=PART-1-1\n')
f.write(f'{nsets.V100 + 1}\n')
f.write(f'*Nset, nset=V011, instance=PART-1-1\n')
f.write(f'{nsets.V011 + 1}\n')
f.write(f'*Nset, nset=V101, instance=PART-1-1\n')
f.write(f'{nsets.V101 + 1}\n')
f.write(f'*Nset, nset=V110, instance=PART-1-1\n')
f.write(f'{nsets.V110 + 1}\n')
f.write(f'*Nset, nset=V111, instance=PART-1-1\n')
f.write(f'{nsets.V111 + 1}\n')
f.write('*Nset, nset=Vertices, instance=PART-1-1\n')
f.write(f'{nsets.V000 + 1}, {nsets.V100 + 1}, {nsets.V010 + 1}, {nsets.V001 + 1}, {nsets.V011 + 1}, '
f'{nsets.V101 + 1}, {nsets.V110 + 1}, {nsets.V111 + 1}\n')
if periodic_bc:
f.write('*Nset, nset=VerticesPeriodic, instance=PART-1-1\n')
f.write(f'{nsets.V001 + 1}, {nsets.V101 + 1}, {nsets.V011 + 1}, {nsets.V000 + 1} \n')
f.write('** 2. EDGES\n')
write_node_set('*Nset, nset=Ex00, instance=PART-1-1\n', nsets.Ex00)
write_node_set('*Nset, nset=Ex01, instance=PART-1-1\n', nsets.Ex01)
write_node_set('*Nset, nset=Ex10, instance=PART-1-1\n', nsets.Ex10)
write_node_set('*Nset, nset=Ex11, instance=PART-1-1\n', nsets.Ex11)
write_node_set('*Nset, nset=E0y0, instance=PART-1-1\n', nsets.E0y0)
write_node_set('*Nset, nset=E0y1, instance=PART-1-1\n', nsets.E0y1)
write_node_set('*Nset, nset=E1y0, instance=PART-1-1\n', nsets.E1y0)
write_node_set('*Nset, nset=E1y1, instance=PART-1-1\n', nsets.E1y1)
write_node_set('*Nset, nset=E00z, instance=PART-1-1\n', nsets.E00z)
write_node_set('*Nset, nset=E01z, instance=PART-1-1\n', nsets.E01z)
write_node_set('*Nset, nset=E10z, instance=PART-1-1\n', nsets.E10z)
write_node_set('*Nset, nset=E11z, instance=PART-1-1\n', nsets.E11z)
f.write('** 3. FACES\n')
write_node_set(f'*Nset, nset=Fxy0, instance=PART-1-1\n', nsets.Fxy0)
write_node_set(f'*Nset, nset=Fxy1, instance=PART-1-1\n', nsets.Fxy1)
write_node_set(f'*Nset, nset=Fx0z, instance=PART-1-1\n', nsets.Fx0z)
write_node_set(f'*Nset, nset=Fx1z, instance=PART-1-1\n', nsets.Fx1z)
write_node_set(f'*Nset, nset=F0yz, instance=PART-1-1\n', nsets.F0yz)
write_node_set(f'*Nset, nset=F1yz, instance=PART-1-1\n', nsets.F1yz)
f.write('** 4. FACES_SETS\n')
write_surface_sets()
f.write('**\n')
f.write('*End Assembly\n')
f.write('**\n')
############################
### Creating Material
############################
f.write('** MATERIALS\n')
f.write('**\n')
# track whether we ever did an include-file inside the loop
did_include = False
for pid in sorted(ph_set):
if files is not None and props_file is not None and files[pid] is None:
if pid == 0 and grain_phase_dict and grain_phase_dict.get(0) == 0:
# Keep the legacy matrix material visible in Abaqus/CAE.
f.write('*Material, name=PHASE0_MAT\n**\n')
continue
f.write('*Material, name=PHASE{}_MAT\n'.format(pid))
if props_file is not None:
include = os.fspath(files[pid])
if any(char in include for char in ('\n', '\r', '"')):
raise ValueError('Invalid include filename.')
f.write(f'*Include, input="{include}"\n')
elif phase_props:
props = phase_props.get(pid, {})
# inline properties as before
if 'damage_init' in props:
di = props['damage_init']
f.write('*Damage Initiation, criterion={}\n'.format(di['criterion']))
f.write(' {}\n'.format(', '.join(map(str, di['values']))))
if 'damage_evol' in props:
de = props['damage_evol']
f.write('*Damage Evolution, type={}\n'.format(de['type']))
f.write(' {}\n'.format(', '.join(map(str, de['values']))))
if 'elastic' in props:
f.write('*Elastic\n')
f.write(' {}\n'.format(', '.join(map(str, props['elastic']))))
if 'plastic' in props:
f.write('*Plastic\n')
for sigma, eps in props['plastic']:
f.write(' {}, {}\n'.format(sigma, eps))
else:
# no inline props → include from file
if dual_phase:
# per‐phase file for each pid
f.write('*Include, input=Material{}.inp\n'.format(pid))
did_include = True
# else: defer the single‐file include until after the loop
f.write('**\n')
# if this wasn’t a dual‐phase run, do one global include once:
if not dual_phase and not did_include and (cp is None or any(cp)):
# strip off last 8 chars (e.g. “_mesh.inp”) and append “mat.inp”
base = file[:-8]
f.write('*Include, input={}mat.inp\n'.format(base))
f.write('**\n')
"""
Previous Material Section
f.write('**__________________________________________________________________\n')
f.write('** MATERIALS\n')
f.write('**\n')
if dual_phase:
for i in ph_set:
f.write('**\n')
f.write('*Material, name=PHASE{}_MAT\n'.format(i))
f.write('**Include, input=Material{}.inp\n'.format(i))
f.write('**\n')
else:
for i in ph_set:
f.write('**\n')
f.write('*Material, name=PHASE{}_MAT\n'.format(i))
f.write('**\n')
f.write('**Include, input={}mat.inp\n'.format(file[0:-8]))
f.write('**\n')
f.write('**__________________________________________________________________')
"""
##########################################################
### Creating non-Periodic or Periodic Boundary Conditions
##########################################################
if apply_bc:
if not periodic_bc :
write_boundary_conditions()
if periodic_bc:
f.write('** BOUNDARY CONDITIONS\n')
f.write('** \n')
f.write('*Boundary \n')
f.write('V001,1 \n')
f.write('V001,2 \n')
f.write('V001,3 \n')
f.write('V101,2 \n')
f.write('V000,1 \n')
f.write('V011,3 \n')
if bc_type.lower() == 'shear':
f.write('V000,3 \n')
f.write('V101,1 \n')
f.write('V011,2 \n')
if loading_direction.lower() == 'xy' :
f.write('V101,3 \n')
f.write('V000,2 \n')
elif loading_direction.lower() == 'xz' :
f.write('V011,1 \n')
f.write('V000,2 \n')
elif loading_direction.lower() == 'yz' :
f.write('V011,1 \n')
f.write('V101,3 \n')
f.write('** \n')
f.write('** \n')
############################
### Creating Step
############################
if (any(cp) if cp is not None else crystal_plasticity): # Using crystal plasticity Umat
f.write('**\n')
f.write('** STEP: Loading\n')
f.write('**\n')
f.write('*Step, name=Loading, nlgeom=YES, inc=500000, unsymm=YES, solver=ITERATIVE\n')
f.write('*Static\n')
f.write('0.001, 1., 1e-10, 0.01\n')
f.write('**\n')
f.write('**\n')
f.write('*CONTROLS, PARAMETER=TIME INCREMENTATION\n')
f.write('35, 50, 9, 50, 28, 5, 12, 45\n')
f.write('**\n')
f.write('*CONTROLS, PARAMETERS=LINE SEARCH\n')
f.write('10\n')
f.write('** Originally that was SOLVER CONTROL\n')
f.write('*SOLVER CONTROL\n')
f.write('1e-5,200,\n')
else: # Using build-in functions
f.write('**\n')
f.write('** STEP: Loading \n')
f.write('** \n')
f.write('*Step, name=Loading, nlgeom=YES, inc=500000 \n')
f.write('*Static \n')
f.write('0.001, 1., 1e-6, 0.02 \n')
f.write('** \n')
f.write('** \n')
#################################
### Creating Load
#################################
if apply_bc:
if not periodic_bc:
write_load()
f.write('** \n')
elif periodic_bc:
write_periodic_load()
#################################
### Creating Fild Output
#################################
f.write('** OUTPUT REQUESTS \n')
f.write('** \n')
f.write('*Restart, write, frequency=0 \n')
f.write('** \n')
f.write('** FIELD OUTPUT: F-Output-1 \n')
f.write('** \n')
f.write('*Output, field \n')
f.write('*Node Output \n')
f.write('CF, COORD, RF, U \n')
f.write('** \n')
f.write('** FIELD OUTPUT: F-Output-2 \n')
f.write('** \n')
f.write('*Element Output, directions=YES \n')
f.write('LE, PE, PEEQ, S, SDEG, SDV \n')
f.write('*Output, history, frequency=0 \n')
f.write('** \n')
f.write('** HISTORY OUTPUT: H-Output-1 \n')
f.write('** \n')
f.write('*Output, history, variable=PRESELECT \n')
f.write('*End Step \n')
print('---->DONE! \n')
return
[docs]
def writeAbaqusMat(
ialloy: Any,
angles: Any,
file=None, path='./',
grain_phase_dict=None,
nsdv=360, props_file=None, crystal_plasticity=None):
"""
Export Euler angles to Abaqus input deck that can be included in the _geom.inp file. If
parameter "grain_phase_dict" is given, the phase number for each grain will be used to select
the corresponding ialloy from a list. If the list ialloy is shorter than the number of phases in
grain_phase_dict, no angles for phases with no corresponding ialloy will be written.
Parameters:
-----------
ialloy : int or list
Identifier, alloy number in ICAMS CP-UMAT: mod_alloys.f
angles : dict or (N, 3)-ndarray
Dict with Euler angles for each grain or array with number of N rows (= number of grains) and
three columns phi1, Phi, phi2
file : str
Filename, optional (default: None)
path : str
Path to save file, option (default: './')
grain_phase_dict: dict
Dict with phase for each grain, optional (default: None)
props_file : str, path-like or list, optional
Shared numeric constants include, or one include per alloy/phase.
Relative paths are resolved against the output material file directory.
Enables the eight-value CP-UMAT header and automatically counts constants.
Files must have eight values per line except the last; comments are allowed.
Files are referenced without modification. None preserves legacy output.
crystal_plasticity : bool, list of bool, or None, optional
Write only grains in CP phases. With multiphase props_file, both options
must be lists with one entry per phase; each True requires an include.
False phases (including None includes) are handled by export2abaqus in
_geom.inp. Grain 0 is reserved for standard plasticity in phase 0.
Omit both options to retain the legacy alloy-based selection.
nsdv : int
Number of state dependant variables, optional (default: 360)
"""
validate_matrix_mapping(grain_phase_dict)
if props_file is not None and ialloy is None:
ialloy = 0
scalar_alloy = not isinstance(ialloy, list)
if scalar_alloy:
ialloy = [ialloy]
nall = len(ialloy)
nphases = max((grain_phase_dict or {}).values(), default=0) + 1
nphases = max(nphases, nall)
if isinstance(crystal_plasticity, list):
nphases = max(nphases, len(crystal_plasticity))
files, cp = _abaqus_phase_options(props_file, crystal_plasticity, nphases)
if props_file is not None and scalar_alloy:
ialloy *= nphases
nall = nphases
if grain_phase_dict and 0 in grain_phase_dict and cp is not None and cp[0]:
raise ValueError('Grain 0 requires standard plasticity in phase 0.')
if type(angles) is not dict:
# converting (N, 3) ndarray to dict
gr_ori_dict = dict()
for igr, ori in enumerate(angles):
gr_ori_dict[igr + 1] = ori
else:
gr_ori_dict = angles
nitem = len(gr_ori_dict.keys())
if file is None:
file = f'abq_px_{nitem}_mat.inp'
path = os.path.normpath(path)
file = os.path.join(path, file)
includes = None
if props_file is not None:
files = list(props_file) if isinstance(props_file, (list, tuple)) else [props_file] * nall
if len(files) != nall:
raise ValueError('props_file must contain one include file per selector.')
includes = []
for pid, include in enumerate(files):
if cp is not None and not cp[pid]:
includes.append(None)
continue
include = os.fspath(include)
if any(char in include for char in ('\n', '\r', '"')):
raise ValueError('Invalid include filename.')
source = os.path.join(os.path.dirname(file), include)
rows = []
with open(source, encoding='utf-8-sig') as stream:
for line in stream:
line = line.strip()
if not line or line.startswith('**'):
continue
values = line.rstrip(',').split(',')
try:
numbers = [float(value.strip().replace('D', 'E').replace('d', 'e'))
for value in values]
except ValueError as exc:
raise ValueError(f'{source}: expected numeric constants only.') from exc
if not all(np.isfinite(numbers)) or not 1 <= len(numbers) <= 8:
raise ValueError(f'{source}: expected 1 to 8 finite constants per line.')
rows.append(numbers)
if not rows or any(len(row) != 8 for row in rows[:-1]):
raise ValueError(f'{source}: use eight constants per line except the final line.')
includes.append((include, 8 + sum(map(len, rows))))
for gid, orientation in gr_ori_dict.items():
pid = (grain_phase_dict or {}).get(gid, 0)
if gid == 0 or pid >= nall or (cp is not None and not cp[pid]):
continue
if orientation is None or np.asarray(orientation).shape != (3,):
raise ValueError(f'Missing or invalid orientation for CP grain {gid}.')
if ialloy[pid] is None:
raise ValueError(f'Missing alloy selector for CP phase {pid}.')
with open(file, 'w') as f:
f.write('**\n')
f.write('** MATERIALS\n')
f.write('**\n')
for igr, ori in gr_ori_dict.items():
if igr==0:
continue # Grain 0 is reserved for matrix phase / unidentified phases
if grain_phase_dict is None:
ip = 0
else:
ip = grain_phase_dict[igr]
if ip > nall - 1:
continue
if cp is not None and not cp[ip]:
continue
f.write('*Material, name=GRAIN{}_MAT\n'.format(igr))
f.write('*Depvar\n')
f.write(' {}\n'.format(nsdv))
if includes is not None:
include, count = includes[ip]
f.write(f'*User Material, constants={count}\n')
f.write(f'{float(ialloy[ip])}, {ori[0]}, {ori[1]}, {ori[2]}, 0., 0., 0., 0.\n')
f.write(f'*Include, input="{include}"\n')
continue
f.write('*User Material, constants=4\n')
f.write('{}, {}, {}, {}\n'.format(float(ialloy[ip]),
ori[0], ori[1], ori[2]))
return
[docs]
def pickle2microstructure(
file: Union[str, os.PathLike[str]],
path: Union[str, os.PathLike[str]] = './') -> Any:
"""
Load a pickled microstructure object from disk.
Parameters
----------
file : str
Filename of the pickled microstructure.
path : str, optional
Directory where the pickle file is stored (default: './').
Returns
-------
pcl : Material
Unpickled microstructure object.
Raises
------
ValueError
If `file` is None.
FileNotFoundError
If the specified pickle file does not exist.
pickle.UnpicklingError
If the file cannot be unpickled (e.g., corrupted or invalid format).
"""
import pickle
if file is None:
raise ValueError('Name for pickled microstructure must be given.')
fname = os.path.join(path, file)
with open(fname, 'rb') as inp:
pcl = pickle.load(inp)
return pcl
[docs]
def import_voxels(
file: Union[str, os.PathLike[str]],
path: Union[str, os.PathLike[str]] = './') -> Any:
"""
Import a voxelized microstructure from a JSON file and reconstruct
a `Microstructure` object for further analysis or simulation.
The voxel file typically contains:
- Grain-wise orientation, phase, and shape information
- Global model metadata (size, periodicity, phase names, etc.)
- Optionally, mesh node coordinates and voxel connectivity
Parameters
----------
file : str
Filename of the voxelized microstructure JSON file to import.
path : str, optional
Path where the JSON file is located (default: './').
Returns
-------
ms : Microstructure
Fully reconstructed microstructure object containing:
- Grain geometry, phases, orientations
- Voxel-level mesh and connectivity
- Simulation box and RVE setup
- Volume fractions and statistics
Raises
------
ValueError
If `file` is None.
FileNotFoundError
If the specified voxel file cannot be found.
KeyError
If required keys (e.g. 'Data', 'Model') are missing in the JSON file.
json.JSONDecodeError
If the file is not a valid JSON.
RuntimeError
If the phase volume fractions do not sum to 1.0 (severe data inconsistency).
Notes
-----
- If `grain_phase_dict` or orientation information is missing,
the function will fall back to single-phase mode and log warnings.
- The mesh is either reconstructed from the JSON (if provided)
or generated using `mesh_creator()`.
- Periodicity and unit information are retained from the model metadata.
Examples
--------
>>> ms = import_voxels('example_vox.json', path='./input')
>>> print(ms.Ngr, 'grains loaded with', ms.nphases, 'phases.')
>>> print(ms.mesh.nodes.shape, 'mesh nodes')
"""
import json
import copy
from .api import Microstructure
from .initializations import RVE_creator, mesh_creator
from .entities import Simulation_Box
if file is None:
raise ValueError('Name for voxel file must be given.')
fname = os.path.join(path, file)
data = json.load(open(fname))
# extract basic model information
sh = tuple(data['Data']['Shape'])
nvox = np.prod(sh)
size = tuple(data['Model']['Size'])
gr_numbers = np.unique(data['Data']['Values'])
grain_keys = np.asarray(gr_numbers, dtype=str)
grains = np.reshape(data['Data']['Values'], sh, order=data['Data']['Order'])
ph_names = data['Model']['Phase_names']
nphases = len(ph_names)
phases = np.zeros(nvox, dtype=int)
grain_dict = dict()
grain_phase_dict = dict()
gr_arr = grains.flatten(order='C')
if 'Grains' in data.keys():
grain_ori_dict = {}
phase_vf = np.zeros(nphases)
ngrain = np.zeros(nphases, dtype=int)
for igr in gr_numbers:
#if igr == 0:
# continue # Grain 0 is reserved for matrix phase / unidentified phases
ind = np.nonzero(gr_arr == igr)[0]
nv = len(ind)
ip = data['Grains'][str(igr)]['Phase']
if not isinstance(ip, int) or not 0 <= ip < nphases:
raise ValueError(f'Invalid phase ID {ip!r} for grain {igr}.')
if igr == 0 and ip != 0:
raise ValueError('GRAIN0 must belong to matrix PHASE0.')
phase_vf[ip] += nv
grain_dict[int(igr)] = ind + 1
grain_phase_dict[int(igr)] = ip
ngrain[ip] += 1
phases[ind] = ip
orientation = data['Grains'][str(igr)].get('Orientation')
if igr != 0 and orientation is not None:
grain_ori_dict[int(igr)] = orientation
grain_ori_dict = grain_ori_dict or None
phase_vf /= nvox
if not np.isclose(np.sum(phase_vf), 1.):
logging.warning(f'Volume fractions do not add up to 1: {phase_vf}')
else:
# no grain-level information in data
grain_ori_dict = None
if nphases > 1:
logging.error('No grain-level information in data file.' +
'Cannot extract phase information or orientations.' +
'Continuing with single phase model.')
nphases = 1
for igr in gr_numbers:
ind = np.nonzero(gr_arr == igr)[0]
grain_dict[int(igr)] = ind + 1
grain_phase_dict[int(igr)] = 0
ngrain = [len(grain_keys)]
phase_vf = [1.]
if 0 in grain_dict and any(gid != 0 and pid == 0 for gid, pid in grain_phase_dict.items()):
raise ValueError('Matrix PHASE0 may only contain GRAIN0.')
# reconstructing microstructure information for RVE
stats_dict = {
'RVE': {'sideX': size[0], 'sideY': size[1], 'sideZ': size[2],
'Nx': sh[0], 'Ny': sh[1], 'Nz': sh[2]},
'Simulation': {'periodicity': data['Model']['Periodicity'],
'output_units': data['Model']['Units']['Length']},
'Phase': {'Name': 'Simulanium', 'Volume fraction': 1.0}
}
# add phase information and construct list of stats_dict's
stats_list = []
for i in range(nphases):
stats_dict['Phase']['Name'] = ph_names[i]
stats_dict['Phase']['Volume fraction'] = phase_vf[i]
phase_stats = copy.deepcopy(stats_dict)
if 0 in grain_dict and i == 0:
phase_stats['Grain type'] = 'Matrix'
stats_list.append(phase_stats)
# Create microstructure object
ms = Microstructure('from_voxels')
ms.name = data['Model']['Material']
ms.Ngr = np.sum(ngrain)
ms.nphases = nphases
ms.descriptor = stats_list
ms.ngrains = ngrain
ms.rve = RVE_creator(stats_list, from_voxels=True)
ms.simbox = Simulation_Box(size)
# initialize voxel structure (= mesh)
ms.mesh = mesh_creator(sh)
ms.mesh.grains = grains
ms.mesh.grain_dict = grain_dict
ms.mesh.grain_ori_dict = grain_ori_dict
ms.mesh.phases = phases.reshape(sh, order='C')
ms.mesh.grain_phase_dict = grain_phase_dict
ms.mesh.ngrains_phase = ngrain
ms.mesh.nphases = nphases
ms.vf_vox = np.asarray(phase_vf)
if 0 in ms.mesh.grain_dict.keys():
porosity = len(ms.mesh.grain_dict[0]) / nvox
ms.precipit = 1.0 - porosity
ms.mesh.prec_vf_voxels = ms.precipit
# import or create mesh
voxel_dict = dict()
vox_centerDict = dict()
if 'Mesh' in data.keys():
nodes = np.array(data['Mesh']['Nodes']['Values'])
for i, el in enumerate(data['Mesh']['Voxels']['Values']):
voxel_dict[i + 1] = el
ind = np.array(el, dtype=int) - 1
vox_centerDict[i + 1] = np.mean(nodes[ind], axis=0)
ms.mesh.voxel_dict = voxel_dict
ms.mesh.vox_center_dict = vox_centerDict
ms.mesh.nodes = nodes
else:
ms.mesh.create_voxels(ms.simbox)
print('\n Voxel structure imported.\n')
return ms
[docs]
def write_stats(
stats: Mapping[str, Any],
file: Union[str, os.PathLike[str]],
path: Union[str, os.PathLike[str]] = './') -> None:
"""
Write microstructure statistical descriptors to a JSON file.
This function serializes the provided statistical data (e.g., phase
fractions, grain sizes, RVE dimensions) and writes it to disk in JSON format.
The `.json` extension is automatically appended if not provided.
Parameters
----------
stats : dict
Dictionary containing statistical descriptors of the microstructure
(e.g., grain count, phase fractions, RVE size, or any other metadata).
file : str
Output filename for the JSON file. The '.json' extension will be appended
automatically if missing.
path : str, optional
Directory path where the JSON file will be written (default: './').
Raises
------
ValueError
If `stats` or `file` is None.
OSError
If the file cannot be written (e.g., due to permission or path issues).
TypeError
If `stats` contains non-serializable objects (invalid for JSON encoding).
Notes
-----
- The function uses Python's built-in :mod:`json` module.
- The output file is overwritten if it already exists.
Examples
--------
>>> stats = {"RVE": {"Nx": 20, "Ny": 20, "Nz": 20}, "Phase": {"Volume fraction": 0.8}}
>>> write_stats(stats, "example_stats")
# Writes './example_stats.json'
"""
import json
if stats is None:
raise ValueError('List or dict with microstructure descriptors must be given.')
if file is None:
raise ValueError('Name for json file with microstructure descriptors must be given.')
if file[-5:].lower() != '.json':
file += '.json'
file = os.path.join(path, file)
with open(file, 'w') as fp:
json.dump(stats, fp)
[docs]
def import_stats(
file: Union[str, os.PathLike[str]],
path: Union[str, os.PathLike[str]] = './') -> Any:
"""
Read microstructure statistical descriptors from a JSON file.
This function loads and parses a JSON file containing statistical
information about a microstructure, such as RVE dimensions, phase fractions,
or grain statistics.
Parameters
----------
file : str
Filename of the JSON file to read. The '.json' extension will be appended
automatically if missing.
path : str, optional
Directory path where the JSON file is located (default: './').
Returns
-------
desc : dict
Dictionary containing the statistical microstructure descriptors.
Raises
------
ValueError
If `file` is None.
FileNotFoundError
If the specified JSON file does not exist.
json.JSONDecodeError
If the file cannot be parsed as valid JSON.
OSError
If the file cannot be read due to permission or I/O errors.
Notes
-----
- The function expects the JSON file to contain key-value pairs describing
microstructure statistics compatible with the output of :func:`write_stats`.
Examples
--------
>>> desc = import_stats("example_stats.json")
>>> print(desc["RVE"]["Nx"])
20
"""
import json
if file is None:
raise ValueError('Name for json file with microstructure descriptors must be given.')
file = os.path.join(path, file)
with open(file, 'r') as inp:
desc = json.load(inp)
return desc