Merge remote-tracking branch 'origin/development' into MiscImprovements

PRIVATE

@@ -1 +1 @@
-Subproject commit ec615d249d39e5d01446b01ab9a5b7e7601340ad
+Subproject commit 6db5f4666fc651b4de3b44ceaed3f2b848170ac9

VERSION

@@ -1 +1 @@
-v2.0.3-1732-gab88ffd2
+v2.0.3-1747-ga2e8e5b1

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),
+              damask.mechanics.Cauchy(table.get(options.stress ).reshape(-1,3,3),
-                                      table.get(options.stress ).reshape(-1,3,3)).reshape(-1,9),
+                                      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),
+              damask.mechanics.PK2(table.get(options.stress ).reshape(-1,3,3),
-                                   table.get(options.stress ).reshape(-1,3,3)).reshape(-1,9),
+                                   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:
+    damask.util.report(scriptName,name)
-    table = damask.ASCIItable(name = name,
-                              buffered = False)
-  except: continue
-  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 = {
+        for i,o in enumerate(['Min','Mid','Max']):
-            'tensor': {'dim': 9, 'shape': [3,3], 'labels':options.tensor, 'column': []},
+            table.add('eigvec{}({})'.format(o,tensor),v[:,:,i],
-          }
+                      scriptID+' '+' '.join(sys.argv[1:]))
-  errors  = []
+        
-  remarks = []
+    table.to_ASCII(sys.stdout if name is None else name)
-  
-  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)

python/damask/__init__.py

@@ -15,6 +15,7 @@ from .config      import Material         # noqa
 from .colormaps   import Colormap, Color  # noqa
 from .orientation import Symmetry, Lattice, Rotation, Orientation # noqa
 from .dadf5       import DADF5            # noqa
+from .dadf5       import DADF5 as Result  # noqa
 
 from .geom        import Geom             # noqa
 from .solver      import Solver           # noqa

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,11 +1,11 @@
 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.
-    
+
     Parameters
     ----------
     F : numpy.array of shape (:,3,3) or (3,3)
@@ -21,67 +21,10 @@ def Cauchy(F,P):
     return symmetric(sigma)
 
 
-def PK2(F,P):
-    """
-    Return 2. Piola-Kirchhoff stress calculated from 1. 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.
-
-    """
-    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
-
-
-def strain_tensor(F,t,m):
-    """
-    Return strain tensor calculated from deformation gradient.
-      
-    For details refer to https://en.wikipedia.org/wiki/Finite_strain_theory and
-    https://de.wikipedia.org/wiki/Verzerrungstensor
-    
-    Parameters
-    ----------
-    F : numpy.array of shape (:,3,3) or (3,3)
-      Deformation gradient.
-    t : {‘V’, ‘U’}
-      Type of the polar decomposition, ‘V’ for left stretch tensor and ‘U’ for right stretch tensor.
-    m : float
-      Order of the strain.
-
-    """
-    F_ = F.reshape((1,3,3)) if F.shape == (3,3) else F
-    if   t == 'V':
-        B   = np.matmul(F_,transpose(F_))
-        w,n = np.linalg.eigh(B)
-    elif t == 'U':
-        C   = np.matmul(transpose(F_),F_)
-        w,n = np.linalg.eigh(C)
-    
-    if   m > 0.0:
-        eps = 1.0/(2.0*abs(m)) * (+ np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n)) 
-                                  - np.broadcast_to(np.eye(3),[F_.shape[0],3,3]))
-    elif m < 0.0:
-        eps = 1.0/(2.0*abs(m)) * (- np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n))
-                                  + np.broadcast_to(np.eye(3),[F_.shape[0],3,3]))
-    else:
-        eps = np.matmul(n,np.einsum('ij,ikj->ijk',0.5*np.log(w),n))
-      
-    return eps.reshape((3,3)) if np.shape(F) == (3,3) else \
-           eps
-
-
 def deviatoric_part(x):
     """
     Return deviatoric part of a tensor.
-      
+
     Parameters
     ----------
     x : numpy.array of shape (:,3,3) or (3,3)
@@ -89,13 +32,151 @@ def deviatoric_part(x):
 
     """
     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)) 
+           x - np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),[x.shape[0],3,3]),spherical_part(x))
+
+
+def eigenvalues(x):
+    """
+    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
+    ----------
+    x : numpy.array of shape (:,3,3) or (3,3)
+      Symmetric tensor of which the eigenvalues are computed.
+
+    """
+    return np.linalg.eigvalsh(symmetric(x))
+
+
+def eigenvectors(x,RHS=False):
+    """
+    Return eigenvectors of a symmetric tensor.
+
+    The eigenvalues are sorted in ascending order of their associated eigenvalues.
+
+    Parameters
+    ----------
+    x : 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(x))
+
+    if RHS:
+        if np.shape(x) == (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(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 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 = eigenvalues(x)
+    return (w[0] - w[2])*0.5 if np.shape(x) == (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
+    ----------
+    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 symmetric(S)
+
+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 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 spherical_part(x,tensor=False):
     """
     Return spherical (hydrostatic) part of a tensor.
-    
+
     Parameters
     ----------
     x : numpy.array of shape (:,3,3) or (3,3)
@@ -113,42 +194,50 @@ def spherical_part(x,tensor=False):
             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):
+def strain_tensor(F,t,m):
     """
-    Return the Mises equivalent of a stress tensor.
+    Return strain tensor calculated from deformation gradient.
-    
+
+    For details refer to https://en.wikipedia.org/wiki/Finite_strain_theory and
+    https://de.wikipedia.org/wiki/Verzerrungstensor
+
     Parameters
     ----------
-    sigma : numpy.array of shape (:,3,3) or (3,3)
+    F : numpy.array of shape (:,3,3) or (3,3)
-      Symmetric stress tensor of which the von Mises equivalent is computed.
+      Deformation gradient.
+    t : {‘V’, ‘U’}
+      Type of the polar decomposition, ‘V’ for left stretch tensor and ‘U’ for right stretch tensor.
+    m : float
+      Order of the strain.
 
     """
-    s = deviatoric_part(sigma)
+    F_ = F.reshape((1,3,3)) if F.shape == (3,3) else F
-    return np.sqrt(3.0/2.0*(np.sum(s**2.0))) if np.shape(sigma) == (3,3) else \
+    if   t == 'V':
-           np.sqrt(3.0/2.0*np.einsum('ijk->i',s**2.0))
+        B   = np.matmul(F_,transpose(F_))
-    
+        w,n = np.linalg.eigh(B)
-    
+    elif t == 'U':
-def Mises_strain(epsilon):
+        C   = np.matmul(transpose(F_),F_)
-    """
+        w,n = np.linalg.eigh(C)
-    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.
 
-    """
+    if   m > 0.0:
-    s = deviatoric_part(epsilon)
+        eps = 1.0/(2.0*abs(m)) * (+ np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n))
-    return np.sqrt(2.0/3.0*(np.sum(s**2.0))) if np.shape(epsilon) == (3,3) else \
+                                  - np.broadcast_to(np.eye(3),[F_.shape[0],3,3]))
-           np.sqrt(2.0/3.0*np.einsum('ijk->i',s**2.0))
+    elif m < 0.0:
+        eps = 1.0/(2.0*abs(m)) * (- np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n))
+                                  + np.broadcast_to(np.eye(3),[F_.shape[0],3,3]))
+    else:
+        eps = np.matmul(n,np.einsum('ij,ikj->ijk',0.5*np.log(w),n))
+
+    return eps.reshape((3,3)) if np.shape(F) == (3,3) else \
+           eps
 
 
 def symmetric(x):
     """
     Return the symmetrized tensor.
-    
+
     Parameters
     ----------
     x : numpy.array of shape (:,3,3) or (3,3)
@@ -158,43 +247,10 @@ def symmetric(x):
     return (x+transpose(x))*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):
     """
     Return the transpose of a tensor.
-      
+
     Parameters
     ----------
     x : numpy.array of shape (:,3,3) or (3,3)
@@ -205,62 +261,23 @@ def transpose(x):
            np.transpose(x,(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):
     """
     Singular value decomposition.
-      
+
     Parameters
     ----------
     x : 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, 
+      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 \
         np.einsum('ijk,ikl->ijl',u,vh)
-    
+
     output = []
     if 'R' in requested:
         output.append(R)
@@ -268,5 +285,22 @@ def __polar_decomposition(x,requested):
         output.append(np.dot(x,R.T) if np.shape(x) == (3,3) else np.einsum('ijk,ilk->ijl',x,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))
-    
+
     return tuple(output)
+
+
+def __Mises(x,s):
+    """
+    Base equation for Mises equivalent of a stres or strain tensor.
+
+    Parameters
+    ----------
+    x : 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(x)
+    return np.sqrt(s*(np.sum(d**2.0))) if np.shape(x) == (3,3) else \
+           np.sqrt(s*np.einsum('ijk->i',d**2.0))

python/damask/util.py

@@ -3,14 +3,18 @@ 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:
     """
     ASCII Colors (Blender code).
-    
+
     https://svn.blender.org/svnroot/bf-blender/trunk/blender/build_files/scons/tools/bcolors.py
     http://stackoverflow.com/questions/287871/print-in-terminal-with-colors-using-python
     """
@@ -36,7 +40,7 @@ class bcolors:
         self.BOLD = ''
         self.UNDERLINE = ''
         self.CROSSOUT = ''
-    
+
 
 # -----------------------------
 def srepr(arg,glue = '\n'):
@@ -159,11 +163,29 @@ def progressBar(iteration, total, prefix='', bar_length=50):
 
   if iteration == total: sys.stderr.write('\n')
   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."""
-  
+
   def __init__(self,message):
     """
     Sets return message.
@@ -175,7 +197,7 @@ class return_message():
 
     """
     self.message = message
-  
+
   def __repr__(self):
     """Return message suitable for interactive shells."""
     return srepr(self.message)

python/tests/test_DADF5.py

@@ -40,7 +40,7 @@ class TestDADF5:
         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')
+        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
@@ -52,8 +52,8 @@ class TestDADF5:
         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),
+        in_memory = mechanics.Cauchy(default.read_dataset(loc['P'],0),
-                                     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)
 
@@ -73,6 +73,54 @@ class TestDADF5:
         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'),
@@ -81,6 +129,24 @@ class TestDADF5:
         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'),
@@ -88,3 +154,30 @@ class TestDADF5:
         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

@@ -2,187 +2,224 @@ import numpy as np
 from damask import mechanics
 
 class TestMechanics:
-   
+
     n = 1000
     c = np.random.randint(n)
-   
+
 
     def test_vectorize_Cauchy(self):
-         P = np.random.random((self.n,3,3))
+        P = np.random.random((self.n,3,3))
-         F = np.random.random((self.n,3,3))
+        F = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.Cauchy(F,P)[self.c],
+        assert np.allclose(mechanics.Cauchy(P,F)[self.c],
-                            mechanics.Cauchy(F[self.c],P[self.c]))
+                           mechanics.Cauchy(P[self.c],F[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))
-
 
     def test_vectorize_deviatoric_part(self):
-         x = np.random.random((self.n,3,3))
+        x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.deviatoric_part(x)[self.c],
+        assert np.allclose(mechanics.deviatoric_part(x)[self.c],
-                            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):
+    def test_vectorize_eigenvectors(self):
-         x = np.random.random((self.n,3,3))
+        x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.spherical_part(x,True)[self.c],
+        assert np.allclose(mechanics.eigenvectors(x)[self.c],
-                            mechanics.spherical_part(x[self.c],True))
+                           mechanics.eigenvectors(x[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_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_left_stretch(self):
-         x = np.random.random((self.n,3,3))
+        x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.left_stretch(x)[self.c],
+        assert np.allclose(mechanics.left_stretch(x)[self.c],
-                            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))
+        x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.right_stretch(x)[self.c],
+        assert np.allclose(mechanics.right_stretch(x)[self.c],
-                            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."""
+        """Ensure Cauchy stress is symmetrized 1. Piola-Kirchhoff stress for no deformation."""
-         P = np.random.random((self.n,3,3))
+        P = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.Cauchy(np.broadcast_to(np.eye(3),(self.n,3,3)),P),
+        assert np.allclose(mechanics.Cauchy(P,np.broadcast_to(np.eye(3),(self.n,3,3))),
-                            mechanics.symmetric(P))
+                           mechanics.symmetric(P))
+
 
     def test_polar_decomposition(self):
-         """F = RU = VR."""
+        """F = RU = VR."""
-         F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.random((self.n,3,3))
+        F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.random((self.n,3,3))
-         R = mechanics.rotational_part(F)
+        R = mechanics.rotational_part(F)
-         V = mechanics.left_stretch(F)
+        V = mechanics.left_stretch(F)
-         U = mechanics.right_stretch(F)
+        U = mechanics.right_stretch(F)
-         assert np.allclose(np.matmul(R,U),
+        assert np.allclose(np.matmul(R,U),
-                            np.matmul(V,R))
+                           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."""
+        """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))
+        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
+        m = np.random.random()*20.0-10.0
-         assert np.allclose(mechanics.strain_tensor(F,'U',m),
+        assert np.allclose(mechanics.strain_tensor(F,'U',m),
-                            mechanics.strain_tensor(F,'V',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."""
+        """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)
+        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
+        m = np.random.random()*5.0-2.5
-         assert np.allclose(np.linalg.det(mechanics.strain_tensor(F,'U',m)),
+        assert np.allclose(np.linalg.det(mechanics.strain_tensor(F,'U',m)),
-                            np.linalg.det(mechanics.strain_tensor(F,'V',m)))
+                           np.linalg.det(mechanics.strain_tensor(F,'V',m)))
 
     def test_strain_tensor_rotation(self):
-         """Ensure that pure rotation results in no strain."""
+        """Ensure that pure rotation results in no strain."""
-         F = mechanics.rotational_part(np.random.random((self.n,3,3)))
+        F = mechanics.rotational_part(np.random.random((self.n,3,3)))
-         t = ['V','U'][np.random.randint(0,2)]
+        t = ['V','U'][np.random.randint(0,2)]
-         m = np.random.random()*2.0 - 1.0
+        m = np.random.random()*2.0 - 1.0
-         assert np.allclose(mechanics.strain_tensor(F,t,m),
+        assert np.allclose(mechanics.strain_tensor(F,t,m),
-                            0.0)
+                           0.0)
-   
-    def test_rotation_determinant(self):
-         """
-         Ensure that the determinant of the rotational part is +- 1.
 
-         Should be +1, but random F might contain a reflection.
+    def test_rotation_determinant(self):
-         """
+        """
-         x = np.random.random((self.n,3,3))
+        Ensure that the determinant of the rotational part is +- 1.
-         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."""
+        """Ensure that full tensor is sum of spherical and deviatoric part."""
-         x = np.random.random((self.n,3,3))
+        x = np.random.random((self.n,3,3))
-         sph = mechanics.spherical_part(x,True)
+        sph = mechanics.spherical_part(x,True)
-         assert np.allclose(sph + mechanics.deviatoric_part(x),
+        assert np.allclose(sph + mechanics.deviatoric_part(x),
-                            x)
+                           x)
 
     def test_deviatoric_Mises(self):
-         """Ensure that Mises equivalent stress depends only on deviatoric part."""
+        """Ensure that Mises equivalent stress depends only on deviatoric part."""
-         x = np.random.random((self.n,3,3))
+        x = np.random.random((self.n,3,3))
-         full = mechanics.Mises_stress(x)
+        full = mechanics.Mises_stress(x)
-         dev  = mechanics.Mises_stress(mechanics.deviatoric_part(x))
+        dev  = mechanics.Mises_stress(mechanics.deviatoric_part(x))
-         assert np.allclose(full,
+        assert np.allclose(full,
-                            dev)
+                           dev)
 
     def test_spherical_mapping(self):
-         """Ensure that mapping to tensor is correct."""
+        """Ensure that mapping to tensor is correct."""
-         x = np.random.random((self.n,3,3))
+        x = np.random.random((self.n,3,3))
-         tensor = mechanics.spherical_part(x,True)
+        tensor = mechanics.spherical_part(x,True)
-         scalar = mechanics.spherical_part(x)
+        scalar = mechanics.spherical_part(x)
-         assert np.allclose(np.linalg.det(tensor),
+        assert np.allclose(np.linalg.det(tensor),
-                            scalar**3.0)
+                           scalar**3.0)
 
     def test_spherical_Mises(self):
-         """Ensure that Mises equivalent strrain of spherical strain is 0."""
+        """Ensure that Mises equivalent strrain of spherical strain is 0."""
-         x = np.random.random((self.n,3,3))
+        x = np.random.random((self.n,3,3))
-         sph = mechanics.spherical_part(x,True)
+        sph = mechanics.spherical_part(x,True)
-         assert np.allclose(mechanics.Mises_strain(sph),
+        assert np.allclose(mechanics.Mises_strain(sph),
-                            0.0)
+                           0.0)
 
     def test_symmetric(self):
-         """Ensure that a symmetric tensor is half of the sum of a tensor and its transpose."""
+        """Ensure that a symmetric tensor is half of the sum of a tensor and its transpose."""
-         x = np.random.random((self.n,3,3))
+        x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.symmetric(x)*2.0,
+        assert np.allclose(mechanics.symmetric(x)*2.0,
-                            mechanics.transpose(x)+x)
+                           mechanics.transpose(x)+x)
 
 
     def test_transpose(self):
-         """Ensure that a symmetric tensor equals its transpose."""
+        """Ensure that a symmetric tensor equals its transpose."""
-         x = mechanics.symmetric(np.random.random((self.n,3,3)))
+        x = mechanics.symmetric(np.random.random((self.n,3,3)))
-         assert np.allclose(mechanics.transpose(x),
+        assert np.allclose(mechanics.transpose(x),
-                            x)
+                           x)
 
 
     def test_Mises(self):
-         """Ensure that equivalent stress is 3/2 of equivalent strain."""
+        """Ensure that equivalent stress is 3/2 of equivalent strain."""
-         x = np.random.random((self.n,3,3))
+        x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.Mises_stress(x)/mechanics.Mises_strain(x),
+        assert np.allclose(mechanics.Mises_stress(x)/mechanics.Mises_strain(x),
-                            1.5)
+                           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/constitutive_plastic_disloUCLA.f90

@@ -6,21 +6,10 @@
 !> @brief crystal plasticity model for bcc metals, especially Tungsten
 !--------------------------------------------------------------------------------------------------
 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,30 +35,30 @@ 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 :: &
+    character(len=pStringLen), allocatable, dimension(:) :: &
-      outputID                                                                                      !< ID of each post result output
+      output
     logical :: &
       dipoleFormation                                                                               !< flag indicating consideration of dipole formation
   end type                                                                                          !< container type for internal constitutive parameters
- 
+
   type :: tDisloUCLAState
     real(pReal), dimension(:,:), pointer :: &
       rho_mob, &
       rho_dip, &
       gamma_sl
   end type tDisloUCLAState
- 
+
   type :: tDisloUCLAdependentState
     real(pReal), dimension(:,:), allocatable :: &
       Lambda_sl, &
@@ -88,41 +77,36 @@ contains
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief module initialization
+!> @brief Perform module initialization.
 !> @details reads in material parameters, allocates arrays, and does sanity checks
 !--------------------------------------------------------------------------------------------------
 module subroutine plastic_disloUCLA_init
- 
+
   integer :: &
     Ninstance, &
     p, i, &
     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_DISLOUCLA_label//' init  -+>>>'; flush(6)
- 
+
-  write(6,'(/,a)')   ' <<<+-  plastic_'//PLASTICITY_DISLOUCLA_label//' init  -+>>>'
+  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) &
     write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
- 
+
   allocate(param(Ninstance))
   allocate(state(Ninstance))
   allocate(dotState(Ninstance))
   allocate(dependentState(Ninstance))
- 
+
- 
+
   do p = 1, size(phase_plasticity)
     if (phase_plasticity(p) /= PLASTICITY_DISLOUCLA_ID) cycle
     associate(prm => param(phase_plasticityInstance(p)), &
@@ -134,9 +118,9 @@ module subroutine plastic_disloUCLA_init
 !--------------------------------------------------------------------------------------------------
 !  optional parameters that need to be defined
     prm%mu = lattice_mu(p)
- 
+
     prm%aTol_rho = config%getFloat('atol_rho')
- 
+
     ! sanity checks
     if (prm%aTol_rho <= 0.0_pReal) extmsg = trim(extmsg)//' atol_rho'
 
@@ -147,7 +131,7 @@ module subroutine plastic_disloUCLA_init
     slipActive: if (prm%sum_N_sl > 0) then
       prm%Schmid = lattice_SchmidMatrix_slip(prm%N_sl,config%getString('lattice_structure'),&
                                              config%getFloat('c/a',defaultVal=0.0_pReal))
- 
+
       if(trim(config%getString('lattice_structure')) == 'bcc') then
         prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
                                                defaultVal = emptyRealArray)
@@ -158,20 +142,20 @@ module subroutine plastic_disloUCLA_init
         prm%nonSchmid_pos  = prm%Schmid
         prm%nonSchmid_neg  = prm%Schmid
       endif
- 
+
       prm%h_sl_sl     = lattice_interaction_SlipBySlip(prm%N_sl, &
                                                        config%getFloats('interaction_slipslip'), &
                                                        config%getString('lattice_structure'))
       prm%forestProjectionEdge = lattice_forestProjection_edge(prm%N_sl,config%getString('lattice_structure'),&
                                                                config%getFloat('c/a',defaultVal=0.0_pReal))
       prm%forestProjectionEdge = transpose(prm%forestProjectionEdge)
-      
+
       prm%rho_mob_0   = config%getFloats('rhoedge0',       requiredSize=size(prm%N_sl))
       prm%rho_dip_0   = config%getFloats('rhoedgedip0',    requiredSize=size(prm%N_sl))
       prm%v0          = config%getFloats('v0',             requiredSize=size(prm%N_sl))
       prm%b_sl        = config%getFloats('slipburgers',    requiredSize=size(prm%N_sl))
       prm%delta_F     = config%getFloats('qedge',          requiredSize=size(prm%N_sl))
- 
+
       prm%i_sl        = config%getFloats('clambdaslip',    requiredSize=size(prm%N_sl))
       prm%tau_0       = config%getFloats('tau_peierls',    requiredSize=size(prm%N_sl))
       prm%p           = config%getFloats('p_slip',         requiredSize=size(prm%N_sl), &
@@ -182,14 +166,14 @@ module subroutine plastic_disloUCLA_init
       prm%w           = config%getFloats('kink_width',     requiredSize=size(prm%N_sl))
       prm%omega       = config%getFloats('omega',          requiredSize=size(prm%N_sl))
       prm%B           = config%getFloats('friction_coeff', requiredSize=size(prm%N_sl))
- 
+
       prm%D                   = config%getFloat('grainsize')
       prm%D_0                 = config%getFloat('d0')
       prm%Q_cl                = config%getFloat('qsd')
       prm%atomicVolume        = config%getFloat('catomicvolume')       * prm%b_sl**3.0_pReal
       prm%D_a                 = config%getFloat('cedgedipmindistance') * prm%b_sl
       prm%dipoleformation     = config%getFloat('dipoleformationfactor') > 0.0_pReal !should be on by default, ToDo: change to /key/-type key
- 
+
       ! expand: family => system
       prm%rho_mob_0      = math_expand(prm%rho_mob_0,      prm%N_sl)
       prm%rho_dip_0      = math_expand(prm%rho_dip_0,      prm%N_sl)
@@ -206,8 +190,8 @@ module subroutine plastic_disloUCLA_init
       prm%i_sl           = math_expand(prm%i_sl,           prm%N_sl)
       prm%atomicVolume   = math_expand(prm%atomicVolume,   prm%N_sl)
       prm%D_a            = math_expand(prm%D_a,            prm%N_sl)
- 
+
- 
+
       ! sanity checks
       if (    prm%D_0            <= 0.0_pReal)  extmsg = trim(extmsg)//' D_0'
       if (    prm%Q_cl           <= 0.0_pReal)  extmsg = trim(extmsg)//' Q_cl'
@@ -219,7 +203,7 @@ module subroutine plastic_disloUCLA_init
       if (any(prm%tau_0          <  0.0_pReal)) extmsg = trim(extmsg)//' tau_0'
       if (any(prm%D_a            <= 0.0_pReal)) extmsg = trim(extmsg)//' cedgedipmindistance or slipb_sl'
       if (any(prm%atomicVolume   <= 0.0_pReal)) extmsg = trim(extmsg)//' catomicvolume or slipb_sl'
- 
+
     else slipActive
       allocate(prm%rho_mob_0(0))
       allocate(prm%rho_dip_0(0))
@@ -232,39 +216,14 @@ module subroutine plastic_disloUCLA_init
 
 !--------------------------------------------------------------------------------------------------
 !  output pararameters
-    outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
+    prm%output = 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
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
     NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
     sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl
     sizeState = sizeDotState
- 
+
     call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0)
 
 !--------------------------------------------------------------------------------------------------
@@ -275,14 +234,14 @@ module subroutine plastic_disloUCLA_init
     stt%rho_mob= spread(prm%rho_mob_0,2,NipcMyPhase)
     dot%rho_mob=>plasticState(p)%dotState(startIndex:endIndex,:)
     plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_rho
- 
+
     startIndex = endIndex + 1
     endIndex   = endIndex + prm%sum_N_sl
     stt%rho_dip=>plasticState(p)%state(startIndex:endIndex,:)
     stt%rho_dip= spread(prm%rho_dip_0,2,NipcMyPhase)
     dot%rho_dip=>plasticState(p)%dotState(startIndex:endIndex,:)
     plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_rho
- 
+
     startIndex = endIndex + 1
     endIndex   = endIndex + prm%sum_N_sl
     stt%gamma_sl=>plasticState(p)%state(startIndex:endIndex,:)
@@ -290,21 +249,21 @@ module subroutine plastic_disloUCLA_init
     plasticState(p)%aTolState(startIndex:endIndex) = 1.0e6_pReal                                    ! Don't use for convergence check
     ! global alias
     plasticState(p)%slipRate        => plasticState(p)%dotState(startIndex:endIndex,:)
- 
+
     allocate(dst%Lambda_sl(prm%sum_N_sl,NipcMyPhase),         source=0.0_pReal)
     allocate(dst%threshold_stress(prm%sum_N_sl,NipcMyPhase),  source=0.0_pReal)
- 
+
     plasticState(p)%state0 = plasticState(p)%state                                                  ! ToDo: this could be done centrally
- 
+
     end associate
- 
+
   enddo
 
 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)
@@ -312,7 +271,7 @@ pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp, &
     Lp                                                                                              !< plastic velocity gradient
   real(pReal), dimension(3,3,3,3), intent(out) :: &
     dLp_dMp                                                                                         !< derivative of Lp with respect to the Mandel stress
- 
+
   real(pReal), dimension(3,3), intent(in) :: &
     Mp                                                                                              !< Mandel stress
   real(pReal),                 intent(in) :: &
@@ -320,18 +279,18 @@ pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp, &
   integer,                     intent(in) :: &
     instance, &
     of
- 
+
   integer :: &
     i,k,l,m,n
   real(pReal), dimension(param(instance)%sum_N_sl) :: &
     dot_gamma_pos,dot_gamma_neg, &
     ddot_gamma_dtau_pos,ddot_gamma_dtau_neg
- 
+
   Lp = 0.0_pReal
   dLp_dMp = 0.0_pReal
- 
+
   associate(prm => param(instance))
- 
+
   call kinetics(Mp,T,instance,of,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg)
   do i = 1, prm%sum_N_sl
     Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%Schmid(1:3,1:3,i)
@@ -340,17 +299,17 @@ pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp, &
                        + ddot_gamma_dtau_pos(i) * prm%Schmid(k,l,i) * prm%nonSchmid_pos(m,n,i) &
                        + ddot_gamma_dtau_neg(i) * prm%Schmid(k,l,i) * prm%nonSchmid_neg(m,n,i)
   enddo
- 
+
   end associate
 
 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)
- 
+
   real(pReal), dimension(3,3),  intent(in) :: &
     Mp                                                                                               !< Mandel stress
   real(pReal),                  intent(in) :: &
@@ -358,7 +317,7 @@ module subroutine plastic_disloUCLA_dotState(Mp,T,instance,of)
   integer,                      intent(in) :: &
     instance, &
     of
- 
+
   real(pReal) :: &
     VacancyDiffusion
   real(pReal), dimension(param(instance)%sum_N_sl) :: &
@@ -369,16 +328,16 @@ module subroutine plastic_disloUCLA_dotState(Mp,T,instance,of)
     dot_rho_dip_formation, &
     dot_rho_dip_climb, &
     dip_distance
- 
+
   associate(prm => param(instance), stt => state(instance),dot => dotState(instance), dst => dependentState(instance))
- 
+
   call kinetics(Mp,T,instance,of,&
                 gdot_pos,gdot_neg, &
                 tau_pos_out = tau_pos,tau_neg_out = tau_neg)
- 
+
   dot%gamma_sl(:,of) = (gdot_pos+gdot_neg)                                                          ! ToDo: needs to be abs
   VacancyDiffusion = prm%D_0*exp(-prm%Q_cl/(kB*T))
- 
+
   where(dEq0(tau_pos))                                                                              ! ToDo: use avg of pos and neg
     dot_rho_dip_formation = 0.0_pReal
     dot_rho_dip_climb     = 0.0_pReal
@@ -393,73 +352,73 @@ module subroutine plastic_disloUCLA_dotState(Mp,T,instance,of)
                   * (1.0_pReal/(dip_distance+prm%D_a))
     dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,of))/(dip_distance-prm%D_a)               ! ToDo: Discuss with Franz: Stress dependency?
   end where
- 
+
   dot%rho_mob(:,of) = abs(dot%gamma_sl(:,of))/(prm%b_sl*dst%Lambda_sl(:,of)) &                      ! multiplication
                     - dot_rho_dip_formation &
                     - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_mob(:,of)*abs(dot%gamma_sl(:,of))        ! Spontaneous annihilation of 2 single edge dislocations
   dot%rho_dip(:,of) = dot_rho_dip_formation &
                     - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_dip(:,of)*abs(dot%gamma_sl(:,of)) &      ! Spontaneous annihilation of a single edge dislocation with a dipole constituent
                     - dot_rho_dip_climb
- 
+
   end associate
 
 end subroutine plastic_disloUCLA_dotState
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief calculates derived quantities from state
+!> @brief Calculate derived quantities from state.
 !--------------------------------------------------------------------------------------------------
 module subroutine plastic_disloUCLA_dependentState(instance,of)
- 
+
   integer,      intent(in) :: &
     instance, &
     of
- 
+
   real(pReal), dimension(param(instance)%sum_N_sl) :: &
     dislocationSpacing
- 
+
   associate(prm => param(instance), stt => state(instance),dst => dependentState(instance))
- 
+
   dislocationSpacing = sqrt(matmul(prm%forestProjectionEdge, &
                                    stt%rho_mob(:,of)+stt%rho_dip(:,of)))
   dst%threshold_stress(:,of) = prm%mu*prm%b_sl &
                              * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,of)+stt%rho_dip(:,of)))
- 
+
   dst%Lambda_sl(:,of) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl)
- 
+
   end associate
 
 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)
 
   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)
+  outputsLoop: do o = 1,size(prm%output)
-    select case(prm%outputID(o))
+    select case(trim(prm%output(o)))
-      case (rho_mob_ID)
+      case('edge_density')                                                                          ! ToDo: should be rho_mob
-        call results_writeDataset(group,stt%rho_mob,'rho_mob',&
+        if(prm%sum_N_sl>0) call results_writeDataset(group,stt%rho_mob,'rho_mob',&
-                                 'mobile dislocation density','1/m²')
+                                                     'mobile dislocation density','1/m²')
-      case (rho_dip_ID)
+      case('dipole_density')                                                                        ! ToDo: should be rho_dip
-        call results_writeDataset(group,stt%rho_dip,'rho_dip',&
+        if(prm%sum_N_sl>0) call results_writeDataset(group,stt%rho_dip,'rho_dip',&
-                                  'dislocation dipole density''1/m²')
+                                                     'dislocation dipole density''1/m²')
-      case (dot_gamma_sl_ID)
+      case('shear_rate_slip')                                                                       ! should be gamma
-        call results_writeDataset(group,stt%gamma_sl,'dot_gamma_sl',&   ! this is not dot!!
+        if(prm%sum_N_sl>0) call results_writeDataset(group,stt%gamma_sl,'dot_gamma_sl',&   ! this is not dot!!
-                                  'plastic shear','1')
+                                                     'plastic shear','1')
-      case (Lambda_sl_ID)
+      case('mfp_slip')                                                                              !ToDo: should be Lambda
-        call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',&
+        if(prm%sum_N_sl>0) call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',&
-                                  'mean free path for slip','m')
+                                                     'mean free path for slip','m')
-      case (tau_pass_ID)
+      case('threshold_stress_slip')                                                                 !ToDo: should be tau_pass
-        call results_writeDataset(group,dst%threshold_stress,'tau_pass',&
+        if(prm%sum_N_sl>0) call results_writeDataset(group,dst%threshold_stress,'tau_pass',&
-                                  'threshold stress for slip','Pa')
+                                                     'threshold stress for slip','Pa')
     end select
   enddo outputsLoop
   end associate
@@ -468,15 +427,15 @@ end subroutine plastic_disloUCLA_results
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Shear rates on slip systems, their derivatives with respect to resolved stress and the
+!> @brief Calculate shear rates on slip systems, their derivatives with respect to resolved
-!  resolved stresss
+!         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
 !--------------------------------------------------------------------------------------------------
 pure subroutine kinetics(Mp,T,instance,of, &
                  dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg,tau_pos_out,tau_neg_out)
- 
+
   real(pReal), dimension(3,3),  intent(in) :: &
     Mp                                                                                              !< Mandel stress
   real(pReal),                  intent(in) :: &
@@ -484,7 +443,7 @@ pure subroutine kinetics(Mp,T,instance,of, &
   integer,                      intent(in) :: &
     instance, &
     of
- 
+
   real(pReal),                  intent(out), dimension(param(instance)%sum_N_sl) :: &
     dot_gamma_pos, &
     dot_gamma_neg
@@ -501,82 +460,82 @@ pure subroutine kinetics(Mp,T,instance,of, &
     t_n, t_k, dtk,dtn, &
     needsGoodName                                                                                   ! ToDo: @Karo: any idea?
   integer :: j
- 
+
   associate(prm => param(instance), stt => state(instance), dst => dependentState(instance))
- 
+
   do j = 1, prm%sum_N_sl
     tau_pos(j) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,j))
     tau_neg(j) = math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,j))
   enddo
- 
+
- 
+
   if (present(tau_pos_out)) tau_pos_out = tau_pos
   if (present(tau_neg_out)) tau_neg_out = tau_neg
- 
+
   associate(BoltzmannRatio  => prm%delta_F/(kB*T), &
             dot_gamma_0     => stt%rho_mob(:,of)*prm%b_sl*prm%v0, &
             effectiveLength => dst%Lambda_sl(:,of) - prm%w)
- 
+
   significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check)
     StressRatio = (abs(tau_pos)-dst%threshold_stress(:,of))/prm%tau_0
     StressRatio_p       = StressRatio** prm%p
     StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal)
     needsGoodName       = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)
-    
+
     t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength)
     t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos)
- 
+
     vel = prm%kink_height/(t_n + t_k)
- 
+
     dot_gamma_pos = dot_gamma_0 * sign(vel,tau_pos) * 0.5_pReal
   else where significantPositiveTau
     dot_gamma_pos = 0.0_pReal
   end where significantPositiveTau
- 
+
   if (present(ddot_gamma_dtau_pos)) then
     significantPositiveTau2: where(abs(tau_pos)-dst%threshold_stress(:,of) > tol_math_check)
       dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) &
           * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_0
       dtk = -1.0_pReal * t_k / tau_pos
-     
+
       dvel = -1.0_pReal * prm%kink_height * (dtk + dtn) / (t_n + t_k)**2.0_pReal
- 
+
       ddot_gamma_dtau_pos = dot_gamma_0 * dvel* 0.5_pReal
     else where significantPositiveTau2
       ddot_gamma_dtau_pos = 0.0_pReal
     end where significantPositiveTau2
   endif
- 
+
   significantNegativeTau: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check)
     StressRatio = (abs(tau_neg)-dst%threshold_stress(:,of))/prm%tau_0
     StressRatio_p       = StressRatio** prm%p
     StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal)
     needsGoodName       = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)
- 
+
     t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength)
     t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos)
- 
+
     vel = prm%kink_height/(t_n + t_k)
- 
+
     dot_gamma_neg = dot_gamma_0 * sign(vel,tau_neg) * 0.5_pReal
   else where significantNegativeTau
     dot_gamma_neg = 0.0_pReal
   end where significantNegativeTau
- 
+
   if (present(ddot_gamma_dtau_neg)) then
     significantNegativeTau2: where(abs(tau_neg)-dst%threshold_stress(:,of) > tol_math_check)
       dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) &
           * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_0
       dtk = -1.0_pReal * t_k / tau_neg
-     
+
       dvel = -1.0_pReal * prm%kink_height * (dtk + dtn) / (t_n + t_k)**2.0_pReal
-   
+
       ddot_gamma_dtau_neg = dot_gamma_0 * dvel * 0.5_pReal
     else where significantNegativeTau2
       ddot_gamma_dtau_neg = 0.0_pReal
     end where significantNegativeTau2
   end if
-  
+
   end associate
   end associate
 

src/constitutive_plastic_isotropic.f90

@@ -8,14 +8,7 @@
 !! untextured polycrystal
 !--------------------------------------------------------------------------------------------------
 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,12 +27,12 @@ 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
     real(pReal), pointer, dimension(:) :: &
       xi, &
@@ -56,50 +49,45 @@ 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  -+>>>'; flush(6)
- 
+
-  write(6,'(/,a)')   ' <<<+-  plastic_'//PLASTICITY_ISOTROPIC_label//' init  -+>>>'
+  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) &
     write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
- 
+
   allocate(param(Ninstance))
   allocate(state(Ninstance))
   allocate(dotState(Ninstance))
- 
+
   do p = 1, size(phase_plasticity)
     if (phase_plasticity(p) /= PLASTICITY_ISOTROPIC_ID) cycle
     associate(prm => param(phase_plasticityInstance(p)), &
               dot => dotState(phase_plasticityInstance(p)), &
               stt => state(phase_plasticityInstance(p)), &
               config => config_phase(p))
- 
+
 #ifdef DEBUG
     if  (p==material_phaseAt(debug_g,debug_e)) &
       prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e)
 #endif
- 
+
     prm%xi_0            = config%getFloat('tau0')
     prm%xi_inf          = config%getFloat('tausat')
     prm%dot_gamma_0     = config%getFloat('gdot0')
@@ -114,9 +102,9 @@ module subroutine plastic_isotropic_init
     prm%a               = config%getFloat('a')
     prm%aTol_xi         = config%getFloat('atol_flowstress',defaultVal=1.0_pReal)
     prm%aTol_gamma      = config%getFloat('atol_shear',     defaultVal=1.0e-6_pReal)
- 
+
     prm%dilatation      = config%keyExists('/dilatation/')
- 
+
 !--------------------------------------------------------------------------------------------------
 !  sanity checks
     extmsg = ''
@@ -128,79 +116,62 @@ module subroutine plastic_isotropic_init
     if (prm%M              <= 0.0_pReal) extmsg = trim(extmsg)//' m'
     if (prm%aTol_xi        <= 0.0_pReal) extmsg = trim(extmsg)//' atol_xi'
     if (prm%aTol_gamma     <= 0.0_pReal) extmsg = trim(extmsg)//' atol_shear'
- 
+
 !--------------------------------------------------------------------------------------------------
 !  exit if any parameter is out of range
     if (extmsg /= '') &
       call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_ISOTROPIC_label//')')
- 
+
 !--------------------------------------------------------------------------------------------------
 !  output pararameters
-    outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
+    prm%output = 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
- 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
     NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
     sizeDotState = size(['xi               ','accumulated_shear'])
     sizeState = sizeDotState
- 
+
     call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0)
- 
+
 !--------------------------------------------------------------------------------------------------
 ! locally defined state aliases and initialization of state0 and aTolState
     stt%xi  => plasticState(p)%state   (1,:)
     stt%xi  = prm%xi_0
     dot%xi  => plasticState(p)%dotState(1,:)
     plasticState(p)%aTolState(1) = prm%aTol_xi
- 
+
     stt%gamma  => plasticState(p)%state   (2,:)
     dot%gamma  => plasticState(p)%dotState(2,:)
     plasticState(p)%aTolState(2) = prm%aTol_gamma
     ! global alias
     plasticState(p)%slipRate        => plasticState(p)%dotState(2:2,:)
- 
+
     plasticState(p)%state0 = plasticState(p)%state                                                  ! ToDo: this could be done centrally
- 
+
     end associate
- 
+
   enddo
 
 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)
- 
+
   real(pReal), dimension(3,3),     intent(out) :: &
     Lp                                                                                              !< plastic velocity gradient
   real(pReal), dimension(3,3,3,3), intent(out) :: &
     dLp_dMp                                                                                         !< derivative of Lp with respect to the Mandel stress
- 
+
   real(pReal), dimension(3,3), intent(in) :: &
     Mp                                                                                              !< Mandel stress
   integer,                     intent(in) :: &
     instance, &
     of
- 
+
   real(pReal), dimension(3,3) :: &
     Mp_dev                                                                                          !< deviatoric part of the Mandel stress
   real(pReal) :: &
@@ -209,16 +180,16 @@ module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
     squarenorm_Mp_dev                                                                               !< square of the norm of the deviatoric part of the Mandel stress
   integer :: &
     k, l, m, n
- 
+
   associate(prm => param(instance), stt => state(instance))
- 
+
   Mp_dev = math_deviatoric33(Mp)
   squarenorm_Mp_dev = math_mul33xx33(Mp_dev,Mp_dev)
   norm_Mp_dev = sqrt(squarenorm_Mp_dev)
- 
+
   if (norm_Mp_dev > 0.0_pReal) then
     dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(of))) **prm%n
- 
+
     Lp = dot_gamma/prm%M * Mp_dev/norm_Mp_dev
 #ifdef DEBUG
     if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 &
@@ -240,42 +211,42 @@ module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
     Lp = 0.0_pReal
     dLp_dMp = 0.0_pReal
   end if
- 
+
   end associate
 
 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)
- 
+
   real(pReal), dimension(3,3), intent(out) :: &
     Li                                                                                              !< inleastic velocity gradient
   real(pReal), dimension(3,3,3,3), intent(out)  :: &
     dLi_dMi                                                                                         !< derivative of Li with respect to Mandel stress
- 
+
   real(pReal), dimension(3,3),   intent(in) :: &
-    Mi                                                                                              !< Mandel stress 
+    Mi                                                                                              !< Mandel stress
   integer,                       intent(in) :: &
     instance, &
     of
- 
+
   real(pReal) :: &
     tr                                                                                              !< trace of spherical part of Mandel stress (= 3 x pressure)
   integer :: &
     k, l, m, n
- 
+
   associate(prm => param(instance), stt => state(instance))
- 
+
   tr=math_trace33(math_spherical33(Mi))
 
   if (prm%dilatation .and. abs(tr) > 0.0_pReal) then                                                 ! no stress or J2 plasticity --> Li and its derivative are zero
     Li = math_I3 &
        * prm%dot_gamma_0/prm%M * (3.0_pReal*prm%M*stt%xi(of))**(-prm%n) &
        * tr * abs(tr)**(prm%n-1.0_pReal)
- 
+
 #ifdef DEBUG
     if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 &
         .and. (of == prm%of_debug .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0)) then
@@ -292,38 +263,38 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of)
     Li      = 0.0_pReal
     dLi_dMi = 0.0_pReal
   endif
- 
+
   end associate
 
  end subroutine plastic_isotropic_LiAndItsTangent
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief calculates the rate of change of microstructure
+!> @brief Calculate the rate of change of microstructure.
 !--------------------------------------------------------------------------------------------------
 module subroutine plastic_isotropic_dotState(Mp,instance,of)
- 
+
   real(pReal), dimension(3,3),  intent(in) :: &
     Mp                                                                                              !< Mandel stress
   integer,                      intent(in) :: &
     instance, &
     of
- 
+
   real(pReal) :: &
     dot_gamma, &                                                                                    !< strainrate
     xi_inf_star, &                                                                                  !< saturation xi
     norm_Mp                                                                                         !< norm of the (deviatoric) Mandel stress
- 
+
   associate(prm => param(instance), stt => state(instance), dot => dotState(instance))
- 
+
   if (prm%dilatation) then
     norm_Mp = sqrt(math_mul33xx33(Mp,Mp))
   else
     norm_Mp = sqrt(math_mul33xx33(math_deviatoric33(Mp),math_deviatoric33(Mp)))
   endif
- 
+
   dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp /(prm%M*stt%xi(of))) **prm%n
- 
+
   if (dot_gamma > 1e-12_pReal) then
     if (dEq0(prm%c_1)) then
       xi_inf_star = prm%xi_inf
@@ -339,28 +310,28 @@ module subroutine plastic_isotropic_dotState(Mp,instance,of)
   else
     dot%xi(of) = 0.0_pReal
   endif
- 
+
   dot%gamma(of) = dot_gamma                                                                         ! ToDo: not really used
- 
+
   end associate
 
 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)
+  outputsLoop: do o = 1,size(prm%output)
-    select case(prm%outputID(o))
+    select case(trim(prm%output(o)))
-      case (xi_ID)
+      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,21 +22,21 @@ 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(:) :: &
+    character(len=pStringLen), allocatable, dimension(:) :: &
-      outputID                                                                                      !< ID of each post result output
+      output
   end type tParameters
- 
+
   type :: tKinehardeningState
     real(pReal), pointer, dimension(:,:) :: &                                                       !< vectors along NipcMyInstance
       crss, &                                                                                       !< critical resolved stress
@@ -59,7 +46,7 @@ submodule(constitutive) plastic_kinehardening
       gamma0, &                                                                                     !< accumulated shear at last switch of stress sense
       accshear                                                                                      !< accumulated (absolute) shear
   end type tKinehardeningState
- 
+
 !--------------------------------------------------------------------------------------------------
 ! containers for parameters and state
   type(tParameters),         allocatable, dimension(:) :: param
@@ -72,37 +59,32 @@ 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  -+>>>'; flush(6)
- 
+
-  write(6,'(/,a)')   ' <<<+-  plastic_'//PLASTICITY_KINEHARDENING_label//' init  -+>>>'
- 
   Ninstance = count(phase_plasticity == PLASTICITY_KINEHARDENING_ID)
   if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
     write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
- 
+
   allocate(param(Ninstance))
   allocate(state(Ninstance))
   allocate(dotState(Ninstance))
   allocate(deltaState(Ninstance))
- 
+
   do p = 1, size(phase_plasticityInstance)
     if (phase_plasticity(p) /= PLASTICITY_KINEHARDENING_ID) cycle
     associate(prm => param(phase_plasticityInstance(p)), &
@@ -110,22 +92,22 @@ module subroutine plastic_kinehardening_init
               dlt => deltaState(phase_plasticityInstance(p)), &
               stt => state(phase_plasticityInstance(p)),&
               config => config_phase(p))
- 
+
 #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
- 
+
 !--------------------------------------------------------------------------------------------------
 !  optional parameters that need to be defined
     prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal)
     prm%aTolShear      = config%getFloat('atol_shear',     defaultVal=1.0e-6_pReal)
- 
+
     ! sanity checks
     if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//' aTolresistance'
     if (prm%aTolShear      <= 0.0_pReal) extmsg = trim(extmsg)//' aTolShear'
- 
+
 !--------------------------------------------------------------------------------------------------
 ! slip related parameters
     prm%Nslip      = config%getInts('nslip',defaultVal=emptyIntArray)
@@ -133,7 +115,7 @@ module subroutine plastic_kinehardening_init
     slipActive: if (prm%totalNslip > 0) then
       prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),&
                                              config%getFloat('c/a',defaultVal=0.0_pReal))
- 
+
       if(trim(config%getString('lattice_structure')) == 'bcc') then
         prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
                                                defaultVal = emptyRealArray)
@@ -146,7 +128,7 @@ module subroutine plastic_kinehardening_init
       prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, &
                                                                 config%getFloats('interaction_slipslip'), &
                                                                 config%getString('lattice_structure'))
- 
+
       prm%crss0    = config%getFloats('crss0',    requiredSize=size(prm%Nslip))
       prm%tau1     = config%getFloats('tau1',     requiredSize=size(prm%Nslip))
       prm%tau1_b   = config%getFloats('tau1_b',   requiredSize=size(prm%Nslip))
@@ -154,10 +136,10 @@ module subroutine plastic_kinehardening_init
       prm%theta1   = config%getFloats('theta1',   requiredSize=size(prm%Nslip))
       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%gdot0    = config%getFloat('gdot0')
-      prm%n = config%getFloat('n_slip')
+      prm%n        = config%getFloat('n_slip')
- 
+
       ! expand: family => system
       prm%crss0    = math_expand(prm%crss0,   prm%Nslip)
       prm%tau1     = math_expand(prm%tau1,    prm%Nslip)
@@ -166,9 +148,9 @@ module subroutine plastic_kinehardening_init
       prm%theta1   = math_expand(prm%theta1,  prm%Nslip)
       prm%theta0_b = math_expand(prm%theta0_b,prm%Nslip)
       prm%theta1_b = math_expand(prm%theta1_b,prm%Nslip)
- 
+
- 
+
- 
+
 !--------------------------------------------------------------------------------------------------
 !  sanity checks
       if (    prm%gdot0  <= 0.0_pReal)   extmsg = trim(extmsg)//' gdot0'
@@ -176,58 +158,29 @@ module subroutine plastic_kinehardening_init
       if (any(prm%crss0  <= 0.0_pReal))  extmsg = trim(extmsg)//' crss0'
       if (any(prm%tau1   <= 0.0_pReal))  extmsg = trim(extmsg)//' tau1'
       if (any(prm%tau1_b <= 0.0_pReal))  extmsg = trim(extmsg)//' tau1_b'
- 
+
       !ToDo: Any sensible checks for theta?
- 
+
     endif slipActive
- 
+
 !--------------------------------------------------------------------------------------------------
 !  exit if any parameter is out of range
     if (extmsg /= '') &
       call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_KINEHARDENING_label//')')
- 
+
 !--------------------------------------------------------------------------------------------------
 !  output pararameters
-    outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
+    prm%output = 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
- 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
     NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
     sizeDotState   = size(['crss     ','crss_back', 'accshear ']) * prm%totalNslip
     sizeDeltaState = size(['sense ',   'chi0  ',    'gamma0'   ]) * prm%totalNslip
     sizeState = sizeDotState + sizeDeltaState
- 
+
     call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState)
- 
+
 !--------------------------------------------------------------------------------------------------
 ! locally defined state aliases and initialization of state0 and aTolState
     startIndex = 1
@@ -236,13 +189,13 @@ module subroutine plastic_kinehardening_init
     stt%crss = spread(prm%crss0, 2, NipcMyPhase)
     dot%crss => plasticState(p)%dotState(startIndex:endIndex,:)
     plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
- 
+
     startIndex = endIndex + 1
     endIndex   = endIndex + prm%totalNslip
     stt%crss_back => plasticState(p)%state   (startIndex:endIndex,:)
     dot%crss_back => plasticState(p)%dotState(startIndex:endIndex,:)
     plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
- 
+
     startIndex = endIndex + 1
     endIndex   = endIndex + prm%totalNslip
     stt%accshear => plasticState(p)%state   (startIndex:endIndex,:)
@@ -250,34 +203,34 @@ module subroutine plastic_kinehardening_init
     plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear
     ! global alias
     plasticState(p)%slipRate        => plasticState(p)%dotState(startIndex:endIndex,:)
- 
+
     o = plasticState(p)%offsetDeltaState
     startIndex = endIndex + 1
     endIndex   = endIndex + prm%totalNslip
     stt%sense => plasticState(p)%state     (startIndex  :endIndex  ,:)
     dlt%sense => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
- 
+
     startIndex = endIndex + 1
     endIndex   = endIndex +  prm%totalNslip
     stt%chi0 => plasticState(p)%state     (startIndex  :endIndex  ,:)
     dlt%chi0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
- 
+
     startIndex = endIndex + 1
     endIndex   = endIndex +  prm%totalNslip
     stt%gamma0 => plasticState(p)%state     (startIndex  :endIndex  ,:)
     dlt%gamma0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
- 
+
     plasticState(p)%state0 = plasticState(p)%state                                                  ! ToDo: this could be done centrally
- 
+
     end associate
- 
+
   enddo
 
 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)
 
@@ -285,26 +238,26 @@ pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,insta
     Lp                                                                                              !< plastic velocity gradient
   real(pReal), dimension(3,3,3,3), intent(out) :: &
     dLp_dMp                                                                                         !< derivative of Lp with respect to the Mandel stress
- 
+
   real(pReal), dimension(3,3), intent(in) :: &
     Mp                                                                                              !< Mandel stress
   integer,               intent(in) :: &
     instance, &
     of
- 
+
   integer :: &
     i,k,l,m,n
   real(pReal), dimension(param(instance)%totalNslip) :: &
     gdot_pos,gdot_neg, &
     dgdot_dtau_pos,dgdot_dtau_neg
- 
+
   Lp = 0.0_pReal
   dLp_dMp = 0.0_pReal
- 
+
   associate(prm => param(instance))
- 
+
   call kinetics(Mp,instance,of,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
- 
+
   do i = 1, prm%totalNslip
     Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%Schmid(1:3,1:3,i)
     forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@@ -312,14 +265,14 @@ pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,insta
                        + dgdot_dtau_pos(i) * prm%Schmid(k,l,i) * prm%nonSchmid_pos(m,n,i) &
                        + dgdot_dtau_neg(i) * prm%Schmid(k,l,i) * prm%nonSchmid_neg(m,n,i)
   enddo
- 
+
   end associate
 
 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)
 
@@ -328,40 +281,40 @@ module subroutine plastic_kinehardening_dotState(Mp,instance,of)
   integer,                      intent(in) :: &
     instance, &
     of
- 
+
   real(pReal) :: &
     sumGamma
   real(pReal), dimension(param(instance)%totalNslip) :: &
     gdot_pos,gdot_neg
- 
+
- 
+
   associate(prm => param(instance), stt => state(instance), dot => dotState(instance))
- 
+
   call kinetics(Mp,instance,of,gdot_pos,gdot_neg)
   dot%accshear(:,of) = abs(gdot_pos+gdot_neg)
   sumGamma = sum(stt%accshear(:,of))
- 
+
-  
+
   dot%crss(:,of) = matmul(prm%interaction_SlipSlip,dot%accshear(:,of)) &
                  * (  prm%theta1 &
                      + (prm%theta0 - prm%theta1 + prm%theta0*prm%theta1*sumGamma/prm%tau1) &
                      * exp(-sumGamma*prm%theta0/prm%tau1) &
                    )
-                     
+
   dot%crss_back(:,of) = stt%sense(:,of)*dot%accshear(:,of) * &
            ( prm%theta1_b + &
              (prm%theta0_b - prm%theta1_b &
                + prm%theta0_b*prm%theta1_b/(prm%tau1_b+stt%chi0(:,of))*(stt%accshear(:,of)-stt%gamma0(:,of))&
              ) *exp(-(stt%accshear(:,of)-stt%gamma0(:,of)) *prm%theta0_b/(prm%tau1_b+stt%chi0(:,of))) &
            )
- 
+
   end associate
 
 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)
 
@@ -370,18 +323,18 @@ module subroutine plastic_kinehardening_deltaState(Mp,instance,of)
   integer,                      intent(in) :: &
     instance, &
     of
- 
+
   real(pReal), dimension(param(instance)%totalNslip) :: &
     gdot_pos,gdot_neg, &
     sense
- 
+
   associate(prm => param(instance), stt => state(instance), dlt => deltaState(instance))
- 
+
   call kinetics(Mp,instance,of,gdot_pos,gdot_neg)
   sense = merge(state(instance)%sense(:,of), &                                                      ! keep existing...
                 sign(1.0_pReal,gdot_pos+gdot_neg), &                                                ! ...or have a defined
                 dEq0(gdot_pos+gdot_neg,1e-10_pReal))                                                ! current sense of shear direction
- 
+
 #ifdef DEBUG
   if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 &
              .and. (of == prm%of_debug &
@@ -390,7 +343,7 @@ module subroutine plastic_kinehardening_deltaState(Mp,instance,of)
     write(6,*) sense,state(instance)%sense(:,of)
   endif
 #endif
- 
+
 !--------------------------------------------------------------------------------------------------
 ! switch in sense of shear?
   where(dNeq(sense,stt%sense(:,of),0.1_pReal))
@@ -402,45 +355,43 @@ module subroutine plastic_kinehardening_deltaState(Mp,instance,of)
     dlt%chi0  (:,of) = 0.0_pReal
     dlt%gamma0(:,of) = 0.0_pReal
   end where
- 
+
   end associate
 
 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)
+  outputsLoop: do o = 1,size(prm%output)
-    select case(prm%outputID(o))
+    select case(trim(prm%output(o)))
-     case (crss_ID)
+     case('resistance')
-       call results_writeDataset(group,stt%crss,'xi_sl', &
+       if(prm%totalNslip>0) call results_writeDataset(group,stt%crss,'xi_sl', &
-                                'resistance against plastic slip','Pa')
+                                                      'resistance against plastic slip','Pa')
-
+     case('backstress')                                                                             ! ToDo: should be 'tau_back'
-     case(crss_back_ID)
+       if(prm%totalNslip>0) call results_writeDataset(group,stt%crss_back,'tau_back', &
-       call results_writeDataset(group,stt%crss_back,'tau_back', &
+                                                      'back stress against plastic slip','Pa')
-                                'back stress against plastic slip','Pa')
+     case ('sense')
-                                
+       if(prm%totalNslip>0) call results_writeDataset(group,stt%sense,'sense_of_shear', &
-     case (sense_ID)
+                                                      'tbd','1')
-       call results_writeDataset(group,stt%sense,'sense_of_shear','tbd','1')
+     case ('chi0')
-
+       if(prm%totalNslip>0) call results_writeDataset(group,stt%chi0,'chi0', &
-     case (chi0_ID)
+                                                      'tbd','Pa')
-       call results_writeDataset(group,stt%chi0,'chi0','tbd','Pa')
+     case ('gamma0')
-       
+       if(prm%totalNslip>0) call results_writeDataset(group,stt%gamma0,'gamma0', &
-     case (gamma0_ID)
+                                                      'tbd','1')
-       call results_writeDataset(group,stt%gamma0,'gamma0','tbd','1')
+     case ('accumulatedshear')
-       
+       if(prm%totalNslip>0) call results_writeDataset(group,stt%accshear,'gamma_sl', &
-     case (accshear_ID)
+                                                      'plastic shear','1')
-       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
+!> @brief Calculate shear rates on slip systems and their derivatives with respect to resolved
-!> @details: Shear rates are calculated only optionally.
+!         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)
@@ -462,44 +414,44 @@ pure subroutine kinetics(Mp,instance,of, &
   integer,                intent(in) :: &
     instance, &
     of
- 
+
   real(pReal),                  intent(out), dimension(param(instance)%totalNslip) :: &
     gdot_pos, &
     gdot_neg
   real(pReal),                  intent(out), optional, dimension(param(instance)%totalNslip) :: &
     dgdot_dtau_pos, &
     dgdot_dtau_neg
- 
+
   real(pReal), dimension(param(instance)%totalNslip) :: &
     tau_pos, &
     tau_neg
   integer :: i
   logical :: nonSchmidActive
- 
+
   associate(prm => param(instance), stt => state(instance))
- 
+
   nonSchmidActive = size(prm%nonSchmidCoeff) > 0
- 
+
   do i = 1, prm%totalNslip
     tau_pos(i) =       math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,of)
     tau_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)) - stt%crss_back(i,of), &
                        0.0_pReal, nonSchmidActive)
   enddo
- 
+
   where(dNeq0(tau_pos))
     gdot_pos = prm%gdot0 * merge(0.5_pReal,1.0_pReal, nonSchmidActive) &                            ! 1/2 if non-Schmid active
              * sign(abs(tau_pos/stt%crss(:,of))**prm%n,  tau_pos)
   else where
     gdot_pos = 0.0_pReal
   end where
- 
+
   where(dNeq0(tau_neg))
     gdot_neg = prm%gdot0 * 0.5_pReal &                                                              ! only used if non-Schmid active, always 1/2
              * sign(abs(tau_neg/stt%crss(:,of))**prm%n,  tau_neg)
   else where
     gdot_neg = 0.0_pReal
   end where
- 
+
   if (present(dgdot_dtau_pos)) then
     where(dNeq0(gdot_pos))
       dgdot_dtau_pos = gdot_pos*prm%n/tau_pos

src/constitutive_plastic_none.f90

@@ -18,19 +18,19 @@ module subroutine plastic_none_init
     Ninstance, &
     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) &
     write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
- 
+
   do p = 1, size(phase_plasticity)
     if (phase_plasticity(p) /= PLASTICITY_NONE_ID) cycle
- 
+
     NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
     call material_allocatePlasticState(p,NipcMyPhase,0,0,0)
- 
+
   enddo
 
 end subroutine plastic_none_init

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)
+    prm%output = 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
 
 !--------------------------------------------------------------------------------------------------
 ! 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)
+  outputsLoop: do o = 1,size(prm%output)
-    select case(prm%outputID(o))
+    select case(trim(prm%output(o)))
-      case (rho_sgl_mob_edg_pos_ID)
+      case('rho_sgl_mob_edg_pos')
-        call results_writeDataset(group,stt%rho_sgl_mob_edg_pos, '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²')
+                                                       'positive mobile edge density','1/m²')
-      case (rho_sgl_imm_edg_pos_ID)
+      case('rho_sgl_imm_edg_pos')
-        call results_writeDataset(group,stt%rho_sgl_imm_edg_pos, '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²')
+                                                       'positive immobile edge density','1/m²')
-      case (rho_sgl_mob_edg_neg_ID)
+      case('rho_sgl_mob_edg_neg')
-        call results_writeDataset(group,stt%rho_sgl_mob_edg_neg, '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²')
+                                                       'negative mobile edge density','1/m²')
-      case (rho_sgl_imm_edg_neg_ID)
+      case('rho_sgl_imm_edg_neg')
-        call results_writeDataset(group,stt%rho_sgl_imm_edg_neg, '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²')
+                                                       'negative immobile edge density','1/m²')
-      case (rho_dip_edg_ID)
+      case('rho_dip_edg')
-        call results_writeDataset(group,stt%rho_dip_edg, 'rho_dip_edg',&
+        if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_dip_edg, 'rho_dip_edg',&
-                                  'edge dipole density','1/m²')
+                                                       'edge dipole density','1/m²')
-      case (rho_sgl_mob_scr_pos_ID)
+      case('rho_sgl_mob_scr_pos')
-        call results_writeDataset(group,stt%rho_sgl_mob_scr_pos, '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²')
+                                                       'positive mobile screw density','1/m²')
-      case (rho_sgl_imm_scr_pos_ID)
+      case('rho_sgl_imm_scr_pos')
-        call results_writeDataset(group,stt%rho_sgl_imm_scr_pos, '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²')
+                                                       'positive immobile screw density','1/m²')
-      case (rho_sgl_mob_scr_neg_ID)
+      case('rho_sgl_mob_scr_neg')
-        call results_writeDataset(group,stt%rho_sgl_mob_scr_neg, '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²')
+                                                       'negative mobile screw density','1/m²')
-      case (rho_sgl_imm_scr_neg_ID)
+      case('rho_sgl_imm_scr_neg')
-        call results_writeDataset(group,stt%rho_sgl_imm_scr_neg, '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²')
+                                                       'negative immobile screw density','1/m²')
-      case (rho_dip_scr_ID)
+      case('rho_dip_scr')
-        call results_writeDataset(group,stt%rho_dip_scr, 'rho_dip_scr',&
+        if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_dip_scr, 'rho_dip_scr',&
-                                  'screw dipole density','1/m²')
+                                                       'screw dipole density','1/m²')
-      case (rho_forest_ID)
+      case('rho_forest')
-        call results_writeDataset(group,stt%rho_forest, 'rho_forest',&
+        if(prm%totalNslip>0) call results_writeDataset(group,stt%rho_forest, 'rho_forest',&
-                                  'forest density','1/m²')
+                                                       'forest density','1/m²')
-      case (v_edg_pos_ID)
+      case('v_edg_pos')
-        call results_writeDataset(group,stt%v_edg_pos, 'v_edg_pos',&
+        if(prm%totalNslip>0) call results_writeDataset(group,stt%v_edg_pos, 'v_edg_pos',&
-                                  'positive edge velocity','m/s')
+                                                       'positive edge velocity','m/s')
-      case (v_edg_neg_ID)
+      case('v_edg_neg')
-        call results_writeDataset(group,stt%v_edg_neg, 'v_edg_neg',&
+        if(prm%totalNslip>0) call results_writeDataset(group,stt%v_edg_neg, 'v_edg_neg',&
-                                  'negative edge velocity','m/s')
+                                                       'negative edge velocity','m/s')
-      case (v_scr_pos_ID)
+      case('v_scr_pos')
-        call results_writeDataset(group,stt%v_scr_pos, 'v_scr_pos',&
+        if(prm%totalNslip>0) call results_writeDataset(group,stt%v_scr_pos, 'v_scr_pos',&
-                                  'positive srew velocity','m/s')
+                                                       'positive srew velocity','m/s')
-      case (v_scr_neg_ID)
+      case('v_scr_neg')
-        call results_writeDataset(group,stt%v_scr_neg, 'v_scr_neg',&
+        if(prm%totalNslip>0) call results_writeDataset(group,stt%v_scr_neg, 'v_scr_neg',&
-                                  'negative screw velocity','m/s')
+                                                       'negative screw velocity','m/s')
-      case(gamma_ID)
+      case('gamma')
-        call results_writeDataset(group,stt%gamma,'gamma',&
+        if(prm%totalNslip>0) call results_writeDataset(group,stt%gamma,'gamma',&
-                                  'plastic shear','1')
+                                                       'plastic shear','1')
-      case (tau_pass_ID)
+      case('tau_pass')
-        call results_writeDataset(group,dst%tau_pass,'tau_pass',&
+        if(prm%totalNslip>0) call results_writeDataset(group,dst%tau_pass,'tau_pass',&
-                                  'passing stress for slip','Pa')
+                                                       '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,13 +44,13 @@ 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(:) :: &
+    character(len=pStringLen), allocatable, dimension(:) :: &
-      outputID                                                                                      !< ID of each post result output
+      output
   end type tParameters
- 
+
   type :: tPhenopowerlawState
     real(pReal), pointer, dimension(:,:) :: &
       xi_slip, &
@@ -71,7 +58,7 @@ submodule(constitutive) plastic_phenopowerlaw
       gamma_slip, &
       gamma_twin
   end type tPhenopowerlawState
- 
+
 !--------------------------------------------------------------------------------------------------
 ! containers for parameters and state
   type(tParameters),         allocatable, dimension(:) :: param
@@ -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,51 +78,46 @@ 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  -+>>>'; flush(6)
- 
+
-  write(6,'(/,a)')   ' <<<+-  plastic_'//PLASTICITY_PHENOPOWERLAW_label//' init  -+>>>'
- 
   Ninstance = count(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)
   if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) &
     write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
- 
+
   allocate(param(Ninstance))
   allocate(state(Ninstance))
   allocate(dotState(Ninstance))
- 
+
   do p = 1, size(phase_plasticity)
     if (phase_plasticity(p) /= PLASTICITY_PHENOPOWERLAW_ID) cycle
     associate(prm => param(phase_plasticityInstance(p)), &
               dot => dotState(phase_plasticityInstance(p)), &
               stt => state(phase_plasticityInstance(p)), &
               config => config_phase(p))
- 
+
 !--------------------------------------------------------------------------------------------------
 !  optional parameters that need to be defined
     prm%c_1            = config%getFloat('twin_c',defaultVal=0.0_pReal)
     prm%c_2            = config%getFloat('twin_b',defaultVal=1.0_pReal)
     prm%c_3            = config%getFloat('twin_e',defaultVal=0.0_pReal)
     prm%c_4            = config%getFloat('twin_d',defaultVal=0.0_pReal)
- 
+
     prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal)
     prm%aTolShear      = config%getFloat('atol_shear',     defaultVal=1.0e-6_pReal)
     prm%aTolTwinfrac   = config%getFloat('atol_twinfrac',  defaultVal=1.0e-6_pReal)
- 
+
     ! sanity checks
     if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//' aTolresistance'
     if (prm%aTolShear      <= 0.0_pReal) extmsg = trim(extmsg)//' aTolShear'
     if (prm%aTolTwinfrac   <= 0.0_pReal) extmsg = trim(extmsg)//' atoltwinfrac'
- 
+
 !--------------------------------------------------------------------------------------------------
 ! slip related parameters
     prm%Nslip      = config%getInts('nslip',defaultVal=emptyIntArray)
@@ -143,7 +125,7 @@ module subroutine plastic_phenopowerlaw_init
     slipActive: if (prm%totalNslip > 0) then
       prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),&
                                                   config%getFloat('c/a',defaultVal=0.0_pReal))
- 
+
       if(trim(config%getString('lattice_structure')) == 'bcc') then
         prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
                                               defaultVal = emptyRealArray)
@@ -157,22 +139,22 @@ module subroutine plastic_phenopowerlaw_init
       prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, &
                                                                 config%getFloats('interaction_slipslip'), &
                                                                 config%getString('lattice_structure'))
- 
+
       prm%xi_slip_0   = config%getFloats('tau0_slip',   requiredSize=size(prm%Nslip))
       prm%xi_slip_sat = config%getFloats('tausat_slip', requiredSize=size(prm%Nslip))
       prm%H_int       = config%getFloats('h_int',       requiredSize=size(prm%Nslip), &
                                          defaultVal=[(0.0_pReal,i=1,size(prm%Nslip))])
- 
+
       prm%gdot0_slip  = config%getFloat('gdot0_slip')
       prm%n_slip      = config%getFloat('n_slip')
       prm%a_slip      = config%getFloat('a_slip')
       prm%h0_SlipSlip = config%getFloat('h0_slipslip')
- 
+
       ! expand: family => system
       prm%xi_slip_0   = math_expand(prm%xi_slip_0,  prm%Nslip)
       prm%xi_slip_sat = math_expand(prm%xi_slip_sat,prm%Nslip)
       prm%H_int       = math_expand(prm%H_int,      prm%Nslip)
- 
+
       ! sanity checks
       if (    prm%gdot0_slip  <= 0.0_pReal)      extmsg = trim(extmsg)//' gdot0_slip'
       if (    prm%a_slip      <= 0.0_pReal)      extmsg = trim(extmsg)//' a_slip'
@@ -183,7 +165,7 @@ module subroutine plastic_phenopowerlaw_init
       allocate(prm%interaction_SlipSlip(0,0))
       allocate(prm%xi_slip_0(0))
     endif slipActive
- 
+
 !--------------------------------------------------------------------------------------------------
 ! twin related parameters
     prm%Ntwin      = config%getInts('ntwin', defaultVal=emptyIntArray)
@@ -196,17 +178,17 @@ module subroutine plastic_phenopowerlaw_init
                                                                 config%getString('lattice_structure'))
       prm%gamma_twin_char      = lattice_characteristicShear_twin(prm%Ntwin,config%getString('lattice_structure'),&
                                                              config%getFloat('c/a'))
- 
+
       prm%xi_twin_0            = config%getFloats('tau0_twin',requiredSize=size(prm%Ntwin))
- 
+
       prm%gdot0_twin           = config%getFloat('gdot0_twin')
       prm%n_twin               = config%getFloat('n_twin')
       prm%spr                  = config%getFloat('s_pr')
       prm%h0_TwinTwin          = config%getFloat('h0_twintwin')
- 
+
       ! expand: family => system
       prm%xi_twin_0   = math_expand(prm%xi_twin_0,  prm%Ntwin)
- 
+
       ! sanity checks
       if (prm%gdot0_twin <= 0.0_pReal)  extmsg = trim(extmsg)//' gdot0_twin'
       if (prm%n_twin     <= 0.0_pReal)  extmsg = trim(extmsg)//' n_twin'
@@ -215,7 +197,7 @@ module subroutine plastic_phenopowerlaw_init
       allocate(prm%xi_twin_0(0))
       allocate(prm%gamma_twin_char(0))
     endif twinActive
- 
+
 !--------------------------------------------------------------------------------------------------
 ! slip-twin related parameters
     slipAndTwinActive: if (prm%totalNslip > 0 .and. prm%totalNtwin > 0) then
@@ -231,63 +213,25 @@ module subroutine plastic_phenopowerlaw_init
       allocate(prm%interaction_TwinSlip(prm%TotalNtwin,prm%TotalNslip))                             ! at least one dimension is 0
       prm%h0_TwinSlip = 0.0_pReal
     endif slipAndTwinActive
- 
+
 !--------------------------------------------------------------------------------------------------
 !  exit if any parameter is out of range
     if (extmsg /= '') &
       call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_PHENOPOWERLAW_label//')')
- 
+
 !--------------------------------------------------------------------------------------------------
 !  output pararameters
-    outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
+    prm%output = 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
- 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
     NipcMyPhase = count(material_phaseAt == p) * discretization_nIP
     sizeDotState = size(['tau_slip  ','gamma_slip']) * prm%totalNslip &
                  + size(['tau_twin  ','gamma_twin']) * prm%totalNtwin
     sizeState = sizeDotState
- 
+
     call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0)
- 
+
 !--------------------------------------------------------------------------------------------------
 ! locally defined state aliases and initialization of state0 and aTolState
     startIndex = 1
@@ -296,14 +240,14 @@ module subroutine plastic_phenopowerlaw_init
     stt%xi_slip = spread(prm%xi_slip_0, 2, NipcMyPhase)
     dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:)
     plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
- 
+
     startIndex = endIndex + 1
     endIndex   = endIndex + prm%totalNtwin
     stt%xi_twin => plasticState(p)%state   (startIndex:endIndex,:)
     stt%xi_twin = spread(prm%xi_twin_0, 2, NipcMyPhase)
     dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:)
     plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
- 
+
     startIndex = endIndex + 1
     endIndex   = endIndex + prm%totalNslip
     stt%gamma_slip => plasticState(p)%state   (startIndex:endIndex,:)
@@ -317,18 +261,18 @@ module subroutine plastic_phenopowerlaw_init
     stt%gamma_twin => plasticState(p)%state   (startIndex:endIndex,:)
     dot%gamma_twin => plasticState(p)%dotState(startIndex:endIndex,:)
     plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear
- 
+
     plasticState(p)%state0 = plasticState(p)%state                                                  ! ToDo: this could be done centrally
- 
+
     end associate
- 
+
   enddo
 
 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
 !--------------------------------------------------------------------------------------------------
@@ -338,13 +282,13 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta
     Lp                                                                                              !< plastic velocity gradient
   real(pReal), dimension(3,3,3,3), intent(out) :: &
     dLp_dMp                                                                                         !< derivative of Lp with respect to the Mandel stress
- 
+
   real(pReal), dimension(3,3), intent(in) :: &
     Mp                                                                                              !< Mandel stress
   integer,               intent(in) :: &
     instance, &
     of
- 
+
   integer :: &
     i,k,l,m,n
   real(pReal), dimension(param(instance)%totalNslip) :: &
@@ -352,12 +296,12 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta
     dgdot_dtauslip_pos,dgdot_dtauslip_neg
   real(pReal), dimension(param(instance)%totalNtwin) :: &
     gdot_twin,dgdot_dtautwin
- 
+
   Lp = 0.0_pReal
   dLp_dMp = 0.0_pReal
- 
+
   associate(prm => param(instance))
- 
+
   call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg)
   slipSystems: do i = 1, prm%totalNslip
     Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i)
@@ -366,7 +310,7 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta
                        + dgdot_dtauslip_pos(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_pos(m,n,i) &
                        + dgdot_dtauslip_neg(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_neg(m,n,i)
   enddo slipSystems
- 
+
   call kinetics_twin(Mp,instance,of,gdot_twin,dgdot_dtautwin)
   twinSystems: do i = 1, prm%totalNtwin
     Lp = Lp + gdot_twin(i)*prm%Schmid_twin(1:3,1:3,i)
@@ -374,14 +318,14 @@ pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,insta
       dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
                        + dgdot_dtautwin(i)*prm%Schmid_twin(k,l,i)*prm%Schmid_twin(m,n,i)
   enddo twinSystems
- 
+
   end associate
 
 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)
 
@@ -390,7 +334,7 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
   integer,                      intent(in) :: &
     instance, &
     of
- 
+
   real(pReal) :: &
     c_SlipSlip,c_TwinSlip,c_TwinTwin, &
     xi_slip_sat_offset,&
@@ -398,9 +342,9 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
   real(pReal), dimension(param(instance)%totalNslip) :: &
     left_SlipSlip,right_SlipSlip, &
     gdot_slip_pos,gdot_slip_neg
- 
+
   associate(prm => param(instance), stt => state(instance), dot => dotState(instance))
- 
+
   sumGamma = sum(stt%gamma_slip(:,of))
   sumF     = sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)
 
@@ -409,26 +353,26 @@ module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
   c_SlipSlip = prm%h0_slipslip * (1.0_pReal + prm%c_1*sumF** prm%c_2)
   c_TwinSlip = prm%h0_TwinSlip * sumGamma**prm%c_3
   c_TwinTwin = prm%h0_TwinTwin * sumF**prm%c_4
- 
+
 !--------------------------------------------------------------------------------------------------
 !  calculate left and right vectors
   left_SlipSlip  = 1.0_pReal + prm%H_int
   xi_slip_sat_offset = prm%spr*sqrt(sumF)
   right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) **prm%a_slip &
                  * sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset))
- 
+
 !--------------------------------------------------------------------------------------------------
 ! shear rates
   call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg)
   dot%gamma_slip(:,of) = abs(gdot_slip_pos+gdot_slip_neg)
   call kinetics_twin(Mp,instance,of,dot%gamma_twin(:,of))
- 
+
 !--------------------------------------------------------------------------------------------------
 ! hardening
   dot%xi_slip(:,of) = c_SlipSlip * left_SlipSlip * &
                       matmul(prm%interaction_SlipSlip,dot%gamma_slip(:,of)*right_SlipSlip) &
                     + matmul(prm%interaction_SlipTwin,dot%gamma_twin(:,of))
- 
+
   dot%xi_twin(:,of) = c_TwinSlip * matmul(prm%interaction_TwinSlip,dot%gamma_slip(:,of)) &
                     + c_TwinTwin * matmul(prm%interaction_TwinTwin,dot%gamma_twin(:,of))
   end associate
@@ -437,33 +381,33 @@ 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)
 
   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)
+  outputsLoop: do o = 1,size(prm%output)
-    select case(prm%outputID(o))
+    select case(trim(prm%output(o)))
+
+      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')
+        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')
 
-      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_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')
-        
     end select
   enddo outputsLoop
   end associate
@@ -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)
@@ -485,44 +430,44 @@ pure subroutine kinetics_slip(Mp,instance,of, &
   integer,                      intent(in) :: &
     instance, &
     of
- 
+
   real(pReal),                  intent(out), dimension(param(instance)%totalNslip) :: &
     gdot_slip_pos, &
     gdot_slip_neg
   real(pReal),                  intent(out), optional, dimension(param(instance)%totalNslip) :: &
     dgdot_dtau_slip_pos, &
     dgdot_dtau_slip_neg
- 
+
   real(pReal), dimension(param(instance)%totalNslip) :: &
     tau_slip_pos, &
     tau_slip_neg
   integer :: i
   logical :: nonSchmidActive
- 
+
   associate(prm => param(instance), stt => state(instance))
- 
+
   nonSchmidActive = size(prm%nonSchmidCoeff) > 0
- 
+
   do i = 1, prm%totalNslip
     tau_slip_pos(i) =       math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i))
     tau_slip_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)), &
                             0.0_pReal, nonSchmidActive)
   enddo
- 
+
   where(dNeq0(tau_slip_pos))
     gdot_slip_pos = prm%gdot0_slip * merge(0.5_pReal,1.0_pReal, nonSchmidActive) &                  ! 1/2 if non-Schmid active
                   * sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_slip,  tau_slip_pos)
   else where
     gdot_slip_pos = 0.0_pReal
   end where
- 
+
   where(dNeq0(tau_slip_neg))
     gdot_slip_neg = prm%gdot0_slip * 0.5_pReal &                                                    ! only used if non-Schmid active, always 1/2
                   * sign(abs(tau_slip_neg/stt%xi_slip(:,of))**prm%n_slip,  tau_slip_neg)
   else where
     gdot_slip_neg = 0.0_pReal
   end where
- 
+
   if (present(dgdot_dtau_slip_pos)) then
     where(dNeq0(gdot_slip_pos))
       dgdot_dtau_slip_pos = gdot_slip_pos*prm%n_slip/tau_slip_pos
@@ -543,9 +488,9 @@ end subroutine kinetics_slip
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Shear rates on twin systems and their derivatives with respect to resolved stress.
+!> @brief Calculate shear rates on twin systems and their derivatives with respect to resolved
-!  twinning is assumed to take place only in untwinned volume.
+!         stress. Twinning is assumed to take place only in untwinned volume.
-!> @details Derivates are calculated only optionally.
+!> @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.
 !--------------------------------------------------------------------------------------------------
@@ -557,29 +502,29 @@ pure subroutine kinetics_twin(Mp,instance,of,&
   integer,                      intent(in) :: &
     instance, &
     of
- 
+
   real(pReal), dimension(param(instance)%totalNtwin), intent(out) :: &
     gdot_twin
   real(pReal), dimension(param(instance)%totalNtwin), intent(out), optional :: &
     dgdot_dtau_twin
- 
+
   real(pReal), dimension(param(instance)%totalNtwin) :: &
     tau_twin
   integer :: i
- 
+
   associate(prm => param(instance), stt => state(instance))
- 
+
   do i = 1, prm%totalNtwin
     tau_twin(i)  = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i))
   enddo
- 
+
   where(tau_twin > 0.0_pReal)
     gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)) &                         ! only twin in untwinned volume fraction
               * prm%gdot0_twin*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_twin
   else where
     gdot_twin = 0.0_pReal
   end where
- 
+
   if (present(dgdot_dtau_twin)) then
     where(dNeq0(gdot_twin))
       dgdot_dtau_twin = gdot_twin*prm%n_twin/tau_twin
@@ -587,7 +532,7 @@ pure subroutine kinetics_twin(Mp,instance,of,&
       dgdot_dtau_twin = 0.0_pReal
     end where
   endif
- 
+
   end associate
 
 end subroutine kinetics_twin

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_subStep(c,i,e) *( crystallite_partionedF (1:3,1:3,c,i,e) &
-                                                                          - crystallite_partionedF0(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_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)), &
+                                                          math_inv33(crystallite_Fp(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_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), &
+        temp_33_1 = transpose(matmul(invFp,invFi))
-                                     crystallite_invFi(1:3,1:3,c,i,e)))
+        temp_33_2 = matmul(crystallite_subF(1:3,1:3,c,i,e),invSubFp0)
-        temp_33_2 = matmul(crystallite_subF(1:3,1:3,c,i,e), &
+        temp_33_3 = matmul(matmul(crystallite_subF(1:3,1:3,c,i,e),invFp), invSubFi0)
-                           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)))
 
         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)), &
+          temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), invFi) &
-                                          crystallite_invFi(1:3,1:3,c,i,e)) &
                                  + 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(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi))
-                                         matmul(temp_3333(1:3,1:3,p,o),crystallite_invFi(1:3,1:3,c,i,e)))
         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_1 = matmul(crystallite_S(1:3,1:3,c,i,e),transpose(invFp))
-        temp_33_2 = matmul(crystallite_invFp(1:3,1:3,c,i,e),temp_33_1)
+        temp_33_2 = matmul(invFp,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_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
 
@@ -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/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   :: &
+    character(len=pStringLen), allocatable, dimension(:) :: &
-      outputID
+      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  -+>>>'
+  write(6,'(/,a)') ' <<<+-  homogenization_'//HOMOGENIZATION_RGC_label//' init  -+>>>'; flush(6)
  
-  write(6,'(/,a)')   ' Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009'
+  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)')   ' 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)') ' 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)')   ' 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,34 +107,8 @@ 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)
+    prm%output = 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
- 
     NofMyHomog = count(material_homogenizationAt == h)
     nIntFaceTot = 3*(  (prm%Nconstituents(1)-1)*prm%Nconstituents(2)*prm%Nconstituents(3) &
                       + prm%Nconstituents(1)*(prm%Nconstituents(2)-1)*prm%Nconstituents(3) &
@@ -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%output)
-  outputsLoop: do o = 1,size(prm%outputID)
+    select case(trim(prm%output(o)))
-    select case(prm%outputID(o))
+      case('constitutivework')
-    
-      case (constitutivework_ID)
         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/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 :: &
+    character(len=pStringLen), allocatable, dimension(:) :: &
-      outputID
+      output
   end type tParameters
   
   type(tparameters),             dimension(:), allocatable :: &
@@ -46,10 +41,9 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine thermal_adiabatic_init
  
-  integer :: maxNinstance,o,h,NofMyHomog   
+  integer :: maxNinstance,h,NofMyHomog   
-  character(len=pStringLen), dimension(:), allocatable :: outputs
  
-  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)
+    prm%output = 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
 
     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%output)
-  outputsLoop: do o = 1,size(prm%outputID)
+    select case(trim(prm%output(o)))
-    select case(prm%outputID(o))
+      case('temperature')                                                                           ! ToDo: should be 'T'
-    
-      case (temperature_ID)
         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 :: &
+    character(len=pStringLen), allocatable, dimension(:) :: &
-      outputID
+      output
   end type tParameters
   
   type(tparameters),             dimension(:), allocatable :: &
@@ -48,10 +42,9 @@ contains
 subroutine thermal_conduction_init
 
   
-  integer :: maxNinstance,o,NofMyHomog,h
+  integer :: maxNinstance,NofMyHomog,h
-  character(len=pStringLen), dimension(:), allocatable :: outputs
  
-  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)
+    prm%output = 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
   
     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%output)
-  outputsLoop: do o = 1,size(prm%outputID)
+    select case(trim(prm%output(o)))
-    select case(prm%outputID(o))
+      case('temperature')                                                                           ! ToDo: should be 'T'
-    
-      case (temperature_ID)
         call results_writeDataset(group,temperature(homog)%p,'T',&
                                   'temperature','K')
     end select