From 2c97a1b3aaf08bce2677ba9d096d6d36b29a7f16 Mon Sep 17 00:00:00 2001 From: AHartmaier Date: Fri, 29 Mar 2024 17:06:01 +0100 Subject: [PATCH] Updates in RVE generation from EBSD maps --- src/kanapy/__init__.py | 7 + src/kanapy/api.py | 26 ++-- src/kanapy/entities.py | 13 +- src/kanapy/initializations.py | 100 +++++++++---- src/kanapy/packing.py | 55 +++---- src/kanapy/plotting.py | 72 ++++++++- src/kanapy/rve_stats.py | 273 +++++++++++++++++++++------------- src/kanapy/textures.py | 35 +++-- src/kanapy/voxelization.py | 3 +- tests/test_packing.py | 3 +- tests/test_texture.py | 8 +- tests/test_voxelization.py | 2 +- 12 files changed, 392 insertions(+), 205 deletions(-) diff --git a/src/kanapy/__init__.py b/src/kanapy/__init__.py index f993aab7..ec4b9e7d 100644 --- a/src/kanapy/__init__.py +++ b/src/kanapy/__init__.py @@ -8,6 +8,7 @@ from kanapy.plotting import plot_voxels_3D, plot_polygons_3D from kanapy.input_output import writeAbaqusMat, pickle2microstructure, import_voxels, \ import_stats, write_stats +from kanapy.rve_stats import find_rot_axis from kanapy.util import ROOT_DIR, MAIN_DIR, MTEX_DIR, poly_scale log_level = 20 # Levels for logging: 10: DEBUG, 20: INFO, 30: WARNING, 40: ERROR @@ -21,6 +22,12 @@ logging.error('Inconsistent installation of Kanapy package. MTEX tools will not be available.') MTEX_AVAIL = False +try: + from kanapy.triple_surface import create_ref_ell + triple_surf = True +except: + triple_surf = False + __author__ = 'Mahesh R.G Prasad, Abhishek Biswas, Golsa Tolooei Eshlaghi, Alexander Hartmaier' __email__ = 'alexander.hartmaier@rub.de' __version__ = version('kanapy') diff --git a/src/kanapy/api.py b/src/kanapy/api.py index 17592a81..cc32f2d4 100644 --- a/src/kanapy/api.py +++ b/src/kanapy/api.py @@ -70,6 +70,8 @@ class Microstructure(object): rve_stats : list List of dictionaries containing statistical information on different RVE types: particles, voxels, polyhedral grains + rve_stats_labels : list + List of strings containing the labels for the RVE types, i.e. Partcls, Voxels, Grains """ def __init__(self, descriptor=None, file=None, name='Microstructure'): @@ -84,6 +86,7 @@ def __init__(self, descriptor=None, file=None, name='Microstructure'): self.simbox = None self.mesh = None self.rve_stats = None + self.rve_stats_labels = None self.from_voxels = False self.ialloy = None @@ -451,10 +454,11 @@ def plot_stats(self, data=None, labels = [] if dual_phase: iphase = ip - print(f'Plotting input & output statistics for phase {ip}') - if data is None or \ - (type(data) is str and 'p' in data.lower()): - # analyse particle statistics + 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.') @@ -465,9 +469,9 @@ def plot_stats(self, data=None, stats_list.append(part_stats) labels.append("Partcls") - if data is None or \ - (type(data) is str and 'v' in data.lower()): - # analyse voxel statistics + # 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.') @@ -477,9 +481,10 @@ def plot_stats(self, data=None, verbose=verbose, save_files=save_files) stats_list.append(vox_stats) labels.append('Voxels') - if data is None or \ - (type(data) is str and 'g' in data.lower()): - # analyse voxel statistics + + # 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.') @@ -489,6 +494,7 @@ def plot_stats(self, data=None, verbose=verbose, save_files=save_files) stats_list.append(grain_stats) labels.append('Grains') + if dual_phase: print(f'\nStatistical microstructure parameters of phase {iphase} in RVE') print('-------------------------------------------------------') diff --git a/src/kanapy/entities.py b/src/kanapy/entities.py index e216f37c..490db93f 100644 --- a/src/kanapy/entities.py +++ b/src/kanapy/entities.py @@ -99,7 +99,7 @@ def __init__(self, iden, x, y, z, a, b, c, quat, phasenum=0, dup=None, points=No self.speedx = 0. self.speedy = 0. self.speedz = 0. - self.rotationMatrixGen() # Initialize roatation matrix for the ellipsoid + self.rotation_matrix = self.rotationMatrixGen() # Initialize roatation matrix for the ellipsoid self.surface_points = self.surfacePointsGen() # Initialize surface points for the ellipsoid self.inside_voxels = [] # List that stores voxels belonging to the ellipsoid self.set_cub() # sets particle cuboid for collision testing with octree boxes @@ -151,6 +151,8 @@ def rotationMatrixGen(self): w, x, y, z = self.quat Nq = w * w + x * x + y * y + z * z + if Nq < FLOAT_EPS: + return np.eye(3) s = 2.0 / Nq X = x * s @@ -166,12 +168,9 @@ def rotationMatrixGen(self): yZ = y * Z zZ = z * Z - if Nq < FLOAT_EPS: - self.rotation_matrix = np.eye(3) - else: - self.rotation_matrix = np.array([[1.0 - (yY + zZ), xY - wZ, xZ + wY], - [xY + wZ, 1.0 - (xX + zZ), yZ - wX], - [xZ - wY, yZ + wX, 1.0 - (xX + yY)]]) + return np.array([[1.0 - (yY + zZ), xY - wZ, xZ + wY], + [xY + wZ, 1.0 - (xX + zZ), yZ - wX], + [xZ - wY, yZ + wX, 1.0 - (xX + yY)]]) def surfacePointsGen(self, nang=20): """ diff --git a/src/kanapy/initializations.py b/src/kanapy/initializations.py index 02fc5f15..0e281692 100644 --- a/src/kanapy/initializations.py +++ b/src/kanapy/initializations.py @@ -1,3 +1,11 @@ +""" +Classes for RVE and mesh initialization + +Authors: Mahesh Prassad, Abhishek Biswas, Alexander Hartmaier +Ruhr University Bochum, Germany + +March 2024 +""" import json import os import numpy as np @@ -80,12 +88,12 @@ def __init__(self, stats_dicts, nsteps=1000, from_voxels=False, poly=None): Attributes ---------- - self.packing_steps = nsteps # max number of steps for packing simulation - self.size = None # tuple of lengths along Cartesian axes - self.dim = None # tuple of number of voxels along Cartesian axes - self.periodic = None # Boolean for periodicity of RVE - self.units = None # Units of RVE dimensions, either "mm" or "um" (micron) - self.nparticles = [] # list of particle numbers for each phase + packing_steps = nsteps # max number of steps for packing simulation + size = None # tuple of lengths along Cartesian axes + dim = None # tuple of number of voxels along Cartesian axes + periodic = None # Boolean for periodicity of RVE + units = None # Units of RVE dimensions, either "mm" or "um" (micron) + nparticles = [] # list of particle numbers for each phase particle_data = [] # list of cits for statistical particle data for each grains phase_names = [] # list of names of phases phase_vf = [] # list of volume fractions of phases @@ -93,18 +101,20 @@ def __init__(self, stats_dicts, nsteps=1000, from_voxels=False, poly=None): Number of alloy in ICAMS CP-UMAT """ - def init_particles(particle_data): + def init_particles(): """ Extract statistical microstructure information from data dictionary and initalize particles for packing accordingly """ def gen_data_basic(pdict): - # Compute the Log-normal PDF & CDF. - frozen_lognorm = lognorm(s=sig_ED, loc=loc_ED, scale=scale_ED) - + """Computes number of particles according to equivalent diameter distribution + and initializes particle ensemble with individual diameters approximating + this lognormal distribution. + """ + # Compute the Log-normal CDF according to parameters for equiv. diameter distrib. xaxis = np.linspace(0.1, 200, 1000) - ycdf = frozen_lognorm.cdf(xaxis) + ycdf = lognorm.cdf(xaxis, s=sig_ED, loc=loc_ED, scale=scale_ED) # Get the mean value for each pair of neighboring points as centers of bins xaxis = np.vstack([xaxis[1:], xaxis[:-1]]).mean(axis=0) @@ -145,7 +155,7 @@ def gen_data_basic(pdict): totalEllipsoids = int(np.sum(num)) # Duplicate the diameter values - eq_Dia = np.repeat(eq_Dia, num) # better calculate num first + eq_Dia = np.repeat(eq_Dia, num) # Raise value error in case the RVE side length is too small to fit grains inside. if len(eq_Dia) == 0: @@ -189,24 +199,34 @@ def gen_data_elong(pdict): # Sample from Normal distribution: It takes mean and std of normal distribution tilt_angle = [] num = pdict['Number'] + iter = 0 while num > 0: tilt = vonmises.rvs(kappa, loc=loc_ori, size=num) index_array = np.where((tilt >= ori_cutoff_min) & (tilt <= ori_cutoff_max)) TA = tilt[index_array].tolist() tilt_angle.extend(TA) num = pdict['Number'] - len(tilt_angle) + iter += 1 + if iter > 10000: + raise StopIteration('Iteration for tilt angles did not converge in 10000 iterations.' + 'Increase cutoff range to assure proper generation of particles.') # Aspect ratio statistics # Sample from lognormal or gamma distribution: # it takes mean, std and scale of the underlying normal distribution finalAR = [] num = pdict['Number'] + iter = 0 while num > 0: ar = lognorm.rvs(sig_AR, loc=loc_AR, scale=scale_AR, size=num) index_array = np.where((ar >= ar_cutoff_min) & (ar <= ar_cutoff_max)) AR = ar[index_array].tolist() finalAR.extend(AR) num = pdict['Number'] - len(finalAR) + iter += 1 + if iter > 10000: + raise StopIteration('Iteration for aspect ratios did not converge in 10000 iterations.' + 'Increase cutoff range to assure proper generation of particles.') finalAR = np.array(finalAR) # Calculate the major, minor axes lengths for particles using: @@ -222,10 +242,9 @@ def gen_data_elong(pdict): pdict['Tilt angle'] = list(tilt_angle) return pdict - if stats["Grain type"] not in ["Elongated", "Equiaxed"]: - raise ValueError('The value for "Grain type" must be either "Equiaxed" or "Elongated".') - - # Attributes for equivalent diameter + # start of particle initialization + # analyze information on grains in "stats" dictionary from outer scope + # 1. Attributes for equivalent diameter if 'mean' in stats['Equivalent diameter'].keys(): # load legacy names to ensure compatibility with previous versions sig, loc, scale, kappa = stat_names(legacy=True) @@ -235,22 +254,28 @@ def gen_data_elong(pdict): sig, loc, scale, kappa = stat_names(legacy=False) sig_ED = stats["Equivalent diameter"][sig] scale_ED = stats["Equivalent diameter"][scale] - loc_ED = stats["Equivalent diameter"][loc] + if loc in stats["Equivalent diameter"].keys(): + loc_ED = stats["Equivalent diameter"][loc] + else: + loc_ED = 0.0 dia_cutoff_min = stats["Equivalent diameter"]["cutoff_min"] dia_cutoff_max = stats["Equivalent diameter"]["cutoff_max"] if dia_cutoff_min / dia_cutoff_max > 0.75: raise ValueError('Min/Max values for cutoffs of equiavalent diameter are too close: ' + f'Max: {dia_cutoff_max}, Min: {dia_cutoff_min}') - # generate dict for particle data + # generate dict for basic particle data: number of particles and equiv. diameters pdict = gen_data_basic(dict({'Type': stats["Grain type"], 'Phase': ip})) - # Additional attributes for elongated grains + # 2. Additional attributes for elongated or freely-defined grains if stats["Grain type"] == "Elongated": # Extract descriptors for aspect ratio distrib. from dict sig_AR = stats["Aspect ratio"][sig] scale_AR = stats["Aspect ratio"][scale] - loc_AR = stats["Aspect ratio"][loc] + if loc in stats["Aspect ratio"].keys(): + loc_AR = stats["Aspect ratio"][loc] + else: + loc_AR = 0.0 ar_cutoff_min = stats["Aspect ratio"]["cutoff_min"] ar_cutoff_max = stats["Aspect ratio"]["cutoff_max"] if ar_cutoff_min / ar_cutoff_max > 0.75: @@ -265,29 +290,40 @@ def gen_data_elong(pdict): if ori_cutoff_min / ori_cutoff_max > 0.75: raise ValueError('Min/Max values for cutoffs of orientation of tilt axis are too close: ' + f'Max: {ori_cutoff_max}, Min: {ori_cutoff_min}') - # Add attributes for elongated particle to dictionary pdict = gen_data_elong(pdict) - particle_data.append(pdict) - self.nparticles.append(pdict['Number']) - return particle_data + elif stats["Grain type"] == "Free": + try: + from kanapy.triple_surface import gen_data_free + except: + raise ModuleNotFoundError('Free grain definitions are not available ' + 'in the public version of Kanapy.') + # Add attributes for elongated particle to dictionary + pdict = gen_data_free(pdict, stats) + return pdict + return pdict + # Start RVE generation if from_voxels: print('Creating an RVE from voxel input') else: print('Creating an RVE based on user defined statistics') - # Extract grain diameter statistics info + # declare attributes of RVE object self.packing_steps = nsteps # max number of steps for packing simulation self.size = None # tuple of lengths along Cartesian axes self.dim = None # tuple of number of voxels along Cartesian axes self.periodic = None # Boolean for periodicity of RVE self.units = None # Units of RVE dimensions, either "mm" or "um" (micron) self.ialloy = None # Number of alloy in ICAMS CP-UMAT - self.nparticles = [] # List of article numbers for each phase - particle_data = [] # list of cits for statistical particle data for each grains phase_names = [] # list of names of phases phase_vf = [] # list of volume fractions of phases ialloy = [] + if from_voxels: + self.particle_data = None + self.nparticles = None + else: + self.particle_data = [] # list of cits for statistical particle data for each grains + self.nparticles = [] # List of article numbers for each phase # extract data from descriptors of individual phases for ip, stats in enumerate(stats_dicts): @@ -351,10 +387,11 @@ def gen_data_elong(pdict): # Extract grains shape attributes to initialize particles if not from_voxels: - particle_data = init_particles(particle_data) - else: - particle_data = None - self.nparticles = None + if stats["Grain type"] not in ["Elongated", "Equiaxed", "Free"]: + raise ValueError('The value for "Grain type" must be either "Equiaxed" or "Elongated".') + part_dict = init_particles() + self.particle_data.append(part_dict) + self.nparticles.append(part_dict['Number']) print(' RVE characteristics:') print(f' RVE side lengths (X, Y, Z) = {self.size} ({self.units})') print(f' Number of voxels (X, Y, Z) = {self.dim}') @@ -367,7 +404,6 @@ def gen_data_elong(pdict): print('\n') self.phase_names = phase_names self.phase_vf = phase_vf - self.particle_data = particle_data nall = len(ialloy) if nall > 0: if nall != len(phase_vf): diff --git a/src/kanapy/packing.py b/src/kanapy/packing.py index 56dc7c33..89664411 100644 --- a/src/kanapy/packing.py +++ b/src/kanapy/packing.py @@ -1,14 +1,13 @@ # -*- coding: utf-8 -*- import random -from tqdm import tqdm import numpy as np - +from tqdm import tqdm from kanapy.input_output import write_dump from kanapy.entities import Ellipsoid, Cuboid, Octree -from kanapy.collisions import collision_routine, collide_detect +from kanapy.collisions import collision_routine -def particle_generator(particle_data, sim_box, periodic, poly): +def particle_generator(particle_data, sim_box, poly): """ Initializes ellipsoids by assigning them random positions and speeds within the simulation box. @@ -18,10 +17,8 @@ def particle_generator(particle_data, sim_box, periodic, poly): :type particle_data: Python dictionary :param sim_box: Simulation box representing RVE. :type sim_box: :class:`entities.Simulation_Box` - :param periodic: Indicates periodicity of RVE. - :type periodic: :class: bool :param poly: Points defining primitive polygon inside ellipsoid (optional). - :type periodic: :class: ndarry(N, 3) + :type poly: :class: ndarry(N, 3) :returns: Ellipsoids for the packing routine :rtype: list @@ -29,8 +26,8 @@ def particle_generator(particle_data, sim_box, periodic, poly): Ellipsoids = [] id_ctr = 0 vell = 0. - # introduce scaling factor to reduce particle overlap for non-periodic box - sf = 0.5 # if periodic else 0.45 + # introduce scaling factor to reduce particle overlap + sf = 0.5 for particle in particle_data: num_particles = particle['Number'] # Total number of ellipsoids for n in range(num_particles): @@ -47,24 +44,28 @@ def particle_generator(particle_data, sim_box, periodic, poly): y = random.uniform(b, sim_box.h - b) z = random.uniform(c, sim_box.d - c) - # Angle represents inclination of Major axis w.r.t positive x-axis - if particle['Type'] == 'Equiaxed': - angle = 0.5 * np.pi + if particle['Type'] == 'Free': + # For free particle definition, quaternion for rotation is given + quat = particle['Quaternion'][n] else: - # Extract the angle - angle = particle['Tilt angle'][n] - # Tilt vector wrt (+ve) x-axis - vec_a = np.array([a * np.cos(angle), a * np.sin(angle), 0.0]) - # Do the cross product to find the quaternion axis - cross_a = np.cross(np.array([1, 0, 0]), vec_a) - # norm of the vector (Magnitude) - norm_cross_a = np.linalg.norm(cross_a, 2) - quat_axis = cross_a / norm_cross_a # normalize the quaternion axis - - # Find the quaternion components - qx, qy, qz = quat_axis * np.sin(angle / 2) - qw = np.cos(angle / 2) - quat = np.array([qw, qx, qy, qz]) + # Angle represents inclination of Major axis w.r.t positive x-axis in xy-plane + if particle['Type'] == 'Equiaxed': + angle = 0.5 * np.pi + else: + # Extract the angle + angle = particle['Tilt angle'][n] + # Tilt vector wrt (+ve) x-axis + vec_a = np.array([a * np.cos(angle), a * np.sin(angle), 0.0]) + # Do the cross product to find the quaternion axis + cross_a = np.cross(np.array([1, 0, 0]), vec_a) + # norm of the vector (Magnitude) + norm_cross_a = np.linalg.norm(cross_a, 2) + quat_axis = cross_a / norm_cross_a # normalize the quaternion axis + + # Find the quaternion components + qx, qy, qz = quat_axis * np.sin(angle / 2) + qw = np.cos(angle / 2) + quat = np.array([qw, qx, qy, qz]) # instance of Ellipsoid class ellipsoid = Ellipsoid(iden, x, y, z, a, b, c, quat, @@ -322,7 +323,7 @@ def packingRoutine(particle_data, periodic, nsteps, sim_box, print(' Creating particles from distribution statistics') # Create instances for particles - Particles = particle_generator(particle_data, sim_box, periodic, poly) + Particles = particle_generator(particle_data, sim_box, poly) # Growth of particle at each time step print(' Particle packing by growth simulation') diff --git a/src/kanapy/plotting.py b/src/kanapy/plotting.py index 582763a5..b1e1adc4 100644 --- a/src/kanapy/plotting.py +++ b/src/kanapy/plotting.py @@ -451,7 +451,10 @@ def plot_init_stats(stats_dict, gs_data=None, ar_data=None, save_files=False): # Extract grain diameter statistics info from input file sd = stats_dict["Equivalent diameter"]["sig"] scale = stats_dict["Equivalent diameter"]["scale"] - loc = stats_dict["Equivalent diameter"]["loc"] + if 'loc' in stats_dict["Equivalent diameter"].keys(): + loc = stats_dict["Equivalent diameter"]["loc"] + else: + loc = 0. dia_cutoff_min = stats_dict["Equivalent diameter"]["cutoff_min"] dia_cutoff_max = stats_dict["Equivalent diameter"]["cutoff_max"] @@ -492,11 +495,14 @@ def plot_init_stats(stats_dict, gs_data=None, ar_data=None, save_files=False): plt.title("Grain equivalent diameter distribution", fontsize=20) plt.legend(fontsize=16) plt.show() - elif stats_dict["Grain type"] == "Elongated": + elif stats_dict["Grain type"] in ["Elongated"]: # Extract grain aspect ratio descriptors from input file sd_AR = stats_dict["Aspect ratio"]["sig"] scale_AR = stats_dict["Aspect ratio"]["scale"] - loc_AR = stats_dict["Aspect ratio"]["loc"] + if 'loc' in stats_dict["Aspect ratio"].keys(): + loc_AR = stats_dict["Aspect ratio"]["loc"] + else: + loc_AR = 0.0 ar_cutoff_min = stats_dict["Aspect ratio"]["cutoff_min"] ar_cutoff_max = stats_dict["Aspect ratio"]["cutoff_max"] @@ -537,6 +543,66 @@ def plot_init_stats(stats_dict, gs_data=None, ar_data=None, save_files=False): ax[1].hist(ar_data, bins=15, density=True, label='experimental data') ax[1].legend(fontsize=14) plt.show() + elif stats_dict["Grain type"] in ["Free"]: + # Compute reference grain aspect ratio descriptors from stats dict + # Identify most likely rotation axis of ellepsoid with free semi-axes + from kanapy.rve_stats import find_rot_axis + semi_axs = stats_dict["Semi axes"]["scale"] + cut_min = stats_dict["Semi axes"]["cutoff_min"] + cut_max = stats_dict["Semi axes"]["cutoff_max"] + irot = find_rot_axis(semi_axs[0], semi_axs[1], semi_axs[2]) + if irot == 0: + ar_scale = 2.0 * semi_axs[0] / (semi_axs[1] + semi_axs[2]) + ar_cutoff_min = 2.0 * cut_min[0] / (cut_min[1] + cut_min[2]) + ar_cutoff_max = 2.0 * cut_max[0] / (cut_max[1] + cut_max[2]) + elif irot == 1: + ar_scale = 2.0 * semi_axs[1] / (semi_axs[0] + semi_axs[2]) + ar_cutoff_min = 2.0 * cut_min[1] / (cut_min[0] + cut_min[2]) + ar_cutoff_max = 2.0 * cut_max[1] / (cut_max[0] + cut_max[2]) + elif irot == 2: + ar_scale = 2.0 * semi_axs[2] / (semi_axs[1] + semi_axs[0]) + ar_cutoff_min = 2.0 * cut_min[2] / (cut_min[1] + cut_min[0]) + ar_cutoff_max = 2.0 * cut_max[2] / (cut_max[1] + cut_max[0]) + else: + raise ValueError('Rotation axis not identified properly') + ar_sig = stats_dict["Semi axes"]["sig"][irot] + + sns.set(color_codes=True) + fig, ax = plt.subplots(2, 1, figsize=(9, 9)) + plt.ion() + + # Plot grain size distribution + ax[0].plot(xaxis, ypdf, linestyle='-', linewidth=3.0) + ax[0].fill_between(xaxis, 0, ypdf, alpha=0.3) + ax[0].set_xlim(left=0.0, right=x_lim) + ax[0].set_xlabel('Equivalent diameter (μm)', fontsize=18) + ax[0].set_ylabel('Density', fontsize=18) + ax[0].tick_params(labelsize=14) + ax[0].axvline(dia_cutoff_min, linestyle='--', linewidth=3.0, + label='Minimum cut-off: {:.2f}'.format(dia_cutoff_min)) + ax[0].axvline(dia_cutoff_max, linestyle='-', linewidth=3.0, + label='Maximum cut-off: {:.2f}'.format(dia_cutoff_max)) + if gs_data is not None: + ax[0].hist(gs_data, bins=80, density=True, label='experimental data') + ax[0].legend(fontsize=14) + + # Plot aspect ratio statistics + # Compute the Log-normal PDF & CDF. + xaxis = np.linspace(0.5 * ar_cutoff_min, 2 * ar_cutoff_max, 500) + ypdf = lognorm.pdf(xaxis, ar_sig, scale=ar_scale) + ax[1].plot(xaxis, ypdf, linestyle='-', linewidth=3.0) + ax[1].fill_between(xaxis, 0, ypdf, alpha=0.3) + ax[1].set_xlabel('Aspect ratio', fontsize=18) + ax[1].set_ylabel('Density', fontsize=18) + ax[1].tick_params(labelsize=14) + ax[1].axvline(ar_cutoff_min, linestyle='--', linewidth=3.0, + label='Minimum cut-off: {:.2f}'.format(ar_cutoff_min)) + ax[1].axvline(ar_cutoff_max, linestyle='-', linewidth=3.0, + label='Maximum cut-off: {:.2f}'.format(ar_cutoff_max)) + if ar_data is not None: + ax[1].hist(ar_data, bins=15, density=True, label='experimental data') + ax[1].legend(fontsize=14) + plt.show() else: raise ValueError('Value for "Grain_type" must be either "Equiaxed" or "Elongated".') diff --git a/src/kanapy/rve_stats.py b/src/kanapy/rve_stats.py index 8e54e5f7..0d6e62f2 100644 --- a/src/kanapy/rve_stats.py +++ b/src/kanapy/rve_stats.py @@ -109,6 +109,72 @@ def pts_in_ellips(Mcomp, pts): return score / len(pts) +def get_diameter(pts): + """ + Get estimate of largest diameter of a set of points. + Parameters + ---------- + pts : (N, dim) ndarray + Point set in dim dimensions + + Returns + ------- + diameter : (dim)-ndarray + + """ + ind0 = np.argmin(pts, axis=0) # index of point with lowest coordinate for each Cartesian axis + ind1 = np.argmax(pts, axis=0) # index of point with highest coordinate for each Cartesian axis + v_min = np.array([pts[i, j] for j, i in enumerate(ind0)]) # min. value for each Cartesian axis + v_max = np.array([pts[i, j] for j, i in enumerate(ind1)]) # max. value for each Cartesian axis + ind_d = np.argmax(v_max - v_min) # Cartesian axis along which largest distance occurs + return pts[ind1[ind_d], :] - pts[ind0[ind_d], :] + +def project_pts(pts, ctr, axis): + """ + Project points (pts) to plane defined via center (ctr) and normal vector (axis). + + Parameters + ---------- + pts : (N, dim) ndarray + Point set in dim dimensions + ctr : (dim)-ndarray + Center point of the projection plane + axis : (dim)-ndarray + Unit vector for plane normal + + Returns + ------- + ppt : (N, dim) ndarray + Points projected to plane with normal axis + """ + dvec = pts - ctr[None, :] # distance vector b/w points and center point + pdist = np.array([np.dot(axis, v) for v in dvec]) + ppt = np.zeros(pts.shape) + for i, p in enumerate(dvec): + ppt[i, :] = p - pdist[i] * axis + return ppt + + +def bbox(pts): + # Approximate smallest rectangular cuboid around points of grains + # to analyse prolate (aspect ratio > 1) and oblate (a.r. < 1) particles correctly + dia = get_diameter(pts) # approx. of largest diameter of grain + ctr = np.mean(pts, axis=0) # approx. center of grain + len_a = np.linalg.norm(dia) # length of largest side + if len_a < 1.e-5: + logging.warning(f'Very small grain at {ctr} with max. diameter = {len_a}') + ax_a = dia / len_a # unit vector along longest side + ppt = project_pts(pts, ctr, ax_a) # project points onto central plane normal to diameter + trans1 = get_diameter(ppt) # largest side transversal to long axis + len_b = np.linalg.norm(trans1) # length of second-largest side + ax_b = trans1 / len_b # unit vector of second axis (normal to diameter) + ax_c = np.cross(ax_a, ax_b) # calculate third orthogonal axes of rectangular cuboid + lpt = project_pts(ppt, np.zeros(3), ax_b) # project points on third axis + pdist = np.array([np.dot(ax_c, v) for v in lpt]) # calculate distance of points on third axis + len_c = np.max(pdist) - np.min(pdist) # get length of shortest side + return 0.5*len_a, 0.5*len_b, 0.5*len_c + + def calc_stats_dict(a, b, c, eqd): arr_a = np.sort(a) arr_b = np.sort(b) @@ -154,6 +220,102 @@ def calc_stats_dict(a, b, c, eqd): return sd +def get_stats_part(part, iphase=None, ax_max=None, + minval=1.e-5, show_plot=True, + verbose=False, save_files=False): + """ Extract statistics about the microstructure from particles. If inner structure is contained + by fitting a 3D ellipsoid to each structure. + + Parameters + ---------- + part + minval + show_plot + verbose + ax_max + + Returns + ------- + + """ + if ax_max is not None: + minval = max(minval, 1. / ax_max ** 2) + cons = ({'type': 'ineq', 'fun': con_fun}) # constraints for minimization + opt = {'maxiter': 200} # options for minimization + mc = np.array([1., 1., 1., 0., 0., 0.]) # start value of matrix for minimization + arr_a = [] + arr_b = [] + arr_c = [] + arr_eqd = [] + for pc in part: + # decide if phase-specific analysis is performed + if iphase is not None and iphase != pc.phasenum: + continue + if pc.inner is not None: + pts = pc.inner.points + rdict = minimize(pts_in_ellips, x0=mc, args=(pts,), method='SLSQP', + constraints=cons, options=opt) + if not rdict['success']: + if verbose: + print(f'Optimization failed for particle {pc.id}') + print(rdict['message']) + continue + eval, evec = np.linalg.eig(arr2mat(rdict['x'])) + if any(eval <= minval): + if verbose: + print(f'Matrix for particle {pc.id} not positive definite or semi-axes too large. ' + f'Eigenvalues: {eval}') + continue + if verbose: + print(f'Optimization succeeded for particle {pc.id} after {rdict["nit"]} iterations.') + print(f'Eigenvalues: {eval}') + print(f'Eigenvectors: {evec}') + # Semi-axes of ellipsoid + ea = 1. / np.sqrt(eval[0]) + eb = 1. / np.sqrt(eval[1]) + ec = 1. / np.sqrt(eval[2]) + eqd = 2.0 * (ea * eb * ec) ** (1.0 / 3.0) + arr_a.append(ea) + arr_b.append(eb) + arr_c.append(ec) + arr_eqd.append(eqd) + else: + eqd = 2.0 * (pc.a * pc.b * pc.c) ** (1.0 / 3.0) + arr_a.append(pc.a) + arr_b.append(pc.b) + arr_c.append(pc.c) + arr_eqd.append(eqd) + + # calculate statistical parameters + part_stats_dict = calc_stats_dict(arr_a, arr_b, arr_c, arr_eqd) + if verbose: + print('\n--------------------------------------------------') + print('Statistical microstructure parameters of particles') + print('--------------------------------------------------') + print('Median lengths of semi-axes of fitted ellipsoids in micron') + print(f'a: {part_stats_dict["a_scale"]:.3f}, b: {part_stats_dict["b_scale"]:.3f}, ' + f'c: {part_stats_dict["c_scale"]:.3f}') + av_std = np.mean([part_stats_dict['a_sig'], part_stats_dict['b_sig'], part_stats_dict['c_sig']]) + print(f'Average standard deviation of semi-axes: {av_std:.4f}') + print('\nAssuming rotational symmetry in grains') + print(f'Rotational axis: {part_stats_dict["ind_rot"]}') + print(f'Median aspect ratio: {part_stats_dict["ar_scale"]:.3f}') + print('\nGrain size') + print(f'Median equivalent grain diameter: {part_stats_dict["eqd_scale"]:.3f} micron') + print(f'Standard deviation of equivalent grain diameter: {part_stats_dict["eqd_sig"]:.4f}') + print('--------------------------------------------------------') + if show_plot: + if part[0].inner is None: + title = 'Particle statistics' + else: + title = 'Statistics of inner particle structures' + if iphase is not None: + title += f' (phase {iphase})' + plot_stats_dict(part_stats_dict, title=title, save_files=save_files) + + return part_stats_dict + + def get_stats_vox(mesh, iphase=None, ax_max=None, minval=1.e-5, show_plot=True, verbose=False, save_files=False): @@ -195,6 +357,9 @@ def get_stats_vox(mesh, iphase=None, ax_max=None, hull = ConvexHull(pts_all) eqd = 2.0 * (gfac * hull.volume) ** (1.0 / 3.0) pts = hull.points[hull.vertices] # outer nodes of grain no. igr + # find bounding box to hull points + ea, eb, ec = bbox(pts) + """ # find best fitting ellipsoid to points rdict = minimize(pts_in_ellips, x0=mc, args=(pts,), method='SLSQP', constraints=cons, options=opt) @@ -216,13 +381,9 @@ def get_stats_vox(mesh, iphase=None, ax_max=None, # Semi-axes of ellipsoid ea = 1. / np.sqrt(eval[0]) eb = 1. / np.sqrt(eval[1]) - ec = 1. / np.sqrt(eval[2]) - arr_a.append(ea) - arr_b.append(eb) - arr_c.append(ec) - arr_eqd.append(eqd) + ec = 1. / np.sqrt(eval[2])""" - """ Plot points on hull with fitted ellipsoid -- only for debugging + """# Plot points on hull with fitted ellipsoid -- only for debugging import matplotlib.pyplot as plt # Points on the outer surface of optimal ellipsoid nang = 100 @@ -248,6 +409,10 @@ def get_stats_vox(mesh, iphase=None, ax_max=None, ax.plot_surface(x, y, z, rstride=4, cstride=4, color=col, linewidth=0) ax.scatter(pts[:, 0], pts[:, 1], pts[:, 2], marker='o') plt.show()""" + arr_a.append(ea) + arr_b.append(eb) + arr_c.append(ec) + arr_eqd.append(eqd) # calculate and print statistical parameters vox_stats_dict = calc_stats_dict(arr_a, arr_b, arr_c, arr_eqd) if verbose: @@ -364,99 +529,3 @@ def get_stats_poly(grains, iphase=None, ax_max=None, plot_stats_dict(poly_stats_dict, title=title, save_files=save_files) return poly_stats_dict - - -def get_stats_part(part, iphase=None, ax_max=None, - minval=1.e-5, show_plot=True, - verbose=False, save_files=False): - """ Extract statistics about the microstructure from particles. If inner structure is contained - by fitting a 3D ellipsoid to each structure. - - Parameters - ---------- - part - minval - show_plot - verbose - ax_max - - Returns - ------- - - """ - if ax_max is not None: - minval = max(minval, 1. / ax_max ** 2) - cons = ({'type': 'ineq', 'fun': con_fun}) # constraints for minimization - opt = {'maxiter': 200} # options for minimization - mc = np.array([1., 1., 1., 0., 0., 0.]) # start value of matrix for minimization - arr_a = [] - arr_b = [] - arr_c = [] - arr_eqd = [] - for pc in part: - # decide if phase-specific analysis is performed - if iphase is not None and iphase != pc.phasenum: - continue - if pc.inner is not None: - pts = pc.inner.points - rdict = minimize(pts_in_ellips, x0=mc, args=(pts,), method='SLSQP', - constraints=cons, options=opt) - if not rdict['success']: - if verbose: - print(f'Optimization failed for particle {pc.id}') - print(rdict['message']) - continue - eval, evec = np.linalg.eig(arr2mat(rdict['x'])) - if any(eval <= minval): - if verbose: - print(f'Matrix for particle {pc.id} not positive definite or semi-axes too large. ' - f'Eigenvalues: {eval}') - continue - if verbose: - print(f'Optimization succeeded for particle {pc.id} after {rdict["nit"]} iterations.') - print(f'Eigenvalues: {eval}') - print(f'Eigenvectors: {evec}') - # Semi-axes of ellipsoid - ea = 1. / np.sqrt(eval[0]) - eb = 1. / np.sqrt(eval[1]) - ec = 1. / np.sqrt(eval[2]) - eqd = 2.0 * (ea * eb * ec) ** (1.0 / 3.0) - arr_a.append(ea) - arr_b.append(eb) - arr_c.append(ec) - arr_eqd.append(eqd) - else: - eqd = 2.0 * (pc.a * pc.b * pc.c) ** (1.0 / 3.0) - arr_a.append(pc.a) - arr_b.append(pc.b) - arr_c.append(pc.c) - arr_eqd.append(eqd) - - # calculate statistical parameters - part_stats_dict = calc_stats_dict(arr_a, arr_b, arr_c, arr_eqd) - if verbose: - print('\n--------------------------------------------------') - print('Statistical microstructure parameters of particles') - print('--------------------------------------------------') - print('Median lengths of semi-axes of fitted ellipsoids in micron') - print(f'a: {part_stats_dict["a_scale"]:.3f}, b: {part_stats_dict["b_scale"]:.3f}, ' - f'c: {part_stats_dict["c_scale"]:.3f}') - av_std = np.mean([part_stats_dict['a_sig'], part_stats_dict['b_sig'], part_stats_dict['c_sig']]) - print(f'Average standard deviation of semi-axes: {av_std:.4f}') - print('\nAssuming rotational symmetry in grains') - print(f'Rotational axis: {part_stats_dict["ind_rot"]}') - print(f'Median aspect ratio: {part_stats_dict["ar_scale"]:.3f}') - print('\nGrain size') - print(f'Median equivalent grain diameter: {part_stats_dict["eqd_scale"]:.3f} micron') - print(f'Standard deviation of equivalent grain diameter: {part_stats_dict["eqd_sig"]:.4f}') - print('--------------------------------------------------------') - if show_plot: - if part[0].inner is None: - title = 'Particle statistics' - else: - title = 'Statistics of inner particle structures' - if iphase is not None: - title += f' (phase {iphase})' - plot_stats_dict(part_stats_dict, title=title, save_files=save_files) - - return part_stats_dict diff --git a/src/kanapy/textures.py b/src/kanapy/textures.py index 3127efbb..9d7def69 100644 --- a/src/kanapy/textures.py +++ b/src/kanapy/textures.py @@ -114,7 +114,7 @@ def __init__(self, fname, matname=None, gs_min=3, vf_min=0.03, plot=True): data['grains'] = eng.eval("grains_w(grains_w.grainSize > {})" .format(gs_min)) - # use MTEX function to analye grains and plot ellipses around grain + # use MTEX function to analyze grains and plot ellipses around grain # centres; calculate orientation, long and short axes of ellipses omega_r, ha, hb = eng.principalComponents(data['grains'], nargout=3) @@ -165,9 +165,11 @@ def __init__(self, fname, matname=None, gs_min=3, vf_min=0.03, plot=True): # grain equivalent diameter deq = 2.0 * np.array(eng.equivalentRadius(data['grains']))[:, 0] - dsig, doffs, dscale = lognorm.fit(deq) # fit log normal distribution - med_eq = lognorm.median(dsig, loc=doffs, scale=dscale) - std_eq = lognorm.std(dsig, loc=doffs, scale=dscale) + # dsig, doffs, dscale = lognorm.fit(deq) # fit log normal distribution + doffs = 0. + deq_log = np.log(deq) + dscale = med_eq = np.exp(np.median(deq_log)) # lognorm.median(dsig, loc=doffs, scale=dscale) + dsig = std_eq = np.std(deq_log) # lognorm.std(dsig, loc=doffs, scale=dscale) data['gs_param'] = np.array([dsig, doffs, dscale]) data['gs_data'] = deq data['gs_moments'] = [med_eq, std_eq] @@ -186,9 +188,11 @@ def __init__(self, fname, matname=None, gs_min=3, vf_min=0.03, plot=True): # grain aspect ratio asp = np.array(eng.aspectRatio(data['grains']))[:, 0] - asig, aoffs, ascale = lognorm.fit(asp) # fit log normal distribution - med_ar = lognorm.median(asig, loc=aoffs, scale=ascale) - std_ar = lognorm.std(asig, loc=aoffs, scale=ascale) + # asig, aoffs, ascale = lognorm.fit(asp) # fit log normal distribution + aoffs = 0. + asp_log = np.log(asp) + ascale = med_ar = np.exp(np.median(asp_log)) # lognorm.median(asig, loc=aoffs, scale=ascale) + asig = std_ar = np.std(asp_log) # lognorm.std(asig, loc=aoffs, scale=ascale) data['ar_param'] = np.array([asig, aoffs, ascale]) data['ar_data'] = asp data['ar_moments'] = [med_ar, std_ar] @@ -207,17 +211,18 @@ def __init__(self, fname, matname=None, gs_min=3, vf_min=0.03, plot=True): # angles of main axis # fit von Mises distribution (circular normal distribution) to data - kappa, oloc, oscale = vonmises.fit(omega) - med_om = vonmises.median(kappa, loc=oloc, scale=oscale) - std_om = vonmises.std(kappa, loc=oloc, scale=oscale) - data['om_param'] = np.array([kappa, oloc, oscale]) - data['om_data'] = omega + omega_p = 2.0*omega - np.pi # rescale angles from [0, pi] to [-pi,pi] for von Mises distr. fit + kappa, oloc, oscale = vonmises.fit(omega_p) + med_om = vonmises.median(kappa, loc=oloc) # scale parameter has no effect on vonmises distribution + std_om = vonmises.std(kappa, loc=oloc) + data['om_param'] = np.array([kappa, oloc]) + data['om_data'] = omega_p data['om_moments'] = [med_om, std_om] if plot: fig, ax = plt.subplots() - x = np.linspace(-np.pi, np.pi, 200) # np.amin(omega), np.amax(omega), 150) - y = vonmises.pdf(x, kappa, loc=oloc, scale=oscale) - ax.plot(x, y, '-r', label='fit') + x = np.linspace(-np.pi, np.pi, 200) # np.amin(omega), np.amax(omega), 150) + y = vonmises.pdf(x, kappa, loc=oloc) + ax.plot(0.5*(x+np.pi), 2*y, '-r', label='fit') ax.hist(omega, bins=40, density=True, label='data') plt.legend() plt.title('Histogram of tilt angles of major axes') diff --git a/src/kanapy/voxelization.py b/src/kanapy/voxelization.py index 4845aac8..b796c6a9 100644 --- a/src/kanapy/voxelization.py +++ b/src/kanapy/voxelization.py @@ -356,8 +356,7 @@ def poly2vox(): vox.append(iv) pa.inside_voxels.append(iv) if gid in mesh.grain_dict.keys(): - for iv in vox: - mesh.grain_dict[gid].append(iv) + mesh.grain_dict[gid].extend(iv) else: mesh.grain_dict[gid] = vox diff --git a/tests/test_packing.py b/tests/test_packing.py index 0c53642f..43fadaa2 100644 --- a/tests/test_packing.py +++ b/tests/test_packing.py @@ -87,8 +87,7 @@ def rot_surf(): def test_particle_generator(par_sim): - periodic = False - ell_list = particle_generator([par_sim[0]], par_sim[1], periodic, None) + ell_list = particle_generator([par_sim[0]], par_sim[1], None) for ell in ell_list: assert isinstance(ell, Ellipsoid) diff --git a/tests/test_texture.py b/tests/test_texture.py index 2a34adad..89342c4e 100644 --- a/tests/test_texture.py +++ b/tests/test_texture.py @@ -25,15 +25,15 @@ def test_mex(): @pytest.mark.skipif(not MTEX_AVAIL, reason="Kanapy is not configured for texture analysis yet!") def test_readEBSD(): from kanapy.textures import EBSDmap - fname = MAIN_DIR + '/tests/ebsd_316L_1000x1000.ang' # name of ang file to be imported + fname = MAIN_DIR + '/tests/ebsd_316L_500x500.ang' # name of ang file to be imported # read EBSD map and evaluate statistics of microstructural features ebsd = EBSDmap(os.path.normpath(fname), plot=False) - assert(len(ebsd.ms_data) == 1) + assert len(ebsd.ms_data) == 1 gs_param = ebsd.ms_data[0]['gs_param'] - assert(np.abs(gs_param[0] - 0.99765477) < 1.e-5) + assert np.abs(gs_param[0] - 0.7177939893510182) < 1.e-5 # get list of orientations for grains in RVE matching the ODF of the EBSD map ori_rve = ebsd.calcORI(20) - assert(np.abs(ori_rve[0, 1] - 0.26179939) < 1.e-5) + assert np.abs(ori_rve[0, 1] - 0.5817764173314431) < 1.e-5 @pytest.mark.skipif(not MTEX_AVAIL, reason="Kanapy is not configured for texture analysis yet!") def test_createORI(): diff --git a/tests/test_voxelization.py b/tests/test_voxelization.py index 1edd8dc9..35ea2648 100644 --- a/tests/test_voxelization.py +++ b/tests/test_voxelization.py @@ -211,7 +211,7 @@ def test_voxelizationRoutine(): mesh = mesh_creator(rve.dim) nphases = 1 mesh.create_voxels(sim_box) - Particles = particle_generator([parDict], sim_box, False, None) + Particles = particle_generator([parDict], sim_box, None) mesh = voxelizationRoutine(Particles, mesh, nphases) assert mesh.vox_center_dict[5][0] == 0.9 assert mesh.voxel_dict[7] == [29, 30, 26, 25, 31, 32, 28, 27]