From 6862016b5984afabdad5cc061c8bed91f59329ae Mon Sep 17 00:00:00 2001 From: AHartmaier Date: Sat, 3 Feb 2024 16:52:36 +0100 Subject: [PATCH] Improved support for dual-phase microstructures --- setup.py | 2 +- src/kanapy/api.py | 13 ++++++++++--- src/kanapy/gb_textureReconstruction.m | 2 +- src/kanapy/import_EBSD.py | 3 +-- src/kanapy/initializations.py | 9 +++++++++ src/kanapy/input_output.py | 14 ++++++++------ src/kanapy/plotting.py | 24 +++++++++++++----------- src/kanapy/textureReconstruction.m | 16 +++++----------- 8 files changed, 48 insertions(+), 35 deletions(-) diff --git a/setup.py b/setup.py index 2b977cef..7b56fb37 100644 --- a/setup.py +++ b/setup.py @@ -48,7 +48,7 @@ setup( name='kanapy', - version='6.1.2', + version='6.1.3', author='Mahesh R.G. Prasad, Abhishek Biswas, Golsa Tolooei Eshlaghi, Napat Vajragupta, Alexander Hartmaier', author_email='alexander.hartmaier@rub.de', classifiers=[ diff --git a/src/kanapy/api.py b/src/kanapy/api.py index 4a7e3cb7..a5fba3ee 100644 --- a/src/kanapy/api.py +++ b/src/kanapy/api.py @@ -19,7 +19,7 @@ from kanapy.grains import calc_polygons, get_stats from kanapy.entities import Simulation_Box -from kanapy.input_output import export2abaqus, writeAbaqusMat +from kanapy.input_output import export2abaqus, writeAbaqusMat, read_dump from kanapy.initializations import RVE_creator, mesh_creator from kanapy.packing import packingRoutine from kanapy.voxelization import voxelizationRoutine @@ -215,7 +215,7 @@ def generate_grains(self): statistical comparison. Final RVE grain volumes and shared grain boundary surface areas info are written out as well.""" - if self.mesh.grains is None: + if self.mesh is None or self.mesh.grains is None: raise ValueError('No information about voxelized microstructure. Run voxelize first.') if self.porosity and 0 in self.mesh.grain_dict.keys(): # in case of porosity, remove irregular grain 0 from analysis @@ -317,7 +317,8 @@ def generate_orientations(self, data, ang=None, omega=None, Nbase=5000, '"random" or "unimodal"') for i, igr in enumerate(self.mesh.grain_dict.keys()): if self.mesh.grain_phase_dict[igr] == ip: - ori_dict[igr] = ori_rve[i, :] + ind = i - ip*self.ngrains[0] + ori_dict[igr] = ori_rve[ind, :] self.mesh.grain_ori_dict = ori_dict return @@ -1043,6 +1044,12 @@ def pckl(self, file=None, path='./'): pickle.dump(self, output, pickle.HIGHEST_PROTOCOL) return + def import_particles(self, file, path='./'): + path = os.path.normpath(path) + file = os.path.join(path, file) + self.simbox, self.particles = read_dump(file) + + """ -------- legacy methods -------- """ diff --git a/src/kanapy/gb_textureReconstruction.m b/src/kanapy/gb_textureReconstruction.m index 4604d33c..c504763f 100644 --- a/src/kanapy/gb_textureReconstruction.m +++ b/src/kanapy/gb_textureReconstruction.m @@ -48,7 +48,7 @@ % switch index index([i1 i2]) = index([i2 i1]); orilist([i1 i2]) = orilist([i2 i1]); -an = angle(orilist'*orilist); +an = angle(orilist'*orilist); %' a_out = an(bin_score_area~=0); [~,loc] = histc(a_out,bins); et = sparse(loc,1:length(loc),area,nbin,length(loc)); diff --git a/src/kanapy/import_EBSD.py b/src/kanapy/import_EBSD.py index 4a9242db..785a97fb 100644 --- a/src/kanapy/import_EBSD.py +++ b/src/kanapy/import_EBSD.py @@ -265,13 +265,12 @@ def calcORI(self, Ng, iphase=0, shared_area=None, nbins=12): Array with Ng Euler angles. """ - ms = self.ms_data[iphase] orired, odfred, ero = \ self.eng.textureReconstruction(Ng, 'orientation', ms['ori'], 'grains', ms['grains'], nargout=3) - if shared_area is None: + if shared_area is None or shared_area == 0: return np.array(self.eng.Euler(orired)) else: orilist, ein, eout, mbin = \ diff --git a/src/kanapy/initializations.py b/src/kanapy/initializations.py index 12c9da5c..dcf1f3a0 100644 --- a/src/kanapy/initializations.py +++ b/src/kanapy/initializations.py @@ -225,6 +225,9 @@ def gen_data_elong(pdict): 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 pdict = gen_data_basic(dict({'Type': stats["Grain type"], 'Phase': ip})) @@ -236,12 +239,18 @@ def gen_data_elong(pdict): loc_AR = stats["Aspect ratio"][loc] 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: + raise ValueError('Min/Max values for cutoffs of aspect ratio are too close: ' + + f'Max: {ar_cutoff_max}, Min: {ar_cutoff_min}') # Extract grain tilt angle statistics info from dict kappa = stats["Tilt angle"][kappa] loc_ori = stats["Tilt angle"][loc] ori_cutoff_min = stats["Tilt angle"]["cutoff_min"] ori_cutoff_max = stats["Tilt angle"]["cutoff_max"] + 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) diff --git a/src/kanapy/input_output.py b/src/kanapy/input_output.py index 9f7aa63d..5ae27131 100644 --- a/src/kanapy/input_output.py +++ b/src/kanapy/input_output.py @@ -8,7 +8,8 @@ def write_dump(Ellipsoids, sim_box): """ - Writes the (.dump) file, which can be read by visualization software OVITO. + Writes the (.dump) files into a sub-directory "dump_files", which can be read by visualization software OVITO + or imported again into Kanapy to avoid the packing simulation. :param Ellipsoids: Contains information of ellipsoids such as its position, axes lengths and tilt angles :type Ellipsoids: list @@ -19,8 +20,8 @@ def write_dump(Ellipsoids, sim_box): """ num_particles = len(Ellipsoids) cwd = os.getcwd() - output_dir = cwd + '/dump_files' # output directory - dump_outfile = output_dir + '/particle' # output dump file + output_dir = os.path.join(cwd, 'dump_files') + dump_outfile = os.path.join(output_dir, 'particle') # output dump file if not os.path.exists(output_dir): os.makedirs(output_dir) @@ -38,8 +39,8 @@ def write_dump(Ellipsoids, sim_box): 'OrientationZ OrientationW\n') for ell in Ellipsoids: qw, qx, qy, qz = ell.quat - f.write('{0} {1} {2} {3} {4} {5} {6} {7} {8} {9} {10}\n'.format( - ell.id, ell.x, ell.y, ell.z, ell.a, ell.b, ell.c, qx, qy, qz, qw)) + f.write('{0} {1} {2} {3} {4} {5} {6} {7} {8} {9} {10} {11}\n'.format( + ell.id, ell.x, ell.y, ell.z, ell.a, ell.b, ell.c, qx, qy, qz, qw, ell.phasenum)) def read_dump(file): @@ -98,8 +99,9 @@ def read_dump(file): a, b, c = values[4], values[5], values[6] # Semi-major length, Semi-minor length 1 & 2 x, y, z = values[1], values[2], values[3] qx, qy, qz, qw = values[7], values[8], values[9], values[10] + ip = int(values[11]) quat = np.array([qw, qx, qy, qz]) - ellipsoid = Ellipsoid(iden, x, y, z, a, b, c, quat) # instance of Ellipsoid class + ellipsoid = Ellipsoid(iden, x, y, z, a, b, c, quat, phasenum=ip) # instance of Ellipsoid class # Find the original particle if the current is duplicate for c in values[0]: diff --git a/src/kanapy/plotting.py b/src/kanapy/plotting.py index adcc8dd2..b1bad3b2 100644 --- a/src/kanapy/plotting.py +++ b/src/kanapy/plotting.py @@ -244,7 +244,7 @@ def plot_output_stats(dataDict, fig, ax = plt.subplots(1, 2, figsize=(15, 9)) # Plot histogram - ax[0].hist(data, density=False, bins=binNum, label=label) + ax[0].hist(data, density=True, bins=binNum, label=label) ax[0].legend(loc="upper right", fontsize=16) ax[0].set_xlabel('Equivalent diameter (μm)', fontsize=18) ax[0].set_ylabel('Frequency', fontsize=18) @@ -307,11 +307,11 @@ def plot_output_stats(dataDict, par_AR = np.sort(np.asarray(dataDict['Particle_Major_diameter']) / np.asarray(dataDict['Particle_Minor_diameter'])) # Concatenate corresponding arrays to compute shared bins - total_AR = np.concatenate([par_AR, grain_AR]) + total_AR = np.concatenate([grain_AR, par_AR]) sig_par, loc_par, sc_par = lognorm.fit(par_AR) par_lognorm = lognorm(sig_par, loc=loc_par, scale=sc_par) - data = [par_AR, grain_AR] - label = ['Particles', 'Grains'] + data = [grain_AR, par_AR] + label = [ 'Grains', 'Particles'] else: total_AR = grain_AR data = [grain_AR] @@ -326,26 +326,28 @@ def plot_output_stats(dataDict, sns.set(color_codes=True) fig, ax = plt.subplots(1, 2, figsize=(15, 9)) # Plot histogram - ax[0].hist(data, density=False, bins=len(shared_AR), label=label) + ax[0].hist(data, density=True, bins=len(shared_AR), label=label) ax[0].legend(loc="upper right", fontsize=16) ax[0].set_xlabel('Aspect ratio', fontsize=18) ax[0].set_ylabel('Frequency', fontsize=18) ax[0].tick_params(labelsize=14) # Plot PDF - if particles and 'Particle_Minor_diameter' in dataDict.keys(): - ypdf1 = par_lognorm.pdf(par_AR) - area = np.trapz(ypdf1, par_AR) - ypdf1 /= area - ax[1].plot(par_AR, ypdf1, linestyle='-', linewidth=3.0, label='Particles') - ax[1].fill_between(par_AR, 0, ypdf1, alpha=0.3) ypdf2 = grain_lognorm.pdf(grain_AR) area = np.trapz(ypdf2, grain_AR) if np.isclose(area, 0.0): + logging.debug('Small area for aspect ratio of grains.') + logging.debug(ypdf2, grain_AR) area = 1.0 ypdf2 /= area ax[1].plot(grain_AR, ypdf2, linestyle='-', linewidth=3.0, label='Grains') ax[1].fill_between(grain_AR, 0, ypdf2, alpha=0.3) + if particles and 'Particle_Minor_diameter' in dataDict.keys(): + ypdf1 = par_lognorm.pdf(par_AR) + area = np.trapz(ypdf1, par_AR) + ypdf1 /= area + ax[1].plot(par_AR, ypdf1, linestyle='-', linewidth=3.0, label='Particles') + ax[1].fill_between(par_AR, 0, ypdf1, alpha=0.3) if ar_param is not None: x0 = np.amin(1.0) x1 = np.amax(grain_AR) diff --git a/src/kanapy/textureReconstruction.m b/src/kanapy/textureReconstruction.m index 383d0046..f8d86b94 100644 --- a/src/kanapy/textureReconstruction.m +++ b/src/kanapy/textureReconstruction.m @@ -32,7 +32,7 @@ % Output: reduced orientation set, ODF and L1 error % %% input fields and checks - +mtex_progress = 1; % Avoid progress output options = {'ebsdMatFile','ebsd','orientation'}; flag = 1; grains = []; @@ -136,17 +136,11 @@ ll=10; hl = 55; - step = 1; - ero=100; - e_mod = []; - lim = 10; - hh=[]; - kappa = psi.halfwidth*180/pi; % initial kernel shape %% tic @@ -218,7 +212,7 @@ oriseq_p = oriseq(fval>1); -ind = num2cell(fval(fval>1)'); +ind = num2cell(fval(fval>1)'); %' pendList = cellfun(@splitMean, oriseq_p, ind, 'UniformOutput', false); ori_f = [ori_f [pendList{:}]]; @@ -227,12 +221,12 @@ hh = [hh h]; er = calcError(odf,odfred); -er = calcError(odf,odfred); +%er = calcError(odf,odfred); if er