Improved support for dual-phase microstructures

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='[email protected]',
     classifiers=[        

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        --------
     """

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));

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 = \

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)

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]:

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)

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<ero
-    orired_f=ori_f;
+    orired_f = ori_f;
     odfred_f = odfred;
     ohw = h;
-    ero=er;
+    ero = er;
 end
 e_mod = [e_mod er];