Merge branch 'development' into python-style

.gitattributes

@@ -8,3 +8,10 @@
 *.jpg binary
 *.hdf5 binary
 *.pdf binary
+
+# ignore files from MSC.Marc in language statistics
+installation/mods_MarcMentat/* linguist-vendored
+src/MarcInclude/*              linguist-vendored
+
+# ignore reference files for tests in language statistics
+python/tests/reference/*       linguist-vendored

.gitlab-ci.yml

@@ -104,7 +104,7 @@ checkout:
     - release
 
 ###################################################################################################
-Pytest:
+Pytest_python:
   stage: python
   script:
     - cd $DAMASKROOT/python
@@ -142,13 +142,6 @@ Pre_General:
     - master
     - release
 
-grid_geometryPacking:
-  stage: preprocessing
-  script: grid_geometryPacking/test.py
-  except:
-    - master
-    - release
-
 ###################################################################################################
 Post_AverageDown:
   stage: postprocessing
@@ -385,6 +378,15 @@ Phenopowerlaw_singleSlip:
     - master
     - release
 
+Pytest_grid:
+  stage: grid
+  script:
+    - cd pytest
+    - pytest
+  except:
+    - master
+    - release
+
 ###################################################################################################
 Marc_compileIfort:
   stage: compileMarc

PRIVATE

@@ -1 +1 @@
-Subproject commit ec615d249d39e5d01446b01ab9a5b7e7601340ad
+Subproject commit 9573ce7bd2c1a7188c1aac5b83aa76d480c2bdb0

VERSION

@@ -1 +1 @@
-v2.0.3-1726-gef4b7437
+v2.0.3-1808-g13c3ce00

examples/ConfigFiles/Phase_Phenopowerlaw_multiField.config

@@ -1,49 +0,0 @@
-[Aluminum]
-elasticity              hooke
-plasticity              phenopowerlaw
-
-(output)                resistance_slip
-(output)                shearrate_slip
-(output)                resolvedstress_slip
-(output)                accumulated_shear_slip
-(output)                resistance_twin
-(output)                shearrate_twin
-(output)                resolvedstress_twin
-(output)                accumulated_shear_twin
-
-lattice_structure       fcc
-Nslip                   12                # per family
-Ntwin                   0                # per family
-
-c11                     106.75e9
-c12                     60.41e9
-c44                     28.34e9
-
-gdot0_slip              0.001
-n_slip                  20
-tau0_slip               31e6                # per family
-tausat_slip             63e6                # per family
-a_slip                  2.25
-gdot0_twin              0.001
-n_twin                  20
-tau0_twin               31e6                # per family
-h0_slipslip             75e6
-interaction_slipslip    1 1 1.4 1.4 1.4 1.4
-atol_resistance         1
-
-(stiffness_degradation) damage
-(stiffness_degradation) porosity
-{./Phase_Damage.config}
-{./Phase_Thermal.config}
-{./Phase_Vacancy.config}
-{./Phase_Porosity.config}
-{./Phase_Hydrogen.config}
-{./Source_Damage_IsoBrittle.config}
-{./Source_Thermal_Dissipation.config}
-{./Source_Vacancy_PhenoPlasticity.config}
-{./Source_Vacancy_Irradiation.config}
-{./Kinematics_Thermal_Expansion.config}
-{./Kinematics_Vacancy_Strain.config}
-{./Kinematics_Hydrogen_Strain.config}
-
-

processing/post/DADF5_postResults.py

@@ -52,15 +52,15 @@ for filename in options.filenames:
     table = damask.Table(np.ones(np.product(results.grid),dtype=int)*int(inc[3:]),{'inc':(1,)})
     table.add('pos',coords.reshape((-1,3)))
 
-    results.set_visible('materialpoints',False)
-    results.set_visible('constituents',  True)
+    results.pick('materialpoints',False)
+    results.pick('constituents',  True)
     for label in options.con:
       x = results.get_dataset_location(label)
       if len(x) != 0:
         table.add(label,results.read_dataset(x,0,plain=True).reshape((results.grid.prod(),-1)))
 
-    results.set_visible('constituents',  False)
-    results.set_visible('materialpoints',True)
+    results.pick('constituents',  False)
+    results.pick('materialpoints',True)
     for label in options.mat:
       x = results.get_dataset_location(label)
       if len(x) != 0:

processing/post/addCauchy.py

@@ -42,8 +42,8 @@ for name in filenames:
 
     table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
     table.add('Cauchy',
-              damask.mechanics.Cauchy(table.get(options.defgrad).reshape(-1,3,3),
-                                      table.get(options.stress ).reshape(-1,3,3)).reshape(-1,9),
+              damask.mechanics.Cauchy(table.get(options.stress ).reshape(-1,3,3),
+                                      table.get(options.defgrad).reshape(-1,3,3)).reshape(-1,9),
               scriptID+' '+' '.join(sys.argv[1:]))
 
     table.to_ASCII(sys.stdout if name is None else name)

processing/post/addPK2.py

@@ -43,8 +43,8 @@ for name in filenames:
     table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
 
     table.add('S',
-              damask.mechanics.PK2(table.get(options.defgrad).reshape(-1,3,3),
-                                   table.get(options.stress ).reshape(-1,3,3)).reshape(-1,9),
+              damask.mechanics.PK2(table.get(options.stress ).reshape(-1,3,3),
+                                   table.get(options.defgrad).reshape(-1,3,3)).reshape(-1,9),
               scriptID+' '+' '.join(sys.argv[1:]))
 
     table.to_ASCII(sys.stdout if name is None else name)

processing/post/addSpectralDecomposition.py

@@ -2,6 +2,7 @@
 
 import os
 import sys
+from io import StringIO
 from optparse import OptionParser
 
 import numpy as np
@@ -33,69 +34,27 @@ parser.add_option('--no-check',
 
 parser.set_defaults(rh = True,
                     )
+
 (options,filenames) = parser.parse_args()
-
-if options.tensor is None:
-  parser.error('no data column specified.')
-
-# --- loop over input files -------------------------------------------------------------------------
-
 if filenames == []: filenames = [None]
 
 for name in filenames:
-  try:
-    table = damask.ASCIItable(name = name,
-                              buffered = False)
-  except: continue
-  damask.util.report(scriptName,name)
+    damask.util.report(scriptName,name)
 
-# ------------------------------------------ read header ------------------------------------------
+    table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
 
-  table.head_read()
+    for tensor in options.tensor:
 
+        t = table.get(tensor).reshape(-1,3,3)
+        (u,v) = np.linalg.eigh(damask.mechanics.symmetric(t))
+        if options.rh: v[np.linalg.det(v) < 0.0,:,2] *= -1.0
 
-# ------------------------------------------ assemble header 1 ------------------------------------
+        for i,o in enumerate(['Min','Mid','Max']):
+            table.add('eigval{}({})'.format(o,tensor),u[:,i],
+                      scriptID+' '+' '.join(sys.argv[1:]))
 
-  items = {
-            'tensor': {'dim': 9, 'shape': [3,3], 'labels':options.tensor, 'column': []},
-          }
-  errors  = []
-  remarks = []
+        for i,o in enumerate(['Min','Mid','Max']):
+            table.add('eigvec{}({})'.format(o,tensor),v[:,:,i],
+                      scriptID+' '+' '.join(sys.argv[1:]))
         
-  for type, data in items.items():
-    for what in data['labels']:
-      dim = table.label_dimension(what)
-      if dim != data['dim']: remarks.append('column {} is not a {}...'.format(what,type))
-      else:
-        items[type]['column'].append(table.label_index(what))
-        for order in ['Min','Mid','Max']:
-          table.labels_append(['eigval{}({})'.format(order,what)])                                         # extend ASCII header with new labels
-        for order in ['Min','Mid','Max']:
-          table.labels_append(['{}_eigvec{}({})'.format(i+1,order,what) for i in range(3)])                # extend ASCII header with new labels
-
-  if remarks != []: damask.util.croak(remarks)
-  if errors  != []:
-    damask.util.croak(errors)
-    table.close(dismiss = True)
-    continue
-
-# ------------------------------------------ assemble header 2 ------------------------------------
-
-  table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
-  table.head_write()
-
-# ------------------------------------------ process data -----------------------------------------
-
-  outputAlive = True
-  while outputAlive and table.data_read():                                                          # read next data line of ASCII table
-    for type, data in items.items():
-      for column in data['column']:
-        (u,v) = np.linalg.eigh(np.array(list(map(float,table.data[column:column+data['dim']]))).reshape(data['shape']))
-        if options.rh and np.dot(np.cross(v[:,0], v[:,1]), v[:,2]) < 0.0 : v[:, 2] *= -1.0          # ensure right-handed eigenvector basis
-        table.data_append(list(u))                                                                  # vector of max,mid,min eigval
-        table.data_append(list(v.transpose().reshape(data['dim'])))                                 # 3x3=9 combo vector of max,mid,min eigvec coordinates
-    outputAlive = table.data_write()                                                                # output processed line in accordance with column labeling
-
-# ------------------------------------------ output finalization -----------------------------------
-
-  table.close()                                                                                     # close input ASCII table (works for stdin)
+    table.to_ASCII(sys.stdout if name is None else name)

processing/post/permuteData.py

@@ -83,7 +83,7 @@ for name in filenames:
        
 # ------------------------------------------ assemble header ---------------------------------------
 
-  randomSeed = int(os.urandom(4).encode('hex'), 16) if options.randomSeed is None else options.randomSeed         # random seed per file
+  randomSeed = int(os.urandom(4).hex(), 16) if options.randomSeed is None else options.randomSeed   # random seed per file
   np.random.seed(randomSeed)
 
   table.info_append([scriptID + '\t' + ' '.join(sys.argv[1:]),

processing/pre/seeds_fromDistribution.py

@@ -1,11 +1,9 @@
 #!/usr/bin/env python3
-# -*- coding: UTF-8 no BOM -*-
 
 import threading,time,os,sys,random
 import numpy as np
 from optparse import OptionParser
 from io import StringIO
-import binascii
 import damask
 
 scriptName = os.path.splitext(os.path.basename(__file__))[0]
@@ -16,9 +14,10 @@ currentSeedsName = None
 
 #---------------------------------------------------------------------------------------------------
 class myThread (threading.Thread):
-  """perturbes seed in seed file, performes Voronoi tessellation, evaluates, and updates best match"""
+  """Perturb seed in seed file, performes Voronoi tessellation, evaluates, and updates best match."""
 
   def __init__(self, threadID):
+    """Threading class with thread ID."""
     threading.Thread.__init__(self)
     self.threadID = threadID
 
@@ -215,20 +214,16 @@ options = parser.parse_args()[0]
 damask.util.report(scriptName,options.seedFile)
 
 if options.randomSeed is None:
-  options.randomSeed = int(binascii.hexlify(os.urandom(4)),16)
+  options.randomSeed = int(os.urandom(4).hex(),16)
 damask.util.croak(options.randomSeed)
 delta = (options.scale/options.grid[0],options.scale/options.grid[1],options.scale/options.grid[2])
 baseFile=os.path.splitext(os.path.basename(options.seedFile))[0]
 points = np.array(options.grid).prod().astype('float')
 
 # ----------- calculate target distribution and bin edges
-targetGeomFile = os.path.splitext(os.path.basename(options.target))[0]+'.geom'
-targetGeomTable = damask.ASCIItable(targetGeomFile,None,labeled=False,readonly=True)
-targetGeomTable.head_read()
-info,devNull =  targetGeomTable.head_getGeom()
-nMicrostructures = info['microstructures']
-targetVolFrac = np.bincount(targetGeomTable.microstructure_read(info['grid']))[1:nMicrostructures+1]/\
-                                                                       float(info['grid'].prod())
+targetGeom = damask.Geom.from_file(os.path.splitext(os.path.basename(options.target))[0]+'.geom')
+nMicrostructures = len(np.unique(targetGeom.microstructure))
+targetVolFrac = np.bincount(targetGeom.microstructure.flatten())/targetGeom.grid.prod().astype(np.float)
 target=[]
 for i in range(1,nMicrostructures+1):
   targetHist,targetBins = np.histogram(targetVolFrac,bins=i) #bin boundaries
@@ -252,13 +247,12 @@ initialGeomVFile = StringIO()
 initialGeomVFile.write(damask.util.execute('geom_fromVoronoiTessellation '+
                                ' -g '+' '.join(list(map(str, options.grid))),bestSeedsVFile)[0])
 initialGeomVFile.seek(0)
-initialGeomTable = damask.ASCIItable(initialGeomVFile,None,labeled=False,readonly=True)
-initialGeomTable.head_read()
-info,devNull =  initialGeomTable.head_getGeom()
+initialGeom = damask.Geom.from_file(initialGeomVFile)
 
-if info['microstructures'] != nMicrostructures: damask.util.croak('error. Microstructure count mismatch')
+if len(np.unique(targetGeom.microstructure)) != nMicrostructures:
+  damask.util.croak('error. Microstructure count mismatch')
 
-initialData = np.bincount(initialGeomTable.microstructure_read(info['grid']))/points
+initialData = np.bincount(initialGeom.microstructure.flatten())/points
 for i in range(nMicrostructures):
   initialHist = np.histogram(initialData,bins=target[i]['bins'])[0]
   target[i]['error']=np.sqrt(np.square(np.array(target[i]['histogram']-initialHist)).sum())
@@ -274,7 +268,7 @@ for i in range(nMicrostructures):
 
 
 if options.maxseeds < 1: 
-  maxSeeds = info['microstructures']
+  maxSeeds = len(np.unique(initialGeom.microstructure))
 else:
   maxSeeds = options.maxseeds
 

python/README

@@ -1,13 +0,0 @@
-DAMASK - The Düsseldorf Advanced Material Simulation Kit
-Visit damask.mpie.de for installation and usage instructions
-
-CONTACT INFORMATION
-
-Max-Planck-Institut für Eisenforschung GmbH
-Max-Planck-Str. 1
-40237 Düsseldorf
-Germany
-
-Email: DAMASK@mpie.de
-https://damask.mpie.de
-https://magit1.mpie.de

python/README

@@ -0,0 +1 @@
+../README

python/damask/__init__.py

@@ -14,9 +14,10 @@ from .asciitable  import ASCIItable       # noqa
 from .config      import Material         # noqa
 from .colormaps   import Colormap, Color  # noqa
 from .rotation    import Rotation         # noqa
-from .lattice     import Symmetry, Lattice # noqa
+from .lattice     import Symmetry, Lattice# noqa
 from .orientation import Orientation      # noqa
-from .dadf5       import DADF5            # noqa
+from .result      import Result           # noqa
+from .result      import Result as DADF5  # noqa
 
 from .geom        import Geom             # noqa
 from .solver      import Solver           # noqa

python/damask/asciitable.py

@@ -16,7 +16,7 @@ class ASCIItable():
   def __init__(self,
                name    = None,
                outname = None,
-               buffered  = False,                                                                   # flush writes
+               buffered  = False,                                                                   # is ignored, only exists for compatibility reasons
                labeled   = True,                                                                    # assume table has labels
                readonly  = False,                                                                   # no reading from file
               ):
@@ -63,7 +63,6 @@ class ASCIItable():
     except AttributeError:
       return str(string)
 
-
 # ------------------------------------------------------------------
   def _quote(self,
              what):
@@ -71,6 +70,7 @@ class ASCIItable():
     return '{quote}{content}{quote}'.format(
              quote   = ('"' if str(what)=='' or re.search(r"\s",str(what)) else ''),
              content = what)
+
 # ------------------------------------------------------------------
   def close(self,
             dismiss = False):
@@ -178,15 +178,11 @@ class ASCIItable():
             'grid':            lambda x: int(x),
             'size':            lambda x: float(x),
             'origin':          lambda x: float(x),
-            'homogenization':  lambda x: int(x),
-            'microstructures': lambda x: int(x),
               }
     info = {
             'grid':            np.zeros(3,'i'),
             'size':            np.zeros(3,'d'),
             'origin':          np.zeros(3,'d'),
-            'homogenization':  0,
-            'microstructures': 0,
            }
     extra_header = []
 
@@ -375,15 +371,6 @@ class ASCIItable():
     self.tags = list(self.__IO__['tags'])                                                           # restore label info found in header (as COPY, not link)
     self.__IO__['labeled'] = len(self.tags) > 0
 
-# ------------------------------------------------------------------
-  def data_skipLines(self,
-                     count):
-    """Wind forward by count number of lines."""
-    for i in range(count):
-      alive = self.data_read()
-
-    return alive
-
 # ------------------------------------------------------------------
   def data_read(self,
                 advance = True,
@@ -473,33 +460,3 @@ class ASCIItable():
         for item in what: self.data_append(item)
       except TypeError:
         self.data += [str(what)]
-
-# ------------------------------------------------------------------
-  def microstructure_read(self,
-                          grid,
-                          type = 'i',
-                          strict = False):
-    """Read microstructure data (from .geom format)."""
-    def datatype(item):
-      return int(item) if type.lower() == 'i' else float(item)
-      
-    N = grid.prod()                                                                          # expected number of microstructure indices in data
-    microstructure = np.zeros(N,type)                                                        # initialize as flat array
-
-    i = 0
-    while i < N and self.data_read():
-      items = self.data
-      if len(items) > 2:
-        if   items[1].lower() == 'of':
-          items = np.ones(datatype(items[0]))*datatype(items[2])
-        elif items[1].lower() == 'to':
-          items = np.linspace(datatype(items[0]),datatype(items[2]),
-                              abs(datatype(items[2])-datatype(items[0]))+1,dtype=int)
-        else:                          items = list(map(datatype,items))
-      else:                            items = list(map(datatype,items))
-
-      s = min(len(items), N-i)                                                              # prevent overflow of microstructure array
-      microstructure[i:i+s] = items[:s]
-      i += len(items)
-
-    return (microstructure, i == N and not self.data_read()) if strict else microstructure  # check for proper point count and end of file

python/damask/environment.py

@@ -8,7 +8,7 @@ class Environment():
     def __init__(self):
         """Read and provide values of DAMASK configuration."""
         self.options = {}
-        self.get_options()
+        self.__get_options()
 
     def relPath(self,relative = '.'):
         return os.path.join(self.rootDir(),relative)
@@ -16,7 +16,7 @@ class Environment():
     def rootDir(self):
         return os.path.normpath(os.path.join(os.path.realpath(__file__),'../../../'))
 
-    def get_options(self):
+    def __get_options(self):
         for item in ['DAMASK_NUM_THREADS',
                      'MSC_ROOT',
                      'MARC_VERSION',

python/damask/mechanics.py

@@ -1,8 +1,8 @@
 import numpy as np
 
-def Cauchy(F,P):
+def Cauchy(P,F):
     """
-    Return Cauchy stress calculated from 1. Piola-Kirchhoff stress and deformation gradient.
+    Return Cauchy stress calculated from first Piola-Kirchhoff stress and deformation gradient.
 
     Resulting tensor is symmetrized as the Cauchy stress needs to be symmetric.
 
@@ -21,23 +21,180 @@ def Cauchy(F,P):
     return symmetric(sigma)
 
 
-def PK2(F,P):
+def deviatoric_part(T):
     """
-    Return 2. Piola-Kirchhoff stress calculated from 1. Piola-Kirchhoff stress and deformation gradient.
+    Return deviatoric part of a tensor.
+
+    Parameters
+    ----------
+    T : numpy.array of shape (:,3,3) or (3,3)
+      Tensor of which the deviatoric part is computed.
+
+    """
+    return T - np.eye(3)*spherical_part(T) if np.shape(T) == (3,3) else \
+           T - np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),[T.shape[0],3,3]),spherical_part(T))
+
+
+def eigenvalues(T_sym):
+    """
+    Return the eigenvalues, i.e. principal components, of a symmetric tensor.
+
+    The eigenvalues are sorted in ascending order, each repeated according to
+    its multiplicity.
+
+    Parameters
+    ----------
+    T_sym : numpy.array of shape (:,3,3) or (3,3)
+      Symmetric tensor of which the eigenvalues are computed.
+
+    """
+    return np.linalg.eigvalsh(symmetric(T_sym))
+
+
+def eigenvectors(T_sym,RHS=False):
+    """
+    Return eigenvectors of a symmetric tensor.
+
+    The eigenvalues are sorted in ascending order of their associated eigenvalues.
+
+    Parameters
+    ----------
+    T_sym : numpy.array of shape (:,3,3) or (3,3)
+      Symmetric tensor of which the eigenvectors are computed.
+    RHS: bool, optional
+      Enforce right-handed coordinate system. Default is False.
+
+    """
+    (u,v) = np.linalg.eigh(symmetric(T_sym))
+
+    if RHS:
+        if np.shape(T_sym) == (3,3):
+            if np.linalg.det(v) < 0.0: v[:,2] *= -1.0
+        else:
+            v[np.linalg.det(v) < 0.0,:,2] *= -1.0
+    return v
+
+
+def left_stretch(T):
+    """
+    Return the left stretch of a tensor.
+
+    Parameters
+    ----------
+    T : numpy.array of shape (:,3,3) or (3,3)
+      Tensor of which the left stretch is computed.
+
+    """
+    return __polar_decomposition(T,'V')[0]
+
+
+def maximum_shear(T_sym):
+    """
+    Return the maximum shear component of a symmetric tensor.
+
+    Parameters
+    ----------
+    T_sym : numpy.array of shape (:,3,3) or (3,3)
+      Symmetric tensor of which the maximum shear is computed.
+
+    """
+    w = eigenvalues(T_sym)
+    return (w[0] - w[2])*0.5 if np.shape(T_sym) == (3,3) else \
+           (w[:,0] - w[:,2])*0.5
+
+
+def Mises_strain(epsilon):
+    """
+    Return the Mises equivalent of a strain tensor.
+
+    Parameters
+    ----------
+    epsilon : numpy.array of shape (:,3,3) or (3,3)
+      Symmetric strain tensor of which the von Mises equivalent is computed.
+
+    """
+    return __Mises(epsilon,2.0/3.0)
+
+
+def Mises_stress(sigma):
+    """
+    Return the Mises equivalent of a stress tensor.
+
+    Parameters
+    ----------
+    sigma : numpy.array of shape (:,3,3) or (3,3)
+      Symmetric stress tensor of which the von Mises equivalent is computed.
+
+    """
+    return __Mises(sigma,3.0/2.0)
+
+
+def PK2(P,F):
+    """
+    Calculate second Piola-Kirchhoff stress from first Piola-Kirchhoff stress and deformation gradient.
 
     Parameters
     ----------
-    F : numpy.array of shape (:,3,3) or (3,3)
-      Deformation gradient.
     P : numpy.array of shape (:,3,3) or (3,3)
       1. Piola-Kirchhoff stress.
+    F : numpy.array of shape (:,3,3) or (3,3)
+      Deformation gradient.
 
     """
     if np.shape(F) == np.shape(P) == (3,3):
         S = np.dot(np.linalg.inv(F),P)
     else:
         S = np.einsum('ijk,ikl->ijl',np.linalg.inv(F),P)
-    return S
+    return symmetric(S)
+
+
+def right_stretch(T):
+    """
+    Return the right stretch of a tensor.
+
+    Parameters
+    ----------
+    T : numpy.array of shape (:,3,3) or (3,3)
+      Tensor of which the right stretch is computed.
+
+    """
+    return __polar_decomposition(T,'U')[0]
+
+
+def rotational_part(T):
+    """
+    Return the rotational part of a tensor.
+
+    Parameters
+    ----------
+    T : numpy.array of shape (:,3,3) or (3,3)
+      Tensor of which the rotational part is computed.
+
+    """
+    return __polar_decomposition(T,'R')[0]
+
+
+def spherical_part(T,tensor=False):
+    """
+    Return spherical (hydrostatic) part of a tensor.
+
+    Parameters
+    ----------
+    T : numpy.array of shape (:,3,3) or (3,3)
+      Tensor of which the hydrostatic part is computed.
+    tensor : bool, optional
+      Map spherical part onto identity tensor. Default is false
+
+    """
+    if T.shape == (3,3):
+        sph = np.trace(T)/3.0
+        return sph if not tensor else np.eye(3)*sph
+    else:
+        sph = np.trace(T,axis1=1,axis2=2)/3.0
+        if not tensor:
+            return sph
+        else:
+            return np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),(T.shape[0],3,3)),sph)
 
 
 def strain_tensor(F,t,m):
@@ -78,195 +235,73 @@ def strain_tensor(F,t,m):
            eps
 
 
-def deviatoric_part(x):
-    """
-    Return deviatoric part of a tensor.
-      
-    Parameters
-    ----------
-    x : numpy.array of shape (:,3,3) or (3,3)
-      Tensor of which the deviatoric part is computed.
-
-    """
-    return x - np.eye(3)*spherical_part(x) if np.shape(x) == (3,3) else \
-           x - np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),[x.shape[0],3,3]),spherical_part(x)) 
-
-
-def spherical_part(x,tensor=False):
-    """
-    Return spherical (hydrostatic) part of a tensor.
-    
-    Parameters
-    ----------
-    x : numpy.array of shape (:,3,3) or (3,3)
-      Tensor of which the hydrostatic part is computed.
-    tensor : bool, optional
-      Map spherical part onto identity tensor. Default is false
-
-    """
-    if x.shape == (3,3):
-        sph = np.trace(x)/3.0
-        return sph if not tensor else np.eye(3)*sph
-    else:
-        sph = np.trace(x,axis1=1,axis2=2)/3.0
-        if not tensor:
-            return sph
-        else:
-            return np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),(x.shape[0],3,3)),sph)
-    
-    
-def Mises_stress(sigma):
-    """
-    Return the Mises equivalent of a stress tensor.
-    
-    Parameters
-    ----------
-    sigma : numpy.array of shape (:,3,3) or (3,3)
-      Symmetric stress tensor of which the von Mises equivalent is computed.
-
-    """
-    s = deviatoric_part(sigma)
-    return np.sqrt(3.0/2.0*(np.sum(s**2.0))) if np.shape(sigma) == (3,3) else \
-           np.sqrt(3.0/2.0*np.einsum('ijk->i',s**2.0))
-    
-    
-def Mises_strain(epsilon):
-    """
-    Return the Mises equivalent of a strain tensor.
-    
-    Parameters
-    ----------
-    epsilon : numpy.array of shape (:,3,3) or (3,3)
-      Symmetric strain tensor of which the von Mises equivalent is computed.
-
-    """
-    s = deviatoric_part(epsilon)
-    return np.sqrt(2.0/3.0*(np.sum(s**2.0))) if np.shape(epsilon) == (3,3) else \
-           np.sqrt(2.0/3.0*np.einsum('ijk->i',s**2.0))
-
-
-def symmetric(x):
+def symmetric(T):
     """
     Return the symmetrized tensor.
 
     Parameters
     ----------
-    x : numpy.array of shape (:,3,3) or (3,3)
+    T : numpy.array of shape (:,3,3) or (3,3)
       Tensor of which the symmetrized values are computed.
 
     """
-    return (x+transpose(x))*0.5
+    return (T+transpose(T))*0.5
 
 
-def maximum_shear(x):
-    """
-    Return the maximum shear component of a symmetric tensor.
-    
-    Parameters
-    ----------
-    x : numpy.array of shape (:,3,3) or (3,3)
-      Symmetric tensor of which the maximum shear is computed.
-
-    """
-    w = np.linalg.eigvalsh(symmetric(x))                                                            # eigenvalues in ascending order
-    return (w[2] - w[0])*0.5 if np.shape(x) == (3,3) else \
-           (w[:,2] - w[:,0])*0.5
-    
-    
-def principal_components(x):
-    """
-    Return the principal components of a symmetric tensor.
-    
-    The principal components (eigenvalues) are sorted in descending order, each repeated according to
-    its multiplicity.
-    
-    Parameters
-    ----------
-    x : numpy.array of shape (:,3,3) or (3,3)
-      Symmetric tensor of which the principal compontents are computed.
-
-    """
-    w = np.linalg.eigvalsh(symmetric(x))                                                            # eigenvalues in ascending order
-    return w[::-1] if np.shape(x) == (3,3) else \
-           w[:,::-1]
-    
-    
-def transpose(x):
+def transpose(T):
     """
     Return the transpose of a tensor.
 
     Parameters
     ----------
-    x : numpy.array of shape (:,3,3) or (3,3)
+    T : numpy.array of shape (:,3,3) or (3,3)
       Tensor of which the transpose is computed.
 
     """
-    return x.T if np.shape(x) == (3,3) else \
-           np.transpose(x,(0,2,1))
+    return T.T if np.shape(T) == (3,3) else \
+           np.transpose(T,(0,2,1))
 
 
-def rotational_part(x):
-    """
-    Return the rotational part of a tensor.
-      
-    Parameters
-    ----------
-    x : numpy.array of shape (:,3,3) or (3,3)
-      Tensor of which the rotational part is computed.
-
-    """
-    return __polar_decomposition(x,'R')[0]
-
-
-def left_stretch(x):
-    """
-    Return the left stretch of a tensor.
-      
-    Parameters
-    ----------
-    x : numpy.array of shape (:,3,3) or (3,3)
-      Tensor of which the left stretch is computed.
-
-    """
-    return __polar_decomposition(x,'V')[0]
-  
-  
-def right_stretch(x):
-    """
-    Return the right stretch of a tensor.
-      
-    Parameters
-    ----------
-    x : numpy.array of shape (:,3,3) or (3,3)
-      Tensor of which the right stretch is computed.
-
-    """
-    return __polar_decomposition(x,'U')[0]
-
-
-def __polar_decomposition(x,requested):
+def __polar_decomposition(T,requested):
     """
     Singular value decomposition.
 
     Parameters
     ----------
-    x : numpy.array of shape (:,3,3) or (3,3)
+    T : numpy.array of shape (:,3,3) or (3,3)
       Tensor of which the singular values are computed.
     requested : iterable of str
       Requested outputs: ‘R’ for the rotation tensor,
       ‘V’ for left stretch tensor and ‘U’ for right stretch tensor.
 
     """
-    u, s, vh = np.linalg.svd(x)
-    R = np.dot(u,vh) if np.shape(x) == (3,3) else \
+    u, s, vh = np.linalg.svd(T)
+    R = np.dot(u,vh) if np.shape(T) == (3,3) else \
         np.einsum('ijk,ikl->ijl',u,vh)
 
     output = []
     if 'R' in requested:
         output.append(R)
     if 'V' in requested:
-        output.append(np.dot(x,R.T) if np.shape(x) == (3,3) else np.einsum('ijk,ilk->ijl',x,R))
+        output.append(np.dot(T,R.T) if np.shape(T) == (3,3) else np.einsum('ijk,ilk->ijl',T,R))
     if 'U' in requested:
-        output.append(np.dot(R.T,x) if np.shape(x) == (3,3) else np.einsum('ikj,ikl->ijl',R,x))
+        output.append(np.dot(R.T,T) if np.shape(T) == (3,3) else np.einsum('ikj,ikl->ijl',R,T))
 
     return tuple(output)
+
+
+def __Mises(T_sym,s):
+    """
+    Base equation for Mises equivalent of a stres or strain tensor.
+
+    Parameters
+    ----------
+    T_sym : numpy.array of shape (:,3,3) or (3,3)
+      Symmetric tensor of which the von Mises equivalent is computed.
+    s : float
+      Scaling factor (2/3 for strain, 3/2 for stress).
+
+    """
+    d = deviatoric_part(T_sym)
+    return np.sqrt(s*(np.sum(d**2.0))) if np.shape(T_sym) == (3,3) else \
+           np.sqrt(s*np.einsum('ijk->i',d**2.0))

python/damask/table.py

@@ -32,7 +32,7 @@ class Table():
         """Label data individually, e.g. v v v ==> 1_v 2_v 3_v."""
         labels = []
         for label,shape in self.shapes.items():
-            size = np.prod(shape)
+            size = int(np.prod(shape))
             labels += ['{}{}'.format('' if size == 1 else '{}_'.format(i+1),label) for i in range(size)]
         self.data.columns = labels
         
@@ -41,14 +41,14 @@ class Table():
         """Label data condensed, e.g. 1_v 2_v 3_v ==> v v v."""
         labels = []
         for label,shape in self.shapes.items():
-            labels += [label] * np.prod(shape)
+            labels += [label] * int(np.prod(shape))
         self.data.columns = labels
 
 
     def __add_comment(self,label,shape,info):
         if info is not None:
             self.comments.append('{}{}: {}'.format(label,
-                                                   ' '+str(shape) if np.prod(shape) > 1 else '',
+                                                   ' '+str(shape) if np.prod(shape,dtype=int) > 1 else '',
                                                    info))
 
     

python/damask/util.py

@@ -3,9 +3,11 @@ import time
 import os
 import subprocess
 import shlex
+from fractions import Fraction
+from functools import reduce
 from optparse import Option
-from queue import Queue
-from threading import Thread
+
+import numpy as np
 
 class bcolors:
     """
@@ -104,6 +106,7 @@ def strikeout(what):
     """Formats string as strikeout."""
     return bcolors.CROSSOUT+srepr(what)+bcolors.ENDC
 
+
 def execute(cmd,
             streamIn = None,
             wd = './',
@@ -207,6 +210,24 @@ def progressBar(iteration, total, prefix='', bar_length=50):
     sys.stderr.flush()
  
 
+def scale_to_coprime(v):
+    """Scale vector to co-prime (relatively prime) integers."""
+    MAX_DENOMINATOR = 1000
+  
+    def get_square_denominator(x):
+        """Denominator of the square of a number."""
+        return Fraction(x ** 2).limit_denominator(MAX_DENOMINATOR).denominator
+  
+    def lcm(a, b):
+        """Least common multiple."""
+        return a * b // np.gcd(a, b)
+  
+    denominators = [int(get_square_denominator(i)) for i in v]
+    s = reduce(lcm, denominators) ** 0.5
+    m = (np.array(v)*s).astype(np.int)
+    return m//reduce(np.gcd,m)
+
+
 class return_message():
     """Object with formatted return message."""
     
@@ -226,56 +247,3 @@ class return_message():
         """Return message suitable for interactive shells."""
         return srepr(self.message)
 
-
-class ThreadPool:
-  """Pool of threads consuming tasks from a queue."""
-
-  class Worker(Thread):
-    """Thread executing tasks from a given tasks queue."""
-
-    def __init__(self, tasks):
-      """Worker for tasks."""
-      Thread.__init__(self)
-      self.tasks = tasks
-      self.daemon = True
-      self.start()
-
-    def run(self):
-      while True:
-        func, args, kargs = self.tasks.get()
-        try:
-          func(*args, **kargs)
-        except Exception as e:
-          # An exception happened in this thread
-          print(e)
-        finally:
-          # Mark this task as done, whether an exception happened or not
-          self.tasks.task_done()
-
-
-  def __init__(self, num_threads):
-    """
-    Thread pool.
-
-    Parameters
-    ----------
-    num_threads : int
-      number of threads
-
-    """
-    self.tasks = Queue(num_threads)
-    for _ in range(num_threads):
-      self.Worker(self.tasks)
-
-  def add_task(self, func, *args, **kargs):
-    """Add a task to the queue."""
-    self.tasks.put((func, args, kargs))
-
-  def map(self, func, args_list):
-    """Add a list of tasks to the queue."""
-    for args in args_list:
-      self.add_task(func, args)
-
-  def wait_completion(self):
-    """Wait for completion of all the tasks in the queue."""
-    self.tasks.join()

python/tests/test_DADF5.py

@@ -1,90 +0,0 @@
-import shutil
-import os
-
-import pytest
-import numpy as np
-
-from damask import DADF5
-from damask import mechanics
-
-@pytest.fixture
-def default(tmp_path,reference_dir):
-    """Small DADF5 file in temp location for modification."""
-    fname = '12grains6x7x8_tensionY.hdf5'
-    shutil.copy(os.path.join(reference_dir,fname),tmp_path)
-    f = DADF5(os.path.join(tmp_path,fname))
-    f.set_by_time(20.0,20.0)
-    return f
-
-@pytest.fixture
-def reference_dir(reference_dir_base):
-    """Directory containing reference results."""
-    return os.path.join(reference_dir_base,'DADF5')
-
-
-class TestDADF5:
-
-    def test_time_increments(self,default):
-        shape = default.read_dataset(default.get_dataset_location('F'),0).shape
-        default.set_by_time(0.0,20.0)
-        for i in default.iter_visible('increments'):
-           assert shape == default.read_dataset(default.get_dataset_location('F'),0).shape
-
-
-    def test_add_absolute(self,default):
-        default.add_absolute('Fe')
-        loc = {'Fe':   default.get_dataset_location('Fe'),
-               '|Fe|': default.get_dataset_location('|Fe|')}
-        in_memory = np.abs(default.read_dataset(loc['Fe'],0))
-        in_file   = default.read_dataset(loc['|Fe|'],0)
-        assert np.allclose(in_memory,in_file)
-
-    def test_add_calculation(self,default):
-        default.add_calculation('2.0*np.abs(#F#)-1.0','x','-','test')
-        loc = {'F':    default.get_dataset_location('F'),
-               'x':    default.get_dataset_location('x')}
-        in_memory = 2.0*np.abs(default.read_dataset(loc['F'],0))-1.0
-        in_file   = default.read_dataset(loc['x'],0)
-        assert np.allclose(in_memory,in_file)
-
-    def test_add_Cauchy(self,default):
-        default.add_Cauchy('P','F')
-        loc = {'F':    default.get_dataset_location('F'),
-               'P':    default.get_dataset_location('P'),
-               'sigma':default.get_dataset_location('sigma')}
-        in_memory = mechanics.Cauchy(default.read_dataset(loc['F'],0),
-                                     default.read_dataset(loc['P'],0))
-        in_file   = default.read_dataset(loc['sigma'],0)
-        assert np.allclose(in_memory,in_file)
-
-    def test_add_determinant(self,default):
-        default.add_determinant('P')
-        loc = {'P':     default.get_dataset_location('P'),
-               'det(P)':default.get_dataset_location('det(P)')}
-        in_memory = np.linalg.det(default.read_dataset(loc['P'],0)).reshape((-1,1))
-        in_file   = default.read_dataset(loc['det(P)'],0)
-        assert np.allclose(in_memory,in_file)
-
-    def test_add_deviator(self,default):
-        default.add_deviator('P')
-        loc = {'P'  :default.get_dataset_location('P'),
-               's_P':default.get_dataset_location('s_P')}
-        in_memory = mechanics.deviatoric_part(default.read_dataset(loc['P'],0))
-        in_file   = default.read_dataset(loc['s_P'],0)
-        assert np.allclose(in_memory,in_file)
-
-    def test_add_norm(self,default):
-        default.add_norm('F',1)
-        loc = {'F':    default.get_dataset_location('F'),
-               '|F|_1':default.get_dataset_location('|F|_1')}
-        in_memory = np.linalg.norm(default.read_dataset(loc['F'],0),ord=1,axis=(1,2),keepdims=True)
-        in_file   = default.read_dataset(loc['|F|_1'],0)
-        assert np.allclose(in_memory,in_file)
-
-    def test_add_spherical(self,default):
-        default.add_spherical('P')
-        loc = {'P':   default.get_dataset_location('P'),
-               'p_P': default.get_dataset_location('p_P')}
-        in_memory = mechanics.spherical_part(default.read_dataset(loc['P'],0)).reshape(-1,1)
-        in_file   = default.read_dataset(loc['p_P'],0)
-        assert np.allclose(in_memory,in_file)

python/tests/test_Result.py

@@ -0,0 +1,183 @@
+import shutil
+import os
+
+import pytest
+import numpy as np
+
+from damask import Result
+from damask import mechanics
+
+@pytest.fixture
+def default(tmp_path,reference_dir):
+    """Small Result file in temp location for modification."""
+    fname = '12grains6x7x8_tensionY.hdf5'
+    shutil.copy(os.path.join(reference_dir,fname),tmp_path)
+    f = Result(os.path.join(tmp_path,fname))
+    f.set_by_time(20.0,20.0)
+    return f
+
+@pytest.fixture
+def reference_dir(reference_dir_base):
+    """Directory containing reference results."""
+    return os.path.join(reference_dir_base,'Result')
+
+
+class TestResult:
+
+    def test_time_increments(self,default):
+        shape = default.read_dataset(default.get_dataset_location('F'),0).shape
+        default.set_by_time(0.0,20.0)
+        for i in default.iter_visible('increments'):
+           assert shape == default.read_dataset(default.get_dataset_location('F'),0).shape
+
+
+    def test_add_absolute(self,default):
+        default.add_absolute('Fe')
+        loc = {'Fe':   default.get_dataset_location('Fe'),
+               '|Fe|': default.get_dataset_location('|Fe|')}
+        in_memory = np.abs(default.read_dataset(loc['Fe'],0))
+        in_file   = default.read_dataset(loc['|Fe|'],0)
+        assert np.allclose(in_memory,in_file)
+
+    def test_add_calculation(self,default):
+        default.add_calculation('x','2.0*np.abs(#F#)-1.0','-','my notes')
+        loc = {'F':    default.get_dataset_location('F'),
+               'x':    default.get_dataset_location('x')}
+        in_memory = 2.0*np.abs(default.read_dataset(loc['F'],0))-1.0
+        in_file   = default.read_dataset(loc['x'],0)
+        assert np.allclose(in_memory,in_file)
+
+    def test_add_Cauchy(self,default):
+        default.add_Cauchy('P','F')
+        loc = {'F':    default.get_dataset_location('F'),
+               'P':    default.get_dataset_location('P'),
+               'sigma':default.get_dataset_location('sigma')}
+        in_memory = mechanics.Cauchy(default.read_dataset(loc['P'],0),
+                                     default.read_dataset(loc['F'],0))
+        in_file   = default.read_dataset(loc['sigma'],0)
+        assert np.allclose(in_memory,in_file)
+
+    def test_add_determinant(self,default):
+        default.add_determinant('P')
+        loc = {'P':     default.get_dataset_location('P'),
+               'det(P)':default.get_dataset_location('det(P)')}
+        in_memory = np.linalg.det(default.read_dataset(loc['P'],0)).reshape((-1,1))
+        in_file   = default.read_dataset(loc['det(P)'],0)
+        assert np.allclose(in_memory,in_file)
+
+    def test_add_deviator(self,default):
+        default.add_deviator('P')
+        loc = {'P'  :default.get_dataset_location('P'),
+               's_P':default.get_dataset_location('s_P')}
+        in_memory = mechanics.deviatoric_part(default.read_dataset(loc['P'],0))
+        in_file   = default.read_dataset(loc['s_P'],0)
+        assert np.allclose(in_memory,in_file)
+
+    def test_add_eigenvalues(self,default):
+        default.add_Cauchy('P','F')
+        default.add_eigenvalues('sigma')
+        loc = {'sigma'        :default.get_dataset_location('sigma'),
+               'lambda(sigma)':default.get_dataset_location('lambda(sigma)')}
+        in_memory = mechanics.eigenvalues(default.read_dataset(loc['sigma'],0))
+        in_file   = default.read_dataset(loc['lambda(sigma)'],0)
+        assert np.allclose(in_memory,in_file)
+
+    def test_add_eigenvectors(self,default):
+        default.add_Cauchy('P','F')
+        default.add_eigenvectors('sigma')
+        loc = {'sigma'   :default.get_dataset_location('sigma'),
+               'v(sigma)':default.get_dataset_location('v(sigma)')}
+        in_memory = mechanics.eigenvectors(default.read_dataset(loc['sigma'],0))
+        in_file   = default.read_dataset(loc['v(sigma)'],0)
+        assert np.allclose(in_memory,in_file)
+
+    def test_add_maximum_shear(self,default):
+        default.add_Cauchy('P','F')
+        default.add_maximum_shear('sigma')
+        loc = {'sigma'           :default.get_dataset_location('sigma'),
+               'max_shear(sigma)':default.get_dataset_location('max_shear(sigma)')}
+        in_memory = mechanics.maximum_shear(default.read_dataset(loc['sigma'],0)).reshape(-1,1)
+        in_file   = default.read_dataset(loc['max_shear(sigma)'],0)
+        assert np.allclose(in_memory,in_file)
+
+    def test_add_Mises_strain(self,default):
+        t = ['V','U'][np.random.randint(0,2)]
+        m = np.random.random()*2.0 - 1.0
+        default.add_strain_tensor('F',t,m)
+        label = 'epsilon_{}^{}(F)'.format(t,m)
+        default.add_Mises(label)
+        loc = {label      :default.get_dataset_location(label),
+               label+'_vM':default.get_dataset_location(label+'_vM')}
+        in_memory = mechanics.Mises_strain(default.read_dataset(loc[label],0)).reshape(-1,1)
+        in_file   = default.read_dataset(loc[label+'_vM'],0)
+        assert np.allclose(in_memory,in_file)
+
+    def test_add_Mises_stress(self,default):
+        default.add_Cauchy('P','F')
+        default.add_Mises('sigma')
+        loc = {'sigma'   :default.get_dataset_location('sigma'),
+               'sigma_vM':default.get_dataset_location('sigma_vM')}
+        in_memory = mechanics.Mises_stress(default.read_dataset(loc['sigma'],0)).reshape(-1,1)
+        in_file   = default.read_dataset(loc['sigma_vM'],0)
+        assert np.allclose(in_memory,in_file)
+
+    def test_add_norm(self,default):
+        default.add_norm('F',1)
+        loc = {'F':    default.get_dataset_location('F'),
+               '|F|_1':default.get_dataset_location('|F|_1')}
+        in_memory = np.linalg.norm(default.read_dataset(loc['F'],0),ord=1,axis=(1,2),keepdims=True)
+        in_file   = default.read_dataset(loc['|F|_1'],0)
+        assert np.allclose(in_memory,in_file)
+
+    def test_add_PK2(self,default):
+        default.add_PK2('P','F')
+        loc = {'F':default.get_dataset_location('F'),
+               'P':default.get_dataset_location('P'),
+               'S':default.get_dataset_location('S')}
+        in_memory = mechanics.PK2(default.read_dataset(loc['P'],0),
+                                  default.read_dataset(loc['F'],0))
+        in_file   = default.read_dataset(loc['S'],0)
+        assert np.allclose(in_memory,in_file)
+
+    def test_add_rotational_part(self,default):
+        default.add_rotational_part('F')
+        loc = {'F':    default.get_dataset_location('F'),
+               'R(F)': default.get_dataset_location('R(F)')}
+        in_memory = mechanics.rotational_part(default.read_dataset(loc['F'],0))
+        in_file   = default.read_dataset(loc['R(F)'],0)
+        assert np.allclose(in_memory,in_file)
+
+    def test_add_spherical(self,default):
+        default.add_spherical('P')
+        loc = {'P':   default.get_dataset_location('P'),
+               'p_P': default.get_dataset_location('p_P')}
+        in_memory = mechanics.spherical_part(default.read_dataset(loc['P'],0)).reshape(-1,1)
+        in_file   = default.read_dataset(loc['p_P'],0)
+        assert np.allclose(in_memory,in_file)
+
+    def test_add_strain(self,default):
+        t = ['V','U'][np.random.randint(0,2)]
+        m = np.random.random()*2.0 - 1.0
+        default.add_strain_tensor('F',t,m)
+        label = 'epsilon_{}^{}(F)'.format(t,m)
+        loc = {'F':   default.get_dataset_location('F'),
+               label: default.get_dataset_location(label)}
+        in_memory = mechanics.strain_tensor(default.read_dataset(loc['F'],0),t,m)
+        in_file   = default.read_dataset(loc[label],0)
+        assert np.allclose(in_memory,in_file)
+
+    def test_add_stretch_right(self,default):
+        default.add_stretch_tensor('F','U')
+        loc = {'F':    default.get_dataset_location('F'),
+               'U(F)': default.get_dataset_location('U(F)')}
+        in_memory = mechanics.right_stretch(default.read_dataset(loc['F'],0))
+        in_file   = default.read_dataset(loc['U(F)'],0)
+        assert np.allclose(in_memory,in_file)
+
+    def test_add_stretch_left(self,default):
+        default.add_stretch_tensor('F','V')
+        loc = {'F':    default.get_dataset_location('F'),
+               'V(F)': default.get_dataset_location('V(F)')}
+        in_memory = mechanics.left_stretch(default.read_dataset(loc['F'],0))
+        in_file   = default.read_dataset(loc['V(F)'],0)
+        assert np.allclose(in_memory,in_file)

python/tests/test_mechanics.py

@@ -8,181 +8,218 @@ class TestMechanics:
 
 
     def test_vectorize_Cauchy(self):
-         P = np.random.random((self.n,3,3))
-         F = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.Cauchy(F,P)[self.c],
-                            mechanics.Cauchy(F[self.c],P[self.c]))
-
-
-    def test_vectorize_strain_tensor(self):
-         F = np.random.random((self.n,3,3))
-         t = ['V','U'][np.random.randint(0,2)]
-         m = np.random.random()*10. -5.0
-         assert np.allclose(mechanics.strain_tensor(F,t,m)[self.c],
-                            mechanics.strain_tensor(F[self.c],t,m))
-
+        P = np.random.random((self.n,3,3))
+        F = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.Cauchy(P,F)[self.c],
+                           mechanics.Cauchy(P[self.c],F[self.c]))
 
     def test_vectorize_deviatoric_part(self):
-         x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.deviatoric_part(x)[self.c],
-                            mechanics.deviatoric_part(x[self.c]))
+        x = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.deviatoric_part(x)[self.c],
+                           mechanics.deviatoric_part(x[self.c]))
 
+    def test_vectorize_eigenvalues(self):
+        x = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.eigenvalues(x)[self.c],
+                           mechanics.eigenvalues(x[self.c]))
 
-    def test_vectorize_spherical_part(self):
-         x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.spherical_part(x,True)[self.c],
-                            mechanics.spherical_part(x[self.c],True))
-
-
-    def test_vectorize_Mises_stress(self):
-         sigma = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.Mises_stress(sigma)[self.c],
-                            mechanics.Mises_stress(sigma[self.c]))
-
-
-    def test_vectorize_Mises_strain(self):
-         epsilon = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.Mises_strain(epsilon)[self.c],
-                            mechanics.Mises_strain(epsilon[self.c]))
-
-
-    def test_vectorize_symmetric(self):
-         x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.symmetric(x)[self.c],
-                            mechanics.symmetric(x[self.c]))
-
-
-    def test_vectorize_maximum_shear(self):
-         x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.maximum_shear(x)[self.c],
-                            mechanics.maximum_shear(x[self.c]))
-
-
-    def test_vectorize_principal_components(self):
-         x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.principal_components(x)[self.c],
-                            mechanics.principal_components(x[self.c]))
-
-
-    def test_vectorize_transpose(self):
-         x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.transpose(x)[self.c],
-                            mechanics.transpose(x[self.c]))
-
-
-    def test_vectorize_rotational_part(self):
-         x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.rotational_part(x)[self.c],
-                            mechanics.rotational_part(x[self.c]))
-
+    def test_vectorize_eigenvectors(self):
+        x = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.eigenvectors(x)[self.c],
+                           mechanics.eigenvectors(x[self.c]))
 
     def test_vectorize_left_stretch(self):
-         x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.left_stretch(x)[self.c],
-                            mechanics.left_stretch(x[self.c]))
+        x = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.left_stretch(x)[self.c],
+                           mechanics.left_stretch(x[self.c]))
 
+    def test_vectorize_maximum_shear(self):
+        x = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.maximum_shear(x)[self.c],
+                           mechanics.maximum_shear(x[self.c]))
+
+    def test_vectorize_Mises_strain(self):
+        epsilon = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.Mises_strain(epsilon)[self.c],
+                           mechanics.Mises_strain(epsilon[self.c]))
+
+    def test_vectorize_Mises_stress(self):
+        sigma = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.Mises_stress(sigma)[self.c],
+                           mechanics.Mises_stress(sigma[self.c]))
+
+    def test_vectorize_PK2(self):
+        F = np.random.random((self.n,3,3))
+        P = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.PK2(P,F)[self.c],
+                           mechanics.PK2(P[self.c],F[self.c]))
 
     def test_vectorize_right_stretch(self):
-         x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.right_stretch(x)[self.c],
-                            mechanics.right_stretch(x[self.c]))
+        x = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.right_stretch(x)[self.c],
+                           mechanics.right_stretch(x[self.c]))
+
+    def test_vectorize_rotational_part(self):
+        x = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.rotational_part(x)[self.c],
+                           mechanics.rotational_part(x[self.c]))
+
+    def test_vectorize_spherical_part(self):
+        x = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.spherical_part(x,True)[self.c],
+                           mechanics.spherical_part(x[self.c],True))
+
+    def test_vectorize_strain_tensor(self):
+        F = np.random.random((self.n,3,3))
+        t = ['V','U'][np.random.randint(0,2)]
+        m = np.random.random()*10. -5.0
+        assert np.allclose(mechanics.strain_tensor(F,t,m)[self.c],
+                           mechanics.strain_tensor(F[self.c],t,m))
+
+    def test_vectorize_symmetric(self):
+        x = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.symmetric(x)[self.c],
+                           mechanics.symmetric(x[self.c]))
+
+    def test_vectorize_transpose(self):
+        x = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.transpose(x)[self.c],
+                           mechanics.transpose(x[self.c]))
 
 
     def test_Cauchy(self):
-         """Ensure Cauchy stress is symmetrized 1. Piola-Kirchhoff stress for no deformation."""
-         P = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.Cauchy(np.broadcast_to(np.eye(3),(self.n,3,3)),P),
-                            mechanics.symmetric(P))
+        """Ensure Cauchy stress is symmetrized 1. Piola-Kirchhoff stress for no deformation."""
+        P = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.Cauchy(P,np.broadcast_to(np.eye(3),(self.n,3,3))),
+                           mechanics.symmetric(P))
+
 
     def test_polar_decomposition(self):
-         """F = RU = VR."""
-         F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.random((self.n,3,3))
-         R = mechanics.rotational_part(F)
-         V = mechanics.left_stretch(F)
-         U = mechanics.right_stretch(F)
-         assert np.allclose(np.matmul(R,U),
-                            np.matmul(V,R))
+        """F = RU = VR."""
+        F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.random((self.n,3,3))
+        R = mechanics.rotational_part(F)
+        V = mechanics.left_stretch(F)
+        U = mechanics.right_stretch(F)
+        assert np.allclose(np.matmul(R,U),
+                           np.matmul(V,R))
+
+
+    def test_PK2(self):
+        """Ensure 2. Piola-Kirchhoff stress is symmetrized 1. Piola-Kirchhoff stress for no deformation."""
+        P = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.PK2(P,np.broadcast_to(np.eye(3),(self.n,3,3))),
+                           mechanics.symmetric(P))
 
 
     def test_strain_tensor_no_rotation(self):
-         """Ensure that left and right stretch give same results for no rotation."""
-         F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.random((self.n,3,3))
-         m = np.random.random()*20.0-10.0
-         assert np.allclose(mechanics.strain_tensor(F,'U',m),
-                            mechanics.strain_tensor(F,'V',m))
+        """Ensure that left and right stretch give same results for no rotation."""
+        F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.random((self.n,3,3))
+        m = np.random.random()*20.0-10.0
+        assert np.allclose(mechanics.strain_tensor(F,'U',m),
+                           mechanics.strain_tensor(F,'V',m))
 
     def test_strain_tensor_rotation_equivalence(self):
-         """Ensure that left and right strain differ only by a rotation."""
-         F = np.broadcast_to(np.eye(3),[self.n,3,3]) + (np.random.random((self.n,3,3))*0.5 - 0.25)
-         m = np.random.random()*5.0-2.5
-         assert np.allclose(np.linalg.det(mechanics.strain_tensor(F,'U',m)),
-                            np.linalg.det(mechanics.strain_tensor(F,'V',m)))
+        """Ensure that left and right strain differ only by a rotation."""
+        F = np.broadcast_to(np.eye(3),[self.n,3,3]) + (np.random.random((self.n,3,3))*0.5 - 0.25)
+        m = np.random.random()*5.0-2.5
+        assert np.allclose(np.linalg.det(mechanics.strain_tensor(F,'U',m)),
+                           np.linalg.det(mechanics.strain_tensor(F,'V',m)))
 
     def test_strain_tensor_rotation(self):
-         """Ensure that pure rotation results in no strain."""
-         F = mechanics.rotational_part(np.random.random((self.n,3,3)))
-         t = ['V','U'][np.random.randint(0,2)]
-         m = np.random.random()*2.0 - 1.0
-         assert np.allclose(mechanics.strain_tensor(F,t,m),
-                            0.0)
+        """Ensure that pure rotation results in no strain."""
+        F = mechanics.rotational_part(np.random.random((self.n,3,3)))
+        t = ['V','U'][np.random.randint(0,2)]
+        m = np.random.random()*2.0 - 1.0
+        assert np.allclose(mechanics.strain_tensor(F,t,m),
+                           0.0)
 
     def test_rotation_determinant(self):
-         """
-         Ensure that the determinant of the rotational part is +- 1.
+        """
+        Ensure that the determinant of the rotational part is +- 1.
 
-         Should be +1, but random F might contain a reflection.
-         """
-         x = np.random.random((self.n,3,3))
-         assert np.allclose(np.abs(np.linalg.det(mechanics.rotational_part(x))),
-                            1.0)
+        Should be +1, but random F might contain a reflection.
+        """
+        x = np.random.random((self.n,3,3))
+        assert np.allclose(np.abs(np.linalg.det(mechanics.rotational_part(x))),
+                           1.0)
 
 
     def test_spherical_deviatoric_part(self):
-         """Ensure that full tensor is sum of spherical and deviatoric part."""
-         x = np.random.random((self.n,3,3))
-         sph = mechanics.spherical_part(x,True)
-         assert np.allclose(sph + mechanics.deviatoric_part(x),
-                            x)
+        """Ensure that full tensor is sum of spherical and deviatoric part."""
+        x = np.random.random((self.n,3,3))
+        sph = mechanics.spherical_part(x,True)
+        assert np.allclose(sph + mechanics.deviatoric_part(x),
+                           x)
 
     def test_deviatoric_Mises(self):
-         """Ensure that Mises equivalent stress depends only on deviatoric part."""
-         x = np.random.random((self.n,3,3))
-         full = mechanics.Mises_stress(x)
-         dev  = mechanics.Mises_stress(mechanics.deviatoric_part(x))
-         assert np.allclose(full,
-                            dev)
+        """Ensure that Mises equivalent stress depends only on deviatoric part."""
+        x = np.random.random((self.n,3,3))
+        full = mechanics.Mises_stress(x)
+        dev  = mechanics.Mises_stress(mechanics.deviatoric_part(x))
+        assert np.allclose(full,
+                           dev)
 
     def test_spherical_mapping(self):
-         """Ensure that mapping to tensor is correct."""
-         x = np.random.random((self.n,3,3))
-         tensor = mechanics.spherical_part(x,True)
-         scalar = mechanics.spherical_part(x)
-         assert np.allclose(np.linalg.det(tensor),
-                            scalar**3.0)
+        """Ensure that mapping to tensor is correct."""
+        x = np.random.random((self.n,3,3))
+        tensor = mechanics.spherical_part(x,True)
+        scalar = mechanics.spherical_part(x)
+        assert np.allclose(np.linalg.det(tensor),
+                           scalar**3.0)
 
     def test_spherical_Mises(self):
-         """Ensure that Mises equivalent strrain of spherical strain is 0."""
-         x = np.random.random((self.n,3,3))
-         sph = mechanics.spherical_part(x,True)
-         assert np.allclose(mechanics.Mises_strain(sph),
-                            0.0)
+        """Ensure that Mises equivalent strrain of spherical strain is 0."""
+        x = np.random.random((self.n,3,3))
+        sph = mechanics.spherical_part(x,True)
+        assert np.allclose(mechanics.Mises_strain(sph),
+                           0.0)
 
     def test_symmetric(self):
-         """Ensure that a symmetric tensor is half of the sum of a tensor and its transpose."""
-         x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.symmetric(x)*2.0,
-                            mechanics.transpose(x)+x)
+        """Ensure that a symmetric tensor is half of the sum of a tensor and its transpose."""
+        x = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.symmetric(x)*2.0,
+                           mechanics.transpose(x)+x)
 
 
     def test_transpose(self):
-         """Ensure that a symmetric tensor equals its transpose."""
-         x = mechanics.symmetric(np.random.random((self.n,3,3)))
-         assert np.allclose(mechanics.transpose(x),
-                            x)
+        """Ensure that a symmetric tensor equals its transpose."""
+        x = mechanics.symmetric(np.random.random((self.n,3,3)))
+        assert np.allclose(mechanics.transpose(x),
+                           x)
 
 
     def test_Mises(self):
-         """Ensure that equivalent stress is 3/2 of equivalent strain."""
-         x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.Mises_stress(x)/mechanics.Mises_strain(x),
-                            1.5)
+        """Ensure that equivalent stress is 3/2 of equivalent strain."""
+        x = np.random.random((self.n,3,3))
+        assert np.allclose(mechanics.Mises_stress(x)/mechanics.Mises_strain(x),
+                           1.5)
+
+
+    def test_eigenvalues(self):
+        """Ensure that the characteristic polynomial can be solved."""
+        A = mechanics.symmetric(np.random.random((self.n,3,3)))
+        lambd = mechanics.eigenvalues(A)
+        s = np.random.randint(self.n)
+        for i in range(3):
+           assert np.allclose(np.linalg.det(A[s]-lambd[s,i]*np.eye(3)),.0)
+
+    def test_eigenvalues_and_vectors(self):
+        """Ensure that eigenvalues and -vectors are the solution to the characteristic polynomial."""
+        A = mechanics.symmetric(np.random.random((self.n,3,3)))
+        lambd = mechanics.eigenvalues(A)
+        x     = mechanics.eigenvectors(A)
+        s = np.random.randint(self.n)
+        for i in range(3):
+           assert np.allclose(np.dot(A[s]-lambd[s,i]*np.eye(3),x[s,:,i]),.0)
+
+    def test_eigenvectors_RHS(self):
+        """Ensure that RHS coordinate system does only change sign of determinant."""
+        A = mechanics.symmetric(np.random.random((self.n,3,3)))
+        LRHS = np.linalg.det(mechanics.eigenvectors(A,RHS=False))
+        RHS  = np.linalg.det(mechanics.eigenvectors(A,RHS=True))
+        assert np.allclose(np.abs(LRHS),RHS)
+
+    def test_spherical_no_shear(self):
+        """Ensure that sherical stress has max shear of 0.0."""
+        A = mechanics.spherical_part(mechanics.symmetric(np.random.random((self.n,3,3))),True)
+        assert np.allclose(mechanics.maximum_shear(A),0.0)

src/IO.f90

@@ -400,6 +400,8 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
       msg = 'number of values does not match'
     case (147)
       msg = 'not supported anymore'
+    case (148)
+      msg = 'Nconstituents mismatch between homogenization and microstructure'
 
 !--------------------------------------------------------------------------------------------------
 ! material error messages and related messages in mesh

src/constitutive_plastic_disloUCLA.f90

@@ -10,17 +10,6 @@ submodule(constitutive) plastic_disloUCLA
   real(pReal), parameter :: &
     kB = 1.38e-23_pReal                                                                             !< Boltzmann constant in J/Kelvin
 
-  enum, bind(c)
-    enumerator :: &
-      undefined_ID, &
-      rho_mob_ID, &
-      rho_dip_ID, &
-      dot_gamma_sl_ID, &
-      gamma_sl_ID, &
-      Lambda_sl_ID, &
-      tau_pass_ID
-  end enum
- 
   type :: tParameters
     real(pReal) :: &
       aTol_rho, &
@@ -28,7 +17,7 @@ submodule(constitutive) plastic_disloUCLA
       mu, &
       D_0, &                                                                                        !< prefactor for self-diffusion coefficient
       Q_cl                                                                                          !< activation energy for dislocation climb
-    real(pReal),                 dimension(:),  allocatable :: &
+    real(pReal),               allocatable, dimension(:) :: &
       rho_mob_0, &                                                                                  !< initial dislocation density
       rho_dip_0, &                                                                                  !< initial dipole density
       b_sl, &                                                                                       !< magnitude of burgers vector [m]
@@ -46,19 +35,19 @@ submodule(constitutive) plastic_disloUCLA
       kink_height, &                                                                                !< height of the kink pair
       w, &                                                                                          !< width of the kink pair
       omega                                                                                         !< attempt frequency for kink pair nucleation
-    real(pReal),                  dimension(:,:),  allocatable :: &
+    real(pReal),               allocatable, dimension(:,:) :: &
       h_sl_sl, &                                                                                    !< slip resistance from slip activity
       forestProjectionEdge
-    real(pReal),                  dimension(:,:,:),  allocatable :: &
+    real(pReal),               allocatable, dimension(:,:,:) :: &
       Schmid, &
       nonSchmid_pos, &
       nonSchmid_neg
     integer :: &
       sum_N_sl                                                                                      !< total number of active slip system
-    integer,                      dimension(:), allocatable :: &
+    integer,                   allocatable, dimension(:) :: &
       N_sl                                                                                          !< number of active slip systems for each family
-    integer(kind(undefined_ID)),  dimension(:),allocatable :: &
-      outputID                                                                                      !< ID of each post result output
+    character(len=pStringLen), allocatable, dimension(:) :: &
+      output
     logical :: &
       dipoleFormation                                                                               !< flag indicating consideration of dipole formation
   end type                                                                                          !< container type for internal constitutive parameters
@@ -88,7 +77,7 @@ contains
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief module initialization
+!> @brief Perform module initialization.
 !> @details reads in material parameters, allocates arrays, and does sanity checks
 !--------------------------------------------------------------------------------------------------
 module subroutine plastic_disloUCLA_init
@@ -100,18 +89,13 @@ module subroutine plastic_disloUCLA_init
     sizeState, sizeDotState, &
     startIndex, endIndex
 
-  integer(kind(undefined_ID)) :: &
-    outputID
- 
   character(len=pStringLen) :: &
     extmsg = ''
-  character(len=pStringLen), dimension(:), allocatable :: &
-    outputs
 
-  write(6,'(/,a)')   ' <<<+-  plastic_'//PLASTICITY_DISLOUCLA_label//' init  -+>>>'
+  write(6,'(/,a)') ' <<<+-  plastic_'//PLASTICITY_DISLOUCLA_label//' init  -+>>>'; flush(6)
 
-  write(6,'(/,a)')   ' Cereceda et al., International Journal of Plasticity 78:242–256, 2016'
-  write(6,'(a)')     ' https://dx.doi.org/10.1016/j.ijplas.2015.09.002'
+  write(6,'(/,a)') ' Cereceda et al., International Journal of Plasticity 78:242–256, 2016'
+  write(6,'(a)')   ' https://dx.doi.org/10.1016/j.ijplas.2015.09.002'
 
   Ninstance = count(phase_plasticity == PLASTICITY_DISLOUCLA_ID)
   if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
@@ -232,32 +216,7 @@ module subroutine plastic_disloUCLA_init
 
 !--------------------------------------------------------------------------------------------------
 !  output pararameters
-    outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
-    allocate(prm%outputID(0))
-    do i=1, size(outputs)
-      outputID = undefined_ID
-      select case(trim(outputs(i)))
- 
-        case ('edge_density')
-          outputID  = merge(rho_mob_ID,undefined_ID,prm%sum_N_sl>0)
-        case ('dipole_density')
-          outputID = merge(rho_dip_ID,undefined_ID,prm%sum_N_sl>0)
-        case ('shear_rate','shearrate','shear_rate_slip','shearrate_slip')
-          outputID = merge(dot_gamma_sl_ID,undefined_ID,prm%sum_N_sl>0)
-        case ('accumulated_shear','accumulatedshear','accumulated_shear_slip')
-          outputID = merge(gamma_sl_ID,undefined_ID,prm%sum_N_sl>0)
-        case ('mfp','mfp_slip')
-          outputID = merge(Lambda_sl_ID,undefined_ID,prm%sum_N_sl>0)
-        case ('threshold_stress','threshold_stress_slip')
-          outputID = merge(tau_pass_ID,undefined_ID,prm%sum_N_sl>0)
- 
-      end select
- 
-      if (outputID /= undefined_ID) then
-        prm%outputID = [prm%outputID, outputID]
-      endif
- 
-    enddo
+    prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
@@ -304,7 +263,7 @@ end subroutine plastic_disloUCLA_init
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief calculates plastic velocity gradient and its tangent
+!> @brief Calculate plastic velocity gradient and its tangent.
 !--------------------------------------------------------------------------------------------------
 pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp, &
                                                   Mp,T,instance,of)
@@ -347,7 +306,7 @@ end subroutine plastic_disloUCLA_LpAndItsTangent
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief calculates the rate of change of microstructure
+!> @brief Calculate the rate of change of microstructure.
 !--------------------------------------------------------------------------------------------------
 module subroutine plastic_disloUCLA_dotState(Mp,T,instance,of)
 
@@ -407,7 +366,7 @@ end subroutine plastic_disloUCLA_dotState
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief calculates derived quantities from state
+!> @brief Calculate derived quantities from state.
 !--------------------------------------------------------------------------------------------------
 module subroutine plastic_disloUCLA_dependentState(instance,of)
 
@@ -433,7 +392,7 @@ end subroutine plastic_disloUCLA_dependentState
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief writes results to HDF5 output file
+!> @brief Write results to HDF5 output file.
 !--------------------------------------------------------------------------------------------------
 module subroutine plastic_disloUCLA_results(instance,group)
 
@@ -443,23 +402,23 @@ module subroutine plastic_disloUCLA_results(instance,group)
   integer :: o
 
   associate(prm => param(instance), stt => state(instance), dst => dependentState(instance))
-  outputsLoop: do o = 1,size(prm%outputID)
-    select case(prm%outputID(o))
-      case (rho_mob_ID)
-        call results_writeDataset(group,stt%rho_mob,'rho_mob',&
-                                 'mobile dislocation density','1/m²')
-      case (rho_dip_ID)
-        call results_writeDataset(group,stt%rho_dip,'rho_dip',&
-                                  'dislocation dipole density''1/m²')
-      case (dot_gamma_sl_ID)
-        call results_writeDataset(group,stt%gamma_sl,'dot_gamma_sl',&   ! this is not dot!!
-                                  'plastic shear','1')
-      case (Lambda_sl_ID)
-        call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',&
-                                  'mean free path for slip','m')
-      case (tau_pass_ID)
-        call results_writeDataset(group,dst%threshold_stress,'tau_pass',&
-                                  'threshold stress for slip','Pa')
+  outputsLoop: do o = 1,size(prm%output)
+    select case(trim(prm%output(o)))
+      case('edge_density')                                                                          ! ToDo: should be rho_mob
+        if(prm%sum_N_sl>0) call results_writeDataset(group,stt%rho_mob,'rho_mob',&
+                                                     'mobile dislocation density','1/m²')
+      case('dipole_density')                                                                        ! ToDo: should be rho_dip
+        if(prm%sum_N_sl>0) call results_writeDataset(group,stt%rho_dip,'rho_dip',&
+                                                     'dislocation dipole density''1/m²')
+      case('shear_rate_slip')                                                                       ! should be gamma
+        if(prm%sum_N_sl>0) call results_writeDataset(group,stt%gamma_sl,'dot_gamma_sl',&   ! this is not dot!!
+                                                     'plastic shear','1')
+      case('mfp_slip')                                                                              !ToDo: should be Lambda
+        if(prm%sum_N_sl>0) call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',&
+                                                     'mean free path for slip','m')
+      case('threshold_stress_slip')                                                                 !ToDo: should be tau_pass
+        if(prm%sum_N_sl>0) call results_writeDataset(group,dst%threshold_stress,'tau_pass',&
+                                                     'threshold stress for slip','Pa')
     end select
   enddo outputsLoop
   end associate
@@ -468,8 +427,8 @@ end subroutine plastic_disloUCLA_results
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Shear rates on slip systems, their derivatives with respect to resolved stress and the
-!  resolved stresss
+!> @brief Calculate shear rates on slip systems, their derivatives with respect to resolved
+!         stress, and the resolved stress.
 !> @details Derivatives and resolved stress are calculated only optionally.
 ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
 ! have the optional arguments at the end

src/constitutive_plastic_dislotwin.f90

@@ -9,27 +9,9 @@
 !--------------------------------------------------------------------------------------------------
 submodule(constitutive) plastic_dislotwin
 
-  real(pReal),                                    parameter :: &
+  real(pReal), parameter :: &
     kB = 1.38e-23_pReal                                                                             !< Boltzmann constant in J/Kelvin
 
-  enum, bind(c)
-    enumerator :: &
-      undefined_ID, &
-      rho_mob_ID, &
-      rho_dip_ID, &
-      dot_gamma_sl_ID, &
-      gamma_sl_ID, &
-      Lambda_sl_ID, &
-      resolved_stress_slip_ID, &
-      tau_pass_ID, &
-      edge_dipole_distance_ID, &
-      f_tw_ID, &
-      Lambda_tw_ID, &
-      resolved_stress_twin_ID, &
-      tau_hat_tw_ID, &
-      f_tr_ID
-  end enum
- 
   type :: tParameters
     real(pReal) :: &
       mu, &
@@ -59,7 +41,7 @@ submodule(constitutive) plastic_dislotwin
       gamma_fcc_hex, &                                                                              !< Free energy difference between austensite and martensite
       i_tr, &                                                                                       !<
       h                                                                                             !< Stack height of hex nucleus
-    real(pReal),                  dimension(:),     allocatable :: & 
+    real(pReal),               allocatable, dimension(:) :: &
       rho_mob_0, &                                                                                  !< initial unipolar dislocation density per slip system
       rho_dip_0, &                                                                                  !< initial dipole dislocation density per slip system
       b_sl, &                                                                                       !< absolute length of burgers vector [m] for each slip system
@@ -79,19 +61,16 @@ submodule(constitutive) plastic_dislotwin
       s, &                                                                                          !< s-exponent in trans nucleation rate
       gamma_char, &                                                                                 !< characteristic shear for twins
       B                                                                                             !< drag coefficient
-    real(pReal),                  dimension(:,:),   allocatable :: & 
+    real(pReal),               allocatable, dimension(:,:) :: &
       h_sl_sl, &                                                                                    !<
       h_sl_tw, &                                                                                    !<
       h_tw_tw, &                                                                                    !<
       h_sl_tr, &                                                                                    !<
-      h_tr_tr                                                                                       !< 
-    integer,                      dimension(:,:),   allocatable :: & 
-      fcc_twinNucleationSlipPair                                                                    ! ToDo: Better name? Is also use for trans
-    real(pReal),                  dimension(:,:),   allocatable :: & 
+      h_tr_tr, &                                                                                    !<
       n0_sl, &                                                                                      !< slip system normal
       forestProjection, &
       C66
-    real(pReal),                  dimension(:,:,:), allocatable :: &
+    real(pReal),               allocatable, dimension(:,:,:) :: &
       P_tr, &
       P_sl, &
       P_tw, &
@@ -101,12 +80,14 @@ submodule(constitutive) plastic_dislotwin
       sum_N_sl, &                                                                                   !< total number of active slip system
       sum_N_tw, &                                                                                   !< total number of active twin system
       sum_N_tr                                                                                      !< total number of active transformation system
-    integer,                      dimension(:),     allocatable :: & 
+    integer,                   allocatable, dimension(:) :: &
       N_sl, &                                                                                       !< number of active slip systems for each family
       N_tw, &                                                                                       !< number of active twin systems for each family
       N_tr                                                                                          !< number of active transformation systems for each family
-    integer(kind(undefined_ID)),  dimension(:),     allocatable :: &
-      outputID                                                                                      !< ID of each post result output
+    integer,                   allocatable, dimension(:,:) :: &
+      fcc_twinNucleationSlipPair                                                                    ! ToDo: Better name? Is also use for trans
+    character(len=pStringLen), allocatable, dimension(:) :: &
+      output
     logical :: &
       ExtendedDislocations, &                                                                       !< consider split into partials for climb calculation
       fccTwinTransNucleation, &                                                                     !< twinning and transformation models are for fcc
@@ -148,7 +129,7 @@ contains
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief module initialization
+!> @brief Perform module initialization.
 !> @details reads in material parameters, allocates arrays, and does sanity checks
 !--------------------------------------------------------------------------------------------------
 module subroutine plastic_dislotwin_init
@@ -156,28 +137,23 @@ module subroutine plastic_dislotwin_init
   integer :: &
     Ninstance, &
     p, i, &
-    NipcMyPhase, outputSize, &
+    NipcMyPhase, &
     sizeState, sizeDotState, &
     startIndex, endIndex
 
-  integer(kind(undefined_ID)) :: &
-    outputID
- 
   character(len=pStringLen) :: &
     extmsg = ''
-  character(len=pStringLen), dimension(:), allocatable :: &
-    outputs
 
-  write(6,'(/,a)')   ' <<<+-  constitutive_'//PLASTICITY_DISLOTWIN_label//' init  -+>>>'
+  write(6,'(/,a)') ' <<<+-  constitutive_'//PLASTICITY_DISLOTWIN_label//' init  -+>>>'; flush(6)
 
-  write(6,'(/,a)')   ' Ma and Roters, Acta Materialia 52(12):3603–3612, 2004'
-  write(6,'(a)')     ' https://doi.org/10.1016/j.actamat.2004.04.012'
+  write(6,'(/,a)') ' Ma and Roters, Acta Materialia 52(12):3603–3612, 2004'
+  write(6,'(a)')   ' https://doi.org/10.1016/j.actamat.2004.04.012'
 
-  write(6,'(/,a)')   ' Roters et al., Computational Materials Science 39:91–95, 2007'
-  write(6,'(a)')     ' https://doi.org/10.1016/j.commatsci.2006.04.014'
+  write(6,'(/,a)') ' Roters et al., Computational Materials Science 39:91–95, 2007'
+  write(6,'(a)')   ' https://doi.org/10.1016/j.commatsci.2006.04.014'
 
-  write(6,'(/,a)')   ' Wong et al., Acta Materialia 118:140–151, 2016'
-  write(6,'(a,/)')   ' https://doi.org/10.1016/j.actamat.2016.07.032'
+  write(6,'(/,a)') ' Wong et al., Acta Materialia 118:140–151, 2016'
+  write(6,'(a,/)') ' https://doi.org/10.1016/j.actamat.2016.07.032'
 
   Ninstance = count(phase_plasticity == PLASTICITY_DISLOTWIN_ID)
 
@@ -206,7 +182,6 @@ module subroutine plastic_dislotwin_init
     prm%nu = lattice_nu(p)
     prm%C66 = lattice_C66(1:6,1:6,p)
 
- 
 !--------------------------------------------------------------------------------------------------
 ! slip related parameters
     prm%N_sl      = config%getInts('nslip',defaultVal=emptyIntArray)
@@ -411,8 +386,6 @@ module subroutine plastic_dislotwin_init
       if (prm%q_sb          <= 0.0_pReal) extmsg = trim(extmsg)//' q_shearband'
     endif
 
- 
- 
     prm%D             = config%getFloat('grainsize')
 
     if (config%keyExists('dipoleformationfactor')) call IO_error(1,ext_msg='use /nodipoleformation/')
@@ -435,48 +408,7 @@ module subroutine plastic_dislotwin_init
         call IO_error(211,el=p,ext_msg='aTol_f_tr ('//PLASTICITY_DISLOTWIN_label//')')
     endif
 
-    outputs = config%getStrings('(output)', defaultVal=emptyStringArray)
-    allocate(prm%outputID(0))
-    do i= 1, size(outputs)
-      outputID = undefined_ID
-      select case(outputs(i))
-        case ('rho_mob')
-          outputID = merge(rho_mob_ID,undefined_ID,prm%sum_N_sl > 0)
-          outputSize = prm%sum_N_sl
-        case ('rho_dip')
-          outputID = merge(rho_dip_ID,undefined_ID,prm%sum_N_sl > 0)
-          outputSize = prm%sum_N_sl
-        case ('gamma_sl')
-          outputID = merge(gamma_sl_ID,undefined_ID,prm%sum_N_sl > 0)
-          outputSize = prm%sum_N_sl
-        case ('lambda_sl')
-          outputID = merge(Lambda_sl_ID,undefined_ID,prm%sum_N_sl > 0)
-          outputSize = prm%sum_N_sl
-        case ('tau_pass')
-          outputID= merge(tau_pass_ID,undefined_ID,prm%sum_N_sl > 0)
-          outputSize = prm%sum_N_sl
- 
-        case ('f_tw')
-          outputID = merge(f_tw_ID,undefined_ID,prm%sum_N_tw >0)
-          outputSize = prm%sum_N_tw
-        case ('lambda_tw')
-          outputID = merge(Lambda_tw_ID,undefined_ID,prm%sum_N_tw >0)
-          outputSize = prm%sum_N_tw
-        case ('tau_hat_tw')
-          outputID = merge(tau_hat_tw_ID,undefined_ID,prm%sum_N_tw >0)
-          outputSize = prm%sum_N_tw
-          
-        case ('f_tr')
-          outputID = f_tr_ID
-          outputSize = prm%sum_N_tr
-         
-      end select
-         
-      if (outputID /= undefined_ID) then
-        prm%outputID = [prm%outputID, outputID]
-      endif
- 
-    enddo
+    prm%output = config%getStrings('(output)', defaultVal=emptyStringArray)
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
@@ -548,7 +480,7 @@ end subroutine plastic_dislotwin_init
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief returns the homogenized elasticity matrix
+!> @brief Return the homogenized elasticity matrix.
 !--------------------------------------------------------------------------------------------------
 module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC)
 
@@ -587,7 +519,7 @@ end function plastic_dislotwin_homogenizedC
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief calculates plastic velocity gradient and its tangent
+!> @brief Calculate plastic velocity gradient and its tangent.
 !--------------------------------------------------------------------------------------------------
 module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
 
@@ -703,7 +635,7 @@ end subroutine plastic_dislotwin_LpAndItsTangent
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief calculates the rate of change of microstructure
+!> @brief Calculate the rate of change of microstructure.
 !--------------------------------------------------------------------------------------------------
 module subroutine plastic_dislotwin_dotState(Mp,T,instance,of)
 
@@ -806,7 +738,7 @@ end subroutine plastic_dislotwin_dotState
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief calculates derived quantities from state
+!> @brief Calculate derived quantities from state.
 !--------------------------------------------------------------------------------------------------
 module subroutine plastic_dislotwin_dependentState(T,instance,of)
 
@@ -898,47 +830,48 @@ end subroutine plastic_dislotwin_dependentState
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief writes results to HDF5 output file
+!> @brief Write results to HDF5 output file.
 !--------------------------------------------------------------------------------------------------
 module subroutine plastic_dislotwin_results(instance,group)
 
   integer,          intent(in) :: instance
   character(len=*), intent(in) :: group
+
   integer :: o
 
-  associate(prm => param(instance), stt => state(instance), dst => dependentState(instance))
-  outputsLoop: do o = 1,size(prm%outputID)
-    select case(prm%outputID(o))
+ associate(prm => param(instance), stt => state(instance), dst => dependentState(instance))
+  outputsLoop: do o = 1,size(prm%output)
+    select case(trim(prm%output(o)))
 
-      case (rho_mob_ID)
-        call results_writeDataset(group,stt%rho_mob,'rho_mob',&
-                                  'mobile dislocation density','1/m²')
-      case (rho_dip_ID)
-        call results_writeDataset(group,stt%rho_dip,'rho_dip',&
-                                  'dislocation dipole density''1/m²')
-      case (gamma_sl_ID)
-        call results_writeDataset(group,stt%gamma_sl,'gamma_sl',&
-                                  'plastic shear','1')
-      case (Lambda_sl_ID)
-        call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',&
-                                  'mean free path for slip','m')
-      case (tau_pass_ID)
-        call results_writeDataset(group,dst%tau_pass,'tau_pass',&
-                                  'passing stress for slip','Pa')
+      case('rho_mob')
+        if(prm%sum_N_sl>0) call results_writeDataset(group,stt%rho_mob,'rho_mob',&
+                                                     'mobile dislocation density','1/m²')
+      case('rho_dip')
+        if(prm%sum_N_sl>0) call results_writeDataset(group,stt%rho_dip,'rho_dip',&
+                                                     'dislocation dipole density''1/m²')
+      case('gamma_sl')
+        if(prm%sum_N_sl>0) call results_writeDataset(group,stt%gamma_sl,'gamma_sl',&
+                                                     'plastic shear','1')
+      case('lambda_sl')
+        if(prm%sum_N_sl>0) call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',&
+                                                     'mean free path for slip','m')
+      case('tau_pass')
+        if(prm%sum_N_sl>0) call results_writeDataset(group,dst%tau_pass,'tau_pass',&
+                                                     'passing stress for slip','Pa')
 
-      case (f_tw_ID)
-        call results_writeDataset(group,stt%f_tw,'f_tw',&
-                                 'twinned volume fraction','m³/m³')
-      case (Lambda_tw_ID)
-        call results_writeDataset(group,dst%Lambda_tw,'Lambda_tw',&
-                                  'mean free path for twinning','m')   
-      case (tau_hat_tw_ID)
-        call results_writeDataset(group,dst%tau_hat_tw,'tau_hat_tw',&
-                                  'threshold stress for twinning','Pa')
+      case('f_tw')
+        if(prm%sum_N_tw>0) call results_writeDataset(group,stt%f_tw,'f_tw',&
+                                                     'twinned volume fraction','m³/m³')
+      case('lambda_tw')
+        if(prm%sum_N_tw>0) call results_writeDataset(group,dst%Lambda_tw,'Lambda_tw',&
+                                                     'mean free path for twinning','m')
+      case('tau_hat_tw')
+        if(prm%sum_N_tw>0) call results_writeDataset(group,dst%tau_hat_tw,'tau_hat_tw',&
+                                                     'threshold stress for twinning','Pa')
 
-      case (f_tr_ID)
-        call results_writeDataset(group,stt%f_tr,'f_tr',&
-                                 'martensite volume fraction','m³/m³')
+      case('f_tr')
+        if(prm%sum_N_tr>0) call results_writeDataset(group,stt%f_tr,'f_tr',&
+                                                     'martensite volume fraction','m³/m³')
 
     end select
   enddo outputsLoop
@@ -948,8 +881,8 @@ end subroutine plastic_dislotwin_results
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Shear rates on slip systems, their derivatives with respect to resolved stress and the
-!  resolved stresss
+!> @brief Calculate shear rates on slip systems, their derivatives with respect to resolved
+!         stress, and the resolved stress.
 !> @details Derivatives and resolved stress are calculated only optionally.
 ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
 ! have the optional arguments at the end
@@ -1025,7 +958,11 @@ end subroutine kinetics_slip
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief calculates shear rates on twin systems
+!> @brief Calculate shear rates on twin systems and their derivatives with respect to resolved
+!         stress.
+!> @details Derivatives are calculated only optionally.
+! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
+! have the optional arguments at the end.
 !--------------------------------------------------------------------------------------------------
 pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,&
                               dot_gamma_twin,ddot_gamma_dtau_twin)
@@ -1090,7 +1027,11 @@ end subroutine kinetics_twin
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief calculates shear rates on twin systems
+!> @brief Calculate shear rates on transformation systems and their derivatives with respect to
+!         resolved stress.
+!> @details Derivatives are calculated only optionally.
+! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
+! have the optional arguments at the end.
 !--------------------------------------------------------------------------------------------------
 pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,&
                               dot_gamma_tr,ddot_gamma_dtau_trans)

src/constitutive_plastic_isotropic.f90

@@ -9,13 +9,6 @@
 !--------------------------------------------------------------------------------------------------
 submodule(constitutive) plastic_isotropic
 
-  enum, bind(c)
-    enumerator :: &
-      undefined_ID, &
-      xi_ID, &
-      dot_gamma_ID
-  end enum
- 
   type :: tParameters
     real(pReal) :: &
       M, &                                                                                          !< Taylor factor
@@ -34,10 +27,10 @@ submodule(constitutive) plastic_isotropic
       aTol_gamma
     integer :: &
       of_debug = 0
-    integer(kind(undefined_ID)), allocatable, dimension(:) :: &
-      outputID
     logical :: &
       dilatation
+    character(len=pStringLen), allocatable, dimension(:) :: &
+      output
   end type tParameters
 
   type :: tIsotropicState
@@ -56,29 +49,24 @@ submodule(constitutive) plastic_isotropic
 contains
 
 !--------------------------------------------------------------------------------------------------
-!> @brief module initialization
+!> @brief Perform module initialization.
 !> @details reads in material parameters, allocates arrays, and does sanity checks
 !--------------------------------------------------------------------------------------------------
 module subroutine plastic_isotropic_init
 
   integer :: &
     Ninstance, &
-    p, i, &
+    p, &
     NipcMyPhase, &
     sizeState, sizeDotState
 
-  integer(kind(undefined_ID)) :: &
-    outputID
- 
   character(len=pStringLen) :: &
     extmsg = ''
-  character(len=pStringLen), dimension(:), allocatable :: &
-    outputs
 
-  write(6,'(/,a)')   ' <<<+-  plastic_'//PLASTICITY_ISOTROPIC_label//' init  -+>>>'
+  write(6,'(/,a)') ' <<<+-  plastic_'//PLASTICITY_ISOTROPIC_label//' init  -+>>>'; flush(6)
 
-  write(6,'(/,a)')   ' Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018'
-  write(6,'(a)')     ' https://doi.org/10.1016/j.scriptamat.2017.09.047'
+  write(6,'(/,a)') ' Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018'
+  write(6,'(a)')   ' https://doi.org/10.1016/j.scriptamat.2017.09.047'
 
   Ninstance = count(phase_plasticity == PLASTICITY_ISOTROPIC_ID)
   if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
@@ -136,24 +124,7 @@ module subroutine plastic_isotropic_init
 
 !--------------------------------------------------------------------------------------------------
 !  output pararameters
-    outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
-    allocate(prm%outputID(0))
-    do i=1, size(outputs)
-      outputID = undefined_ID
-      select case(outputs(i))
- 
-        case ('flowstress')
-          outputID = xi_ID
-        case ('strainrate')
-          outputID = dot_gamma_ID
- 
-      end select
- 
-      if (outputID /= undefined_ID) then
-        prm%outputID = [prm%outputID, outputID]
-     endif
- 
-    enddo
+    prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
@@ -186,7 +157,7 @@ end subroutine plastic_isotropic_init
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief calculates plastic velocity gradient and its tangent
+!> @brief Calculate plastic velocity gradient and its tangent.
 !--------------------------------------------------------------------------------------------------
 module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
 
@@ -247,7 +218,7 @@ end subroutine plastic_isotropic_LpAndItsTangent
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief calculates plastic velocity gradient and its tangent
+!> @brief Calculate inelastic velocity gradient and its tangent.
 !--------------------------------------------------------------------------------------------------
 module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of)
 
@@ -299,7 +270,7 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of)
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief calculates the rate of change of microstructure
+!> @brief Calculate the rate of change of microstructure.
 !--------------------------------------------------------------------------------------------------
 module subroutine plastic_isotropic_dotState(Mp,instance,of)
 
@@ -348,19 +319,19 @@ end subroutine plastic_isotropic_dotState
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief writes results to HDF5 output file
+!> @brief Write results to HDF5 output file.
 !--------------------------------------------------------------------------------------------------
 module subroutine plastic_isotropic_results(instance,group)
 
-  integer, intent(in) :: instance
+  integer,          intent(in) :: instance
   character(len=*), intent(in) :: group
 
   integer :: o
 
   associate(prm => param(instance), stt => state(instance))
-  outputsLoop: do o = 1,size(prm%outputID)
-    select case(prm%outputID(o))
-      case (xi_ID)
+  outputsLoop: do o = 1,size(prm%output)
+    select case(trim(prm%output(o)))
+      case ('flowstress')                                                                           ! ToDo: should be 'xi'
         call results_writeDataset(group,stt%xi,'xi','resistance against plastic flow','Pa')
     end select
   enddo outputsLoop

src/constitutive_plastic_kinehardening.f90

@@ -7,26 +7,13 @@
 !--------------------------------------------------------------------------------------------------
 submodule(constitutive) plastic_kinehardening
 
-  enum, bind(c)
-    enumerator :: &
-      undefined_ID, &
-      crss_ID, &                                                                                    !< critical resolved stress
-      crss_back_ID, &                                                                               !< critical resolved back stress
-      sense_ID, &                                                                                   !< sense of acting shear stress (-1 or +1)
-      chi0_ID, &                                                                                    !< backstress at last switch of stress sense (positive?)
-      gamma0_ID, &                                                                                  !< accumulated shear at last switch of stress sense (at current switch?)
-      accshear_ID, &
-      shearrate_ID, &
-      resolvedstress_ID
-  end enum
- 
   type :: tParameters
     real(pReal) :: &
       gdot0, &                                                                                      !< reference shear strain rate for slip
       n, &                                                                                          !< stress exponent for slip
       aTolResistance, &
       aTolShear
-    real(pReal), allocatable, dimension(:) :: &
+    real(pReal),              allocatable, dimension(:) :: &
       crss0, &                                                                                      !< initial critical shear stress for slip
       theta0, &                                                                                     !< initial hardening rate of forward stress for each slip
       theta1, &                                                                                     !< asymptotic hardening rate of forward stress for each slip
@@ -35,19 +22,19 @@ submodule(constitutive) plastic_kinehardening
       tau1, &
       tau1_b, &
       nonSchmidCoeff
-    real(pReal), allocatable, dimension(:,:) :: &
+    real(pReal),              allocatable, dimension(:,:) :: &
       interaction_slipslip                                                                          !< slip resistance from slip activity
-    real(pReal), allocatable, dimension(:,:,:) :: &
+    real(pReal),              allocatable, dimension(:,:,:) :: &
       Schmid, &
       nonSchmid_pos, &
       nonSchmid_neg
     integer :: &
       totalNslip, &                                                                                 !< total number of active slip system
       of_debug = 0
-    integer,     allocatable, dimension(:) :: &
+    integer,                   allocatable, dimension(:) :: &
       Nslip                                                                                         !< number of active slip systems for each family
-    integer(kind(undefined_ID)), allocatable, dimension(:) :: &
-      outputID                                                                                      !< ID of each post result output
+    character(len=pStringLen), allocatable, dimension(:) :: &
+      output
   end type tParameters
 
   type :: tKinehardeningState
@@ -72,27 +59,22 @@ contains
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief module initialization
+!> @brief Perform module initialization.
 !> @details reads in material parameters, allocates arrays, and does sanity checks
 !--------------------------------------------------------------------------------------------------
 module subroutine plastic_kinehardening_init
 
   integer :: &
     Ninstance, &
-    p, i, o, &
+    p, o, &
     NipcMyPhase, &
     sizeState, sizeDeltaState, sizeDotState, &
     startIndex, endIndex
 
-  integer(kind(undefined_ID)) :: &
-    outputID
- 
   character(len=pStringLen) :: &
     extmsg = ''
-  character(len=pStringLen), dimension(:), allocatable :: &
-    outputs
 
-  write(6,'(/,a)')   ' <<<+-  plastic_'//PLASTICITY_KINEHARDENING_label//' init  -+>>>'
+  write(6,'(/,a)') ' <<<+-  plastic_'//PLASTICITY_KINEHARDENING_label//' init  -+>>>'; flush(6)
 
   Ninstance = count(phase_plasticity == PLASTICITY_KINEHARDENING_ID)
   if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
@@ -113,7 +95,7 @@ module subroutine plastic_kinehardening_init
 
 #ifdef DEBUG
     if  (p==material_phaseAt(debug_g,debug_e)) then
-       prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e)
+      prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e)
     endif
 #endif
 
@@ -155,8 +137,8 @@ module subroutine plastic_kinehardening_init
       prm%theta0_b = config%getFloats('theta0_b', requiredSize=size(prm%Nslip))
       prm%theta1_b = config%getFloats('theta1_b', requiredSize=size(prm%Nslip))
 
-      prm%gdot0  = config%getFloat('gdot0')
-      prm%n = config%getFloat('n_slip')
+      prm%gdot0    = config%getFloat('gdot0')
+      prm%n        = config%getFloat('n_slip')
 
       ! expand: family => system
       prm%crss0    = math_expand(prm%crss0,   prm%Nslip)
@@ -188,36 +170,7 @@ module subroutine plastic_kinehardening_init
 
 !--------------------------------------------------------------------------------------------------
 !  output pararameters
-    outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
-    allocate(prm%outputID(0))
-    do i=1, size(outputs)
-      outputID = undefined_ID
-      select case(outputs(i))
- 
-        case ('resistance')
-          outputID = merge(crss_ID,undefined_ID,prm%totalNslip>0)
-        case ('accumulatedshear')
-          outputID = merge(accshear_ID,undefined_ID,prm%totalNslip>0)
-        case ('shearrate')
-          outputID = merge(shearrate_ID,undefined_ID,prm%totalNslip>0)
-        case ('resolvedstress')
-          outputID = merge(resolvedstress_ID,undefined_ID,prm%totalNslip>0)
-        case ('backstress')
-          outputID = merge(crss_back_ID,undefined_ID,prm%totalNslip>0)
-        case ('sense')
-          outputID = merge(sense_ID,undefined_ID,prm%totalNslip>0)
-        case ('chi0')
-          outputID = merge(chi0_ID,undefined_ID,prm%totalNslip>0)
-        case ('gamma0')
-          outputID = merge(gamma0_ID,undefined_ID,prm%totalNslip>0)
- 
-      end select
- 
-      if (outputID /= undefined_ID) then
-        prm%outputID = [prm%outputID , outputID]
-      endif
- 
-    enddo
+    prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
@@ -277,7 +230,7 @@ end subroutine plastic_kinehardening_init
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief calculates plastic velocity gradient and its tangent
+!> @brief Calculate plastic velocity gradient and its tangent.
 !--------------------------------------------------------------------------------------------------
 pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
 
@@ -319,7 +272,7 @@ end subroutine plastic_kinehardening_LpAndItsTangent
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief calculates the rate of change of microstructure
+!> @brief Calculate the rate of change of microstructure.
 !--------------------------------------------------------------------------------------------------
 module subroutine plastic_kinehardening_dotState(Mp,instance,of)
 
@@ -361,7 +314,7 @@ end subroutine plastic_kinehardening_dotState
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief calculates (instantaneous) incremental change of microstructure
+!> @brief Calculate (instantaneous) incremental change of microstructure.
 !--------------------------------------------------------------------------------------------------
 module subroutine plastic_kinehardening_deltaState(Mp,instance,of)
 
@@ -409,38 +362,36 @@ end subroutine plastic_kinehardening_deltaState
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief writes results to HDF5 output file
+!> @brief Write results to HDF5 output file.
 !--------------------------------------------------------------------------------------------------
 module subroutine plastic_kinehardening_results(instance,group)
 
   integer,          intent(in) :: instance
   character(len=*), intent(in) :: group
+
   integer :: o
 
   associate(prm => param(instance), stt => state(instance))
-  outputsLoop: do o = 1,size(prm%outputID)
-    select case(prm%outputID(o))
-     case (crss_ID)
-       call results_writeDataset(group,stt%crss,'xi_sl', &
-                                'resistance against plastic slip','Pa')
-
-     case(crss_back_ID)
-       call results_writeDataset(group,stt%crss_back,'tau_back', &
-                                'back stress against plastic slip','Pa')
-                                
-     case (sense_ID)
-       call results_writeDataset(group,stt%sense,'sense_of_shear','tbd','1')
-
-     case (chi0_ID)
-       call results_writeDataset(group,stt%chi0,'chi0','tbd','Pa')
-       
-     case (gamma0_ID)
-       call results_writeDataset(group,stt%gamma0,'gamma0','tbd','1')
-       
-     case (accshear_ID)
-       call results_writeDataset(group,stt%accshear,'gamma_sl', &
-                                  'plastic shear','1')
-       
+  outputsLoop: do o = 1,size(prm%output)
+    select case(trim(prm%output(o)))
+     case('resistance')
+       if(prm%totalNslip>0) call results_writeDataset(group,stt%crss,'xi_sl', &
+                                                      'resistance against plastic slip','Pa')
+     case('backstress')                                                                             ! ToDo: should be 'tau_back'
+       if(prm%totalNslip>0) call results_writeDataset(group,stt%crss_back,'tau_back', &
+                                                      'back stress against plastic slip','Pa')
+     case ('sense')
+       if(prm%totalNslip>0) call results_writeDataset(group,stt%sense,'sense_of_shear', &
+                                                      'tbd','1')
+     case ('chi0')
+       if(prm%totalNslip>0) call results_writeDataset(group,stt%chi0,'chi0', &
+                                                      'tbd','Pa')
+     case ('gamma0')
+       if(prm%totalNslip>0) call results_writeDataset(group,stt%gamma0,'gamma0', &
+                                                      'tbd','1')
+     case ('accumulatedshear')
+       if(prm%totalNslip>0) call results_writeDataset(group,stt%accshear,'gamma_sl', &
+                                                      'plastic shear','1')
     end select
   enddo outputsLoop
   end associate
@@ -449,10 +400,11 @@ end subroutine plastic_kinehardening_results
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief calculates shear rates on slip systems and derivatives with respect to resolved stress
-!> @details: Shear rates are calculated only optionally.
+!> @brief Calculate shear rates on slip systems and their derivatives with respect to resolved
+!         stress.
+!> @details: Derivatives are calculated only optionally.
 ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
-! have the optional arguments at the end
+! have the optional arguments at the end.
 !--------------------------------------------------------------------------------------------------
 pure subroutine kinetics(Mp,instance,of, &
                          gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)

src/constitutive_plastic_none.f90

@@ -19,7 +19,7 @@ module subroutine plastic_none_init
     p, &
     NipcMyPhase
 
-  write(6,'(/,a)')   ' <<<+-  plastic_'//PLASTICITY_NONE_label//' init  -+>>>'
+  write(6,'(/,a)') ' <<<+-  plastic_'//PLASTICITY_NONE_label//' init  -+>>>'; flush(6)
 
   Ninstance = count(phase_plasticity == PLASTICITY_NONE_ID)
   if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &

src/constitutive_plastic_nonlocal.f90

@@ -48,32 +48,6 @@ submodule(constitutive) plastic_nonlocal
   real(pReal), dimension(:,:,:,:,:,:), allocatable :: &
     compatibility                                                                                   !< slip system compatibility between me and my neighbors
 
-  enum, bind(c)
-    enumerator :: &
-      undefined_ID, &
-      rho_sgl_mob_edg_pos_ID, &
-      rho_sgl_mob_edg_neg_ID, &
-      rho_sgl_mob_scr_pos_ID, &
-      rho_sgl_mob_scr_neg_ID, &
-      rho_sgl_imm_edg_pos_ID, &
-      rho_sgl_imm_edg_neg_ID, &
-      rho_sgl_imm_scr_pos_ID, &
-      rho_sgl_imm_scr_neg_ID, &
-      rho_dip_edg_ID, &
-      rho_dip_scr_ID, &
-      rho_forest_ID, &
-      resolvedstress_back_ID, &
-      tau_pass_ID, &
-      rho_dot_sgl_ID, &
-      rho_dot_sgl_mobile_ID, &
-      rho_dot_dip_ID, &
-      v_edg_pos_ID, &
-      v_edg_neg_ID, &
-      v_scr_pos_ID, &
-      v_scr_neg_ID, &
-      gamma_ID
-  end enum
-
   type :: tParameters                                                                               !< container type for internal constitutive parameters
     real(pReal) :: &
       atomicVolume, &                                                                               !< atomic volume
@@ -135,14 +109,12 @@ submodule(constitutive) plastic_nonlocal
     integer, dimension(:) ,allocatable :: &
       Nslip,&
       colinearSystem                                                                                !< colinear system to the active slip system (only valid for fcc!)
-
+    character(len=pStringLen), allocatable, dimension(:) :: &
+      output
     logical :: &
       shortRangeStressCorrection, &                                                                 !< flag indicating the use of the short range stress correction by a excess density gradient term
       probabilisticMultiplication
 
-    integer(kind(undefined_ID)),  dimension(:), allocatable :: &
-      outputID                                                                                      !< ID of each post result output
-
   end type tParameters
 
   type :: tNonlocalMicrostructure
@@ -198,22 +170,19 @@ module subroutine plastic_nonlocal_init
   integer :: &
     sizeState, sizeDotState,sizeDependentState, sizeDeltaState, &
     maxNinstances, &
-    p, i, &
+    p, &
     l, &
     s1, s2, &
     s, &
     t, &
     c
 
-  integer(kind(undefined_ID)) :: &
-    outputID
   character(len=pStringLen) :: &
     extmsg    = '', &
     structure
-  character(len=pStringLen), dimension(:), allocatable :: outputs
   integer :: NofMyPhase
 
-  write(6,'(/,a)') ' <<<+-  constitutive_'//PLASTICITY_NONLOCAL_label//' init  -+>>>'
+  write(6,'(/,a)') ' <<<+-  constitutive_'//PLASTICITY_NONLOCAL_label//' init  -+>>>'; flush(6)
 
   write(6,'(/,a)') ' Reuber et al., Acta Materialia 71:333–348, 2014'
   write(6,'(a)')   ' https://doi.org/10.1016/j.actamat.2014.03.012'
@@ -407,60 +376,7 @@ module subroutine plastic_nonlocal_init
 
     endif slipActive
 
-    outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
-    allocate(prm%outputID(0))
-    do i=1, size(outputs)
-      outputID = undefined_ID
-      select case(trim(outputs(i)))
-        case ('rho_sgl_mob_edg_pos')
-          outputID = merge(rho_sgl_mob_edg_pos_ID,undefined_ID,prm%totalNslip>0)
-        case ('rho_sgl_mob_edg_neg')
-          outputID = merge(rho_sgl_mob_edg_neg_ID,undefined_ID,prm%totalNslip>0)
-        case ('rho_sgl_mob_scr_pos')
-          outputID = merge(rho_sgl_mob_scr_pos_ID,undefined_ID,prm%totalNslip>0)
-        case ('rho_sgl_mob_scr_neg')
-          outputID = merge(rho_sgl_mob_scr_neg_ID,undefined_ID,prm%totalNslip>0)
-        case ('rho_sgl_imm_edg_pos')
-          outputID = merge(rho_sgl_imm_edg_pos_ID,undefined_ID,prm%totalNslip>0)
-        case ('rho_sgl_imm_edg_neg')
-          outputID = merge(rho_sgl_imm_edg_neg_ID,undefined_ID,prm%totalNslip>0)
-        case ('rho_sgl_imm_scr_pos')
-          outputID = merge(rho_sgl_imm_scr_pos_ID,undefined_ID,prm%totalNslip>0)
-        case ('rho_sgl_imm_scr_neg')
-          outputID = merge(rho_sgl_imm_scr_neg_ID,undefined_ID,prm%totalNslip>0)
-        case ('rho_dip_edg')
-          outputID = merge(rho_dip_edg_ID,undefined_ID,prm%totalNslip>0)
-        case ('rho_dip_scr')
-          outputID = merge(rho_dip_scr_ID,undefined_ID,prm%totalNslip>0)
-        case ('rho_forest')
-          outputID = merge(rho_forest_ID,undefined_ID,prm%totalNslip>0)
-        case ('resolvedstress_back')
-          outputID = merge(resolvedstress_back_ID,undefined_ID,prm%totalNslip>0)
-        case ('tau_pass')
-          outputID = merge(tau_pass_ID,undefined_ID,prm%totalNslip>0)
-        case ('rho_dot_sgl')
-          outputID = merge(rho_dot_sgl_ID,undefined_ID,prm%totalNslip>0)
-        case ('rho_dot_sgl_mobile')
-          outputID = merge(rho_dot_sgl_mobile_ID,undefined_ID,prm%totalNslip>0)
-        case ('rho_dot_dip')
-          outputID = merge(rho_dot_dip_ID,undefined_ID,prm%totalNslip>0)
-        case ('v_edg_pos')
-          outputID = merge(v_edg_pos_ID,undefined_ID,prm%totalNslip>0)
-        case ('v_edg_neg')
-          outputID = merge(v_edg_neg_ID,undefined_ID,prm%totalNslip>0)
-        case ('v_scr_pos')
-          outputID = merge(v_scr_pos_ID,undefined_ID,prm%totalNslip>0)
-        case ('v_scr_neg')
-          outputID = merge(v_scr_neg_ID,undefined_ID,prm%totalNslip>0)
-        case ('gamma')
-          outputID = merge(gamma_ID,undefined_ID,prm%totalNslip>0)
-      end select
-
-      if (outputID /= undefined_ID) then
-        prm%outputID = [prm%outputID , outputID]
-      endif
-
-    enddo
+    prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
@@ -561,7 +477,7 @@ module subroutine plastic_nonlocal_init
         stt%v_scr_neg  => plasticState(p)%state             (15*prm%totalNslip + 1:16*prm%totalNslip ,1:NofMyPhase)
 
     allocate(dst%tau_pass(prm%totalNslip,NofMyPhase),source=0.0_pReal)
-    allocate(dst%tau_Back(prm%totalNslip,NofMyPhase),     source=0.0_pReal)
+    allocate(dst%tau_back(prm%totalNslip,NofMyPhase),source=0.0_pReal)
     end associate
 
     if (NofMyPhase > 0) call stateInit(p,NofMyPhase)
@@ -1968,62 +1884,63 @@ module subroutine plastic_nonlocal_results(instance,group)
 
   integer,         intent(in) :: instance
   character(len=*),intent(in) :: group
+
   integer :: o
 
   associate(prm => param(instance),dst => microstructure(instance),stt=>state(instance))
-  outputsLoop: do o = 1,size(prm%outputID)
-    select case(prm%outputID(o))
-      case (rho_sgl_mob_edg_pos_ID)
-        call results_writeDataset(group,stt%rho_sgl_mob_edg_pos, 'rho_sgl_mob_edg_pos', &
-                                  'positive mobile edge density','1/m²')
-      case (rho_sgl_imm_edg_pos_ID)
-        call results_writeDataset(group,stt%rho_sgl_imm_edg_pos, 'rho_sgl_imm_edg_pos',&
-                                  'positive immobile edge density','1/m²')
-      case (rho_sgl_mob_edg_neg_ID)
-        call results_writeDataset(group,stt%rho_sgl_mob_edg_neg, 'rho_sgl_mob_edg_neg',&
-                                  'negative mobile edge density','1/m²')
-      case (rho_sgl_imm_edg_neg_ID)
-        call results_writeDataset(group,stt%rho_sgl_imm_edg_neg, 'rho_sgl_imm_edg_neg',&
-                                  'negative immobile edge density','1/m²')
-      case (rho_dip_edg_ID)
-        call results_writeDataset(group,stt%rho_dip_edg, 'rho_dip_edg',&
-                                  'edge dipole density','1/m²')
-      case (rho_sgl_mob_scr_pos_ID)
-        call results_writeDataset(group,stt%rho_sgl_mob_scr_pos, 'rho_sgl_mob_scr_pos',&
-                                  'positive mobile screw density','1/m²')
-      case (rho_sgl_imm_scr_pos_ID)
-        call results_writeDataset(group,stt%rho_sgl_imm_scr_pos, 'rho_sgl_imm_scr_pos',&
-                                  'positive immobile screw density','1/m²')
-      case (rho_sgl_mob_scr_neg_ID)
-        call results_writeDataset(group,stt%rho_sgl_mob_scr_neg, 'rho_sgl_mob_scr_neg',&
-                                  'negative mobile screw density','1/m²')
-      case (rho_sgl_imm_scr_neg_ID)
-        call results_writeDataset(group,stt%rho_sgl_imm_scr_neg, 'rho_sgl_imm_scr_neg',&
-                                  'negative immobile screw density','1/m²')
-      case (rho_dip_scr_ID)
-        call results_writeDataset(group,stt%rho_dip_scr, 'rho_dip_scr',&
-                                  'screw dipole density','1/m²')
-      case (rho_forest_ID)
-        call results_writeDataset(group,stt%rho_forest, 'rho_forest',&
-                                  'forest density','1/m²')
-      case (v_edg_pos_ID)
-        call results_writeDataset(group,stt%v_edg_pos, 'v_edg_pos',&
-                                  'positive edge velocity','m/s')
-      case (v_edg_neg_ID)
-        call results_writeDataset(group,stt%v_edg_neg, 'v_edg_neg',&
-                                  'negative edge velocity','m/s')
-      case (v_scr_pos_ID)
-        call results_writeDataset(group,stt%v_scr_pos, 'v_scr_pos',&
-                                  'positive srew velocity','m/s')
-      case (v_scr_neg_ID)
-        call results_writeDataset(group,stt%v_scr_neg, 'v_scr_neg',&
-                                  'negative screw velocity','m/s')
-      case(gamma_ID)
-        call results_writeDataset(group,stt%gamma,'gamma',&
-                                  'plastic shear','1')
-      case (tau_pass_ID)
-        call results_writeDataset(group,dst%tau_pass,'tau_pass',&
-                                  'passing stress for slip','Pa')
+  outputsLoop: do o = 1,size(prm%output)
+    select case(trim(prm%output(o)))
+      case('rho_sgl_mob_edg_pos')
+        if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_mob_edg_pos, 'rho_sgl_mob_edg_pos', &
+                                                       'positive mobile edge density','1/m²')
+      case('rho_sgl_imm_edg_pos')
+        if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_imm_edg_pos, 'rho_sgl_imm_edg_pos',&
+                                                       'positive immobile edge density','1/m²')
+      case('rho_sgl_mob_edg_neg')
+        if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_mob_edg_neg, 'rho_sgl_mob_edg_neg',&
+                                                       'negative mobile edge density','1/m²')
+      case('rho_sgl_imm_edg_neg')
+        if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_imm_edg_neg, 'rho_sgl_imm_edg_neg',&
+                                                       'negative immobile edge density','1/m²')
+      case('rho_dip_edg')
+        if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_dip_edg, 'rho_dip_edg',&
+                                                       'edge dipole density','1/m²')
+      case('rho_sgl_mob_scr_pos')
+        if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_mob_scr_pos, 'rho_sgl_mob_scr_pos',&
+                                                       'positive mobile screw density','1/m²')
+      case('rho_sgl_imm_scr_pos')
+        if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_imm_scr_pos, 'rho_sgl_imm_scr_pos',&
+                                                       'positive immobile screw density','1/m²')
+      case('rho_sgl_mob_scr_neg')
+        if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_mob_scr_neg, 'rho_sgl_mob_scr_neg',&
+                                                       'negative mobile screw density','1/m²')
+      case('rho_sgl_imm_scr_neg')
+        if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_sgl_imm_scr_neg, 'rho_sgl_imm_scr_neg',&
+                                                       'negative immobile screw density','1/m²')
+      case('rho_dip_scr')
+        if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_dip_scr, 'rho_dip_scr',&
+                                                       'screw dipole density','1/m²')
+      case('rho_forest')
+        if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_forest, 'rho_forest',&
+                                                       'forest density','1/m²')
+      case('v_edg_pos')
+        if(prm%totalNslip>0) call results_writeDataset(group,stt%v_edg_pos, 'v_edg_pos',&
+                                                       'positive edge velocity','m/s')
+      case('v_edg_neg')
+        if(prm%totalNslip>0) call results_writeDataset(group,stt%v_edg_neg, 'v_edg_neg',&
+                                                       'negative edge velocity','m/s')
+      case('v_scr_pos')
+        if(prm%totalNslip>0) call results_writeDataset(group,stt%v_scr_pos, 'v_scr_pos',&
+                                                       'positive srew velocity','m/s')
+      case('v_scr_neg')
+        if(prm%totalNslip>0) call results_writeDataset(group,stt%v_scr_neg, 'v_scr_neg',&
+                                                       'negative screw velocity','m/s')
+      case('gamma')
+        if(prm%totalNslip>0) call results_writeDataset(group,stt%gamma,'gamma',&
+                                                       'plastic shear','1')
+      case('tau_pass')
+        if(prm%totalNslip>0) call results_writeDataset(group,dst%tau_pass,'tau_pass',&
+                                                       'passing stress for slip','Pa')
     end select
   enddo outputsLoop
   end associate

src/constitutive_plastic_phenopowerlaw.f90

@@ -6,19 +6,6 @@
 !--------------------------------------------------------------------------------------------------
 submodule(constitutive) plastic_phenopowerlaw
 
-  enum, bind(c)
-    enumerator :: &
-      undefined_ID, &
-      resistance_slip_ID, &
-      accumulatedshear_slip_ID, &
-      shearrate_slip_ID, &
-      resolvedstress_slip_ID, &
-      resistance_twin_ID, &
-      accumulatedshear_twin_ID, &
-      shearrate_twin_ID, &
-      resolvedstress_twin_ID
-  end enum
- 
   type :: tParameters
     real(pReal) :: &
       gdot0_slip, &                                                                                 !< reference shear strain rate for slip
@@ -37,19 +24,19 @@ submodule(constitutive) plastic_phenopowerlaw
       aTolResistance, &                                                                             !< absolute tolerance for integration of xi
       aTolShear, &                                                                                  !< absolute tolerance for integration of gamma
       aTolTwinfrac                                                                                  !< absolute tolerance for integration of f
-    real(pReal), allocatable, dimension(:) :: &
+    real(pReal),               allocatable, dimension(:) :: &
       xi_slip_0, &                                                                                  !< initial critical shear stress for slip
       xi_twin_0, &                                                                                  !< initial critical shear stress for twin
       xi_slip_sat, &                                                                                !< maximum critical shear stress for slip
       nonSchmidCoeff, &
       H_int, &                                                                                      !< per family hardening activity (optional)
       gamma_twin_char                                                                               !< characteristic shear for twins
-    real(pReal), allocatable, dimension(:,:) :: &
+    real(pReal),               allocatable, dimension(:,:) :: &
       interaction_SlipSlip, &                                                                       !< slip resistance from slip activity
       interaction_SlipTwin, &                                                                       !< slip resistance from twin activity
       interaction_TwinSlip, &                                                                       !< twin resistance from slip activity
       interaction_TwinTwin                                                                          !< twin resistance from twin activity
-    real(pReal),  allocatable, dimension(:,:,:) :: &
+    real(pReal),               allocatable, dimension(:,:,:) :: &
       Schmid_slip, &
       Schmid_twin, &
       nonSchmid_pos, &
@@ -57,11 +44,11 @@ submodule(constitutive) plastic_phenopowerlaw
     integer :: &
       totalNslip, &                                                                                 !< total number of active slip system
       totalNtwin                                                                                    !< total number of active twin systems
-    integer,      allocatable, dimension(:) :: &
+    integer,                   allocatable, dimension(:) :: &
       Nslip, &                                                                                      !< number of active slip systems for each family
       Ntwin                                                                                         !< number of active twin systems for each family
-    integer(kind(undefined_ID)), allocatable, dimension(:) :: &
-      outputID                                                                                      !< ID of each post result output
+    character(len=pStringLen), allocatable, dimension(:) :: &
+      output
   end type tParameters
 
   type :: tPhenopowerlawState
@@ -83,7 +70,7 @@ contains
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief module initialization
+!> @brief Perform module initialization.
 !> @details reads in material parameters, allocates arrays, and does sanity checks
 !--------------------------------------------------------------------------------------------------
 module subroutine plastic_phenopowerlaw_init
@@ -91,19 +78,14 @@ module subroutine plastic_phenopowerlaw_init
   integer :: &
     Ninstance, &
     p, i, &
-    NipcMyPhase, outputSize, &
+    NipcMyPhase, &
     sizeState, sizeDotState, &
     startIndex, endIndex
 
-  integer(kind(undefined_ID)) :: &
-    outputID
- 
   character(len=pStringLen) :: &
     extmsg = ''
-  character(len=pStringLen), dimension(:), allocatable :: &
-    outputs
 
-  write(6,'(/,a)')   ' <<<+-  plastic_'//PLASTICITY_PHENOPOWERLAW_label//' init  -+>>>'
+  write(6,'(/,a)') ' <<<+-  plastic_'//PLASTICITY_PHENOPOWERLAW_label//' init  -+>>>'; flush(6)
 
   Ninstance = count(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)
   if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
@@ -239,45 +221,7 @@ module subroutine plastic_phenopowerlaw_init
 
 !--------------------------------------------------------------------------------------------------
 !  output pararameters
-    outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
-    allocate(prm%outputID(0))
-    do i=1, size(outputs)
-      outputID = undefined_ID
-      select case(outputs(i))
- 
-        case ('resistance_slip')
-          outputID = merge(resistance_slip_ID,undefined_ID,prm%totalNslip>0)
-          outputSize = prm%totalNslip
-        case ('accumulatedshear_slip')
-          outputID = merge(accumulatedshear_slip_ID,undefined_ID,prm%totalNslip>0)
-          outputSize = prm%totalNslip
-        case ('shearrate_slip')
-          outputID = merge(shearrate_slip_ID,undefined_ID,prm%totalNslip>0)
-          outputSize = prm%totalNslip
-        case ('resolvedstress_slip')
-          outputID = merge(resolvedstress_slip_ID,undefined_ID,prm%totalNslip>0)
-          outputSize = prm%totalNslip
- 
-        case ('resistance_twin')
-          outputID = merge(resistance_twin_ID,undefined_ID,prm%totalNtwin>0)
-          outputSize = prm%totalNtwin
-        case ('accumulatedshear_twin')
-          outputID = merge(accumulatedshear_twin_ID,undefined_ID,prm%totalNtwin>0)
-          outputSize = prm%totalNtwin
-        case ('shearrate_twin')
-          outputID = merge(shearrate_twin_ID,undefined_ID,prm%totalNtwin>0)
-          outputSize = prm%totalNtwin
-        case ('resolvedstress_twin')
-          outputID = merge(resolvedstress_twin_ID,undefined_ID,prm%totalNtwin>0)
-          outputSize = prm%totalNtwin
- 
-      end select
- 
-      if (outputID /= undefined_ID) then
-        prm%outputID = [prm%outputID, outputID]
-      endif
- 
-    enddo
+    prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
@@ -328,7 +272,7 @@ end subroutine plastic_phenopowerlaw_init
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief calculates plastic velocity gradient and its tangent
+!> @brief Calculate plastic velocity gradient and its tangent.
 !> @details asummes that deformation by dislocation glide affects twinned and untwinned volume
 !  equally (Taylor assumption). Twinning happens only in untwinned volume
 !--------------------------------------------------------------------------------------------------
@@ -381,7 +325,7 @@ end subroutine plastic_phenopowerlaw_LpAndItsTangent
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief calculates the rate of change of microstructure
+!> @brief Calculate the rate of change of microstructure.
 !--------------------------------------------------------------------------------------------------
 module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
 
@@ -437,7 +381,7 @@ end subroutine plastic_phenopowerlaw_dotState
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief writes results to HDF5 output file
+!> @brief Write results to HDF5 output file.
 !--------------------------------------------------------------------------------------------------
 module subroutine plastic_phenopowerlaw_results(instance,group)
 
@@ -447,22 +391,22 @@ module subroutine plastic_phenopowerlaw_results(instance,group)
   integer :: o
 
   associate(prm => param(instance), stt => state(instance))
-  outputsLoop: do o = 1,size(prm%outputID)
-    select case(prm%outputID(o))
+  outputsLoop: do o = 1,size(prm%output)
+    select case(trim(prm%output(o)))
 
-      case (resistance_slip_ID)
-        call results_writeDataset(group,stt%xi_slip,   'xi_sl', &
-                                  'resistance against plastic slip','Pa')
-      case (accumulatedshear_slip_ID)
-        call results_writeDataset(group,stt%gamma_slip,'gamma_sl', &
-                                  'plastic shear','1')
+      case('resistance_slip')
+        if(prm%totalNslip>0) call results_writeDataset(group,stt%xi_slip,   'xi_sl', &
+                                                       'resistance against plastic slip','Pa')
+      case('accumulatedshear_slip')
+        if(prm%totalNslip>0) call results_writeDataset(group,stt%gamma_slip,'gamma_sl', &
+                                                       'plastic shear','1')
 
-      case (resistance_twin_ID)
-        call results_writeDataset(group,stt%xi_twin,   'xi_tw', &
-                                  'resistance against twinning','Pa')
-      case (accumulatedshear_twin_ID)
-        call results_writeDataset(group,stt%gamma_twin,'gamma_tw', &
-                                  'twinning shear','1')
+      case('resistance_twin')
+        if(prm%totalNtwin>0) call results_writeDataset(group,stt%xi_twin,   'xi_tw', &
+                                                       'resistance against twinning','Pa')
+      case('accumulatedshear_twin')
+        if(prm%totalNtwin>0) call results_writeDataset(group,stt%gamma_twin,'gamma_tw', &
+                                                       'twinning shear','1')
 
     end select
   enddo outputsLoop
@@ -472,10 +416,11 @@ end subroutine plastic_phenopowerlaw_results
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Shear rates on slip systems and their derivatives with respect to resolved stress
+!> @brief Calculate shear rates on slip systems and their derivatives with respect to resolved
+!         stress.
 !> @details Derivatives are calculated only optionally.
 ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
-! have the optional arguments at the end
+! have the optional arguments at the end.
 !--------------------------------------------------------------------------------------------------
 pure subroutine kinetics_slip(Mp,instance,of, &
                               gdot_slip_pos,gdot_slip_neg,dgdot_dtau_slip_pos,dgdot_dtau_slip_neg)
@@ -543,9 +488,9 @@ end subroutine kinetics_slip
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Shear rates on twin systems and their derivatives with respect to resolved stress.
-!  twinning is assumed to take place only in untwinned volume.
-!> @details Derivates are calculated only optionally.
+!> @brief Calculate shear rates on twin systems and their derivatives with respect to resolved
+!         stress. Twinning is assumed to take place only in untwinned volume.
+!> @details Derivatives are calculated only optionally.
 ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
 ! have the optional arguments at the end.
 !--------------------------------------------------------------------------------------------------

src/crystallite.f90

@@ -57,9 +57,7 @@ module crystallite
     crystallite_Li0, &                                                                              !< intermediate velocitiy grad at start of FE inc
     crystallite_partionedLi0                                                                        !< intermediate velocity grad at start of homog inc
   real(pReal),                dimension(:,:,:,:,:),    allocatable :: &
-    crystallite_invFp, &                                                                            !< inverse of current plastic def grad (end of converged time step)
     crystallite_subFp0,&                                                                            !< plastic def grad at start of crystallite inc
-    crystallite_invFi, &                                                                            !< inverse of current intermediate def grad (end of converged time step)
     crystallite_subFi0,&                                                                            !< intermediate def grad at start of crystallite inc
     crystallite_subF,  &                                                                            !< def grad to be reached at end of crystallite inc
     crystallite_subF0, &                                                                            !< def grad at start of crystallite inc
@@ -145,12 +143,10 @@ subroutine crystallite_init
   allocate(crystallite_partionedFp0(3,3,cMax,iMax,eMax),      source=0.0_pReal)
   allocate(crystallite_subFp0(3,3,cMax,iMax,eMax),            source=0.0_pReal)
   allocate(crystallite_Fp(3,3,cMax,iMax,eMax),                source=0.0_pReal)
-  allocate(crystallite_invFp(3,3,cMax,iMax,eMax),             source=0.0_pReal)
   allocate(crystallite_Fi0(3,3,cMax,iMax,eMax),               source=0.0_pReal)
   allocate(crystallite_partionedFi0(3,3,cMax,iMax,eMax),      source=0.0_pReal)
   allocate(crystallite_subFi0(3,3,cMax,iMax,eMax),            source=0.0_pReal)
   allocate(crystallite_Fi(3,3,cMax,iMax,eMax),                source=0.0_pReal)
-  allocate(crystallite_invFi(3,3,cMax,iMax,eMax),             source=0.0_pReal)
   allocate(crystallite_Fe(3,3,cMax,iMax,eMax),                source=0.0_pReal)
   allocate(crystallite_Lp0(3,3,cMax,iMax,eMax),               source=0.0_pReal)
   allocate(crystallite_partionedLp0(3,3,cMax,iMax,eMax),      source=0.0_pReal)
@@ -408,9 +404,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
           else
             crystallite_subStep(c,i,e)       = num%subStepSizeCryst * crystallite_subStep(c,i,e)
             crystallite_Fp   (1:3,1:3,c,i,e) =            crystallite_subFp0(1:3,1:3,c,i,e)
-            crystallite_invFp(1:3,1:3,c,i,e) = math_inv33(crystallite_Fp    (1:3,1:3,c,i,e))
             crystallite_Fi   (1:3,1:3,c,i,e) =            crystallite_subFi0(1:3,1:3,c,i,e)
-            crystallite_invFi(1:3,1:3,c,i,e) = math_inv33(crystallite_Fi    (1:3,1:3,c,i,e))
             crystallite_S    (1:3,1:3,c,i,e) =            crystallite_S0    (1:3,1:3,c,i,e)
             if (crystallite_subStep(c,i,e) < 1.0_pReal) then                                        ! actual (not initial) cutback
               crystallite_Lp (1:3,1:3,c,i,e) =            crystallite_subLp0(1:3,1:3,c,i,e)
@@ -431,11 +425,11 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
 !  prepare for integration
           if (crystallite_todo(c,i,e)) then
             crystallite_subF(1:3,1:3,c,i,e) = crystallite_subF0(1:3,1:3,c,i,e) &
-                                            + crystallite_subStep(c,i,e) * (crystallite_partionedF (1:3,1:3,c,i,e) &
-                                                                          - crystallite_partionedF0(1:3,1:3,c,i,e))
-            crystallite_Fe(1:3,1:3,c,i,e) = matmul(matmul(crystallite_subF (1:3,1:3,c,i,e), &
-                                                          crystallite_invFp(1:3,1:3,c,i,e)), &
-                                                   crystallite_invFi(1:3,1:3,c,i,e))
+                                            + crystallite_subStep(c,i,e) *( crystallite_partionedF (1:3,1:3,c,i,e) &
+                                                                           -crystallite_partionedF0(1:3,1:3,c,i,e))
+            crystallite_Fe(1:3,1:3,c,i,e) = matmul(matmul(crystallite_subF(1:3,1:3,c,i,e), &
+                                                          math_inv33(crystallite_Fp(1:3,1:3,c,i,e))), &
+                                                   math_inv33(crystallite_Fi(1:3,1:3,c,i,e)))
             crystallite_subdt(c,i,e) = crystallite_subStep(c,i,e) * crystallite_dt(c,i,e)
             crystallite_converged(c,i,e) = .false.
           endif
@@ -477,7 +471,9 @@ subroutine crystallite_stressTangent
     o, &
     p
 
-  real(pReal), dimension(3,3)     ::   temp_33_1, devNull,invSubFi0, temp_33_2, temp_33_3, temp_33_4
+  real(pReal), dimension(3,3)     ::   devNull, &
+                                       invSubFp0,invSubFi0,invFp,invFi, &
+                                       temp_33_1, temp_33_2, temp_33_3, temp_33_4
   real(pReal), dimension(3,3,3,3) ::   dSdFe, &
                                        dSdF, &
                                        dSdFi, &
@@ -493,7 +489,8 @@ subroutine crystallite_stressTangent
   real(pReal), dimension(9,9)::        temp_99
   logical :: error
 
-  !$OMP PARALLEL DO PRIVATE(dSdF,dSdFe,dSdFi,dLpdS,dLpdFi,dFpinvdF,dLidS,dLidFi,dFidS,invSubFi0,o,p, &
+  !$OMP PARALLEL DO PRIVATE(dSdF,dSdFe,dSdFi,dLpdS,dLpdFi,dFpinvdF,dLidS,dLidFi,dFidS,o,p, &
+  !$OMP                     invSubFp0,invSubFi0,invFp,invFi, &
   !$OMP                     rhs_3333,lhs_3333,temp_99,temp_33_1,temp_33_2,temp_33_3,temp_33_4,temp_3333,error)
   elementLooping: do e = FEsolving_execElem(1),FEsolving_execElem(2)
     do i = FEsolving_execIP(1),FEsolving_execIP(2)
@@ -507,16 +504,20 @@ subroutine crystallite_stressTangent
                                            crystallite_Fi(1:3,1:3,c,i,e), &
                                            c,i,e)
 
+        invFp = math_inv33(crystallite_Fp(1:3,1:3,c,i,e))
+        invFi = math_inv33(crystallite_Fi(1:3,1:3,c,i,e))
+        invSubFp0 = math_inv33(crystallite_subFp0(1:3,1:3,c,i,e))
+        invSubFi0 = math_inv33(crystallite_subFi0(1:3,1:3,c,i,e))
+
         if (sum(abs(dLidS)) < tol_math_check) then
           dFidS = 0.0_pReal
         else
-          invSubFi0 = math_inv33(crystallite_subFi0(1:3,1:3,c,i,e))
           lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal
           do o=1,3; do p=1,3
             lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) &
                                   + crystallite_subdt(c,i,e)*matmul(invSubFi0,dLidFi(1:3,1:3,o,p))
             lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) &
-                                  + crystallite_invFi(1:3,1:3,c,i,e)*crystallite_invFi(p,o,c,i,e)
+                                  + invFi*invFi(p,o)
             rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) &
                                   - crystallite_subdt(c,i,e)*matmul(invSubFi0,dLidS(1:3,1:3,o,p))
           enddo; enddo
@@ -538,18 +539,13 @@ subroutine crystallite_stressTangent
 
 !--------------------------------------------------------------------------------------------------
 ! calculate dSdF
-        temp_33_1 = transpose(matmul(crystallite_invFp(1:3,1:3,c,i,e), &
-                                     crystallite_invFi(1:3,1:3,c,i,e)))
-        temp_33_2 = matmul(crystallite_subF(1:3,1:3,c,i,e), &
-                           math_inv33(crystallite_subFp0(1:3,1:3,c,i,e)))
-        temp_33_3 = matmul(matmul(crystallite_subF(1:3,1:3,c,i,e), &
-                                  crystallite_invFp (1:3,1:3,c,i,e)), &
-                           math_inv33(crystallite_subFi0(1:3,1:3,c,i,e)))
+        temp_33_1 = transpose(matmul(invFp,invFi))
+        temp_33_2 = matmul(crystallite_subF(1:3,1:3,c,i,e),invSubFp0)
+        temp_33_3 = matmul(matmul(crystallite_subF(1:3,1:3,c,i,e),invFp), invSubFi0)
 
         do o=1,3; do p=1,3
           rhs_3333(p,o,1:3,1:3)  = matmul(dSdFe(p,o,1:3,1:3),temp_33_1)
-          temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), &
-                                          crystallite_invFi(1:3,1:3,c,i,e)) &
+          temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), invFi) &
                                  + matmul(temp_33_3,dLidS(1:3,1:3,p,o))
         enddo; enddo
         lhs_3333 = crystallite_subdt(c,i,e)*math_mul3333xx3333(dSdFe,temp_3333) &
@@ -569,15 +565,14 @@ subroutine crystallite_stressTangent
         temp_3333 = math_mul3333xx3333(dLpdS,dSdF)
         do o=1,3; do p=1,3
           dFpinvdF(1:3,1:3,p,o) = -crystallite_subdt(c,i,e) &
-                                * matmul(math_inv33(crystallite_subFp0(1:3,1:3,c,i,e)), &
-                                         matmul(temp_3333(1:3,1:3,p,o),crystallite_invFi(1:3,1:3,c,i,e)))
+                                * matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi))
         enddo; enddo
 
 !--------------------------------------------------------------------------------------------------
 ! assemble dPdF
-        temp_33_1 = matmul(crystallite_S(1:3,1:3,c,i,e),transpose(crystallite_invFp(1:3,1:3,c,i,e)))
-        temp_33_2 = matmul(crystallite_invFp(1:3,1:3,c,i,e),temp_33_1)
-        temp_33_3 = matmul(crystallite_subF(1:3,1:3,c,i,e),crystallite_invFp(1:3,1:3,c,i,e))
+        temp_33_1 = matmul(crystallite_S(1:3,1:3,c,i,e),transpose(invFp))
+        temp_33_2 = matmul(invFp,temp_33_1)
+        temp_33_3 = matmul(crystallite_subF(1:3,1:3,c,i,e),invFp)
         temp_33_4 = matmul(temp_33_3,crystallite_S(1:3,1:3,c,i,e))
 
         crystallite_dPdF(1:3,1:3,1:3,1:3,c,i,e) = 0.0_pReal
@@ -589,7 +584,7 @@ subroutine crystallite_stressTangent
                                               + matmul(matmul(crystallite_subF(1:3,1:3,c,i,e), &
                                                        dFpinvdF(1:3,1:3,p,o)),temp_33_1) &
                                               + matmul(matmul(temp_33_3,dSdF(1:3,1:3,p,o)), &
-                                                       transpose(crystallite_invFp(1:3,1:3,c,i,e))) &
+                                                       transpose(invFp)) &
                                               + matmul(temp_33_4,transpose(dFpinvdF(1:3,1:3,p,o)))
         enddo; enddo
 
@@ -734,7 +729,7 @@ subroutine crystallite_results
     real(pReal), allocatable, dimension(:,:,:) :: select_tensors
     integer :: e,i,c,j
 
-    allocate(select_tensors(3,3,count(material_phaseAt==instance)*homogenization_maxNgrains*discretization_nIP))
+    allocate(select_tensors(3,3,count(material_phaseAt==instance)*discretization_nIP))
 
     j=0
     do e = 1, size(material_phaseAt,2)
@@ -793,28 +788,28 @@ logical function integrateStress(ipc,ip,el,timeFraction)
 
   real(pReal), dimension(3,3)::       F, &                                                          ! deformation gradient at end of timestep
                                       Fp_new, &                                                     ! plastic deformation gradient at end of timestep
-                                      Fe_new, &                                                     ! elastic deformation gradient at end of timestep
                                       invFp_new, &                                                  ! inverse of Fp_new
-                                      Fi_new, &                                                     ! gradient of intermediate deformation stages
-                                      invFi_new, &
                                       invFp_current, &                                              ! inverse of Fp_current
-                                      invFi_current, &                                              ! inverse of Fp_current
                                       Lpguess, &                                                    ! current guess for plastic velocity gradient
                                       Lpguess_old, &                                                ! known last good guess for plastic velocity gradient
                                       Lp_constitutive, &                                            ! plastic velocity gradient resulting from constitutive law
                                       residuumLp, &                                                 ! current residuum of plastic velocity gradient
                                       residuumLp_old, &                                             ! last residuum of plastic velocity gradient
                                       deltaLp, &                                                    ! direction of next guess
+                                      Fi_new, &                                                     ! gradient of intermediate deformation stages
+                                      invFi_new, &
+                                      invFi_current, &                                              ! inverse of Fi_current
                                       Liguess, &                                                    ! current guess for intermediate velocity gradient
                                       Liguess_old, &                                                ! known last good guess for intermediate velocity gradient
                                       Li_constitutive, &                                            ! intermediate velocity gradient resulting from constitutive law
                                       residuumLi, &                                                 ! current residuum of intermediate velocity gradient
                                       residuumLi_old, &                                             ! last residuum of intermediate velocity gradient
                                       deltaLi, &                                                    ! direction of next guess
+                                      Fe, &                                                         ! elastic deformation gradient
+                                      Fe_new, &
                                       S, &                                                          ! 2nd Piola-Kirchhoff Stress in plastic (lattice) configuration
                                       A, &
                                       B, &
-                                      Fe, &                                                         ! elastic deformation gradient
                                       temp_33
   real(pReal), dimension(9) ::        temp_9                                                        ! needed for matrix inversion by LAPACK
   integer,     dimension(9) ::        devNull_9                                                     ! needed for matrix inversion by LAPACK
@@ -993,16 +988,13 @@ logical function integrateStress(ipc,ip,el,timeFraction)
   Fe_new = matmul(matmul(F,invFp_new),invFi_new)
 
   integrateStress = .true.
-  crystallite_P    (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new))) ! ToDo: We propably do not need to store P!
+  crystallite_P    (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new)))
   crystallite_S    (1:3,1:3,ipc,ip,el) = S
   crystallite_Lp   (1:3,1:3,ipc,ip,el) = Lpguess
   crystallite_Li   (1:3,1:3,ipc,ip,el) = Liguess
   crystallite_Fp   (1:3,1:3,ipc,ip,el) = Fp_new
   crystallite_Fi   (1:3,1:3,ipc,ip,el) = Fi_new
   crystallite_Fe   (1:3,1:3,ipc,ip,el) = Fe_new
-  crystallite_invFp(1:3,1:3,ipc,ip,el) = invFp_new
-  crystallite_invFi(1:3,1:3,ipc,ip,el) = invFi_new
-
 
 end function integrateStress
 

src/grid/grid_mech_FEM.f90

@@ -475,7 +475,7 @@ subroutine formResidual(da_local,x_local, &
     write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter+1, '≤', itmax
     if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
       write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
-              ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotTensor2(F_aim,active=.true.))
+              ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
     write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
               ' deformation gradient aim       =', transpose(F_aim)
     flush(6)
@@ -491,7 +491,7 @@ subroutine formResidual(da_local,x_local, &
       x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk)
     enddo; enddo; enddo
     ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1
-    F(1:3,1:3,ii,jj,kk) = params%rotation_BC%rotTensor2(F_aim,active=.true.) + transpose(matmul(BMat,x_elem)) 
+    F(1:3,1:3,ii,jj,kk) = params%rotation_BC%rotate(F_aim,active=.true.) + transpose(matmul(BMat,x_elem))
   enddo; enddo; enddo
   call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr)
 

src/grid/grid_mech_spectral_basic.f90

@@ -297,7 +297,7 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo
 
     Fdot = utilities_calculateRate(guess, &
                                    F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, &
-                                   rotation_BC%rotTensor2(F_aimDot,active=.true.))
+                                   rotation_BC%rotate(F_aimDot,active=.true.))
     F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3])
 
     materialpoint_F0 = reshape(F, [3,3,1,product(grid(1:2))*grid3])
@@ -307,7 +307,7 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo
 ! update average and local deformation gradients
   F_aim = F_aim_lastInc + F_aimDot * timeinc
   F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, &                                       ! estimate of F at end of time+timeinc that matches rotated F_aim on average
-              rotation_BC%rotTensor2(F_aim,active=.true.)),[9,grid(1),grid(2),grid3])
+              rotation_BC%rotate(F_aim,active=.true.)),[9,grid(1),grid(2),grid3])
   call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
 
 end subroutine grid_mech_spectral_basic_forward
@@ -441,7 +441,7 @@ subroutine formResidual(in, F, &
     write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
     if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
       write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
-              ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotTensor2(F_aim,active=.true.))
+              ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
     write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
               ' deformation gradient aim       =', transpose(F_aim)
     flush(6)
@@ -466,7 +466,7 @@ subroutine formResidual(in, F, &
   tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residuum                                  ! store fPK field for subsequent FFT forward transform
   call utilities_FFTtensorForward                                                                   ! FFT forward of global "tensorField_real"
   err_div = Utilities_divergenceRMS()                                                               ! divRMS of tensorField_fourier for later use
-  call utilities_fourierGammaConvolution(params%rotation_BC%rotTensor2(deltaF_aim,active=.true.))   ! convolution of Gamma and tensorField_fourier
+  call utilities_fourierGammaConvolution(params%rotation_BC%rotate(deltaF_aim,active=.true.))       ! convolution of Gamma and tensorField_fourier
   call utilities_FFTtensorBackward                                                                  ! FFT backward of global tensorField_fourier
 
 !--------------------------------------------------------------------------------------------------

src/grid/grid_mech_spectral_polarisation.f90

@@ -321,10 +321,10 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
 
     Fdot     = utilities_calculateRate(guess, &
                                        F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, &
-                                       rotation_BC%rotTensor2(F_aimDot,active=.true.))
+                                       rotation_BC%rotate(F_aimDot,active=.true.))
     F_tauDot = utilities_calculateRate(guess, &
                                        F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid3]), timeinc_old, &
-                                       rotation_BC%rotTensor2(F_aimDot,active=.true.))
+                                       rotation_BC%rotate(F_aimDot,active=.true.))
     F_lastInc     = reshape(F,    [3,3,grid(1),grid(2),grid3])
     F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3])
 
@@ -335,7 +335,7 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
 ! update average and local deformation gradients
   F_aim = F_aim_lastInc + F_aimDot * timeinc
   F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, &                                      ! estimate of F at end of time+timeinc that matches rotated F_aim on average
-                                     rotation_BC%rotTensor2(F_aim,active=.true.)),&
+                                     rotation_BC%rotate(F_aim,active=.true.)),&
               [9,grid(1),grid(2),grid3])
   if (guess) then
      F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), &
@@ -510,7 +510,7 @@ subroutine formResidual(in, FandF_tau, &
     write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
     if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
       write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
-              ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotTensor2(F_aim,active=.true.))
+              ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
     write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
               ' deformation gradient aim       =', transpose(F_aim)
     flush(6)
@@ -529,7 +529,7 @@ subroutine formResidual(in, FandF_tau, &
 !--------------------------------------------------------------------------------------------------
 ! doing convolution in Fourier space
   call utilities_FFTtensorForward
-  call utilities_fourierGammaConvolution(params%rotation_BC%rotTensor2(polarBeta*F_aim,active=.true.)) 
+  call utilities_fourierGammaConvolution(params%rotation_BC%rotate(polarBeta*F_aim,active=.true.))
   call utilities_FFTtensorBackward
 
 !--------------------------------------------------------------------------------------------------
@@ -547,7 +547,7 @@ subroutine formResidual(in, FandF_tau, &
 ! stress BC handling
   F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC)))                                  ! S = 0.0 for no bc
   err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim &
-                                                 -params%rotation_BC%rotTensor2(F_av)) + &
+                                                 -params%rotation_BC%rotate(F_av)) + &
                                  params%stress_mask  * (P_av-params%stress_BC)))                    ! mask = 0.0 for no bc
 ! calculate divergence
   tensorField_real = 0.0_pReal

src/grid/spectral_utilities.f90

@@ -688,9 +688,10 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
   type(rotation), intent(in)                     :: rot_BC                                          !< rotation of load frame
   logical,        intent(in), dimension(3,3)     :: mask_stress                                     !< mask of stress BC
 
-  integer :: j, k, m, n
-  logical, dimension(9) :: mask_stressVector
-  real(pReal), dimension(9,9) :: temp99_Real
+  integer :: i, j
+  logical, dimension(9)   :: mask_stressVector
+  logical, dimension(9,9) :: mask
+  real(pReal), dimension(9,9) :: temp99_real
   integer :: size_reduced = 0
   real(pReal),              dimension(:,:), allocatable ::  &
     s_reduced, &                                                                                    !< reduced compliance matrix (depending on number of stress BC)
@@ -701,11 +702,8 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
 
   mask_stressVector = reshape(transpose(mask_stress), [9])
   size_reduced = count(mask_stressVector)
-  if(size_reduced > 0 )then
-    allocate (c_reduced(size_reduced,size_reduced), source =0.0_pReal)
-    allocate (s_reduced(size_reduced,size_reduced), source =0.0_pReal)
-    allocate (sTimesC(size_reduced,size_reduced),   source =0.0_pReal)
-    temp99_Real = math_3333to99(rot_BC%rotTensor4(C))
+  if(size_reduced > 0) then
+    temp99_real = math_3333to99(rot_BC%rotate(C))
 
     if(debugGeneral) then
       write(6,'(/,a)') ' ... updating masked compliance ............................................'
@@ -713,42 +711,21 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
                                                    transpose(temp99_Real)*1.0e-9_pReal
       flush(6)
     endif
-    k = 0                                                                                           ! calculate reduced stiffness
-    do n = 1,9
-      if(mask_stressVector(n)) then
-        k = k + 1
-        j = 0
-        do m = 1,9
-          if(mask_stressVector(m)) then
-            j = j + 1
-            c_reduced(k,j) = temp99_Real(n,m)
-    endif; enddo; endif; enddo
 
+    do i = 1,9; do j = 1,9
+      mask(i,j) = mask_stressVector(i) .and. mask_stressVector(j)
+    enddo; enddo
+    c_reduced = reshape(pack(temp99_Real,mask),[size_reduced,size_reduced])
+
+    allocate(s_reduced,mold = c_reduced)
     call math_invert(s_reduced, errmatinv, c_reduced)                                               ! invert reduced stiffness
     if (any(IEEE_is_NaN(s_reduced))) errmatinv = .true.
     if (errmatinv) call IO_error(error_ID=400,ext_msg='utilities_maskedCompliance')
-    temp99_Real = 0.0_pReal                                                                         ! fill up compliance with zeros
-     k = 0
-     do n = 1,9
-       if(mask_stressVector(n)) then
-         k = k + 1
-         j = 0
-         do m = 1,9
-           if(mask_stressVector(m)) then
-             j = j + 1
-             temp99_Real(n,m) = s_reduced(k,j)
-    endif; enddo; endif; enddo
 
 !--------------------------------------------------------------------------------------------------
 ! check if inversion was successful
     sTimesC = matmul(c_reduced,s_reduced)
-    do m=1, size_reduced
-      do n=1, size_reduced
-        errmatinv = errmatinv &
-               .or. (m==n .and. abs(sTimesC(m,n)-1.0_pReal) > 1.0e-12_pReal) &                      ! diagonal elements of S*C should be 1
-               .or. (m/=n .and. abs(sTimesC(m,n))           > 1.0e-12_pReal)                        ! off-diagonal elements of S*C should be 0
-      enddo
-    enddo
+    errmatinv = errmatinv .or. any(dNeq(sTimesC,math_identity2nd(size_reduced),1.0e-12_pReal))
     if (debugGeneral .or. errmatinv) then
       write(formatString, '(i2)') size_reduced
       formatString = '(/,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))'
@@ -757,15 +734,18 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
       write(6,trim(formatString),advance='no') ' S (load) ', transpose(s_reduced)
       if(errmatinv) call IO_error(error_ID=400,ext_msg='utilities_maskedCompliance')
     endif
+    temp99_real = reshape(unpack(reshape(s_reduced,[size_reduced**2]),reshape(mask,[81]),0.0_pReal),[9,9])
   else
     temp99_real = 0.0_pReal
   endif
+
+  utilities_maskedCompliance = math_99to3333(temp99_Real)
+
   if(debugGeneral) then
     write(6,'(/,a,/,9(9(2x,f10.5,1x)/),/)',advance='no') &
       ' Masked Compliance (load) * GPa =', transpose(temp99_Real)*1.0e9_pReal
     flush(6)
   endif
-  utilities_maskedCompliance = math_99to3333(temp99_Real)
 
 end function utilities_maskedCompliance
 
@@ -862,7 +842,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
   write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',&
                                                       transpose(P_av)*1.e-6_pReal
   if(present(rotation_BC)) &
-    P_av = rotation_BC%rotTensor2(P_av)
+    P_av = rotation_BC%rotate(P_av)
   write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress       / MPa =',&
                                                       transpose(P_av)*1.e-6_pReal
   flush(6)

src/homogenization.f90

@@ -146,9 +146,7 @@ subroutine homogenization_init
 !--------------------------------------------------------------------------------------------------
 ! allocate and initialize global variables
   allocate(materialpoint_dPdF(3,3,3,3,discretization_nIP,discretization_nElem),       source=0.0_pReal)
-  allocate(materialpoint_F0(3,3,discretization_nIP,discretization_nElem),             source=0.0_pReal)
   materialpoint_F0 = spread(spread(math_I3,3,discretization_nIP),4,discretization_nElem)            ! initialize to identity
-  allocate(materialpoint_F(3,3,discretization_nIP,discretization_nElem),              source=0.0_pReal)
   materialpoint_F = materialpoint_F0                                                                ! initialize to identity
   allocate(materialpoint_subF0(3,3,discretization_nIP,discretization_nElem),          source=0.0_pReal)
   allocate(materialpoint_subF(3,3,discretization_nIP,discretization_nElem),           source=0.0_pReal)
@@ -333,12 +331,10 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
             !$OMP FLUSH(terminallyIll)
             if (.not. terminallyIll) then                                                           ! so first signals terminally ill...
               !$OMP CRITICAL (write2out)
-                write(6,*) 'Integration point ', i,' at element ', e, ' terminally ill'
+              write(6,*) 'Integration point ', i,' at element ', e, ' terminally ill'
               !$OMP END CRITICAL (write2out)
             endif
-            !$OMP CRITICAL (setTerminallyIll)
-              terminallyIll = .true.                                                                ! ...and kills all others
-            !$OMP END CRITICAL (setTerminallyIll)
+            terminallyIll = .true.                                                                  ! ...and kills all others
           else                                                                                      ! cutback makes sense
             materialpoint_subStep(i,e) = subStepSizeHomog * materialpoint_subStep(i,e)              ! crystallite had severe trouble, so do a significant cutback
 

src/homogenization_mech_RGC.f90

@@ -9,17 +9,6 @@
 submodule(homogenization) homogenization_mech_RGC
   use rotations
 
-  enum, bind(c) 
-    enumerator :: &
-      undefined_ID, &
-      constitutivework_ID, &
-      penaltyenergy_ID, &
-      volumediscrepancy_ID, &
-      averagerelaxrate_ID,&
-      maximumrelaxrate_ID,&
-      magnitudemismatch_ID
-  end enum
-  
   type :: tParameters
     integer, dimension(:), allocatable :: &
       Nconstituents
@@ -31,8 +20,8 @@ submodule(homogenization) homogenization_mech_RGC
       angles
     integer :: &
       of_debug = 0
-    integer(kind(undefined_ID)),  dimension(:),   allocatable   :: &
-      outputID
+    character(len=pStringLen), allocatable, dimension(:) :: &
+      output
   end type tParameters
  
   type :: tRGCstate
@@ -71,23 +60,18 @@ module subroutine mech_RGC_init
 
   integer :: &
     Ninstance, &
-    h, i, &
+    h, &
     NofMyHomog, &
     sizeState, nIntFaceTot
  
-  integer(kind(undefined_ID)) :: &
-    outputID
   
-  character(len=pStringLen), dimension(:), allocatable :: &
-    outputs
+  write(6,'(/,a)') ' <<<+-  homogenization_'//HOMOGENIZATION_RGC_label//' init  -+>>>'; flush(6)
  
-  write(6,'(/,a)')   ' <<<+-  homogenization_'//HOMOGENIZATION_RGC_label//' init  -+>>>'
+  write(6,'(/,a)') ' Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009'
+  write(6,'(a)')   ' https://doi.org/10.1007/s12289-009-0619-1'
  
-  write(6,'(/,a)')   ' Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009'
-  write(6,'(a)')     ' https://doi.org/10.1007/s12289-009-0619-1'
- 
-  write(6,'(/,a)')   ' Tjahjanto et al., Modelling and Simulation in Materials Science and Engineering 18:015006, 2010'
-  write(6,'(a)')     ' https://doi.org/10.1088/0965-0393/18/1/015006'
+  write(6,'(/,a)') ' Tjahjanto et al., Modelling and Simulation in Materials Science and Engineering 18:015006, 2010'
+  write(6,'(a)')   ' https://doi.org/10.1088/0965-0393/18/1/015006'
  
   Ninstance = count(homogenization_type == HOMOGENIZATION_RGC_ID)
   if (iand(debug_level(debug_HOMOGENIZATION),debug_levelBasic) /= 0) &
@@ -123,33 +107,7 @@ module subroutine mech_RGC_init
     prm%dAlpha  = config%getFloats('grainsize',         requiredSize=3)
     prm%angles  = config%getFloats('clusterorientation',requiredSize=3)
  
-    outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
-    allocate(prm%outputID(0))
- 
-    do i=1, size(outputs)
-      outputID = undefined_ID
-      select case(outputs(i))
-      
-        case('constitutivework')
-          outputID = constitutivework_ID
-        case('penaltyenergy')
-          outputID = penaltyenergy_ID
-        case('volumediscrepancy')
-          outputID = volumediscrepancy_ID
-        case('averagerelaxrate')
-          outputID = averagerelaxrate_ID
-        case('maximumrelaxrate')
-          outputID = maximumrelaxrate_ID
-        case('magnitudemismatch')
-          outputID = magnitudemismatch_ID
- 
-      end select
-      
-      if (outputID /= undefined_ID) then
-        prm%outputID = [prm%outputID , outputID]
-      endif
-      
-    enddo
+    prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
 
     NofMyHomog = count(material_homogenizationAt == h)
     nIntFaceTot = 3*(  (prm%Nconstituents(1)-1)*prm%Nconstituents(2)*prm%Nconstituents(3) &
@@ -711,7 +669,7 @@ module procedure mech_RGC_updateState
        nDef = 0.0_pReal
        do i = 1,3; do j = 1,3
          do k = 1,3; do l = 1,3
-           nDef(i,j) = nDef(i,j) - nVect(k)*gDef(i,l)*math_civita(j,k,l)                            ! compute the interface mismatch tensor from the jump of deformation gradient
+           nDef(i,j) = nDef(i,j) - nVect(k)*gDef(i,l)*math_LeviCivita(j,k,l)                        ! compute the interface mismatch tensor from the jump of deformation gradient
          enddo; enddo
          nDefNorm = nDefNorm + nDef(i,j)**2.0_pReal                                                 ! compute the norm of the mismatch tensor
        enddo; enddo
@@ -731,7 +689,7 @@ module procedure mech_RGC_updateState
          rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*bgGrain + muGNghb*bgGNghb)*prm%xiAlpha &
                                                 *surfCorr(abs(intFace(1)))/prm%dAlpha(abs(intFace(1))) &
                                                 *cosh(prm%ciAlpha*nDefNorm) &
-                                                *0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_civita(k,l,j) &
+                                                *0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_LeviCivita(k,l,j) &
                                                 *tanh(nDefNorm/xSmoo_RGC)
        enddo; enddo;enddo; enddo
      enddo interfaceLoop
@@ -934,26 +892,24 @@ module subroutine mech_RGC_results(instance,group)
   integer :: o
   
   associate(stt => state(instance), dst => dependentState(instance), prm => param(instance))
-
-  outputsLoop: do o = 1,size(prm%outputID)
-    select case(prm%outputID(o))
-    
-      case (constitutivework_ID)
+  outputsLoop: do o = 1,size(prm%output)
+    select case(trim(prm%output(o)))
+      case('constitutivework')
         call results_writeDataset(group,stt%work,'W',&
                                   'work density','J/m³')
-      case (magnitudemismatch_ID)
+      case('magnitudemismatch')
         call results_writeDataset(group,dst%mismatch,'N',&
                                   'average mismatch tensor','1')
-      case (penaltyenergy_ID)
+      case('penaltyenergy')
         call results_writeDataset(group,stt%penaltyEnergy,'R',&
                                   'mismatch penalty density','J/m³')
-      case (volumediscrepancy_ID)
+      case('volumediscrepancy')
         call results_writeDataset(group,dst%volumeDiscrepancy,'Delta_V',&
                                   'volume discrepancy','m³')
-      case (maximumrelaxrate_ID)
+      case('maximumrelaxrate')
         call results_writeDataset(group,dst%relaxationrate_max,'max_alpha_dot',&
                                   'maximum relaxation rate','m/s')
-      case (averagerelaxrate_ID)
+      case('averagerelaxrate')
         call results_writeDataset(group,dst%relaxationrate_avg,'avg_alpha_dot',&
                                   'average relaxation rate','m/s')
     end select

src/lattice.f90

@@ -524,7 +524,7 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine lattice_init
 
-  integer :: Nphases, p
+  integer :: Nphases, p,i
   character(len=pStringLen) :: &
     tag  = ''
   real(pReal) :: CoverA
@@ -550,27 +550,11 @@ subroutine lattice_init
   allocate(lattice_mu(Nphases), source=0.0_pReal)
   allocate(lattice_nu(Nphases), source=0.0_pReal)
 
- 
   allocate(lattice_Scleavage(3,3,3,lattice_maxNcleavage,Nphases),source=0.0_pReal)
   allocate(lattice_NcleavageSystem(lattice_maxNcleavageFamily,Nphases),source=0)
 
-  do p = 1, size(config_phase)
-    tag = config_phase(p)%getString('lattice_structure')
-    select case(trim(tag(1:3)))
-      case('iso')
-        lattice_structure(p) = LATTICE_iso_ID
-      case('fcc')
-        lattice_structure(p) = LATTICE_fcc_ID
-      case('bcc')
-        lattice_structure(p) = LATTICE_bcc_ID
-      case('hex')
-        lattice_structure(p) = LATTICE_hex_ID
-      case('bct')
-        lattice_structure(p) = LATTICE_bct_ID
-      case('ort')
-        lattice_structure(p) = LATTICE_ort_ID
-    end select
 
+  do p = 1, size(config_phase)
 
     lattice_C66(1,1,p) = config_phase(p)%getFloat('c11',defaultVal=0.0_pReal)
     lattice_C66(1,2,p) = config_phase(p)%getFloat('c12',defaultVal=0.0_pReal)
@@ -582,9 +566,44 @@ subroutine lattice_init
     lattice_C66(5,5,p) = config_phase(p)%getFloat('c55',defaultVal=0.0_pReal)
     lattice_C66(6,6,p) = config_phase(p)%getFloat('c66',defaultVal=0.0_pReal)
 
- 
     CoverA = config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal)
 
+    tag = config_phase(p)%getString('lattice_structure')
+    select case(tag(1:3))
+      case('iso')
+        lattice_structure(p) = LATTICE_iso_ID
+      case('fcc')
+        lattice_structure(p) = LATTICE_fcc_ID
+      case('bcc')
+        lattice_structure(p) = LATTICE_bcc_ID
+      case('hex')
+        if(CoverA < 1.0_pReal .or. CoverA > 2.0_pReal) call IO_error(131,el=p)
+        lattice_structure(p) = LATTICE_hex_ID
+      case('bct')
+        if(CoverA > 2.0_pReal) call IO_error(131,el=p)
+        lattice_structure(p) = LATTICE_bct_ID
+      case('ort')
+        lattice_structure(p) = LATTICE_ort_ID
+      case default
+        call IO_error(130,ext_msg='lattice_init')
+    end select
+
+    lattice_C66(1:6,1:6,p) = lattice_symmetrizeC66(lattice_structure(p),lattice_C66(1:6,1:6,p))
+
+    lattice_mu(p) = 0.2_pReal *(lattice_C66(1,1,p) -lattice_C66(1,2,p) +3.0_pReal*lattice_C66(4,4,p)) ! (C11iso-C12iso)/2 with C11iso=(3*C11+2*C12+4*C44)/5 and C12iso=(C11+4*C12-2*C44)/5
+    lattice_nu(p) = (          lattice_C66(1,1,p) +4.0_pReal*lattice_C66(1,2,p) -2.0_pReal*lattice_C66(4,4,p)) &
+                  / (4.0_pReal*lattice_C66(1,1,p) +6.0_pReal*lattice_C66(1,2,p) +2.0_pReal*lattice_C66(4,4,p))! C12iso/(C11iso+C12iso) with C11iso=(3*C11+2*C12+4*C44)/5 and C12iso=(C11+4*C12-2*C44)/5
+
+    lattice_C3333(1:3,1:3,1:3,1:3,p) = math_Voigt66to3333(lattice_C66(1:6,1:6,p))                   ! Literature data is Voigt
+    lattice_C66(1:6,1:6,p)           = math_sym3333to66(lattice_C3333(1:3,1:3,1:3,1:3,p))           ! DAMASK uses Mandel-weighting
+
+    do i = 1, 6
+      if (abs(lattice_C66(i,i,p))<tol_math_check) &
+        call IO_error(135,el=i,ip=p,ext_msg='matrix diagonal "el"ement of phase "ip"')
+    enddo
+
+
+    ! should not be part of lattice
     lattice_thermalConductivity33(1,1,p) = config_phase(p)%getFloat('thermal_conductivity11',defaultVal=0.0_pReal)
     lattice_thermalConductivity33(2,2,p) = config_phase(p)%getFloat('thermal_conductivity22',defaultVal=0.0_pReal)
     lattice_thermalConductivity33(3,3,p) = config_phase(p)%getFloat('thermal_conductivity33',defaultVal=0.0_pReal)
@@ -604,10 +623,6 @@ subroutine lattice_init
     lattice_DamageDiffusion33(3,3,p) = config_phase(p)%getFloat( 'damage_diffusion33',defaultVal=0.0_pReal)
     lattice_DamageMobility(p) = config_phase(p)%getFloat( 'damage_mobility',defaultVal=0.0_pReal)
 
-    if ((CoverA < 1.0_pReal .or. CoverA > 2.0_pReal) &
-        .and. lattice_structure(p) == LATTICE_hex_ID) call IO_error(131,el=p)                       ! checking physical significance of c/a
-    if ((CoverA > 2.0_pReal) &
-        .and. lattice_structure(p) == LATTICE_bct_ID) call IO_error(131,el=p)                       ! checking physical significance of c/a
     call lattice_initializeStructure(p, CoverA)
   enddo
 
@@ -624,27 +639,7 @@ subroutine lattice_initializeStructure(myPhase,CoverA)
     CoverA
 
   integer :: &
-   i,  &
-    myNcleavage
- 
-  lattice_C66(1:6,1:6,myPhase) = lattice_symmetrizeC66(lattice_structure(myPhase),&
-                                                       lattice_C66(1:6,1:6,myPhase))
- 
-  lattice_mu(myPhase) = 0.2_pReal *(  lattice_C66(1,1,myPhase) &
-                                    - lattice_C66(1,2,myPhase) &
-                                    + 3.0_pReal*lattice_C66(4,4,myPhase))                           ! (C11iso-C12iso)/2 with C11iso=(3*C11+2*C12+4*C44)/5 and C12iso=(C11+4*C12-2*C44)/5
-  lattice_nu(myPhase) = (  lattice_C66(1,1,myPhase) &
-                         + 4.0_pReal*lattice_C66(1,2,myPhase) &
-                         - 2.0_pReal*lattice_C66(4,4,myPhase)) &
-                                                              /(  4.0_pReal*lattice_C66(1,1,myPhase) &
-                                                                + 6.0_pReal*lattice_C66(1,2,myPhase) &
-                                                                + 2.0_pReal*lattice_C66(4,4,myPhase))! C12iso/(C11iso+C12iso) with C11iso=(3*C11+2*C12+4*C44)/5 and C12iso=(C11+4*C12-2*C44)/5
-  lattice_C3333(1:3,1:3,1:3,1:3,myPhase) = math_Voigt66to3333(lattice_C66(1:6,1:6,myPhase))         ! Literature data is Voigt
-  lattice_C66(1:6,1:6,myPhase) = math_sym3333to66(lattice_C3333(1:3,1:3,1:3,1:3,myPhase))           ! DAMASK uses Mandel-weighting
-  do i = 1, 6
-    if (abs(lattice_C66(i,i,myPhase))<tol_math_check) &
-      call IO_error(135,el=i,ip=myPhase,ext_msg='matrix diagonal "el"ement of phase "ip"')
-  enddo
+   i
 
   do i = 1,3
     lattice_thermalExpansion33 (1:3,1:3,i,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),&
@@ -655,51 +650,34 @@ subroutine lattice_initializeStructure(myPhase,CoverA)
                                                                           lattice_thermalConductivity33  (1:3,1:3,myPhase))
   lattice_DamageDiffusion33      (1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),&
                                                                           lattice_DamageDiffusion33      (1:3,1:3,myPhase))
-  myNcleavage   = 0
 
   select case(lattice_structure(myPhase))
 
     case (LATTICE_fcc_ID)
-      myNcleavage = lattice_fcc_Ncleavage
       lattice_NcleavageSystem(1:2,myPhase)  = lattice_fcc_NcleavageSystem
- 
-      lattice_Scleavage(1:3,1:3,1:3,1:myNcleavage,myPhase) = &
+      lattice_Scleavage(1:3,1:3,1:3,1:lattice_fcc_Ncleavage,myPhase) = &
       lattice_SchmidMatrix_cleavage(lattice_fcc_ncleavageSystem,'fcc',covera)
 
     case (LATTICE_bcc_ID)
-      myNcleavage = lattice_bcc_Ncleavage
       lattice_NcleavageSystem(1:2,myPhase)  = lattice_bcc_NcleavageSystem
- 
-      lattice_Scleavage(1:3,1:3,1:3,1:myNcleavage,myPhase) = &
+      lattice_Scleavage(1:3,1:3,1:3,1:lattice_bcc_Ncleavage,myPhase) = &
       lattice_SchmidMatrix_cleavage(lattice_bcc_ncleavagesystem,'bcc',covera)
 
     case (LATTICE_hex_ID)
-      myNcleavage = lattice_hex_Ncleavage
       lattice_NcleavageSystem(1:1,myPhase) = lattice_hex_NcleavageSystem
- 
-      lattice_Scleavage(1:3,1:3,1:3,1:myNcleavage,myPhase) = &
+      lattice_Scleavage(1:3,1:3,1:3,1:lattice_hex_Ncleavage,myPhase) = &
       lattice_SchmidMatrix_cleavage(lattice_hex_ncleavagesystem,'hex',covera)
 
-    case (LATTICE_bct_ID)
- 
     case (LATTICE_ort_ID)
-      myNcleavage   = lattice_ort_Ncleavage
       lattice_NcleavageSystem(1:3,myPhase)  = lattice_ort_NcleavageSystem
- 
-      lattice_Scleavage(1:3,1:3,1:3,1:myNcleavage,myPhase) = &
+      lattice_Scleavage(1:3,1:3,1:3,1:lattice_ort_Ncleavage,myPhase) = &
       lattice_SchmidMatrix_cleavage(lattice_ort_NcleavageSystem,'ort',covera)
 
     case (LATTICE_iso_ID)
-      myNcleavage   = lattice_iso_Ncleavage
       lattice_NcleavageSystem(1:1,myPhase)  = lattice_iso_NcleavageSystem
- 
-      lattice_Scleavage(1:3,1:3,1:3,1:myNcleavage,myPhase) = &
+      lattice_Scleavage(1:3,1:3,1:3,1:lattice_iso_Ncleavage,myPhase) = &
       lattice_SchmidMatrix_cleavage(lattice_iso_NcleavageSystem,'iso',covera)
 
-!--------------------------------------------------------------------------------------------------
-! something went wrong
-    case default
-      call IO_error(130,ext_msg='lattice_initializeStructure')
   end select
 
 end subroutine lattice_initializeStructure

src/material.f90

@@ -268,6 +268,7 @@ subroutine material_init
     if(microstructure_Nconstituents(m) < 1) &
          call IO_error(151,m)
   enddo
+  if(homogenization_maxNgrains > size(microstructure_phase,1)) call IO_error(148)
  
   debugOut: if (iand(myDebug,debug_levelExtensive) /= 0) then
     write(6,'(/,a,/)') ' MATERIAL configuration'
@@ -290,6 +291,7 @@ subroutine material_init
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 ! new mappings
+
   allocate(material_phaseAt(homogenization_maxNgrains,discretization_nElem),                   source=0)
   allocate(material_texture(homogenization_maxNgrains,discretization_nIP,discretization_nElem),source=0)   !this is only needed by plasticity nonlocal
   allocate(material_orientation0(homogenization_maxNgrains,discretization_nIP,discretization_nElem))
@@ -298,15 +300,24 @@ subroutine material_init
      do i = 1, discretization_nIP
        myMicro = discretization_microstructureAt(e)
        do c = 1, homogenization_Ngrains(discretization_homogenizationAt(e))
-         material_phaseAt(c,e)        = microstructure_phase(c,myMicro)
-         material_texture(c,i,e)      = microstructure_texture(c,myMicro)
-         material_orientation0(c,i,e) = texture_orientation(material_texture(c,i,e))
+         if(microstructure_phase(c,myMicro) > 0) then
+           material_phaseAt(c,e)        = microstructure_phase(c,myMicro)
+         else
+           call IO_error(150,ext_msg='phase')
+         endif
+         if(microstructure_texture(c,myMicro) > 0) then
+           material_texture(c,i,e)      = microstructure_texture(c,myMicro)
+           material_orientation0(c,i,e) = texture_orientation(material_texture(c,i,e))
+         else
+           call IO_error(150,ext_msg='texture')
+         endif
        enddo
      enddo
    enddo
  
   deallocate(microstructure_phase)
   deallocate(microstructure_texture)
+  deallocate(texture_orientation)
  
   
   allocate(material_homogenizationAt,source=discretization_homogenizationAt)
@@ -464,7 +475,7 @@ subroutine material_parseMicrostructure
   real(pReal), dimension(:,:), allocatable :: &
     microstructure_fraction                                                                         !< vol fraction of each constituent in microstructure
   integer :: &
-    microstructure_maxNconstituents                                                                 !< max number of constituents in any phase
+    maxNconstituents                                                                                !< max number of constituents in any phase
  
   allocate(microstructure_Nconstituents(size(config_microstructure)), source=0)
  
@@ -475,10 +486,10 @@ subroutine material_parseMicrostructure
     microstructure_Nconstituents(m) =  config_microstructure(m)%countKeys('(constituent)')
   enddo
  
-  microstructure_maxNconstituents = maxval(microstructure_Nconstituents)
-  allocate(microstructure_phase   (microstructure_maxNconstituents,size(config_microstructure)),source=0)
-  allocate(microstructure_texture (microstructure_maxNconstituents,size(config_microstructure)),source=0)
-  allocate(microstructure_fraction(microstructure_maxNconstituents,size(config_microstructure)),source=0.0_pReal)
+  maxNconstituents = maxval(microstructure_Nconstituents)
+  allocate(microstructure_phase   (maxNconstituents,size(config_microstructure)),source=0)
+  allocate(microstructure_texture (maxNconstituents,size(config_microstructure)),source=0)
+  allocate(microstructure_fraction(maxNconstituents,size(config_microstructure)),source=0.0_pReal)
  
   allocate(strings(1))                                                                              ! Intel 16.0 Bug
   do m=1, size(config_microstructure)

src/math.f90

@@ -73,6 +73,11 @@ module math
       3,3  &
       ],[2,9])                                                                                      !< arrangement in Plain notation
 
+
+  interface math_mul33xx33
+    module procedure math_tensordot
+  end interface math_mul33xx33
+
 !---------------------------------------------------------------------------------------------------
  private :: &
    unitTest
@@ -266,31 +271,30 @@ end function math_identity4th
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief permutation tensor e_ijk used for computing cross product of two tensors
+!> @brief permutation tensor e_ijk
 ! e_ijk =  1 if even permutation of ijk
 ! e_ijk = -1 if odd permutation of ijk
 ! e_ijk =  0 otherwise
 !--------------------------------------------------------------------------------------------------
-real(pReal) pure function math_civita(i,j,k)
+real(pReal) pure function math_LeviCivita(i,j,k)
 
   integer, intent(in) :: i,j,k
 
-  math_civita = 0.0_pReal
-  if (((i == 1).and.(j == 2).and.(k == 3)) .or. &
-      ((i == 2).and.(j == 3).and.(k == 1)) .or. &
-      ((i == 3).and.(j == 1).and.(k == 2))) math_civita = 1.0_pReal
-  if (((i == 1).and.(j == 3).and.(k == 2)) .or. &
-      ((i == 2).and.(j == 1).and.(k == 3)) .or. &
-      ((i == 3).and.(j == 2).and.(k == 1))) math_civita = -1.0_pReal
+  if     (all([i,j,k] == [1,2,3]) .or. all([i,j,k] == [2,3,1]) .or. all([i,j,k] == [3,1,2])) then
+    math_LeviCivita = +1.0_pReal
+  elseif (all([i,j,k] == [3,2,1]) .or. all([i,j,k] == [2,1,3]) .or. all([i,j,k] == [1,3,2])) then
+    math_LeviCivita = -1.0_pReal
+  else
+    math_LeviCivita =  0.0_pReal
+  endif
 
-end function math_civita
+end function math_LeviCivita
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief kronecker delta function d_ij
 ! d_ij = 1 if i = j
 ! d_ij = 0 otherwise
-! inspired by http://fortraninacworld.blogspot.de/2012/12/ternary-operator.html
 !--------------------------------------------------------------------------------------------------
 real(pReal) pure function math_delta(i,j)
 
@@ -317,7 +321,7 @@ end function math_cross
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief outer product A \otimes B of arbitrary sized vectors A and B
+!> @brief outer product of arbitrary sized vectors (A ⊗ B / i,j)
 !--------------------------------------------------------------------------------------------------
 pure function math_outer(A,B)
 
@@ -333,7 +337,7 @@ end function math_outer
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief outer product A \otimes B of arbitrary sized vectors A and B
+!> @brief inner product of arbitrary sized vectors (A · B / i,i)
 !--------------------------------------------------------------------------------------------------
 real(pReal) pure function math_inner(A,B)
 
@@ -346,24 +350,19 @@ end function math_inner
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief matrix multiplication 33xx33 = 1 (double contraction --> ij * ij)
+!> @brief double contraction of 3x3 matrices (A : B / ij,ij)
 !--------------------------------------------------------------------------------------------------
-real(pReal) pure function math_mul33xx33(A,B)
+real(pReal) pure function math_tensordot(A,B)
 
   real(pReal), dimension(3,3), intent(in) :: A,B
-  integer :: i,j
-  real(pReal), dimension(3,3) ::             C
 
-  do i=1,3; do j=1,3
-    C(i,j) = A(i,j) * B(i,j)
-  enddo; enddo
-  math_mul33xx33 = sum(C)
+  math_tensordot = sum(A*B)
 
-end function math_mul33xx33
+end function math_tensordot
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief matrix multiplication 3333x33 = 33 (double contraction --> ijkl *kl = ij)
+!> @brief matrix double contraction 3333x33 = 33 (ijkl,kl)
 !--------------------------------------------------------------------------------------------------
 pure function math_mul3333xx33(A,B)
 
@@ -380,7 +379,7 @@ end function math_mul3333xx33
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief matrix multiplication 3333x3333 = 3333 (ijkl *klmn = ijmn)
+!> @brief matrix multiplication 3333x3333 = 3333 (ijkl,klmn)
 !--------------------------------------------------------------------------------------------------
 pure function math_mul3333xx3333(A,B)
 
@@ -402,8 +401,9 @@ end function math_mul3333xx3333
 pure function math_exp33(A,n)
 
   real(pReal), dimension(3,3), intent(in) :: A
-  integer, intent(in), optional :: n
+  integer,                     intent(in), optional :: n
   real(pReal), dimension(3,3) :: B, math_exp33
+
   real(pReal) :: invFac
   integer     :: n_,i
 
@@ -433,10 +433,11 @@ end function math_exp33
 !--------------------------------------------------------------------------------------------------
 pure function math_inv33(A)
 
-  real(pReal),dimension(3,3),intent(in)  :: A
-  real(pReal)                :: DetA
-  real(pReal),dimension(3,3) :: math_inv33
-  logical                    :: error
+  real(pReal), dimension(3,3), intent(in) :: A
+  real(pReal), dimension(3,3) :: math_inv33
+
+  real(pReal) :: DetA
+  logical     :: error
 
   call math_invert33(math_inv33,DetA,error,A)
   if(error) math_inv33 = 0.0_pReal
@@ -451,10 +452,10 @@ end function math_inv33
 !--------------------------------------------------------------------------------------------------
 pure subroutine math_invert33(InvA, DetA, error, A)
 
-  logical, intent(out) :: error
-  real(pReal),dimension(3,3),intent(in)  :: A
-  real(pReal),dimension(3,3),intent(out) :: InvA
-  real(pReal), intent(out) :: DetA
+  real(pReal), dimension(3,3), intent(out) :: InvA
+  real(pReal),                 intent(out) :: DetA
+  logical,                     intent(out) :: error
+  real(pReal), dimension(3,3), intent(in)  :: A
 
   InvA(1,1) =  A(2,2) * A(3,3) - A(2,3) * A(3,2)
   InvA(2,1) = -A(2,1) * A(3,3) + A(2,3) * A(3,1)
@@ -1307,6 +1308,7 @@ subroutine unitTest
     sort_out_  = reshape([-1,-1,  +1,+5,  +5,+6,  +3,-2],[2,4])
 
   integer,     dimension(5) :: range_out_ = [1,2,3,4,5]
+  integer,     dimension(3) :: ijk
 
   real(pReal)                 :: det
   real(pReal), dimension(3)   :: v3_1,v3_2,v3_3,v3_4
@@ -1415,13 +1417,23 @@ subroutine unitTest
     call IO_error(0,ext_msg='math_invert(t99)')
 
   if(any(dNeq(math_clip([4.0_pReal,9.0_pReal],5.0_pReal,6.5_pReal),[5.0_pReal,6.5_pReal]))) &
-     call IO_error(0,ext_msg='math_clip')
+    call IO_error(0,ext_msg='math_clip')
 
   if(math_factorial(10) /= 3628800) &
-     call IO_error(0,ext_msg='math_factorial')
+    call IO_error(0,ext_msg='math_factorial')
 
   if(math_binomial(49,6) /= 13983816) &
-     call IO_error(0,ext_msg='math_binomial')
+    call IO_error(0,ext_msg='math_binomial')
+
+  ijk = cshift([1,2,3],int(r*1.0e2_pReal))
+  if(dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),+1.0_pReal)) &
+    call IO_error(0,ext_msg='math_LeviCivita(even)')
+  ijk = cshift([3,2,1],int(r*2.0e2_pReal))
+  if(dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),-1.0_pReal)) &
+    call IO_error(0,ext_msg='math_LeviCivita(odd)')
+  ijk = cshift([2,2,1],int(r*2.0e2_pReal))
+  if(dNeq0(math_LeviCivita(ijk(1),ijk(2),ijk(3))))&
+    call IO_error(0,ext_msg='math_LeviCivita')
 
 end subroutine unitTest
 

src/results.f90

@@ -485,15 +485,15 @@ end subroutine results_writeScalarDataset_rotation
 !--------------------------------------------------------------------------------------------------
 !> @brief adds the unique mapping from spatial position and constituent ID to results
 !--------------------------------------------------------------------------------------------------
-subroutine results_mapping_constituent(phaseAt,memberAt,label)
+subroutine results_mapping_constituent(phaseAt,memberAtLocal,label)
     
   integer,          dimension(:,:),   intent(in) :: phaseAt                                         !< phase section at (constituent,element)
-  integer,          dimension(:,:,:), intent(in) :: memberAt                                        !< phase member at (constituent,IP,element)
-  character(len=*), dimension(:),     intent(in) :: label                                           !< label of each phase section
+  integer,                   dimension(:,:,:), intent(in) :: memberAtLocal                          !< phase member at (constituent,IP,element)
+  character(len=pStringLen), dimension(:),     intent(in) :: label                                  !< label of each phase section
   
-  integer, dimension(size(memberAt,1),size(memberAt,2),size(memberAt,3)) :: &
-    phaseAt_perIP, &
-    memberAt_total
+  integer, dimension(size(memberAtLocal,1),size(memberAtLocal,2),size(memberAtLocal,3)) :: &
+    phaseAtMaterialpoint, &
+    memberAtGlobal
   integer, dimension(size(label),0:worldsize-1) :: memberOffset                                     !< offset in member counting per process
   integer, dimension(0:worldsize-1)             :: writeSize                                        !< amount of data written per process
   integer(HSIZE_T), dimension(2) :: &
@@ -543,10 +543,10 @@ subroutine results_mapping_constituent(phaseAt,memberAt,label)
   call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, ierr)
   memberOffset = 0
   do i=1, size(label)
-    memberOffset(i,worldrank) = count(phaseAt == i)*size(memberAt,2)                                ! number of points/instance of this process
+    memberOffset(i,worldrank) = count(phaseAt == i)*size(memberAtLocal,2)                                ! number of points/instance of this process
   enddo
   writeSize = 0
-  writeSize(worldrank) = size(memberAt(1,:,:))                                                      ! total number of points by this process
+  writeSize(worldrank) = size(memberAtLocal(1,:,:))                                                      ! total number of points by this process
 
 !--------------------------------------------------------------------------------------------------
 ! MPI settings and communication
@@ -578,14 +578,14 @@ subroutine results_mapping_constituent(phaseAt,memberAt,label)
 
 !---------------------------------------------------------------------------------------------------
 ! expand phaseAt to consider IPs (is not stored per IP)
-  do i = 1, size(phaseAt_perIP,2)
-    phaseAt_perIP(:,i,:) = phaseAt
+  do i = 1, size(phaseAtMaterialpoint,2)
+    phaseAtMaterialpoint(:,i,:) = phaseAt
   enddo
   
 !---------------------------------------------------------------------------------------------------
 ! renumber member from my process to all processes
   do i = 1, size(label)
-    where(phaseAt_perIP == i) memberAt_total = memberAt + sum(memberOffset(i,0:worldrank-1)) -1     ! convert to 0-based
+    where(phaseAtMaterialpoint == i) memberAtGlobal = memberAtLocal + sum(memberOffset(i,0:worldrank-1)) -1     ! convert to 0-based
   enddo
 
 !--------------------------------------------------------------------------------------------------
@@ -596,10 +596,10 @@ subroutine results_mapping_constituent(phaseAt,memberAt,label)
   call h5dcreate_f(loc_id, 'constituent', dtype_id, filespace_id, dset_id, ierr)
   if (ierr < 0) call IO_error(1,ext_msg='results_mapping_constituent: h5dcreate_f')
   
-  call h5dwrite_f(dset_id, name_id, reshape(label(pack(phaseAt_perIP,.true.)),myShape), &
+  call h5dwrite_f(dset_id, name_id, reshape(label(pack(phaseAtMaterialpoint,.true.)),myShape), &
                   myShape, ierr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
   if (ierr < 0) call IO_error(1,ext_msg='results_mapping_constituent: h5dwrite_f/name_id')
-  call h5dwrite_f(dset_id, position_id, reshape(pack(memberAt_total,.true.),myShape), &
+  call h5dwrite_f(dset_id, position_id, reshape(pack(memberAtGlobal,.true.),myShape), &
                   myShape, ierr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
   if (ierr < 0) call IO_error(1,ext_msg='results_mapping_constituent: h5dwrite_f/position_id')
 
@@ -620,15 +620,15 @@ end subroutine results_mapping_constituent
 !--------------------------------------------------------------------------------------------------
 !> @brief adds the unique mapping from spatial position and constituent ID to results
 !--------------------------------------------------------------------------------------------------
-subroutine results_mapping_materialpoint(homogenizationAt,memberAt,label)
+subroutine results_mapping_materialpoint(homogenizationAt,memberAtLocal,label)
     
   integer,          dimension(:),   intent(in) :: homogenizationAt                                  !< homogenization section at (element)
-  integer,          dimension(:,:), intent(in) :: memberAt                                          !< homogenization member at (IP,element)
-  character(len=*), dimension(:),   intent(in) :: label                                             !< label of each homogenization section
+  integer,                   dimension(:,:), intent(in) :: memberAtLocal                            !< homogenization member at (IP,element)
+  character(len=pStringLen), dimension(:),   intent(in) :: label                                    !< label of each homogenization section
   
-  integer, dimension(size(memberAt,1),size(memberAt,2)) :: &
-    homogenizationAt_perIP, &
-    memberAt_total
+  integer, dimension(size(memberAtLocal,1),size(memberAtLocal,2)) :: &
+    homogenizationAtMaterialpoint, &
+    memberAtGlobal
   integer, dimension(size(label),0:worldsize-1) :: memberOffset                                     !< offset in member counting per process
   integer, dimension(0:worldsize-1)             :: writeSize                                        !< amount of data written per process
   integer(HSIZE_T), dimension(1) :: &
@@ -678,10 +678,10 @@ subroutine results_mapping_materialpoint(homogenizationAt,memberAt,label)
   call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, ierr)
   memberOffset = 0
   do i=1, size(label)
-    memberOffset(i,worldrank) = count(homogenizationAt == i)*size(memberAt,1)                       ! number of points/instance of this process
+    memberOffset(i,worldrank) = count(homogenizationAt == i)*size(memberAtLocal,1)                       ! number of points/instance of this process
   enddo
   writeSize = 0
-  writeSize(worldrank) = size(memberAt)                                                             ! total number of points by this process
+  writeSize(worldrank) = size(memberAtLocal)                                                             ! total number of points by this process
 
 !--------------------------------------------------------------------------------------------------
 ! MPI settings and communication
@@ -713,14 +713,14 @@ subroutine results_mapping_materialpoint(homogenizationAt,memberAt,label)
 
 !---------------------------------------------------------------------------------------------------
 ! expand phaseAt to consider IPs (is not stored per IP)
-  do i = 1, size(homogenizationAt_perIP,1)
-    homogenizationAt_perIP(i,:) = homogenizationAt
+  do i = 1, size(homogenizationAtMaterialpoint,1)
+    homogenizationAtMaterialpoint(i,:) = homogenizationAt
   enddo
   
 !---------------------------------------------------------------------------------------------------
 ! renumber member from my process to all processes
   do i = 1, size(label)
-    where(homogenizationAt_perIP == i) memberAt_total = memberAt + sum(memberOffset(i,0:worldrank-1)) - 1  ! convert to 0-based
+    where(homogenizationAtMaterialpoint == i) memberAtGlobal = memberAtLocal + sum(memberOffset(i,0:worldrank-1)) - 1  ! convert to 0-based
   enddo
 
 !--------------------------------------------------------------------------------------------------
@@ -731,10 +731,10 @@ subroutine results_mapping_materialpoint(homogenizationAt,memberAt,label)
   call h5dcreate_f(loc_id, 'materialpoint', dtype_id, filespace_id, dset_id, ierr)
   if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5dcreate_f')
   
-  call h5dwrite_f(dset_id, name_id, reshape(label(pack(homogenizationAt_perIP,.true.)),myShape), &
+  call h5dwrite_f(dset_id, name_id, reshape(label(pack(homogenizationAtMaterialpoint,.true.)),myShape), &
                   myShape, ierr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
   if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5dwrite_f/name_id')
-  call h5dwrite_f(dset_id, position_id, reshape(pack(memberAt_total,.true.),myShape), &
+  call h5dwrite_f(dset_id, position_id, reshape(pack(memberAtGlobal,.true.),myShape), &
                   myShape, ierr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
   if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5dwrite_f/position_id')
 

src/thermal_adiabatic.f90

@@ -16,14 +16,9 @@ module thermal_adiabatic
   implicit none
   private
 
-  enum, bind(c) 
-    enumerator :: undefined_ID, &
-                  temperature_ID
-  end enum
-  
   type :: tParameters
-    integer(kind(undefined_ID)), dimension(:), allocatable :: &
-      outputID
+    character(len=pStringLen), allocatable, dimension(:) :: &
+      output
   end type tParameters
   
   type(tparameters),             dimension(:), allocatable :: &
@@ -46,10 +41,9 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine thermal_adiabatic_init
  
-  integer :: maxNinstance,o,h,NofMyHomog   
-  character(len=pStringLen), dimension(:), allocatable :: outputs
+  integer :: maxNinstance,h,NofMyHomog   
  
-  write(6,'(/,a)')   ' <<<+-  thermal_'//THERMAL_ADIABATIC_label//' init  -+>>>'; flush(6)
+  write(6,'(/,a)') ' <<<+-  thermal_'//THERMAL_ADIABATIC_label//' init  -+>>>'; flush(6)
   
   maxNinstance = count(thermal_type == THERMAL_adiabatic_ID)
   if (maxNinstance == 0) return
@@ -60,15 +54,7 @@ subroutine thermal_adiabatic_init
     if (thermal_type(h) /= THERMAL_adiabatic_ID) cycle
     associate(prm => param(thermal_typeInstance(h)),config => config_homogenization(h))
               
-    outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
-    allocate(prm%outputID(0))
-
-    do o=1, size(outputs)
-      select case(outputs(o))
-        case('temperature')
-          prm%outputID = [prm%outputID, temperature_ID]
-      end select
-    enddo
+    prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
 
     NofMyHomog=count(material_homogenizationAt==h)
     thermalState(h)%sizeState = 1
@@ -76,7 +62,6 @@ subroutine thermal_adiabatic_init
     allocate(thermalState(h)%subState0(1,NofMyHomog), source=thermal_initialT(h))
     allocate(thermalState(h)%state    (1,NofMyHomog), source=thermal_initialT(h))
  
-    nullify(thermalMapping(h)%p)
     thermalMapping(h)%p => material_homogenizationMemberAt
     deallocate(temperature(h)%p)
     temperature(h)%p => thermalState(h)%state(1,:)
@@ -246,14 +231,13 @@ subroutine thermal_adiabatic_results(homog,group)
 
   integer,          intent(in) :: homog
   character(len=*), intent(in) :: group
+
   integer :: o
   
   associate(prm => param(damage_typeInstance(homog)))
-
-  outputsLoop: do o = 1,size(prm%outputID)
-    select case(prm%outputID(o))
-    
-      case (temperature_ID)
+  outputsLoop: do o = 1,size(prm%output)
+    select case(trim(prm%output(o)))
+      case('temperature')                                                                           ! ToDo: should be 'T'
         call results_writeDataset(group,temperature(homog)%p,'T',&
                                   'temperature','K')
     end select

src/thermal_conduction.f90

@@ -15,15 +15,9 @@ module thermal_conduction
   implicit none
   private
 
-  enum, bind(c) 
-    enumerator :: &
-      undefined_ID, &
-      temperature_ID
-  end enum
-  
   type :: tParameters
-    integer(kind(undefined_ID)), dimension(:), allocatable :: &
-      outputID
+    character(len=pStringLen), allocatable, dimension(:) :: &
+      output
   end type tParameters
   
   type(tparameters),             dimension(:), allocatable :: &
@@ -48,10 +42,9 @@ contains
 subroutine thermal_conduction_init
 
   
-  integer :: maxNinstance,o,NofMyHomog,h
-  character(len=pStringLen), dimension(:), allocatable :: outputs
+  integer :: maxNinstance,NofMyHomog,h
  
-  write(6,'(/,a)')   ' <<<+-  thermal_'//THERMAL_CONDUCTION_label//' init  -+>>>'; flush(6)
+  write(6,'(/,a)') ' <<<+-  thermal_'//THERMAL_CONDUCTION_label//' init  -+>>>'; flush(6)
   
   maxNinstance = count(thermal_type == THERMAL_conduction_ID)
   if (maxNinstance == 0) return
@@ -62,15 +55,7 @@ subroutine thermal_conduction_init
     if (thermal_type(h) /= THERMAL_conduction_ID) cycle
     associate(prm => param(thermal_typeInstance(h)),config => config_homogenization(h))
               
-    outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
-    allocate(prm%outputID(0))
-
-    do o=1, size(outputs)
-      select case(outputs(o))
-        case('temperature')
-          prm%outputID = [prm%outputID, temperature_ID]
-      end select
-    enddo
+    prm%output = config%getStrings('(output)',defaultVal=emptyStringArray)
   
     NofMyHomog=count(material_homogenizationAt==h)
     thermalState(h)%sizeState = 0
@@ -78,7 +63,6 @@ subroutine thermal_conduction_init
     allocate(thermalState(h)%subState0(0,NofMyHomog))
     allocate(thermalState(h)%state    (0,NofMyHomog))
  
-    nullify(thermalMapping(h)%p)
     thermalMapping(h)%p => material_homogenizationMemberAt
     deallocate(temperature    (h)%p)
     allocate  (temperature    (h)%p(NofMyHomog), source=thermal_initialT(h))
@@ -259,14 +243,13 @@ subroutine thermal_conduction_results(homog,group)
 
   integer,          intent(in) :: homog
   character(len=*), intent(in) :: group
+
   integer :: o
   
   associate(prm => param(damage_typeInstance(homog)))
-
-  outputsLoop: do o = 1,size(prm%outputID)
-    select case(prm%outputID(o))
-    
-      case (temperature_ID)
+  outputsLoop: do o = 1,size(prm%output)
+    select case(trim(prm%output(o)))
+      case('temperature')                                                                           ! ToDo: should be 'T'
         call results_writeDataset(group,temperature(homog)%p,'T',&
                                   'temperature','K')
     end select