Improved support for export of dual-phase microstructures

setup.py

@@ -48,7 +48,7 @@
 
 setup(
     name='kanapy',
-    version='6.1.4',
+    version='6.1.5',
     author='Mahesh R.G. Prasad, Abhishek Biswas, Golsa Tolooei Eshlaghi, Napat Vajragupta, Alexander Hartmaier',
     author_email='[email protected]',
     classifiers=[        

src/kanapy/api.py

@@ -64,7 +64,7 @@ class Microstructure(object):
                 nvox, phases, prec_vf_voxels, vox_center_dict, voxel_dict
         geometry : dict
             Dictionary of grain geometries;
-            Dict keys: "Ngrains", "Vertices", "Points", "Simplices", "Facets", "Grains", "GBnodes", GBarea" "GBfaces"
+            Dict keys: "Vertices", "Points", "Simplices", "Facets", "Grains", "GBnodes", GBarea" "GBfaces"
             "Grains" : dictionary with key grain_number  
             Keys:"Vertices", "Points", "Center", "Simplices", "Volume", "Area", "Phase", "eqDia", "majDia", "minDia"
         res_data : list
@@ -252,11 +252,11 @@ def generate_grains(self):
         self.geometry: dict = \
             calc_polygons(self.rve, self.mesh)  # updates RVE_data
         # verify that geometry['Grains'] and mesh.grain_dict are consistent
-        if np.any(self.geometry['Ngrains'] != self.ngrains):
+        """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.')
+                    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 = []
@@ -290,7 +290,8 @@ def generate_grains(self):
                 print(f'{ip}: {self.rve.phase_names[ip]} ({vf.round(1)}%)')
         # analyse grains w.r.t. statistical descriptors
         self.res_data = \
-            get_stats(self.rve.particle_data, self.geometry, self.rve.units, nphases)
+            get_stats(self.rve.particle_data, self.geometry, self.rve.units,
+                      nphases, self.mesh.ngrains_phase)
         if empty_vox is not None:
             # add removed grain again
             self.mesh.grain_dict[0] = empty_vox
@@ -481,8 +482,36 @@ def plot_slice(self, cut='xy', data=None, pos=None, fname=None,
     """
 
     def write_abq(self, nodes=None, file=None, path='./', voxel_dict=None, grain_dict=None,
-                  dual_phase=False, thermal=False, units=None, ialloy=None):
-        """ Writes out the Abaqus (.inp) file for the generated RVE."""
+                  dual_phase=False, thermal=False, units=None,
+                  ialloy=None, nsdv=200):
+        """
+        Writes out the Abaqus deck (.inp file) for the generated RVE. The parameter nodes should be
+        a string indicating if voxel ("v") or smoothened ("s") mesh should be written. It can also
+        provide an array of nodal positions. If dual_phase is true, the
+        generated deck contains plain material definitions for each phase. Material parameters must
+        be specified by the user. If ialloy is provided, the generated deck material definitions
+        for each grain. For dual phase structures to be used with crystal plasticity, ialloy
+        can be a list with all required material definitions. If the list ialloy is shorted than the
+        number of phases in the RVE, plain material definitions for the remaining phases will
+        be included in the input deck.
+
+        Parameters
+        ----------
+        nodes
+        file
+        path
+        voxel_dict
+        grain_dict
+        dual_phase
+        thermal
+        units
+        ialloy
+        nsdv
+
+        Returns
+        -------
+
+        """
         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')
@@ -537,6 +566,13 @@ def write_abq(self, nodes=None, file=None, path='./', voxel_dict=None, grain_dic
             nct = f'abq_px_{len(grain_dict)}'
         if ialloy is None:
             ialloy = self.rve.ialloy
+        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 = self.mesh.grain_phase_dict
+        else:
+            grpd = None
         if file is None:
             if self.name == 'Microstructure':
                 file = nct + ntag + '_geom.inp'
@@ -546,11 +582,14 @@ def write_abq(self, nodes=None, file=None, path='./', voxel_dict=None, grain_dic
         file = os.path.join(path, file)
         export2abaqus(nodes, file, grain_dict, voxel_dict,
                       units=units, gb_area=faces,
-                      dual_phase=dual_phase, thermal=thermal)
+                      dual_phase=dual_phase,
+                      ialloy=ialloy, grain_phase_dict=grpd,
+                      thermal=thermal)
         # if orientation exists also write material file with Euler angles
         if 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')
+                           file=file[0:-8] + 'mat.inp',
+                           grain_phase_dict=grpd, nsdv=nsdv)
         return file
 
     def write_abq_ori(self, ialloy=None, ori=None, file=None, path='./', nsdv=200):
@@ -1025,7 +1064,7 @@ def write_voxels(self, angles=None, script_name=None, file=None, path='./',
             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_dict.keys():
+                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()):

src/kanapy/grains.py

@@ -2,7 +2,6 @@
 import logging
 import numpy as np
 from scipy.spatial import ConvexHull, Delaunay
-from scipy.spatial.distance import euclidean
 from tqdm import tqdm
 
 
@@ -230,7 +229,6 @@ def project_pts(pts, ctr, axis):
 
     # create dicts for GB facets, including fake facets at surfaces
     geometry = dict()
-    geometry['Ngrains'] = mesh.ngrains_phase
     grain_facesDict = dict()  # {Grain: faces}
     gb_vox_dict = dict()
     for i in range(1, Ng_max + 7):
@@ -593,7 +591,7 @@ def project_pts(pts, ctr, axis):
     return geometry
 
 
-def get_stats(particle_data, geometry, units, nphases):
+def get_stats(particle_data, geometry, units, nphases, ngrains):
     """
     Compare the geometries of particles used for packing and the resulting
     grains.
@@ -637,7 +635,7 @@ def get_stats(particle_data, geometry, units, nphases):
     output_data_list = []
     for ip in range(nphases):
         # Create dictionaries to store the data generated
-        output_data = {'Number': geometry['Ngrains'][ip],
+        output_data = {'Number': ngrains[ip],
                        'Unit_scale': units,
                        'Grain_Equivalent_diameter': np.array(grain_eqDia[ip]),
                        'Grain_Major_diameter': np.array(grain_majDia[ip]),

src/kanapy/initializations.py

@@ -224,7 +224,7 @@ 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:
+            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
@@ -238,7 +238,7 @@ 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:
+                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}')
 
@@ -247,7 +247,7 @@ def gen_data_elong(pdict):
                 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:
+                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}')
 
@@ -267,11 +267,12 @@ def gen_data_elong(pdict):
         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.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 = []
 
         # extract data from descriptors of individual phases
         for ip, stats in enumerate(stats_dicts):
@@ -297,13 +298,17 @@ def gen_data_elong(pdict):
                                         'Using first value.')
                 # Extract Alloy number for ICAMS CP-UMAT
                 if "ialloy" in stats['RVE'].keys():
-                    self.ialloy = stats['RVE']['ialloy']
+                    ialloy.append(stats['RVE']['ialloy'])
             elif ip == 0:
                 raise ValueError('RVE properties must be specified in descriptors for first phase.')
 
             # Extract other simulation attributes, must be specified for phase 0
             if "Simulation" in stats.keys():
-                periodic = bool(stats["Simulation"]["periodicity"])
+                peri = stats["Simulation"]["periodicity"]
+                if type(peri) is bool:
+                    periodic = peri
+                else:
+                    periodic = True if peri == 'True' else False
                 if self.periodic is None:
                     self.periodic = periodic
                 elif self.periodic != periodic:
@@ -348,6 +353,11 @@ def gen_data_elong(pdict):
         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):
+                logging.warning(f'{nall} values of "ialloy" provided, but only {len(phase_vf)} phases defined.')
+            self.ialloy = ialloy
         return
 
 
@@ -515,7 +525,10 @@ def set_stats(grains, ar=None, omega=None, deq_min=None, deq_max=None,
     # number of voxels
     nx = ny = nz = voxels  # number of voxels in each direction
     # specify RVE info
-    pbc = 'True' if periodicity else 'False'
+    if type(periodicity) is bool:
+        pbc = periodicity
+    else:
+        pbc = True if periodicity == 'True' else False
 
     # check grain type
     # create the corresponding dict with statistical grain geometry information

src/kanapy/input_output.py

@@ -119,10 +119,17 @@ def read_dump(file):
 
 
 def export2abaqus(nodes, file, grain_dict, voxel_dict, units='mm',
-                  gb_area=None, dual_phase=False, thermal=False):
+                  gb_area=None, dual_phase=False, thermal=False,
+                  ialloy=None, grain_phase_dict=None):
     r"""
     Creates an ABAQUS input file with microstructure morphology information
-    in the form of nodes, elements and element sets.
+    in the form of nodes, elements and element sets. If "dual_phase" is true,
+    element sets with phase numbers will be defined and assigned to materials
+    "PHASE_{phase_id}MAT" plain material definitions for phases will be included.
+    Otherwise, it will be assumed that each grain refers to a material
+    "GRAIN_{}grain_id}MAT. In this case, a "_mat.inp" file with the same name
+    trunc will be included, in which the alloy number and Euler angles for each
+    grain must be defined.
 
     .. note:: The nodal coordinates are written out in units of 1 mm or 1 :math:`\mu` m scale, as requested by the
                  user in the input file.
@@ -142,9 +149,15 @@ def write_grain_sets():
                     f.write('%d\n' % el)
         # Create sections
         for k in grain_dict.keys():
-            f.write(
-                '*Solid Section, elset=GRAIN{0}_SET, material=GRAIN{1}_MAT\n'
-                .format(k, k))
+            if grain_phase_dict is None or grain_phase_dict[k] < nall:
+                f.write(
+                    '*Solid Section, elset=GRAIN{0}_SET, material=GRAIN{1}_MAT\n'
+                    .format(k, k))
+            else:
+                f.write(
+                    '*Solid Section, elset=GRAIN{0}_SET, material=PHASE{1}_MAT\n'
+                    .format(k, grain_phase_dict[k]))
+                ph_set.add(grain_phase_dict[k])
         return
 
     def write_phase_sets():
@@ -164,6 +177,7 @@ def write_phase_sets():
             f.write(
                 '*Solid Section, elset=PHASE{0}_SET, material=PHASE{1}_MAT\n'
                 .format(k, k))
+            ph_set.add(k)
         return
 
     print('')
@@ -180,6 +194,11 @@ def write_phase_sets():
         print('Using tet element type SFM3D4.')
         eltype = 'SFM3D4'
 
+    # select material definition
+    if type(ialloy) is not list:
+        ialloy = [ialloy]
+    nall = len(ialloy)
+    ph_set = set()
     # Factor used to generate nodal coordinates in 'mm' or 'um' scale
     if units == 'mm':
         scale_fact = 1.e-3
@@ -262,25 +281,38 @@ def write_phase_sets():
         f.write('** MATERIALS\n')
         f.write('**\n')
         if dual_phase:
-            for i in range(1, len(grain_dict.keys())):
+            # declare plane materials for Abaqus standard materials for each phase
+            for i in ph_set:
                 f.write('**\n')
                 f.write('*Material, name=PHASE{}_MAT\n'.format(i))
                 f.write('**Include, input=Material{}.inp\n'.format(i))
                 f.write('**\n')
         else:
+            # declare plane materials for Abaqus standard materials
+            for i in ph_set:
+                f.write('**\n')
+                f.write('*Material, name=PHASE{}_MAT\n'.format(i))
+            # include file with definition for CP parameters for each grain
+            f.write('**\n')
             f.write('*Include, input={}mat.inp\n'.format(file[0:-8]))
             f.write('**\n')
     print('---->DONE!\n')
     return
 
 
-def writeAbaqusMat(ialloy, angles, file=None, path='./', nsdv=200):
+def writeAbaqusMat(ialloy, angles,
+                   file=None, path='./',
+                   grain_phase_dict=None,
+                   nsdv=200):
     """
-    Export Euler angles to Abaqus input deck that can be included in the _geom.inp file.
+    Export Euler angles to Abaqus input deck that can be included in the _geom.inp file. If
+    parameter "grain_phase_dict" is given, the phase number for each grain will be used to select
+    the corresponding ialloy from a list. If the list ialloy is shorter than the number of phases in
+    grain_phase_dict, no angles for phases with no corresponding ialloy will be written.
 
     Parameters:
     -----------
-    ialloy : int
+    ialloy : int or list
         Identifier, alloy number in ICAMS CP-UMAT: mod_alloys.f
     angles : dict or (N, 3)-ndarray
         Dict with Euler angles for each grain or array with number of N rows (= number of grains) and
@@ -289,17 +321,21 @@ def writeAbaqusMat(ialloy, angles, file=None, path='./', nsdv=200):
         Filename, optional (default: None)
     path : str
         Path to save file, option (default: './')
+    grain_phase_dict: dict
+        Dict with phase for each grain, optional (default: None)
     nsdv : int
         Number of state dependant variables, optional (default: 200)
     """
+    if type(ialloy) is not list:
+        ialloy = [ialloy]
+    nall = len(ialloy)
     if type(angles) is not dict:
         # converting (N, 3) ndarray to dict
         gr_ori_dict = dict()
         for igr, ori in enumerate(angles):
             gr_ori_dict[igr+1] = ori
     else:
         gr_ori_dict = angles
-
     nitem = len(gr_ori_dict.keys())
     if file is None:
         file = f'abq_px_{nitem}_mat.inp'
@@ -310,11 +346,17 @@ def writeAbaqusMat(ialloy, angles, file=None, path='./', nsdv=200):
         f.write('** MATERIALS\n')
         f.write('**\n')
         for igr, ori in gr_ori_dict.items():
+            if grain_phase_dict is None:
+                ip = 0
+            else:
+                ip = grain_phase_dict[igr]
+                if ip > nall - 1:
+                    continue
             f.write('*Material, name=GRAIN{}_MAT\n'.format(igr))
             f.write('*Depvar\n')
             f.write('    {}\n'.format(nsdv))
             f.write('*User Material, constants=4\n')
-            f.write('{}, {}, {}, {}\n'.format(float(ialloy),
+            f.write('{}, {}, {}, {}\n'.format(float(ialloy[ip]),
                                               ori[0], ori[1], ori[2]))
     return