""" Module defining class Microstructure that contains the necessary
methods and attributes to analyze experimental microstructures in form
of EBSD maps to generate statistical descriptors for 3D microstructures, and
to create synthetic RVE that fulfill the required statistical microstructure
descriptors.
The methods of the class Microstructure for an API that can be used to generate
Python workflows.
Authors: Alexander Hartmaier, Golsa Tolooei Eshlghi, Abhishek Biswas
Institution: ICAMS, Ruhr University Bochum
"""
import os
import json
import logging
import platform
import hashlib
import orix
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Delaunay
from typing import Dict, Any, List, Optional, Union, Mapping
from importlib.metadata import version as pkg_version
from datetime import datetime
from .units import normalize_length_unit, length_scale_from_um
from .grains import calc_polygons
from .entities import Simulation_Box
from .input_output import export2abaqus, writeAbaqusMat, read_dump, _abaqus_phase_options
from .initializations import RVE_creator, mesh_creator, normalize_phase_descriptors, validate_matrix_mapping
from .packing import packingRoutine
from .voxelization import voxelizationRoutine
from .smoothingGB import smoothingRoutine
from .rve_stats import get_stats_vox, get_stats_part, get_stats_poly
from .plotting import plot_init_stats, plot_voxels_3D, plot_ellipsoids_3D, \
plot_polygons_3D, plot_output_stats, plot_particles_3D
[docs]
class Microstructure(object):
"""
Define a class for creating and managing synthetic microstructures
This class provides tools to define, generate, and analyze synthetic
microstructures composed of one or multiple phases. It integrates particle
packing, voxelization, and grain geometry generation for RVE (Representative
Volume Element) modeling.
Parameters
----------
descriptor : dict, list of dict, str, or None, optional
List of dictionaries describing the microstructure of each phase.
Dictionary keys typically include "Grains type", "Equivalent diameter",
"Aspect ratio", "Tilt Angle", "RVE", and "Simulation".
file : str or os.PathLike or None, optional
Path to a JSON or other input file to initialize the microstructure.
name : str, default='Microstructure'
Name of the microstructure
Attributes
----------
name : str
Name of the microstructure
nphases : int or None
Number of phases in the microstructure
ngrains : ndarray or None
Array of grain counts in each phase
nparticles : list or None
List of the number of particles in each phase
descriptor : list of dict or None
Descriptor used to define microstructure if provided
precipit : float or None
Indicates the presence of precipitates, pores, or secondary particles
in a continuous matrix. If float, specifies their volume fraction
from_voxels : bool
True if the microstructure is imported from a voxel file
particles : list or None
List of particle objects containing geometric information
rve : RVE_creator or None
Object containing RVE information including mesh dimensions,
particle count, periodicity, and phase volume fractions
simbox : Simulation_Box or None
Object defining the geometric boundaries of the RVE simulation domain
mesh : mesh_creator or None
Object storing voxelized mesh data including grain assignments,
voxel connectivity, and smoothed node coordinates
geometry : dict or None
Dictionary of grain geometries. Keys may include "Vertices", "Points",
"Simplices", "Facets", "Grains", "GBnodes", "GBarea", "GBfaces"
rve_stats : list of dict or None
Statistical information for different RVE representations,
including particles, voxels, and polyhedral grains
rve_stats_labels : list of str or None
Labels corresponding to the types of RVEs analyzed, e.g., "Particles",
"Voxels", "Grains"
vf_vox : ndarray or None
Phase volume fractions obtained from voxelized structures
ialloy : int or None
Alloy identifier used for material library mapping
Notes
-----
- The class can initialize from either a statistical descriptor or a JSON data file
- When only one phase is specified with a volume fraction < 1.0,
a matrix phase is automatically added with a complementary fraction (1 - vf)
- Kanapy is tested for up to two phases; using more may yield unpredictable results
"""
def __init__(
self,
descriptor: Optional[Union[Mapping[str, Any], List[Mapping[str, Any]], str]] = None,
file: Optional[Union[str, os.PathLike[str]]] = None,
name: str = 'Microstructure') -> None:
self.name = name
self.nphases = None
self.ngrains = None
self.nparticles = None
self.precipit = None
self.rve = None
self.particles = None
self.geometry = None
self.simbox = None
self.mesh = None
self.rve_stats = None
self.rve_stats_labels = None
self.from_voxels = False
self.ialloy = None
self.vf_vox = None
if descriptor is None:
if file is None:
raise ValueError('Please provide either a dictionary with statistics or an data file name')
# Open the user data statistics file and read the data
try:
with open(os.path.normpath(file)) as json_file:
self.descriptor = json.load(json_file)
except Exception as e:
logging.error(f'An unexpected exception occurred: {e}')
raise FileNotFoundError("File: '{}' does not exist in the current working directory!\n".format(file))
elif descriptor == 'from_voxels':
self.from_voxels = True
return
else:
self.descriptor = descriptor
if file is not None:
logging.warning('Input descriptor and file are given. Only descriptor will be used.')
self.descriptor = normalize_phase_descriptors(self.descriptor)
self.nphases = len(self.descriptor)
if any(d.get('Grain type') == 'Matrix' for d in self.descriptor):
self.precipit = 1.0 - self.descriptor[0]['Phase']['Volume fraction']
return
"""
-------- Routines for user interface --------
"""
[docs]
def init_RVE(
self,
descriptor: Optional[Union[Mapping[str, Any], List[Mapping[str, Any]]]] = None,
nsteps: int = 1000) -> None:
"""
Initialize the Representative Volume Element (RVE) of the microstructure
This method creates the voxel- or particle-based RVE using the provided
phase descriptor(s). It sets up the mesh dimensions, particle distribution,
and simulation box geometry. The RVE object is constructed using the
`RVE_creator` class, and key attributes such as `nparticles` and `simbox`
are updated accordingly.
Parameters
----------
descriptor : dict, list of dict, or None, optional, default=None
Description of the microstructure phases. Each dictionary specifies
parameters such as grain type, equivalent diameter, aspect ratio,
and tilt angle. If `None`, the class attribute `self.descriptor` is used.
nsteps : int, optional, default=1000
Number of optimization or relaxation steps for RVE generation.
Default is 1000.
Notes
-----
- Assumes that `RVE_creator` and `Simulation_Box` are available and correctly configured.
- A ``Grain type: Matrix`` descriptor reserves phase 0 for GRAIN0.
A single grain phase with fraction below one adds an implicit matrix.
Attributes Updated
------------------
rve : RVE_creator
The created RVE object containing voxel/particle information, size,
periodicity, and phase volume fractions.
nparticles : list
Number of particles per phase after RVE generation.
simbox : Simulation_Box
Object containing the geometric boundaries of the RVE domain.
precipit : float or None
Total fraction of grain-bearing phases when a matrix is present.
The matrix fraction is ``1.0 - precipit``.
Returns
-------
None
Updates class attributes with the initialized RVE and simulation box.
"""
if descriptor is None:
descriptor = self.descriptor
descriptor = normalize_phase_descriptors(descriptor, validate=True)
self.rve = RVE_creator(descriptor, nsteps=nsteps)
self.descriptor = descriptor
self.nphases = len(descriptor)
self.precipit = (1.0 - self.rve.phase_vf[0]
if self.rve.matrix_phase is not None else None)
self.nparticles = self.rve.nparticles
# store geometry in simbox object
self.simbox = Simulation_Box(self.rve.size)
[docs]
def pack(
self,
particle_data: Any = None,
k_rep: float = 0.0,
k_att: float = 0.0,
fill_factor: Optional[float] = None,
poly: Any = None,
save_files: bool = False,
verbose: bool = False) -> None:
"""
Pack particles into the simulation box according to the RVE settings.
Parameters
----------
particle_data : array-like or None, optional, default=None
Particle information. If None, uses ``self.rve.particle_data``.
k_rep : float, optional, default=0.0
Repulsion coefficient between particles. Default is 0.0.
k_att : float, optional, default=0.0
Attraction coefficient between particles. Default is 0.0.
fill_factor : float or None, optional, default=None
Fraction of simulation box to fill. If None and ``self.precipit`` is
set, the effective value is 1.0.
poly : object or None, optional, default=None
Additional packing options for polyhedral particles.
save_files : bool, optional, default=False
If True, saves packed particle data to files.
verbose : bool, optional, default=False
If True, prints progress and warnings during packing.
Returns
-------
None
Updates ``self.particles`` and ``self.simbox`` in place.
Raises
------
ValueError
If no particle data is available.
AttributeError
If RVE or simulation box is not initialized before packing.
Notes
-----
If self.precipit > 0.65, particle overlaps may occur.
Requires self.rve and self.simbox to be initialized before calling.
"""
if particle_data is None:
particle_data = self.rve.particle_data
if particle_data is None:
raise ValueError('No particle_data in pack. Run create_RVE first.')
if fill_factor is None and self.precipit is not None:
fill_factor = 1.0 # pack to full volume fraction defined in particles
print(f'Sparse particles (precipitates/pores): '
f'Packing up to particle volume fraction of {(100 * self.precipit):.1f}%.')
if self.precipit > 0.65:
print('Overlap of particles will occur since volume fraction > 65%')
self.particles, self.simbox = \
packingRoutine(particle_data, self.rve.periodic,
self.rve.packing_steps, self.simbox,
k_rep=k_rep, k_att=k_att, fill_factor=fill_factor,
poly=poly, save_files=save_files, verbose=verbose)
[docs]
def voxelize(
self,
particles: Optional[List[Any]] = None,
dim: Optional[tuple[int, int, int]] = None) -> None:
"""
Generate the RVE by assigning voxels to grains.
Parameters
----------
particles : list of object or None, optional, default=None
List of particle objects to voxelize. If None, uses ``self.particles``.
dim : tuple of int or None, optional, default=None
3-tuple specifying the number of voxels in each spatial direction.
If None, uses `self.rve.dim`.
Returns
-------
None
The voxelized mesh is stored in `self.mesh`. Phase voxel fractions are
stored in `self.vf_vox` and printed to the console.
Raises
------
ValueError
- If `particles` is None and `self.particles` is not available.
- If `dim` is not a 3-tuple when specified.
AttributeError
If `self.rve` or `self.simbox` is not initialized before voxelization.
Notes
-----
- Requires `self.particles` to be initialized by `pack()`.
- Updates `self.mesh`, `self.ngrains`, `self.Ngr`, and `self.vf_vox`.
- Logs warnings if phase volume fractions do not sum to 1.
- Removes polyhedral grain geometries (`self.geometry`) after re-meshing to avoid inconsistencies.
"""
if particles is None:
particles = self.particles
if particles is None:
raise ValueError('No particles in voxelize. Run pack first.')
if dim is None:
dim = self.rve.dim
else:
if len(dim) != 3 or type(dim) is not tuple:
raise ValueError(f'"dim" must be a 3-tuple of the voxel numbers in each direction, not {dim}.')
self.rve.dim = dim
# initialize voxel structure (= mesh)
self.mesh = mesh_creator(dim)
self.mesh.nphases = self.nphases
self.mesh.create_voxels(self.simbox)
self.mesh = \
voxelizationRoutine(particles, self.mesh, self.nphases, prec_vf=self.precipit)
if np.any(self.nparticles != self.mesh.ngrains_phase):
logging.info(f'Number of grains per phase changed from {self.nparticles} to ' +
f'{list(self.mesh.ngrains_phase)} during voxelization.')
self.ngrains = self.mesh.ngrains_phase
self.Ngr = np.sum(self.mesh.ngrains_phase, dtype=int)
# extract volume fractions from voxelized grains
if self.nphases > 1:
self.vf_vox = np.zeros(self.nphases)
vox_count = np.zeros(self.nphases, dtype=int)
for igr, ip in self.mesh.grain_phase_dict.items():
vox_count[ip] += len(self.mesh.grain_dict[igr])
print('Volume fractions of phases in voxel structure:')
vt = 0.
for ip in range(self.nphases):
vf_act = vox_count[ip] / self.mesh.nvox
self.vf_vox[ip] = vf_act
vt += vf_act
print(f'{ip}: {self.rve.phase_names[ip]} ({(vf_act * 100):.3f}%)')
if not np.isclose(vt, 1.0):
logging.warning(f'Volume fractions of phases in voxels do not add up to 1. Value: {vt}')
else:
self.vf_vox = np.ones(1)
# remove grain information if it already exists to avoid inconsistencies
if self.geometry is not None:
logging.info('Removing polyhedral grain geometries and statistical data after re-meshing.')
self.geometry = None
[docs]
def smoothen(
self,
nodes_v: Any = None,
voxel_dict: Optional[Mapping[Any, Any]] = None,
grain_dict: Optional[Mapping[Any, Any]] = None) -> None:
"""
Generate smoothed grain boundaries from a voxelated mesh.
Parameters
----------
nodes_v : array-like or None, optional, default=None
Mesh node coordinates. If None, uses `self.mesh.nodes`.
voxel_dict : dict or None, optional, default=None
Dictionary of voxels in the mesh. If None, uses `self.mesh.voxel_dict`.
grain_dict : dict or None, optional, default=None
Dictionary mapping grains to their voxels. If None, uses `self.mesh.grain_dict`.
Returns
-------
None
The smoothed nodes are stored in `self.mesh.nodes_smooth`. Grain boundary
faces are optionally stored in `self.geometry['GBfaces']`.
Raises
------
ValueError
If `nodes_v` is None and `self.mesh.nodes` is not available.
AttributeError
If `self.mesh` is not initialized before calling `smoothen()`.
Notes
-----
- Requires `self.mesh` to be initialized by `voxelize()`.
- Updates `self.mesh.nodes_smooth` with smoothed node coordinates.
- Adds 'GBfaces' entry to `self.geometry` if it exists.
"""
if nodes_v is None:
nodes_v = self.mesh.nodes
if nodes_v is None:
raise ValueError('No nodes_v in smoothen. Run voxelize first.')
if voxel_dict is None:
voxel_dict = self.mesh.voxel_dict
if grain_dict is None:
grain_dict = self.mesh.grain_dict
self.mesh.nodes_smooth, grain_facesDict = \
smoothingRoutine(nodes_v, voxel_dict, grain_dict)
if isinstance(self.geometry, dict):
self.geometry['GBfaces'] = grain_facesDict
[docs]
def generate_grains(self) -> None:
"""
Calculate and store polyhedral grain geometry, including particle- and grain-diameter attributes,
for statistical comparison
Notes
-----
- Requires `self.mesh` to be initialized by `voxelize`.
- Updates `self.geometry` with polyhedral grain volumes and shared grain boundary (GB) areas.
- If `self.precipit` is True, irregular grain 0 is temporarily removed from analysis.
- Logs warnings if grains are not represented in the polyhedral geometry.
- Prints volume fractions of each phase in the polyhedral geometry.
Parameters
----------
None
Returns
-------
None
The calculated grain geometry is stored in `self.geometry`. Phase volume fractions
are printed to the console.
Raises
------
ValueError
- If `self.mesh` or `self.mesh.grains` is None (i.e., voxelized microstructure not available).
"""
if self.mesh is None or self.mesh.grains is None:
raise ValueError('No information about voxelized microstructure. Run voxelize first.')
if self.precipit and 0 in self.mesh.grain_dict.keys():
# in case of precipit, remove irregular grain 0 from analysis
empty_vox = self.mesh.grain_dict.pop(0)
grain_store = self.mesh.grain_phase_dict.pop(0)
else:
empty_vox = None
grain_store = None
try:
self.geometry = calc_polygons(self.rve, self.mesh)
finally:
if empty_vox is not None:
self.mesh.grain_dict[0] = empty_vox
self.mesh.grain_phase_dict[0] = grain_store
# verify that geometry['Grains'] and mesh.grain_dict are consistent
"""if np.any(self.geometry['Ngrains'] != self.ngrains):
logging.warning(f'Only facets for {self.geometry["Ngrains"]} created, but {self.Ngr} grains in voxels.')
for igr in self.mesh.grain_dict.keys():
if igr not in self.geometry['Grains'].keys():
logging.warning(f'Grain: {igr} not in geometry. Be aware when creating GB textures.')"""
# verify that geometry['GBarea'] is consistent with geometry['Grains']
gba = self.geometry['GBarea']
ind = []
igr = []
for i, gblist in enumerate(gba):
if not gblist[0] in self.geometry['Grains'].keys():
ind.append(i)
igr.append(gblist[0])
continue
if not gblist[1] in self.geometry['Grains'].keys():
ind.append(i)
igr.append(gblist[1])
if len(ind) > 0:
logging.warning(f'{len(ind)} grains are not represented in polyhedral geometry.')
# logging.warning('Consider increasing the number of voxels, as grains appear to be very irregular.')
"""ind.reverse()
igr.reverse()
for j, i in enumerate(ind):
logging.warning(f'Removing {gba[i]} from GBarea as grain {igr[j]} does not exist.')
gba.pop(i)
self.geometry['GBarea'] = gba"""
# extract volume fractions from polyhedral grains
if self.nphases > 1:
ph_vol = np.zeros(self.nphases)
for igr, grd in self.geometry['Grains'].items():
ip = grd['Phase']
ph_vol[ip] += grd['Volume']
print('Volume fractions of phases in polyhedral geometry:')
for ip in range(self.nphases):
if self.precipit is not None and ip == 0:
continue # Matrix GRAIN0 has no reconstructed grain geometry.
vf = 100.0 * ph_vol[ip] / np.prod(self.rve.size)
print(f'{ip}: {self.rve.phase_names[ip]} ({vf.round(1)}%)')
[docs]
def generate_orientations(
self,
data: Any,
ang: Optional[float] = None,
omega: Optional[float] = None,
Nbase: int = 5000,
hist: Any = None,
shared_area: Any = None,
iphase: Optional[int] = None,
verbose: bool = False,
**kwargs: Any) -> None:
"""
Generate orientations for grains in a representative volume element (RVE)
to achieve a desired crystallographic texture
This method assigns grain orientations to achieve a specified texture. The
input can be an EBSDmap object or a string defining the type of orientation set
(random or unimodal). For unimodal textures, `ang` and `omega` must be specified.
Generated orientations are stored in `self.mesh.grain_ori_dict`, and
`self.mesh.texture` is updated accordingly.
Parameters
----------
data : EBSDmap or str
Input source for orientations. Can be an `EBSDmap` object or a string
specifying the type of orientation set:
- 'random' or 'rnd' : generate a random orientation set
- 'unimodal', 'uni_mod', or 'uni_modal' : generate a unimodal orientation set
ang : float or None, optional, default=None
Orientation angle for unimodal texture (required if `data` is unimodal).
omega : float or None, optional, default=None
Kernel halfwidth for unimodal texture (required if `data` is unimodal).
Nbase : int, default=5000
Number of base orientations used in random or unimodal generation.
hist : array_like or None, optional, default=None
Histogram for the grain orientations, used to weight orientations.
shared_area : array_like, float, or None, optional, default=None
Shared grain boundary area for weighted orientation generation.
iphase : int or None, optional, default=None
Phase index to update, preserving orientations in other phases. If None,
all grain-bearing phases are processed. Matrix GRAIN0 is always excluded.
verbose : bool, default=False
If True, prints additional information during orientation generation.
ebsd_phase_map : dict, optional
Passed through kwargs; maps canonical Kanapy phase IDs to EBSD phase
indices. By default, a matrix occupies phase 0 and grain phase i
uses EBSD phase i-1; without a matrix, indices are unchanged.
Returns
-------
None
The generated orientations are stored in `self.mesh.grain_ori_dict` and
`self.mesh.texture` is updated to reflect the type of generated texture.
Raises
------
ValueError
- If grain geometry is not defined (i.e., `self.mesh.grains` is None).
- If `data` is unimodal but `ang` or `omega` are not provided.
- If `data` is neither `EBSDmap` nor a recognized string option.
- If histogram is provided but GB areas are not defined (and `shared_area` is None).
Examples
--------
>>> # Generate random orientations for all grains
>>> rve.generate_orientations('random')
>>> # Generate unimodal orientations with specified angle and halfwidth
>>> rve.generate_orientations('unimodal', ang=30.0, omega=10.0)
>>> # Generate orientations for a specific phase with histogram weighting
>>> rve.generate_orientations('random', hist=hist_array, iphase=1)
>>> # Generate orientations using an EBSDmap object
>>> rve.generate_orientations(ebsd_map_obj)
"""
from kanapy import __backend__
if __backend__ == 'mtex':
from kanapy_mtex.texture import EBSDmap, createOrisetRandom, createOriset
logging.info('Using MTEX library to read EBSD maps and generate orientations.')
MTEX = True
else:
from kanapy.texture import EBSDmap, createOrisetRandom, createOriset
logging.info('Using ORIX library to read EBSD maps and generate orientations.')
MTEX = False
if self.mesh.grains is None:
raise ValueError('Grain geometry is not defined. Run voxelize first.')
if shared_area is None:
if hist is None:
gba = None
else:
if self.geometry is None:
raise ValueError('If histogram for GB texture is provided, GB areas must be defined.\n' +
'Run generate_grains() first, to calculate GB areas.')
gba = self.geometry['GBarea']
else:
if shared_area == 0:
gba = None
else:
gba = shared_area
if iphase is not None and not 0 <= iphase < len(self.ngrains):
raise ValueError('iphase must identify an existing phase.')
ori_dict = dict(self.mesh.grain_ori_dict or {}) if iphase is not None else {}
ori_dict.pop(0, None) # GRAIN0 is the unoriented matrix, never a CP grain.
for ip in range(len(self.ngrains)):
if iphase is not None and iphase != ip:
continue
phase_grains = [gid for gid in self.mesh.grain_dict
if gid != 0 and self.mesh.grain_phase_dict[gid] == ip]
ngr = len(phase_grains)
if ngr == 0:
continue
if isinstance(data, EBSDmap):
if iphase is None or iphase == ip:
if gba is not None and not MTEX:
logging.warning('Shared GB area (option: gba) is currently only available in kanapy-mtex.\n'
'This option will be ignored.')
gba = None
allowed = ["res_low", "res_high", "res_step", "lim", "hw_init"]
ori_kwargs = {k: v for k, v in kwargs.items() if k in allowed}
source_phase = kwargs.get("ebsd_phase_map", {}).get(
ip, ip - 1 if getattr(self, "precipit", None) is not None else ip)
ori_rve = data.calcORI(ngr, iphase=source_phase, shared_area=gba, verbose=verbose, **ori_kwargs)
self.mesh.texture = "ODF"
elif isinstance(data, str):
if data.lower() in ['random', 'rnd']:
ori_rve = createOrisetRandom(ngr, Nbase=Nbase, hist=hist, shared_area=gba)
self.mesh.texture = "Random"
elif data.lower() in ['unimodal', 'uni_mod', 'uni_modal']:
if ang is None or omega is None:
raise ValueError('To generate orientation sets of type "unimodal" angle "ang" and kernel' +
'halfwidth "omega" are required.')
ori_rve = createOriset(ngr, ang, omega, hist=hist, shared_area=gba, verbose=verbose)
self.mesh.texture = "Unimodal"
else:
self.mesh.texture = None
raise ValueError('Argument to generate grain orientation must be either of type EBSDmap or ' +
'"random" or "unimodal"')
for ind, igr in enumerate(phase_grains):
ori_dict[igr] = ori_rve[ind, :]
self.mesh.grain_ori_dict = ori_dict
return
"""
-------- Plotting methods --------
"""
[docs]
def plot_ellipsoids(
self,
cmap: str = 'prism',
dual_phase: Optional[bool] = None,
phases: bool = False) -> None:
"""
Generate a 3D plot of particles in the RVE
This function visualizes ellipsoidal particles in the Representative Volume Element (RVE).
Particles can be colored according to phase, and a custom colormap can be used.
Note that `dual_phase` is deprecated; use `phases` instead.
Parameters
----------
cmap : str, optional, default='prism'
Colormap used for plotting particles. Default is 'prism'.
dual_phase : bool or None, optional, default=None
Deprecated parameter for indicating dual-phase visualization.
Use `phases` instead. Default is None.
phases : bool, optional, default=False
If True, color particles according to their phase. Default is False.
Notes
-----
- Requires `self.particles` to be initialized by `pack`.
- Automatically calculates aspect ratios of the RVE for proper 3D plotting.
- Prints a warning if `dual_phase` is used and maps it to `phases`.
Returns
-------
None
Displays a 3D plot of ellipsoidal particles using `plot_ellipsoids_3D`.
Raises
------
ValueError
- If `self.particles` is None.
Examples
--------
>>> # Simple 3D plot with default colormap
>>> rve.plot_ellipsoids()
>>> # Color particles according to phase
>>> rve.plot_ellipsoids(phases=True)
>>> # Use a custom colormap
>>> rve.plot_ellipsoids(cmap='viridis')
"""
if dual_phase is not None:
print('Use of "dual_phase" is depracted. Use parameter "phases" instead.')
phases = dual_phase
if self.particles is None:
raise ValueError('No particle to plot. Run pack first.')
hmin = min(self.rve.size)
asp_arr = [int(self.rve.size[0] / hmin),
int(self.rve.size[1] / hmin),
int(self.rve.size[2] / hmin)]
plot_ellipsoids_3D(self.particles, cmap=cmap, phases=phases, asp_arr=asp_arr)
[docs]
def plot_particles(
self,
cmap: str = 'prism',
dual_phase: Optional[bool] = None,
phases: bool = False,
plot_hull: bool = True) -> None:
"""
Generate a 3D plot of particles in the RVE.
This function visualizes particles in the Representative Volume Element (RVE) in 3D.
Particles can be colored according to phase, plotted with their convex hulls, and
displayed using a custom colormap.
Parameters
----------
cmap : str, optional, default='prism'
Colormap used for plotting particles. Default is 'prism'.
dual_phase : bool or None, optional, default=None
Deprecated parameter for indicating dual-phase visualization.
Use `phases` instead. Default is None.
phases : bool, optional, default=False
If True, color particles according to their phase. Default is False.
plot_hull : bool, optional, default=True
If True, plot the convex hull (inner polygon) of each particle. Default is True.
Notes
-----
- Requires `self.particles` to be initialized by `pack`.
- Automatically calculates aspect ratios of the RVE for proper 3D plotting.
- Prints a warning if `dual_phase` is used and maps it to `phases`.
Returns
-------
None
Displays a 3D plot of particles using `plot_particles_3D`.
Raises
------
ValueError
- If `self.particles` is None.
- If particles do not have inner polygons.
Examples
--------
>>> # Simple 3D plot with default colormap and convex hulls
>>> rve.plot_particles()
>>> # Color particles according to phase
>>> rve.plot_particles(phases=True)
>>> # Plot particles without convex hulls and use a custom colormap
>>> rve.plot_particles(cmap='viridis', plot_hull=False)
"""
if dual_phase is not None:
print('Use of "dual_phase" is depracted. Use parameter "phases" instead.')
phases = dual_phase
if self.particles is None:
raise ValueError('No particle to plot. Run pack first.')
if self.particles[0].inner is None:
raise ValueError('Ellipsoids without inner polygon cannot be plotted.')
hmin = min(self.rve.size)
asp_arr = [int(self.rve.size[0] / hmin),
int(self.rve.size[1] / hmin),
int(self.rve.size[2] / hmin)]
plot_particles_3D(self.particles, cmap=cmap,
phases=phases, plot_hull=plot_hull, asp_arr=asp_arr)
[docs]
def plot_voxels(
self,
sliced: bool = False,
dual_phase: Optional[bool] = None,
phases: bool = False,
cmap: str = 'prism',
ori: Any = None,
color_key: int = 0,
silent: bool = False) -> Optional[Any]:
"""
Generate a 3D visualization of the voxelized RVE structure
This function visualizes the voxel-based microstructure of the RVE. It supports
coloring by phase, grain ID, or crystallographic orientation, and can optionally
render a sliced view of the 3D voxel mesh.
Parameters
----------
sliced : bool, optional, default=False
If True, generates a sliced view of the voxel mesh to visualize the internal structure.
Default is False.
dual_phase : bool or None, optional, default=None
Deprecated parameter for dual-phase visualization. Use `phases` instead.
Default is None.
phases : bool, optional, default=False
If True, color voxels by phase instead of grain ID. Default is False.
cmap : str, optional, default='prism'
Name of the matplotlib colormap used for rendering. Default is 'prism'.
ori : array-like, bool, or None, optional, default=None
Array of grain orientations, or True to use `self.mesh.grain_ori_dict` for coloring
via inverse pole figure (IPF) mapping. Default is None.
color_key : int, optional, default=0
Selects the color mapping for orientations:
- 0: iphHSVKey
- 1: BungeColorKey
- 2: ipfHKLKey
Default is 0.
silent : bool, optional, default=False
If True, suppresses figure display and returns the matplotlib figure object instead.
Default is False.
Returns
-------
fig : matplotlib.figure.Figure or None
Returns the figure object if `silent=True`. Otherwise, displays the 3D plot and returns None.
Raises
------
ValueError
If `self.mesh.grains` is None, indicating that voxelization has not been performed.
Notes
-----
- Requires that the voxel mesh (`self.mesh`) has been generated by `voxelize()`.
- If `ori` is provided, orientation-based coloring is applied using `get_ipf_colors()`.
- Prints a warning if `dual_phase` is used, since it is deprecated and replaced by `phases`.
- Automatically computes aspect ratios of the RVE to ensure accurate 3D visualization.
Examples
--------
>>> # Simple 3D plot of voxelized RVE
>>> rve.plot_voxels()
>>> # Plot a sliced view of the voxel mesh
>>> rve.plot_voxels(sliced=True)
>>> # Color voxels by phase and use a custom colormap
>>> rve.plot_voxels(phases=True, cmap='viridis')
>>> # Use orientation-based coloring via IPF mapping
>>> rve.plot_voxels(ori=True, color_key=2)
>>> # Suppress figure display and get the matplotlib Figure object
>>> fig = rve.plot_voxels(silent=True)
"""
if dual_phase is not None:
print('Use of "dual_phase" is depracted. Use parameter "phases" instead.')
phases = dual_phase
if self.mesh.grains is None:
raise ValueError('No voxels or elements to plot. Run voxelize first.')
elif phases:
data = self.mesh.phases
else:
data = self.mesh.grains
if ori is not None:
from kanapy import __backend__
if __backend__ == "mtex":
from kanapy_mtex.texture import get_ipf_colors
else:
from kanapy.texture import get_ipf_colors
if isinstance(ori, bool) and ori:
ori = np.array([val for val in self.mesh.grain_ori_dict.values()])
clist = get_ipf_colors(ori, color_key)
else:
clist = None
hmin = min(self.rve.size)
asp_arr = [int(self.rve.size[0] / hmin),
int(self.rve.size[1] / hmin),
int(self.rve.size[2] / hmin)]
fig = plot_voxels_3D(data, sliced=sliced,
phases=phases, cmap=cmap, clist=clist,
silent=silent, asp_arr=asp_arr)
if silent:
return fig
[docs]
def plot_grains(
self,
geometry: Optional[Mapping[str, Any]] = None,
cmap: str = 'prism',
alpha: float = 0.4,
ec: Optional[List[float]] = None,
dual_phase: Optional[bool] = None,
phases: bool = False) -> None:
"""
Plot the polygonalized microstructure of the RVE in 3D.
This function visualizes the polygonal (polyhedral) grain geometry of the
microstructure, typically after `generate_grains()` has been executed.
Users can customize the colormap, transparency, and edge color. It also
supports coloring grains by phase.
Parameters
----------
geometry : dict or None, optional, default=None
Dictionary containing the polygonal grain geometries. If None, uses
`self.geometry`.
cmap : str, optional, default='prism'
Matplotlib colormap name for rendering grain colors. Default is 'prism'.
alpha : float, optional, default=0.4
Transparency level of the grain surfaces (0 = fully transparent,
1 = fully opaque). Default is 0.4.
ec : list of float or None, optional, default=None
Edge color specified as an RGBA list, e.g. `[0.5, 0.5, 0.5, 0.1]`.
Default is `[0.5, 0.5, 0.5, 0.1]`.
dual_phase : bool or None, optional, default=None
Deprecated parameter for dual-phase coloring. Use `phases` instead.
Default is None.
phases : bool, optional, default=False
If True, color grains by phase rather than by grain ID. Default is False.
Notes
-----
- Requires polygonal grain data generated by `generate_grains()`.
- Automatically computes the aspect ratio of the RVE for correct 3D scaling.
- Visualization is handled by `plot_polygons_3D()`.
- The parameter `dual_phase` is deprecated; prefer using `phases` instead.
Returns
-------
None
Displays a 3D plot of polygonal grains.
Raises
------
ValueError
If no polygonal geometry is available (i.e., `self.geometry` is None).
Examples
--------
>>> # Simple 3D plot of polygonal grains with default settings
>>> rve.plot_grains()
>>> # Color grains by phase
>>> rve.plot_grains(phases=True)
>>> # Use a custom colormap and set transparency
>>> rve.plot_grains(cmap='viridis', alpha=0.6)
>>> # Specify edge color
>>> rve.plot_grains(ec=[0.2, 0.2, 0.2, 0.3])
"""
if ec is None:
ec = [0.5, 0.5, 0.5, 0.1]
if dual_phase is not None:
print('Use of "dual_phase" is depracted. Use parameter "phases" instead.')
phases = dual_phase
if geometry is None:
geometry = self.geometry
if geometry is None:
raise ValueError('No polygons for grains defined. Run generate_grains() first')
hmin = min(self.rve.size)
asp_arr = [int(self.rve.size[0] / hmin),
int(self.rve.size[1] / hmin),
int(self.rve.size[2] / hmin)]
plot_polygons_3D(geometry, cmap=cmap, alpha=alpha, ec=ec,
phases=phases, asp_arr=asp_arr)
[docs]
def plot_stats(
self,
data: Optional[str] = None,
gs_data: Any = None,
gs_param: Any = None,
ar_data: Any = None,
ar_param: Any = None,
dual_phase: Optional[bool] = None,
phases: bool = False,
save_files: bool = False,
show_all: bool = False,
verbose: bool = False,
silent: bool = False,
enhanced_plot: bool = False) -> Optional[List[Any]]:
"""
Plot particle, voxel, and grain diameter statistics for comparison
This method analyzes the microstructure at different representation levels
(particles, voxels, and polyhedral grains) and plots corresponding statistical
distributions. It can optionally handle multiphase materials, save plots,
and display detailed information about geometric parameters.
Parameters
----------
data : str or None, optional, default=None
Specifies which type of data to analyze and plot:
- 'p' : particles
- 'v' : voxels
- 'g' : grains
If None, all available data types are analyzed.
gs_data : list or array-like or None, optional, default=None
Grain size data for comparison with simulation results.
gs_param : list or array-like or None, optional, default=None
Parameters for fitting grain size distributions.
ar_data : list or array-like or None, optional, default=None
Aspect ratio data for comparison with simulation results.
ar_param : list or array-like or None, optional, default=None
Parameters for fitting aspect ratio distributions.
dual_phase : bool or None, optional, default=None
Deprecated. Use `phases` instead.
phases : bool, optional, default=False
If True, perform separate statistical analysis for each phase.
save_files : bool, optional, default=False
If True, save generated plots and statistical results to files.
show_all : bool, optional, default=False
If True, display all generated plots interactively.
verbose : bool, optional, default=False
If True, print detailed numerical results during analysis.
silent : bool, optional, default=False
If True, suppresses console output and returns figures directly.
enhanced_plot : bool, optional, default=False
If True, use enhanced plot styling (automatically enabled when `silent=True`).
Returns
-------
flist : list of matplotlib.figure.Figure or None
List of generated figure objects if `silent=True`, otherwise None.
Raises
------
None
This function does not raise exceptions if data is missing. Instead, it logs an
error message and returns immediately.
Notes
-----
- Computes and compares geometric parameters such as principal axis lengths
(a, b, c), aspect ratios, rotation axes, and equivalent diameters.
- Results are stored in `self.rve_stats` and labeled in `self.rve_stats_labels`.
- For multiphase structures, statistics are calculated per phase.
- `dual_phase` is deprecated; prefer using `phases`.
Examples
--------
>>> # Plot all statistics with default settings
>>> rve.plot_stats()
>>> # Plot only particle statistics and show all plots interactively
>>> rve.plot_stats(data='p', show_all=True)
>>> # Perform multiphase statistical analysis and save figures
>>> rve.plot_stats(phases=True, save_files=True)
>>> # Suppress console output and get figure objects
>>> figs = rve.plot_stats(silent=True)
"""
if dual_phase is not None:
print('Use of "dual_phase" is depracted. Use parameter "phases" instead.')
phases = dual_phase
if silent:
verbose = False
show_all = False
enhanced_plot = True
ax_max = np.prod(self.rve.size) ** (1 / 3)
if phases:
phase_ids = [ip for ip in range(self.nphases)
if not (self.precipit is not None and ip == 0)]
if self.mesh is not None:
present = {pid for gid, pid in self.mesh.grain_phase_dict.items() if gid != 0}
phase_ids = [ip for ip in phase_ids if ip in present]
else:
phase_ids = [None]
nphases = len(phase_ids)
# Comparison lists may include a matrix slot; retain canonical phase IDs.
gs_data, gs_param, ar_data, ar_param = [
[values[ip] for ip in phase_ids]
if phases and isinstance(values, list) and len(values) == self.nphases else values
for values in (gs_data, gs_param, ar_data, ar_param)
]
if not (isinstance(gs_data, list) and len(gs_data) == nphases):
gs_data = [gs_data] * nphases
if not (isinstance(gs_param, list) and len(gs_param) == nphases):
gs_param = [gs_param] * nphases
if not (isinstance(ar_data, list) and len(ar_data) == nphases):
ar_data = [ar_data] * nphases
if not (isinstance(ar_param, list) and len(ar_param) == len(gs_param)):
ar_param = [ar_param] * nphases
"""
gs_data = [ebsd.ms_data[i]['gs_data'] for i in range(nphases)]
gs_param = [ebsd.ms_data[i]['gs_param'] for i in range(nphases)]
ar_data = [ebsd.ms_data[i]['ar_data'] for i in range(nphases)]
ar_param = [ebsd.ms_data[i]['ar_param'] for i in range(nphases)]
"""
iphase = None
flist = []
for slot, ip in enumerate(phase_ids):
stats_list = []
labels = []
if phases:
iphase = ip
print(f'Plotting statistical information for phase {ip}')
# Analyze and plot particles statistics
if (data is None and self.particles is not None) or \
(type(data) is str and 'p' in data.lower()):
if self.particles is None:
logging.error('Particle statistics requested, but no particles defined. '
'Run "pack()" first.')
return
part_stats = get_stats_part(self.particles, iphase=iphase, ax_max=ax_max,
show_plot=show_all,
verbose=verbose, save_files=save_files)
stats_list.append(part_stats)
labels.append("Partcls")
# Analyze and plot statistics of voxel structure in RVE
if (data is None and self.mesh is not None) or \
(type(data) is str and 'v' in data.lower()):
if self.mesh is None:
logging.error('Voxel statistics requested, but no voxel mesh defined. '
'Run "voxelize()" first.')
return
vox_stats = get_stats_vox(self.mesh, iphase=iphase, ax_max=ax_max,
show_plot=show_all,
verbose=verbose, save_files=save_files)
stats_list.append(vox_stats)
labels.append('Voxels')
# Analyze and plot statistics of polyhedral grains in RVE
if (data is None and self.geometry is not None) or \
(type(data) is str and 'g' in data.lower()):
if self.geometry is None:
logging.error('Geometry statistics requested, but no polyhedral grains defined. '
'Run "generate_grains()" first.')
return
grain_stats = get_stats_poly(self.geometry['Grains'], iphase=iphase, ax_max=ax_max,
show_plot=show_all, phase_dict=self.mesh.grain_phase_dict,
verbose=verbose, save_files=save_files)
stats_list.append(grain_stats)
labels.append('Grains')
if phases:
print(f'\nStatistical microstructure parameters of phase {iphase} in RVE')
print('-------------------------------------------------------')
else:
print('\nStatistical microstructure parameters of RVE')
print('--------------------------------------------')
print(f'Type\t| a (µm) \t| b (µm) \t| c (µm) \t| std.dev\t| rot.axis\t| asp.ratio\t| std.dev\t|'
f' equ.dia. (µm)\t| std.dev')
for i, sd in enumerate(stats_list):
av_std = np.mean([sd['a_sig'], sd['b_sig'], sd['c_sig']])
print(f'{labels[i]}\t| {sd["a_scale"]:.3f}\t| {sd["b_scale"]:.3f}\t| {sd["c_scale"]:.3f}\t| '
f'{av_std:.4f}\t| {sd["ind_rot"]} \t| {sd["ar_scale"]:.3f}\t| {sd["ar_sig"]:.4f}\t| '
f'{sd["eqd_scale"]:.3f} \t| {sd["eqd_sig"]:.4f}')
self.rve_stats = stats_list
self.rve_stats_labels = labels
fig = plot_output_stats(stats_list, labels, iphase=iphase,
gs_data=gs_data[slot], gs_param=gs_param[slot],
ar_data=ar_data[slot], ar_param=ar_param[slot],
save_files=save_files, silent=silent,
enhanced_plot=enhanced_plot)
flist.append(fig)
if silent:
return flist
[docs]
def plot_stats_init(
self,
descriptor: Optional[Union[Mapping[str, Any], List[Mapping[str, Any]]]] = None,
gs_data: Any = None,
ar_data: Any = None,
porous: bool = False,
get_res: bool = False,
show_res: bool = False,
save_files: bool = False,
silent: bool = False,
return_descriptors: bool = False) -> Optional[tuple[List[Any], List[Dict[str, Any]]]]:
"""
Plot initial statistical microstructure descriptors of RVE and optionally return computed descriptors
This method analyzes the initial microstructure defined by input descriptors
(or the class attribute `self.descriptor`) and plots the distributions of
equivalent diameters and aspect ratios. For elongated grains, it prints a
summary of input and output statistics. Optionally, statistical descriptors
can be returned for further processing.
Parameters
----------
descriptor : dict, list of dict, or None, optional, default=None
Microstructure phase descriptor(s). If None, uses `self.descriptor`.
gs_data : list or array-like or None, optional, default=None
Grain size data for comparison with initial statistics.
ar_data : list or array-like or None, optional, default=None
Aspect ratio data for comparison with initial statistics.
porous : bool, optional, default=False
If True, only the first phase is considered (e.g., for porous structures).
get_res : bool, optional, default=False
If True, computes statistical descriptors from the voxelized structure.
show_res : bool, optional, default=False
If True, prints detailed statistical results to the console.
save_files : bool, optional, default=False
If True, saves generated plots to files.
silent : bool, optional, default=False
If True, suppresses console output and returns figures directly.
return_descriptors : bool, optional, default=False
If True, returns computed statistical descriptors along with figures.
Returns
-------
tuple of (list, list) or None
Returns ``(figures, descriptors)`` when ``silent`` or
``return_descriptors`` is True; otherwise returns None.
Notes
-----
- Requires voxelized microstructure (`self.mesh`) for computing descriptors.
- Computes statistics for equivalent diameter, aspect ratio, principal axes,
and rotation axes when applicable.
- Can handle multiphase or single-phase microstructures.
Examples
--------
>>> # Plot initial statistics with default settings
>>> rve.plot_stats_init()
>>> # Plot only the first phase (porous) and save figures
>>> rve.plot_stats_init(porous=True, save_files=True)
>>> # Compute descriptors and get figures and descriptor objects
>>> figs, descs = rve.plot_stats_init(get_res=True, return_descriptors=True)
>>> # Suppress console output and return figure objects
>>> figs = rve.plot_stats_init(silent=True)
"""
def analyze_voxels(
ip: int,
des: Mapping[str, Any]) -> tuple[List[float], List[float], Dict[str, Any]]:
"""
Compute voxel-based statistical descriptors for a given phase of the RVE
This function calculates equivalent diameter, aspect ratio, principal axes, and rotation axis
statistics from the voxelized microstructure. It prints input vs output statistics if the phase
is elongated and returns formatted arrays for grain size, aspect ratio, and a dictionary of descriptors
Parameters
----------
ip : int
Index of the phase to analyze
des : dict
Descriptor dictionary for the phase containing grain type and target statistics
Returns
-------
gsp : list
Summary of equivalent diameter statistics for plotting
arp : list
Summary of aspect ratio statistics for plotting
statistical_descriptors : dict
Detailed statistics including mean, std, min, max of eqd, ar, axes, and rotation axis
"""
if self.mesh is None:
raise ValueError('show_res is True, but no voxels have been defined. Run voxelize first.')
vox_stats = get_stats_vox(self.mesh, iphase=ip, show_plot=show_res)
gsp = [vox_stats['eqd_sig'], 0.0, vox_stats['eqd_scale'],
min(vox_stats['eqd']), max(vox_stats['eqd'])]
arp = [vox_stats['ar_sig'], 0.0, vox_stats['ar_scale'],
min(vox_stats['ar']), max(vox_stats['ar'])]
if 'Grain type' in des.keys() and des['Grain type'] == 'Elongated':
if nel > 1:
print(f'\nStatistical microstructure parameters of phase {ip} in RVE')
print('-------------------------------------------------------')
else:
print('\nStatistical microstructure parameters of RVE')
print('--------------------------------------------')
print(f'Type\t| a (µm) \t| b (µm) \t| c (µm) \t| std.dev\t| rot.axis\t| asp.ratio\t| std.dev\t|'
f' equ.dia. (µm)\t| std.dev')
av_std = np.mean([vox_stats['a_sig'], vox_stats['b_sig'], vox_stats['c_sig']])
print(f'Input\t| - \t| - \t| - \t| - \t| - \t| '
f'{des["Aspect ratio"]["scale"]:.3f}\t| {des["Aspect ratio"]["sig"]:.4f}\t| '
f'{des["Equivalent diameter"]["scale"]:.3f} \t| {des["Equivalent diameter"]["sig"]:.4f}')
print(f'Output\t| {vox_stats["a_scale"]:.3f}\t| {vox_stats["b_scale"]:.3f}\t| '
f'{vox_stats["c_scale"]:.3f}\t| {av_std:.4f}\t| {vox_stats["ind_rot"]} \t| '
f'{vox_stats["ar_scale"]:.3f}\t| {vox_stats["ar_sig"]:.4f}\t| '
f'{vox_stats["eqd_scale"]:.3f} \t| {vox_stats["eqd_sig"]:.4f}')
statistical_descriptors = {
'eqd': {
'mean': float(vox_stats['eqd_scale']),
'std': float(vox_stats['eqd_sig']),
'min': float(min(vox_stats['eqd'])),
'max': float(max(vox_stats['eqd'])),
},
'ar': {
'mean': float(vox_stats['ar_scale']),
'std': float(vox_stats['ar_sig']),
'min': float(min(vox_stats['ar'])),
'max': float(max(vox_stats['ar'])),
},
'axes': {
'a': float(vox_stats['a_scale']),
'b': float(vox_stats['b_scale']),
'c': float(vox_stats['c_scale']),
'a_std': float(vox_stats.get('a_sig', np.nan)),
'b_std': float(vox_stats.get('b_sig', np.nan)),
'c_std': float(vox_stats.get('c_sig', np.nan)),
},
'rotation_axis': vox_stats.get('ind_rot', None),
}
return gsp, arp, statistical_descriptors
if show_res: get_res = True
if silent: show_res = False
if descriptor is None: descriptor = self.descriptor
descriptor = normalize_phase_descriptors(descriptor)
if porous: descriptor = descriptor[0:1]
nel = len(descriptor)
if not (isinstance(gs_data, list) and len(gs_data) == nel): gs_data = [gs_data] * nel
if not (isinstance(ar_data, list) and len(ar_data) == nel): ar_data = [ar_data] * nel
flist, descs = [] , []
for ip, des in enumerate(descriptor):
if des.get("Grain type") == "Matrix":
continue
gsp = arp = None
statistical_descriptors = None
if get_res:
gsp, arp, statistical_descriptors = analyze_voxels(ip, des)
fig = plot_init_stats(des, gs_data=gs_data[ip], ar_data=ar_data[ip],
gs_param=gsp, ar_param=arp,
save_files=save_files, silent=silent)
flist.append(fig)
if return_descriptors:
descs.append({'phase': ip, **(statistical_descriptors or {})})
if return_descriptors or silent: return flist, descs
[docs]
def plot_slice(
self,
cut: str = 'xy',
data: Optional[str] = None,
pos: Optional[Union[str, float, int]] = None,
fname: Optional[Union[str, os.PathLike[str]]] = None,
dual_phase: bool = False,
save_files: bool = False) -> None:
"""
Plot a 2D slice through the microstructure.
The function visualizes a cross-section of the microstructure. If a polygonalized
microstructure is available, it will be used as the plotting basis; otherwise,
or if `data='voxels'`, the voxelized microstructure will be plotted. This method
internally calls `output_ang` with plotting enabled and file writing disabled.
Parameters
----------
cut : str, optional, default='xy'
The cutting plane of the slice. Options are 'xy', 'xz', or 'yz'.
Default is 'xy'.
data : str or None, optional, default=None
Data basis for plotting. Options are 'voxels' or 'poly'. Default is None.
pos : str, float, int, or None, optional, default=None
Position of the slice, either as an absolute value or as one of
'top', 'bottom', 'left', 'right'. Default is None.
fname : str or os.PathLike or None, optional, default=None
Filename to save the figure as a PDF. Default is None.
dual_phase : bool, optional, default=False
If True, enable dual-phase visualization. Default is False.
save_files : bool, optional, default=False
If True, the figure will be saved to disk. Default is False.
Returns
-------
None
Examples
--------
>>> micro.plot_slice(cut='xz', data='poly', pos='top', save_files=True)
"""
self.output_ang(cut=cut, data=data, plot=True, save_files=False,
pos=pos, fname=fname, dual_phase=dual_phase,
save_plot=save_files)
"""
-------- Import/Export methods --------
"""
[docs]
def write_abq(
self,
nodes: Any = None,
file: Optional[Union[str, os.PathLike[str]]] = None,
path: Union[str, os.PathLike[str]] = './',
voxel_dict: Optional[Mapping[Any, Any]] = None,
grain_dict: Optional[Mapping[Any, Any]] = None,
dual_phase: bool = False,
thermal: bool = False,
units: Optional[str] = None,
ialloy: Any = None,
nsdv: int = 360,
crystal_plasticity: Any = None,
phase_props: Optional[Mapping[str, Any]] = None,
boundary_conditions: Optional[Dict[str, Any]] = None,
props_file: Any = None) -> str:
"""
Write the Abaqus input deck (.inp) for the generated RVE
This method generates an Abaqus input file for the current RVE. It supports
voxelized or smoothened meshes, dual-phase materials, crystal plasticity,
thermal analysis, custom material properties, and optional boundary conditions.
Material definitions and mesh data are automatically handled based on provided
arguments or class attributes.
Parameters
----------
nodes : str, array-like, or None, optional, default=None
Defines the mesh to write:
- 'voxels', 'v' : use voxelized mesh
- 'smooth', 's' : use smoothened mesh
- array-like : explicit nodal coordinates
Default is None, automatically selecting available mesh.
file : str or os.PathLike or None, optional, default=None
Filename for the Abaqus input deck. Default is auto-generated.
path : str or os.PathLike, optional, default='./'
Directory path to save the input deck. Default is './'.
voxel_dict : dict or None, optional, default=None
Dictionary with voxel information. Default is `self.mesh.voxel_dict`.
grain_dict : dict or None, optional, default=None
Dictionary mapping grain IDs to nodes. Default is `self.mesh.grain_dict`.
dual_phase : bool, optional, default=False
If True, group elements into phase sets for standard plasticity.
Keep False for grain-wise CP materials, including multiphase RVEs.
thermal : bool, optional, default=False
If True, include thermal material definitions. Default is False.
units : str or None, optional, default=None
Output length units: 'µm', 'um', 'mm', or 'm'. Default is `self.rve.units`.
ialloy : list, object, or None, optional, default=None
Material definitions for each phase. Default is `self.rve.ialloy`.
props_file : str, path-like, list, or None, optional
For multiple phases, both props_file and crystal_plasticity must be
lists of length nphases, indexed by zero-based phase ID. For CP phases,
supply numeric shared constants; each grain's definition in _mat.inp
includes them after the selector, Euler angles and four zeros.
For standard J2 phases, supply an Abaqus material-property include,
written under PHASE{id}_MAT in _geom.inp. A None entry requires False
for crystal_plasticity and writes only the solid-section reference;
define that material yourself in CAE or by editing _geom.inp.
For grain 0, an empty PHASE0_MAT definition is also written so the
assigned material is visible in CAE.
Relative include paths are resolved against path. Omit for legacy output.
nsdv : int, optional
Number of state variables per integration point (default: 360).
crystal_plasticity : bool, list of bool, or None, optional
Select CP per phase. Multiple phases with props_file require a list.
CP requires grain-wise sets (dual_phase=False) and grain orientations.
Grain 0 is reserved for standard plasticity in phase 0. With a
single-phase props_file, None retains the historical CP default;
pass False explicitly for standard plasticity.
phase_props : dict or None, optional, default=None
Additional phase-specific material properties.
boundary_conditions : dict or None, optional, default=None
Dictionary specifying boundary conditions to write. Default is None.
Returns
-------
file : str
Full path to the generated Abaqus input file.
Raises
------
ValueError
- If no voxelized or smoothened mesh is available when required.
- If invalid `nodes` argument is provided.
- If units are not 'µm', 'um', 'mm', or 'm'.
- If the list `ialloy` is longer than the number of phases in the RVE.
- If periodic boundary conditions are requested but the RVE is non-periodic.
Notes
-----
- Automatically selects voxel or smoothened mesh if `nodes` is None.
- Handles dual-phase structures and crystal plasticity input.
- Writes additional material files if orientations are available and `ialloy` is provided.
- Visualization or mesh generation must be performed before calling this function.
Examples
--------
Write voxelized RVE input deck with default settings
abq_file = rve.write_abq(nodes='voxels')
# Write smoothened RVE deck for dual-phase material
abq_file = rve.write_abq(nodes='smooth', dual_phase=True, ialloy=alloy_list)
"""
if nodes is None:
if self.mesh.nodes_smooth is not None and 'GBarea' in self.geometry.keys():
logging.warning('\nWarning: No argument "nodes" is given, will write smoothened structure')
nodes = self.mesh.nodes_smooth
faces = self.geometry['GBarea']
ntag = '_smooth'
elif self.mesh.nodes is not None:
logging.warning('\nWarning: No argument "nodes" is given, will write voxelized structure')
nodes = self.mesh.nodes
faces = None
ntag = '_voxels'
else:
raise ValueError('No information about voxelized microstructure. Run voxelize first.')
elif type(nodes) is not str:
faces = None
ntag = '_voxels'
elif nodes.lower() in ['smooth', 's']:
if self.mesh.nodes_smooth is not None and 'GBarea' in self.geometry.keys():
nodes = self.mesh.nodes_smooth
faces = self.geometry['GBarea'] # use tet elements for smoothened structure
ntag = '_smooth'
else:
raise ValueError('No information about smoothed microstructure. Run smoothen first.')
elif nodes.lower() in ['voxels', 'v', 'voxel']:
if self.mesh.nodes is not None:
nodes = self.mesh.nodes
faces = None # use brick elements for voxel structure
ntag = '_voxels'
else:
raise ValueError('No information about voxelized microstructure. Run voxelize first.')
else:
raise ValueError('Wrong value for parameter "nodes". Must be either "smooth" ' +
f'or "voxels", not {nodes}')
if voxel_dict is None:
voxel_dict = self.mesh.voxel_dict
if units is None:
units = self.rve.units
units = normalize_length_unit(units)
validate_matrix_mapping(getattr(self.mesh, "grain_phase_dict", None))
if dual_phase:
nct = 'abq_dual_phase'
if grain_dict is None:
grain_dict = dict()
for i in range(self.nphases):
grain_dict[i] = list()
for igr, ip in self.mesh.grain_phase_dict.items():
grain_dict[ip] = np.concatenate(
[grain_dict[ip], self.mesh.grain_dict[igr]])
else:
if grain_dict is None:
grain_dict = self.mesh.grain_dict
nct = f'abq_px_{len(grain_dict)}'
if ialloy is None:
ialloy = self.rve.ialloy
files, cp = _abaqus_phase_options(props_file, crystal_plasticity, self.nphases)
if cp is not None:
if dual_phase and any(cp):
raise ValueError('Crystal plasticity requires grain-wise sets (dual_phase=False).')
if 0 in self.mesh.grain_dict and (cp[0]):
raise ValueError('Grain 0 requires standard plasticity in phase 0.')
if any(cp) and self.mesh.grain_ori_dict is None:
raise ValueError('Crystal plasticity requires grain orientations. Run generate_orientations first.')
if props_file is not None:
if ialloy is None:
ialloy = 0
if not isinstance(ialloy, list):
ialloy = [ialloy] * self.nphases
if len(ialloy) != self.nphases:
raise ValueError('ialloy must contain one selector per phase with props_file.')
if type(ialloy) is list and len(ialloy) > self.nphases:
raise ValueError('List of values in ialloy is larger than number of phases in RVE.' +
f'({len(ialloy)} > {self.nphases})')
if self.nphases > 1:
grpd = dict(self.mesh.grain_phase_dict)
else:
grpd = None
if cp is not None and any(cp):
missing = [gid for gid in grain_dict if gid != 0
and cp[(grpd or {}).get(gid, 0)]
and (self.mesh.grain_ori_dict.get(gid) is None
or np.asarray(self.mesh.grain_ori_dict[gid]).shape != (3,))]
if missing:
raise ValueError(f'Missing orientations for CP grains: {missing}')
if boundary_conditions is None:
boundary_conditions = {
"apply_bc": False,
"periodic_bc": self.rve.periodic,
"type_bc": "displacement",
"components_bc": ['*', '*', '*', '*', '*', '*'],
}
elif isinstance(boundary_conditions, dict):
# check consistency of PBC with RVE-type
periodic = self.rve.periodic
periodicBC = boundary_conditions["periodic_bc"]
if not periodic and periodicBC:
raise ValueError("Periodic boundary conditions cannot be applied to a non-periodic RVE.")
else:
raise TypeError('Boundary conditions must be a dictionary.')
if file is None:
if self.name == 'Microstructure':
file = nct + ntag + '_geom.inp'
else:
file = self.name + ntag + '_geom.inp'
path = os.path.normpath(path)
file = os.path.join(path, file)
export2abaqus(nodes, file, grain_dict, voxel_dict,
units=units, gb_area=faces,
dual_phase=dual_phase,
ialloy=ialloy, grain_phase_dict=grpd,
thermal=thermal,
crystal_plasticity=cp if cp is not None else crystal_plasticity,
props_file=props_file,
phase_props=phase_props,
boundary_conditions=boundary_conditions)
# if orientations exist and ialloy is defined also write material file with Euler angles
if (cp is None or any(cp)) and not (self.mesh.grain_ori_dict is None or ialloy is None):
writeAbaqusMat(ialloy, self.mesh.grain_ori_dict,
file=file[0:-8] + 'mat.inp',
grain_phase_dict=grpd, nsdv=nsdv, props_file=props_file,
crystal_plasticity=cp)
return file
[docs]
def write_abq_ori(
self,
ialloy: Any = None,
props_file: Any = None,
ori: Any = None,
file: Optional[Union[str, os.PathLike[str]]] = None,
path: Union[str, os.PathLike[str]] = './',
nsdv: int = 360,
crystal_plasticity=None) -> None:
"""
Write Abaqus material input file using grain orientations for crystal plasticity
This function generates an Abaqus `.inp` file containing material definitions
for crystal plasticity, using the provided grain orientations and material IDs.
If `ialloy` or `ori` are not specified, the method uses the RVE's default values.
Parameters
----------
ialloy : list, object, or None, optional, default=None
Material identifiers for each phase. Default is `self.rve.ialloy`.
props_file : str, path-like, list, or None, optional, default=None
Optional CP-UMAT include file or one include file per phase.
crystal_plasticity : list of bool or None, optional
Per-phase CP flags when using include files; PHASE0 must be False
when GRAIN0 is present. Lists use canonical phase IDs.
ori : dict, array-like, or None, optional, default=None
Dictionary or array of grain orientations. Default is `self.mesh.grain_ori_dict`.
file : str or os.PathLike or None, optional, default=None
Name of the output Abaqus material file. Default is auto-generated based on RVE name.
path : str or os.PathLike, optional, default='./'
Directory path where the file will be saved. Default is './'.
nsdv : int, optional, default=360
Number of state variables per integration point for crystal plasticity. Default is 360.
Returns
-------
None
Writes the Abaqus material input file and does not return its path.
Raises
------
ValueError
- If `ialloy` is not provided and `self.rve.ialloy` is None.
- If `ori` is not provided and `self.mesh.grain_ori_dict` is None.
Notes
-----
- Calls `writeAbaqusMat` internally to generate the material definitions.
- Assumes the RVE and grain orientations have been initialized prior to calling.
Examples
--------
>>> # Write material file using default RVE orientations and material numbers
>>> rve.write_abq_ori()
>>> # Specify custom material IDs and orientations
>>> rve.write_abq_ori(ialloy=alloy_list, ori=ori_dict, file='custom_mat.inp')
"""
if ialloy is None:
ialloy = self.rve.ialloy
if ialloy is None:
raise ValueError('Value of material number in ICAMS CP-UMAT (ialloy) not defined.')
if ori is None:
ori = self.mesh.grain_ori_dict
if ori is None:
raise ValueError('No orientations present. Run "generate_orientations" first.')
if file is None:
if self.name == 'Microstructure':
file = f'abq_px_{self.Ngr}_mat.inp'
else:
file = self.name + '_mat.inp'
path = os.path.normpath(path)
file = os.path.join(path, file)
writeAbaqusMat(ialloy, ori, props_file=props_file, file=file, nsdv=nsdv,
grain_phase_dict=self.mesh.grain_phase_dict,
crystal_plasticity=crystal_plasticity)
[docs]
def output_neper(self) -> None:
"""
Write particle position and weight files for Neper tessellation
This function generates two text files required for tessellation in Neper:
- `sphere_positions.txt` containing the x, y, z coordinates of each particle
- `sphere_weights.txt` containing the weight (radius) of each particle
The files are written in the current working directory. Each particle in
`self.particles` is included in the output.
Raises
------
ValueError
If `self.particles` is None, indicating that particles have not been packed yet.
Notes
-----
- Requires that particles are initialized via the `pack()` method.
- The particle dictionary is built internally but not returned.
Examples
--------
>>> # Write Neper input files for current particles
>>> rve.output_neper()
"""
# write_position_weights(timestep)
if self.particles is None:
raise ValueError('No particle to plot. Run pack first.')
print('')
print('Writing position and weights files for NEPER', end="")
par_dict = dict()
for pa in self.particles:
x, y, z = pa.x, pa.y, pa.z
a = pa.a
par_dict[pa] = [x, y, z, a]
with open('sphere_positions.txt', 'w') as fd:
for key, value in par_dict.items():
fd.write('{0} {1} {2}\n'.format(value[0], value[1], value[2]))
with open('sphere_weights.txt', 'w') as fd:
for key, value in par_dict.items():
fd.write('{0}\n'.format(value[3]))
print('---->DONE!\n')
[docs]
def output_ang(
self,
ori: Any = None,
cut: str = 'xy',
data: Optional[str] = None,
plot: bool = True,
cs: Optional[str] = None,
pos: Optional[Union[str, float, int]] = None,
fname: Optional[Union[str, os.PathLike[str]]] = None,
matname: str = 'XXXX',
save_files: bool = True,
dual_phase: bool = False,
save_plot: bool = False) -> Optional[str]:
"""
Convert microstructure orientations into a .ang file, optionally plotting a slice
This function generates a .ang file mimicking an EBSD map. If polygonalized
microstructure is available, it will be used; otherwise, or if `data='voxels'`,
the voxelized microstructure will be exported. If no orientations are provided,
each grain receives a random Euler angle. The function can also plot a 2D slice
of the orientation map along a specified cutting plane.
Parameters
----------
ori : array-like of shape (self.Ngr, 3), dict, or None, optional, default=None
Euler angles of grains. If None, random angles are generated.
cut : {'xy', 'xz', 'yz'}, optional, default='xy'
Plane along which to take a slice. Default is 'xy'.
data : {'voxels', 'poly'}, or None, optional, default=None
Basis for generating the ANG file. Default is None, automatically selected.
plot : bool, optional, default=True
Whether to display a plot of the slice. Default is True.
cs : str or None, optional, default=None
Crystal symmetry. Default is None. Currently inactive.
pos : float, int, str, or None, optional, default=None
Slice position as a numerical value or one of 'top', 'bottom', 'left', 'right'.
Default is None (uses top/right edge).
fname : str or os.PathLike or None, optional, default=None
Filename for the output ANG file. The value is replaced by an
automatically generated filename when ``save_files=True``.
matname : str, optional, default='XXXX'
Material name to write in the ANG file. Default is 'XXXX'.
save_files : bool, optional, default=True
Whether to save the ANG file. Default is True.
dual_phase : bool, optional, default=False
Whether to output dual-phase information. Default is False.
save_plot : bool, optional, default=False
Whether to save the plot as a PDF. Default is False.
Returns
-------
str or None
Name of the generated ANG file when ``save_files=True``; otherwise
returns the input ``fname`` value.
Raises
------
ValueError
- If `pos` is not a float or one of 'top', 'bottom', 'left', 'right'.
- If `cut` is not 'xy', 'xz', or 'yz'.
- If neither polygonal nor voxel microstructure is available.
Notes
-----
- Uses Delaunay triangulation to determine which grain each pixel belongs to
when exporting polygons.
- Slice plotting uses a colormap (default 'prism') and includes optional
dual-phase coloring.
Examples
--------
>>> # Export an XY slice with default settings
>>> fname = rve.output_ang(cut='xy')
>>> # Export a slice using voxel data and a specific Euler angles array
>>> fname = rve.output_ang(ori=euler_array, data='voxels', cut='xz', pos=10)
>>> # Plot and save the slice as PDF
>>> fname = rve.output_ang(save_plot=True)
"""
if type(ori) is dict:
ori = np.array([val for val in ori.values()])
cut = cut.lower()
if type(pos) is str:
pos = pos.lower()
botlist = ['bottom', 'bot', 'left', 'b', 'l']
toplist = ['top', 'right', 't', 'r']
if cut == 'xy':
sizeX = self.rve.size[0]
sizeY = self.rve.size[1]
(sx, sy, sz) = np.divide(self.rve.size, self.rve.dim)
ix = np.arange(self.rve.dim[0])
iy = np.arange(self.rve.dim[1])
if pos is None or pos in toplist:
iz = self.rve.dim[2] - 1
elif pos in botlist:
iz = 0
elif type(pos) is float or type(pos) is int:
iz = int(pos / sz)
else:
raise ValueError('"pos" must be either float or "top", "bottom", "left" or "right"')
if pos is None:
pos = int(iz * sz)
xl = r'x ($\mu$m)'
yl = r'y ($\mu$m)'
title = r'XY slice at z={} $\mu$m'.format(round(iz * sz, 1))
elif cut == 'xz':
sizeX = self.rve.size[0]
sizeY = self.rve.size[2]
vox_res = np.divide(self.rve.size, self.rve.dim)
sx = vox_res[0]
sy = vox_res[1]
sz = vox_res[2]
ix = np.arange(self.rve.dim[0])
iy = np.arange(self.rve.dim[2])
if pos is None or pos in toplist:
iz = self.rve.dim[1] - 1
elif pos in botlist:
iz = 0
elif type(pos) is float or type(pos) is int:
iz = int(pos / sy)
else:
raise ValueError('"pos" must be either float or "top", "bottom", "left" or "right"')
if pos is None:
pos = int(iz * sz)
xl = r'x ($\mu$m)'
yl = r'z ($\mu$m)'
title = r'XZ slice at y={} $\mu$m'.format(round(iz * sz, 1))
elif cut == 'yz':
sizeX = self.rve.size[1]
sizeY = self.rve.size[2]
vox_res = np.divide(self.rve.size, self.rve.dim)
sx = vox_res[0]
sy = vox_res[1]
sz = vox_res[2]
ix = np.arange(self.rve.dim[1])
iy = np.arange(self.rve.dim[2])
if pos is None or pos in toplist:
iz = self.rve.dim[0] - 1
elif pos in botlist:
iz = 0
elif type(pos) is float or type(pos) is int:
iz = int(pos / sx)
else:
raise ValueError('"pos" must be either float or "top", "bottom", "left" or "right"')
if pos is None:
pos = int(iz * sz)
xl = r'y ($\mu$m)'
yl = r'z ($\mu$m)'
title = r'YZ slice at x={} $\mu$m'.format(round(iz * sz, 1))
else:
raise ValueError('"cut" must bei either "xy", "xz" or "yz".')
# ANG file header
head = ['# TEM_PIXperUM 1.000000\n',
'# x-star 0.000000\n',
'# y-star 0.000000\n',
'# z-star 0.000000\n',
'# WorkingDistance 0.000000\n',
'#\n',
'# Phase 0\n',
'# MaterialName {}\n'.format(matname),
'# Formula\n',
'# Info\n',
'# Symmetry m-3m\n',
'# LatticeConstants 4.050 4.050 4.050 90.000 90.000 90.000\n',
'# NumberFamilies 0\n',
'# ElasticConstants 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000\n',
'# ElasticConstants 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000\n',
'# ElasticConstants 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000\n',
'# ElasticConstants 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000\n',
'# ElasticConstants 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000\n',
'# ElasticConstants 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000\n',
'# Categories0 0 0 0 0\n',
'# \n',
'# GRID: SqrGrid\n',
'# XSTEP: {}\n'.format(round(sx, 6)),
'# YSTEP: {}\n'.format(round(sy, 6)),
'# NCOLS_ODD: {}\n'.format(ix),
'# NCOLS_EVEN: {}\n'.format(ix),
'# NROWS: {}\n'.format(iy),
'#\n',
'# OPERATOR: Administrator\n',
'#\n',
'# SAMPLEID:\n',
'#\n',
'# SCANID:\n',
'#\n'
]
# determine whether polygons or voxels shall be exported
if data is None:
if 'Grains' in self.geometry.keys():
data = 'poly'
elif self.mesh.voxels is None:
raise ValueError('Neither polygons nor voxels for grains are present.\
\nRun voxelize and generate_grains first.')
else:
data = 'voxels'
elif data != 'voxels' and data != 'poly':
raise ValueError('"data" must be either "voxels" or "poly".')
if data == 'voxels':
title += ' (Voxels)'
if cut == 'xy':
g_slice = np.array(self.mesh.grains[:, :, iz], dtype=int)
elif cut == 'xz':
g_slice = np.array(self.mesh.grains[:, iz, :], dtype=int)
else:
g_slice = np.array(self.mesh.grains[iz, :, :], dtype=int)
if dual_phase:
if cut == 'xy':
g_slice_phase = np.array(self.mesh.phases[:, :, iz], dtype=int)
elif cut == 'xz':
g_slice_phase = np.array(self.mesh.phases[:, iz, :], dtype=int)
else:
g_slice_phase = np.array(self.mesh.phases[iz, :, :], dtype=int)
else:
title += ' (Polygons)'
xv, yv = np.meshgrid(ix * sx, iy * sy, indexing='ij')
grain_slice = np.ones(len(ix) * len(iy), dtype=int)
if cut == 'xy':
mesh_slice = np.array([xv.flatten(), yv.flatten(), grain_slice * iz * sz]).T
elif cut == 'xz':
mesh_slice = np.array([xv.flatten(), grain_slice * iz * sz, yv.flatten()]).T
else:
mesh_slice = np.array([grain_slice * iz * sz, xv.flatten(), yv.flatten()]).T
grain_slice = np.zeros(len(ix) * len(iy), dtype=int)
for igr in self.geometry['Grains'].keys():
pts = self.geometry['Grains'][igr]['Points']
try:
tri = Delaunay(pts)
i = tri.find_simplex(mesh_slice)
ind = np.nonzero(i >= 0)[0]
grain_slice[ind] = igr
except Exception as e:
logging.error(f'An unexpected exception occurred: {e}')
logging.error('Grain #{} has no convex hull (Nvertices: {})'
.format(igr, len(pts)))
if np.any(grain_slice == 0):
ind = np.nonzero(grain_slice == 0)[0]
logging.error('Incomplete slicing for {} pixels in {} slice at {}.'
.format(len(ind), cut, pos))
g_slice = grain_slice.reshape(xv.shape)
if save_files:
if ori is None:
ori = np.zeros((self.Ngr, 3))
ori[:, 0] = np.random.rand(self.Ngr) * 2 * np.pi
ori[:, 1] = np.random.rand(self.Ngr) * 0.5 * np.pi
ori[:, 2] = np.random.rand(self.Ngr) * 0.5 * np.pi
# write data to ang file
fname = '{0}_slice_{1}_{2}.ang'.format(cut.upper(), pos, data)
with open(fname, 'w') as f:
f.writelines(head)
for j in iy:
for i in ix:
p1 = ori[g_slice[j, i] - 1, 0]
P = ori[g_slice[j, i] - 1, 1]
p2 = ori[g_slice[j, i] - 1, 2]
f.write(' {0} {1} {2} {3} {4} 0.0 0.000 0 1 0.000\n'
.format(round(p1, 5), round(P, 5), round(p2, 5),
round(sizeX - i * sx, 5), round(sizeY - j * sy, 5)))
if plot:
# plot grains on slice
# cmap = plt.cm.get_cmap('gist_rainbow')
cmap = plt.cm.get_cmap('prism')
fig, ax = plt.subplots(1)
ax.grid(False)
ax.imshow(g_slice, cmap=cmap, interpolation='none',
extent=[0, sizeX, 0, sizeY])
ax.set(xlabel=xl, ylabel=yl)
ax.set_title(title)
if save_plot:
plt.savefig(fname[:-4] + '.pdf', format='pdf', dpi=300)
plt.show()
if dual_phase:
fig, ax = plt.subplots(1)
ax.grid(False)
ax.imshow(g_slice_phase, cmap=cmap, interpolation='none',
extent=[0, sizeX, 0, sizeY])
ax.set(xlabel=xl, ylabel=yl)
ax.set_title(title)
if save_plot:
plt.savefig(fname[:-4] + '.pdf', format='pdf', dpi=300)
plt.show()
return fname
[docs]
def write_stl(
self,
data: str = 'grains',
file: Optional[Union[str, os.PathLike[str]]] = None,
path: Union[str, os.PathLike[str]] = './',
phases: bool = False,
phase_num: Optional[int] = None) -> None:
"""
Export grains or particles as STL files representing convex polyhedra
This function writes STL files with triangular facets for grains or particles.
Each facet is written in standard STL ASCII format:
```
solid name
facet normal n1 n2 n3
outer loop
vertex p1x p1y p1z
vertex p2x p2y p2z
vertex p3x p3y p3z
endloop
endfacet
endsolid name
```
Parameters
----------
data : {'grains', 'particles'}, str, optional, default='grains'
Determines whether to export grains or particles. Default is 'grains'.
file : str or os.PathLike or None, optional, default=None
Filename for the STL file. Default is generated from `self.name`.
path : str or os.PathLike, optional, default='./'
Directory to save the STL file. Default is './'.
phases : bool, optional, default=False
If True, export only grains of a specific phase. Default is False.
phase_num : int or None, optional, default=None
Phase number to export if `phases=True`. Required in that case.
Returns
-------
None
Raises
------
ValueError
- If `phases=True` but `phase_num` is not provided.
Notes
-----
- Facet normals are automatically computed using cross products of vertices.
- Acute or irregular facets trigger warnings.
- For particle export, each particle must have a valid inner polyhedron.
Examples
--------
>>> # Export all grains to STL
>>> rve.write_stl()
>>> # Export particles to STL
>>> rve.write_stl(data='particles')
>>> # Export only phase 1 grains
>>> rve.write_stl(phases=True, phase_num=1)
"""
def write_facet(nv: Any, pts: Any, ft: int) -> None:
"""
Write a single triangular facet to the STL file with a normalized normal vector
Checks for degenerate or nearly zero-length normals and issues warnings
if the facet is acute or irregular. Writes the facet in ASCII STL format.
Parameters
----------
nv : ndarray
Initial normal vector of the facet.
pts : ndarray
3x3 array of vertex coordinates defining the facet.
ft : int
Facet index, used for logging warnings.
"""
if np.linalg.norm(nv) < 1.e-5:
logging.warning(f'Acute facet detected. Facet: {ft}')
nv = np.cross(pts[1] - pts[0], pts[2] - pts[1])
if np.linalg.norm(nv) < 1.e-5:
logging.warning(f'Irregular facet detected. Facet: {ft}')
nv /= np.linalg.norm(nv)
f.write(" facet normal {} {} {}\n"
.format(nv[0], nv[1], nv[2]))
f.write(" outer loop\n")
f.write(" vertex {} {} {}\n"
.format(pts[0, 0], pts[0, 1], pts[0, 2]))
f.write(" vertex {} {} {}\n"
.format(pts[1, 0], pts[1, 1], pts[1, 2]))
f.write(" vertex {} {} {}\n"
.format(pts[2, 0], pts[2, 1], pts[2, 2]))
f.write(" endloop\n")
f.write(" endfacet\n")
def write_grains() -> None:
"""
Write all grain facets of the microstructure to the STL file
Iterates over all facets defined in `self.geometry['Facets']`, computes
the facet normal, and calls `write_facet` to output each triangular facet.
"""
for ft in self.geometry['Facets']:
pts = self.geometry['Points'][ft]
nv = np.cross(pts[1] - pts[0], pts[2] - pts[0]) # facet normal
write_facet(nv, pts, ft)
def write_phases(ip: int) -> None:
"""
Write facets of grains belonging to a specific phase to the STL file
Parameters
----------
ip : int
Phase number for which the grain facets should be exported.
Notes
-----
Iterates over all grains in `self.geometry['Grains']` and writes only
the facets of grains whose 'Phase' matches `ip`. Each facet is output
using the `write_facet` function.
"""
for grain in self.geometry['Grains'].values():
if grain['Phase'] == ip:
for ft in grain['Simplices']:
pts = self.geometry['Points'][ft]
nv = np.cross(pts[1] - pts[0], pts[2] - pts[0]) # facet normal
write_facet(nv, pts, ft)
def write_particles() -> None:
"""
Write facets of all particle convex hulls to the STL file
Notes
-----
Iterates over all particles in `self.particles` and exports each facet
from the particle's inner convex hull using the `write_facet` function.
"""
for pa in self.particles:
for ft in pa.inner.convex_hull:
pts = pa.inner.points[ft]
nv = np.cross(pts[1] - pts[0], pts[2] - pts[0]) # facet normal
write_facet(nv, pts, ft)
if file is None:
if self.name == 'Microstructure':
file = 'px_{}grains.stl'.format(self.Ngr)
else:
file = self.name + '.stl'
path = os.path.normpath(path)
file = os.path.join(path, file)
with open(file, 'w') as f:
f.write("solid {}\n".format(self.name))
if data in ['particles', 'pa', 'p']:
if self.particles[0].inner is None:
logging.error("Particles don't contain inner polyhedron, cannot write STL file.")
else:
for pa in self.particles:
pa.sync_poly()
write_particles()
elif data in ['grains', 'gr', 'g']:
if phases:
if phase_num is None:
raise ValueError('Phase-specific output requested, but no phase number specified.')
write_phases(phase_num)
else:
write_grains()
else:
raise ValueError(f"Invalid data type specified for STL export, must be either 'particles' or 'grains', not {data}.")
f.write("endsolid\n")
return
[docs]
def write_centers(
self,
file: Optional[Union[str, os.PathLike[str]]] = None,
path: Union[str, os.PathLike[str]] = './',
grains: Optional[Mapping[Any, Any]] = None) -> None:
"""
Write the center positions of grains to a CSV file
Parameters
----------
file : str or os.PathLike or None, optional, default=None
Filename for the CSV file. Default is auto-generated based on `self.name`.
path : str or os.PathLike, optional, default='./'
Directory path to save the CSV file. Default is './'.
grains : dict or None, optional, default=None
Dictionary of grain data. Default is `self.geometry['Grains']`.
Returns
-------
None
Notes
-----
- Each line in the CSV file corresponds to a grain center in the format: x, y, z.
- Polyhedral grains without simplices are skipped.
"""
if file is None:
if self.name == 'Microstructure':
file = 'px_{}grains_centroid.csv'.format(self.Ngr)
else:
file = self.name + '_centroid.csv'
path = os.path.normpath(path)
file = os.path.join(path, file)
if grains is None:
grains = self.geometry['Grains']
with open(file, 'w') as f:
for gr in grains.values():
# if polyhedral grain has no simplices, center should not be written!!!
ctr = gr['Center']
f.write('{}, {}, {}\n'.format(ctr[0], ctr[1], ctr[2]))
return
[docs]
def write_ori(
self,
angles: Any = None,
file: Optional[Union[str, os.PathLike[str]]] = None,
path: Union[str, os.PathLike[str]] = './') -> None:
"""
Write grain orientations (Euler angles) to a CSV file
This method exports the orientations of grains in the microstructure. If no angles
are provided, the method uses the stored grain orientations from the mesh. The CSV
file contains one row per grain with three values representing the Euler angles
(phi1, Phi, phi2).
Parameters
----------
angles : array-like of shape (N, 3) or None, optional, default=None
Euler angles for each grain. If None, the stored orientations are used.
file : str or os.PathLike or None, optional, default=None
Name of the output CSV file. If None, a default name based on the microstructure
name and number of grains is used.
path : str or os.PathLike, optional, default='./'
Directory path where the CSV file will be saved. Default is './'.
Returns
-------
None
Writes the CSV file and does not return its path.
Raises
------
ValueError
If `angles` is None and no grain orientations are stored in `self.mesh.grain_ori_dict`.
Notes
-----
- Each row in the CSV corresponds to one grain.
- Values are written in the order: phi1, Phi, phi2 (in radians).
Examples
--------
>>> # Export stored orientations to default file
>>> rve.write_ori()
>>> # Export provided orientations to a specific file
>>> rve.write_ori(angles=my_angles, file='grain_orientations.csv', path='./data')
"""
if file is None:
if self.name == 'Microstructure':
file = 'px_{}grains_ori.csv'.format(self.Ngr)
else:
file = self.name + '_ori.csv'
path = os.path.normpath(path)
file = os.path.join(path, file)
if angles is None:
if self.mesh.grain_ori_dict is None:
raise ValueError('No grain orientations given or stored.')
angles = [val for val in self.mesh.grain_ori_dict.values()]
with open(file, 'w') as f:
for ori in angles:
f.write('{}, {}, {}\n'.format(ori[0], ori[1], ori[2]))
return
[docs]
def write_voxels(
self,
angles: Any = None,
script_name: Optional[str] = None,
file: Optional[Union[str, os.PathLike[str]]] = None,
path: Union[str, os.PathLike[str]] = './',
mesh: bool = True,
source: Any = None,
system: bool = False) -> None:
"""
Write voxel structure of the microstructure to a JSON file
This function exports the voxel-based representation of the microstructure,
including grain numbers, optional orientations, and mesh data. Metadata such
as owner, system information, and script name are included.
Parameters
----------
angles : np.ndarray or None, optional, default=None
Array of Euler angles for grain orientations. If None, stored orientations
are used if available. Default is None.
script_name : str or None, optional, default=None
Name of the script used to generate the microstructure. Default is the
current script file.
file : str or os.PathLike or None, optional, default=None
Name of the output JSON file. Default is based on microstructure name and
number of grains.
path : str or os.PathLike, optional, default='./'
Directory path to save the JSON file. Default is current directory './'.
mesh : bool, optional, default=True
Whether to include mesh nodes and voxel connectivity in the output. Default
is True.
source : Any, optional, default=None
Optional input or source information to include in the JSON metadata.
system : bool, optional, default=False
Whether to include system information in the JSON metadata. Default is False.
Returns
-------
None
The function writes a JSON file and does not return any value.
"""
import platform
import getpass
from datetime import date
from kanapy import __version__
if script_name is None:
script_name = __file__
if file is None:
if self.name == 'Microstructure':
file = f'px_{self.Ngr}grains_voxels.json'
else:
file = self.name + '_voxels.json'
path = os.path.normpath(path)
file = os.path.join(path, file)
print(f'Writing voxel information of microstructure to {file}.')
# metadata
today = str(date.today()) # date
owner = getpass.getuser() # username
sys_info = platform.uname() # system information
# output dict
structure = {
"Info": {
"Owner": owner,
"Institution": "ICAMS, Ruhr University Bochum, Germany",
"Date": today,
"Description": "Voxels of microstructure",
"Method": "Synthetic microstructure generator Kanapy",
},
"Model": {
"Creator": "kanapy",
"Version": __version__,
"Repository": "https://github.com/ICAMS/Kanapy.git",
"Input": source,
"Script": script_name,
"Material": self.name,
"Phase_names": self.rve.phase_names,
"Size": [int(val) for val in self.rve.size],
"Periodicity": str(self.rve.periodic),
"Units": {
'Length': self.rve.units,
},
},
"Data": {
"Description": 'Grain numbers per voxel',
"Type": 'int',
"Shape": self.rve.dim,
"Order": 'C',
"Values": [int(val) for val in self.mesh.grains.flatten()],
},
"Grains": {
"Description": "Grain-related data",
"Orientation": "Euler-Bunge angle",
"Phase": "Phase number"
},
}
if system:
structure["Info"]["System"] = {
"sysname": sys_info[0],
"nodename": sys_info[1],
"release": sys_info[2],
"version": sys_info[3],
"machine": sys_info[4],
}
for igr in self.mesh.grain_dict.keys():
structure["Grains"][int(igr)] = {
"Phase": int(self.mesh.grain_phase_dict[igr])
}
if angles is None:
if self.mesh.grain_ori_dict is None:
logging.info('No angles for grains are given. Writing only geometry of RVE.')
else:
for igr in self.mesh.grain_ori_dict.keys():
structure["Grains"][igr]["Orientation"] = list(self.mesh.grain_ori_dict[igr])
else:
for i, igr in enumerate(self.mesh.grain_dict.keys()):
structure["Grains"][igr]["Orientation"] = list(angles[i, :])
if mesh:
structure['Mesh'] = {
"Nodes": {
"Description": 'Nodal coordinates',
"Type": 'float',
"Shape": self.mesh.nodes.shape,
"Values": [list(val) for val in self.mesh.nodes],
},
"Voxels": {
"Description": 'Node list per voxel',
"Type": 'int',
"Shape": (len(self.mesh.voxel_dict.keys()), 8),
"Values": [val for val in self.mesh.voxel_dict.values()],
}
}
# write file
with open(file, 'w') as fp:
json.dump(structure, fp)
return
[docs]
def create_microstructure_identifier(
self,
microstructure_step: Mapping[str, Any],
hash_length: int = 8,
sig_figs: int = 6,
) -> str:
"""
Create a deterministic identifier for one microstructure snapshot.
The identifier is based on the microstructure content: grid, grains,
voxels, phase IDs, orientations, voxel-grain mapping, voxel indices,
voxel coordinates, and voxel volumes. It is written as the prefix
``"S_"`` followed by the leading ``hash_length`` hexadecimal characters
of a SHA-256 digest, per the microstructure_state_id convention.
Only undeformed, simulation-ready snapshots should be given an
identifier -- this function itself is content-only and does not inspect
``microstructure_step["grid"]["status"]``; the caller is responsible for
only calling it on undeformed/regridded snapshots, never on intermediate
deformed states.
Parameters
----------
microstructure_step : Mapping[str, Any]
One microstructure snapshot dictionary.
hash_length : int, optional, default=8
Number of hexadecimal characters taken from the SHA-256 digest
(after the ``"S_"`` prefix). Default is 8.
sig_figs : int, optional, default=6
Number of significant figures each float is rounded to before
hashing. Default is 6.
Returns
-------
str
Deterministic identifier, for example ``"S_bb1f711a"``.
Notes
-----
Rounding is done to `sig_figs` **significant figures**, not decimal
places. This schema stores lengths/volumes in SI base units (meters,
cubic meters), where grain/voxel volumes are ~1e-17-1e-19 and voxel
centroid coordinates are ~1e-7-1e-5. A fixed-decimal-places round (e.g.
``round(value, 6)``) collapses all such values to exactly ``0.0``,
making the hash blind to grain volume entirely and aliasing many
distinct voxel positions onto the same rounded value -- confirmed
empirically against a real MiMeDO snapshot (all 24 distinct grain
volumes and many of 8000 distinct voxel centroids collapsed together).
Rounding to significant figures instead keeps every field distinguishable
regardless of its physical unit scale, so the hash reproduces the same
value from the same content every time, and only from that content, as
required by a content-derived identifier.
"""
if not isinstance(microstructure_step, Mapping):
raise TypeError("microstructure_step must be a dictionary-like mapping.")
if hash_length <= 0:
raise ValueError("hash_length must be a positive integer.")
def _is_numpy_array(value: Any) -> bool:
return hasattr(value, "shape") and hasattr(value, "tolist")
def _is_numpy_scalar(value: Any) -> bool:
return hasattr(value, "item") and not _is_numpy_array(value)
def _is_empty(value: Any) -> bool:
"""
NumPy-safe emptiness check.
Do not compare arrays using value == [] or value == ''.
"""
if value is None:
return True
if isinstance(value, str):
return value == ""
if _is_numpy_array(value):
return value.size == 0
if isinstance(value, Mapping):
return len(value) == 0
if isinstance(value, (list, tuple, set)):
return len(value) == 0
return False
def _round_sig(value: float, sig: int) -> float:
"""Round to `sig` significant figures, magnitude-independent."""
if value == 0.0:
return 0.0
return float(f"{value:.{sig}g}")
def _make_json_safe(value: Any) -> Any:
"""
Convert NumPy/Python objects into deterministic JSON-safe values.
"""
# NumPy arrays first
if _is_numpy_array(value):
return _make_json_safe(value.tolist())
# NumPy scalars
if _is_numpy_scalar(value):
try:
return _make_json_safe(value.item())
except Exception:
pass
if isinstance(value, float):
return _round_sig(value, sig_figs)
if isinstance(value, bool):
return value
if isinstance(value, int):
return value
if isinstance(value, str):
return value
if isinstance(value, Mapping):
cleaned = {}
for key, item in value.items():
if key in ("id", "microstructure_state_id"):
continue
safe_item = _make_json_safe(item)
if not _is_empty(safe_item):
cleaned[str(key)] = safe_item
return cleaned
if isinstance(value, (list, tuple, set)):
cleaned = []
for item in value:
safe_item = _make_json_safe(item)
if not _is_empty(safe_item):
cleaned.append(safe_item)
return cleaned
return value
def _safe_int(value: Any) -> int:
try:
return int(value)
except (TypeError, ValueError):
return 0
grains = microstructure_step.get("grains", [])
voxels = microstructure_step.get("voxels", [])
phase_ids = sorted(
{
grain.get("phase_id")
for grain in grains
if isinstance(grain, Mapping) and grain.get("phase_id") is not None
}
)
grain_payload = [
{
"grain_id": grain.get("grain_id"),
"phase_id": grain.get("phase_id"),
"grain_volume": grain.get("grain_volume"),
"orientation": grain.get("orientation"),
}
for grain in grains
if isinstance(grain, Mapping)
]
voxel_payload = [
{
"voxel_id": voxel.get("voxel_id"),
"grain_id": voxel.get("grain_id"),
"centroid_coordinates": voxel.get("centroid_coordinates"),
"voxel_index": voxel.get("voxel_index"),
"voxel_volume": voxel.get("voxel_volume"),
"orientation": voxel.get("orientation"),
}
for voxel in voxels
if isinstance(voxel, Mapping)
]
payload = {
"grid": microstructure_step.get("grid"),
"grain_count": len(grain_payload),
"voxel_count": len(voxel_payload),
"phase_count": len(phase_ids),
"phase_ids": phase_ids,
"grains": sorted(
grain_payload,
key=lambda item: _safe_int(item.get("grain_id")),
),
"voxels": sorted(
voxel_payload,
key=lambda item: _safe_int(item.get("voxel_id")),
),
}
canonical_payload = json.dumps(
_make_json_safe(payload),
sort_keys=True,
separators=(",", ":"),
ensure_ascii=True,
allow_nan=False,
)
digest = hashlib.sha256(canonical_payload.encode("utf-8")).hexdigest()[:hash_length]
return f"S_{digest}"
[docs]
def write_data(
self,
user_metadata: Optional[Dict[str, Any]] = None,
boundary_condition: Optional[Dict[str, Any]] = None,
phases: Optional[Union[Dict[str, Any], List[Dict[str, Any]]]] = None,
interactive: bool = True,
structured: bool = True,
ialloy: int = 0,
length_unit: str = 'µm') -> Dict[str, Any]:
"""
Generate a JSON-compatible data schema containing user, system, and job-specific elements
This function builds a comprehensive data dictionary for a microstructure simulation,
including metadata, boundary conditions, phase information, voxel data, and grain data. The exported initial voxel records include both grain_id and phase_id.
It supports interactive prompts for missing information, structured/unstructured meshes,
and unit scaling.
Parameters
----------
user_metadata : dict or None, optional, default=None
Prefilled metadata fields; required if interactive=False. Default is None.
boundary_condition : dict or None, optional, default=None
Dictionary specifying mechanical or thermal boundary conditions. Default is None.
phases : dict, list of dicts, or None, optional, default=None
Phase-specific information; if None, default material library and ialloy are used.
interactive : bool, optional, default=True
If True, prompt user for missing information. Default is True.
structured : bool, optional, default=True
Whether the mesh is structured. Default is True.
ialloy : int, optional, default=0
Alloy index for selecting material properties from the built-in library. Default is 0.
length_unit : {'µm', 'um', 'mm', 'm'}, optional, default='µm'
Unit for length scaling.
Returns
-------
dict
A nested dictionary containing the complete data schema for the simulation.
Raises
------
ValueError
If interactive=False and required metadata fields or boundary conditions are missing,
or if length_unit is invalid
TypeError
If `phases` is not a dict or a list of dicts when provided
"""
length_scale = length_scale_from_um(length_unit)
# Material library definitions (pulled from mod_alloys.f)
material_library = {
1: { # Aluminum
'ialloy': 1,
'material_identifier': 'Aluminum',
'elastic_model_name': 'Anisotropic Elasticity',
'elastic_parameters': {
'C11': 247000.0, 'C12': 147000.0, 'C44': 125000.0
},
'plastic_model_name': 'Crystal Plasticity',
'plastic_parameters': {
'reference_shear_rate': 1e-6, # shrt0
'initial_critical_resolved_shear_stress': 20.0, # crss0
'saturated_slip_resistance': 1500.0, # crsss
'strain_rate_sensitivity_exponent': 20.0, # pwfl
'reference_hardening_rate': 60.0, # hdrt0
'hardening_exponent': 2.25 # pwhd
}
},
2: { # Copper
'ialloy': 2,
'material_identifier': 'Copper',
'elastic_model_name': 'Anisotropic Elasticity',
'elastic_parameters': {
'C11': 170000.0, 'C12': 124000.0, 'C44': 75000.0
},
'plastic_model_name': 'Crystal Plasticity',
'plastic_parameters': {
'reference_shear_rate': 0.001, # shrt0
'initial_critical_resolved_shear_stress': 16.0, # crss0
'saturated_slip_resistance': 148.0, # crsss
'strain_rate_sensitivity_exponent': 83.0, # pwfl
'reference_hardening_rate': 250.0, # hdrt0
'hardening_exponent': 2.25 # pwhd
}
},
3: { # Ferrite
'ialloy': 3,
'material_identifier': 'Ferrite',
'elastic_model_name': 'Anisotropic Elasticity',
'elastic_parameters': {
'C11': 230000.0, 'C12': 135000.0, 'C44': 116000.0
},
'plastic_model_name': 'Crystal Plasticity',
'plastic_parameters': {
'reference_shear_rate': 1e-6,
'initial_critical_resolved_shear_stress': 25.0,
'saturated_slip_resistance': 1600.0,
'strain_rate_sensitivity_exponent': 18.0,
'reference_hardening_rate': 70.0,
'hardening_exponent': 2.0
}
},
4: { # Austenite
'ialloy': 4,
'material_identifier': 'Austenite',
'elastic_model_name': 'Anisotropic Elasticity',
'elastic_parameters': {
'C11': 190000.0, 'C12': 130000.0, 'C44': 115000.0
},
'plastic_model_name': 'Crystal Plasticity',
'plastic_parameters': {
'reference_shear_rate': 1e-6,
'initial_critical_resolved_shear_stress': 15.0,
'saturated_slip_resistance': 1400.0,
'strain_rate_sensitivity_exponent': 22.0,
'reference_hardening_rate': 65.0,
'hardening_exponent': 2.2
}
},
5: { # Superalloy
'ialloy': 5,
'material_identifier': 'Superalloy',
'elastic_model_name': 'Anisotropic Elasticity',
'elastic_parameters': {
'C11': 260000.0, 'C12': 150000.0, 'C44': 120000.0
},
'plastic_model_name': 'Crystal Plasticity',
'plastic_parameters': {
'reference_shear_rate': 1e-6,
'initial_critical_resolved_shear_stress': 30.0,
'saturated_slip_resistance': 1700.0,
'strain_rate_sensitivity_exponent': 25.0,
'reference_hardening_rate': 80.0,
'hardening_exponent': 2.8
}
},
6: { # Nickel
'ialloy': 6,
'material_identifier': 'Nickel',
'elastic_model_name': 'Anisotropic Elasticity',
'elastic_parameters': {
'C11': 246000.0, 'C12': 147000.0, 'C44': 124000.0
},
'plastic_model_name': 'Crystal Plasticity',
'plastic_parameters': {
'reference_shear_rate': 1e-6,
'initial_critical_resolved_shear_stress': 18.0,
'saturated_slip_resistance': 1450.0,
'strain_rate_sensitivity_exponent': 20.0,
'reference_hardening_rate': 75.0,
'hardening_exponent': 2.4
}
}
}
# Define required fields
required_fields = [
'title',
'creator', 'creator_ORCID', 'creator_affiliation', 'creator_institute', 'creator_group',
'contributor', 'contributor_ORCID', 'contributor_affiliation', 'contributor_institute', 'contributor_group',
'date', 'shared_with', 'description', 'rights', 'rights_holder', 'funder_name', 'fund_identifier',
'publisher', 'relation', 'keywords'
]
def prompt_list(field_name: str) -> List[str]:
"""
Prompt the user to enter a comma-separated list of values
This function requests user input for a given field name and splits
the input string into a list of stripped strings. Returns an empty
list if the input is blank.
Parameters
----------
field_name : str
Name of the field to prompt the user for
Returns
-------
list of str
List of values entered by the user, with leading/trailing whitespace removed
"""
vals = input(f"Enter comma-separated {field_name}: ").strip()
return [v.strip() for v in vals.split(',')] if vals else []
# Gather metadata
if interactive and user_metadata is None:
use: Dict[str, Any] = {}
# identifier
ident = input("Identifier (leave blank to auto-generate): ").strip()
if not ident:
now = datetime.utcnow().isoformat()
ident = hashlib.sha256(now.encode()).hexdigest()[:8]
# use['identifier'] = ident
use['title'] = input("Title: ").strip()
# creator fields
use['creator'] = prompt_list('creator names (e.g. Last, First)')
use['creator_ORCID'] = prompt_list('creator ORCID(s)')
use['creator_affiliation'] = prompt_list('creator affiliations')
use['creator_institute'] = prompt_list('creator institutes')
use['creator_group'] = prompt_list('creator groups')
# contributor fields
use['contributor'] = prompt_list('contributor names')
use['contributor_ORCID'] = prompt_list('contributor ORCID(s)')
use['contributor_affiliation'] = prompt_list('contributor affiliations')
use['contributor_institute'] = prompt_list('contributor institutes')
use['contributor_group'] = prompt_list('contributor groups')
use['date'] = input("Date (YYYY-MM-DD): ").strip() or datetime.utcnow().strftime('%Y-%m-%d')
# shared_with
shared: List[Dict[str, str]] = []
print("Enter shared_with access entries. Valid types: c, u, g, all. Blank to stop.")
while True:
atype = input(" access_type: ").strip()
if not atype:
break
shared.append({'access_type': atype})
use['shared_with'] = shared
use['description'] = input("Description: ").strip()
use['rights'] = input("Rights (e.g. Creative Commons Attribution 4.0 International): ").strip()
use['rights_holder'] = prompt_list('rights_holder')
# other fields
use['funder_name'] = input("Funder name: ").strip()
use['fund_identifier'] = input("Fund identifier: ").strip()
use['publisher'] = input("Publisher: ").strip()
use['relation'] = prompt_list('relation (DOI or URL)')
use['keywords'] = prompt_list('keywords')
else:
if user_metadata is None or len(user_metadata) == 0:
raise ValueError("user_metadata must be provided when interactive is False.")
use = user_metadata.copy()
# Validate presence of required fields
missing = [f for f in required_fields if f not in use]
if missing:
raise ValueError(f"Missing required metadata fields: {', '.join(missing)}")
ig = {
'RVE_size': [float(v) * length_scale for v in self.rve.size],
'RVE_continuity': self.rve.periodic,
'discretization_type': 'Structured' if structured else 'Unstructured',
'discretization_unit_size': [(float(s) / float(d)) * length_scale for s, d in
zip(self.rve.size, self.rve.dim)],
'discretization_count': int(self.mesh.nvox),
'origin': {
'software': 'kanapy',
'software_version': pkg_version('kanapy'),
'system': platform.system(),
'system_version': platform.version()
}
}
# ─── Show vertex‐diagram for BC reference ────────────────────────────────
try:
# locate the package root, two levels up from this file
script_dir = os.path.dirname(__file__)
project_root = os.path.abspath(os.path.join(script_dir, os.pardir, os.pardir))
img_path = os.path.join(project_root, 'docs', 'figs', 'RVE', 'Vertices.png')
if os.path.exists(img_path):
# first try with Pillow
try:
from PIL import Image
Image.open(img_path).show()
except ImportError:
# fallback to the system browser/viewer
import webbrowser
webbrowser.open(f'file://{img_path}')
except (FileNotFoundError, OSError) as e:
logging.warning(f"[Warning] Unable to open Vertices.png: {e}")
# Boundary conditions
if boundary_condition is not None and len(boundary_condition) > 0:
# Ensure 'type' is present
if 'mechanical_BC' in boundary_condition:
# Determine mechanical vs thermal by key presence
mech_list = boundary_condition['mechanical_BC']
# Normalize to list
if not isinstance(mech_list, list):
mech_list = [mech_list]
job_bc = {'mechanical_BC': mech_list}
elif 'thermal_BC' in boundary_condition:
th_list = boundary_condition['thermal_BC']
if not isinstance(th_list, list):
th_list = [th_list]
job_bc = {'thermal_BC': th_list}
else:
# fallback if wrong keys provided
if interactive:
# ask user which type to populate
bc_choice = input("Boundary condition key not found. Enter 'mechanical' or 'thermal': ").strip()
if bc_choice == 'mechanical':
job_bc = {'mechanical_BC': []}
else:
job_bc = {'thermal_BC': []}
else:
raise ValueError(
"boundary_condition dict must include 'mechanical_BC' or 'thermal_BC' key when interactive=False.")
elif interactive: # Interactive entry of multiple mechanical BCs
mech_entries: List[Dict[str, Any]] = []
bc_type = input("Boundary condition type ('mechanical' or 'thermal'): ").strip()
if bc_type == 'mechanical':
while True:
print("Define a mechanical BC (leave vertex_list blank to stop):")
vertex_list = prompt_list('vertex_list')
if not vertex_list:
break
constraints = input("Constraints xyz (e.g. 'free,fixed,loaded'): ").split(',')
loading_type = input("Loading type (force/displacement/stress/strain/none): ").strip()
loading_mode = input("Loading mode (static/cyclic/monotonic/intermittent): ").strip()
# count how many of those constraints are “loaded”:
num = sum(1 for c in constraints if c.strip().lower() == 'loaded')
loads = []
for i in range(num):
print(f"Load entry #{i + 1}:")
mag = float(input(" magnitude: "))
freq = float(input(" frequency: "))
dur = float(input(" duration: "))
R = float(input(" R: "))
loads.append({'magnitude': mag, 'frequency': freq, 'duration': dur, 'R': R})
mech_entries.append({
'vertex_list': vertex_list,
'constraints': constraints,
'loading_type': loading_type,
'loading_mode': loading_mode,
'applied_load': loads
})
job_bc = {'mechanical_BC': mech_entries}
else:
job_bc = {'thermal_BC': []}
else:
job_bc = {}
# Phase data
phase_list = []
if phases is not None and len(phases) > 0: # user provided a dict or list of dicts
if isinstance(phases, dict):
phase_list = [phases]
elif isinstance(phases, list):
phase_list = phases
else:
raise TypeError("`phases` must be a dict or list of dicts.")
else: # fallback: use ialloy + material_library
if not ialloy or ialloy not in material_library:
if interactive:
ialloy = int(input(f"Choose ialloy from {list(material_library.keys())}: "))
else:
raise ValueError(
f"No phases provided and invalid ialloy. "
f"Valid ialloy values: {list(material_library.keys())}"
)
for idx in range(self.nphases):
mat = material_library[ialloy]
pe = mat['elastic_parameters']
pp = mat['plastic_parameters']
cs = self.rve.ms_data[idx]["cs"]
phase_entry = {
"id": idx,
"phase_name": self.rve.phase_names[idx],
"volume_fraction": self.rve.phase_vf[idx],
"lattice_structure": None,
"constitutive_model": {
"$schema": "http://json-schema.org/draft-04/schema#",
"elastic_model_name": mat['elastic_model_name'],
"elastic_parameters": pe,
"plastic_model_name": mat['plastic_model_name'],
"plastic_parameters": pp
},
"orientation": {
"euler_angles": {
"Phi1": None,
"Phi": None,
"Phi2": None,
},
"grain_count": self.mesh.ngrains_phase[idx],
"orientation_identifier": None,
"texture_type": None,
"group_symmetry": cs.name,
"software": "orix",
"software_version": orix.__version__,
}
}
phase_list.append(phase_entry)
# ─── pull Mesh + RVE into locals ─────────────────────────────────────────
# ─── pull Mesh + RVE into locals ─────────────────────────────────────────
grain_phase_dict = getattr(self.mesh, "grain_phase_dict", None)
if grain_phase_dict is None:
grain_phase_dict = {}
grain_ori_dict = getattr(self.mesh, "grain_ori_dict", None)
vox_center_dict = getattr(self.mesh, "vox_center_dict", None)
if vox_center_dict is None:
vox_center_dict = {}
grain_to_voxels = getattr(self.mesh, "grain_dict", None)
if grain_to_voxels is None:
grain_to_voxels = {}
rve_size = getattr(self.rve, "size", None)
if rve_size is None:
rve_size = [0, 0, 0]
else:
rve_size = list(rve_size)
rve_dim = getattr(self.rve, "dim", None)
if rve_dim is None:
rve_dim = [1, 1, 1]
else:
rve_dim = list(rve_dim)
# Is orientation available?
include_orientation = isinstance(grain_ori_dict, dict) and len(grain_ori_dict) > 0
# ─── compute one‐voxel volume ─────────────────────────────────────────────
unit_sizes = []
for s, d in zip(rve_size, rve_dim):
try:
unit_sizes.append((float(s) / float(d)) * length_scale if float(d) != 0 else 0.0)
except Exception:
unit_sizes.append(0.0)
voxel_volume = (unit_sizes[0] if len(unit_sizes) > 0 else 0.0) \
* (unit_sizes[1] if len(unit_sizes) > 1 else 0.0) \
* (unit_sizes[2] if len(unit_sizes) > 2 else 0.0)
# grid size per axis
Nx, Ny, Nz = int(rve_dim[0]), int(rve_dim[1]), int(rve_dim[2])
# infer origin from centers if available: origin = min(center) - 0.5*unit
if len(vox_center_dict) > 0:
xs = [float(c[0]) * length_scale for c in vox_center_dict.values()]
ys = [float(c[1]) * length_scale for c in vox_center_dict.values()]
zs = [float(c[2]) * length_scale for c in vox_center_dict.values()]
ox = (min(xs) if xs else 0.0) - 0.5 * unit_sizes[0]
oy = (min(ys) if ys else 0.0) - 0.5 * unit_sizes[1]
oz = (min(zs) if zs else 0.0) - 0.5 * unit_sizes[2]
else:
ox = oy = oz = 0.0
origin = [ox, oy, oz]
# ─── precompute voxel→grain lookup ───────────────────────────────────────
voxel_to_grain = {
vid: gid
for gid, vids in grain_to_voxels.items()
for vid in vids
}
# ─── build time-0 grid dict ────────────────────────────────────────────
grid_t0 = {
"status": "undeformed",
"grid_size": [float(v) * length_scale for v in self.rve.size],
"grid_spacing": [(float(s) / float(d)) * length_scale for s, d in zip(self.rve.size, self.rve.dim)],
}
# ─── build time‐0 grain dict ────────────────────────────────────────────
grains_t0 = []
for gid in grain_phase_dict.keys():
entry = {
"grain_id": gid,
"phase_id": grain_phase_dict.get(gid),
"grain_volume": len(grain_to_voxels.get(gid, [])) * voxel_volume,
}
if include_orientation:
ori = grain_ori_dict.get(gid)
if ori is not None:
entry["orientation"] = list(ori)
grains_t0.append(entry)
grains_t0_sorted = sorted(grains_t0, key=lambda d: int(d["grain_id"]))
# ─── Build time‐0 voxel dictionary ────────────────────────────────────────
def _coord_to_index_1based(c: float, o: float, d: float) -> int:
"""
Convert a coordinate to a 1-based voxel index
This function computes the voxel index corresponding to a coordinate `c`,
given the origin `o` and voxel size `d`. The result is rounded and adjusted
to be 1-based, robust to small floating-point errors.
Parameters
----------
c : float
Coordinate value
o : float
Origin coordinate of the voxel grid
d : float
Size of a voxel along the corresponding axis
Returns
-------
int
1-based voxel index corresponding to the coordinate
"""
# i = round((c - o)/d + 0.5), robust to tiny float noise
return int(round((float(c) - float(o)) / float(d) + 0.5))
def _clamp(v: float, lo: float, hi: float) -> float:
"""
Clamp a value between a lower and upper bound
This function ensures that the input value `v` lies within the range
[`lo`, `hi`]. Values below `lo` are set to `lo`, and values above `hi`
are set to `hi`.
Parameters
----------
v : float
Value to be clamped
lo : float
Lower bound
hi : float
Upper bound
Returns
-------
float
The clamped value within the specified bounds
"""
return max(lo, min(hi, v))
voxels_t0 = []
for vid, gid in voxel_to_grain.items():
cx, cy, cz = vox_center_dict.get(vid, (0.0, 0.0, 0.0))
centroid = [float(cx) * length_scale,
float(cy) * length_scale,
float(cz) * length_scale]
# Compute 1-based voxel indices from centers
ix = _coord_to_index_1based(centroid[0], origin[0], unit_sizes[0])
iy = _coord_to_index_1based(centroid[1], origin[1], unit_sizes[1])
iz = _coord_to_index_1based(centroid[2], origin[2], unit_sizes[2])
# Clamp to valid range [1..N*]
ix = _clamp(ix, 1, Nx)
iy = _clamp(iy, 1, Ny)
iz = _clamp(iz, 1, Nz)
# Phase is stored at grain level in grain_phase_dict.
# Since each voxel belongs to one grain, the voxel phase is inherited
# directly from its parent grain. This keeps phase_id available at both
# grain and voxel level in the exported MiMeDat object.
phase_id = grain_phase_dict.get(gid)
phase_id = int(phase_id) if phase_id is not None else None
entry = {
"voxel_id": int(vid),
"grain_id": int(gid),
"phase_id": phase_id,
"centroid_coordinates": centroid, # scaled for output
"voxel_index": [int(ix), int(iy), int(iz)], # 1-based indices
"voxel_volume": float(voxel_volume),
}
if include_orientation:
ori = grain_ori_dict.get(gid)
if ori is not None:
entry["orientation"] = list(ori)
voxels_t0.append(entry)
voxels_t0_sorted = sorted(voxels_t0, key=lambda d: int(d["voxel_id"]))
# ─── wrap under the time‐step keys ────────────────────────────────────────
microstructure_t0 = {
"time_point": 0,
"grid": grid_t0,
"grains": grains_t0_sorted,
"voxels": voxels_t0_sorted,
}
microstructure_t0["microstructure_state_id"] = self.create_microstructure_identifier(microstructure_t0)
time_steps = [
microstructure_t0,
# …etc…
]
# Assemble final structure with placeholders
data = {
**use, # expand user-specific dict entries directly
"software": "",
"software_version": "",
"system": "",
"system_version": "",
"processor_specifications": "",
"input_path": "",
"results_path": "",
**ig, # expand initial geometry dict entries directly
"global_temperature": 298,
**job_bc, # expand boundary condition dict entries directly
"phase": phase_list,
"microstructure": time_steps # Time-level storage: time-frame data of voxels and grains
}
return data
[docs]
def pckl(
self,
file: Optional[Union[str, os.PathLike[str]]] = None,
path: Union[str, os.PathLike[str]] = './') -> None:
"""
Write the microstructure into a pickle file
This function serializes the microstructure object using pickle, allowing
complex structures to be stored and later reloaded.
Parameters
----------
file : str or os.PathLike or None, optional, default=None
File name for the pickled microstructure. If None, the filename defaults
to the microstructure name plus '.pckl'.
path : str or os.PathLike, optional, default='./'
Directory path to save the pickle file. Default is current directory './'.
Returns
-------
None
"""
import pickle
if file is None:
if self.name == 'Microstructure':
file = 'px_{}grains_microstructure.pckl'.format(self.Ngr)
else:
file = self.name + '_microstructure.pckl'
path = os.path.normpath(path)
file = os.path.join(path, file)
with open(file, 'wb') as output:
pickle.dump(self, output, pickle.HIGHEST_PROTOCOL)
return
[docs]
def import_particles(
self,
file: Union[str, os.PathLike[str]],
path: Union[str, os.PathLike[str]] = './') -> None:
"""
Import particle data from a dump file
This function reads particle information from a specified dump file and
stores it in the microstructure object, updating the simulation box and
particle list.
Parameters
----------
file : str or os.PathLike
Name of the input dump file containing particle data
path : str or os.PathLike, optional, default='./'
Directory path where the dump file is located. Default is current directory './'
Returns
-------
None
Updates the object's simbox and particles attributes
"""
path = os.path.normpath(path)
file = os.path.join(path, file)
self.simbox, self.particles = read_dump(file)
"""
-------- legacy methods --------
"""
[docs]
def init_stats(
self,
descriptor: Any = None,
gs_data: Any = None,
ar_data: Any = None,
porous: bool = False,
save_files: bool = False) -> None:
"""
Legacy function that redirects to plot_stats_init
This function is deprecated and serves as a wrapper for `plot_stats_init`.
A warning is logged advising users to use `plot_stats_init` directly.
Parameters
----------
descriptor : Any, optional, default=None
Descriptor data passed to `plot_stats_init`. Default is None.
gs_data : Any, optional, default=None
Grain size data passed to `plot_stats_init`. Default is None.
ar_data : Any, optional, default=None
Aspect ratio data passed to `plot_stats_init`. Default is None.
porous : bool, optional, default=False
Deprecated parameter, ignored. Default is False.
save_files : bool, optional, default=False
Whether to save generated statistics files. Default is False.
Returns
-------
None
Calls `plot_stats_init` and does not return any value
"""
logging.warning('"init_stats" is a legacy function and will be depracted, please use "plot_stats_init()".')
self.plot_stats_init(descriptor, gs_data=gs_data, ar_data=ar_data, save_files=save_files)
[docs]
def output_abq(
self,
nodes: Any = None,
name: Optional[Union[str, os.PathLike[str]]] = None,
voxel_dict: Optional[Mapping[Any, Any]] = None,
grain_dict: Optional[Mapping[Any, Any]] = None,
faces: Any = None,
dual_phase: bool = False,
thermal: bool = False,
units: Optional[str] = None) -> None:
"""
Legacy function that redirects to write_abq
This function is deprecated and serves as a wrapper for `write_abq`.
Warnings are logged advising users to use `write_abq` directly. The `faces`
parameter is determined automatically if provided.
Parameters
----------
nodes : array-like or None, optional, default=None
Nodal coordinates to include in the output. Default is None.
name : str or os.PathLike or None, optional, default=None
File name for the Abaqus input file. Default is None.
voxel_dict : dict or None, optional, default=None
Voxel connectivity information. Default is None.
grain_dict : dict or None, optional, default=None
Grain-to-voxel mapping. Default is None.
faces : Any, optional, default=None
Deprecated parameter, automatically determined. Default is None.
dual_phase : bool, optional, default=False
Whether to include dual-phase data. Default is False.
thermal : bool, optional, default=False
Whether to include thermal data. Default is False.
units : str or None, optional, default=None
Units for the Abaqus model.
Returns
-------
None
Calls `write_abq` and does not return any value
"""
logging.warning('"output_abq" is a legacy function and will be depracted, please use "write_abq()".')
if faces is not None:
logging.warning('Parameter "faces" will be determined automatically.')
self.write_abq(nodes=nodes, file=name, voxel_dict=voxel_dict, grain_dict=grain_dict,
dual_phase=dual_phase, thermal=thermal, units=units)