Merge branch 'development' into fix_homogenization_restart

VERSION

@@ -1 +1 @@
-v3.0.0-alpha2-530-g0d0226f70
+v3.0.0-alpha2-602-ge2d4ab427

examples/ConfigFiles/Microstructure_DP_Steel.config

@@ -1,4 +0,0 @@
-[DP_Steel]
-crystallite 1
-(constituent)	phase	1	texture	1 fraction	0.82
-(constituent)	phase	2	texture	2	fraction	0.18 

examples/ConfigFiles/Phase_Dislotwin_TWIP-Steel-FeMnC.config

@@ -1,64 +0,0 @@
-[TWIP_Steel_FeMnC]
-
-elasticity          hooke
-plasticity          dislotwin
-
-(output)           rho_mob
-(output)           rho_dip
-(output)           gamma_sl
-(output)           lambda_sl
-(output)           tau_pass
-(output)           f_tw
-(output)           lambda_tw
-(output)           tau_hat_tw
-(output)           f_tr
-
-
-### Material parameters ###
-lattice_structure	fcc
-C11                   175.0e9       # From Music et al. Applied Physics Letters 91, 191904 (2007)
-C12                   115.0e9
-C44                   135.0e9
-grainsize             2.0e-5        # Average grain size [m]
-SolidSolutionStrength 1.5e8         # Strength due to elements in solid solution
-
-### Dislocation glide parameters ###
-Nslip               12
-slipburgers         2.56e-10        # Burgers vector of slip system [m]
-rhoedgedip0         1.0             # Initial dislocation density [m/m**3]
-rhoedge0            1.0e12          # Initial dislocation density [m/m**3]
-v0                  1.0e-4          # Initial glide velocity [m/s]
-Qedge               3.7e-19         # Activation energy for dislocation glide [J]
-p_slip              1.0             # p-exponent in glide velocity
-q_slip              1.0             # q-exponent in glide velocity
-
-# hardening of glide
-CLambdaSlip         10.0            # Adj. parameter controlling dislocation mean free path
-D0                  4.0e-5          # Vacancy diffusion prefactor [m**2/s]
-Qsd                 4.5e-19         # Activation energy for climb [J]
-Catomicvolume       1.0             # Adj. parameter controlling the atomic volume [in b^3]
-Cedgedipmindistance 1.0             # Adj. parameter controlling the minimum dipole distance [in b]
-interactionSlipSlip 0.122 0.122 0.625 0.07 0.137 0.122 # Interaction coefficients (Kubin et al. 2008)
-
-### Shearband parameters ###
-shearbandresistance 180e6
-shearbandvelocity   0e-4            # set to zero to turn shear banding of
-QedgePerSbSystem    3.7e-19         # Activation energy for shear banding [J]
-p_shearband         1.0             # p-exponent in glide velocity
-q_shearband         1.0             # q-exponent in glide velocity
-
-### Twinning parameters ###
-Ntwin               12
-twinburgers         1.47e-10        # Burgers vector of twin system [m]
-twinsize            5.0e-8          # Twin stack mean thickness [m]	
-L0_twin             442.0           # Length of twin nuclei in Burgers vectors
-maxtwinfraction     1.0             # Maximum admissible twin volume fraction
-xc_twin             1.0e-9          # critical distance for formation of twin nucleus
-VcrossSlip          1.67e-29        # cross slip volume
-r_twin              10.0            # r-exponent in twin formation probability
-Cmfptwin            1.0             # Adj. parameter controlling twin mean free path
-Cthresholdtwin      1.0             # Adj. parameter controlling twin threshold stress
-interactionSlipTwin 0.0 1.0 1.0     # Dislocation-Twin interaction coefficients
-interactionTwinTwin 0.0 1.0         # Twin-Twin interaction coefficients
-SFE_0K              -0.0396         # stacking fault energy at zero K; TWIP steel: -0.0526; Cu: -0.0396
-dSFE_dT             0.0002          # temperature dependance of stacking fault energy

examples/ConfigFiles/Phase_Dislotwin_TWIP-Steel-FeMnC.yaml

@@ -0,0 +1,41 @@
+TWIP_Steel_FeMnC:
+  lattice: cF
+  mechanics:
+    elasticity: {type: hooke, C_11: 175.0e9, C_12: 115.0e9, C_44: 135.0e9}
+    plasticity:
+      type: dislotwin
+      output: [rho_mob, rho_dip, gamma_sl, Lambda_sl, tau_pass, f_tw, Lambda_tw, tau_hat_tw, f_tr]
+      D: 2.0e-5
+      N_sl: [12]
+      b_sl: [2.56e-10]
+      rho_mob_0: [1.0e12]
+      rho_dip_0: [1.0]
+      v_0: [1.0e4]
+      Q_s: [3.7e-19]
+      p_sl: [1.0]
+      q_sl: [1.0]
+      tau_0: [1.5e8]
+      i_sl: [10.0]                 # Adj. parameter controlling dislocation mean free path
+      D_0: 4.0e-5                  # Vacancy diffusion prefactor  / m^2/s
+      D_a: 1.0                     # minimum dipole distance / b
+      Q_cl: 4.5e-19                # Activation energy for climb / J
+      h_sl_sl: [0.122, 0.122, 0.625, 0.07, 0.137, 0.122]  # Interaction coefficients (Kubin et al. 2008)
+# shear band parameters
+      xi_sb: 180.0e6
+      Q_sb: 3.7e-19
+      p_sb: 1.0
+      q_sb: 1.0
+      v_sb: 0.0                    # set to 0, to turn it off
+# twinning parameters
+      N_tw: [12]
+      b_tw: [1.47e-10]             # Burgers vector length of twin system / b
+      t_tw: [5.0e-8]               # Twin stack mean thickness / m
+      L_tw: 442.0                  # Length of twin nuclei / b
+      x_c_tw: 1.0e-9               # critical distance for formation of twin nucleus / m
+      V_cs:  1.67e-29              # cross slip volume / m^3
+      p_tw: [10.0]                 # r-exponent in twin formation probability
+      i_tw: 1.0                    # Adj. parameter controlling twin mean free path
+      h_sl_tw: [0.0, 1.0, 1.0]     # dislocation-twin interaction coefficients
+      h_tw_tw: [0.0, 1.0]          # twin-twin interaction coefficients
+      Gamma_sf_0K: -0.0396         # stacking fault energy / J/m^2 at zero K; TWIP steel: -0.0526; Cu: -0.0396
+      dGamma_sf_dT: 0.0002         # temperature dependence / J/(m^2 K) of stacking fault energy

examples/ConfigFiles/Phase_Dislotwin_Tungsten.config

@@ -1,36 +0,0 @@
-[Tungsten]
-
-elasticity          hooke
-plasticity          dislotwin
-
-### Material parameters ###
-lattice_structure     bcc
-C11                   523.0e9       # From Marinica et al. Journal of Physics: Condensed Matter(2013)
-C12                   202.0e9
-C44                   161.0e9
-
-grainsize             2.0e-5        # Average grain size [m]
-SolidSolutionStrength 1.5e8         # Strength due to elements in solid solution
-
-### Dislocation glide parameters ###
-#per family
-Nslip               12 
-slipburgers         2.72e-10        # Burgers vector of slip system [m]
-rhoedge0            1.0e12          # Initial edge dislocation density [m/m**3]
-rhoedgedip0         1.0             # Initial edged dipole dislocation density [m/m**3]
-v0                  1.0e-4          # Initial glide velocity [m/s]
-Qedge               2.725e-19       # Activation energy for dislocation glide [J]
-p_slip              0.78            # p-exponent in glide velocity
-q_slip              1.58            # q-exponent in glide velocity
-tau_peierls         2.03e9          # peierls stress (for bcc)
-dipoleformationfactor 0             # to have hardening due to dipole formation off
-
-#hardening
-CLambdaSlip         10.0            # Adj. parameter controlling dislocation mean free path
-D0                  4.0e-5          # Vacancy diffusion prefactor [m**2/s]
-Qsd                 4.5e-19         # Activation energy for climb [J]
-Catomicvolume       1.0             # Adj. parameter controlling the atomic volume [in b]
-Cedgedipmindistance 1.0             # Adj. parameter controlling the minimum dipole distance [in b]
-interaction_slipslip    1 1 1.4 1.4 1.4 1.4
-
-

examples/ConfigFiles/Phase_Dislotwin_Tungsten.yaml

@@ -0,0 +1,21 @@
+Tungsten:
+  lattice: cI
+  mechanics:
+    elasticity: {type: hooke, C_11: 523.0e9, C_12: 202.0e9, C_44: 161.0e9} # Marinica et al. Journal of Physics: Condensed Matter(2013)
+    plasticity:
+      type: dislotwin
+      D: 2.0e-5              # Average grain size / m
+      N_sl: [12]
+      b_sl: [2.72e-10]       # Burgers vector length of slip families / m
+      rho_mob_0: [1.0e12]
+      rho_dip_0: [1.0]
+      v_0: [1.0e4]           # Initial glide velocity / m/s
+      Q_s: [2.725e-19]       # Activation energy for dislocation glide / J
+      p_sl: [0.78]           # p-exponent in glide velocity
+      q_sl: [1.58]           # q-exponent in glide velocity
+      tau_0: [1.5e8]         # solid solution strength / Pa
+      i_sl: [10.0]           # Adj. parameter controlling dislocation mean free path
+      D_0: 4.0e-5            # Vacancy diffusion prefactor / m^2/s
+      D_a: 1.0               # minimum dipole distance / b
+      Q_cl: 4.5e-19          # Activation energy for climb / J
+      h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4]

examples/ConfigFiles/Phase_Hydrogen.config

@@ -1,3 +0,0 @@
-hydrogenflux_diffusion11        1.0
-hydrogenflux_mobility11         1.0
-hydrogenVolume                  1e-28

examples/ConfigFiles/Phase_Phenopowerlaw_cpTi.yaml

@@ -1,13 +1,14 @@
 # M. Levy, Handbook of Elastic Properties of Solids, Liquids, and Gases (2001)
 # C. Zambaldi, "Orientation informed nanoindentation of a-titanium: Indentation pileup in hexagonal metals deforming by prismatic slip", J. Mater. Res., Vol. 27, No. 1, Jan 14, 2012
-Ti-alpha:
+# Better use values from L. Wang, Z. Zheng, H. Phukan, P. Kenesei, J.-S. Park, J. Lind, R.M. Suter, T.R. Bieler, Direct measurement of critical resolved shear stress of prismatic and basal slip in polycrystalline Ti using high energy X-ray diffraction microscopy, Acta Mater 2017
+cpTi:
   lattice: hP
   c/a: 1.587
   mechanics:
     output: [F, P, F_e, F_p, L_p, O]
     elasticity: {C_11: 160.0e9, C_12: 90.0e9, C_13: 66.0e9, C_33: 181.7e9, C_44: 46.5e9, type: hooke}
     plasticity:
-      N_sl: [3, 3, 0, 0, 12]
+      N_sl: [3, 3, 0, 6, 12]
       a_sl: 2.0
       dot_gamma_0_sl: 0.001
       h_0_sl_sl: 200e6
@@ -15,5 +16,5 @@ Ti-alpha:
       n_sl: 20
       output: [gamma_sl]
       type: phenopowerlaw
-      xi_0_sl: [349e6, 150e6, 0, 0, 1107e6]
+      xi_0_sl: [0.15e9, 0.09e9, 0, 0.20e9, 0.25e9]
-      xi_inf_sl: [568e6, 1502e6, 0, 0, 3420e6]
+      xi_inf_sl: [0.24e9, 0.5e9, 0, 0.6e9, 0.8e9]

python/damask/_config.py

@@ -1,5 +1,6 @@
 import copy
 from io import StringIO
+from collections.abc import Iterable
 import abc
 
 import numpy as np
@@ -44,6 +45,42 @@ class Config(dict):
     copy = __copy__
 
 
+    def __or__(self,other):
+        """
+        Update configuration with contents of other.
+
+        Parameters
+        ----------
+        other : damask.Config or dict
+            Key-value pairs that update self.
+
+        """
+        duplicate = self.copy()
+        duplicate.update(other)
+        return duplicate
+
+
+    def __ior__(self,other):
+        """Update configuration with contents of other."""
+        return self.__or__(other)
+
+
+    def delete(self,keys):
+        """
+        Remove configuration keys.
+
+        Parameters
+        ----------
+        keys : iterable or scalar
+            Label of the key(s) to remove.
+
+        """
+        duplicate = self.copy()
+        for k in keys if isinstance(keys, Iterable) and not isinstance(keys, str) else [keys]:
+            del duplicate[k]
+        return duplicate
+
+
     @classmethod
     def load(cls,fname):
         """
@@ -99,30 +136,6 @@ class Config(dict):
             fhandle.write(yaml.dump(self,Dumper=NiceDumper,**kwargs))
 
 
-    def add(self,d):
-        """
-        Add dictionary.
-
-        d : dict
-            Dictionary to append.
-        """
-        duplicate = self.copy()
-        duplicate.update(d)
-        return duplicate
-
-
-    def delete(self,key):
-        """
-        Delete item.
-
-        key : str or scalar
-            Label of the key to remove.
-        """
-        duplicate = self.copy()
-        del duplicate[key]
-        return duplicate
-
-
     @property
     @abc.abstractmethod
     def is_complete(self):

python/damask/_result.py

@@ -105,7 +105,7 @@ class Result:
         self.view('increments',visible_increments)
 
         in_between = '' if len(visible_increments) < 3 else \
-                     ''.join([f'\n{inc}\n  ...\n' for inc in visible_increments[1:-2]])
+                     ''.join([f'\n{inc}\n  ...\n' for inc in visible_increments[1:-1]])
 
         return util.srepr(first + in_between + last)
 
@@ -119,10 +119,10 @@ class Result:
         action : str
             Select from 'set', 'add', and 'del'.
         what : str
-            Attribute to change (must be in self.visible).
+            Attribute to change (must be from self.visible).
         datasets : list of str or bool
-            Name of datasets as list, supports ? and * wildcards.
+            Name of datasets as list; supports ? and * wildcards.
-            True is equivalent to [*], False is equivalent to []
+            True is equivalent to [*], False is equivalent to [].
 
         """
         def natural_sort(key):
@@ -200,7 +200,7 @@ class Result:
         self._allow_modification = True
 
     def disallow_modification(self):
-        """Disllow to overwrite existing data (default case)."""
+        """Disallow to overwrite existing data (default case)."""
         self._allow_modification = False
 
 
@@ -272,10 +272,10 @@ class Result:
         Parameters
         ----------
         what : str
-            attribute to change (must be from self.visible)
+            Attribute to change (must be from self.visible).
         datasets : list of str or bool
-            name of datasets as list, supports ? and * wildcards.
+            Name of datasets as list; supports ? and * wildcards.
-            True is equivalent to [*], False is equivalent to []
+            True is equivalent to [*], False is equivalent to [].
 
         """
         self._manage_view('set',what,datasets)
@@ -288,10 +288,10 @@ class Result:
         Parameters
         ----------
         what : str
-            attribute to change (must be from self.visible)
+            Attribute to change (must be from self.visible).
         datasets : list of str or bool
-            name of datasets as list, supports ? and * wildcards.
+            Name of datasets as list; supports ? and * wildcards.
-            True is equivalent to [*], False is equivalent to []
+            True is equivalent to [*], False is equivalent to [].
 
         """
         self._manage_view('add',what,datasets)
@@ -304,10 +304,10 @@ class Result:
         Parameters
         ----------
         what : str
-            attribute to change (must be from self.visible)
+            Attribute to change (must be from self.visible).
         datasets : list of str or bool
-            name of datasets as list, supports ? and * wildcards.
+            Name of datasets as list; supports ? and * wildcards.
-            True is equivalent to [*], False is equivalent to []
+            True is equivalent to [*], False is equivalent to [].
 
         """
         self._manage_view('del',what,datasets)
@@ -315,14 +315,14 @@ class Result:
 
     def rename(self,name_old,name_new):
         """
-        Rename datasets.
+        Rename dataset.
 
         Parameters
         ----------
         name_old : str
-            name of the datasets to be renamed
+            Name of the dataset to be renamed.
         name_new : str
-            new name of the datasets
+            New name of the dataset.
 
         """
         if self._allow_modification:
@@ -353,13 +353,13 @@ class Result:
         ----------
           datasets : iterable or str
           constituent : int
-              Constituent to consider for phase data
+              Constituent to consider for phase data.
           tagged : bool
-              tag Table.column name with '#constituent'
+              Tag Table.column name with '#constituent'.
-              defaults to False
+              Defaults to False.
           split : bool
-              split Table by increment and return dictionary of Tables
+              Split Table by increment and return dictionary of Tables.
-              defaults to True
+              Defaults to True.
 
         """
         sets = datasets if hasattr(datasets,'__iter__') and not isinstance(datasets,str) else \
@@ -371,7 +371,7 @@ class Result:
         with h5py.File(self.fname,'r') as f:
             for dataset in sets:
                 for group in self.groups_with_datasets(dataset):
-                    path = os.path.join(group,dataset)
+                    path = '/'.join([group,dataset])
                     inc,prop,name,cat,item = (path.split('/') + ['']*5)[:5]
                     key = '/'.join([prop,name+tag])
                     if key not in inGeom:
@@ -388,15 +388,15 @@ class Result:
                                    np.nan,
                                    dtype=np.dtype(f[path]))
                     data[inGeom[key]] = (f[path] if len(shape)>1 else np.expand_dims(f[path],1))[inData[key]]
-                    path = (os.path.join(*([prop,name]+([cat] if cat else [])+([item] if item else []))) if split else path)+tag
+                    path = ('/'.join([prop,name]+([cat] if cat else [])+([item] if item else [])) if split else path)+tag
                     if split:
                         try:
-                            tbl[inc].add(path,data)
+                            tbl[inc] = tbl[inc].add(path,data)
                         except KeyError:
                             tbl[inc] = Table(data.reshape(self.N_materialpoints,-1),{path:data.shape[1:]})
                     else:
                         try:
-                            tbl.add(path,data)
+                            tbl = tbl.add(path,data)
                         except AttributeError:
                             tbl = Table(data.reshape(self.N_materialpoints,-1),{path:data.shape[1:]})
 
@@ -415,7 +415,7 @@ class Result:
         are considered as they contain user-relevant data.
         Single strings will be treated as list with one entry.
 
-        Wild card matching is allowed, but the number of arguments need to fit.
+        Wild card matching is allowed, but the number of arguments needs to fit.
 
         Parameters
         ----------

python/damask/_rotation.py

@@ -475,15 +475,17 @@ class Rotation:
 
         Parameters
         ----------
-        vector : bool, optional
+        compact : bool, optional
-            Return as actual Rodrigues-Frank vector, i.e. axis
+            Return as actual Rodrigues-Frank vector,
-            and angle argument are not separated.
+            i.e. axis and angle argument are not separated.
 
         Returns
         -------
-        rho : numpy.ndarray of shape (...,4) unless vector == True:
+        rho : numpy.ndarray of shape (...,4) containing 
-            numpy.ndarray of shape (...,3)
+            [n_1, n_2, n_3, tan(ω/2)], ǀnǀ = 1 and ω ∈ [0,π]
-            Rodrigues-Frank vector: [n_1, n_2, n_3, tan(ω/2)], ǀnǀ = 1 and ω ∈ [0,π].
+            unless compact == True:
+            numpy.ndarray of shape (...,3) containing 
+            tan(ω/2) [n_1, n_2, n_3], ω ∈ [0,π].
 
         """
         ro = Rotation._qu2ro(self.quaternion)

python/damask/_table.py

@@ -189,6 +189,11 @@ class Table:
         label : str
             Column label.
 
+        Returns
+        -------
+        data : numpy.ndarray
+            Array of column data.
+
         """
         if re.match(r'[0-9]*?_',label):
             idx,key = label.split('_',1)
@@ -212,6 +217,11 @@ class Table:
         info : str, optional
             Human-readable information about the new data.
 
+        Returns
+        -------
+        table : Table
+            Updated table.
+
         """
         dup = self.copy()
         dup._add_comment(label,data.shape[1:],info)
@@ -238,6 +248,11 @@ class Table:
         info : str, optional
             Human-readable information about the modified data.
 
+        Returns
+        -------
+        table : Table
+            Updated table.
+
         """
         dup = self.copy()
         dup._add_comment(label,data.shape[1:],info)
@@ -261,6 +276,11 @@ class Table:
         label : str
             Column label.
 
+        Returns
+        -------
+        table : Table
+            Updated table.
+
         """
         dup = self.copy()
         dup.data.drop(columns=label,inplace=True)
@@ -279,6 +299,11 @@ class Table:
         label_new : str or iterable of str
             New column label(s).
 
+        Returns
+        -------
+        table : Table
+            Updated table.
+
         """
         dup = self.copy()
         columns = dict(zip([old] if isinstance(old,str) else old,
@@ -300,6 +325,11 @@ class Table:
         ascending : bool or list, optional
             Set sort order.
 
+        Returns
+        -------
+        table : Table
+            Updated table.
+
         """
         dup = self.copy()
         dup._label_discrete()
@@ -320,6 +350,11 @@ class Table:
         other : Table
             Table to append.
 
+        Returns
+        -------
+        table : Table
+            Concatenated table.
+
         """
         if self.shapes != other.shapes or not self.data.columns.equals(other.data.columns):
             raise KeyError('Labels or shapes or order do not match')
@@ -340,6 +375,11 @@ class Table:
         other : Table
             Table to join.
 
+        Returns
+        -------
+        table : Table
+            Joined table.
+
         """
         if set(self.shapes) & set(other.shapes) or self.data.shape[0] != other.data.shape[0]:
             raise KeyError('Dublicated keys or row count mismatch')

python/damask/util.py

@@ -133,6 +133,8 @@ def execute(cmd,
     stdout = stdout.decode('utf-8').replace('\x08','')
     stderr = stderr.decode('utf-8').replace('\x08','')
     if process.returncode != 0:
+        print(stdout)
+        print(stderr)
         raise RuntimeError(f"'{cmd}' failed with returncode {process.returncode}")
     return stdout, stderr
 
@@ -193,7 +195,7 @@ def scale_to_coprime(v):
     return m
 
 
-def project_stereographic(vector,normalize=False):
+def project_stereographic(vector,direction='z',normalize=True,keepdims=False):
     """
     Apply stereographic projection to vector.
 
@@ -201,18 +203,37 @@ def project_stereographic(vector,normalize=False):
     ----------
     vector : numpy.ndarray of shape (...,3)
         Vector coordinates to be projected.
+    direction : str
+        Projection direction 'x', 'y', or 'z'.
+        Defaults to 'z'.
     normalize : bool
-        Ensure unit length for vector. Defaults to False.
+        Ensure unit length of input vector. Defaults to True.
+    keepdims : bool
+        Maintain three-dimensional output coordinates.
+        Default two-dimensional output uses right-handed frame spanned by
+        the next and next-next axis relative to the projection direction,
+        e.g. x-y when projecting along z and z-x when projecting along y.
 
     Returns
     -------
-    coordinates : numpy.ndarray of shape (...,2)
+    coordinates : numpy.ndarray of shape (...,2 | 3)
         Projected coordinates.
 
+    Examples
+    --------
+    >>> project_stereographic(np.ones(3))
+        [0.3660254, 0.3660254]
+    >>> project_stereographic(np.ones(3),direction='x',normalize=False,keepdims=True)
+        [0, 0.5, 0.5]
+    >>> project_stereographic([0,1,1],direction='y',normalize=True,keepdims=False)
+        [0.41421356, 0]
+
     """
-    v_ = vector/np.linalg.norm(vector,axis=-1,keepdims=True) if normalize else vector
+    shift = 'zyx'.index(direction)
-    return np.block([v_[...,:2]/(1+np.abs(v_[...,2:3])),
+    v_ = np.roll(vector/np.linalg.norm(vector,axis=-1,keepdims=True) if normalize else vector,
-                     np.zeros_like(v_[...,2:3])])
+                 shift,axis=-1)
+    return np.roll(np.block([v_[...,:2]/(1+np.abs(v_[...,2:3])),np.zeros_like(v_[...,2:3])]),
+                   -shift if keepdims else 0,axis=-1)[...,:3 if keepdims else 2]
 
 
 def execution_stamp(class_name,function_name=None):
@@ -418,7 +439,7 @@ class _ProgressBar:
             bar = '█' * filled_length + '░' * (self.bar_length - filled_length)
             delta_time = datetime.datetime.now() - self.start_time
             remaining_time = (self.total - (iteration+1)) * delta_time / (iteration+1)
-            remaining_time -= datetime.timedelta(microseconds=remaining_time.microseconds)           # remove μs
+            remaining_time -= datetime.timedelta(microseconds=remaining_time.microseconds)          # remove μs
             sys.stderr.write(f'\r{self.prefix} {bar} {fraction:>4.0%} ETA {remaining_time}')
             sys.stderr.flush()
 

python/tests/test_Config.py

@@ -23,8 +23,17 @@ class TestConfig:
             assert Config.load(f) == config
 
     def test_add_remove(self):
+        dummy = {'hello':'world','foo':'bar'}
         config = Config()
-        assert config.add({'hello':'world'}).delete('hello') == config
+        config |= dummy
+        assert config == Config() | dummy
+        config = config.delete(dummy)
+        assert config == Config()
+        assert (config |        dummy ).delete(        'hello'            ) == config | {'foo':'bar'}
+        assert (config |        dummy ).delete([       'hello',  'foo'   ]) == config
+        assert (config | Config(dummy)).delete({       'hello':1,'foo':2 }) == config
+        assert (config | Config(dummy)).delete(Config({'hello':1        })) == config | {'foo':'bar'}
+
 
     def test_repr(self,tmp_path):
         config = Config()

python/tests/test_util.py

@@ -49,17 +49,18 @@ class TestUtil:
         dist_sampled = np.histogram(centers[selected],bins)[0]/N_samples*np.sum(dist)
         assert np.sqrt(((dist - dist_sampled) ** 2).mean()) < .025 and selected.shape[0]==N_samples
 
-    @pytest.mark.parametrize('point,normalize,answer',
+    @pytest.mark.parametrize('point,direction,normalize,keepdims,answer',
                              [
-                              ([1,0,0],False,[1,0,0]),
+                              ([1,0,0],'z',False,True, [1,0,0]),
-                              ([1,0,0],True, [1,0,0]),
+                              ([1,0,0],'z',True, False,[1,0]),
-                              ([0,1,1],False,[0,0.5,0]),
+                              ([0,1,1],'z',False,True, [0,0.5,0]),
-                              ([0,1,1],True, [0,0.41421356,0]),
+                              ([0,1,1],'y',True, False,[0.41421356,0]),
-                              ([1,1,1],False,[0.5,0.5,0]),
+                              ([1,1,0],'x',False,False,[0.5,0]),
-                              ([1,1,1],True, [0.3660254, 0.3660254, 0]),
+                              ([1,1,1],'y',True, True, [0.3660254, 0,0.3660254]),
                              ])
-    def test_project_stereographic(self,point,normalize,answer):
+    def test_project_stereographic(self,point,direction,normalize,keepdims,answer):
-        assert np.allclose(util.project_stereographic(np.array(point),normalize=normalize),answer)
+        assert np.allclose(util.project_stereographic(np.array(point),direction=direction,
+                                                      normalize=normalize,keepdims=keepdims),answer)
 
     @pytest.mark.parametrize('fro,to,mode,answer',
                              [

src/CPFEM.f90

@@ -178,11 +178,11 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
 
   if (iand(mode, CPFEM_AGERESULTS) /= 0_pInt) call CPFEM_forward
 
-    chosenThermal1: select case (thermal_type(material_homogenizationAt(elCP)))
+    !chosenThermal1: select case (thermal_type(material_homogenizationAt(elCP)))
-      !case (THERMAL_conduction_ID) chosenThermal1
+    !  case (THERMAL_conduction_ID) chosenThermal1
-      !  temperature(material_homogenizationAt(elCP))%p(material_homogenizationMemberAt(ip,elCP)) = &
+    !    temperature(material_homogenizationAt(elCP))%p(material_homogenizationMemberAt(ip,elCP)) = &
-      !    temperature_inp
+    !      temperature_inp
-    end select chosenThermal1
+    !end select chosenThermal1
     homogenization_F0(1:3,1:3,ma) = ffn
     homogenization_F(1:3,1:3,ma) = ffn1
 

src/YAML_types.f90

@@ -72,7 +72,7 @@ module YAML_types
       getKey                      => tNode_getKey_byIndex
     procedure :: &
       contains                    => tNode_contains
- 
+
     generic :: &
       get           => tNode_get_byIndex, &
                        tNode_get_byKey
@@ -157,7 +157,7 @@ module YAML_types
     emptyDict
   type(tList), target, public :: &
     emptyList
-  
+
   abstract interface
 
     recursive function asFormattedString(self,indent)
@@ -179,7 +179,7 @@ module YAML_types
 
   public :: &
     YAML_types_init, &
-    output_asStrings, &                                        !ToDo: Hack for GNU. Remove later 
+    output_asStrings, &                                        !ToDo: Hack for GNU. Remove later
     assignment(=)
 
 contains
@@ -207,11 +207,11 @@ subroutine selfTest
   select type(s1)
     class is(tScalar)
       s1 = '1'
-      if(s1%asInt() /= 1)              error stop 'tScalar_asInt'
+      if (s1%asInt() /= 1)              error stop 'tScalar_asInt'
-      if(dNeq(s1%asFloat(),1.0_pReal)) error stop 'tScalar_asFloat'
+      if (dNeq(s1%asFloat(),1.0_pReal)) error stop 'tScalar_asFloat'
       s1 = 'true'
-      if(.not. s1%asBool())            error stop 'tScalar_asBool'
+      if (.not. s1%asBool())            error stop 'tScalar_asBool'
-      if(s1%asString() /= 'true')      error stop 'tScalar_asString'
+      if (s1%asString() /= 'true')      error stop 'tScalar_asString'
   end select
 
   block
@@ -232,18 +232,18 @@ subroutine selfTest
         call l1%append(s1)
         call l1%append(s2)
         n => l1
-        if(any(l1%asInts() /= [2,3]))                      error stop 'tList_asInts'
+        if (any(l1%asInts() /= [2,3]))                      error stop 'tList_asInts'
-        if(any(dNeq(l1%asFloats(),[2.0_pReal,3.0_pReal]))) error stop 'tList_asFloats'
+        if (any(dNeq(l1%asFloats(),[2.0_pReal,3.0_pReal]))) error stop 'tList_asFloats'
-        if(n%get_asInt(1) /= 2)                            error stop 'byIndex_asInt'
+        if (n%get_asInt(1) /= 2)                            error stop 'byIndex_asInt'
-        if(dNeq(n%get_asFloat(2),3.0_pReal))               error stop 'byIndex_asFloat'
+        if (dNeq(n%get_asFloat(2),3.0_pReal))               error stop 'byIndex_asFloat'
     endselect
 
     allocate(tList::l2)
     select type(l2)
       class is(tList)
         call l2%append(l1)
-        if(any(l2%get_asInts(1) /= [2,3]))                      error stop 'byIndex_asInts'
+        if (any(l2%get_asInts(1) /= [2,3]))                      error stop 'byIndex_asInts'
-        if(any(dNeq(l2%get_asFloats(1),[2.0_pReal,3.0_pReal]))) error stop 'byIndex_asFloats'
+        if (any(dNeq(l2%get_asFloats(1),[2.0_pReal,3.0_pReal]))) error stop 'byIndex_asFloats'
         n => l2
     end select
     deallocate(n)
@@ -265,10 +265,10 @@ subroutine selfTest
      call l1%append(s2)
      n => l1
 
-     if(any(l1%asBools() .neqv. [.true., .false.])) error stop 'tList_asBools'
+     if (any(l1%asBools() .neqv. [.true., .false.])) error stop 'tList_asBools'
-     if(any(l1%asStrings() /=   ['true ','False'])) error stop 'tList_asStrings'
+     if (any(l1%asStrings() /=   ['true ','False'])) error stop 'tList_asStrings'
-     if(n%get_asBool(2))                            error stop 'byIndex_asBool'
+     if (n%get_asBool(2))                            error stop 'byIndex_asBool'
-     if(n%get_asString(1) /= 'true')                error stop 'byIndex_asString'
+     if (n%get_asString(1) /= 'true')                error stop 'byIndex_asString'
   end block
 
 end subroutine selfTest
@@ -418,7 +418,7 @@ function tNode_get_byIndex(self,i) result(node)
   integer :: j
 
   self_ => self%asList()
-  if(i < 1 .or. i > self_%length) call IO_error(150,ext_msg='tNode_get_byIndex')
+  if (i < 1 .or. i > self_%length) call IO_error(150,ext_msg='tNode_get_byIndex')
 
   j = 1
   item => self_%first
@@ -599,7 +599,7 @@ function tNode_getKey_byIndex(self,i)  result(key)
   dict => self%asDict()
   item => dict%first
   do j = 1, dict%length
-    if(j == i) then
+    if (j == i) then
       key = item%key
       exit
     else
@@ -613,7 +613,7 @@ end function tNode_getKey_byIndex
 !-------------------------------------------------------------------------------------------------
 !> @brief Checks if a given key/item is present in the dict/list
 !-------------------------------------------------------------------------------------------------
-function tNode_contains(self,k)  result(exists)   
+function tNode_contains(self,k)  result(exists)
 
   class(tNode),     intent(in), target  :: self
   character(len=*), intent(in)          :: k
@@ -624,18 +624,18 @@ function tNode_contains(self,k)  result(exists)
   type(tDict), pointer :: dict
 
   exists = .false.
-  if(self%isDict()) then
+  if (self%isDict()) then
     dict => self%asDict()
     do j=1, dict%length
-      if(dict%getKey(j) == k) then
+      if (dict%getKey(j) == k) then
         exists = .true.
         return
       endif
     enddo
-  elseif(self%isList()) then
+  elseif (self%isList()) then
     list => self%asList()
-    do j =1, list%length
+    do j=1, list%length
-      if(list%get_asString(j) == k) then
+      if (list%get_asString(j) == k) then
         exists = .true.
         return
       endif
@@ -663,8 +663,8 @@ function tNode_get_byKey(self,k,defaultVal) result(node)
   logical :: found
 
   found = present(defaultVal)
-  if(found) node => defaultVal
+  if (found) node => defaultVal
-  
+
   self_ => self%asDict()
 
   j = 1
@@ -677,11 +677,11 @@ function tNode_get_byKey(self,k,defaultVal) result(node)
     item => item%next
     j = j + 1
   enddo
-  
+
   if (.not. found) then
     call IO_error(143,ext_msg=k)
   else
-    if(associated(item)) node => item%node
+    if (associated(item)) node => item%node
   endif
 
 end function tNode_get_byKey
@@ -700,11 +700,11 @@ function tNode_get_byKey_asFloat(self,k,defaultVal) result(nodeAsFloat)
   class(tNode),  pointer :: node
   type(tScalar), pointer :: scalar
 
-  if(self%contains(k)) then
+  if (self%contains(k)) then
     node => self%get(k)
     scalar => node%asScalar()
     nodeAsFloat = scalar%asFloat()
-  elseif(present(defaultVal)) then
+  elseif (present(defaultVal)) then
     nodeAsFloat = defaultVal
   else
     call IO_error(143,ext_msg=k)
@@ -726,11 +726,11 @@ function tNode_get_byKey_asInt(self,k,defaultVal) result(nodeAsInt)
   class(tNode),  pointer :: node
   type(tScalar), pointer :: scalar
 
-  if(self%contains(k)) then
+  if (self%contains(k)) then
     node => self%get(k)
     scalar => node%asScalar()
     nodeAsInt = scalar%asInt()
-  elseif(present(defaultVal)) then
+  elseif (present(defaultVal)) then
     nodeAsInt = defaultVal
   else
     call IO_error(143,ext_msg=k)
@@ -752,11 +752,11 @@ function tNode_get_byKey_asBool(self,k,defaultVal) result(nodeAsBool)
   class(tNode),  pointer :: node
   type(tScalar), pointer :: scalar
 
-  if(self%contains(k)) then
+  if (self%contains(k)) then
     node => self%get(k)
     scalar => node%asScalar()
     nodeAsBool = scalar%asBool()
-  elseif(present(defaultVal)) then
+  elseif (present(defaultVal)) then
     nodeAsBool = defaultVal
   else
     call IO_error(143,ext_msg=k)
@@ -778,11 +778,11 @@ function tNode_get_byKey_asString(self,k,defaultVal) result(nodeAsString)
   class(tNode),  pointer :: node
   type(tScalar), pointer :: scalar
 
-  if(self%contains(k)) then
+  if (self%contains(k)) then
     node => self%get(k)
     scalar => node%asScalar()
     nodeAsString = scalar%asString()
-  elseif(present(defaultVal)) then
+  elseif (present(defaultVal)) then
     nodeAsString = defaultVal
   else
     call IO_error(143,ext_msg=k)
@@ -806,18 +806,18 @@ function tNode_get_byKey_asFloats(self,k,defaultVal,requiredSize) result(nodeAsF
   class(tNode), pointer :: node
   type(tList),  pointer :: list
 
-  if(self%contains(k)) then
+  if (self%contains(k)) then
     node => self%get(k)
     list => node%asList()
     nodeAsFloats = list%asFloats()
-  elseif(present(defaultVal)) then
+  elseif (present(defaultVal)) then
     nodeAsFloats = defaultVal
   else
     call IO_error(143,ext_msg=k)
   endif
 
-  if(present(requiredSize)) then
+  if (present(requiredSize)) then
-    if(requiredSize /= size(nodeAsFloats)) call IO_error(146,ext_msg=k)
+    if (requiredSize /= size(nodeAsFloats)) call IO_error(146,ext_msg=k)
   endif
 
 end function tNode_get_byKey_asFloats
@@ -837,18 +837,18 @@ function tNode_get_byKey_asInts(self,k,defaultVal,requiredSize) result(nodeAsInt
   class(tNode), pointer :: node
   type(tList),  pointer :: list
 
-  if(self%contains(k)) then
+  if (self%contains(k)) then
     node => self%get(k)
     list => node%asList()
     nodeAsInts = list%asInts()
-  elseif(present(defaultVal)) then
+  elseif (present(defaultVal)) then
     nodeAsInts = defaultVal
   else
     call IO_error(143,ext_msg=k)
   endif
 
-  if(present(requiredSize)) then
+  if (present(requiredSize)) then
-    if(requiredSize /= size(nodeAsInts)) call IO_error(146,ext_msg=k)
+    if (requiredSize /= size(nodeAsInts)) call IO_error(146,ext_msg=k)
   endif
 
 end function tNode_get_byKey_asInts
@@ -867,11 +867,11 @@ function tNode_get_byKey_asBools(self,k,defaultVal) result(nodeAsBools)
   class(tNode), pointer :: node
   type(tList),  pointer :: list
 
-  if(self%contains(k)) then
+  if (self%contains(k)) then
     node => self%get(k)
     list => node%asList()
     nodeAsBools = list%asBools()
-  elseif(present(defaultVal)) then
+  elseif (present(defaultVal)) then
     nodeAsBools = defaultVal
   else
     call IO_error(143,ext_msg=k)
@@ -893,11 +893,11 @@ function tNode_get_byKey_asStrings(self,k,defaultVal) result(nodeAsStrings)
   class(tNode), pointer :: node
   type(tList),  pointer :: list
 
-  if(self%contains(k)) then
+  if (self%contains(k)) then
     node => self%get(k)
     list => node%asList()
     nodeAsStrings = list%asStrings()
-  elseif(present(defaultVal)) then
+  elseif (present(defaultVal)) then
     nodeAsStrings = defaultVal
   else
     call IO_error(143,ext_msg=k)
@@ -925,7 +925,7 @@ function output_asStrings(self)  result(output)                   !ToDo: SR: Rem
 
 
 end function output_asStrings
- 
+
 
 !--------------------------------------------------------------------------------------------------
 !> @brief Returns the index of a key in a dictionary
@@ -944,7 +944,7 @@ function tNode_get_byKey_asIndex(self,key)  result(keyIndex)
   item => dict%first
   keyIndex = -1
   do i = 1, dict%length
-    if(key == item%key) then
+    if (key == item%key) then
       keyIndex = i
       exit
     else
@@ -952,9 +952,9 @@ function tNode_get_byKey_asIndex(self,key)  result(keyIndex)
     endif
   enddo
 
-  if(keyIndex == -1) call IO_error(140,ext_msg=key)
+  if (keyIndex == -1) call IO_error(140,ext_msg=key)
+
 
- 
 end function tNode_get_byKey_asIndex
 
 
@@ -985,7 +985,7 @@ recursive function tList_asFormattedString(self,indent) result(str)
   integer :: i, indent_
 
   str = ''
-  if(present(indent)) then
+  if (present(indent)) then
     indent_ = indent
   else
     indent_ = 0
@@ -993,7 +993,7 @@ recursive function tList_asFormattedString(self,indent) result(str)
 
   item => self%first
   do i = 1, self%length
-    if(i /= 1) str = str//repeat(' ',indent_)
+    if (i /= 1) str = str//repeat(' ',indent_)
     str = str//'- '//item%node%asFormattedString(indent_+2)
     item => item%next
   end do
@@ -1014,7 +1014,7 @@ recursive function tDict_asFormattedString(self,indent) result(str)
   integer :: i, indent_
 
   str = ''
-  if(present(indent)) then
+  if (present(indent)) then
     indent_ = indent
   else
     indent_ = 0
@@ -1022,7 +1022,7 @@ recursive function tDict_asFormattedString(self,indent) result(str)
 
   item => self%first
   do i = 1, self%length
-    if(i /= 1) str = str//repeat(' ',indent_)
+    if (i /= 1) str = str//repeat(' ',indent_)
     select type(node_1 =>item%node)
       class is(tScalar)
         str = str//trim(item%key)//': '//item%node%asFormattedString(indent_+len_trim(item%key)+2)
@@ -1270,7 +1270,7 @@ recursive subroutine tItem_finalize(self)
   type(tItem),intent(inout) :: self
 
   deallocate(self%node)
-  if(associated(self%next)) deallocate(self%next)
+  if (associated(self%next)) deallocate(self%next)
 
 end subroutine tItem_finalize
 

src/grid/DAMASK_grid.f90

@@ -37,9 +37,10 @@ program DAMASK_grid
                                 f_out, &                                                            !< frequency of result writes
                                 f_restart                                                           !< frequency of restart writes
     logical ::                  estimate_rate                                                       !< follow trajectory of former loadcase
-    integer(kind(FIELD_UNDEFINED_ID)), allocatable :: ID(:)
   end type tLoadCase
 
+  integer(kind(FIELD_UNDEFINED_ID)), allocatable :: ID(:)
+
 !--------------------------------------------------------------------------------------------------
 ! variables related to information from load case and geom file
   real(pReal), dimension(9) :: temp_valueVector                                                     !< temporarily from loadcase file when reading in tensors (initialize to 0.0)
@@ -53,6 +54,7 @@ program DAMASK_grid
   integer, parameter :: &
     subStepFactor = 2                                                                               !< for each substep, divide the last time increment by 2.0
   real(pReal) :: &
+    T_0 = 300.0_pReal, &
     time = 0.0_pReal, &                                                                             !< elapsed time
     time0 = 0.0_pReal, &                                                                            !< begin of interval
     timeinc = 1.0_pReal, &                                                                          !< current time interval
@@ -78,8 +80,7 @@ program DAMASK_grid
     maxCutBack, &                                                                                   !< max number of cut backs
     stagItMax                                                                                       !< max number of field level staggered iterations
   character(len=pStringLen) :: &
-    incInfo, &
+    incInfo
-    loadcase_string
 
   type(tLoadCase), allocatable, dimension(:) :: loadCases                                           !< array of all load cases
   type(tSolutionState), allocatable, dimension(:) :: solres
@@ -98,10 +99,13 @@ program DAMASK_grid
     quit
   class (tNode), pointer :: &
     num_grid, &
-    debug_grid, &                                                                                   ! pointer to grid debug options
     config_load, &
     load_steps, &
     load_step, &
+    solver, &
+    initial_conditions, &
+    ic_thermal, &
+    thermal, &
     step_bc, &
     step_mech, &
     step_discretization, &
@@ -112,17 +116,11 @@ program DAMASK_grid
 ! init DAMASK (all modules)
 
   call CPFEM_initAll
-  print'(/,a)',   ' <<<+-  DAMASK_spectral init  -+>>>'; flush(IO_STDOUT)
+  print'(/,a)',   ' <<<+-  DAMASK_grid init  -+>>>'; flush(IO_STDOUT)
 
   print*, 'Shanthraj et al., Handbook of Mechanics of Materials, 2019'
   print*, 'https://doi.org/10.1007/978-981-10-6855-3_80'
 
-!--------------------------------------------------------------------------------------------------
-! initialize field solver information
-  nActiveFields = 1
-  if (any(thermal_type  == THERMAL_conduction_ID  )) nActiveFields = nActiveFields + 1
-  if (any(damage_type   == DAMAGE_nonlocal_ID     )) nActiveFields = nActiveFields + 1
-  allocate(solres(nActiveFields))
 
 !-------------------------------------------------------------------------------------------------
 ! reading field paramters from numerics file and do sanity checks
@@ -133,19 +131,22 @@ program DAMASK_grid
   if (stagItMax < 0)    call IO_error(301,ext_msg='maxStaggeredIter')
   if (maxCutBack < 0)   call IO_error(301,ext_msg='maxCutBack')
 
+  config_load => YAML_parse_file(trim(interface_loadFile))
+  solver => config_load%get('solver')
+
 !--------------------------------------------------------------------------------------------------
 ! assign mechanics solver depending on selected type
 
-  debug_grid => config_debug%get('grid',defaultVal=emptyList)
+  nActiveFields = 1
-  select case (trim(num_grid%get_asString('solver', defaultVal = 'Basic')))
+  select case (solver%get_asString('mechanical'))
-    case ('Basic')
+    case ('spectral_basic')
       mechanical_init         => grid_mechanical_spectral_basic_init
       mechanical_forward      => grid_mechanical_spectral_basic_forward
       mechanical_solution     => grid_mechanical_spectral_basic_solution
       mechanical_updateCoords => grid_mechanical_spectral_basic_updateCoords
       mechanical_restartWrite => grid_mechanical_spectral_basic_restartWrite
 
-    case ('Polarisation')
+    case ('spectral_polarization')
       mechanical_init         => grid_mechanical_spectral_polarisation_init
       mechanical_forward      => grid_mechanical_spectral_polarisation_forward
       mechanical_solution     => grid_mechanical_spectral_polarisation_solution
@@ -160,32 +161,36 @@ program DAMASK_grid
       mechanical_restartWrite => grid_mechanical_FEM_restartWrite
 
     case default
-      call IO_error(error_ID = 891, ext_msg = trim(num_grid%get_asString('solver')))
+      call IO_error(error_ID = 891, ext_msg = trim(solver%get_asString('mechanical')))
 
   end select
 
+!--------------------------------------------------------------------------------------------------
+! initialize field solver information
+  if (solver%get_asString('thermal',defaultVal = 'n/a') == 'spectral') nActiveFields = nActiveFields + 1
+  if (solver%get_asString('damage', defaultVal = 'n/a') == 'spectral') nActiveFields = nActiveFields + 1
+
+  allocate(solres(nActiveFields))
+  allocate(    ID(nActiveFields))
+
+  field = 1
+  ID(field) = FIELD_MECH_ID                                                                         ! mechanical active by default
+  thermalActive: if (solver%get_asString('thermal',defaultVal = 'n/a') == 'spectral') then
+    field = field + 1
+    ID(field) = FIELD_THERMAL_ID
+  endif thermalActive
+  damageActive: if (solver%get_asString('damage',defaultVal = 'n/a') == 'spectral') then
+    field = field + 1
+    ID(field) = FIELD_DAMAGE_ID
+  endif damageActive
+
 
 !--------------------------------------------------------------------------------------------------
-! reading information from load case file and to sanity checks
-  config_load => YAML_parse_file(trim(interface_loadFile))
-
   load_steps => config_load%get('loadstep')
   allocate(loadCases(load_steps%length))                                                            ! array of load cases
 
   do l = 1, load_steps%length
 
-    allocate(loadCases(l)%ID(nActiveFields))
-    field = 1
-    loadCases(l)%ID(field) = FIELD_MECH_ID                                                          ! mechanical active by default
-    thermalActive: if (any(thermal_type == THERMAL_conduction_ID)) then
-      field = field + 1
-      loadCases(l)%ID(field) = FIELD_THERMAL_ID
-    endif thermalActive
-    damageActive: if (any(damage_type == DAMAGE_nonlocal_ID)) then
-      field = field + 1
-      loadCases(l)%ID(field) = FIELD_DAMAGE_ID
-    endif damageActive
-
     load_step => load_steps%get(l)
     step_bc   => load_step%get('boundary_conditions')
     step_mech => step_bc%get('mechanical')
@@ -220,19 +225,17 @@ program DAMASK_grid
     if (.not. allocated(loadCases(l)%deformation%myType)) call IO_error(error_ID=837,ext_msg = 'L/dot_F/F missing')
 
     step_discretization => load_step%get('discretization')
-    if(.not. step_discretization%contains('t')) call IO_error(error_ID=837,ext_msg = 't missing')
+    if (.not. step_discretization%contains('t')) call IO_error(error_ID=837,ext_msg = 't missing')
-    if(.not. step_discretization%contains('N')) call IO_error(error_ID=837,ext_msg = 'N missing')
+    if (.not. step_discretization%contains('N')) call IO_error(error_ID=837,ext_msg = 'N missing')
     loadCases(l)%t         = step_discretization%get_asFloat('t')
     loadCases(l)%N         = step_discretization%get_asInt  ('N')
     loadCases(l)%r         = step_discretization%get_asFloat('r',         defaultVal= 1.0_pReal)
 
     loadCases(l)%f_restart = load_step%get_asInt('f_restart', defaultVal=huge(0))
     loadCases(l)%f_out     = load_step%get_asInt('f_out',     defaultVal=1)
-    loadCases(l)%estimate_rate = (load_step%get_asBool('estimate_rate',defaultVal=.true.) .and. & 
+    loadCases(l)%estimate_rate = (load_step%get_asBool('estimate_rate',defaultVal=.true.) .and. l>1)
-                                                       merge(.true.,.false.,l > 1))
 
     reportAndCheck: if (worldrank == 0) then
-      write (loadcase_string, '(i0)' ) l
       print'(/,a,i0)', ' load case: ', l
       print*, ' estimate_rate:', loadCases(l)%estimate_rate
       if (loadCases(l)%deformation%myType == 'L') then
@@ -286,13 +289,13 @@ program DAMASK_grid
       else
         print'(a,f0.3)', '  r: ', loadCases(l)%r
       endif
-      print'(a,f0.3)', '  t: ', loadCases(l)%t
+      print'(a,f0.3)',   '  t: ', loadCases(l)%t
-      print'(a,i0)',   '  N: ', loadCases(l)%N
+      print'(a,i0)',     '  N: ', loadCases(l)%N
-      print'(a,i0)',   '  f_out: ', loadCases(l)%f_out
+      print'(a,i0)',     '  f_out: ', loadCases(l)%f_out
       if (loadCases(l)%f_restart < huge(0)) &
-        print'(a,i0)', '  f_restart: ', loadCases(l)%f_restart
+        print'(a,i0)',   '  f_restart: ', loadCases(l)%f_restart
 
-      if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string)                 ! exit with error message
+      if (errorID > 0) call IO_error(error_ID = errorID, el = l)
 
     endif reportAndCheck
   enddo
@@ -301,12 +304,14 @@ program DAMASK_grid
 ! doing initialization depending on active solvers
   call spectral_Utilities_init
   do field = 1, nActiveFields
-    select case (loadCases(1)%ID(field))
+    select case (ID(field))
       case(FIELD_MECH_ID)
         call mechanical_init
 
       case(FIELD_THERMAL_ID)
-        call grid_thermal_spectral_init
+        initial_conditions => config_load%get('initial_conditions',defaultVal=emptyDict)
+        thermal => initial_conditions%get('thermal',defaultVal=emptyDict)
+        call grid_thermal_spectral_init(thermal%get_asFloat('T',defaultVal = T_0))
 
       case(FIELD_DAMAGE_ID)
         call grid_damage_spectral_init
@@ -377,7 +382,7 @@ program DAMASK_grid
 !--------------------------------------------------------------------------------------------------
 ! forward fields
           do field = 1, nActiveFields
-            select case(loadCases(l)%ID(field))
+            select case(ID(field))
               case(FIELD_MECH_ID)
                 call mechanical_forward (&
                         cutBack,guess,timeinc,timeIncOld,remainingLoadCaseTime, &
@@ -397,7 +402,7 @@ program DAMASK_grid
           stagIterate = .true.
           do while (stagIterate)
             do field = 1, nActiveFields
-              select case(loadCases(l)%ID(field))
+              select case(ID(field))
                 case(FIELD_MECH_ID)
                   solres(field) = mechanical_solution(incInfo)
                 case(FIELD_THERMAL_ID)

src/grid/grid_damage_spectral.f90

@@ -160,7 +160,7 @@ function grid_damage_spectral_solution(timeinc) result(solution)
 
   real(pReal), intent(in) :: &
     timeinc                                                                                         !< increment in time for current solution
-  integer :: i, j, k, cell
+  integer :: i, j, k, ce
   type(tSolutionState) :: solution
   PetscInt  :: devNull
   PetscReal :: phi_min, phi_max, stagNorm, solnNorm
@@ -193,10 +193,10 @@ function grid_damage_spectral_solution(timeinc) result(solution)
 
 !--------------------------------------------------------------------------------------------------
 ! updating damage state
-  cell = 0
+  ce = 0
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
-    cell = cell + 1
+    ce = ce + 1
-    call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),1,cell)
+    call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),ce)
   enddo; enddo; enddo
 
   call VecMin(solution_vec,devNull,phi_min,ierr); CHKERRQ(ierr)
@@ -216,7 +216,7 @@ end function grid_damage_spectral_solution
 subroutine grid_damage_spectral_forward(cutBack)
 
   logical, intent(in) :: cutBack
-  integer                               :: i, j, k, cell
+  integer                               :: i, j, k, ce
   DM :: dm_local
   PetscScalar,  dimension(:,:,:), pointer     :: x_scal
   PetscErrorCode                              :: ierr
@@ -226,14 +226,14 @@ subroutine grid_damage_spectral_forward(cutBack)
     phi_stagInc = phi_lastInc
 !--------------------------------------------------------------------------------------------------
 ! reverting damage field state
-    cell = 0
+    ce = 0
     call SNESGetDM(damage_snes,dm_local,ierr); CHKERRQ(ierr)
     call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr)                       !< get the data out of PETSc to work with
     x_scal(xstart:xend,ystart:yend,zstart:zend) = phi_current
     call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr)
     do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
-      cell = cell + 1
+      ce = ce + 1
-      call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),1,cell)
+      call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),ce)
     enddo; enddo; enddo
   else
     phi_lastInc = phi_current
@@ -258,7 +258,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
     f_scal
   PetscObject :: dummy
   PetscErrorCode :: ierr
-  integer :: i, j, k, cell
+  integer :: i, j, k, ce
   real(pReal)   :: phiDot, dPhiDot_dPhi, mobility
 
   phi_current = x_scal
@@ -269,20 +269,20 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
   call utilities_FFTscalarForward
   call utilities_fourierScalarGradient                                                              !< calculate gradient of damage field
   call utilities_FFTvectorBackward
-  cell = 0
+  ce = 0
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
-    cell = cell + 1
+    ce = ce + 1
-    vectorField_real(1:3,i,j,k) = matmul(damage_nonlocal_getDiffusion(1,cell) - K_ref, &
+    vectorField_real(1:3,i,j,k) = matmul(damage_nonlocal_getDiffusion(ce) - K_ref, &
                                          vectorField_real(1:3,i,j,k))
   enddo; enddo; enddo
   call utilities_FFTvectorForward
   call utilities_fourierVectorDivergence                                                            !< calculate damage divergence in fourier field
   call utilities_FFTscalarBackward
-  cell = 0
+  ce = 0
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
-    cell = cell + 1
+    ce = ce + 1
-    call damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi_current(i,j,k), 1, cell)
+    call damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi_current(i,j,k),ce)
-    mobility = damage_nonlocal_getMobility(1,cell)
+    mobility = damage_nonlocal_getMobility(ce)
     scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + phiDot) &
                             + mobility*(phi_lastInc(i,j,k) - phi_current(i,j,k)) &
                             + mu_ref*phi_current(i,j,k)
@@ -310,15 +310,15 @@ end subroutine formResidual
 !--------------------------------------------------------------------------------------------------
 subroutine updateReference
 
-  integer :: i,j,k,cell,ierr
+  integer :: i,j,k,ce,ierr
 
-  cell = 0
+  ce = 0
   K_ref = 0.0_pReal
   mu_ref = 0.0_pReal
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
-    cell = cell + 1
+    ce = ce + 1
-    K_ref  = K_ref  + damage_nonlocal_getDiffusion(1,cell)
+    K_ref  = K_ref  + damage_nonlocal_getDiffusion(ce)
-    mu_ref = mu_ref + damage_nonlocal_getMobility(1,cell)
+    mu_ref = mu_ref + damage_nonlocal_getMobility(ce)
   enddo; enddo; enddo
   K_ref = K_ref*wgt
   call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)

src/grid/grid_mech_spectral_polarisation.f90

@@ -116,7 +116,7 @@ subroutine grid_mechanical_spectral_polarisation_init
     num_grid, &
     debug_grid
 
-  print'(/,a)', ' <<<+-  grid_mechanical_spectral_polarisation init  -+>>>'; flush(IO_STDOUT)
+  print'(/,a)', ' <<<+-  grid_mechanical_spectral_polarization init  -+>>>'; flush(IO_STDOUT)
 
   print*, 'Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
   print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006'

src/grid/grid_thermal_spectral.f90

@@ -61,7 +61,9 @@ contains
 !> @brief allocates all neccessary fields and fills them with data
 ! ToDo: Restart not implemented
 !--------------------------------------------------------------------------------------------------
-subroutine grid_thermal_spectral_init
+subroutine grid_thermal_spectral_init(T_0)
+
+  real(pReal), intent(in) :: T_0
 
   PetscInt, dimension(0:worldsize-1) :: localK
   integer :: i, j, k, ce
@@ -131,9 +133,10 @@ subroutine grid_thermal_spectral_init
   ce = 0
   do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
     ce = ce + 1
-    T_current(i,j,k) = homogenization_thermal_T(ce)
+    T_current(i,j,k) = T_0
     T_lastInc(i,j,k) = T_current(i,j,k)
     T_stagInc(i,j,k) = T_current(i,j,k)
+    call homogenization_thermal_setField(T_0,0.0_pReal,ce)
   enddo; enddo; enddo
   call DMDAVecGetArrayF90(thermal_grid,solution_vec,x_scal,ierr); CHKERRQ(ierr)                     !< get the data out of PETSc to work with
   x_scal(xstart:xend,ystart:yend,zstart:zend) = T_current
@@ -268,7 +271,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
   ce = 0
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
     ce = ce + 1
-    vectorField_real(1:3,i,j,k) = matmul(thermal_conduction_getConductivity(1,ce) - K_ref, &
+    vectorField_real(1:3,i,j,k) = matmul(thermal_conduction_getConductivity(ce) - K_ref, &
                                          vectorField_real(1:3,i,j,k))
   enddo; enddo; enddo
   call utilities_FFTvectorForward
@@ -277,7 +280,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
   ce = 0
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
     ce = ce + 1
-    call thermal_conduction_getSource(Tdot, 1,ce)
+    call thermal_conduction_getSource(Tdot,1,ce)
     scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + Tdot) &
                             + thermal_conduction_getMassDensity (ce)* &
                               thermal_conduction_getSpecificHeat(ce)*(T_lastInc(i,j,k)  - &
@@ -310,7 +313,7 @@ subroutine updateReference
   mu_ref = 0.0_pReal
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
     ce = ce + 1
-    K_ref  = K_ref  + thermal_conduction_getConductivity(1,ce)
+    K_ref  = K_ref  + thermal_conduction_getConductivity(ce)
     mu_ref = mu_ref + thermal_conduction_getMassDensity(ce)* thermal_conduction_getSpecificHeat(ce)
   enddo; enddo; enddo
   K_ref = K_ref*wgt

src/homogenization.f90

@@ -35,19 +35,11 @@ module homogenization
     homogState, &
     damageState_h
 
-  integer, dimension(:), allocatable, public, protected :: &
+  integer(kind(THERMAL_isothermal_ID)),       dimension(:),   allocatable :: &
-    homogenization_typeInstance, &                                                                  !< instance of particular type of each homogenization
-    thermal_typeInstance, &                                                                         !< instance of particular type of each thermal transport
-    damage_typeInstance                                                                             !< instance of particular type of each nonlocal damage
-
-  real(pReal), dimension(:), allocatable, public, protected :: &
-    thermal_initialT
-
-  integer(kind(THERMAL_isothermal_ID)),       dimension(:),   allocatable, public, protected :: &
     thermal_type                                                                                    !< thermal transport model
-  integer(kind(DAMAGE_none_ID)),              dimension(:),   allocatable, public, protected :: &
+  integer(kind(DAMAGE_none_ID)),              dimension(:),   allocatable :: &
     damage_type                                                                                     !< nonlocal damage model
-  integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:),   allocatable, public, protected :: &
+  integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:),   allocatable :: &
     homogenization_type                                                                             !< type of each homogenization
 
   type, private :: tNumerics_damage
@@ -95,12 +87,11 @@ module homogenization
     module subroutine damage_init
     end subroutine damage_init
 
-    module subroutine mechanical_partition(subF,ip,el)
+    module subroutine mechanical_partition(subF,ce)
       real(pReal), intent(in), dimension(3,3) :: &
         subF
       integer,     intent(in) :: &
-        ip, &                                                                                       !< integration point
+        ce
-        el                                                                                          !< element number
     end subroutine mechanical_partition
 
     module subroutine thermal_partition(ce)
@@ -115,37 +106,31 @@ module homogenization
       integer, intent(in) :: ip,el
     end subroutine thermal_homogenize
 
-    module subroutine mechanical_homogenize(dt,ip,el)
+    module subroutine mechanical_homogenize(dt,ce)
      real(pReal), intent(in) :: dt
      integer, intent(in) :: &
-       ip, &                                                                                        !< integration point
+       ce                                                                                           !< cell
-       el                                                                                           !< element number
     end subroutine mechanical_homogenize
 
-    module subroutine mechanical_results(group_base,h)
+    module subroutine mechanical_results(group_base,ho)
       character(len=*), intent(in) :: group_base
-      integer, intent(in)          :: h
+      integer, intent(in)          :: ho
     end subroutine mechanical_results
 
-    module function mechanical_updateState(subdt,subF,ip,el) result(doneAndHappy)
+    module function mechanical_updateState(subdt,subF,ce) result(doneAndHappy)
       real(pReal), intent(in) :: &
-        subdt                                                                                           !< current time step
+        subdt                                                                                       !< current time step
       real(pReal), intent(in), dimension(3,3) :: &
         subF
       integer,     intent(in) :: &
-        ip, &                                                                                           !< integration point
+        ce                                                                                          !< cell
-        el                                                                                              !< element number
       logical, dimension(2) :: doneAndHappy
     end function mechanical_updateState
 
 
-    module function thermal_conduction_getConductivity(ip,el) result(K)
+    module function thermal_conduction_getConductivity(ce) result(K)
-
+      integer, intent(in) :: ce
-      integer, intent(in) :: &
-        ip, &                                                                                           !< integration point number
-        el                                                                                              !< element number
       real(pReal), dimension(3,3) :: K
-
     end function thermal_conduction_getConductivity
 
     module function thermal_conduction_getSpecificHeat(ce) result(c_P)
@@ -173,48 +158,37 @@ module homogenization
       real(pReal) :: T
     end function homogenization_thermal_T
 
-    module subroutine thermal_conduction_getSource(Tdot, ip,el)
+    module subroutine thermal_conduction_getSource(Tdot, ip, el)
       integer, intent(in) :: &
-        ip, &                                                                                           !< integration point number
+        ip, &
-        el                                                                                              !< element number
+        el
       real(pReal), intent(out) :: Tdot
     end subroutine thermal_conduction_getSource
 
-   module function damage_nonlocal_getMobility(ip,el) result(M)
+    module function damage_nonlocal_getMobility(ce) result(M)
-    integer, intent(in) :: &
+      integer, intent(in) :: ce
-      ip, &                                                                                           !< integration point number
+      real(pReal) :: M
-      el                                                                                              !< element number
-    real(pReal) :: M
     end function damage_nonlocal_getMobility
 
-    module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el)
+    module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ce)
-
+      integer, intent(in) :: ce
-      integer, intent(in) :: &
-        ip, &                                                                                           !< integration point number
-        el                                                                                              !< element number
       real(pReal),   intent(in) :: &
         phi
       real(pReal) :: &
         phiDot, dPhiDot_dPhi
     end subroutine damage_nonlocal_getSourceAndItsTangent
 
-
+    module subroutine damage_nonlocal_putNonLocalDamage(phi,ce)
-    module subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el)
+      integer, intent(in) :: ce
-
-      integer, intent(in) :: &
-        ip, &                                                                                           !< integration point number
-        el                                                                                              !< element number
       real(pReal),   intent(in) :: &
         phi
-
     end subroutine damage_nonlocal_putNonLocalDamage
 
-    module subroutine damage_nonlocal_results(homog,group)
+    module subroutine damage_nonlocal_results(ho,group)
-
+      integer,          intent(in) :: ho
-      integer,          intent(in) :: homog
       character(len=*), intent(in) :: group
-
     end subroutine damage_nonlocal_results
+
   end interface
 
   public ::  &
@@ -257,21 +231,18 @@ subroutine homogenization_init()
 
   allocate(homogState      (size(material_name_homogenization)))
   allocate(damageState_h   (size(material_name_homogenization)))
-  call material_parseHomogenization
+  call material_parseHomogenization()
-
 
   num_homog        => config_numerics%get('homogenization',defaultVal=emptyDict)
   num_homogGeneric => num_homog%get('generic',defaultVal=emptyDict)
 
-  num%nMPstate          = num_homogGeneric%get_asInt  ('nMPstate',     defaultVal=10)
+  num%nMPstate  = num_homogGeneric%get_asInt('nMPstate',defaultVal=10)
-  if (num%nMPstate < 1)                   call IO_error(301,ext_msg='nMPstate')
+  if (num%nMPstate < 1) call IO_error(301,ext_msg='nMPstate')
-
 
   call mechanical_init(num_homog)
   call thermal_init()
   call damage_init()
-
+  call damage_nonlocal_init()
-  call damage_nonlocal_init
 
 
 end subroutine homogenization_init
@@ -318,7 +289,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
 
 
         if (.not. doneAndHappy(1)) then
-          call mechanical_partition(homogenization_F(1:3,1:3,ce),ip,el)
+          call mechanical_partition(homogenization_F(1:3,1:3,ce),ce)
           converged = .true.
           do co = 1, myNgrains
             converged = converged .and. crystallite_stress(dt,co,ip,el)
@@ -327,7 +298,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
           if (.not. converged) then
             doneAndHappy = [.true.,.false.]
           else
-            doneAndHappy = mechanical_updateState(dt,homogenization_F(1:3,1:3,ce),ip,el)
+            doneAndHappy = mechanical_updateState(dt,homogenization_F(1:3,1:3,ce),ce)
             converged = all(doneAndHappy)
           endif
         endif
@@ -341,7 +312,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
   enddo
   !$OMP END DO
 
-  if (.not. terminallyIll ) then
+  if (.not. terminallyIll) then
     !$OMP DO PRIVATE(ho,ph,ce)
     do el = FEsolving_execElem(1),FEsolving_execElem(2)
       if (terminallyIll) continue
@@ -355,21 +326,22 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
             if (.not. terminallyIll) &                                                           ! so first signals terminally ill...
               print*, ' Integration point ', ip,' at element ', el, ' terminally ill'
             terminallyIll = .true.                                                                  ! ...and kills all others
-         endif
+          endif
-         call thermal_homogenize(ip,el)
         enddo
+        call thermal_homogenize(ip,el)
       enddo
     enddo
     !$OMP END DO
 
-    !$OMP DO PRIVATE(ho)
+    !$OMP DO PRIVATE(ho,ce)
     elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2)
       ho = material_homogenizationAt(el)
       IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2)
+        ce = (el-1)*discretization_nIPs + ip
         do co = 1, homogenization_Nconstituents(ho)
           call crystallite_orientations(co,ip,el)
         enddo
-        call mechanical_homogenize(dt,ip,el)
+        call mechanical_homogenize(dt,ce)
       enddo IpLooping3
     enddo elementLooping3
     !$OMP END DO
@@ -527,26 +499,25 @@ end subroutine damage_nonlocal_init
 !--------------------------------------------------------------------------------------------------
 !> @brief returns homogenized non local damage diffusion tensor in reference configuration
 !--------------------------------------------------------------------------------------------------
-function damage_nonlocal_getDiffusion(ip,el)
+function damage_nonlocal_getDiffusion(ce)
 
-  integer, intent(in) :: &
+  integer, intent(in) :: ce
-    ip, &                                                                                           !< integration point number
-    el                                                                                              !< element number
   real(pReal), dimension(3,3) :: &
     damage_nonlocal_getDiffusion
   integer :: &
-    homog, &
+    ho, &
-    grain
+    co
 
-  homog  = material_homogenizationAt(el)
+  ho  = material_homogenizationAt2(ce)
   damage_nonlocal_getDiffusion = 0.0_pReal
-  do grain = 1, homogenization_Nconstituents(homog)
+
+  do co = 1, homogenization_Nconstituents(ho)
     damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + &
-      crystallite_push33ToRef(grain,ip,el,lattice_D(1:3,1:3,material_phaseAt(grain,el)))
+      crystallite_push33ToRef(co,ce,lattice_D(1:3,1:3,material_phaseAt2(co,ce)))
   enddo
 
   damage_nonlocal_getDiffusion = &
-    num_damage%charLength**2*damage_nonlocal_getDiffusion/real(homogenization_Nconstituents(homog),pReal)
+    num_damage%charLength**2*damage_nonlocal_getDiffusion/real(homogenization_Nconstituents(ho),pReal)
 
 end function damage_nonlocal_getDiffusion
 
@@ -568,13 +539,9 @@ subroutine material_parseHomogenization
 
   material_homogenization => config_material%get('homogenization')
 
-  allocate(homogenization_type(size(material_name_homogenization)),           source=HOMOGENIZATION_undefined_ID)
+  allocate(homogenization_type(size(material_name_homogenization)), source=HOMOGENIZATION_undefined_ID)
-  allocate(thermal_type(size(material_name_homogenization)),                  source=THERMAL_isothermal_ID)
+  allocate(thermal_type(size(material_name_homogenization)),        source=THERMAL_isothermal_ID)
-  allocate(damage_type (size(material_name_homogenization)),                  source=DAMAGE_none_ID)
+  allocate(damage_type (size(material_name_homogenization)),        source=DAMAGE_none_ID)
-  allocate(homogenization_typeInstance(size(material_name_homogenization)),   source=0)
-  allocate(thermal_typeInstance(size(material_name_homogenization)),          source=0)
-  allocate(damage_typeInstance(size(material_name_homogenization)),           source=0)
-  allocate(thermal_initialT(size(material_name_homogenization)),              source=300.0_pReal)
 
   do h=1, size(material_name_homogenization)
     homog => material_homogenization%get(h)
@@ -590,28 +557,20 @@ subroutine material_parseHomogenization
         call IO_error(500,ext_msg=homogMech%get_asString('type'))
     end select
 
-    homogenization_typeInstance(h) = count(homogenization_type==homogenization_type(h))
+    if (homog%contains('thermal')) then
-
-    if(homog%contains('thermal')) then
       homogThermal => homog%get('thermal')
-        thermal_initialT(h) =  homogThermal%get_asFloat('T_0',defaultVal=300.0_pReal)
-
         select case (homogThermal%get_asString('type'))
-          case('isothermal')
+          case('pass')
-            thermal_type(h) = THERMAL_isothermal_ID
-          case('conduction')
             thermal_type(h) = THERMAL_conduction_ID
           case default
             call IO_error(500,ext_msg=homogThermal%get_asString('type'))
         end select
     endif
 
-    if(homog%contains('damage')) then
+    if (homog%contains('damage')) then
       homogDamage => homog%get('damage')
         select case (homogDamage%get_asString('type'))
-          case('none')
+          case('pass')
-            damage_type(h) = DAMAGE_none_ID
-          case('nonlocal')
             damage_type(h) = DAMAGE_nonlocal_ID
           case default
             call IO_error(500,ext_msg=homogDamage%get_asString('type'))
@@ -619,12 +578,6 @@ subroutine material_parseHomogenization
     endif
   enddo
 
-  do h=1, size(material_name_homogenization)
-    homogenization_typeInstance(h)  = count(homogenization_type(1:h) == homogenization_type(h))
-    thermal_typeInstance(h)         = count(thermal_type       (1:h) == thermal_type       (h))
-    damage_typeInstance(h)          = count(damage_type        (1:h) == damage_type        (h))
-  enddo
-
 end subroutine material_parseHomogenization
 
 

src/homogenization_damage.f90

@@ -85,22 +85,20 @@ end subroutine damage_partition
 !--------------------------------------------------------------------------------------------------
 !> @brief Returns homogenized nonlocal damage mobility
 !--------------------------------------------------------------------------------------------------
-module function damage_nonlocal_getMobility(ip,el) result(M)
+module function damage_nonlocal_getMobility(ce) result(M)
 
-  integer, intent(in) :: &
+  integer, intent(in) :: ce
-    ip, &                                                                                           !< integration point number
-    el                                                                                              !< element number
   integer :: &
     co
   real(pReal) :: M
 
   M = 0.0_pReal
 
-  do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
+  do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce))
-    M = M + lattice_M(material_phaseAt(co,el))
+    M = M + lattice_M(material_phaseAt2(co,ce))
   enddo
 
-  M = M/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
+  M = M/real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal)
 
 end function damage_nonlocal_getMobility
 
@@ -108,11 +106,9 @@ end function damage_nonlocal_getMobility
 !--------------------------------------------------------------------------------------------------
 !> @brief  calculates homogenized damage driving forces
 !--------------------------------------------------------------------------------------------------
-module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el)
+module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ce)
 
-  integer, intent(in) :: &
+  integer, intent(in) :: ce
-    ip, &                                                                                           !< integration point number
-    el                                                                                              !< element number
   real(pReal),   intent(in) :: &
     phi
   real(pReal) :: &
@@ -121,9 +117,9 @@ module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, p
   phiDot = 0.0_pReal
   dPhiDot_dPhi = 0.0_pReal
 
-  call phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el)
+  call phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce)
-  phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
+  phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal)
-  dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
+  dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal)
 
 end subroutine damage_nonlocal_getSourceAndItsTangent
 
@@ -131,20 +127,18 @@ end subroutine damage_nonlocal_getSourceAndItsTangent
 !--------------------------------------------------------------------------------------------------
 !> @brief updated nonlocal damage field with solution from damage phase field PDE
 !--------------------------------------------------------------------------------------------------
-module subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el)
+module subroutine damage_nonlocal_putNonLocalDamage(phi,ce)
 
-  integer, intent(in) :: &
+  integer, intent(in) :: ce
-    ip, &                                                                                           !< integration point number
-    el                                                                                              !< element number
   real(pReal),   intent(in) :: &
     phi
   integer :: &
-    homog, &
+    ho, &
-    offset
+    me
 
-  homog  = material_homogenizationAt(el)
+  ho = material_homogenizationAt2(ce)
-  offset = material_homogenizationMemberAt(ip,el)
+  me = material_homogenizationMemberAt2(ce)
-  damagestate_h(homog)%state(1,offset) = phi
+  damagestate_h(ho)%state(1,me) = phi
 
 end subroutine damage_nonlocal_putNonLocalDamage
 
@@ -152,18 +146,18 @@ end subroutine damage_nonlocal_putNonLocalDamage
 !--------------------------------------------------------------------------------------------------
 !> @brief writes results to HDF5 output file
 !--------------------------------------------------------------------------------------------------
-module subroutine damage_nonlocal_results(homog,group)
+module subroutine damage_nonlocal_results(ho,group)
 
-  integer,          intent(in) :: homog
+  integer,          intent(in) :: ho
   character(len=*), intent(in) :: group
 
   integer :: o
 
-  associate(prm => param(damage_typeInstance(homog)))
+  associate(prm => param(ho))
   outputsLoop: do o = 1,size(prm%output)
     select case(prm%output(o))
       case ('phi')
-        call results_writeDataset(group,damagestate_h(homog)%state(1,:),prm%output(o),&
+        call results_writeDataset(group,damagestate_h(ho)%state(1,:),prm%output(o),&
                                   'damage indicator','-')
     end select
   enddo outputsLoop

src/homogenization_mechanical.f90

@@ -24,35 +24,34 @@ submodule(homogenization) mechanical
       real(pReal),   dimension (3,3),   intent(in)  :: avgF                                         !< average deformation gradient at material point
     end subroutine mechanical_isostrain_partitionDeformation
 
-    module subroutine mechanical_RGC_partitionDeformation(F,avgF,instance,of)
+    module subroutine mechanical_RGC_partitionDeformation(F,avgF,ce)
       real(pReal),   dimension (:,:,:), intent(out) :: F                                            !< partitioned deformation gradient
       real(pReal),   dimension (3,3),   intent(in)  :: avgF                                         !< average deformation gradient at material point
       integer,                          intent(in)  :: &
-        instance, &
+        ce
-        of
     end subroutine mechanical_RGC_partitionDeformation
 
 
-    module subroutine mechanical_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance)
+    module subroutine mechanical_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho)
       real(pReal),   dimension (3,3),       intent(out) :: avgP                                     !< average stress at material point
       real(pReal),   dimension (3,3,3,3),   intent(out) :: dAvgPdAvgF                               !< average stiffness at material point
 
       real(pReal),   dimension (:,:,:),     intent(in)  :: P                                        !< partitioned stresses
       real(pReal),   dimension (:,:,:,:,:), intent(in)  :: dPdF                                     !< partitioned stiffnesses
-      integer,                              intent(in)  :: instance
+      integer,                              intent(in)  :: ho
     end subroutine mechanical_isostrain_averageStressAndItsTangent
 
-    module subroutine mechanical_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance)
+    module subroutine mechanical_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho)
       real(pReal),   dimension (3,3),       intent(out) :: avgP                                     !< average stress at material point
       real(pReal),   dimension (3,3,3,3),   intent(out) :: dAvgPdAvgF                               !< average stiffness at material point
 
       real(pReal),   dimension (:,:,:),     intent(in)  :: P                                        !< partitioned stresses
       real(pReal),   dimension (:,:,:,:,:), intent(in)  :: dPdF                                     !< partitioned stiffnesses
-      integer,                              intent(in)  :: instance
+      integer,                              intent(in)  :: ho
     end subroutine mechanical_RGC_averageStressAndItsTangent
 
 
-    module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAndHappy)
+    module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
       logical, dimension(2) :: doneAndHappy
       real(pReal), dimension(:,:,:),     intent(in)    :: &
         P,&                                                                                         !< partitioned stresses
@@ -61,13 +60,12 @@ submodule(homogenization) mechanical
       real(pReal), dimension(3,3),       intent(in) :: avgF                                         !< average F
       real(pReal),                       intent(in) :: dt                                           !< time increment
       integer,                           intent(in) :: &
-        ip, &                                                                                       !< integration point number
+        ce                                                                                          !< cell
-        el                                                                                          !< element number
     end function mechanical_RGC_updateState
 
 
-    module subroutine mechanical_RGC_results(instance,group)
+    module subroutine mechanical_RGC_results(ho,group)
-      integer,          intent(in) :: instance                                                      !< homogenization instance
+      integer,          intent(in) :: ho                                                            !< homogenization type
       character(len=*), intent(in) :: group                                                         !< group name in HDF5 file
     end subroutine mechanical_RGC_results
 
@@ -104,19 +102,18 @@ end subroutine mechanical_init
 !--------------------------------------------------------------------------------------------------
 !> @brief Partition F onto the individual constituents.
 !--------------------------------------------------------------------------------------------------
-module subroutine mechanical_partition(subF,ip,el)
+module subroutine mechanical_partition(subF,ce)
 
   real(pReal), intent(in), dimension(3,3) :: &
     subF
   integer,     intent(in) :: &
-    ip, &                                                                                           !< integration point
+    ce
-    el                                                                                              !< element number
 
   integer :: co
-  real(pReal), dimension (3,3,homogenization_Nconstituents(material_homogenizationAt(el))) :: Fs
+  real(pReal), dimension (3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) :: Fs
 
 
-  chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
+  chosenHomogenization: select case(homogenization_type(material_homogenizationAt2(ce)))
 
     case (HOMOGENIZATION_NONE_ID) chosenHomogenization
       Fs(1:3,1:3,1) = subF
@@ -125,12 +122,12 @@ module subroutine mechanical_partition(subF,ip,el)
       call mechanical_isostrain_partitionDeformation(Fs,subF)
 
     case (HOMOGENIZATION_RGC_ID) chosenHomogenization
-      call mechanical_RGC_partitionDeformation(Fs,subF,ip,el)
+      call mechanical_RGC_partitionDeformation(Fs,subF,ce)
 
   end select chosenHomogenization
 
-  do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
+  do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce))
-    call phase_mechanical_setF(Fs(1:3,1:3,co),co,ip,el)
+    call phase_mechanical_setF(Fs(1:3,1:3,co),co,ce)
   enddo
 
 
@@ -140,46 +137,43 @@ end subroutine mechanical_partition
 !--------------------------------------------------------------------------------------------------
 !> @brief Average P and dPdF from the individual constituents.
 !--------------------------------------------------------------------------------------------------
-module subroutine mechanical_homogenize(dt,ip,el)
+module subroutine mechanical_homogenize(dt,ce)
 
   real(pReal), intent(in) :: dt
-  integer, intent(in) :: &
+  integer, intent(in) :: ce
-       ip, &                                                                                        !< integration point
-       el                                                                                           !< element number
 
-  integer :: co,ce
+  integer :: co
-  real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
+  real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt2(ce)))
-  real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
+  real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationAt2(ce)))
 
 
-  ce = (el-1)* discretization_nIPs + ip
+  chosenHomogenization: select case(homogenization_type(material_homogenizationAt2(ce)))
-  chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
 
     case (HOMOGENIZATION_NONE_ID) chosenHomogenization
-        homogenization_P(1:3,1:3,ce)            = phase_mechanical_getP(1,ip,el)
+        homogenization_P(1:3,1:3,ce)            = phase_mechanical_getP(1,ce)
-        homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(dt,1,ip,el)
+        homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(dt,1,ce)
 
     case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
-      do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
+      do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce))
-        dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ip,el)
+        dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ce)
-        Ps(:,:,co) = phase_mechanical_getP(co,ip,el)
+        Ps(:,:,co) = phase_mechanical_getP(co,ce)
       enddo
       call mechanical_isostrain_averageStressAndItsTangent(&
         homogenization_P(1:3,1:3,ce), &
         homogenization_dPdF(1:3,1:3,1:3,1:3,ce),&
         Ps,dPdFs, &
-        homogenization_typeInstance(material_homogenizationAt(el)))
+        material_homogenizationAt2(ce))
 
     case (HOMOGENIZATION_RGC_ID) chosenHomogenization
-      do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
+      do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce))
-        dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ip,el)
+        dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ce)
-        Ps(:,:,co) = phase_mechanical_getP(co,ip,el)
+        Ps(:,:,co) = phase_mechanical_getP(co,ce)
       enddo
       call mechanical_RGC_averageStressAndItsTangent(&
         homogenization_P(1:3,1:3,ce), &
         homogenization_dPdF(1:3,1:3,1:3,1:3,ce),&
         Ps,dPdFs, &
-        homogenization_typeInstance(material_homogenizationAt(el)))
+        material_homogenizationAt2(ce))
 
   end select chosenHomogenization
 
@@ -190,30 +184,29 @@ end subroutine mechanical_homogenize
 !> @brief update the internal state of the homogenization scheme and tell whether "done" and
 !> "happy" with result
 !--------------------------------------------------------------------------------------------------
-module function mechanical_updateState(subdt,subF,ip,el) result(doneAndHappy)
+module function mechanical_updateState(subdt,subF,ce) result(doneAndHappy)
 
   real(pReal), intent(in) :: &
     subdt                                                                                           !< current time step
   real(pReal), intent(in), dimension(3,3) :: &
     subF
   integer,     intent(in) :: &
-    ip, &                                                                                           !< integration point
+    ce
-    el                                                                                              !< element number
   logical, dimension(2) :: doneAndHappy
 
   integer :: co
-  real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
+  real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt2(ce)))
-  real(pReal) :: Fs(3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
+  real(pReal) :: Fs(3,3,homogenization_Nconstituents(material_homogenizationAt2(ce)))
-  real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
+  real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationAt2(ce)))
 
 
-  if (homogenization_type(material_homogenizationAt(el)) == HOMOGENIZATION_RGC_ID) then
+  if (homogenization_type(material_homogenizationAt2(ce)) == HOMOGENIZATION_RGC_ID) then
-      do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
+      do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce))
-        dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(subdt,co,ip,el)
+        dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(subdt,co,ce)
-        Fs(:,:,co)        = phase_mechanical_getF(co,ip,el)
+        Fs(:,:,co)        = phase_mechanical_getF(co,ce)
-        Ps(:,:,co)        = phase_mechanical_getP(co,ip,el)
+        Ps(:,:,co)        = phase_mechanical_getP(co,ce)
       enddo
-      doneAndHappy = mechanical_RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ip,el)
+      doneAndHappy = mechanical_RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ce)
   else
     doneAndHappy = .true.
   endif
@@ -224,20 +217,20 @@ end function mechanical_updateState
 !--------------------------------------------------------------------------------------------------
 !> @brief Write results to file.
 !--------------------------------------------------------------------------------------------------
-module subroutine mechanical_results(group_base,h)
+module subroutine mechanical_results(group_base,ho)
 
   character(len=*), intent(in) :: group_base
-  integer, intent(in)          :: h
+  integer, intent(in)          :: ho
 
   character(len=:), allocatable :: group
 
   group = trim(group_base)//'/mech'
   call results_closeGroup(results_addGroup(group))
 
-  select case(homogenization_type(h))
+  select case(homogenization_type(ho))
 
     case(HOMOGENIZATION_rgc_ID)
-      call mechanical_RGC_results(homogenization_typeInstance(h),group)
+      call mechanical_RGC_results(ho,group)
 
   end select
 

src/homogenization_mechanical_RGC.f90

@@ -77,8 +77,7 @@ module subroutine mechanical_RGC_init(num_homogMech)
     num_homogMech                                                                                   !< pointer to mechanical homogenization numerics data
 
   integer :: &
-    Ninstances, &
+    ho, &
-    h, &
     Nmaterialpoints, &
     sizeState, nIntFaceTot
 
@@ -90,8 +89,7 @@ module subroutine mechanical_RGC_init(num_homogMech)
 
   print'(/,a)', ' <<<+-  homogenization:mechanical:RGC init  -+>>>'
 
-  Ninstances = count(homogenization_type == HOMOGENIZATION_RGC_ID)
+  print'(a,i2)', ' # instances: ',count(homogenization_type == HOMOGENIZATION_RGC_ID); flush(IO_STDOUT)
-  print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
 
   print*, 'Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009'
   print*, 'https://doi.org/10.1007/s12289-009-0619-1'//IO_EOL
@@ -101,10 +99,11 @@ module subroutine mechanical_RGC_init(num_homogMech)
 
 
 
-  allocate(param(Ninstances))
+  material_homogenization => config_material%get('homogenization')
-  allocate(state(Ninstances))
+  allocate(param(material_homogenization%length))
-  allocate(state0(Ninstances))
+  allocate(state(material_homogenization%length))
-  allocate(dependentState(Ninstances))
+  allocate(state0(material_homogenization%length))
+  allocate(dependentState(material_homogenization%length))
 
   num_RGC => num_homogMech%get('RGC',defaultVal=emptyDict)
 
@@ -137,15 +136,14 @@ module subroutine mechanical_RGC_init(num_homogMech)
   if (num%volDiscrPow <= 0.0_pReal)  call IO_error(301,ext_msg='volDiscrPw_RGC')
 
 
-  material_homogenization => config_material%get('homogenization')
+  do ho = 1, size(homogenization_type)
-  do h = 1, size(homogenization_type)
+    if (homogenization_type(ho) /= HOMOGENIZATION_RGC_ID) cycle
-    if (homogenization_type(h) /= HOMOGENIZATION_RGC_ID) cycle
+    homog => material_homogenization%get(ho)
-    homog => material_homogenization%get(h)
     homogMech => homog%get('mechanics')
-    associate(prm => param(homogenization_typeInstance(h)), &
+    associate(prm => param(ho), &
-              stt => state(homogenization_typeInstance(h)), &
+              stt => state(ho), &
-              st0 => state0(homogenization_typeInstance(h)), &
+              st0 => state0(ho), &
-              dst => dependentState(homogenization_typeInstance(h)))
+              dst => dependentState(ho))
 
 #if defined (__GFORTRAN__)
     prm%output = output_asStrings(homogMech)
@@ -154,7 +152,7 @@ module subroutine mechanical_RGC_init(num_homogMech)
 #endif
 
     prm%N_constituents = homogMech%get_asInts('cluster_size',requiredSize=3)
-    if (homogenization_Nconstituents(h) /= product(prm%N_constituents)) &
+    if (homogenization_Nconstituents(ho) /= product(prm%N_constituents)) &
       call IO_error(211,ext_msg='N_constituents (mechanical_RGC)')
 
     prm%xi_alpha = homogMech%get_asFloat('xi_alpha')
@@ -163,18 +161,18 @@ module subroutine mechanical_RGC_init(num_homogMech)
     prm%D_alpha  = homogMech%get_asFloats('D_alpha', requiredSize=3)
     prm%a_g      = homogMech%get_asFloats('a_g',     requiredSize=3)
 
-    Nmaterialpoints = count(material_homogenizationAt == h)
+    Nmaterialpoints = count(material_homogenizationAt == ho)
     nIntFaceTot = 3*(  (prm%N_constituents(1)-1)*prm%N_constituents(2)*prm%N_constituents(3) &
                       + prm%N_constituents(1)*(prm%N_constituents(2)-1)*prm%N_constituents(3) &
                       + prm%N_constituents(1)*prm%N_constituents(2)*(prm%N_constituents(3)-1))
     sizeState = nIntFaceTot
 
-    homogState(h)%sizeState = sizeState
+    homogState(ho)%sizeState = sizeState
-    allocate(homogState(h)%state0   (sizeState,Nmaterialpoints), source=0.0_pReal)
+    allocate(homogState(ho)%state0   (sizeState,Nmaterialpoints), source=0.0_pReal)
-    allocate(homogState(h)%state    (sizeState,Nmaterialpoints), source=0.0_pReal)
+    allocate(homogState(ho)%state    (sizeState,Nmaterialpoints), source=0.0_pReal)
 
-    stt%relaxationVector   => homogState(h)%state(1:nIntFaceTot,:)
+    stt%relaxationVector   => homogState(ho)%state(1:nIntFaceTot,:)
-    st0%relaxationVector   => homogState(h)%state0(1:nIntFaceTot,:)
+    st0%relaxationVector   => homogState(ho)%state0(1:nIntFaceTot,:)
 
     allocate(dst%volumeDiscrepancy(   Nmaterialpoints), source=0.0_pReal)
     allocate(dst%relaxationRate_avg(  Nmaterialpoints), source=0.0_pReal)
@@ -183,7 +181,7 @@ module subroutine mechanical_RGC_init(num_homogMech)
 
 !--------------------------------------------------------------------------------------------------
 ! assigning cluster orientations
-    dependentState(homogenization_typeInstance(h))%orientation = spread(eu2om(prm%a_g*inRad),3,Nmaterialpoints)
+    dependentState(ho)%orientation = spread(eu2om(prm%a_g*inRad),3,Nmaterialpoints)
     !dst%orientation = spread(eu2om(prm%a_g*inRad),3,Nmaterialpoints) ifort version 18.0.1 crashes (for whatever reason)
 
     end associate
@@ -196,22 +194,23 @@ end subroutine mechanical_RGC_init
 !--------------------------------------------------------------------------------------------------
 !> @brief partitions the deformation gradient onto the constituents
 !--------------------------------------------------------------------------------------------------
-module subroutine mechanical_RGC_partitionDeformation(F,avgF,instance,of)
+module subroutine mechanical_RGC_partitionDeformation(F,avgF,ce)
 
   real(pReal),   dimension (:,:,:), intent(out) :: F                                                !< partitioned F  per grain
 
   real(pReal),   dimension (3,3),   intent(in)  :: avgF                                             !< averaged F
   integer,                          intent(in)  :: &
-    instance, &
+    ce
-    of
 
   real(pReal), dimension(3) :: aVect,nVect
   integer,     dimension(4) :: intFace
   integer,     dimension(3) :: iGrain3
-  integer ::  iGrain,iFace,i,j
+  integer ::  iGrain,iFace,i,j,ho,me
-
-  associate(prm => param(instance))
 
+  associate(prm => param(material_homogenizationAt2(ce)))
+ 
+  ho = material_homogenizationAt2(ce)
+  me = material_homogenizationMemberAt2(ce)
 !--------------------------------------------------------------------------------------------------
 ! compute the deformation gradient of individual grains due to relaxations
   F = 0.0_pReal
@@ -219,8 +218,8 @@ module subroutine mechanical_RGC_partitionDeformation(F,avgF,instance,of)
     iGrain3 = grain1to3(iGrain,prm%N_constituents)
     do iFace = 1,6
       intFace = getInterface(iFace,iGrain3)                                                         ! identifying 6 interfaces of each grain
-      aVect = relaxationVector(intFace,instance,of)                                                 ! get the relaxation vectors for each interface from global relaxation vector array
+      aVect = relaxationVector(intFace,ho,me)                                                       ! get the relaxation vectors for each interface from global relaxation vector array
-      nVect = interfaceNormal(intFace,instance,of)
+      nVect = interfaceNormal(intFace,ho,me)
       forall (i=1:3,j=1:3) &
         F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j)                                           ! calculating deformation relaxations due to interface relaxation
     enddo
@@ -236,7 +235,7 @@ end subroutine mechanical_RGC_partitionDeformation
 !> @brief update the internal state of the homogenization scheme and tell whether "done" and
 ! "happy" with result
 !--------------------------------------------------------------------------------------------------
-module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAndHappy)
+module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
       logical, dimension(2) :: doneAndHappy
       real(pReal), dimension(:,:,:),     intent(in)    :: &
         P,&                                                                                         !< partitioned stresses
@@ -245,12 +244,11 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn
       real(pReal), dimension(3,3),       intent(in) :: avgF                                         !< average F
       real(pReal),                       intent(in) :: dt                                           !< time increment
       integer,                           intent(in) :: &
-        ip, &                                                                                       !< integration point number
+        ce                                                                                          !< cell
-        el                                                                                          !< element number
 
   integer, dimension(4) :: intFaceN,intFaceP,faceID
   integer, dimension(3) :: nGDim,iGr3N,iGr3P
-  integer :: instance,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ipert,iGrain,nGrain, of
+  integer :: ho,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ipert,iGrain,nGrain, me
   real(pReal), dimension(3,3,size(P,3)) :: R,pF,pR,D,pD
   real(pReal), dimension(3,size(P,3))   :: NN,devNull
   real(pReal), dimension(3)             :: normP,normN,mornP,mornN
@@ -264,10 +262,10 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn
     return
   endif zeroTimeStep
 
-  instance  = homogenization_typeInstance(material_homogenizationAt(el))
+  ho  = material_homogenizationAt2(ce)
-  of = material_homogenizationMemberAt(ip,el)
 
-  associate(stt => state(instance), st0 => state0(instance), dst => dependentState(instance), prm => param(instance))
+  me = material_homogenizationMemberAt2(ce)
+  associate(stt => state(ho), st0 => state0(ho), dst => dependentState(ho), prm => param(ho))
 
 !--------------------------------------------------------------------------------------------------
 ! get the dimension of the cluster (grains and interfaces)
@@ -281,36 +279,36 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn
 ! allocate the size of the global relaxation arrays/jacobian matrices depending on the size of the cluster
   allocate(resid(3*nIntFaceTot), source=0.0_pReal)
   allocate(tract(nIntFaceTot,3), source=0.0_pReal)
-  relax  = stt%relaxationVector(:,of)
+  relax  = stt%relaxationVector(:,me)
-  drelax = stt%relaxationVector(:,of) - st0%relaxationVector(:,of)
+  drelax = stt%relaxationVector(:,me) - st0%relaxationVector(:,me)
 
 !--------------------------------------------------------------------------------------------------
 ! computing interface mismatch and stress penalty tensor for all interfaces of all grains
-  call stressPenalty(R,NN,avgF,F,ip,el,instance,of)
+  call stressPenalty(R,NN,avgF,F,ho,me)
 
 !--------------------------------------------------------------------------------------------------
 ! calculating volume discrepancy and stress penalty related to overall volume discrepancy
-  call volumePenalty(D,dst%volumeDiscrepancy(of),avgF,F,nGrain,instance,of)
+  call volumePenalty(D,dst%volumeDiscrepancy(me),avgF,F,nGrain)
 
 !------------------------------------------------------------------------------------------------
 ! computing the residual stress from the balance of traction at all (interior) interfaces
   do iNum = 1,nIntFaceTot
-    faceID = interface1to4(iNum,param(instance)%N_constituents)                                     ! identifying the interface ID in local coordinate system (4-dimensional index)
+    faceID = interface1to4(iNum,param(ho)%N_constituents)                                           ! identifying the interface ID in local coordinate system (4-dimensional index)
 
 !--------------------------------------------------------------------------------------------------
 ! identify the left/bottom/back grain (-|N)
     iGr3N = faceID(2:4)                                                                             ! identifying the grain ID in local coordinate system (3-dimensional index)
-    iGrN = grain3to1(iGr3N,param(instance)%N_constituents)                                          ! translate the local grain ID into global coordinate system (1-dimensional index)
+    iGrN = grain3to1(iGr3N,param(ho)%N_constituents)                                                ! translate the local grain ID into global coordinate system (1-dimensional index)
     intFaceN = getInterface(2*faceID(1),iGr3N)
-    normN = interfaceNormal(intFaceN,instance,of)
+    normN = interfaceNormal(intFaceN,ho,me)
 
 !--------------------------------------------------------------------------------------------------
 ! identify the right/up/front grain (+|P)
     iGr3P = iGr3N
     iGr3P(faceID(1)) = iGr3N(faceID(1))+1                                                           ! identifying the grain ID in local coordinate system (3-dimensional index)
-    iGrP = grain3to1(iGr3P,param(instance)%N_constituents)                                          ! translate the local grain ID into global coordinate system (1-dimensional index)
+    iGrP = grain3to1(iGr3P,param(ho)%N_constituents)                                                ! translate the local grain ID into global coordinate system (1-dimensional index)
     intFaceP = getInterface(2*faceID(1)-1,iGr3P)
-    normP = interfaceNormal(intFaceP,instance,of)
+    normP = interfaceNormal(intFaceP,ho,me)
 
 !--------------------------------------------------------------------------------------------------
 ! compute the residual of traction at the interface (in local system, 4-dimensional index)
@@ -338,9 +336,9 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn
   if (residMax < num%rtol*stresMax .or. residMax < num%atol) then
     doneAndHappy = .true.
 
-    dst%mismatch(1:3,of)       = sum(NN,2)/real(nGrain,pReal)
+    dst%mismatch(1:3,me)       = sum(NN,2)/real(nGrain,pReal)
-    dst%relaxationRate_avg(of) = sum(abs(drelax))/dt/real(3*nIntFaceTot,pReal)
+    dst%relaxationRate_avg(me) = sum(abs(drelax))/dt/real(3*nIntFaceTot,pReal)
-    dst%relaxationRate_max(of) = maxval(abs(drelax))/dt
+    dst%relaxationRate_max(me) = maxval(abs(drelax))/dt
 
     return
 
@@ -359,18 +357,18 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn
 ! ... of the constitutive stress tangent, assembled from dPdF or material constitutive model "smatrix"
   allocate(smatrix(3*nIntFaceTot,3*nIntFaceTot), source=0.0_pReal)
   do iNum = 1,nIntFaceTot
-    faceID = interface1to4(iNum,param(instance)%N_constituents)                                     ! assembling of local dPdF into global Jacobian matrix
+    faceID = interface1to4(iNum,param(ho)%N_constituents)                                           ! assembling of local dPdF into global Jacobian matrix
 
 !--------------------------------------------------------------------------------------------------
 ! identify the left/bottom/back grain (-|N)
     iGr3N = faceID(2:4)                                                                             ! identifying the grain ID in local coordinate sytem
-    iGrN = grain3to1(iGr3N,param(instance)%N_constituents)                                          ! translate into global grain ID
+    iGrN = grain3to1(iGr3N,param(ho)%N_constituents)                                                ! translate into global grain ID
     intFaceN = getInterface(2*faceID(1),iGr3N)                                                      ! identifying the connecting interface in local coordinate system
-    normN = interfaceNormal(intFaceN,instance,of)
+    normN = interfaceNormal(intFaceN,ho,me)
     do iFace = 1,6
       intFaceN = getInterface(iFace,iGr3N)                                                          ! identifying all interfaces that influence relaxation of the above interface
-      mornN = interfaceNormal(intFaceN,instance,of)
+      mornN = interfaceNormal(intFaceN,ho,me)
-      iMun = interface4to1(intFaceN,param(instance)%N_constituents)                                 ! translate the interfaces ID into local 4-dimensional index
+      iMun = interface4to1(intFaceN,param(ho)%N_constituents)                                       ! translate the interfaces ID into local 4-dimensional index
       if (iMun > 0) then                                                                            ! get the corresponding tangent
         do i=1,3; do j=1,3; do k=1,3; do l=1,3
           smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) &
@@ -385,13 +383,13 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn
 ! identify the right/up/front grain (+|P)
     iGr3P = iGr3N
     iGr3P(faceID(1)) = iGr3N(faceID(1))+1                                                           ! identifying the grain ID in local coordinate sytem
-    iGrP = grain3to1(iGr3P,param(instance)%N_constituents)                                          ! translate into global grain ID
+    iGrP = grain3to1(iGr3P,param(ho)%N_constituents)                                                ! translate into global grain ID
     intFaceP = getInterface(2*faceID(1)-1,iGr3P)                                                    ! identifying the connecting interface in local coordinate system
-    normP = interfaceNormal(intFaceP,instance,of)
+    normP = interfaceNormal(intFaceP,ho,me)
     do iFace = 1,6
       intFaceP = getInterface(iFace,iGr3P)                                                          ! identifying all interfaces that influence relaxation of the above interface
-      mornP = interfaceNormal(intFaceP,instance,of)
+      mornP = interfaceNormal(intFaceP,ho,me)
-      iMun = interface4to1(intFaceP,param(instance)%N_constituents)                                 ! translate the interfaces ID into local 4-dimensional index
+      iMun = interface4to1(intFaceP,param(ho)%N_constituents)                                       ! translate the interfaces ID into local 4-dimensional index
       if (iMun > 0) then                                                                            ! get the corresponding tangent
         do i=1,3; do j=1,3; do k=1,3; do l=1,3
           smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) &
@@ -411,31 +409,31 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn
   do ipert = 1,3*nIntFaceTot
     p_relax = relax
     p_relax(ipert) = relax(ipert) + num%pPert                                                       ! perturb the relaxation vector
-    stt%relaxationVector(:,of) = p_relax
+    stt%relaxationVector(:,me) = p_relax
-    call grainDeformation(pF,avgF,instance,of)                                                      ! rain deformation from perturbed state
+    call grainDeformation(pF,avgF,ho,me)                                                            ! rain deformation from perturbed state
-    call stressPenalty(pR,DevNull,      avgF,pF,ip,el,instance,of)                                  ! stress penalty due to interface mismatch from perturbed state
+    call stressPenalty(pR,DevNull,      avgF,pF,ho,me)                                              ! stress penalty due to interface mismatch from perturbed state
-    call volumePenalty(pD,devNull(1,1), avgF,pF,nGrain,instance,of)                                 ! stress penalty due to volume discrepancy from perturbed state
+    call volumePenalty(pD,devNull(1,1), avgF,pF,nGrain)                                             ! stress penalty due to volume discrepancy from perturbed state
 
 !--------------------------------------------------------------------------------------------------
 ! computing the global stress residual array from the perturbed state
     p_resid = 0.0_pReal
     do iNum = 1,nIntFaceTot
-      faceID = interface1to4(iNum,param(instance)%N_constituents)                                   ! identifying the interface ID in local coordinate system (4-dimensional index)
+      faceID = interface1to4(iNum,param(ho)%N_constituents)                                         ! identifying the interface ID in local coordinate system (4-dimensional index)
 
 !--------------------------------------------------------------------------------------------------
 ! identify the left/bottom/back grain (-|N)
       iGr3N = faceID(2:4)                                                                           ! identify the grain ID in local coordinate system (3-dimensional index)
-      iGrN = grain3to1(iGr3N,param(instance)%N_constituents)                                        ! translate the local grain ID into global coordinate system (1-dimensional index)
+      iGrN = grain3to1(iGr3N,param(ho)%N_constituents)                                              ! translate the local grain ID into global coordinate system (1-dimensional index)
       intFaceN = getInterface(2*faceID(1),iGr3N)                                                    ! identify the interface ID of the grain
-      normN = interfaceNormal(intFaceN,instance,of)
+      normN = interfaceNormal(intFaceN,ho,me)
 
 !--------------------------------------------------------------------------------------------------
 ! identify the right/up/front grain (+|P)
       iGr3P = iGr3N
       iGr3P(faceID(1)) = iGr3N(faceID(1))+1                                                         ! identify the grain ID in local coordinate system (3-dimensional index)
-      iGrP = grain3to1(iGr3P,param(instance)%N_constituents)                                        ! translate the local grain ID into global coordinate system (1-dimensional index)
+      iGrP = grain3to1(iGr3P,param(ho)%N_constituents)                                              ! translate the local grain ID into global coordinate system (1-dimensional index)
       intFaceP = getInterface(2*faceID(1)-1,iGr3P)                                                  ! identify the interface ID of the grain
-      normP = interfaceNormal(intFaceP,instance,of)
+      normP = interfaceNormal(intFaceP,ho,me)
 
 !--------------------------------------------------------------------------------------------------
 ! compute the residual stress (contribution of mismatch and volume penalties) from perturbed state
@@ -475,11 +473,11 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn
   do i = 1,3*nIntFaceTot;do j = 1,3*nIntFaceTot
     drelax(i) = drelax(i) - jnverse(i,j)*resid(j)                                                   ! Calculate the correction for the state variable
   enddo; enddo
-  stt%relaxationVector(:,of) = relax + drelax                                                       ! Updateing the state variable for the next iteration
+  stt%relaxationVector(:,me) = relax + drelax                                                       ! Updateing the state variable for the next iteration
   if (any(abs(drelax) > num%maxdRelax)) then                                                        ! Forcing cutback when the incremental change of relaxation vector becomes too large
     doneAndHappy = [.true.,.false.]
     !$OMP CRITICAL (write2out)
-    print'(a,i3,a,i3,a)',' RGC_updateState: ip ',ip,' | el ',el,' enforces cutback'
+    print'(a,i3,a,i3,a)',' RGC_updateState: enforces cutback'
     print'(a,e15.8)',' due to large relaxation change = ',maxval(abs(drelax))
     flush(IO_STDOUT)
     !$OMP END CRITICAL (write2out)
@@ -491,14 +489,14 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn
   !------------------------------------------------------------------------------------------------
   !> @brief calculate stress-like penalty due to deformation mismatch
   !------------------------------------------------------------------------------------------------
-  subroutine stressPenalty(rPen,nMis,avgF,fDef,ip,el,instance,of)
+  subroutine stressPenalty(rPen,nMis,avgF,fDef,ho,me)
 
     real(pReal),   dimension (:,:,:), intent(out) :: rPen                                           !< stress-like penalty
     real(pReal),   dimension (:,:),   intent(out) :: nMis                                           !< total amount of mismatch
 
     real(pReal),   dimension (:,:,:), intent(in)  :: fDef                                           !< deformation gradients
     real(pReal),   dimension (3,3),   intent(in)  :: avgF                                           !< initial effective stretch tensor
-    integer,                          intent(in)  :: ip,el,instance,of
+    integer,                          intent(in)  :: ho, me
 
     integer, dimension (4)   :: intFace
     integer, dimension (3)   :: iGrain3,iGNghb3,nGDim
@@ -510,7 +508,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn
       nDefToler = 1.0e-10_pReal, &
       b = 2.5e-10_pReal                                                                             ! Length of Burgers vector
 
-    nGDim = param(instance)%N_constituents
+    nGDim = param(ho)%N_constituents
     rPen = 0.0_pReal
     nMis = 0.0_pReal
 
@@ -518,27 +516,26 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn
     ! get the correction factor the modulus of penalty stress representing the evolution of area of
     ! the interfaces due to deformations
 
-    surfCorr = surfaceCorrection(avgF,instance,of)
+    surfCorr = surfaceCorrection(avgF,ho,me)
-
-    associate(prm => param(instance))
 
+    associate(prm => param(ho))
 
    !-----------------------------------------------------------------------------------------------
    ! computing the mismatch and penalty stress tensor of all grains
    grainLoop: do iGrain = 1,product(prm%N_constituents)
-     muGrain = equivalentMu(iGrain,ip,el)
+     muGrain = equivalentMu(iGrain,ce)
      iGrain3 = grain1to3(iGrain,prm%N_constituents)                                                 ! get the grain ID in local 3-dimensional index (x,y,z)-position
 
      interfaceLoop: do iFace = 1,6
        intFace = getInterface(iFace,iGrain3)                                                        ! get the 4-dimensional index of the interface in local numbering system of the grain
-       nVect = interfaceNormal(intFace,instance,of)
+       nVect = interfaceNormal(intFace,ho,me)
        iGNghb3 = iGrain3                                                                            ! identify the neighboring grain across the interface
        iGNghb3(abs(intFace(1))) = iGNghb3(abs(intFace(1))) &
                                 + int(real(intFace(1),pReal)/real(abs(intFace(1)),pReal))
        where(iGNghb3 < 1)    iGNghb3 = nGDim
        where(iGNghb3 >nGDim) iGNghb3 = 1
        iGNghb  = grain3to1(iGNghb3,prm%N_constituents)                                              ! get the ID of the neighboring grain
-       muGNghb = equivalentMu(iGNghb,ip,el)
+       muGNghb = equivalentMu(iGNghb,ce)
        gDef = 0.5_pReal*(fDef(1:3,1:3,iGNghb) - fDef(1:3,1:3,iGrain))                               ! difference/jump in deformation gradeint across the neighbor
 
        !-------------------------------------------------------------------------------------------
@@ -577,7 +574,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn
   !------------------------------------------------------------------------------------------------
   !> @brief calculate stress-like penalty due to volume discrepancy
   !------------------------------------------------------------------------------------------------
-  subroutine volumePenalty(vPen,vDiscrep,fAvg,fDef,nGrain,instance,of)
+  subroutine volumePenalty(vPen,vDiscrep,fAvg,fDef,nGrain)
 
     real(pReal), dimension (:,:,:), intent(out) :: vPen                                             ! stress-like penalty due to volume
     real(pReal),                    intent(out) :: vDiscrep                                         ! total volume discrepancy
@@ -585,9 +582,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn
     real(pReal), dimension (:,:,:), intent(in)  :: fDef                                             ! deformation gradients
     real(pReal), dimension (3,3),   intent(in)  :: fAvg                                             ! overall deformation gradient
     integer,                        intent(in) :: &
-      Ngrain, &
+      Ngrain
-      instance, &
-      of
 
     real(pReal), dimension(size(vPen,3)) :: gVol
     integer :: i
@@ -617,14 +612,14 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn
   !> @brief compute the correction factor accouted for surface evolution (area change) due to
   ! deformation
   !--------------------------------------------------------------------------------------------------
-  function surfaceCorrection(avgF,instance,of)
+  function surfaceCorrection(avgF,ho,me)
 
     real(pReal), dimension(3)               :: surfaceCorrection
 
     real(pReal), dimension(3,3), intent(in) :: avgF                                                 !< average F
     integer,                     intent(in) :: &
-      instance, &
+      ho, &
-      of
+      me
     real(pReal), dimension(3,3)             :: invC
     real(pReal), dimension(3)               :: nVect
     real(pReal)  :: detF
@@ -635,7 +630,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn
 
     surfaceCorrection = 0.0_pReal
     do iBase = 1,3
-      nVect = interfaceNormal([iBase,1,1,1],instance,of)
+      nVect = interfaceNormal([iBase,1,1,1],ho,me)
       do i = 1,3; do j = 1,3
         surfaceCorrection(iBase) = surfaceCorrection(iBase) + invC(i,j)*nVect(i)*nVect(j)           ! compute the component of (the inverse of) the stretch in the direction of the normal
       enddo; enddo
@@ -648,17 +643,16 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn
   !-------------------------------------------------------------------------------------------------
   !> @brief compute the equivalent shear and bulk moduli from the elasticity tensor
   !-------------------------------------------------------------------------------------------------
-  real(pReal) function equivalentMu(grainID,ip,el)
+  real(pReal) function equivalentMu(grainID,ce)
 
     integer, intent(in)    :: &
       grainID,&
-      ip, &                                                                                         !< integration point number
+      ce
-      el                                                                                            !< element number
 
     real(pReal), dimension(6,6) :: C
 
 
-    C = phase_homogenizedC(material_phaseAt(grainID,el),material_phaseMemberAt(grainID,ip,el))
+    C = phase_homogenizedC(material_phaseAt2(grainID,ce),material_phaseMemberAt2(grainID,ce))
     equivalentMu = lattice_equivalent_mu(C,'voigt')
 
   end function equivalentMu
@@ -668,14 +662,14 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn
   !> @brief calculating the grain deformation gradient (the same with
   ! homogenization_RGC_partitionDeformation, but used only for perturbation scheme)
   !-------------------------------------------------------------------------------------------------
-  subroutine grainDeformation(F, avgF, instance, of)
+  subroutine grainDeformation(F, avgF, ho, me)
 
     real(pReal),   dimension(:,:,:), intent(out) :: F                                               !< partitioned F  per grain
 
     real(pReal),   dimension(:,:),   intent(in)  :: avgF                                            !< averaged F
     integer,                          intent(in)  :: &
-      instance, &
+      ho, &
-      of
+      me
 
     real(pReal),   dimension(3) :: aVect,nVect
     integer,       dimension(4) :: intFace
@@ -685,15 +679,15 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ip,el) result(doneAn
     !-----------------------------------------------------------------------------------------------
     ! compute the deformation gradient of individual grains due to relaxations
 
-    associate(prm => param(instance))
+    associate (prm => param(ho))
 
     F = 0.0_pReal
     do iGrain = 1,product(prm%N_constituents)
       iGrain3 = grain1to3(iGrain,prm%N_constituents)
       do iFace = 1,6
         intFace = getInterface(iFace,iGrain3)
-        aVect   = relaxationVector(intFace,instance,of)
+        aVect   = relaxationVector(intFace,ho,me)
-        nVect   = interfaceNormal(intFace,instance,of)
+        nVect   = interfaceNormal(intFace,ho,me)
         forall (i=1:3,j=1:3) &
           F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j)                                         ! effective relaxations
       enddo
@@ -710,17 +704,17 @@ end function mechanical_RGC_updateState
 !--------------------------------------------------------------------------------------------------
 !> @brief derive average stress and stiffness from constituent quantities
 !--------------------------------------------------------------------------------------------------
-module subroutine mechanical_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance)
+module subroutine mechanical_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho)
 
   real(pReal), dimension (3,3),        intent(out) :: avgP                                          !< average stress at material point
   real(pReal), dimension (3,3,3,3),    intent(out) :: dAvgPdAvgF                                    !< average stiffness at material point
 
   real(pReal), dimension (:,:,:),      intent(in)  :: P                                             !< partitioned stresses
   real(pReal), dimension (:,:,:,:,:),  intent(in)  :: dPdF                                          !< partitioned stiffnesses
-  integer,                             intent(in)  :: instance
+  integer,                             intent(in)  :: ho
 
-  avgP       = sum(P,3)   /real(product(param(instance)%N_constituents),pReal)
+  avgP       = sum(P,3)   /real(product(param(ho)%N_constituents),pReal)
-  dAvgPdAvgF = sum(dPdF,5)/real(product(param(instance)%N_constituents),pReal)
+  dAvgPdAvgF = sum(dPdF,5)/real(product(param(ho)%N_constituents),pReal)
 
 end subroutine mechanical_RGC_averageStressAndItsTangent
 
@@ -728,14 +722,14 @@ end subroutine mechanical_RGC_averageStressAndItsTangent
 !--------------------------------------------------------------------------------------------------
 !> @brief writes results to HDF5 output file
 !--------------------------------------------------------------------------------------------------
-module subroutine mechanical_RGC_results(instance,group)
+module subroutine mechanical_RGC_results(ho,group)
 
-  integer,          intent(in) :: instance
+  integer,          intent(in) :: ho
   character(len=*), intent(in) :: group
 
   integer :: o
 
-  associate(stt => state(instance), dst => dependentState(instance), prm => param(instance))
+  associate(stt => state(ho), dst => dependentState(ho), prm => param(ho))
   outputsLoop: do o = 1,size(prm%output)
     select case(trim(prm%output(o)))
       case('M')
@@ -760,11 +754,11 @@ end subroutine mechanical_RGC_results
 !--------------------------------------------------------------------------------------------------
 !> @brief collect relaxation vectors of an interface
 !--------------------------------------------------------------------------------------------------
-pure function relaxationVector(intFace,instance,of)
+pure function relaxationVector(intFace,ho,me)
 
   real(pReal), dimension (3)            :: relaxationVector
 
-  integer,                   intent(in) :: instance,of
+  integer,                   intent(in) :: ho,me
   integer,     dimension(4), intent(in) :: intFace                                                  !< set of interface ID in 4D array (normal and position)
 
   integer :: iNum
@@ -772,29 +766,35 @@ pure function relaxationVector(intFace,instance,of)
 !--------------------------------------------------------------------------------------------------
 ! collect the interface relaxation vector from the global state array
 
-  iNum = interface4to1(intFace,param(instance)%N_constituents)                                      ! identify the position of the interface in global state array
+  associate (prm => param(ho), &
+             stt => state(ho))
+
+  iNum = interface4to1(intFace,prm%N_constituents)                                                  ! identify the position of the interface in global state array
   if (iNum > 0) then
-    relaxationVector = state(instance)%relaxationVector((3*iNum-2):(3*iNum),of)
+    relaxationVector = stt%relaxationVector((3*iNum-2):(3*iNum),me)
   else
     relaxationVector = 0.0_pReal
   endif
 
+  end associate
+
 end function relaxationVector
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief identify the normal of an interface
 !--------------------------------------------------------------------------------------------------
-pure function interfaceNormal(intFace,instance,of)
+pure function interfaceNormal(intFace,ho,me)
 
   real(pReal), dimension(3)             :: interfaceNormal
 
   integer,     dimension(4), intent(in) :: intFace                                                  !< interface ID in 4D array (normal and position)
   integer,                   intent(in) :: &
-    instance, &
+    ho, &
-    of
+    me
 
   integer :: nPos
+  associate (dst => dependentState(ho))
 
 !--------------------------------------------------------------------------------------------------
 ! get the normal of the interface, identified from the value of intFace(1)
@@ -802,7 +802,9 @@ pure function interfaceNormal(intFace,instance,of)
   nPos = abs(intFace(1))                                                                            ! identify the position of the interface in global state array
   interfaceNormal(nPos) = real(intFace(1)/abs(intFace(1)),pReal)                                    ! get the normal vector w.r.t. cluster axis
 
-  interfaceNormal = matmul(dependentState(instance)%orientation(1:3,1:3,of),interfaceNormal)        ! map the normal vector into sample coordinate system (basis)
+  interfaceNormal = matmul(dst%orientation(1:3,1:3,me),interfaceNormal)                             ! map the normal vector into sample coordinate system (basis)
+  
+  end associate
 
 end function interfaceNormal
 

src/homogenization_mechanical_isostrain.f90

@@ -29,7 +29,6 @@ contains
 module subroutine mechanical_isostrain_init
 
   integer :: &
-    Ninstances, &
     h, &
     Nmaterialpoints
   class(tNode), pointer :: &
@@ -39,17 +38,16 @@ module subroutine mechanical_isostrain_init
 
   print'(/,a)', ' <<<+-  homogenization:mechanical:isostrain init  -+>>>'
 
-  Ninstances = count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)
+  print'(a,i2)', ' # instances: ',count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID); flush(IO_STDOUT)
-  print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
-
-  allocate(param(Ninstances))                                                                        ! one container of parameters per instance
 
   material_homogenization => config_material%get('homogenization')
+  allocate(param(material_homogenization%length))                                                   ! one container of parameters per homog
+
   do h = 1, size(homogenization_type)
     if (homogenization_type(h) /= HOMOGENIZATION_ISOSTRAIN_ID) cycle
     homog => material_homogenization%get(h)
     homogMech => homog%get('mechanics')
-    associate(prm => param(homogenization_typeInstance(h)))
+    associate(prm => param(h))
 
     prm%N_constituents = homogenization_Nconstituents(h)
     select case(homogMech%get_asString('mapping',defaultVal = 'sum'))
@@ -90,16 +88,16 @@ end subroutine mechanical_isostrain_partitionDeformation
 !--------------------------------------------------------------------------------------------------
 !> @brief derive average stress and stiffness from constituent quantities
 !--------------------------------------------------------------------------------------------------
-module subroutine mechanical_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance)
+module subroutine mechanical_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho)
 
   real(pReal),   dimension (3,3),       intent(out) :: avgP                                         !< average stress at material point
   real(pReal),   dimension (3,3,3,3),   intent(out) :: dAvgPdAvgF                                   !< average stiffness at material point
 
   real(pReal),   dimension (:,:,:),     intent(in)  :: P                                            !< partitioned stresses
   real(pReal),   dimension (:,:,:,:,:), intent(in)  :: dPdF                                         !< partitioned stiffnesses
-  integer,                              intent(in)  :: instance
+  integer,                              intent(in)  :: ho
 
-  associate(prm => param(instance))
+  associate(prm => param(ho))
 
   select case (prm%mapping)
     case (parallel_ID)

src/homogenization_thermal.f90

@@ -44,7 +44,7 @@ module subroutine thermal_init()
   allocate(current(configHomogenizations%length))
 
   do ho = 1, configHomogenizations%length
-    allocate(current(ho)%T(count(material_homogenizationAt2==ho)), source=thermal_initialT(ho))
+    allocate(current(ho)%T(count(material_homogenizationAt2==ho)), source=300.0_pReal)
     allocate(current(ho)%dot_T(count(material_homogenizationAt2==ho)), source=0.0_pReal)
     configHomogenization => configHomogenizations%get(ho)
     associate(prm => param(ho))
@@ -99,24 +99,21 @@ end subroutine thermal_homogenize
 !--------------------------------------------------------------------------------------------------
 !> @brief return homogenized thermal conductivity in reference configuration
 !--------------------------------------------------------------------------------------------------
-module function thermal_conduction_getConductivity(ip,el) result(K)
+module function thermal_conduction_getConductivity(ce) result(K)
 
-  integer, intent(in) :: &
+  integer, intent(in) :: ce
-    ip, &                                                                                           !< integration point number
-    el                                                                                              !< element number
   real(pReal), dimension(3,3) :: K
 
   integer :: &
     co
 
-
   K = 0.0_pReal
 
-  do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
+  do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce))
-    K = K + crystallite_push33ToRef(co,ip,el,lattice_K(:,:,material_phaseAt(co,el)))
+    K = K + crystallite_push33ToRef(co,ce,lattice_K(:,:,material_phaseAt2(co,ce)))
   enddo
 
-  K = K / real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
+  K = K / real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal)
 
 end function thermal_conduction_getConductivity
 
@@ -219,11 +216,11 @@ end function homogenization_thermal_T
 !--------------------------------------------------------------------------------------------------
 !> @brief return heat generation rate
 !--------------------------------------------------------------------------------------------------
-module subroutine thermal_conduction_getSource(Tdot, ip,el)
+module subroutine thermal_conduction_getSource(Tdot, ip, el)
 
   integer, intent(in) :: &
-    ip, &                                                                                           !< integration point number
+    ip, &
-    el                                                                                              !< element number
+    el
   real(pReal), intent(out) :: &
     Tdot
 

src/material.f90

@@ -17,7 +17,7 @@ module material
   private
 
   integer, dimension(:), allocatable, public, protected :: &
-    homogenization_Nconstituents                                                                 !< number of grains in each homogenization
+    homogenization_Nconstituents                                                                    !< number of grains in each homogenization
 
   character(len=:), public, protected, allocatable, dimension(:) :: &
     material_name_phase, &                                                                          !< name of each phase
@@ -30,7 +30,7 @@ module material
     material_homogenizationAt, &                                                                    !< homogenization ID of each element
     material_homogenizationAt2, &                                                                   !< per cell
     material_homogenizationMemberAt2                                                                !< cell
-  integer, dimension(:,:),   allocatable, public, protected :: &                                    ! (ip,elem)
+  integer, dimension(:,:),   allocatable :: &                                                       ! (ip,elem)
     material_homogenizationMemberAt                                                                 !< position of the element within its homogenization instance
   integer, dimension(:,:),   allocatable, public, protected :: &                                    ! (constituent,elem)
     material_phaseAt, &                                                                             !< phase ID of each element
@@ -39,9 +39,6 @@ module material
   integer, dimension(:,:,:), allocatable, public, protected :: &                                    ! (constituent,IP,elem)
     material_phaseMemberAt                                                                          !< position of the element within its phase instance
 
-  type(Rotation), dimension(:,:,:), allocatable, public, protected :: &
-    material_orientation0                                                                           !< initial orientation of each grain,IP,element
-
   public :: &
     material_init
 
@@ -125,8 +122,6 @@ subroutine material_parseMaterial
   allocate(material_phaseAt2(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0)
   allocate(material_phaseMemberAt2(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0)
 
-  allocate(material_orientation0(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems))
-
   do el = 1, discretization_Nelems
     material     => materials%get(discretization_materialAt(el))
     constituents => material%get('constituents')
@@ -135,7 +130,7 @@ subroutine material_parseMaterial
     do ip = 1, discretization_nIPs
       ce = (el-1)*discretization_nIPs + ip
       counterHomogenization(material_homogenizationAt(el)) = counterHomogenization(material_homogenizationAt(el)) + 1
-      material_homogenizationMemberAt(ip,el)                = counterHomogenization(material_homogenizationAt(el))
+      material_homogenizationMemberAt(ip,el)               = counterHomogenization(material_homogenizationAt(el))
       material_homogenizationAt2(ce)       = material_homogenizationAt(el)
       material_homogenizationMemberAt2(ce) = material_homogenizationMemberAt(ip,el)
     enddo
@@ -153,7 +148,6 @@ subroutine material_parseMaterial
 
         material_phaseAt2(co,ce)       = material_phaseAt(co,el)
         material_phaseMemberAt2(co,ce) = material_phaseMemberAt(co,ip,el)
-        call material_orientation0(co,ip,el)%fromQuaternion(constituent%get_asFloats('O',requiredSize=4)) ! should be done in crystallite
       enddo
 
     enddo

src/math.f90

@@ -1133,6 +1133,7 @@ real(pReal) pure function math_areaTriangle(v1,v2,v3)
 
   real(pReal), dimension (3), intent(in) :: v1,v2,v3
 
+
   math_areaTriangle = 0.5_pReal * norm2(math_cross(v1-v2,v1-v3))
 
 end function math_areaTriangle
@@ -1147,11 +1148,13 @@ real(pReal) pure elemental function math_clip(a, left, right)
   real(pReal), intent(in) :: a
   real(pReal), intent(in), optional :: left, right
 
+
   math_clip = a
   if (present(left))  math_clip = max(left,math_clip)
   if (present(right)) math_clip = min(right,math_clip)
-  if (present(left) .and. present(right)) &
+  if (present(left) .and. present(right)) then
-    math_clip = merge (IEEE_value(1.0_pReal,IEEE_quiet_NaN),math_clip, left>right)
+    if(left>right) error stop 'left > right'
+  endif
 
 end function math_clip
 
@@ -1182,6 +1185,7 @@ subroutine selfTest
   integer                     :: d
   logical                     :: e
 
+
   if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,3.0_pReal,3.0_pReal,3.0_pReal] - &
               math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1,2,3,0])) > tol_math_check)) &
     error stop 'math_expand [1,2,3] by [1,2,3,0] => [1,2,2,3,3,3]'

src/phase.f90

@@ -19,6 +19,9 @@ module phase
   implicit none
   private
 
+  type(Rotation), dimension(:,:,:), allocatable :: &
+    material_orientation0                                                                           !< initial orientation of each grain,IP,element
+
   type(rotation),            dimension(:,:,:),        allocatable :: &
     crystallite_orientation                                                                         !< current orientation
 
@@ -77,8 +80,8 @@ module phase
   interface
 
 ! == cleaned:begin =================================================================================
-    module subroutine mechanical_init(phases)
+    module subroutine mechanical_init(materials,phases)
-      class(tNode), pointer :: phases
+      class(tNode), pointer :: materials,phases
     end subroutine mechanical_init
 
     module subroutine damage_init
@@ -117,12 +120,11 @@ module phase
     end subroutine mechanical_restore
 
 
-    module function phase_mechanical_dPdF(dt,co,ip,el) result(dPdF)
+    module function phase_mechanical_dPdF(dt,co,ce) result(dPdF)
       real(pReal), intent(in) :: dt
       integer, intent(in) :: &
         co, &                                                                                       !< counter in constituent loop
-        ip, &                                                                                       !< counter in integration point loop
+        ce
-        el                                                                                          !< counter in element loop
       real(pReal), dimension(3,3,3,3) :: dPdF
     end function phase_mechanical_dPdF
 
@@ -147,8 +149,8 @@ module phase
       real(pReal), dimension(3,3) :: L_p
     end function mechanical_L_p
 
-    module function phase_mechanical_getF(co,ip,el) result(F)
+    module function phase_mechanical_getF(co,ce) result(F)
-      integer, intent(in) :: co, ip, el
+      integer, intent(in) :: co, ce
       real(pReal), dimension(3,3) :: F
     end function phase_mechanical_getF
 
@@ -157,8 +159,8 @@ module phase
       real(pReal), dimension(3,3) :: F_e
     end function mechanical_F_e
 
-    module function phase_mechanical_getP(co,ip,el) result(P)
+    module function phase_mechanical_getP(co,ce) result(P)
-      integer, intent(in) :: co, ip, el
+      integer, intent(in) :: co, ce
       real(pReal), dimension(3,3) :: P
     end function phase_mechanical_getP
 
@@ -183,9 +185,9 @@ module phase
     end function damage_phi
 
 
-    module subroutine phase_mechanical_setF(F,co,ip,el)
+    module subroutine phase_mechanical_setF(F,co,ce)
       real(pReal), dimension(3,3), intent(in) :: F
-      integer, intent(in) :: co, ip, el
+      integer, intent(in) :: co, ce
     end subroutine phase_mechanical_setF
 
     module subroutine phase_thermal_setField(T,dot_T, co,ce)
@@ -229,10 +231,8 @@ module phase
     end function phase_homogenizedC
 
 
-    module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el)
+    module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce)
-      integer, intent(in) :: &
+      integer, intent(in) :: ce
-        ip, &                                                                                       !< integration point number
-        el                                                                                          !< element number
       real(pReal), intent(in) :: &
         phi                                                                                         !< damage parameter
       real(pReal), intent(inout) :: &
@@ -342,6 +342,7 @@ subroutine phase_init
     so                                                                                               !< counter in source loop
   class (tNode), pointer :: &
     debug_constitutive, &
+    materials, &
     phases
 
 
@@ -356,9 +357,10 @@ subroutine phase_init
   debugConstitutive%grain      =  config_debug%get_asInt('grain',defaultVal = 1)
 
 
+  materials => config_material%get('material')
   phases => config_material%get('phase')
 
-  call mechanical_init(phases)
+  call mechanical_init(materials,phases)
   call damage_init
   call thermal_init(phases)
 
@@ -624,19 +626,20 @@ end subroutine crystallite_orientations
 !--------------------------------------------------------------------------------------------------
 !> @brief Map 2nd order tensor to reference config
 !--------------------------------------------------------------------------------------------------
-function crystallite_push33ToRef(co,ip,el, tensor33)
+function crystallite_push33ToRef(co,ce, tensor33)
 
   real(pReal), dimension(3,3), intent(in) :: tensor33
   integer, intent(in):: &
-    el, &
+    co, &
-    ip, &
+    ce
-    co
   real(pReal), dimension(3,3) :: crystallite_push33ToRef
 
   real(pReal), dimension(3,3)             :: T
+  integer :: ph, me
 
-
+  ph = material_phaseAt2(co,ce)
-  T = matmul(material_orientation0(co,ip,el)%asMatrix(),transpose(math_inv33(phase_mechanical_getF(co,ip,el)))) ! ToDo: initial orientation correct?
+  me = material_phaseMemberAt2(co,ce)
+  T = matmul(material_orientation0(co,ph,me)%asMatrix(),transpose(math_inv33(phase_mechanical_getF(co,ce)))) ! ToDo: initial orientation correct?
 
   crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T))
 

src/phase_damage.f90

@@ -133,7 +133,7 @@ module subroutine damage_init
 
   integer :: &
     ph, &                                                                                           !< counter in phase loop
-    Nconstituents
+    Nmembers
   class(tNode), pointer :: &
    phases, &
    phase, &
@@ -151,10 +151,10 @@ module subroutine damage_init
 
   do ph = 1,phases%length
 
-    Nconstituents = count(material_phaseAt2 == ph)
+    Nmembers = count(material_phaseAt2 == ph)
 
-    allocate(current(ph)%phi(Nconstituents),source=1.0_pReal)
+    allocate(current(ph)%phi(Nmembers),source=1.0_pReal)
-    allocate(current(ph)%d_phi_d_dot_phi(Nconstituents),source=0.0_pReal)
+    allocate(current(ph)%d_phi_d_dot_phi(Nmembers),source=0.0_pReal)
 
     phase => phases%get(ph)
     sources => phase%get('damage',defaultVal=emptyList)
@@ -179,11 +179,9 @@ end subroutine damage_init
 !----------------------------------------------------------------------------------------------
 !< @brief returns local part of nonlocal damage driving force
 !----------------------------------------------------------------------------------------------
-module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el)
+module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce)
 
-  integer, intent(in) :: &
+  integer, intent(in) :: ce
-    ip, &                                                                                           !< integration point number
-    el                                                                                              !< element number
   real(pReal), intent(in) :: &
     phi                                                                                             !< damage parameter
   real(pReal), intent(inout) :: &
@@ -201,9 +199,9 @@ module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi,
    phiDot = 0.0_pReal
    dPhiDot_dPhi = 0.0_pReal
 
-   do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
+   do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce))
-     ph = material_phaseAt(co,el)
+     ph = material_phaseAt2(co,ce)
-     me = material_phasememberAt(co,ip,el)
+     me = material_phasememberAt2(co,ce)
 
        select case(phase_source(ph))
          case (DAMAGE_ISOBRITTLE_ID)

src/phase_damage_anisobrittle.f90

@@ -40,7 +40,7 @@ module function anisobrittle_init() result(mySources)
     phase, &
     sources, &
     src
-  integer :: Nconstituents,p
+  integer :: Nmembers,p
   integer, dimension(:), allocatable :: N_cl
   character(len=pStringLen) :: extmsg = ''
 
@@ -92,8 +92,8 @@ module function anisobrittle_init() result(mySources)
         if (any(prm%g_crit <  0.0_pReal)) extmsg = trim(extmsg)//' g_crit'
         if (any(prm%s_crit <  0.0_pReal)) extmsg = trim(extmsg)//' s_crit'
 
-        Nconstituents = count(material_phaseAt==p) * discretization_nIPs
+        Nmembers = count(material_phaseAt==p) * discretization_nIPs
-        call phase_allocateState(damageState(p),Nconstituents,1,1,0)
+        call phase_allocateState(damageState(p),Nmembers,1,1,0)
         damageState(p)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal)
         if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol'
 

src/phase_damage_anisoductile.f90

@@ -35,7 +35,7 @@ module function anisoductile_init() result(mySources)
     pl, &
     sources, &
     src
-  integer :: Ninstances,Nconstituents,p
+  integer :: Ninstances,Nmembers,p
   integer, dimension(:), allocatable :: N_sl
   character(len=pStringLen) :: extmsg = ''
 
@@ -78,8 +78,8 @@ module function anisoductile_init() result(mySources)
         if (prm%q              <= 0.0_pReal)  extmsg = trim(extmsg)//' q'
         if (any(prm%gamma_crit <  0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit'
 
-        Nconstituents=count(material_phaseAt2==p)
+        Nmembers=count(material_phaseAt2==p)
-        call phase_allocateState(damageState(p),Nconstituents,1,1,0)
+        call phase_allocateState(damageState(p),Nmembers,1,1,0)
         damageState(p)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal)
         if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol'
 

src/phase_damage_isobrittle.f90

@@ -31,7 +31,7 @@ module function isobrittle_init() result(mySources)
     phase, &
     sources, &
     src
-  integer :: Nconstituents,p
+  integer :: Nmembers,p
   character(len=pStringLen) :: extmsg = ''
 
 
@@ -64,8 +64,8 @@ module function isobrittle_init() result(mySources)
         ! sanity checks
         if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit'
 
-        Nconstituents = count(material_phaseAt2==p)
+        Nmembers = count(material_phaseAt2==p)
-        call phase_allocateState(damageState(p),Nconstituents,1,1,1)
+        call phase_allocateState(damageState(p),Nmembers,1,1,1)
         damageState(p)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal)
         if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol'
 

src/phase_damage_isoductile.f90

@@ -33,7 +33,7 @@ module function isoductile_init() result(mySources)
     phase, &
     sources, &
     src
-  integer :: Ninstances,Nconstituents,p
+  integer :: Ninstances,Nmembers,p
   character(len=pStringLen) :: extmsg = ''
 
 
@@ -68,8 +68,8 @@ module function isoductile_init() result(mySources)
         if (prm%q          <= 0.0_pReal) extmsg = trim(extmsg)//' q'
         if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit'
 
-        Nconstituents=count(material_phaseAt2==p)
+        Nmembers=count(material_phaseAt2==p)
-        call phase_allocateState(damageState(p),Nconstituents,1,1,0)
+        call phase_allocateState(damageState(p),Nmembers,1,1,0)
         damageState(p)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal)
         if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol'
 

src/phase_mechanical.f90

@@ -177,21 +177,26 @@ contains
 !> @brief Initialize mechanical field related constitutive models
 !> @details Initialize elasticity, plasticity and stiffness degradation models.
 !--------------------------------------------------------------------------------------------------
-module subroutine mechanical_init(phases)
+module subroutine mechanical_init(materials,phases)
 
   class(tNode), pointer :: &
+    materials, &
     phases
 
   integer :: &
     el, &
     ip, &
     co, &
+    ce, &
     ph, &
     me, &
     stiffDegradationCtr, &
-    Nconstituents
+    Nmembers
   class(tNode), pointer :: &
     num_crystallite, &
+    material, &
+    constituents, &
+    constituent, &
     phase, &
     mech, &
     elastic, &
@@ -221,23 +226,25 @@ module subroutine mechanical_init(phases)
   allocate(phase_mechanical_P(phases%length))
   allocate(phase_mechanical_S0(phases%length))
 
-  do ph = 1, phases%length
+  allocate(material_orientation0(homogenization_maxNconstituents,phases%length,maxVal(material_phaseMemberAt)))
-    Nconstituents = count(material_phaseAt == ph) * discretization_nIPs
 
-    allocate(phase_mechanical_Fi(ph)%data(3,3,Nconstituents))
+  do ph = 1, phases%length
-    allocate(phase_mechanical_Fe(ph)%data(3,3,Nconstituents))
+    Nmembers = count(material_phaseAt == ph) * discretization_nIPs
-    allocate(phase_mechanical_Fi0(ph)%data(3,3,Nconstituents))
+
-    allocate(phase_mechanical_Fp(ph)%data(3,3,Nconstituents))
+    allocate(phase_mechanical_Fi(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_Fp0(ph)%data(3,3,Nconstituents))
+    allocate(phase_mechanical_Fe(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_Li(ph)%data(3,3,Nconstituents))
+    allocate(phase_mechanical_Fi0(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_Li0(ph)%data(3,3,Nconstituents))
+    allocate(phase_mechanical_Fp(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_Lp0(ph)%data(3,3,Nconstituents))
+    allocate(phase_mechanical_Fp0(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_Lp(ph)%data(3,3,Nconstituents))
+    allocate(phase_mechanical_Li(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_S(ph)%data(3,3,Nconstituents),source=0.0_pReal)
+    allocate(phase_mechanical_Li0(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_P(ph)%data(3,3,Nconstituents),source=0.0_pReal)
+    allocate(phase_mechanical_Lp0(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_S0(ph)%data(3,3,Nconstituents),source=0.0_pReal)
+    allocate(phase_mechanical_Lp(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_F(ph)%data(3,3,Nconstituents))
+    allocate(phase_mechanical_S(ph)%data(3,3,Nmembers),source=0.0_pReal)
-    allocate(phase_mechanical_F0(ph)%data(3,3,Nconstituents))
+    allocate(phase_mechanical_P(ph)%data(3,3,Nmembers),source=0.0_pReal)
+    allocate(phase_mechanical_S0(ph)%data(3,3,Nmembers),source=0.0_pReal)
+    allocate(phase_mechanical_F(ph)%data(3,3,Nmembers))
+    allocate(phase_mechanical_F0(ph)%data(3,3,Nmembers))
 
     phase   => phases%get(ph)
     mech    => phase%get('mechanics')
@@ -271,14 +278,19 @@ module subroutine mechanical_init(phases)
     enddo
   endif
 
-  !$OMP PARALLEL DO PRIVATE(ph,me)
+
   do el = 1, size(material_phaseMemberAt,3); do ip = 1, size(material_phaseMemberAt,2)
     do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
+      material     => materials%get(discretization_materialAt(el))
+      constituents => material%get('constituents')
+      constituent => constituents%get(co)
 
       ph = material_phaseAt(co,el)
       me = material_phaseMemberAt(co,ip,el)
 
-      phase_mechanical_Fp0(ph)%data(1:3,1:3,me) = material_orientation0(co,ip,el)%asMatrix()                      ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005)
+      call material_orientation0(co,ph,me)%fromQuaternion(constituent%get_asFloats('O',requiredSize=4))
+
+      phase_mechanical_Fp0(ph)%data(1:3,1:3,me) = material_orientation0(co,ph,me)%asMatrix()                         ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005)
       phase_mechanical_Fp0(ph)%data(1:3,1:3,me) = phase_mechanical_Fp0(ph)%data(1:3,1:3,me) &
                                                 / math_det33(phase_mechanical_Fp0(ph)%data(1:3,1:3,me))**(1.0_pReal/3.0_pReal)
       phase_mechanical_Fi0(ph)%data(1:3,1:3,me) = math_I3
@@ -292,7 +304,6 @@ module subroutine mechanical_init(phases)
 
     enddo
   enddo; enddo
-  !$OMP END PARALLEL DO
 
 
 ! initialize plasticity
@@ -342,12 +353,11 @@ end subroutine mechanical_init
 !> the elastic and intermediate deformation gradients using Hooke's law
 !--------------------------------------------------------------------------------------------------
 subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
-                                              Fe, Fi, co, ip, el)
+                                              Fe, Fi, ph, me)
 
   integer, intent(in) :: &
-    co, &                                                                                          !< component-ID of integration point
+    ph, &
-    ip, &                                                                                           !< integration point
+    me
-    el                                                                                              !< element
   real(pReal),   intent(in),  dimension(3,3) :: &
     Fe, &                                                                                           !< elastic deformation gradient
     Fi                                                                                              !< intermediate deformation gradient
@@ -360,17 +370,15 @@ subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
   real(pReal), dimension(3,3) :: E
   real(pReal), dimension(3,3,3,3) :: C
   integer :: &
-    ho, &                                                                                           !< homogenization
+    d, &                                                                                            !< counter in degradation loop
-    d, &                                                                                               !< counter in degradation loop
+    i, j
-    i, j, ph, me
 
-  ho = material_homogenizationAt(el)
+  C = math_66toSym3333(phase_homogenizedC(ph,me))
-  C = math_66toSym3333(phase_homogenizedC(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el)))
 
-  DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phaseAt(co,el))
+  DegradationLoop: do d = 1, phase_NstiffnessDegradations(ph)
-    degradationType: select case(phase_stiffnessDegradation(d,material_phaseAt(co,el)))
+    degradationType: select case(phase_stiffnessDegradation(d,ph))
       case (STIFFNESS_DEGRADATION_damage_ID) degradationType
-        C = C * phase_damage_get_phi(co,ip,el)**2
+        C = C * damage_phi(ph,me)**2
     end select degradationType
   enddo DegradationLoop
 
@@ -528,7 +536,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
       B  = math_I3 - Delta_t*Lpguess
       Fe = matmul(matmul(A,B), invFi_new)
       call phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
-                                        Fe, Fi_new, co, ip, el)
+                                        Fe, Fi_new, ph, me)
 
       call plastic_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, &
                                          S, Fi_new, ph,me)
@@ -1250,13 +1258,12 @@ end subroutine mechanical_restore
 !--------------------------------------------------------------------------------------------------
 !> @brief Calculate tangent (dPdF).
 !--------------------------------------------------------------------------------------------------
-module function phase_mechanical_dPdF(dt,co,ip,el) result(dPdF)
+module function phase_mechanical_dPdF(dt,co,ce) result(dPdF)
 
   real(pReal), intent(in) :: dt
   integer, intent(in) :: &
     co, &                                                                                            !< counter in constituent loop
-    ip, &                                                                                            !< counter in integration point loop
+    ce
-    el                                                                                               !< counter in element loop
   real(pReal), dimension(3,3,3,3) :: dPdF
 
   integer :: &
@@ -1281,12 +1288,12 @@ module function phase_mechanical_dPdF(dt,co,ip,el) result(dPdF)
   logical :: error
 
 
-  ph = material_phaseAt(co,el)
+  ph = material_phaseAt2(co,ce)
-  me = material_phaseMemberAt(co,ip,el)
+  me = material_phaseMemberAt2(co,ce)
 
   call phase_hooke_SandItsTangents(devNull,dSdFe,dSdFi, &
                                           phase_mechanical_Fe(ph)%data(1:3,1:3,me), &
-                                          phase_mechanical_Fi(ph)%data(1:3,1:3,me),co,ip,el)
+                                          phase_mechanical_Fi(ph)%data(1:3,1:3,me),ph,me)
   call phase_LiAndItsTangents(devNull,dLidS,dLidFi, &
                                      phase_mechanical_S(ph)%data(1:3,1:3,me), &
                                      phase_mechanical_Fi(ph)%data(1:3,1:3,me), &
@@ -1311,7 +1318,7 @@ module function phase_mechanical_dPdF(dt,co,ip,el) result(dPdF)
     enddo; enddo
     call math_invert(temp_99,error,math_3333to99(lhs_3333))
     if (error) then
-      call IO_warning(warning_ID=600,el=el,ip=ip,g=co, &
+      call IO_warning(warning_ID=600, &
                       ext_msg='inversion error in analytic tangent calculation')
       dFidS = 0.0_pReal
     else
@@ -1341,7 +1348,7 @@ module function phase_mechanical_dPdF(dt,co,ip,el) result(dPdF)
 
   call math_invert(temp_99,error,math_eye(9)+math_3333to99(lhs_3333))
   if (error) then
-    call IO_warning(warning_ID=600,el=el,ip=ip,g=co, &
+    call IO_warning(warning_ID=600, &
                     ext_msg='inversion error in analytic tangent calculation')
     dSdF = rhs_3333
   else
@@ -1440,13 +1447,13 @@ end function mechanical_L_p
 !----------------------------------------------------------------------------------------------
 !< @brief Get deformation gradient (for use by homogenization)
 !----------------------------------------------------------------------------------------------
-module function phase_mechanical_getF(co,ip,el) result(F)
+module function phase_mechanical_getF(co,ce) result(F)
 
-  integer, intent(in) :: co, ip, el
+  integer, intent(in) :: co, ce
   real(pReal), dimension(3,3) :: F
 
 
-  F = phase_mechanical_F(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el))
+  F = phase_mechanical_F(material_phaseAt2(co,ce))%data(1:3,1:3,material_phaseMemberAt2(co,ce))
 
 end function phase_mechanical_getF
 
@@ -1469,25 +1476,25 @@ end function mechanical_F_e
 !----------------------------------------------------------------------------------------------
 !< @brief Get second Piola-Kichhoff stress (for use by homogenization)
 !----------------------------------------------------------------------------------------------
-module function phase_mechanical_getP(co,ip,el) result(P)
+module function phase_mechanical_getP(co,ce) result(P)
 
-  integer, intent(in) :: co, ip, el
+  integer, intent(in) :: co, ce
   real(pReal), dimension(3,3) :: P
 
 
-  P = phase_mechanical_P(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el))
+  P = phase_mechanical_P(material_phaseAt2(co,ce))%data(1:3,1:3,material_phaseMemberAt2(co,ce))
 
 end function phase_mechanical_getP
 
 
 ! setter for homogenization
-module subroutine phase_mechanical_setF(F,co,ip,el)
+module subroutine phase_mechanical_setF(F,co,ce)
 
   real(pReal), dimension(3,3), intent(in) :: F
-  integer, intent(in) :: co, ip, el
+  integer, intent(in) :: co, ce
 
 
-  phase_mechanical_F(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) = F
+  phase_mechanical_F(material_phaseAt2(co,ce))%data(1:3,1:3,material_phaseMemberAt2(co,ce)) = F
 
 end subroutine phase_mechanical_setF
 

src/phase_mechanical_plastic_dislotungsten.f90

@@ -79,7 +79,7 @@ module function plastic_dislotungsten_init() result(myPlasticity)
   logical, dimension(:), allocatable :: myPlasticity
   integer :: &
     ph, i, &
-    Nconstituents, &
+    Nmembers, &
     sizeState, sizeDotState, &
     startIndex, endIndex
   integer,    dimension(:), allocatable :: &
@@ -220,18 +220,18 @@ module function plastic_dislotungsten_init() result(myPlasticity)
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
-    Nconstituents = count(material_phaseAt2 == ph)
+    Nmembers = count(material_phaseAt2 == ph)
     sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl
     sizeState = sizeDotState
 
-    call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0)
+    call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
 
 !--------------------------------------------------------------------------------------------------
 ! state aliases and initialization
     startIndex = 1
     endIndex   = prm%sum_N_sl
     stt%rho_mob => plasticState(ph)%state(startIndex:endIndex,:)
-    stt%rho_mob =  spread(rho_mob_0,2,Nconstituents)
+    stt%rho_mob =  spread(rho_mob_0,2,Nmembers)
     dot%rho_mob => plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
     if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho'
@@ -239,7 +239,7 @@ module function plastic_dislotungsten_init() result(myPlasticity)
     startIndex = endIndex + 1
     endIndex   = endIndex + prm%sum_N_sl
     stt%rho_dip => plasticState(ph)%state(startIndex:endIndex,:)
-    stt%rho_dip =  spread(rho_dip_0,2,Nconstituents)
+    stt%rho_dip =  spread(rho_dip_0,2,Nmembers)
     dot%rho_dip => plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
 
@@ -251,8 +251,8 @@ module function plastic_dislotungsten_init() result(myPlasticity)
     ! global alias
     plasticState(ph)%slipRate        => plasticState(ph)%dotState(startIndex:endIndex,:)
 
-    allocate(dst%Lambda_sl(prm%sum_N_sl,Nconstituents),         source=0.0_pReal)
+    allocate(dst%Lambda_sl(prm%sum_N_sl,Nmembers),         source=0.0_pReal)
-    allocate(dst%threshold_stress(prm%sum_N_sl,Nconstituents),  source=0.0_pReal)
+    allocate(dst%threshold_stress(prm%sum_N_sl,Nmembers),  source=0.0_pReal)
 
     plasticState(ph)%state0 = plasticState(ph)%state                                                ! ToDo: this could be done centrally
 

src/phase_mechanical_plastic_dislotwin.f90

@@ -24,7 +24,6 @@ submodule(phase:plastic) dislotwin
       q_sb                = 1.0_pReal, &                                                            !< q-exponent in shear band velocity
       D_a                 = 1.0_pReal, &                                                            !< adjustment parameter to calculate minimum dipole distance
       i_tw                = 1.0_pReal, &                                                            !< adjustment parameter to calculate MFP for twinning
-      tau_0               = 1.0_pReal, &                                                            !< strength due to elements in solid solution
       L_tw                = 1.0_pReal, &                                                            !< Length of twin nuclei in Burgers vectors
       L_tr                = 1.0_pReal, &                                                            !< Length of trans nuclei in Burgers vectors
       x_c_tw              = 1.0_pReal, &                                                            !< critical distance for formation of twin nucleus
@@ -53,6 +52,7 @@ submodule(phase:plastic) dislotwin
       q, &                                                                                          !< q-exponent in glide velocity
       r, &                                                                                          !< r-exponent in twin nucleation rate
       s, &                                                                                          !< s-exponent in trans nucleation rate
+      tau_0, &                                                                                      !< strength due to elements in solid solution
       gamma_char, &                                                                                 !< characteristic shear for twins
       B                                                                                             !< drag coefficient
     real(pReal),               allocatable, dimension(:,:) :: &
@@ -81,7 +81,7 @@ submodule(phase:plastic) dislotwin
     logical :: &
       ExtendedDislocations, &                                                                       !< consider split into partials for climb calculation
       fccTwinTransNucleation, &                                                                     !< twinning and transformation models are for fcc
-      dipoleFormation                                                                               !< flag indicating consideration of dipole formation
+      omitDipoles                                                                                   !< flag controlling consideration of dipole formation
   end type                                                                                          !< container type for internal constitutive parameters
 
   type :: tDislotwinState
@@ -127,7 +127,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
   logical, dimension(:), allocatable :: myPlasticity
   integer :: &
     ph, i, &
-    Nconstituents, &
+    Nmembers, &
     sizeState, sizeDotState, &
     startIndex, endIndex
   integer,     dimension(:), allocatable :: &
@@ -213,10 +213,10 @@ module function plastic_dislotwin_init() result(myPlasticity)
       prm%i_sl                 = pl%get_asFloats('i_sl',        requiredSize=size(N_sl))
       prm%p                    = pl%get_asFloats('p_sl',        requiredSize=size(N_sl))
       prm%q                    = pl%get_asFloats('q_sl',        requiredSize=size(N_sl))
+      prm%tau_0                = pl%get_asFloats('tau_0',       requiredSize=size(N_sl))
       prm%B                    = pl%get_asFloats('B',           requiredSize=size(N_sl), &
                                                   defaultVal=[(0.0_pReal, i=1,size(N_sl))])
 
-      prm%tau_0                = pl%get_asFloat('tau_0')
       prm%D_a                  = pl%get_asFloat('D_a')
       prm%D_0                  = pl%get_asFloat('D_0')
       prm%Q_cl                 = pl%get_asFloat('Q_cl')
@@ -226,7 +226,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
         prm%dGamma_sf_dT       = pl%get_asFloat('dGamma_sf_dT')
       endif
 
-      prm%dipoleformation = .not. pl%get_asBool('no_dipole_formation',defaultVal = .false.)
+      prm%omitDipoles = pl%get_asBool('omit_dipoles',defaultVal = .false.)
 
       ! multiplication factor according to crystal structure (nearest neighbors bcc vs fcc/hex)
       ! details: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981
@@ -242,6 +242,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
       prm%i_sl         = math_expand(prm%i_sl,        N_sl)
       prm%p            = math_expand(prm%p,           N_sl)
       prm%q            = math_expand(prm%q,           N_sl)
+      prm%tau_0        = math_expand(prm%tau_0,       N_sl)
       prm%B            = math_expand(prm%B,           N_sl)
 
       ! sanity checks
@@ -405,21 +406,21 @@ module function plastic_dislotwin_init() result(myPlasticity)
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
-    Nconstituents  = count(material_phaseAt2 == ph)
+    Nmembers  = count(material_phaseAt2 == ph)
     sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl &
                  + size(['f_tw'])                           * prm%sum_N_tw &
                  + size(['f_tr'])                           * prm%sum_N_tr
     sizeState = sizeDotState
 
 
-    call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0)
+    call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
 
 !--------------------------------------------------------------------------------------------------
 ! locally defined state aliases and initialization of state0 and atol
     startIndex = 1
     endIndex   = prm%sum_N_sl
     stt%rho_mob=>plasticState(ph)%state(startIndex:endIndex,:)
-    stt%rho_mob= spread(rho_mob_0,2,Nconstituents)
+    stt%rho_mob= spread(rho_mob_0,2,Nmembers)
     dot%rho_mob=>plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
     if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho'
@@ -427,7 +428,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
     startIndex = endIndex + 1
     endIndex   = endIndex + prm%sum_N_sl
     stt%rho_dip=>plasticState(ph)%state(startIndex:endIndex,:)
-    stt%rho_dip= spread(rho_dip_0,2,Nconstituents)
+    stt%rho_dip= spread(rho_dip_0,2,Nmembers)
     dot%rho_dip=>plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
 
@@ -443,28 +444,28 @@ module function plastic_dislotwin_init() result(myPlasticity)
     endIndex   = endIndex + prm%sum_N_tw
     stt%f_tw=>plasticState(ph)%state(startIndex:endIndex,:)
     dot%f_tw=>plasticState(ph)%dotState(startIndex:endIndex,:)
-    plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('f_twin',defaultVal=1.0e-7_pReal)
+    plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_f_tw',defaultVal=1.0e-7_pReal)
-    if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' f_twin'
+    if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_f_tw'
 
     startIndex = endIndex + 1
     endIndex   = endIndex + prm%sum_N_tr
     stt%f_tr=>plasticState(ph)%state(startIndex:endIndex,:)
     dot%f_tr=>plasticState(ph)%dotState(startIndex:endIndex,:)
-    plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('f_trans',defaultVal=1.0e-6_pReal)
+    plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_f_tr',defaultVal=1.0e-6_pReal)
-    if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' f_trans'
+    if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_f_tr'
 
-    allocate(dst%Lambda_sl             (prm%sum_N_sl,Nconstituents),source=0.0_pReal)
+    allocate(dst%Lambda_sl             (prm%sum_N_sl,Nmembers),source=0.0_pReal)
-    allocate(dst%tau_pass              (prm%sum_N_sl,Nconstituents),source=0.0_pReal)
+    allocate(dst%tau_pass              (prm%sum_N_sl,Nmembers),source=0.0_pReal)
 
-    allocate(dst%Lambda_tw             (prm%sum_N_tw,Nconstituents),source=0.0_pReal)
+    allocate(dst%Lambda_tw             (prm%sum_N_tw,Nmembers),source=0.0_pReal)
-    allocate(dst%tau_hat_tw            (prm%sum_N_tw,Nconstituents),source=0.0_pReal)
+    allocate(dst%tau_hat_tw            (prm%sum_N_tw,Nmembers),source=0.0_pReal)
-    allocate(dst%tau_r_tw              (prm%sum_N_tw,Nconstituents),source=0.0_pReal)
+    allocate(dst%tau_r_tw              (prm%sum_N_tw,Nmembers),source=0.0_pReal)
-    allocate(dst%V_tw                  (prm%sum_N_tw,Nconstituents),source=0.0_pReal)
+    allocate(dst%V_tw                  (prm%sum_N_tw,Nmembers),source=0.0_pReal)
 
-    allocate(dst%Lambda_tr             (prm%sum_N_tr,Nconstituents),source=0.0_pReal)
+    allocate(dst%Lambda_tr             (prm%sum_N_tr,Nmembers),source=0.0_pReal)
-    allocate(dst%tau_hat_tr            (prm%sum_N_tr,Nconstituents),source=0.0_pReal)
+    allocate(dst%tau_hat_tr            (prm%sum_N_tr,Nmembers),source=0.0_pReal)
-    allocate(dst%tau_r_tr              (prm%sum_N_tr,Nconstituents),source=0.0_pReal)
+    allocate(dst%tau_r_tr              (prm%sum_N_tr,Nmembers),source=0.0_pReal)
-    allocate(dst%V_tr                  (prm%sum_N_tr,Nconstituents),source=0.0_pReal)
+    allocate(dst%V_tr                  (prm%sum_N_tr,Nmembers),source=0.0_pReal)
 
     plasticState(ph)%state0 = plasticState(ph)%state                                                ! ToDo: this could be done centrally
 
@@ -535,9 +536,9 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
     dot_gamma_sl,ddot_gamma_dtau_slip
   real(pReal), dimension(param(ph)%sum_N_tw) :: &
-    dot_gamma_twin,ddot_gamma_dtau_twin
+    dot_gamma_tw,ddot_gamma_dtau_tw
   real(pReal), dimension(param(ph)%sum_N_tr) :: &
-    dot_gamma_tr,ddot_gamma_dtau_trans
+    dot_gamma_tr,ddot_gamma_dtau_tr
   real(pReal):: dot_gamma_sb
   real(pReal), dimension(3,3) :: eigVectors, P_sb
   real(pReal), dimension(3)   :: eigValues
@@ -578,20 +579,20 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,me)
                        + ddot_gamma_dtau_slip(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i)
   enddo slipContribution
 
-  call kinetics_twin(Mp,T,dot_gamma_sl,ph,me,dot_gamma_twin,ddot_gamma_dtau_twin)
+  call kinetics_twin(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tw,ddot_gamma_dtau_tw)
   twinContibution: do i = 1, prm%sum_N_tw
-    Lp = Lp + dot_gamma_twin(i)*prm%P_tw(1:3,1:3,i)
+    Lp = Lp + dot_gamma_tw(i)*prm%P_tw(1:3,1:3,i)
     forall (k=1:3,l=1:3,m=1:3,n=1:3) &
       dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
-                       + ddot_gamma_dtau_twin(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
+                       + ddot_gamma_dtau_tw(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
   enddo twinContibution
 
-  call kinetics_trans(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tr,ddot_gamma_dtau_trans)
+  call kinetics_trans(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tr,ddot_gamma_dtau_tr)
   transContibution: do i = 1, prm%sum_N_tr
     Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i)
     forall (k=1:3,l=1:3,m=1:3,n=1:3) &
       dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
-                       + ddot_gamma_dtau_trans(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i)
+                       + ddot_gamma_dtau_tr(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i)
   enddo transContibution
 
   Lp      = Lp      * f_unrotated
@@ -646,7 +647,6 @@ module subroutine dislotwin_dotState(Mp,T,ph,me)
     f_unrotated, &
     rho_dip_distance, &
     v_cl, &                                                                                         !< climb velocity
-    Gamma, &                                                                                        !< stacking fault energy
     tau, &
     sigma_cl, &                                                                                     !< climb stress
     b_d                                                                                             !< ratio of Burgers vector to stacking fault width
@@ -656,7 +656,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,me)
     rho_dip_distance_min, &
     dot_gamma_sl
   real(pReal), dimension(param(ph)%sum_N_tw) :: &
-    dot_gamma_twin
+    dot_gamma_tw
   real(pReal), dimension(param(ph)%sum_N_tr) :: &
     dot_gamma_tr
 
@@ -675,7 +675,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,me)
   slipState: do i = 1, prm%sum_N_sl
     tau = math_tensordot(Mp,prm%P_sl(1:3,1:3,i))
 
-    significantSlipStress: if (dEq0(tau)) then
+    significantSlipStress: if (dEq0(tau) .or. prm%omitDipoles) then
       dot_rho_dip_formation(i) = 0.0_pReal
       dot_rho_dip_climb(i) = 0.0_pReal
     else significantSlipStress
@@ -683,24 +683,18 @@ module subroutine dislotwin_dotState(Mp,T,ph,me)
       rho_dip_distance = math_clip(rho_dip_distance, right = dst%Lambda_sl(i,me))
       rho_dip_distance = math_clip(rho_dip_distance, left  = rho_dip_distance_min(i))
 
-      if (prm%dipoleFormation) then
+      dot_rho_dip_formation(i) = 2.0_pReal*(rho_dip_distance-rho_dip_distance_min(i))/prm%b_sl(i) &
-        dot_rho_dip_formation(i) = 2.0_pReal*(rho_dip_distance-rho_dip_distance_min(i))/prm%b_sl(i) &
+                               * stt%rho_mob(i,me)*abs(dot_gamma_sl(i))
-                                 * stt%rho_mob(i,me)*abs(dot_gamma_sl(i))
-      else
-        dot_rho_dip_formation(i) = 0.0_pReal
-      endif
 
       if (dEq(rho_dip_distance,rho_dip_distance_min(i))) then
         dot_rho_dip_climb(i) = 0.0_pReal
       else
-      !@details: Refer: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981
+        ! Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981
         sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i)))
-        if (prm%ExtendedDislocations) then
+        b_d = merge(24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu) &
-          Gamma = prm%Gamma_sf_0K + prm%dGamma_sf_dT * T
+                      * (prm%Gamma_sf_0K + prm%dGamma_sf_dT * T) / (prm%mu*prm%b_sl(i)), &
-          b_d = 24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu)* Gamma/(prm%mu*prm%b_sl(i))
+                    1.0_pReal, &
-        else
+                    prm%ExtendedDislocations)
-          b_d = 1.0_pReal
-        endif
         v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(kB*T)) &
              * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal)
 
@@ -718,8 +712,8 @@ module subroutine dislotwin_dotState(Mp,T,ph,me)
                     - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,me)*abs(dot_gamma_sl) &
                     - dot_rho_dip_climb
 
-  call kinetics_twin(Mp,T,dot_gamma_sl,ph,me,dot_gamma_twin)
+  call kinetics_twin(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tw)
-  dot%f_tw(:,me) = f_unrotated*dot_gamma_twin/prm%gamma_char
+  dot%f_tw(:,me) = f_unrotated*dot_gamma_tw/prm%gamma_char
 
   call kinetics_trans(Mp,T,dot_gamma_sl,ph,me,dot_gamma_tr)
   dot%f_tr(:,me) = f_unrotated*dot_gamma_tr
@@ -741,11 +735,9 @@ module subroutine dislotwin_dependentState(T,ph,me)
     T
 
   real(pReal) :: &
-    sumf_twin,Gamma,sumf_trans
+    sumf_tw,Gamma,sumf_tr
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
-    inv_lambda_sl_sl, &                                                                             !< 1/mean free distance between 2 forest dislocations seen by a moving dislocation
+    inv_lambda_sl
-    inv_lambda_sl_tw, &                                                                             !< 1/mean free distance between 2 twin stacks from different systems seen by a moving dislocation
-    inv_lambda_sl_tr                                                                                !< 1/mean free distance between 2 martensite lamellar from different systems seen by a moving dislocation
   real(pReal), dimension(param(ph)%sum_N_tw) :: &
     inv_lambda_tw_tw, &                                                                             !< 1/mean free distance between 2 twin stacks from different systems seen by a growing twin
     f_over_t_tw
@@ -760,38 +752,27 @@ module subroutine dislotwin_dependentState(T,ph,me)
             stt => state(ph),&
             dst => dependentState(ph))
 
-  sumf_twin  = sum(stt%f_tw(1:prm%sum_N_tw,me))
+  sumf_tw  = sum(stt%f_tw(1:prm%sum_N_tw,me))
-  sumf_trans = sum(stt%f_tr(1:prm%sum_N_tr,me))
+  sumf_tr = sum(stt%f_tr(1:prm%sum_N_tr,me))
 
   Gamma = prm%Gamma_sf_0K + prm%dGamma_sf_dT * T
 
   !* rescaled volume fraction for topology
   f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,me)/prm%t_tw                                                ! this is per system ...
-  f_over_t_tr = sumf_trans/prm%t_tr                                                                 ! but this not
+  f_over_t_tr = sumf_tr/prm%t_tr                                                                    ! but this not
                                                                                                     ! ToDo ...Physically correct, but naming could be adjusted
 
-  inv_lambda_sl_sl = sqrt(matmul(prm%forestProjection, &
+  inv_lambda_sl = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,me)+stt%rho_dip(:,me)))/prm%i_sl
-                                 stt%rho_mob(:,me)+stt%rho_dip(:,me)))/prm%i_sl
-
   if (prm%sum_N_tw > 0 .and. prm%sum_N_sl > 0) &
-    inv_lambda_sl_tw = matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_twin)
+    inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_tw)
-
-  inv_lambda_tw_tw = matmul(prm%h_tw_tw,f_over_t_tw)/(1.0_pReal-sumf_twin)
-
   if (prm%sum_N_tr > 0 .and. prm%sum_N_sl > 0) &
-    inv_lambda_sl_tr = matmul(prm%h_sl_tr,f_over_t_tr)/(1.0_pReal-sumf_trans)
+    inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tr,f_over_t_tr)/(1.0_pReal-sumf_tr)
-
+  dst%Lambda_sl(:,me) = prm%D / (1.0_pReal+prm%D*inv_lambda_sl)
-  inv_lambda_tr_tr = matmul(prm%h_tr_tr,f_over_t_tr)/(1.0_pReal-sumf_trans)
-
-  if ((prm%sum_N_tw > 0) .or. (prm%sum_N_tr > 0)) then                                              ! ToDo: better logic needed here
-    dst%Lambda_sl(:,me) = prm%D &
-                        / (1.0_pReal+prm%D*(inv_lambda_sl_sl + inv_lambda_sl_tw + inv_lambda_sl_tr))
-  else
-    dst%Lambda_sl(:,me) = prm%D &
-                        / (1.0_pReal+prm%D*inv_lambda_sl_sl) !!!!!! correct?
-  endif
 
+  inv_lambda_tw_tw = matmul(prm%h_tw_tw,f_over_t_tw)/(1.0_pReal-sumf_tw)
   dst%Lambda_tw(:,me) = prm%i_tw*prm%D/(1.0_pReal+prm%D*inv_lambda_tw_tw)
+
+  inv_lambda_tr_tr = matmul(prm%h_tr_tr,f_over_t_tr)/(1.0_pReal-sumf_tr)
   dst%Lambda_tr(:,me) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr)
 
   !* threshold stress for dislocation motion
@@ -957,7 +938,7 @@ end subroutine kinetics_slip
 ! have the optional arguments at the end.
 !--------------------------------------------------------------------------------------------------
 pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,me,&
-                              dot_gamma_twin,ddot_gamma_dtau_twin)
+                              dot_gamma_tw,ddot_gamma_dtau_tw)
 
   real(pReal), dimension(3,3),  intent(in) :: &
     Mp                                                                                              !< Mandel stress
@@ -970,9 +951,9 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,me,&
     dot_gamma_sl
 
   real(pReal), dimension(param(ph)%sum_N_tw), intent(out) :: &
-    dot_gamma_twin
+    dot_gamma_tw
   real(pReal), dimension(param(ph)%sum_N_tw), optional, intent(out) :: &
-    ddot_gamma_dtau_twin
+    ddot_gamma_dtau_tw
 
   real, dimension(param(ph)%sum_N_tw) :: &
     tau, &
@@ -1004,16 +985,16 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,me,&
 
   significantStress: where(tau > tol_math_check)
     StressRatio_r   = (dst%tau_hat_tw(:,me)/tau)**prm%r
-    dot_gamma_twin  = prm%gamma_char * dst%V_tw(:,me) * Ndot0*exp(-StressRatio_r)
+    dot_gamma_tw    = prm%gamma_char * dst%V_tw(:,me) * Ndot0*exp(-StressRatio_r)
-    ddot_gamma_dtau = (dot_gamma_twin*prm%r/tau)*StressRatio_r
+    ddot_gamma_dtau = (dot_gamma_tw*prm%r/tau)*StressRatio_r
   else where significantStress
-    dot_gamma_twin  = 0.0_pReal
+    dot_gamma_tw    = 0.0_pReal
     ddot_gamma_dtau = 0.0_pReal
   end where significantStress
 
   end associate
 
-  if(present(ddot_gamma_dtau_twin)) ddot_gamma_dtau_twin = ddot_gamma_dtau
+  if(present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw = ddot_gamma_dtau
 
 end subroutine kinetics_twin
 
@@ -1026,7 +1007,7 @@ end subroutine kinetics_twin
 ! have the optional arguments at the end.
 !--------------------------------------------------------------------------------------------------
 pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,me,&
-                              dot_gamma_tr,ddot_gamma_dtau_trans)
+                              dot_gamma_tr,ddot_gamma_dtau_tr)
 
   real(pReal), dimension(3,3),  intent(in) :: &
     Mp                                                                                              !< Mandel stress
@@ -1041,7 +1022,7 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,me,&
   real(pReal), dimension(param(ph)%sum_N_tr), intent(out) :: &
     dot_gamma_tr
   real(pReal), dimension(param(ph)%sum_N_tr), optional, intent(out) :: &
-    ddot_gamma_dtau_trans
+    ddot_gamma_dtau_tr
 
   real, dimension(param(ph)%sum_N_tr) :: &
     tau, &
@@ -1081,7 +1062,7 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,me,&
 
   end associate
 
-  if(present(ddot_gamma_dtau_trans)) ddot_gamma_dtau_trans = ddot_gamma_dtau
+  if(present(ddot_gamma_dtau_tr)) ddot_gamma_dtau_tr = ddot_gamma_dtau
 
 end subroutine kinetics_trans
 

src/phase_mechanical_plastic_isotropic.f90

@@ -52,7 +52,7 @@ module function plastic_isotropic_init() result(myPlasticity)
   logical, dimension(:), allocatable :: myPlasticity
   integer :: &
     ph, &
-    Nconstituents, &
+    Nmembers, &
     sizeState, sizeDotState
   real(pReal) :: &
     xi_0                                                                                            !< initial critical stress
@@ -119,11 +119,11 @@ module function plastic_isotropic_init() result(myPlasticity)
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
-    Nconstituents = count(material_phaseAt2 == ph)
+    Nmembers = count(material_phaseAt2 == ph)
     sizeDotState = size(['xi   ','gamma'])
     sizeState = sizeDotState
 
-    call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0)
+    call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
 
 !--------------------------------------------------------------------------------------------------
 ! state aliases and initialization

src/phase_mechanical_plastic_kinehardening.f90

@@ -62,7 +62,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
   logical, dimension(:), allocatable :: myPlasticity
   integer :: &
     ph, o,  &
-    Nconstituents, &
+    Nmembers, &
     sizeState, sizeDeltaState, sizeDotState, &
     startIndex, endIndex
   integer,     dimension(:), allocatable :: &
@@ -165,19 +165,19 @@ module function plastic_kinehardening_init() result(myPlasticity)
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
-    Nconstituents = count(material_phaseAt2 == ph)
+    Nmembers = count(material_phaseAt2 == ph)
     sizeDotState   = size(['crss     ','crss_back', 'accshear ']) * prm%sum_N_sl !ToDo: adjust names like in material.yaml
     sizeDeltaState = size(['sense ',   'chi0  ',    'gamma0'   ]) * prm%sum_N_sl !ToDo: adjust names like in material.yaml
     sizeState = sizeDotState + sizeDeltaState
 
-    call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,sizeDeltaState)
+    call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState)
 
 !--------------------------------------------------------------------------------------------------
 ! state aliases and initialization
     startIndex = 1
     endIndex   = prm%sum_N_sl
     stt%crss => plasticState(ph)%state   (startIndex:endIndex,:)
-    stt%crss = spread(xi_0, 2, Nconstituents)
+    stt%crss = spread(xi_0, 2, Nmembers)
     dot%crss => plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
     if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'

src/phase_mechanical_plastic_nonlocal.f90

@@ -177,7 +177,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
   integer :: &
     Ninstances, &
     ph, &
-    Nconstituents, &
+    Nmembers, &
     sizeState, sizeDotState, sizeDependentState, sizeDeltaState, &
     s1, s2, &
     s, t, l
@@ -398,7 +398,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
-    Nconstituents  = count(material_phaseAt2 == ph)
+    Nmembers  = count(material_phaseAt2 == ph)
     sizeDotState = size([   'rhoSglEdgePosMobile   ','rhoSglEdgeNegMobile   ', &
                             'rhoSglScrewPosMobile  ','rhoSglScrewNegMobile  ', &
                             'rhoSglEdgePosImmobile ','rhoSglEdgeNegImmobile ', &
@@ -412,9 +412,9 @@ module function plastic_nonlocal_init() result(myPlasticity)
                                 'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl     !< other dependent state variables that are not updated by microstructure
     sizeDeltaState            = sizeDotState
 
-    call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,sizeDeltaState)
+    call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState)
 
-    allocate(geom(ph)%V_0(Nconstituents))
+    allocate(geom(ph)%V_0(Nmembers))
     call storeGeometry(ph)
 
     plasticState(ph)%nonlocal         = pl%get_asBool('nonlocal')
@@ -486,26 +486,26 @@ module function plastic_nonlocal_init() result(myPlasticity)
         dot%rho_dip_scr => plasticState(ph)%dotState               (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
         del%rho_dip_scr => plasticState(ph)%deltaState             (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
 
-    stt%gamma => plasticState(ph)%state                      (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
+    stt%gamma => plasticState(ph)%state                      (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nmembers)
-    dot%gamma => plasticState(ph)%dotState                   (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
+    dot%gamma => plasticState(ph)%dotState                   (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nmembers)
-    del%gamma => plasticState(ph)%deltaState                 (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
+    del%gamma => plasticState(ph)%deltaState                 (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nmembers)
     plasticState(ph)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl )  = pl%get_asFloat('atol_gamma', defaultVal = 1.0e-2_pReal)
     if(any(plasticState(ph)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl) < 0.0_pReal)) &
       extmsg = trim(extmsg)//' atol_gamma'
-    plasticState(ph)%slipRate => plasticState(ph)%dotState    (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
+    plasticState(ph)%slipRate => plasticState(ph)%dotState    (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nmembers)
 
-    stt%rho_forest => plasticState(ph)%state                 (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:Nconstituents)
+    stt%rho_forest => plasticState(ph)%state                 (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:Nmembers)
-    stt%v          => plasticState(ph)%state                 (12*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nconstituents)
+    stt%v          => plasticState(ph)%state                 (12*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nmembers)
-        stt%v_edg_pos  => plasticState(ph)%state             (12*prm%sum_N_sl + 1:13*prm%sum_N_sl,1:Nconstituents)
+        stt%v_edg_pos  => plasticState(ph)%state             (12*prm%sum_N_sl + 1:13*prm%sum_N_sl,1:Nmembers)
-        stt%v_edg_neg  => plasticState(ph)%state             (13*prm%sum_N_sl + 1:14*prm%sum_N_sl,1:Nconstituents)
+        stt%v_edg_neg  => plasticState(ph)%state             (13*prm%sum_N_sl + 1:14*prm%sum_N_sl,1:Nmembers)
-        stt%v_scr_pos  => plasticState(ph)%state             (14*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:Nconstituents)
+        stt%v_scr_pos  => plasticState(ph)%state             (14*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:Nmembers)
-        stt%v_scr_neg  => plasticState(ph)%state             (15*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nconstituents)
+        stt%v_scr_neg  => plasticState(ph)%state             (15*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nmembers)
 
-    allocate(dst%tau_pass(prm%sum_N_sl,Nconstituents),source=0.0_pReal)
+    allocate(dst%tau_pass(prm%sum_N_sl,Nmembers),source=0.0_pReal)
-    allocate(dst%tau_back(prm%sum_N_sl,Nconstituents),source=0.0_pReal)
+    allocate(dst%tau_back(prm%sum_N_sl,Nmembers),source=0.0_pReal)
     end associate
 
-    if (Nconstituents > 0) call stateInit(ini,ph,Nconstituents)
+    if (Nmembers > 0) call stateInit(ini,ph,Nmembers)
     plasticState(ph)%state0 = plasticState(ph)%state
 
 !--------------------------------------------------------------------------------------------------
@@ -527,7 +527,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
     if(.not. myPlasticity(ph)) cycle
 
     phase => phases%get(ph)
-    Nconstituents = count(material_phaseAt2 == ph)
+    Nmembers = count(material_phaseAt2 == ph)
     l = 0
     do t = 1,4
       do s = 1,param(ph)%sum_N_sl
@@ -1403,8 +1403,10 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
 
   integer :: &
     n, &                                                                                            ! neighbor index
+    me, &
     neighbor_e, &                                                                                   ! element index of my neighbor
     neighbor_i, &                                                                                   ! integration point index of my neighbor
+    neighbor_me, &
     neighbor_phase, &
     ns, &                                                                                           ! number of active slip systems
     s1, &                                                                                           ! slip system index (me)
@@ -1422,6 +1424,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
   associate(prm => param(ph))
   ns = prm%sum_N_sl
 
+  me = material_phaseMemberAt(1,i,e)
   !*** start out fully compatible
   my_compatibility = 0.0_pReal
   forall(s1 = 1:ns) my_compatibility(:,s1,s1,:) = 1.0_pReal
@@ -1429,7 +1432,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
   neighbors: do n = 1,nIPneighbors
     neighbor_e = IPneighborhood(1,n,i,e)
     neighbor_i = IPneighborhood(2,n,i,e)
-
+    neighbor_me = material_phaseMemberAt(1,neighbor_i,neighbor_e)
     neighbor_phase = material_phaseAt(1,neighbor_e)
 
     if (neighbor_e <= 0 .or. neighbor_i <= 0) then
@@ -1447,8 +1450,8 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
     elseif (prm%chi_GB >= 0.0_pReal) then
       !* GRAIN BOUNDARY !
       !* fixed transmissivity for adjacent ips with different texture (only if explicitly given in material.config)
-      if (any(dNeq(material_orientation0(1,i,e)%asQuaternion(), &
+      if (any(dNeq(material_orientation0(1,ph,me)%asQuaternion(), &
-                   material_orientation0(1,neighbor_i,neighbor_e)%asQuaternion())) .and. &
+                   material_orientation0(1,neighbor_phase,neighbor_me)%asQuaternion())) .and. &
           (.not. phase_localPlasticity(neighbor_phase))) &
         forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB)
     else
@@ -1576,13 +1579,13 @@ end subroutine plastic_nonlocal_results
 !--------------------------------------------------------------------------------------------------
 !> @brief populates the initial dislocation density
 !--------------------------------------------------------------------------------------------------
-subroutine stateInit(ini,phase,Nconstituents)
+subroutine stateInit(ini,phase,Nmembers)
 
   type(tInitialParameters) :: &
     ini
   integer,intent(in) :: &
     phase, &
-    Nconstituents
+    Nmembers
   integer :: &
     i, &
     e, &
@@ -1599,7 +1602,7 @@ subroutine stateInit(ini,phase,Nconstituents)
     totalVolume, &
     densityBinning, &
     minimumIpVolume
-  real(pReal), dimension(Nconstituents) :: &
+  real(pReal), dimension(Nmembers) :: &
     volume
 
 
@@ -1619,13 +1622,13 @@ subroutine stateInit(ini,phase,Nconstituents)
     meanDensity = 0.0_pReal
     do while(meanDensity < ini%random_rho_u)
       call random_number(rnd)
-      phasemember = nint(rnd(1)*real(Nconstituents,pReal) + 0.5_pReal)
+      phasemember = nint(rnd(1)*real(Nmembers,pReal) + 0.5_pReal)
       s           = nint(rnd(2)*real(sum(ini%N_sl),pReal)*4.0_pReal + 0.5_pReal)
       meanDensity = meanDensity + densityBinning * volume(phasemember) / totalVolume
       stt%rhoSglMobile(s,phasemember) = densityBinning
     enddo
   else                                ! homogeneous distribution with noise
-    do e = 1, Nconstituents
+    do e = 1, Nmembers
       do f = 1,size(ini%N_sl,1)
         from = 1 + sum(ini%N_sl(1:f-1))
         upto = sum(ini%N_sl(1:f))
@@ -1806,7 +1809,7 @@ pure function getRho0(ph,me)
     getRho0(:,mob) = max(getRho0(:,mob),0.0_pReal)
     getRho0(:,dip) = max(getRho0(:,dip),0.0_pReal)
 
-    where(abs(getRho0) < max(prm%rho_min/geom(ph)%V_0(me)**(2.0_pReal/3.0_pReal),prm%rho_significant)) &
+    where (abs(getRho0) < max(prm%rho_min/geom(ph)%V_0(me)**(2.0_pReal/3.0_pReal),prm%rho_significant)) &
       getRho0 = 0.0_pReal
 
   end associate
@@ -1819,16 +1822,13 @@ subroutine storeGeometry(ph)
   integer, intent(in) :: ph
 
   integer :: ip, el, ce, co
+  real(pReal), dimension(:), allocatable :: V
 
-  ce = 0
+
-  do el = 1, size(material_homogenizationMemberAt,2)
+  V = reshape(IPvolume,[product(shape(IPvolume))])
-    do ip = 1, size(material_homogenizationMemberAt,1)
+  do ce = 1, size(material_homogenizationMemberAt2,1)
-      ce = ce + 1
+    do co = 1, homogenization_maxNconstituents
-      do co = 1, homogenization_maxNconstituents
+      if (material_phaseAt2(co,ce) == ph) geom(ph)%V_0(material_phaseMemberAt2(co,ce)) = V(ce)
-        if(material_phaseAt2(co,ce) == ph) then
-           geom(ph)%V_0(material_phaseMemberAt2(co,ce)) = IPvolume(ip,el)
-        endif
-      enddo
     enddo
   enddo
 

src/phase_mechanical_plastic_phenopowerlaw.f90

@@ -71,7 +71,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
   logical, dimension(:), allocatable :: myPlasticity
   integer :: &
     ph, i, &
-    Nconstituents, &
+    Nmembers, &
     sizeState, sizeDotState, &
     startIndex, endIndex
   integer,     dimension(:), allocatable :: &
@@ -223,20 +223,20 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
-    Nconstituents = count(material_phaseAt2 == ph)
+    Nmembers = count(material_phaseAt2 == ph)
     sizeDotState = size(['xi_sl   ','gamma_sl']) * prm%sum_N_sl &
                  + size(['xi_tw   ','gamma_tw']) * prm%sum_N_tw
     sizeState = sizeDotState
 
 
-    call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0)
+    call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
 
 !--------------------------------------------------------------------------------------------------
 ! state aliases and initialization
     startIndex = 1
     endIndex   = prm%sum_N_sl
     stt%xi_slip => plasticState(ph)%state   (startIndex:endIndex,:)
-    stt%xi_slip =  spread(xi_0_sl, 2, Nconstituents)
+    stt%xi_slip =  spread(xi_0_sl, 2, Nmembers)
     dot%xi_slip => plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
     if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
@@ -244,7 +244,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
     startIndex = endIndex + 1
     endIndex   = endIndex + prm%sum_N_tw
     stt%xi_twin => plasticState(ph)%state   (startIndex:endIndex,:)
-    stt%xi_twin =  spread(xi_0_tw, 2, Nconstituents)
+    stt%xi_twin =  spread(xi_0_tw, 2, Nmembers)
     dot%xi_twin => plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
     if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'

src/phase_thermal.f90

@@ -15,13 +15,13 @@ submodule(phase) thermal
     THERMAL_EXTERNALHEAT_ID
   end enum
 
-  type :: tDataContainer
+  type :: tDataContainer             ! ?? not very telling name. Better: "fieldQuantities" ??
     real(pReal), dimension(:), allocatable :: T, dot_T
   end type tDataContainer
   integer(kind(THERMAL_UNDEFINED_ID)),  dimension(:,:), allocatable :: &
     thermal_source
 
-  type(tDataContainer), dimension(:), allocatable :: current
+  type(tDataContainer), dimension(:), allocatable :: current          ! ?? not very telling name. Better: "field" ??
 
   integer :: thermal_source_maxSizeDotState
 
@@ -78,41 +78,36 @@ module subroutine thermal_init(phases)
 
   integer :: &
     ph, so, &
-    Nconstituents
+    Nmembers
 
 
   print'(/,a)', ' <<<+-  phase:thermal init  -+>>>'
 
   allocate(current(phases%length))
 
-  allocate(thermalState (phases%length))
+  allocate(thermalState(phases%length))
   allocate(thermal_Nsources(phases%length),source = 0)
 
   do ph = 1, phases%length
-
+    Nmembers = count(material_phaseAt2 == ph)
-    Nconstituents = count(material_phaseAt2 == ph)
+    allocate(current(ph)%T(Nmembers),source=300.0_pReal)
-
+    allocate(current(ph)%dot_T(Nmembers),source=0.0_pReal)
-    allocate(current(ph)%T(Nconstituents),source=300.0_pReal)
-    allocate(current(ph)%dot_T(Nconstituents),source=0.0_pReal)
     phase => phases%get(ph)
-    if(phase%contains('thermal')) then
+    thermal => phase%get('thermal',defaultVal=emptyDict)
-      thermal => phase%get('thermal')
+    sources => thermal%get('source',defaultVal=emptyList)
-      sources => thermal%get('source',defaultVal=emptyList)
+    thermal_Nsources(ph) = sources%length
-
-      thermal_Nsources(ph) = sources%length
-    endif
     allocate(thermalstate(ph)%p(thermal_Nsources(ph)))
   enddo
 
   allocate(thermal_source(maxval(thermal_Nsources),phases%length), source = THERMAL_UNDEFINED_ID)
 
-  if(maxval(thermal_Nsources) /= 0) then
+  if (maxval(thermal_Nsources) /= 0) then
     where(dissipation_init (maxval(thermal_Nsources))) thermal_source = THERMAL_DISSIPATION_ID
     where(externalheat_init(maxval(thermal_Nsources))) thermal_source = THERMAL_EXTERNALHEAT_ID
   endif
 
   thermal_source_maxSizeDotState = 0
-  PhaseLoop2:do ph = 1,phases%length
+  do ph = 1,phases%length
 
     do so = 1,thermal_Nsources(ph)
       thermalState(ph)%p(so)%state  = thermalState(ph)%p(so)%state0
@@ -120,7 +115,7 @@ module subroutine thermal_init(phases)
 
     thermal_source_maxSizeDotState  = max(thermal_source_maxSizeDotState, &
                                           maxval(thermalState(ph)%p%sizeDotState))
-  enddo PhaseLoop2
+  enddo
 
 end subroutine thermal_init
 
@@ -145,18 +140,17 @@ module subroutine phase_thermal_getRate(TDot, ph,me)
   do so = 1, thermal_Nsources(ph)
    select case(thermal_source(so,ph))
      case (THERMAL_DISSIPATION_ID)
-      call dissipation_getRate(my_Tdot, ph,me)
+       call dissipation_getRate(my_Tdot, ph,me)
 
      case (THERMAL_EXTERNALHEAT_ID)
-      call externalheat_getRate(my_Tdot, ph,me)
+       call externalheat_getRate(my_Tdot, ph,me)
 
      case default
-      my_Tdot = 0.0_pReal
+       my_Tdot = 0.0_pReal
    end select
    Tdot = Tdot + my_Tdot
   enddo
 
-
 end subroutine phase_thermal_getRate
 
 
@@ -185,7 +179,7 @@ function phase_thermal_collectDotState(ph,me) result(broken)
 end function phase_thermal_collectDotState
 
 
-module function thermal_stress(Delta_t,ph,me) result(converged_)
+module function thermal_stress(Delta_t,ph,me) result(converged_)           ! ?? why is this called "stress" when it seems closer to "updateState" ??
 
   real(pReal), intent(in) :: Delta_t
   integer, intent(in) :: ph, me
@@ -212,7 +206,7 @@ function integrateThermalState(Delta_t, ph,me) result(broken)
     sizeDotState
 
   broken = phase_thermal_collectDotState(ph,me)
-  if(broken) return
+  if (broken) return
 
   do so = 1, thermal_Nsources(ph)
     sizeDotState = thermalState(ph)%p(so)%sizeDotState
@@ -301,14 +295,12 @@ function thermal_active(source_label,src_length)  result(active_source)
   allocate(active_source(src_length,phases%length), source = .false. )
   do p = 1, phases%length
     phase => phases%get(p)
-    if (phase%contains('thermal')) then
+    thermal =>  phase%get('thermal',defaultVal=emptyDict)
-      thermal =>  phase%get('thermal',defaultVal=emptyList)
+    sources =>  thermal%get('source',defaultVal=emptyList)
-      sources =>  thermal%get('source',defaultVal=emptyList)
+    do s = 1, sources%length
-      do s = 1, sources%length
+      src => sources%get(s)
-        src => sources%get(s)
+      active_source(s,p) = src%get_asString('type') == source_label
-        if(src%get_asString('type') == source_label) active_source(s,p) = .true.
+    enddo
-      enddo
-    endif
   enddo
 
 

src/phase_thermal_dissipation.f90

@@ -31,7 +31,7 @@ module function dissipation_init(source_length) result(mySources)
     phase, &
     sources, thermal, &
     src
-  integer :: so,Nconstituents,ph
+  integer :: so,Nmembers,ph
 
 
   mySources = thermal_active('dissipation',source_length)
@@ -54,8 +54,8 @@ module function dissipation_init(source_length) result(mySources)
           src => sources%get(so)
 
           prm%kappa = src%get_asFloat('kappa')
-          Nconstituents = count(material_phaseAt2 == ph)
+          Nmembers = count(material_phaseAt2 == ph)
-          call phase_allocateState(thermalState(ph)%p(so),Nconstituents,0,0,0)
+          call phase_allocateState(thermalState(ph)%p(so),Nmembers,0,0,0)
 
         end associate
       endif

src/phase_thermal_externalheat.f90

@@ -38,7 +38,7 @@ module function externalheat_init(source_length) result(mySources)
     phase, &
     sources, thermal, &
     src
-  integer :: so,Nconstituents,ph
+  integer :: so,Nmembers,ph
 
 
   mySources = thermal_active('externalheat',source_length)
@@ -67,8 +67,8 @@ module function externalheat_init(source_length) result(mySources)
 
           prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n))
 
-          Nconstituents = count(material_phaseAt2 == ph)
+          Nmembers = count(material_phaseAt2 == ph)
-          call phase_allocateState(thermalState(ph)%p(so),Nconstituents,1,1,0)
+          call phase_allocateState(thermalState(ph)%p(so),Nmembers,1,1,0)
         end associate
       endif
     enddo