Merge branch 'development' into 32_NewStyleNonlocal-3

.gitlab-ci.yml

@@ -482,24 +482,40 @@ AbaqusStd:
   script:
     - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel $Doxygen
     - $DAMASKROOT/PRIVATE/documenting/runDoxygen.sh $DAMASKROOT abaqus
-  only:
+  except:
-    - development
+    - master
+    - release
 
 Marc:
   stage: createDocumentation
   script:
     - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel $Doxygen
     - $DAMASKROOT/PRIVATE/documenting/runDoxygen.sh $DAMASKROOT marc
-  only:
+  except:
-    - development
+    - master
+    - release
 
 Spectral:
   stage: createDocumentation
   script:
     - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel $Doxygen
     - $DAMASKROOT/PRIVATE/documenting/runDoxygen.sh $DAMASKROOT spectral
-  only:
+  except:
-    - development
+    - master
+    - release
+
+Processing:
+  stage: createDocumentation
+  script:
+    - cd $DAMASKROOT/processing/pre
+    - rm 3DRVEfrom2Dang.py abq_addUserOutput.py marc_addUserOutput.py
+    - $DAMASKROOT/PRIVATE/documenting/scriptHelpToWiki.py --debug *.py
+    - cd $DAMASKROOT/processing/post
+    - rm marc_to_vtk.py vtk2ang.py
+    - $DAMASKROOT/PRIVATE/documenting/scriptHelpToWiki.py --debug *.py
+  except:
+    - master
+    - release
 
 ##################################################################################################
 backupData:
@@ -508,11 +524,10 @@ backupData:
     - cd $TESTROOT/performance                                                                      # location of new runtime results
     - git commit -am"${CI_PIPELINE_ID}_${CI_COMMIT_SHA}"
     - mkdir $BACKUP/${CI_PIPELINE_ID}_${CI_COMMIT_SHA}
-    - cp $TESTROOT/performance/time.txt $BACKUP/${CI_PIPELINE_ID}_${CI_COMMIT_SHA}/
     - mv $TESTROOT/performance/time.png $BACKUP/${CI_PIPELINE_ID}_${CI_COMMIT_SHA}/
-    - cp $TESTROOT/performance/memory.txt $BACKUP/${CI_PIPELINE_ID}_${CI_COMMIT_SHA}/
     - mv $TESTROOT/performance/memory.png $BACKUP/${CI_PIPELINE_ID}_${CI_COMMIT_SHA}/
     - mv $DAMASKROOT/PRIVATE/documenting/DAMASK_* $BACKUP/${CI_PIPELINE_ID}_${CI_COMMIT_SHA}/
+    - mv $DAMASKROOT/processing $BACKUP/${CI_PIPELINE_ID}_${CI_COMMIT_SHA}/
   only:
     - development
 

PRIVATE

@@ -1 +1 @@
-Subproject commit 18ba1ba6a5e9ba446dc9311acf2acf2781614db1
+Subproject commit ddb0dae72af9012cca45d9fa5665da41815e88f7

VERSION

@@ -1 +1 @@
-v2.0.2-1789-g524bfb8c
+v2.0.2-1826-gd2a9f55a

processing/post/addAPS34IDEstrainCoords.py

@@ -14,20 +14,15 @@ scriptID   = ' '.join([scriptName,damask.version])
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Transform X,Y,Z,F APS BeamLine 34 coordinates to x,y,z APS strain coordinates.
 
 """, version = scriptID)
 
-parser.add_option('-f',
+parser.add_option('-f','--frame',dest='frame', metavar='string',
-                  '--frame',
+                                 help='label of APS X,Y,Z coords')
-                  dest='frame',
+parser.add_option('--depth',     dest='depth',          metavar='string',
-                  metavar='string',
+                                 help='depth')
-                  help='APS X,Y,Z coords')
-parser.add_option('--depth',
-                  dest='depth',
-                  metavar='string',
-                  help='depth')
 
 (options,filenames) = parser.parse_args()
 

processing/post/addCalculation.py

@@ -18,7 +18,7 @@ def listify(x):
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Add or alter column(s) with derived values according to user-defined arithmetic operation between column(s).
 Column labels are tagged by '#label#' in formulas. Use ';' for ',' in functions.
 Numpy is available as 'np'.
@@ -41,10 +41,7 @@ parser.add_option('-f','--formula',
 
 parser.add_option('-c','--condition',
                   dest   = 'condition', metavar='string',
-                  help   = 'condition to alter existing column data')
+                  help   = 'condition to alter existing column data (optional)')
-
-parser.set_defaults(condition = None,
-                   )
 
 (options,filenames) = parser.parse_args()
 
@@ -80,7 +77,7 @@ for name in filenames:
     condition = options.condition                                                                   # copy per file, since might be altered inline
     breaker = False
   
-    for position,(all,marker,column) in enumerate(set(re.findall(r'#(([s]#)?(.+?))#',condition))):              # find three groups
+    for position,(all,marker,column) in enumerate(set(re.findall(r'#(([s]#)?(.+?))#',condition))):  # find three groups
       idx = table.label_index(column)
       dim = table.label_dimension(column)
       if idx < 0 and column not in specials:
@@ -89,15 +86,15 @@ for name in filenames:
       else:
         if column in specials:
           replacement = 'specials["{}"]'.format(column)
-        elif dim == 1:                                                                                # scalar input
+        elif dim == 1:                                                                              # scalar input
           replacement = '{}(table.data[{}])'.format({  '':'float',
-                                                        's#':'str'}[marker],idx)                      # take float or string value of data column
+                                                        's#':'str'}[marker],idx)                    # take float or string value of data column
-        elif dim > 1:                                                                                 # multidimensional input (vector, tensor, etc.)
+        elif dim > 1:                                                                               # multidimensional input (vector, tensor, etc.)
-          replacement = 'np.array(table.data[{}:{}],dtype=float)'.format(idx,idx+dim)                 # use (flat) array representation
+          replacement = 'np.array(table.data[{}:{}],dtype=float)'.format(idx,idx+dim)               # use (flat) array representation
        
         condition = condition.replace('#'+all+'#',replacement)
     
-    if breaker: continue                                                                              # found mistake in condition evaluation --> next file
+    if breaker: continue                                                                            # found mistake in condition evaluation --> next file
 
 # ------------------------------------------ build formulas ----------------------------------------
 

processing/post/addCauchy.py

@@ -13,8 +13,8 @@ scriptID   = ' '.join([scriptName,damask.version])
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
-Add column(s) containing Cauchy stress based on given column(s) of deformation gradient and first Piola--Kirchhoff stress.
+Add column containing Cauchy stress based on deformation gradient and first Piola--Kirchhoff stress.
 
 """, version = scriptID)
 

processing/post/addCompatibilityMismatch.py

@@ -209,7 +209,7 @@ def shapeMismatch(size,F,nodes,centres):
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options file[s]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Add column(s) containing the shape and volume mismatch resulting from given deformation gradient.
 Operates on periodic three-dimensional x,y,z-ordered data sets.
 

processing/post/addCumulative.py

@@ -13,7 +13,7 @@ scriptID   = ' '.join([scriptName,damask.version])
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Add cumulative (sum of first to current row) values for given label(s).
 """, version = scriptID)
 
@@ -22,12 +22,9 @@ parser.add_option('-l','--label',
                   action = 'extend', metavar = '<string LIST>',
                   help = 'columns to cumulate')
 
-parser.set_defaults(label = [],
-                   )
-                    
 (options,filenames) = parser.parse_args()
 
-if len(options.label) == 0:
+if options.label is None:
   parser.error('no data column(s) specified.')
 
 # --- loop over input files -------------------------------------------------------------------------

processing/post/addCurl.py

@@ -49,14 +49,14 @@ def curlFFT(geomdim,field):
 
   curl_fourier = np.einsum(einsums[n],e,k_s,field_fourier)*TWOPIIMG
 
-  return np.fft.irfftn(curl_fourier,s=shapeFFT,axes=(0,1,2)).reshape([N,n])
+  return np.fft.irfftn(curl_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n])
 
 
 # --------------------------------------------------------------------
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [ASCIItable(s)]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Add column(s) containing curl of requested column(s).
 Operates on periodic ordered three-dimensional data sets
 of vector and tensor fields.

processing/post/addDerivative.py

@@ -34,7 +34,7 @@ def derivative(coordinates,what):
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Add column(s) containing numerical derivative of requested column(s) with respect to given coordinates.
 
 """, version = scriptID)

processing/post/addDeterminant.py

@@ -20,7 +20,7 @@ def determinant(m):
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Add column(s) containing determinant of requested tensor column(s).
 
 """, version = scriptID)

processing/post/addDeviator.py

@@ -23,7 +23,7 @@ def deviator(m,spherical = False):
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(2)]', description = """
 Add column(s) containing deviator of requested tensor column(s).
 
 """, version = scriptID)

processing/post/addDisplacement.py

@@ -87,7 +87,7 @@ def displacementFluctFFT(F,grid,size,nodal=False,transformed=False):
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog [options] [ASCIItable(s)]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Add displacments resulting from deformation gradient field.
 Operates on periodic three-dimensional x,y,z-ordered data sets.
 Outputs at cell centers or cell nodes (into separate file).
@@ -111,7 +111,6 @@ parser.add_option('--nodal',
 
 parser.set_defaults(defgrad = 'f',
                     pos     = 'pos',
-                    nodal   = False,
                    )
 
 (options,filenames) = parser.parse_args()

processing/post/addDivergence.py

@@ -45,7 +45,7 @@ def divFFT(geomdim,field):
 
   div_fourier = np.einsum(einsums[n],k_s,field_fourier)*TWOPIIMG
 
-  return np.fft.irfftn(div_fourier,s=shapeFFT,axes=(0,1,2)).reshape([N,n//3])
+  return np.fft.irfftn(div_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,n//3])
 
 
 # --------------------------------------------------------------------

processing/post/addEhkl.py

@@ -30,7 +30,7 @@ def E_hkl(stiffness,vec):   # stiffness = (c11,c12,c44)
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Add column(s) containing directional stiffness based on given cubic stiffness values C11, C12, and C44 in consecutive columns.
 
 """, version = scriptID)

processing/post/addEuclideanDistance.py

@@ -83,7 +83,7 @@ neighborhoods = {
                                       ])
                 }
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Add column(s) containing Euclidean distance to grain structural features: boundaries, triple lines, and quadruple points.
 
 """, version = scriptID)

processing/post/addGaussian.py

@@ -15,7 +15,7 @@ scriptID   = ' '.join([scriptName,damask.version])
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [ASCIItable(s)]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog option [ASCIItable(s)]', description = """
 Add column(s) containing Gaussian filtered values of requested column(s).
 Operates on periodic and non-periodic ordered three-dimensional data sets.
 For details see scipy.ndimage documentation.
@@ -34,12 +34,12 @@ parser.add_option('-o','--order',
                   dest = 'order',
                   type = int,
                   metavar = 'int',
-                  help = 'order of the filter')
+                  help = 'order of the filter [%default]')
 parser.add_option('--sigma',
                   dest = 'sigma',
                   type = float,
                   metavar = 'float',
-                  help = 'standard deviation')
+                  help = 'standard deviation [%default]')
 parser.add_option('--periodic',
                   dest = 'periodic',
                   action = 'store_true',
@@ -50,7 +50,6 @@ parser.add_option('--periodic',
 parser.set_defaults(pos = 'pos',
                     order = 0,
                     sigma = 1,
-                    periodic = False,
                    )
 
 (options,filenames) = parser.parse_args()

processing/post/addGradient.py

@@ -45,14 +45,14 @@ def gradFFT(geomdim,field):
   k_s = np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3).astype('c16')                           
   grad_fourier = np.einsum(einsums[n],field_fourier,k_s)*TWOPIIMG
 
-  return np.fft.irfftn(grad_fourier,s=shapeFFT,axes=(0,1,2)).reshape([N,3*n])
+  return np.fft.irfftn(grad_fourier,axes=(0,1,2),s=shapeFFT).reshape([N,3*n])
 
 
 # --------------------------------------------------------------------
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [ASCIItable(s)]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog option [ASCIItable(s)]', description = """
 Add column(s) containing gradient of requested column(s).
 Operates on periodic ordered three-dimensional data sets
 of vector and scalar fields.

processing/post/addGrainID.py

@@ -28,9 +28,9 @@ parser.add_option('-d',
                   help = 'disorientation threshold in degrees [%default]')
 parser.add_option('-s',
                   '--symmetry',
-                  dest = 'symmetry',
+                  dest = 'symmetry', type = 'choice', choices = damask.Symmetry.lattices[1:],
                   metavar = 'string',
-                  help = 'crystal symmetry [%default]')
+                  help = 'crystal symmetry [%default] {{{}}} '.format(', '.join(damask.Symmetry.lattices[1:])))
 parser.add_option('-o',
                   '--orientation',
                   dest = 'quaternion',
@@ -49,7 +49,7 @@ parser.add_option('--quiet',
 parser.set_defaults(disorientation = 5,
                     verbose = True,
                     quaternion = 'orientation',
-                    symmetry = 'cubic',
+                    symmetry = damask.Symmetry.lattices[-1],
                     pos      = 'pos',
                    )
 

processing/post/addIPFcolor.py

@@ -13,7 +13,7 @@ scriptID   = ' '.join([scriptName,damask.version])
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Add RGB color value corresponding to TSL-OIM scheme for inverse pole figures.
 
 """, version = scriptID)

processing/post/addIndexed.py

@@ -13,7 +13,7 @@ scriptID   = ' '.join([scriptName,damask.version])
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Add data in column(s) of mapped ASCIItable selected from the row indexed by the value in a mapping column.
 Row numbers start at 1.
 

processing/post/addInfo.py

@@ -13,7 +13,7 @@ scriptID   = ' '.join([scriptName,damask.version])
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options file[s]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Add info lines to ASCIItable header.
 
 """, version = scriptID)
@@ -23,11 +23,12 @@ parser.add_option('-i',
                   dest = 'info', action = 'extend', metavar = '<string LIST>',
                   help    = 'items to add')
 
-parser.set_defaults(info = [],
-                   )
 
 (options,filenames) = parser.parse_args()
 
+if options.info is None:
+  parser.error('no info specified.')
+
 # --- loop over input files ------------------------------------------------------------------------
 
 if filenames == []: filenames = [None]

processing/post/addLinked.py

@@ -13,7 +13,7 @@ scriptID   = ' '.join([scriptName,damask.version])
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Add data of selected column(s) from (first) row of linked ASCIItable that shares the linking column value.
 
 """, version = scriptID)

processing/post/addMises.py

@@ -24,7 +24,7 @@ def Mises(what,tensor):
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Add vonMises equivalent values for symmetric part of requested strains and/or stresses.
 
 """, version = scriptID)
@@ -38,13 +38,9 @@ parser.add_option('-s','--stress',
                   action = 'extend', metavar = '<string LIST>',
                   help = 'heading(s) of columns containing stress tensors')
 
-parser.set_defaults(strain = [],
-                    stress = [],
-                   )
-
 (options,filenames) = parser.parse_args()
 
-if len(options.stress+options.strain) == 0:
+if options.stress is None and options.strain is None:
   parser.error('no data column specified...')
 
 # --- loop over input files -------------------------------------------------------------------------

processing/post/addNorm.py

@@ -9,6 +9,7 @@ scriptName = os.path.splitext(os.path.basename(__file__))[0]
 scriptID   = ' '.join([scriptName,damask.version])
 
 # definition of element-wise p-norms for matrices
+# ToDo: better use numpy.linalg.norm
 
 def norm(which,object):
 
@@ -18,12 +19,14 @@ def norm(which,object):
     return math.sqrt(sum([x*x for x in object]))
   elif which == 'Max':                                                                              # p = inf
     return max(map(abs, object))
+  else:
+    return -1
 
 # --------------------------------------------------------------------
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Add column(s) containing norm of requested column(s) being either vectors or tensors.
 
 """, version = scriptID)
@@ -43,6 +46,8 @@ parser.set_defaults(norm = 'frobenius',
 
 (options,filenames) = parser.parse_args()
 
+if options.norm.lower() not in normChoices:
+  parser.error('invalid norm ({}) specified.'.format(options.norm))
 if options.label is None:
   parser.error('no data column specified.')
 
@@ -74,7 +79,7 @@ for name in filenames:
     else:
       dims.append(dim)
       columns.append(table.label_index(what))
-      table.labels_append('norm{}({})'.format(options.norm.capitalize(),what))                    # extend ASCII header with new labels
+      table.labels_append('norm{}({})'.format(options.norm.capitalize(),what))                      # extend ASCII header with new labels
 
   if remarks != []: damask.util.croak(remarks)
   if errors  != []:

processing/post/addOrientations.py

@@ -38,7 +38,7 @@ def check_matrix(M):
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Add quaternion and/or Bunge Euler angle representation of crystal lattice orientation.
 Orientation is given by quaternion, Euler angles, rotation matrix, or crystal frame coordinates
 (i.e. component vectors of rotation matrix).
@@ -68,12 +68,12 @@ parser.add_option('-R',
                   '--labrotation',
                   dest='labrotation',
                   type = 'float', nargs = 4, metavar = ' '.join(['float']*4),
-                  help = 'angle and axis of additional lab frame rotation')
+                  help = 'angle and axis of additional lab frame rotation [%default]')
 parser.add_option('-r',
                   '--crystalrotation',
                   dest='crystalrotation',
                   type = 'float', nargs = 4, metavar = ' '.join(['float']*4),
-                  help = 'angle and axis of additional crystal frame rotation')
+                  help = 'angle and axis of additional crystal frame rotation [%default]')
 parser.add_option('--eulers',
                   dest = 'eulers',
                   metavar = 'string',
@@ -106,7 +106,6 @@ parser.add_option('-z',
 parser.set_defaults(output = [],
                     labrotation     = (0.,1.,1.,1.),                                                # no rotation about 1,1,1
                     crystalrotation = (0.,1.,1.,1.),                                                # no rotation about 1,1,1
-                    degrees = False,
                    )
 
 (options, filenames) = parser.parse_args()

processing/post/addPK2.py

@@ -13,7 +13,7 @@ scriptID   = ' '.join([scriptName,damask.version])
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Add column(s) containing Second Piola--Kirchhoff stress based on given column(s) of deformation
 gradient and first Piola--Kirchhoff stress.
 

processing/post/addPole.py

@@ -13,7 +13,7 @@ scriptID   = ' '.join([scriptName,damask.version])
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Add coordinates of stereographic projection of given direction (pole) in crystal frame.
 
 """, version = scriptID)
@@ -35,7 +35,6 @@ parser.add_option('-o',
 
 parser.set_defaults(pole = (1.0,0.0,0.0),
                     quaternion = 'orientation',
-                    polar   = False,
                    )
 
 (options, filenames) = parser.parse_args()

processing/post/addSchmidfactors.py

@@ -103,7 +103,7 @@ slipSystems = {
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Add columns listing Schmid factors (and optional trace vector of selected system) for given Euler angles.
 
 """, version = scriptID)
@@ -115,7 +115,7 @@ parser.add_option('-l',
                   help = 'type of lattice structure [%default] {}'.format(latticeChoices))
 parser.add_option('--covera',
                   dest = 'CoverA', type = 'float', metavar = 'float',
-                  help = 'C over A ratio for hexagonal systems')
+                  help = 'C over A ratio for hexagonal systems [%default]')
 parser.add_option('-f',
                   '--force',
                   dest = 'force',

processing/post/addSpectralDecomposition.py

@@ -13,7 +13,7 @@ scriptID   = ' '.join([scriptName,damask.version])
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Add column(s) containing eigenvalues and eigenvectors of requested symmetric tensor column(s).
 
 """, version = scriptID)

processing/post/addStrainTensors.py

@@ -25,7 +25,7 @@ def operator(stretch,strain,eigenvalues):
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Add column(s) containing given strains based on given stretches of requested deformation gradient column(s).
 
 """, version = scriptID)
@@ -56,16 +56,15 @@ parser.add_option('-f','--defgrad',
                   metavar = '<string LIST>',
                   help = 'heading(s) of columns containing deformation tensor values [%default]')
 
-parser.set_defaults(right       = False,
+parser.set_defaults(
-                    left        = False,
-                    logarithmic = False,
-                    biot        = False,
-                    green       = False,
                     defgrad     = ['f'],
                    )
 
 (options,filenames) = parser.parse_args()
 
+if len(options.defgrad) > 1:
+  options.defgrad = options.defgrad[1:]
+
 stretches = []
 strains = []
 

processing/post/addTable.py

@@ -12,7 +12,7 @@ scriptID   = ' '.join([scriptName,damask.version])
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Append data of ASCIItable(s).
 
 """, version = scriptID)
@@ -24,6 +24,10 @@ parser.add_option('-a', '--add','--table',
 
 (options,filenames) = parser.parse_args()
 
+if options.table is None:
+  parser.error('no table specified.')
+
+
 # --- loop over input files -------------------------------------------------------------------------
 
 if filenames == []: filenames = [None]

processing/post/averageDown.py

@@ -14,7 +14,7 @@ scriptID   = ' '.join([scriptName,damask.version])
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Average each data block of size 'packing' into single values thus reducing the former grid to grid/packing.
 
 """, version = scriptID)
@@ -34,16 +34,14 @@ parser.add_option('--shift',
 parser.add_option('-g', '--grid',
                   dest = 'grid',
                   type = 'int', nargs = 3, metavar = 'int int int',
-                  help = 'grid in x,y,z [autodetect]')
+                  help = 'grid in x,y,z (optional)')
 parser.add_option('-s', '--size',
                   dest = 'size',
                   type = 'float', nargs = 3, metavar = 'float float float',
-                  help = 'size in x,y,z [autodetect]')
+                  help = 'size in x,y,z (optional)')
 parser.set_defaults(pos     = 'pos',
                     packing = (2,2,2),
                     shift   = (0,0,0),
-                    grid    = (0,0,0),
-                    size    = (0.0,0.0,0.0),
                    )
 
 (options,filenames) = parser.parse_args()
@@ -92,7 +90,7 @@ for name in filenames:
 
   table.data_readArray()
 
-  if (any(options.grid) == 0 or any(options.size) == 0.0):
+  if (options.grid is None or options.size is None):
     grid,size = damask.util.coordGridAndSize(table.data[:,table.label_indexrange(options.pos)])
   else:
     grid   = np.array(options.grid,'i')

processing/post/binXY.py

@@ -13,7 +13,7 @@ scriptID   = ' '.join([scriptName,damask.version])
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Produces a binned grid of two columns from an ASCIItable, i.e. a two-dimensional probability density map.
 
 """, version = scriptID)
@@ -37,15 +37,15 @@ parser.add_option('-t','--type',
 parser.add_option('-x','--xrange',
                   dest = 'xrange',
                   type = 'float', nargs = 2, metavar = 'float float',
-                  help = 'min max value in x direction [autodetect]')
+                  help = 'min max limits in x direction (optional)')
 parser.add_option('-y','--yrange',
                   dest = 'yrange',
                   type = 'float', nargs = 2, metavar = 'float float',
-                  help = 'min max value in y direction [autodetect]')
+                  help = 'min max limits in y direction (optional)')
 parser.add_option('-z','--zrange',
                   dest = 'zrange',
                   type = 'float', nargs = 2, metavar = 'float float',
-                  help = 'min max value in z direction [autodetect]')
+                  help = 'min max limits in z direction (optional)')
 parser.add_option('-i','--invert',
                   dest = 'invert',
                   action = 'store_true',
@@ -64,9 +64,6 @@ parser.set_defaults(bins = (10,10),
                     xrange = (0.0,0.0),
                     yrange = (0.0,0.0),
                     zrange = (0.0,0.0),
-                    invert = False,
-                    normRow = False,
-                    normCol = False,
                    )
 
 (options,filenames) = parser.parse_args()

processing/post/blowUp.py

@@ -13,7 +13,7 @@ scriptID   = ' '.join([scriptName,damask.version])
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Blows up each value to a surrounding data block of size 'packing' thus increasing the former resolution
 to resolution*packing.
 
@@ -27,10 +27,10 @@ parser.add_option('-p','--packing',
                   help = 'dimension of packed group [%default]')
 parser.add_option('-g','--grid',
                   dest = 'resolution', type = 'int', nargs = 3, metavar = 'int int int',
-                  help = 'resolution in x,y,z [autodetect]')
+                  help = 'grid in x,y,z (optional)')
 parser.add_option('-s','--size',
                   dest = 'dimension', type = 'float', nargs = 3, metavar = 'int int int',
-                  help = 'dimension in x,y,z [autodetect]')
+                  help = 'size in x,y,z (optional)')
 parser.set_defaults(pos  = 'pos',
                     packing = (2,2,2),
                     grid    = (0,0,0),

processing/post/filterTable.py

@@ -30,7 +30,7 @@ def sortingList(labels,whitelistitems):
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Filter rows according to condition and columns by either white or black listing.
 
 Examples:

processing/post/groupTable.py

@@ -20,7 +20,7 @@ scriptID   = ' '.join([scriptName,damask.version])
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Apply a user-specified function to condense into a single row all those rows for which columns 'label' have identical values.
 Output table will contain as many rows as there are different (unique) values in the grouping column(s).
 Periodic domain averaging of coordinate values is supported.

processing/post/permuteData.py

@@ -13,7 +13,7 @@ scriptID   = ' '.join([scriptName,damask.version])
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Permute all values in given column(s).
 
 """, version = scriptID)

processing/post/reLabel.py

@@ -12,7 +12,7 @@ scriptID   = ' '.join([scriptName,damask.version])
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog [options] dfile[s]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Rename scalar, vectorial, and/or tensorial data header labels.
 
 """, version = scriptID)

processing/post/rotateData.py

@@ -13,7 +13,7 @@ scriptID   = ' '.join([scriptName,damask.version])
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Rotate vector and/or tensor column data by given angle around given axis.
 
 """, version = scriptID)
@@ -29,7 +29,7 @@ parser.add_option('-r', '--rotation',
 parser.add_option('--degrees',
                   dest = 'degrees',
                   action = 'store_true',
-                  help = 'angles are given in degrees [%default]')
+                  help = 'angles are given in degrees')
 
 parser.set_defaults(rotation = (0.,1.,1.,1.),                                                       # no rotation about 1,1,1
                     degrees = False,

processing/post/scaleData.py

@@ -13,7 +13,7 @@ scriptID   = ' '.join([scriptName,damask.version])
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Uniformly scale column values by given factor.
 
 """, version = scriptID)

processing/post/shiftData.py

@@ -13,7 +13,7 @@ scriptID   = ' '.join([scriptName,damask.version])
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Uniformly shift column values by given offset.
 
 """, version = scriptID)

processing/post/sortTable.py

@@ -13,7 +13,7 @@ scriptID   = ' '.join([scriptName,damask.version])
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Sort rows by given (or all) column label(s).
 
 Examples:

processing/post/viewTable.py

@@ -12,7 +12,7 @@ scriptID   = ' '.join([scriptName,damask.version])
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(usage='%prog [options] [file[s]]', description = """
+parser = OptionParser(usage='%prog options [ASCIItable(s)]', description = """
 Show components of given ASCIItable(s).
 
 """, version = scriptID)

processing/post/vtk_addGridData.py

@@ -17,7 +17,7 @@ scriptID   = ' '.join([scriptName,damask.version])
 msg = "Add scalars, vectors, and/or an RGB tuple from"
 msg += "an ASCIItable to existing VTK grid (.vtr/.vtk/.vtu)."
 parser = OptionParser(option_class=damask.extendableOption,
-                      usage='%prog options [file[s]]',
+                      usage='%prog options [ASCIItable(s)]',
                       description = msg,
                       version = scriptID)
 
@@ -172,8 +172,7 @@ for name in filenames:
 
   writer.SetDataModeToBinary()
   writer.SetCompressorTypeToZLib()
-  if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(rGrid)
+  writer.SetInputData(rGrid)
-  else:                          writer.SetInputData(rGrid)
   writer.Write()
 
 # ------------------------------------------ render result ---------------------------------------

processing/post/vtk_addPointcloudData.py

@@ -15,7 +15,7 @@ scriptID   = ' '.join([scriptName,damask.version])
 # --------------------------------------------------------------------
 
 parser = OptionParser(option_class=damask.extendableOption,
-                      usage='%prog options [file[s]]',
+                      usage='%prog options [ASCIItable(s)]',
                       description = """Add scalar and RGB tuples from ASCIItable to existing VTK point cloud (.vtp).""",
                       version = scriptID)
 
@@ -46,8 +46,6 @@ parser.add_option('-c', '--color',   dest='color', action='extend',
 parser.set_defaults(data = [],
                     tensor = [],
                     color = [],
-                    inplace = False,
-                    render = False,
 )
 
 (options, filenames) = parser.parse_args()
@@ -151,14 +149,12 @@ for name in filenames:
 # ------------------------------------------ output result ---------------------------------------
 
   Polydata.Modified()
-  if vtk.VTK_MAJOR_VERSION <= 5: Polydata.Update()
 
   writer = vtk.vtkXMLPolyDataWriter()
   writer.SetDataModeToBinary()
   writer.SetCompressorTypeToZLib()
   writer.SetFileName(os.path.splitext(options.vtk)[0]+('.vtp' if options.inplace else '_added.vtp'))
-  if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(Polydata)
+  writer.SetInputData(Polydata)
-  else:                          writer.SetInputData(Polydata)
   writer.Write()
 
 # ------------------------------------------ render result ---------------------------------------

processing/post/vtk_pointCloud.py

@@ -13,7 +13,7 @@ scriptID   = ' '.join([scriptName,damask.version])
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Produce a VTK point cloud dataset based on coordinates given in an ASCIItable.
 
 """, version = scriptID)
@@ -78,7 +78,6 @@ for name in filenames:
   Polydata.SetPoints(Points)
   Polydata.SetVerts(Vertices)
   Polydata.Modified()
-  if vtk.VTK_MAJOR_VERSION <= 5: Polydata.Update()
  
 # ------------------------------------------ output result ---------------------------------------  
 
@@ -94,8 +93,8 @@ for name in filenames:
     writer.SetHeader('# powered by '+scriptID)
     writer.WriteToOutputStringOn()
   
-  if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(Polydata)
+  
-  else:                          writer.SetInputData(Polydata)
+  writer.SetInputData(Polydata)
 
   writer.Write()
 

processing/post/vtk_rectilinearGrid.py

@@ -13,7 +13,7 @@ scriptID   = ' '.join([scriptName,damask.version])
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
 Create regular voxel grid from points in an ASCIItable.
 
 """, version = scriptID)
@@ -125,8 +125,7 @@ for name in filenames:
     writer.SetHeader('# powered by '+scriptID)
     writer.WriteToOutputStringOn()
   
-  if vtk.VTK_MAJOR_VERSION <= 5: writer.SetInput(rGrid)
+  writer.SetInputData(rGrid)
-  else:                          writer.SetInputData(rGrid)
 
   writer.Write()
 

processing/pre/geom_addPrimitive.py

@@ -25,7 +25,7 @@ mappings = {
         'microstructures': lambda x: int(x),
           }
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [geomfile(s)]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog option [geomfile(s)]', description = """
 Positions a geometric object within the (three-dimensional) canvas of a spectral geometry description.
 Depending on the sign of the dimension parameters, these objects can be boxes, cylinders, or ellipsoids.
 

processing/pre/geom_clean.py

@@ -18,7 +18,7 @@ def mostFrequent(arr):
 #                                MAIN
 #--------------------------------------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [geomfile(s)]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [geomfile(s)]', description = """
 Smooth geometry by selecting most frequent microstructure index within given stencil at each location.
 
 """, version=scriptID)

processing/pre/seeds_fromRandom.py

@@ -28,7 +28,7 @@ def kdtree_search(cloud, queryPoints):
 #                                MAIN
 # --------------------------------------------------------------------
 
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog [options]', description = """
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog options', description = """
 Distribute given number of points randomly within (a fraction of) the three-dimensional cube [0.0,0.0,0.0]--[1.0,1.0,1.0].
 Reports positions with random crystal orientations in seeds file format to STDOUT.
 

src/constitutive.f90

@@ -160,16 +160,16 @@ subroutine constitutive_init()
  call IO_checkAndRewind(FILEUNIT)
  if (any(phase_source == SOURCE_thermal_dissipation_ID))     call source_thermal_dissipation_init(FILEUNIT)
  if (any(phase_source == SOURCE_thermal_externalheat_ID))    call source_thermal_externalheat_init(FILEUNIT)
- if (any(phase_source == SOURCE_damage_isoBrittle_ID))       call source_damage_isoBrittle_init(FILEUNIT)
+ if (any(phase_source == SOURCE_damage_isoBrittle_ID))       call source_damage_isoBrittle_init
- if (any(phase_source == SOURCE_damage_isoDuctile_ID))       call source_damage_isoDuctile_init(FILEUNIT)
+ if (any(phase_source == SOURCE_damage_isoDuctile_ID))       call source_damage_isoDuctile_init
- if (any(phase_source == SOURCE_damage_anisoBrittle_ID))     call source_damage_anisoBrittle_init(FILEUNIT)
+ if (any(phase_source == SOURCE_damage_anisoBrittle_ID))     call source_damage_anisoBrittle_init
- if (any(phase_source == SOURCE_damage_anisoDuctile_ID))     call source_damage_anisoDuctile_init(FILEUNIT)
+ if (any(phase_source == SOURCE_damage_anisoDuctile_ID))     call source_damage_anisoDuctile_init
-
+ 
 !--------------------------------------------------------------------------------------------------
 ! parse kinematic mechanisms from config file
  call IO_checkAndRewind(FILEUNIT)
- if (any(phase_kinematics == KINEMATICS_cleavage_opening_ID))  call kinematics_cleavage_opening_init(FILEUNIT)
+ if (any(phase_kinematics == KINEMATICS_cleavage_opening_ID))  call kinematics_cleavage_opening_init
- if (any(phase_kinematics == KINEMATICS_slipplane_opening_ID)) call kinematics_slipplane_opening_init(FILEUNIT)
+ if (any(phase_kinematics == KINEMATICS_slipplane_opening_ID)) call kinematics_slipplane_opening_init
  if (any(phase_kinematics == KINEMATICS_thermal_expansion_ID)) call kinematics_thermal_expansion_init(FILEUNIT)
  close(FILEUNIT)
 
@@ -602,9 +602,9 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, e
  KinematicsLoop: do k = 1_pInt, phase_Nkinematics(material_phase(ipc,ip,el))
    kinematicsType: select case (phase_kinematics(k,material_phase(ipc,ip,el)))
      case (KINEMATICS_cleavage_opening_ID) kinematicsType
-       call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S6, ipc, ip, el)
+       call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, math_6toSym33(S6), ipc, ip, el)
      case (KINEMATICS_slipplane_opening_ID) kinematicsType
-       call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S6, ipc, ip, el)
+       call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, math_6toSym33(S6), ipc, ip, el)
      case (KINEMATICS_thermal_expansion_ID) kinematicsType
        call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, ipc, ip, el)
      case default kinematicsType
@@ -891,7 +891,7 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac
    sourceType: select case (phase_source(s,material_phase(ipc,ip,el)))
 
      case (SOURCE_damage_anisoBrittle_ID) sourceType
-       call source_damage_anisoBrittle_dotState (S6, ipc, ip, el) !< correct stress?
+       call source_damage_anisoBrittle_dotState (math_6toSym33(S6), ipc, ip, el) !< correct stress?
 
      case (SOURCE_damage_isoDuctile_ID) sourceType
        call source_damage_isoDuctile_dotState   (         ipc, ip, el)
@@ -1112,16 +1112,18 @@ function constitutive_postResults(S, Fi, FeArray, ipc, ip, el)
  SourceLoop: do i = 1_pInt, phase_Nsources(material_phase(ipc,ip,el))
    startPos = endPos + 1_pInt
    endPos = endPos + sourceState(material_phase(ipc,ip,el))%p(i)%sizePostResults
+   of = phasememberAt(ipc,ip,el)
    sourceType: select case (phase_source(i,material_phase(ipc,ip,el)))
      case (SOURCE_damage_isoBrittle_ID) sourceType
-       constitutive_postResults(startPos:endPos) = source_damage_isoBrittle_postResults(ipc, ip, el)
+       constitutive_postResults(startPos:endPos) = source_damage_isoBrittle_postResults(material_phase(ipc,ip,el),of)
      case (SOURCE_damage_isoDuctile_ID) sourceType
-       constitutive_postResults(startPos:endPos) = source_damage_isoDuctile_postResults(ipc, ip, el)
+       constitutive_postResults(startPos:endPos) = source_damage_isoDuctile_postResults(material_phase(ipc,ip,el),of)
      case (SOURCE_damage_anisoBrittle_ID) sourceType
-       constitutive_postResults(startPos:endPos) = source_damage_anisoBrittle_postResults(ipc, ip, el)
+       constitutive_postResults(startPos:endPos) = source_damage_anisoBrittle_postResults(material_phase(ipc,ip,el),of)
      case (SOURCE_damage_anisoDuctile_ID) sourceType
-       constitutive_postResults(startPos:endPos) = source_damage_anisoDuctile_postResults(ipc, ip, el)
+       constitutive_postResults(startPos:endPos) = source_damage_anisoDuctile_postResults(material_phase(ipc,ip,el),of)
    end select sourceType
+
  enddo SourceLoop
 
 end function constitutive_postResults

src/damage_local.f90

@@ -225,6 +225,7 @@ subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el
    homogenization_Ngrains, &
    mappingHomogenization, &
    phaseAt, &
+   phasememberAt, &
    phase_source, &
    phase_Nsources, &
    SOURCE_damage_isoBrittle_ID, &
@@ -249,7 +250,8 @@ subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el
  integer(pInt) :: &
    phase, &
    grain, &
-   source
+   source, &
+   constituent
  real(pReal) :: &
    phiDot, dPhiDot_dPhi, localphiDot, dLocalphiDot_dPhi  
 
@@ -257,19 +259,20 @@ subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el
  dPhiDot_dPhi = 0.0_pReal
  do grain = 1, homogenization_Ngrains(mappingHomogenization(2,ip,el))
    phase = phaseAt(grain,ip,el)
+   constituent = phasememberAt(grain,ip,el)
    do source = 1, phase_Nsources(phase)
      select case(phase_source(source,phase))                                                   
        case (SOURCE_damage_isoBrittle_ID)
-        call source_damage_isobrittle_getRateAndItsTangent  (localphiDot, dLocalphiDot_dPhi, phi, grain, ip,  el)
+        call source_damage_isobrittle_getRateAndItsTangent  (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
 
        case (SOURCE_damage_isoDuctile_ID)
-        call source_damage_isoductile_getRateAndItsTangent  (localphiDot, dLocalphiDot_dPhi, phi, grain, ip,  el)
+        call source_damage_isoductile_getRateAndItsTangent  (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
 
        case (SOURCE_damage_anisoBrittle_ID)
-        call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, grain, ip,  el)
+        call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
 
        case (SOURCE_damage_anisoDuctile_ID)
-        call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, grain, ip,  el)
+        call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
 
        case default
         localphiDot = 0.0_pReal

src/damage_nonlocal.f90

@@ -186,6 +186,7 @@ subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip,
    homogenization_Ngrains, &
    mappingHomogenization, &
    phaseAt, &
+   phasememberAt, &
    phase_source, &
    phase_Nsources, &
    SOURCE_damage_isoBrittle_ID, &
@@ -210,7 +211,8 @@ subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip,
  integer(pInt) :: &
    phase, &
    grain, &
-   source
+   source, &
+   constituent
  real(pReal) :: &
    phiDot, dPhiDot_dPhi, localphiDot, dLocalphiDot_dPhi  
 
@@ -218,19 +220,20 @@ subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip,
  dPhiDot_dPhi = 0.0_pReal
  do grain = 1, homogenization_Ngrains(mappingHomogenization(2,ip,el))
    phase = phaseAt(grain,ip,el)
+   constituent = phasememberAt(grain,ip,el)
    do source = 1_pInt, phase_Nsources(phase)
      select case(phase_source(source,phase))                                                   
        case (SOURCE_damage_isoBrittle_ID)
-        call source_damage_isobrittle_getRateAndItsTangent  (localphiDot, dLocalphiDot_dPhi, phi, grain, ip,  el)
+        call source_damage_isobrittle_getRateAndItsTangent  (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
 
        case (SOURCE_damage_isoDuctile_ID)
-        call source_damage_isoductile_getRateAndItsTangent  (localphiDot, dLocalphiDot_dPhi, phi, grain, ip,  el)
+        call source_damage_isoductile_getRateAndItsTangent  (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
 
        case (SOURCE_damage_anisoBrittle_ID)
-        call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, grain, ip,  el)
+        call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
 
        case (SOURCE_damage_anisoDuctile_ID)
-        call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, grain, ip,  el)
+        call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
 
        case default
         localphiDot = 0.0_pReal

src/kinematics_cleavage_opening.f90

@@ -11,20 +11,22 @@ module kinematics_cleavage_opening
 
  implicit none
  private
- integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
+ integer(pInt), dimension(:), allocatable, private :: kinematics_cleavage_opening_instance
-   kinematics_cleavage_opening_sizePostResults, &                                                                !< cumulative size of post results
-   kinematics_cleavage_opening_offset, &                                                                         !< which kinematics is my current damage mechanism?
-   kinematics_cleavage_opening_instance                                                                          !< instance of damage kinematics mechanism
 
- integer(pInt),                       dimension(:,:),         allocatable, target, public  :: &
+ type, private :: tParameters                                                                       !< container type for internal constitutive parameters
-   kinematics_cleavage_opening_sizePostResult                                                                    !< size of each post result output
+   integer(pInt) :: &
-
+     totalNcleavage
- character(len=64),                   dimension(:,:),         allocatable, target, public  :: &
+   integer(pInt), dimension(:),   allocatable :: &
-   kinematics_cleavage_opening_output                                                                            !< name of each post result output
+     Ncleavage                                                                                      !< active number of cleavage systems per family
-   
+   real(pReal) :: &
- integer(pInt),                       dimension(:),           allocatable, target, public  :: &
+     sdot0, &
-   kinematics_cleavage_opening_Noutput                                                                           !< number of outputs per instance of this damage 
+     n
+   real(pReal),   dimension(:),   allocatable :: &
+     critDisp, &
+     critLoad
+ end type
 
+! Begin Deprecated
  integer(pInt),                       dimension(:),           allocatable,         private :: &
    kinematics_cleavage_opening_totalNcleavage                                                                    !< total number of cleavage systems
    
@@ -38,6 +40,7 @@ module kinematics_cleavage_opening
  real(pReal),                         dimension(:,:),         allocatable,         private :: &
    kinematics_cleavage_opening_critDisp, &
    kinematics_cleavage_opening_critLoad
+! End Deprecated
 
  public :: &
    kinematics_cleavage_opening_init, &
@@ -50,7 +53,7 @@ contains
 !> @brief module initialization
 !> @details reads in material parameters, allocates arrays, and does sanity checks
 !--------------------------------------------------------------------------------------------------
-subroutine kinematics_cleavage_opening_init(fileUnit)
+subroutine kinematics_cleavage_opening_init()
 #if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
  use, intrinsic :: iso_fortran_env, only: &
    compiler_version, &
@@ -60,41 +63,25 @@ subroutine kinematics_cleavage_opening_init(fileUnit)
    debug_level,&
    debug_constitutive,&
    debug_levelBasic
+ use config, only: &
+   config_phase
  use IO, only: &
-   IO_read, &
-   IO_lc, &
-   IO_getTag, &
-   IO_isBlank, &
-   IO_stringPos, &
-   IO_stringValue, &
-   IO_floatValue, &
-   IO_intValue, &
    IO_warning, &
    IO_error, &
-   IO_timeStamp, &
+   IO_timeStamp
-   IO_EOF
  use material, only: &
    phase_kinematics, &
-   phase_Nkinematics, &
-   phase_Noutput, &
    KINEMATICS_cleavage_opening_label, &
    KINEMATICS_cleavage_opening_ID
- use config, only: &
-   material_Nphase, &
-   MATERIAL_partPhase
  use lattice, only: &
    lattice_maxNcleavageFamily, &
    lattice_NcleavageSystem
 
  implicit none
- integer(pInt), intent(in) :: fileUnit
+ integer(pInt), allocatable, dimension(:) :: tempInt
+ real(pReal), allocatable, dimension(:) :: tempFloat
 
- integer(pInt), allocatable, dimension(:) :: chunkPos
+ integer(pInt) :: maxNinstance,p,instance,kinematics
- integer(pInt) :: maxNinstance,phase,instance,kinematics
- integer(pInt) :: Nchunks_CleavageFamilies = 0_pInt, j   
- character(len=65536) :: &
-   tag  = '', &
-   line = ''
 
  write(6,'(/,a)')   ' <<<+-  kinematics_'//KINEMATICS_cleavage_opening_LABEL//' init  -+>>>'
  write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
@@ -106,21 +93,11 @@ subroutine kinematics_cleavage_opening_init(fileUnit)
  if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
    write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
  
- allocate(kinematics_cleavage_opening_offset(material_Nphase), source=0_pInt)
+ allocate(kinematics_cleavage_opening_instance(size(config_phase)), source=0_pInt)
- allocate(kinematics_cleavage_opening_instance(material_Nphase), source=0_pInt)
+ do p = 1_pInt, size(config_phase)
- do phase = 1, material_Nphase
+   kinematics_cleavage_opening_instance(p) = count(phase_kinematics(:,1:p) == kinematics_cleavage_opening_ID) ! ToDo: count correct?
-   kinematics_cleavage_opening_instance(phase) = count(phase_kinematics(:,1:phase) == kinematics_cleavage_opening_ID)
-   do kinematics = 1, phase_Nkinematics(phase)
-     if (phase_kinematics(kinematics,phase) == kinematics_cleavage_opening_ID) &
-       kinematics_cleavage_opening_offset(phase) = kinematics
-   enddo    
  enddo
    
- allocate(kinematics_cleavage_opening_sizePostResults(maxNinstance),                      source=0_pInt)
- allocate(kinematics_cleavage_opening_sizePostResult(maxval(phase_Noutput),maxNinstance), source=0_pInt)
- allocate(kinematics_cleavage_opening_output(maxval(phase_Noutput),maxNinstance))
-          kinematics_cleavage_opening_output = ''
- allocate(kinematics_cleavage_opening_Noutput(maxNinstance),                              source=0_pInt) 
  allocate(kinematics_cleavage_opening_critDisp(lattice_maxNcleavageFamily,maxNinstance),  source=0.0_pReal) 
  allocate(kinematics_cleavage_opening_critLoad(lattice_maxNcleavageFamily,maxNinstance),  source=0.0_pReal) 
  allocate(kinematics_cleavage_opening_Ncleavage(lattice_maxNcleavageFamily,maxNinstance), source=0_pInt)
@@ -128,90 +105,51 @@ subroutine kinematics_cleavage_opening_init(fileUnit)
  allocate(kinematics_cleavage_opening_sdot_0(maxNinstance),                               source=0.0_pReal) 
  allocate(kinematics_cleavage_opening_N(maxNinstance),                                    source=0.0_pReal) 
 
- rewind(fileUnit)
+ do p = 1_pInt, size(config_phase)
- phase = 0_pInt
+   if (all(phase_kinematics(:,p) /= KINEMATICS_cleavage_opening_ID)) cycle
- do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase)         ! wind forward to <phase>
+   instance = kinematics_cleavage_opening_instance(p)
-   line = IO_read(fileUnit)
+   kinematics_cleavage_opening_sdot_0(instance) = config_phase(p)%getFloat('anisobrittle_sdot0')
- enddo
+   kinematics_cleavage_opening_N(instance)      = config_phase(p)%getFloat('anisobrittle_ratesensitivity')
- 
+   tempInt = config_phase(p)%getInts('ncleavage')     
- parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
+   kinematics_cleavage_opening_Ncleavage(1:size(tempInt),instance) = tempInt
-   line = IO_read(fileUnit)
-   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
-   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
-     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
-     exit                                                                                           
-   endif   
-   if (IO_getTag(line,'[',']') /= '') then                                                          ! next phase section
-     phase = phase + 1_pInt                                                                         ! advance phase section counter
-     cycle                                                                                          ! skip to next line
-   endif
-   if (phase > 0_pInt ) then; if (any(phase_kinematics(:,phase) == KINEMATICS_cleavage_opening_ID)) then         ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
-     instance = kinematics_cleavage_opening_instance(phase)                                                         ! which instance of my damage is present phase
-     chunkPos = IO_stringPos(line)
-     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
-     select case(tag)
-       case ('anisobrittle_sdot0')
-         kinematics_cleavage_opening_sdot_0(instance) = IO_floatValue(line,chunkPos,2_pInt)
-         
-       case ('anisobrittle_ratesensitivity')
-         kinematics_cleavage_opening_N(instance) = IO_floatValue(line,chunkPos,2_pInt)
-         
-       case ('ncleavage')  !
-         Nchunks_CleavageFamilies = chunkPos(1) - 1_pInt
-         do j = 1_pInt, Nchunks_CleavageFamilies
-           kinematics_cleavage_opening_Ncleavage(j,instance) = IO_intValue(line,chunkPos,1_pInt+j)
-         enddo
 
-       case ('anisobrittle_criticaldisplacement')
+   tempFloat = config_phase(p)%getFloats('anisobrittle_criticaldisplacement',requiredSize=size(tempInt))
-         do j = 1_pInt, Nchunks_CleavageFamilies
+   kinematics_cleavage_opening_critDisp(1:size(tempInt),instance) = tempFloat
-           kinematics_cleavage_opening_critDisp(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j)
-         enddo
 
-       case ('anisobrittle_criticalload')
+   tempFloat = config_phase(p)%getFloats('anisobrittle_criticalload',requiredSize=size(tempInt))
-         do j = 1_pInt, Nchunks_CleavageFamilies
+   kinematics_cleavage_opening_critLoad(1:size(tempInt),instance) = tempFloat
-           kinematics_cleavage_opening_critLoad(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j)
-         enddo
 
-     end select
+   kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance) = &
-   endif; endif
+       min(lattice_NcleavageSystem(1:lattice_maxNcleavageFamily,p),&                            ! limit active cleavage systems per family to min of available and requested
- enddo parsingFile
-
-!--------------------------------------------------------------------------------------------------
-!  sanity checks
- sanityChecks: do phase = 1_pInt, material_Nphase   
-   myPhase: if (any(phase_kinematics(:,phase) == KINEMATICS_cleavage_opening_ID)) then
-     instance = kinematics_cleavage_opening_instance(phase)
-     kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance) = &
-       min(lattice_NcleavageSystem(1:lattice_maxNcleavageFamily,phase),&                            ! limit active cleavage systems per family to min of available and requested
            kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance))
      kinematics_cleavage_opening_totalNcleavage(instance)  = sum(kinematics_cleavage_opening_Ncleavage(:,instance)) ! how many cleavage systems altogether
      if (kinematics_cleavage_opening_sdot_0(instance) <= 0.0_pReal) &
        call IO_error(211_pInt,el=instance,ext_msg='sdot_0 ('//KINEMATICS_cleavage_opening_LABEL//')')
-     if (any(kinematics_cleavage_opening_critDisp(1:Nchunks_CleavageFamilies,instance) < 0.0_pReal)) &
+     if (any(kinematics_cleavage_opening_critDisp(1:size(tempInt),instance) < 0.0_pReal)) &
        call IO_error(211_pInt,el=instance,ext_msg='critical_displacement ('//KINEMATICS_cleavage_opening_LABEL//')')
-     if (any(kinematics_cleavage_opening_critLoad(1:Nchunks_CleavageFamilies,instance) < 0.0_pReal)) &
+     if (any(kinematics_cleavage_opening_critLoad(1:size(tempInt),instance) < 0.0_pReal)) &
        call IO_error(211_pInt,el=instance,ext_msg='critical_load ('//KINEMATICS_cleavage_opening_LABEL//')')
      if (kinematics_cleavage_opening_N(instance) <= 0.0_pReal) &
        call IO_error(211_pInt,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_cleavage_opening_LABEL//')')
-   endif myPhase
+ enddo
- enddo sanityChecks
  
 end subroutine kinematics_cleavage_opening_init
 
 !--------------------------------------------------------------------------------------------------
 !> @brief  contains the constitutive equation for calculating the velocity gradient  
 !--------------------------------------------------------------------------------------------------
-subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar3333, Tstar_v, ipc, ip, el)
+subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el)
  use prec, only: &
    tol_math_check
+ use math, only: &
+   math_mul33xx33
  use material, only: &
-   phaseAt, phasememberAt, &
+   material_phase, &
    material_homog, &
    damage, &
    damageMapping
  use lattice, only: &
    lattice_Scleavage, &
-   lattice_Scleavage_v, &
    lattice_maxNcleavageFamily, &
    lattice_NcleavageSystem
  
@@ -220,36 +158,33 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar3333, Tstar
    ipc, &                                                                                           !< grain number
    ip, &                                                                                            !< integration point number
    el                                                                                               !< element number
- real(pReal),   intent(in),  dimension(6) :: &
+ real(pReal),   intent(in),  dimension(3,3) :: &
-   Tstar_v                                                                                          !< 2nd Piola-Kirchhoff stress
+   S
  real(pReal),   intent(out), dimension(3,3) :: &
    Ld                                                                                               !< damage velocity gradient
  real(pReal),   intent(out), dimension(3,3,3,3) :: &
-   dLd_dTstar3333                                                                                   !< derivative of Ld with respect to Tstar (4th-order tensor)
+   dLd_dTstar                                                                                       !< derivative of Ld with respect to Tstar (4th-order tensor)
  integer(pInt) :: &
-   phase, &
+   instance, phase, &
-   constituent, &
-   instance, &
    homog, damageOffset, &
    f, i, index_myFamily, k, l, m, n
  real(pReal) :: &
    traction_d, traction_t, traction_n, traction_crit, &
    udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt
 
- phase = phaseAt(ipc,ip,el)
+ phase = material_phase(ipc,ip,el)
- constituent = phasememberAt(ipc,ip,el)
  instance = kinematics_cleavage_opening_instance(phase)
  homog = material_homog(ip,el)
  damageOffset = damageMapping(homog)%p(ip,el)
  
  Ld = 0.0_pReal
- dLd_dTstar3333 = 0.0_pReal
+ dLd_dTstar = 0.0_pReal
  do f = 1_pInt,lattice_maxNcleavageFamily
    index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase))                                   ! at which index starts my family
    do i = 1_pInt,kinematics_cleavage_opening_Ncleavage(f,instance)                                            ! process each (active) cleavage system in family
-     traction_d    = dot_product(Tstar_v,lattice_Scleavage_v(1:6,1,index_myFamily+i,phase))
+     traction_d    = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase))
-     traction_t    = dot_product(Tstar_v,lattice_Scleavage_v(1:6,2,index_myFamily+i,phase))
+     traction_t    = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase))
-     traction_n    = dot_product(Tstar_v,lattice_Scleavage_v(1:6,3,index_myFamily+i,phase))
+     traction_n    = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase))
      traction_crit = kinematics_cleavage_opening_critLoad(f,instance)* &
                      damage(homog)%p(damageOffset)*damage(homog)%p(damageOffset)
      udotd = &
@@ -261,7 +196,7 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar3333, Tstar
        dudotd_dt = sign(1.0_pReal,traction_d)*udotd*kinematics_cleavage_opening_N(instance)/ &
                    max(0.0_pReal, abs(traction_d) - traction_crit)
        forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
-         dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + &
+         dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
            dudotd_dt*lattice_Scleavage(k,l,1,index_myFamily+i,phase)* &
                      lattice_Scleavage(m,n,1,index_myFamily+i,phase)
      endif                
@@ -275,7 +210,7 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar3333, Tstar
        dudott_dt = sign(1.0_pReal,traction_t)*udott*kinematics_cleavage_opening_N(instance)/ &
                    max(0.0_pReal, abs(traction_t) - traction_crit)  
        forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
-         dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + &
+         dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
            dudott_dt*lattice_Scleavage(k,l,2,index_myFamily+i,phase)* &
                      lattice_Scleavage(m,n,2,index_myFamily+i,phase)
      endif                
@@ -289,11 +224,10 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar3333, Tstar
        dudotn_dt = sign(1.0_pReal,traction_n)*udotn*kinematics_cleavage_opening_N(instance)/ &
                    max(0.0_pReal, abs(traction_n) - traction_crit)
        forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
-         dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + &
+         dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
            dudotn_dt*lattice_Scleavage(k,l,3,index_myFamily+i,phase)* &
                      lattice_Scleavage(m,n,3,index_myFamily+i,phase)
      endif                
-
    enddo
  enddo
  

src/kinematics_slipplane_opening.f90

@@ -11,20 +11,22 @@ module kinematics_slipplane_opening
 
  implicit none
  private
- integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
+ integer(pInt), dimension(:), allocatable, private :: kinematics_slipplane_opening_instance
-   kinematics_slipplane_opening_sizePostResults, &                                                            !< cumulative size of post results
-   kinematics_slipplane_opening_offset, &                                                                     !< which kinematics is my current damage mechanism?
-   kinematics_slipplane_opening_instance                                                                      !< instance of damage kinematics mechanism
-
- integer(pInt),                       dimension(:,:),         allocatable, target, public  :: &
-   kinematics_slipplane_opening_sizePostResult                                                                !< size of each post result output
-
- character(len=64),                   dimension(:,:),         allocatable, target, public  :: &
-   kinematics_slipplane_opening_output                                                                        !< name of each post result output
-   
- integer(pInt),                       dimension(:),           allocatable, target, public  :: &
-   kinematics_slipplane_opening_Noutput                                                                       !< number of outputs per instance of this damage 
    
+ type, private :: tParameters                                                                       !< container type for internal constitutive parameters
+   integer(pInt) :: &
+     totalNslip
+   integer(pInt), dimension(:),   allocatable :: &
+     Nslip                                                                                          !< active number of slip systems per family
+   real(pReal) :: &
+     sdot0, &
+     n
+   real(pReal),   dimension(:),   allocatable :: &
+     critDisp, &
+     critPlasticStrain
+ end type
+
+! Begin Deprecated
  integer(pInt),                       dimension(:),           allocatable,         private :: &
    kinematics_slipplane_opening_totalNslip                                                                    !< total number of slip systems
 
@@ -38,6 +40,7 @@ module kinematics_slipplane_opening
  real(pReal),                         dimension(:,:),         allocatable,         private :: &
    kinematics_slipplane_opening_critPlasticStrain, &
    kinematics_slipplane_opening_critLoad
+! End Deprecated
    
  public :: &
    kinematics_slipplane_opening_init, &
@@ -50,7 +53,7 @@ contains
 !> @brief module initialization
 !> @details reads in material parameters, allocates arrays, and does sanity checks
 !--------------------------------------------------------------------------------------------------
-subroutine kinematics_slipplane_opening_init(fileUnit)
+subroutine kinematics_slipplane_opening_init()
 #if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
  use, intrinsic :: iso_fortran_env, only: &
    compiler_version, &
@@ -60,41 +63,25 @@ subroutine kinematics_slipplane_opening_init(fileUnit)
    debug_level,&
    debug_constitutive,&
    debug_levelBasic
+ use config, only: &
+   config_phase
  use IO, only: &
-   IO_read, &
-   IO_lc, &
-   IO_getTag, &
-   IO_isBlank, &
-   IO_stringPos, &
-   IO_stringValue, &
-   IO_floatValue, &
-   IO_intValue, &
    IO_warning, &
    IO_error, &
-   IO_timeStamp, &
+   IO_timeStamp
-   IO_EOF
  use material, only: &
    phase_kinematics, &
-   phase_Nkinematics, &
-   phase_Noutput, &
    KINEMATICS_slipplane_opening_label, &
    KINEMATICS_slipplane_opening_ID
- use config, only: &
-   material_Nphase, &
-   MATERIAL_partPhase
  use lattice, only: &
    lattice_maxNslipFamily, &
    lattice_NslipSystem
 
  implicit none
- integer(pInt), intent(in) :: fileUnit
+ integer(pInt), allocatable, dimension(:) :: tempInt
+ real(pReal), allocatable, dimension(:) :: tempFloat
 
- integer(pInt), allocatable, dimension(:) :: chunkPos
+ integer(pInt) :: maxNinstance,p,instance,kinematics
- integer(pInt) :: maxNinstance,phase,instance,kinematics
- integer(pInt) :: Nchunks_SlipFamilies = 0_pInt, j   
- character(len=65536) :: &
-   tag  = '', &
-   line = ''
 
  write(6,'(/,a)')   ' <<<+-  kinematics_'//KINEMATICS_slipplane_opening_LABEL//' init  -+>>>'
  write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
@@ -106,21 +93,11 @@ subroutine kinematics_slipplane_opening_init(fileUnit)
  if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
    write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
  
- allocate(kinematics_slipplane_opening_offset(material_Nphase), source=0_pInt)
+ allocate(kinematics_slipplane_opening_instance(size(config_phase)), source=0_pInt)
- allocate(kinematics_slipplane_opening_instance(material_Nphase), source=0_pInt)
+ do p = 1_pInt, size(config_phase)
- do phase = 1, material_Nphase
+   kinematics_slipplane_opening_instance(p) = count(phase_kinematics(:,1:p) == kinematics_slipplane_opening_ID) ! ToDo: count correct?
-   kinematics_slipplane_opening_instance(phase) = count(phase_kinematics(:,1:phase) == kinematics_slipplane_opening_ID)
-   do kinematics = 1, phase_Nkinematics(phase)
-     if (phase_kinematics(kinematics,phase) == kinematics_slipplane_opening_ID) &
-       kinematics_slipplane_opening_offset(phase) = kinematics
-   enddo    
  enddo
    
- allocate(kinematics_slipplane_opening_sizePostResults(maxNinstance),                     source=0_pInt)
- allocate(kinematics_slipplane_opening_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
- allocate(kinematics_slipplane_opening_output(maxval(phase_Noutput),maxNinstance))
-          kinematics_slipplane_opening_output = ''
- allocate(kinematics_slipplane_opening_Noutput(maxNinstance),                             source=0_pInt) 
  allocate(kinematics_slipplane_opening_critLoad(lattice_maxNslipFamily,maxNinstance),     source=0.0_pReal) 
  allocate(kinematics_slipplane_opening_critPlasticStrain(lattice_maxNslipFamily,maxNinstance),source=0.0_pReal) 
  allocate(kinematics_slipplane_opening_Nslip(lattice_maxNslipFamily,maxNinstance),        source=0_pInt)
@@ -128,61 +105,22 @@ subroutine kinematics_slipplane_opening_init(fileUnit)
  allocate(kinematics_slipplane_opening_N(maxNinstance),                                   source=0.0_pReal) 
  allocate(kinematics_slipplane_opening_sdot_0(maxNinstance),                              source=0.0_pReal) 
 
- rewind(fileUnit)
+ do p = 1_pInt, size(config_phase)
- phase = 0_pInt
+   if (all(phase_kinematics(:,p) /= KINEMATICS_slipplane_opening_ID)) cycle
- do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase)         ! wind forward to <phase>
+   instance = kinematics_slipplane_opening_instance(p)
-   line = IO_read(fileUnit)
+   kinematics_slipplane_opening_sdot_0(instance) = config_phase(p)%getFloat('anisoductile_sdot0')
- enddo
+   kinematics_slipplane_opening_N(instance)      = config_phase(p)%getFloat('anisoductile_ratesensitivity')
- 
+   tempInt = config_phase(p)%getInts('ncleavage')     
- parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
+   kinematics_slipplane_opening_Nslip(1:size(tempInt),instance) = tempInt
-   line = IO_read(fileUnit)
-   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
-   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
-     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
-     exit                                                                                           
-   endif   
-   if (IO_getTag(line,'[',']') /= '') then                                                          ! next phase section
-     phase = phase + 1_pInt                                                                         ! advance phase section counter
-     cycle                                                                                          ! skip to next line
-   endif
-   if (phase > 0_pInt ) then; if (any(phase_kinematics(:,phase) == KINEMATICS_slipplane_opening_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
-     instance = kinematics_slipplane_opening_instance(phase)                                        ! which instance of my damage is present phase
-     chunkPos = IO_stringPos(line)
-     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                              ! extract key
-     select case(tag)
-       case ('nslip')  !
-         Nchunks_SlipFamilies = chunkPos(1) - 1_pInt
-         do j = 1_pInt, Nchunks_SlipFamilies
-           kinematics_slipplane_opening_Nslip(j,instance) = IO_intValue(line,chunkPos,1_pInt+j)
-         enddo
 
-       case ('anisoductile_sdot0')
+   tempFloat = config_phase(p)%getFloats('anisoductile_criticalplasticstrain',requiredSize=size(tempInt))
-         kinematics_slipplane_opening_sdot_0(instance) = IO_floatValue(line,chunkPos,2_pInt)
+   kinematics_slipplane_opening_critPlasticStrain(1:size(tempInt),instance) = tempFloat
-         
-       case ('anisoductile_criticalplasticstrain')
-         do j = 1_pInt, Nchunks_SlipFamilies
-           kinematics_slipplane_opening_critPlasticStrain(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j)
-         enddo
-         
-       case ('anisoductile_ratesensitivity')
-         kinematics_slipplane_opening_N(instance) = IO_floatValue(line,chunkPos,2_pInt)
 
-       case ('anisoductile_criticalload')
+   tempFloat = config_phase(p)%getFloats('anisoductile_criticalload',requiredSize=size(tempInt))
-         do j = 1_pInt, Nchunks_SlipFamilies
+   kinematics_slipplane_opening_critLoad(1:size(tempInt),instance) = tempFloat
-           kinematics_slipplane_opening_critLoad(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j)
-         enddo
-         
-     end select
-   endif; endif
- enddo parsingFile
 
-!--------------------------------------------------------------------------------------------------
-!  sanity checks
- sanityChecks: do phase = 1_pInt, material_Nphase   
-   myPhase: if (any(phase_kinematics(:,phase) == KINEMATICS_slipplane_opening_ID)) then
-     instance = kinematics_slipplane_opening_instance(phase)
      kinematics_slipplane_opening_Nslip(1:lattice_maxNslipFamily,instance) = &
-       min(lattice_NslipSystem(1:lattice_maxNslipFamily,phase),&                                    ! limit active cleavage systems per family to min of available and requested
+       min(lattice_NslipSystem(1:lattice_maxNslipFamily,p),&                                    ! limit active cleavage systems per family to min of available and requested
            kinematics_slipplane_opening_Nslip(1:lattice_maxNslipFamily,instance))
          kinematics_slipplane_opening_totalNslip(instance) = sum(kinematics_slipplane_opening_Nslip(:,instance))
      if (kinematics_slipplane_opening_sdot_0(instance) <= 0.0_pReal) &
@@ -191,18 +129,18 @@ subroutine kinematics_slipplane_opening_init(fileUnit)
        call IO_error(211_pInt,el=instance,ext_msg='criticaPlasticStrain ('//KINEMATICS_slipplane_opening_LABEL//')')
      if (kinematics_slipplane_opening_N(instance) <= 0.0_pReal) &
        call IO_error(211_pInt,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_slipplane_opening_LABEL//')')
-   endif myPhase
+ enddo
- enddo sanityChecks
   
-
 end subroutine kinematics_slipplane_opening_init
 
 !--------------------------------------------------------------------------------------------------
 !> @brief  contains the constitutive equation for calculating the velocity gradient  
 !--------------------------------------------------------------------------------------------------
-subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar3333, Tstar_v, ipc, ip, el)
+subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el)
  use prec, only: &
    tol_math_check
+ use math, only: &
+   math_mul33xx33
  use lattice, only: &
    lattice_maxNslipFamily, &
    lattice_NslipSystem, &
@@ -210,53 +148,41 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar3333, Tsta
    lattice_st, &
    lattice_sn
  use material, only: &
-   phaseAt, phasememberAt, &
+   material_phase, &
    material_homog, &
    damage, &
    damageMapping
  use math, only: &
-   math_Plain3333to99, &
+   math_tensorproduct33
-   math_I3, &
-   math_identity4th, &
-   math_symmetric33, &
-   math_Mandel33to6, &
-   math_tensorproduct33, &
-   math_det33, &
-   math_mul33x33
  
  implicit none
  integer(pInt), intent(in) :: &
    ipc, &                                                                                           !< grain number
    ip, &                                                                                            !< integration point number
    el                                                                                               !< element number
- real(pReal),   intent(in),  dimension(6) :: &
+ real(pReal),   intent(in),  dimension(3,3) :: &
-   Tstar_v                                                                                          !< 2nd Piola-Kirchhoff stress
+   S
  real(pReal),   intent(out), dimension(3,3) :: &
    Ld                                                                                               !< damage velocity gradient
  real(pReal),   intent(out), dimension(3,3,3,3) :: &
-   dLd_dTstar3333                                                                                   !< derivative of Ld with respect to Tstar (4th-order tensor)
+   dLd_dTstar                                                                                       !< derivative of Ld with respect to Tstar (4th-order tensor)
  real(pReal),   dimension(3,3) :: &
    projection_d, projection_t, projection_n                                                         !< projection modes 3x3 tensor
- real(pReal),   dimension(6) :: &
-   projection_d_v, projection_t_v, projection_n_v                                                   !< projection modes 3x3 vector
  integer(pInt) :: &
-   phase, &
+   instance, phase, &
-   constituent, &
-   instance, &
    homog, damageOffset, &
    f, i, index_myFamily, k, l, m, n
  real(pReal) :: &
    traction_d, traction_t, traction_n, traction_crit, &
    udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt
    
- phase = phaseAt(ipc,ip,el)
+ phase = material_phase(ipc,ip,el)
- constituent = phasememberAt(ipc,ip,el)
  instance = kinematics_slipplane_opening_instance(phase)
  homog = material_homog(ip,el)
  damageOffset = damageMapping(homog)%p(ip,el)
  
  Ld = 0.0_pReal
- dLd_dTstar3333 = 0.0_pReal
+ dLd_dTstar = 0.0_pReal
  do f = 1_pInt,lattice_maxNslipFamily
    index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,phase))                                      ! at which index starts my family
    do i = 1_pInt,kinematics_slipplane_opening_Nslip(f,instance)                                              ! process each (active) slip system in family
@@ -267,13 +193,10 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar3333, Tsta
      projection_n = math_tensorproduct33(lattice_sn(1:3,index_myFamily+i,phase),&
                                        lattice_sn(1:3,index_myFamily+i,phase))
 
-     projection_d_v(1:6) = math_Mandel33to6(math_symmetric33(projection_d(1:3,1:3)))
-     projection_t_v(1:6) = math_Mandel33to6(math_symmetric33(projection_t(1:3,1:3)))
-     projection_n_v(1:6) = math_Mandel33to6(math_symmetric33(projection_n(1:3,1:3)))
    
-     traction_d    = dot_product(Tstar_v,projection_d_v(1:6))
+     traction_d    = math_mul33xx33(S,projection_d)
-     traction_t    = dot_product(Tstar_v,projection_t_v(1:6))
+     traction_t    = math_mul33xx33(S,projection_t)
-     traction_n    = dot_product(Tstar_v,projection_n_v(1:6))
+     traction_n    = math_mul33xx33(S,projection_n)
      
      traction_crit = kinematics_slipplane_opening_critLoad(f,instance)* &
                      damage(homog)%p(damageOffset)                                                        ! degrading critical load carrying capacity by damage 
@@ -287,7 +210,7 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar3333, Tsta
        Ld = Ld + udotd*projection_d
        dudotd_dt = udotd*kinematics_slipplane_opening_N(instance)/traction_d
        forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
-         dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + &
+         dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
            dudotd_dt*projection_d(k,l)*projection_d(m,n)
      endif                
  
@@ -300,9 +223,10 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar3333, Tsta
        Ld = Ld + udott*projection_t
        dudott_dt = udott*kinematics_slipplane_opening_N(instance)/traction_t
        forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
-         dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + &
+         dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
            dudott_dt*projection_t(k,l)*projection_t(m,n)
      endif
+
      udotn = &
        kinematics_slipplane_opening_sdot_0(instance)* &
        (max(0.0_pReal,traction_n)/traction_crit - &
@@ -311,7 +235,7 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar3333, Tsta
        Ld = Ld + udotn*projection_n
        dudotn_dt = udotn*kinematics_slipplane_opening_N(instance)/traction_n
        forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
-         dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + &
+         dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) + &
            dudotn_dt*projection_n(k,l)*projection_n(m,n)
      endif 
    enddo

src/lattice.f90

@@ -26,9 +26,6 @@ module lattice
  real(pReal), allocatable, dimension(:,:,:,:,:), protected, public :: &
    lattice_Scleavage                                                                                !< Schmid matrices for cleavage systems
 
- real(pReal), allocatable, dimension(:,:,:,:), protected, public :: &
-   lattice_Scleavage_v                                                                              !< Mandel notation of lattice_Scleavege
-
  real(pReal), allocatable, dimension(:,:,:), protected, public :: &
    lattice_sn, &                                                                                    !< normal direction of slip system
    lattice_st, &                                                                                    !< sd x sn
@@ -569,7 +566,6 @@ subroutine lattice_init
  allocate(lattice_NslipSystem(lattice_maxNslipFamily,Nphases),source=0_pInt)
 
  allocate(lattice_Scleavage(3,3,3,lattice_maxNslip,Nphases),source=0.0_pReal)
- allocate(lattice_Scleavage_v(6,3,lattice_maxNslip,Nphases),source=0.0_pReal)
  allocate(lattice_NcleavageSystem(lattice_maxNcleavageFamily,Nphases),source=0_pInt)
 
  allocate(CoverA(Nphases),source=0.0_pReal)
@@ -657,8 +653,6 @@ subroutine lattice_initializeStructure(myPhase,CoverA)
   tol_math_check
  use math, only: &
    math_mul33x33, &
-   math_symmetric33, &
-   math_sym33to6, &
    math_sym3333to66, &
    math_Voigt66to3333, &
    math_crossproduct
@@ -805,13 +799,6 @@ subroutine lattice_initializeStructure(myPhase,CoverA)
    lattice_st(1:3,i,myPhase) = math_crossproduct(lattice_sd(1:3,i,myPhase),lattice_sn(1:3,i,myPhase))
  enddo
 
- do i = 1_pInt,myNcleavage                                                                          ! store slip system vectors and Schmid matrix for my structure
-   do j = 1_pInt,3_pInt
-     lattice_Scleavage_v(1:6,j,i,myPhase) = &
-       math_sym33to6(math_symmetric33(lattice_Scleavage(1:3,1:3,j,i,myPhase)))
-   enddo
- enddo
-
 end subroutine lattice_initializeStructure
 
 

src/material.f90

@@ -235,6 +235,7 @@ module material
  public :: &
    material_init, &
    material_allocatePlasticState, &
+   material_allocateSourceState, &
    ELASTICITY_hooke_ID ,&
    PLASTICITY_none_ID, &
    PLASTICITY_isotropic_ID, &
@@ -964,6 +965,49 @@ subroutine material_allocatePlasticState(phase,NofMyPhase,&
 end subroutine material_allocatePlasticState
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief allocates the source state of a phase
+!--------------------------------------------------------------------------------------------------
+subroutine material_allocateSourceState(phase,of,NofMyPhase,&
+                                        sizeState,sizeDotState,sizeDeltaState)
+ use numerics, only: &
+   numerics_integrator2 => numerics_integrator                                                      ! compatibility hack
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   phase, &
+   of, &
+   NofMyPhase, &
+   sizeState, sizeDotState,sizeDeltaState
+ integer(pInt) :: numerics_integrator                                                               ! compatibility hack
+ numerics_integrator = numerics_integrator2(1)                                                      ! compatibility hack
+
+ sourceState(phase)%p(of)%sizeState        = sizeState
+ sourceState(phase)%p(of)%sizeDotState     = sizeDotState
+ sourceState(phase)%p(of)%sizeDeltaState   = sizeDeltaState
+ plasticState(phase)%offsetDeltaState = sizeState-sizeDeltaState                                    ! deltaState occupies latter part of state by definition
+
+ allocate(sourceState(phase)%p(of)%aTolState           (sizeState),                source=0.0_pReal)
+ allocate(sourceState(phase)%p(of)%state0              (sizeState,NofMyPhase),     source=0.0_pReal)
+ allocate(sourceState(phase)%p(of)%partionedState0     (sizeState,NofMyPhase),     source=0.0_pReal)
+ allocate(sourceState(phase)%p(of)%subState0           (sizeState,NofMyPhase),     source=0.0_pReal)
+ allocate(sourceState(phase)%p(of)%state               (sizeState,NofMyPhase),     source=0.0_pReal)
+
+ allocate(sourceState(phase)%p(of)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
+ if (numerics_integrator == 1_pInt) then
+   allocate(sourceState(phase)%p(of)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)
+   allocate(sourceState(phase)%p(of)%previousDotState2 (sizeDotState,NofMyPhase),  source=0.0_pReal)
+ endif
+ if (numerics_integrator == 4_pInt) &
+   allocate(sourceState(phase)%p(of)%RK4dotState       (sizeDotState,NofMyPhase),  source=0.0_pReal)
+ if (numerics_integrator == 5_pInt) &
+   allocate(sourceState(phase)%p(of)%RKCK45dotState  (6,sizeDotState,NofMyPhase),  source=0.0_pReal)
+
+ allocate(plasticState(phase)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
+
+end subroutine material_allocateSourceState
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief populates the grains
 !> @details populates the grains by identifying active microstructure/homogenization pairs,

src/source_damage_anisoBrittle.f90

@@ -12,7 +12,6 @@ module source_damage_anisoBrittle
  implicit none
  private
  integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
-   source_damage_anisoBrittle_sizePostResults, &                                                                !< cumulative size of post results
    source_damage_anisoBrittle_offset, &                                                                         !< which source is my current source mechanism?
    source_damage_anisoBrittle_instance                                                                          !< instance of source mechanism
 
@@ -22,31 +21,32 @@ module source_damage_anisoBrittle
  character(len=64),                   dimension(:,:),         allocatable, target, public  :: &
    source_damage_anisoBrittle_output                                                                            !< name of each post result output
    
- integer(pInt),                       dimension(:),           allocatable, target, public  :: &
-   source_damage_anisoBrittle_Noutput                                                                           !< number of outputs per instance of this source 
-
- integer(pInt),                       dimension(:),           allocatable,         private :: &
-   source_damage_anisoBrittle_totalNcleavage                                                                    !< total number of cleavage systems
-   
  integer(pInt),                       dimension(:,:),         allocatable,         private :: &
    source_damage_anisoBrittle_Ncleavage                                                                         !< number of cleavage systems per family
-   
- real(pReal),                         dimension(:),           allocatable,         private :: &
-   source_damage_anisoBrittle_aTol, &
-   source_damage_anisoBrittle_sdot_0, &
-   source_damage_anisoBrittle_N
-
- real(pReal),                         dimension(:,:),         allocatable,         private :: &
-   source_damage_anisoBrittle_critDisp, &
-   source_damage_anisoBrittle_critLoad
 
  enum, bind(c) 
    enumerator :: undefined_ID, &
                  damage_drivingforce_ID
  end enum                                                
- 
+
- integer(kind(undefined_ID)),         dimension(:,:),         allocatable,          private :: & 
+
-   source_damage_anisoBrittle_outputID                                                                  !< ID of each post result output
+ type, private :: tParameters                                                                       !< container type for internal constitutive parameters
+   real(pReal) :: &
+     aTol, &
+     sdot_0, &
+     N
+   real(pReal), dimension(:), allocatable :: &
+     critDisp, &
+     critLoad
+   integer(pInt) :: &
+     totalNcleavage
+   integer(pInt), dimension(:), allocatable :: &
+     Ncleavage
+   integer(kind(undefined_ID)), allocatable, dimension(:) :: &
+     outputID                                                                                       !< ID of each post result output
+ end type tParameters
+
+ type(tParameters), dimension(:), allocatable, private :: param                                     !< containers of constitutive parameters (len Ninstance)
 
 
  public :: &
@@ -62,30 +62,24 @@ contains
 !> @brief module initialization
 !> @details reads in material parameters, allocates arrays, and does sanity checks
 !--------------------------------------------------------------------------------------------------
-subroutine source_damage_anisoBrittle_init(fileUnit)
+subroutine source_damage_anisoBrittle_init
 #if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
  use, intrinsic :: iso_fortran_env, only: &
    compiler_version, &
    compiler_options
 #endif
+ use prec, only: &
+   pStringLen
  use debug, only: &
    debug_level,&
    debug_constitutive,&
    debug_levelBasic
  use IO, only: &
-   IO_read, &
+   IO_error
-   IO_lc, &
+    use math, only: &
-   IO_getTag, &
+   math_expand
-   IO_isBlank, &
-   IO_stringPos, &
-   IO_stringValue, &
-   IO_floatValue, &
-   IO_intValue, &
-   IO_warning, &
-   IO_error, &
-   IO_timeStamp, &
-   IO_EOF
  use material, only: &
+   material_allocateSourceState, &
    phase_source, &
    phase_Nsources, &
    phase_Noutput, &
@@ -94,35 +88,34 @@ subroutine source_damage_anisoBrittle_init(fileUnit)
    material_phase, &
    sourceState
  use config, only: &
+   config_phase, &
    material_Nphase, &
    MATERIAL_partPhase
- use numerics,only: &
-   numerics_integrator
  use lattice, only: &
-   lattice_maxNcleavageFamily, &
+   lattice_maxNcleavageFamily
-   lattice_NcleavageSystem
 
  implicit none
- integer(pInt), intent(in) :: fileUnit
 
- integer(pInt), allocatable, dimension(:) :: chunkPos
+ integer(pInt) :: Ninstance,phase,instance,source,sourceOffset
- integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,source,sourceOffset,o
+ integer(pInt) :: NofMyPhase,p   ,i
- integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
+ integer(pInt),          dimension(0), parameter :: emptyIntArray    = [integer(pInt)::]
- integer(pInt) :: NofMyPhase   
+ character(len=65536),   dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
- integer(pInt) :: Nchunks_CleavageFamilies = 0_pInt, j   
+ integer(kind(undefined_ID)) :: &
- character(len=65536) :: &
+   outputID
-   tag  = '', &
-   line = ''
 
- write(6,'(/,a)')   ' <<<+-  source_'//SOURCE_damage_anisoBrittle_LABEL//' init  -+>>>'
+ character(len=pStringLen) :: &
- write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
+   extmsg = ''
+ character(len=65536), dimension(:), allocatable :: &
+   outputs
+
+ write(6,'(/,a)')   ' <<<+-  source_'//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//' init  -+>>>'
 #include "compilation_info.f90"
 
- maxNinstance = int(count(phase_source == SOURCE_damage_anisoBrittle_ID),pInt)
+ Ninstance = int(count(phase_source == SOURCE_damage_anisoBrittle_ID),pInt)
- if (maxNinstance == 0_pInt) return
+ if (Ninstance == 0_pInt) return
  
  if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
-   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
+   write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
  
  allocate(source_damage_anisoBrittle_offset(material_Nphase), source=0_pInt)
  allocate(source_damage_anisoBrittle_instance(material_Nphase), source=0_pInt)
@@ -133,160 +126,89 @@ subroutine source_damage_anisoBrittle_init(fileUnit)
        source_damage_anisoBrittle_offset(phase) = source
    enddo    
  enddo
-   
- allocate(source_damage_anisoBrittle_sizePostResults(maxNinstance),                      source=0_pInt)
- allocate(source_damage_anisoBrittle_sizePostResult(maxval(phase_Noutput),maxNinstance), source=0_pInt)
- allocate(source_damage_anisoBrittle_output(maxval(phase_Noutput),maxNinstance))
-          source_damage_anisoBrittle_output = ''
- allocate(source_damage_anisoBrittle_outputID(maxval(phase_Noutput),maxNinstance),       source=undefined_ID)
- allocate(source_damage_anisoBrittle_Noutput(maxNinstance),                              source=0_pInt) 
- allocate(source_damage_anisoBrittle_critDisp(lattice_maxNcleavageFamily,maxNinstance),  source=0.0_pReal) 
- allocate(source_damage_anisoBrittle_critLoad(lattice_maxNcleavageFamily,maxNinstance),  source=0.0_pReal) 
- allocate(source_damage_anisoBrittle_Ncleavage(lattice_maxNcleavageFamily,maxNinstance), source=0_pInt)
- allocate(source_damage_anisoBrittle_totalNcleavage(maxNinstance),                       source=0_pInt)
- allocate(source_damage_anisoBrittle_aTol(maxNinstance),                                 source=0.0_pReal) 
- allocate(source_damage_anisoBrittle_sdot_0(maxNinstance),                               source=0.0_pReal) 
- allocate(source_damage_anisoBrittle_N(maxNinstance),                                    source=0.0_pReal) 
-
- rewind(fileUnit)
- phase = 0_pInt
- do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase)         ! wind forward to <phase>
-   line = IO_read(fileUnit)
- enddo
  
- parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
+ allocate(source_damage_anisoBrittle_sizePostResult(maxval(phase_Noutput),Ninstance), source=0_pInt)
-   line = IO_read(fileUnit)
+ allocate(source_damage_anisoBrittle_output(maxval(phase_Noutput),Ninstance))
-   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
+          source_damage_anisoBrittle_output = ''
-   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
-     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
-     exit                                                                                           
-   endif   
-   if (IO_getTag(line,'[',']') /= '') then                                                          ! next phase section
-     phase = phase + 1_pInt                                                                         ! advance phase section counter
-     cycle                                                                                          ! skip to next line
-   endif
-   if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_damage_anisoBrittle_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
-     instance = source_damage_anisoBrittle_instance(phase)                                          ! which instance of my damage is present phase
-     chunkPos = IO_stringPos(line)
-     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
-     select case(tag)
-       case ('(output)')
-         select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
-           case ('anisobrittle_drivingforce')
-             source_damage_anisoBrittle_Noutput(instance) = source_damage_anisoBrittle_Noutput(instance) + 1_pInt
-             source_damage_anisoBrittle_outputID(source_damage_anisoBrittle_Noutput(instance),instance) = damage_drivingforce_ID
-             source_damage_anisoBrittle_output(source_damage_anisoBrittle_Noutput(instance),instance) = &
-                                                       IO_lc(IO_stringValue(line,chunkPos,2_pInt))
-          end select
 
-       case ('anisobrittle_atol')
+ allocate(source_damage_anisoBrittle_Ncleavage(lattice_maxNcleavageFamily,Ninstance), source=0_pInt)
-         source_damage_anisoBrittle_aTol(instance) = IO_floatValue(line,chunkPos,2_pInt)
-         
-       case ('anisobrittle_sdot0')
-         source_damage_anisoBrittle_sdot_0(instance) = IO_floatValue(line,chunkPos,2_pInt)
-         
-       case ('anisobrittle_ratesensitivity')
-         source_damage_anisoBrittle_N(instance) = IO_floatValue(line,chunkPos,2_pInt)
-         
-       case ('ncleavage')  !
-         Nchunks_CleavageFamilies = chunkPos(1) - 1_pInt
-         do j = 1_pInt, Nchunks_CleavageFamilies
-           source_damage_anisoBrittle_Ncleavage(j,instance) = IO_intValue(line,chunkPos,1_pInt+j)
-         enddo
 
-       case ('anisobrittle_criticaldisplacement')
+ allocate(param(Ninstance))
-         do j = 1_pInt, Nchunks_CleavageFamilies
+ 
-           source_damage_anisoBrittle_critDisp(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j)
+ do p=1, size(config_phase)
-         enddo
+   if (all(phase_source(:,p) /= SOURCE_DAMAGE_ANISOBRITTLE_ID)) cycle
+   associate(prm => param(source_damage_anisoBrittle_instance(p)), &
+             config => config_phase(p))
+             
+   prm%aTol      = config%getFloat('anisobrittle_atol',defaultVal = 1.0e-3_pReal)
 
-       case ('anisobrittle_criticalload')
+   prm%N         = config%getFloat('anisobrittle_ratesensitivity')
-         do j = 1_pInt, Nchunks_CleavageFamilies
+   prm%sdot_0    = config%getFloat('anisobrittle_sdot0')
-           source_damage_anisoBrittle_critLoad(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j)
+   
-         enddo
+   ! sanity checks
+   if (prm%aTol      < 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_atol'
+   
+   if (prm%N        <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_ratesensitivity'
+   if (prm%sdot_0   <= 0.0_pReal) extmsg = trim(extmsg)//' anisobrittle_sdot0'
+   
+   prm%Ncleavage = config%getInts('ncleavage',defaultVal=emptyIntArray)
+
+   prm%critDisp = config%getFloats('anisobrittle_criticaldisplacement',requiredSize=size(prm%Ncleavage))
+   prm%critLoad = config%getFloats('anisobrittle_criticalload',        requiredSize=size(prm%Ncleavage))
+
+     ! expand: family => system
+     prm%critDisp  = math_expand(prm%critDisp, prm%Ncleavage)
+     prm%critLoad  = math_expand(prm%critLoad, prm%Ncleavage)
+     
+     if (any(prm%critLoad < 0.0_pReal))     extmsg = trim(extmsg)//' anisobrittle_criticalload'
+    if (any(prm%critDisp < 0.0_pReal))     extmsg = trim(extmsg)//' anisobrittle_criticaldisplacement'  
+!--------------------------------------------------------------------------------------------------
+!  exit if any parameter is out of range
+   if (extmsg /= '') &
+     call IO_error(211_pInt,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISOBRITTLE_LABEL//')')
+
+!--------------------------------------------------------------------------------------------------
+!  output pararameters
+   outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
+   allocate(prm%outputID(0))
+   do i=1_pInt, size(outputs)
+     outputID = undefined_ID
+     select case(outputs(i))
+     
+       case ('anisobrittle_drivingforce')
+         source_damage_anisoBrittle_sizePostResult(i,source_damage_anisoBrittle_instance(p)) = 1_pInt
+         source_damage_anisoBrittle_output(i,source_damage_anisoBrittle_instance(p)) = outputs(i)
+         prm%outputID = [prm%outputID, damage_drivingforce_ID]
 
      end select
-   endif; endif
- enddo parsingFile
 
-!--------------------------------------------------------------------------------------------------
+   enddo
-!  sanity checks
+
- sanityChecks: do phase = 1_pInt, material_Nphase  
+   end associate
-   myPhase: if (any(phase_source(:,phase) == SOURCE_damage_anisoBrittle_ID)) then
+   
-     instance = source_damage_anisoBrittle_instance(phase)
+   phase = p
-     source_damage_anisoBrittle_Ncleavage(1:lattice_maxNcleavageFamily,instance) = &
+   NofMyPhase=count(material_phase==phase)
-       min(lattice_NcleavageSystem(1:lattice_maxNcleavageFamily,phase),&                            ! limit active cleavage systems per family to min of available and requested
+   instance = source_damage_anisoBrittle_instance(phase)
-           source_damage_anisoBrittle_Ncleavage(1:lattice_maxNcleavageFamily,instance))
+   sourceOffset = source_damage_anisoBrittle_offset(phase)
-     source_damage_anisoBrittle_totalNcleavage(instance)  = sum(source_damage_anisoBrittle_Ncleavage(:,instance)) ! how many cleavage systems altogether
+
-     if (source_damage_anisoBrittle_aTol(instance) < 0.0_pReal) &
+
-       source_damage_anisoBrittle_aTol(instance) = 1.0e-3_pReal                                     ! default absolute tolerance 1e-3
+   call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1_pInt,1_pInt,0_pInt)
-     if (source_damage_anisoBrittle_sdot_0(instance) <= 0.0_pReal) &
+   sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_anisoBrittle_sizePostResult(:,instance))
-       call IO_error(211_pInt,el=instance,ext_msg='sdot_0 ('//SOURCE_damage_anisoBrittle_LABEL//')')
+   sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol
-     if (any(source_damage_anisoBrittle_critDisp(1:Nchunks_CleavageFamilies,instance) < 0.0_pReal)) &
+
-       call IO_error(211_pInt,el=instance,ext_msg='critical_displacement ('//SOURCE_damage_anisoBrittle_LABEL//')')
+
-     if (any(source_damage_anisoBrittle_critLoad(1:Nchunks_CleavageFamilies,instance) < 0.0_pReal)) &
+   source_damage_anisoBrittle_Ncleavage(1:size(param(instance)%Ncleavage),instance) = param(instance)%Ncleavage
-       call IO_error(211_pInt,el=instance,ext_msg='critical_load ('//SOURCE_damage_anisoBrittle_LABEL//')')
+ enddo
-     if (source_damage_anisoBrittle_N(instance) <= 0.0_pReal) &
+
-       call IO_error(211_pInt,el=instance,ext_msg='rate_sensitivity ('//SOURCE_damage_anisoBrittle_LABEL//')')
-   endif myPhase
- enddo sanityChecks
  
- initializeInstances: do phase = 1_pInt, material_Nphase
-   if (any(phase_source(:,phase) == SOURCE_damage_anisoBrittle_ID)) then
-     NofMyPhase=count(material_phase==phase)
-     instance = source_damage_anisoBrittle_instance(phase)
-     sourceOffset = source_damage_anisoBrittle_offset(phase)
-
-!--------------------------------------------------------------------------------------------------
-!  Determine size of postResults array
-     outputsLoop: do o = 1_pInt,source_damage_anisoBrittle_Noutput(instance)
-       select case(source_damage_anisoBrittle_outputID(o,instance))
-         case(damage_drivingforce_ID)
-           mySize = 1_pInt
-       end select
- 
-       if (mySize > 0_pInt) then  ! any meaningful output found
-          source_damage_anisoBrittle_sizePostResult(o,instance) = mySize
-          source_damage_anisoBrittle_sizePostResults(instance)  = source_damage_anisoBrittle_sizePostResults(instance) + mySize
-       endif
-     enddo outputsLoop
-
-!--------------------------------------------------------------------------------------------------
-! Determine size of state array
-     sizeDotState              = 1_pInt
-     sizeDeltaState            = 0_pInt
-     sizeState                 = 1_pInt
-
-     sourceState(phase)%p(sourceOffset)%sizeState      = sizeState
-     sourceState(phase)%p(sourceOffset)%sizeDotState   = sizeDotState
-     sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState
-     sourceState(phase)%p(sourceOffset)%sizePostResults = source_damage_anisoBrittle_sizePostResults(instance)
-     allocate(sourceState(phase)%p(sourceOffset)%aTolState           (sizeState),                &
-              source=source_damage_anisoBrittle_aTol(instance))
-     allocate(sourceState(phase)%p(sourceOffset)%state0              (sizeState,NofMyPhase),     source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%partionedState0     (sizeState,NofMyPhase),     source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%subState0           (sizeState,NofMyPhase),     source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%state               (sizeState,NofMyPhase),     source=0.0_pReal)
-
-     allocate(sourceState(phase)%p(sourceOffset)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
-     if (any(numerics_integrator == 1_pInt)) then
-       allocate(sourceState(phase)%p(sourceOffset)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)
-       allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     endif
-     if (any(numerics_integrator == 4_pInt)) &
-       allocate(sourceState(phase)%p(sourceOffset)%RK4dotState       (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     if (any(numerics_integrator == 5_pInt)) &
-       allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState    (6,sizeDotState,NofMyPhase),source=0.0_pReal)
-
-   endif
- 
- enddo initializeInstances
 end subroutine source_damage_anisoBrittle_init
 
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates derived quantities from state
 !--------------------------------------------------------------------------------------------------
-subroutine source_damage_anisoBrittle_dotState(Tstar_v, ipc, ip, el)
+subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el)
+ use math, only: &
+   math_mul33xx33
  use material, only: &
    phaseAt, phasememberAt, &
    sourceState, &
@@ -294,7 +216,7 @@ subroutine source_damage_anisoBrittle_dotState(Tstar_v, ipc, ip, el)
    damage, &
    damageMapping
  use lattice, only: &
-   lattice_Scleavage_v, &
+   lattice_Scleavage, &
    lattice_maxNcleavageFamily, &
    lattice_NcleavageSystem
 
@@ -303,8 +225,8 @@ subroutine source_damage_anisoBrittle_dotState(Tstar_v, ipc, ip, el)
    ipc, &                                                                                           !< component-ID of integration point
    ip, &                                                                                            !< integration point
    el                                                                                               !< element
- real(pReal),  intent(in), dimension(6) :: &
+ real(pReal),  intent(in), dimension(3,3) :: &
-   Tstar_v                                                                                          !< 2nd Piola Kirchhoff stress tensor (Mandel)
+   S
  integer(pInt) :: &
    phase, &
    constituent, &
@@ -312,7 +234,7 @@ subroutine source_damage_anisoBrittle_dotState(Tstar_v, ipc, ip, el)
    sourceOffset, &
    damageOffset, &
    homog, &
-   f, i, index_myFamily
+   f, i, index_myFamily, index
  real(pReal) :: &
    traction_d, traction_t, traction_n, traction_crit
 
@@ -324,23 +246,26 @@ subroutine source_damage_anisoBrittle_dotState(Tstar_v, ipc, ip, el)
  damageOffset = damageMapping(homog)%p(ip,el)
  
  sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal
+ 
+ index = 1_pInt
  do f = 1_pInt,lattice_maxNcleavageFamily
-   index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase))                                    ! at which index starts my family
+   index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase))                                  ! at which index starts my family
-   do i = 1_pInt,source_damage_anisoBrittle_Ncleavage(f,instance)                                     ! process each (active) cleavage system in family
+   do i = 1_pInt,source_damage_anisoBrittle_Ncleavage(f,instance)                                   ! process each (active) cleavage system in family
-     traction_d    = dot_product(Tstar_v,lattice_Scleavage_v(1:6,1,index_myFamily+i,phase))
+     traction_d    = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase))
-     traction_t    = dot_product(Tstar_v,lattice_Scleavage_v(1:6,2,index_myFamily+i,phase))
+     traction_t    = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase))
-     traction_n    = dot_product(Tstar_v,lattice_Scleavage_v(1:6,3,index_myFamily+i,phase))
+     traction_n    = math_mul33xx33(S,lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase))
      
-     traction_crit = source_damage_anisoBrittle_critLoad(f,instance)* &
+     traction_crit = param(instance)%critLoad(index)* &
                      damage(homog)%p(damageOffset)*damage(homog)%p(damageOffset)
      sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = &
        sourceState(phase)%p(sourceOffset)%dotState(1,constituent) + &
-       source_damage_anisoBrittle_sdot_0(instance)* &
+       param(instance)%sdot_0* &
-       ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**source_damage_anisoBrittle_N(instance) + &
+       ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**param(instance)%N + &
-        (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**source_damage_anisoBrittle_N(instance) + &
+        (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**param(instance)%N + &
-        (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**source_damage_anisoBrittle_N(instance))/ &
+        (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**param(instance)%N)/ &
-       source_damage_anisoBrittle_critDisp(f,instance)
+       param(instance)%critDisp(index)
 
+   index = index + 1_pInt
    enddo
  enddo
 
@@ -349,30 +274,26 @@ end subroutine source_damage_anisoBrittle_dotState
 !--------------------------------------------------------------------------------------------------
 !> @brief returns local part of nonlocal damage driving force
 !--------------------------------------------------------------------------------------------------
-subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ipc, ip,  el)
+subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
  use material, only: &
-   phaseAt, phasememberAt, &
    sourceState
 
  implicit none
  integer(pInt), intent(in) :: &
-   ipc, &                                                                                           !< component-ID of integration point
+   phase, &
-   ip, &                                                                                            !< integration point
+   constituent
-   el                                                                                               !< element
  real(pReal),  intent(in) :: &
    phi
  real(pReal),  intent(out) :: &
    localphiDot, &
    dLocalphiDot_dPhi
  integer(pInt) :: &
-   phase, constituent, sourceOffset
+   sourceOffset
 
- phase = phaseAt(ipc,ip,el)
- constituent = phasememberAt(ipc,ip,el)
  sourceOffset = source_damage_anisoBrittle_offset(phase)
  
- localphiDot = 1.0_pReal - &
+ localphiDot = 1.0_pReal &
-               sourceState(phase)%p(sourceOffset)%state(1,constituent)*phi
+             - sourceState(phase)%p(sourceOffset)%state(1,constituent)*phi
  
  dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent)
  
@@ -381,33 +302,28 @@ end subroutine source_damage_anisobrittle_getRateAndItsTangent
 !--------------------------------------------------------------------------------------------------
 !> @brief return array of local damage results
 !--------------------------------------------------------------------------------------------------
-function source_damage_anisoBrittle_postResults(ipc,ip,el)
+function source_damage_anisoBrittle_postResults(phase, constituent)
  use material, only: &
-   phaseAt, phasememberAt, &
    sourceState
 
  implicit none
- integer(pInt),              intent(in) :: &
+ integer(pInt), intent(in) :: &
-   ipc, &                                                                                           !< component-ID of integration point
+   phase, &
-   ip, &                                                                                            !< integration point
+   constituent
-   el                                                                                               !< element
+ real(pReal), dimension(sum(source_damage_anisoBrittle_sizePostResult(:, &
- real(pReal), dimension(source_damage_anisoBrittle_sizePostResults( &
+                          source_damage_anisoBrittle_instance(phase)))) :: &
-                          source_damage_anisoBrittle_instance(phaseAt(ipc,ip,el)))) :: &
    source_damage_anisoBrittle_postResults
 
  integer(pInt) :: &
-   instance, phase, constituent, sourceOffset, o, c
+   instance, sourceOffset, o, c
    
- phase = phaseAt(ipc,ip,el)
- constituent = phasememberAt(ipc,ip,el)
  instance = source_damage_anisoBrittle_instance(phase)
  sourceOffset = source_damage_anisoBrittle_offset(phase)
 
  c = 0_pInt
- source_damage_anisoBrittle_postResults = 0.0_pReal
 
- do o = 1_pInt,source_damage_anisoBrittle_Noutput(instance)
+ do o = 1_pInt,size(param(instance)%outputID)
-    select case(source_damage_anisoBrittle_outputID(o,instance))
+    select case(param(instance)%outputID(o))
       case (damage_drivingforce_ID)
         source_damage_anisoBrittle_postResults(c+1_pInt) = &
           sourceState(phase)%p(sourceOffset)%state(1,constituent)

src/source_damage_anisoDuctile.f90

@@ -12,7 +12,6 @@ module source_damage_anisoDuctile
  implicit none
  private
  integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
-   source_damage_anisoDuctile_sizePostResults, &                                                                !< cumulative size of post results
    source_damage_anisoDuctile_offset, &                                                                         !< which source is my current damage mechanism?
    source_damage_anisoDuctile_instance                                                                          !< instance of damage source mechanism
 
@@ -22,35 +21,31 @@ module source_damage_anisoDuctile
  character(len=64),                   dimension(:,:),         allocatable, target, public  :: &
    source_damage_anisoDuctile_output                                                                            !< name of each post result output
    
- integer(pInt),                       dimension(:),           allocatable, target, public  :: &
-   source_damage_anisoDuctile_Noutput                                                                           !< number of outputs per instance of this damage 
-   
- integer(pInt),                       dimension(:),           allocatable,         private :: &
-   source_damage_anisoDuctile_totalNslip                                                                    !< total number of slip systems
 
  integer(pInt),                       dimension(:,:),         allocatable,         private :: &
    source_damage_anisoDuctile_Nslip                                                                         !< number of slip systems per family
    
- real(pReal),                         dimension(:),           allocatable,         private :: &
-   source_damage_anisoDuctile_aTol
-
- real(pReal),                         dimension(:,:),         allocatable,         private :: &
-   source_damage_anisoDuctile_critPlasticStrain
-
- real(pReal),                         dimension(:),           allocatable,         private :: &
-   source_damage_anisoDuctile_sdot_0, &
-   source_damage_anisoDuctile_N
-
- real(pReal),                         dimension(:,:),         allocatable,         private :: &
-   source_damage_anisoDuctile_critLoad
-   
  enum, bind(c) 
    enumerator :: undefined_ID, &
                  damage_drivingforce_ID
  end enum 
- 
+
- integer(kind(undefined_ID)),         dimension(:,:),         allocatable,          private :: & 
+
-   source_damage_anisoDuctile_outputID                                                                  !< ID of each post result output
+ type, private :: tParameters                                                                       !< container type for internal constitutive parameters
+   real(pReal) :: &
+     aTol, &
+     N
+   real(pReal), dimension(:), allocatable :: &
+     critPlasticStrain
+   integer(pInt) :: &
+     totalNslip
+   integer(pInt), dimension(:), allocatable :: &
+     Nslip
+   integer(kind(undefined_ID)), allocatable, dimension(:) :: &
+     outputID
+ end type tParameters
+
+ type(tParameters), dimension(:), allocatable, private :: param                                     !< containers of constitutive parameters (len Ninstance)
 
 
  public :: &
@@ -66,30 +61,24 @@ contains
 !> @brief module initialization
 !> @details reads in material parameters, allocates arrays, and does sanity checks
 !--------------------------------------------------------------------------------------------------
-subroutine source_damage_anisoDuctile_init(fileUnit)
+subroutine source_damage_anisoDuctile_init
 #if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
  use, intrinsic :: iso_fortran_env, only: &
    compiler_version, &
    compiler_options
 #endif
+ use prec, only: &
+   pStringLen
  use debug, only: &
    debug_level,&
    debug_constitutive,&
    debug_levelBasic
  use IO, only: &
-   IO_read, &
+   IO_error
-   IO_lc, &
+    use math, only: &
-   IO_getTag, &
+   math_expand
-   IO_isBlank, &
-   IO_stringPos, &
-   IO_stringValue, &
-   IO_floatValue, &
-   IO_intValue, &
-   IO_warning, &
-   IO_error, &
-   IO_timeStamp, &
-   IO_EOF
  use material, only: &
+   material_allocateSourceState, &
    phase_source, &
    phase_Nsources, &
    phase_Noutput, &
@@ -98,35 +87,35 @@ subroutine source_damage_anisoDuctile_init(fileUnit)
    material_phase, &  
    sourceState
  use config, only: &
+   config_phase, &
    material_Nphase, &
    MATERIAL_partPhase
- use numerics,only: &
-   numerics_integrator
  use lattice, only: &
-   lattice_maxNslipFamily, &
+   lattice_maxNslipFamily
-   lattice_NslipSystem
+   
-
  implicit none
- integer(pInt), intent(in) :: fileUnit
 
- integer(pInt), allocatable, dimension(:) :: chunkPos
+ integer(pInt) :: Ninstance,phase,instance,source,sourceOffset
- integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,source,sourceOffset,o
+ integer(pInt) :: NofMyPhase,p ,i
- integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
- integer(pInt) :: NofMyPhase   
- integer(pInt) :: Nchunks_SlipFamilies = 0_pInt, j   
- character(len=65536) :: &
-   tag  = '', &
-   line = ''
 
- write(6,'(/,a)')   ' <<<+-  source_'//SOURCE_damage_anisoDuctile_LABEL//' init  -+>>>'
+ integer(pInt),          dimension(0), parameter :: emptyIntArray    = [integer(pInt)::]
- write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
+ character(len=65536),   dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
+ integer(kind(undefined_ID)) :: &
+   outputID
+
+ character(len=pStringLen) :: &
+   extmsg = ''
+ character(len=65536), dimension(:), allocatable :: &
+   outputs
+
+ write(6,'(/,a)')   ' <<<+-  source_'//SOURCE_DAMAGE_ANISODUCTILE_LABEL//' init  -+>>>'
 #include "compilation_info.f90"
 
- maxNinstance = int(count(phase_source == SOURCE_damage_anisoDuctile_ID),pInt)
+ Ninstance = int(count(phase_source == SOURCE_damage_anisoDuctile_ID),pInt)
- if (maxNinstance == 0_pInt) return
+ if (Ninstance == 0_pInt) return
  
  if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
-   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
+   write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
  
  allocate(source_damage_anisoDuctile_offset(material_Nphase), source=0_pInt)
  allocate(source_damage_anisoDuctile_instance(material_Nphase), source=0_pInt)
@@ -138,151 +127,75 @@ subroutine source_damage_anisoDuctile_init(fileUnit)
    enddo    
  enddo
    
- allocate(source_damage_anisoDuctile_sizePostResults(maxNinstance),                     source=0_pInt)
+ allocate(source_damage_anisoDuctile_sizePostResult(maxval(phase_Noutput),Ninstance),source=0_pInt)
- allocate(source_damage_anisoDuctile_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
+ allocate(source_damage_anisoDuctile_output(maxval(phase_Noutput),Ninstance))
- allocate(source_damage_anisoDuctile_output(maxval(phase_Noutput),maxNinstance))
           source_damage_anisoDuctile_output = ''
- allocate(source_damage_anisoDuctile_outputID(maxval(phase_Noutput),maxNinstance),      source=undefined_ID)
- allocate(source_damage_anisoDuctile_Noutput(maxNinstance),                             source=0_pInt) 
- allocate(source_damage_anisoDuctile_critLoad(lattice_maxNslipFamily,maxNinstance), source=0.0_pReal) 
- allocate(source_damage_anisoDuctile_critPlasticStrain(lattice_maxNslipFamily,maxNinstance),source=0.0_pReal) 
- allocate(source_damage_anisoDuctile_Nslip(lattice_maxNslipFamily,maxNinstance),        source=0_pInt)
- allocate(source_damage_anisoDuctile_totalNslip(maxNinstance),                          source=0_pInt)
- allocate(source_damage_anisoDuctile_N(maxNinstance),                                   source=0.0_pReal) 
- allocate(source_damage_anisoDuctile_sdot_0(maxNinstance),                              source=0.0_pReal) 
- allocate(source_damage_anisoDuctile_aTol(maxNinstance),                                source=0.0_pReal) 
 
- rewind(fileUnit)
+ allocate(source_damage_anisoDuctile_Nslip(lattice_maxNslipFamily,Ninstance),        source=0_pInt)
- phase = 0_pInt
+
- do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase)         ! wind forward to <phase>
+ allocate(param(Ninstance))
-   line = IO_read(fileUnit)
- enddo
  
- parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
+ do p=1, size(config_phase)
-   line = IO_read(fileUnit)
+   if (all(phase_source(:,p) /= SOURCE_DAMAGE_ANISODUCTILE_ID)) cycle
-   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
+   associate(prm => param(source_damage_anisoDuctile_instance(p)), &
-   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
+             config => config_phase(p))
-     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
+             
-     exit                                                                                           
+   prm%aTol   = config%getFloat('anisoductile_atol',defaultVal = 1.0e-3_pReal)
-   endif   
-   if (IO_getTag(line,'[',']') /= '') then                                                          ! next phase section
-     phase = phase + 1_pInt                                                                         ! advance phase section counter
-     cycle                                                                                          ! skip to next line
-   endif
-   if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_damage_anisoDuctile_ID)) then         ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
-     instance = source_damage_anisoDuctile_instance(phase)                                                         ! which instance of my damage is present phase
-     chunkPos = IO_stringPos(line)
-     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
-     select case(tag)
-       case ('(output)')
-         select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
-           case ('anisoductile_drivingforce')
-             source_damage_anisoDuctile_Noutput(instance) = source_damage_anisoDuctile_Noutput(instance) + 1_pInt
-             source_damage_anisoDuctile_outputID(source_damage_anisoDuctile_Noutput(instance),instance) = damage_drivingforce_ID
-             source_damage_anisoDuctile_output(source_damage_anisoDuctile_Noutput(instance),instance) = &
-                                                       IO_lc(IO_stringValue(line,chunkPos,2_pInt))
-          end select
 
-       case ('anisoductile_atol')
+   prm%N      = config%getFloat('anisoductile_ratesensitivity')
-         source_damage_anisoDuctile_aTol(instance) = IO_floatValue(line,chunkPos,2_pInt)
+   
-         
+   ! sanity checks
-       case ('nslip')  !
+   if (prm%aTol                 < 0.0_pReal) extmsg = trim(extmsg)//' anisoductile_atol'
-         Nchunks_SlipFamilies = chunkPos(1) - 1_pInt
+   
-         do j = 1_pInt, Nchunks_SlipFamilies
+   if (prm%N                   <= 0.0_pReal) extmsg = trim(extmsg)//' anisoductile_ratesensitivity'
-           source_damage_anisoDuctile_Nslip(j,instance) = IO_intValue(line,chunkPos,1_pInt+j)
+   
-         enddo
+   prm%Nslip  = config%getInts('nslip',defaultVal=emptyIntArray)
+   
+   prm%critPlasticStrain = config%getFloats('anisoductile_criticalplasticstrain',requiredSize=size(prm%Nslip))
 
-       case ('anisoductile_sdot0')
+     ! expand: family => system
-         source_damage_anisoDuctile_sdot_0(instance) = IO_floatValue(line,chunkPos,2_pInt)
+     prm%critPlasticStrain   = math_expand(prm%critPlasticStrain,  prm%Nslip)
-         
+     
-       case ('anisoductile_criticalplasticstrain')
+          if (any(prm%critPlasticStrain < 0.0_pReal))     extmsg = trim(extmsg)//' anisoductile_criticalplasticstrain'
-         do j = 1_pInt, Nchunks_SlipFamilies
+   
-           source_damage_anisoDuctile_critPlasticStrain(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j)
+!--------------------------------------------------------------------------------------------------
-         enddo
+!  exit if any parameter is out of range
-         
+   if (extmsg /= '') &
-       case ('anisoductile_ratesensitivity')
+     call IO_error(211_pInt,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ANISODUCTILE_LABEL//')')
-         source_damage_anisoDuctile_N(instance) = IO_floatValue(line,chunkPos,2_pInt)
+
+!--------------------------------------------------------------------------------------------------
+!  output pararameters
+   outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
+   allocate(prm%outputID(0))
+   do i=1_pInt, size(outputs)
+     outputID = undefined_ID
+     select case(outputs(i))
+     
+       case ('anisoductile_drivingforce')
+         source_damage_anisoDuctile_sizePostResult(i,source_damage_anisoDuctile_instance(p)) = 1_pInt
+         source_damage_anisoDuctile_output(i,source_damage_anisoDuctile_instance(p)) = outputs(i)
+         prm%outputID = [prm%outputID, damage_drivingforce_ID]
 
-       case ('anisoductile_criticalload')
-         do j = 1_pInt, Nchunks_SlipFamilies
-           source_damage_anisoDuctile_critLoad(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j)
-         enddo
-         
      end select
-   endif; endif
- enddo parsingFile
 
-!--------------------------------------------------------------------------------------------------
+   enddo
-!  sanity checks
+
- sanityChecks: do phase = 1_pInt, size(phase_source)   
+   end associate
-   myPhase: if (any(phase_source(:,phase) == SOURCE_damage_anisoDuctile_ID)) then
+   
-     instance = source_damage_anisoDuctile_instance(phase)
+   phase = p
-     source_damage_anisoDuctile_Nslip(1:lattice_maxNslipFamily,instance) = &
+   
-       min(lattice_NslipSystem(1:lattice_maxNslipFamily,phase),&                                    ! limit active cleavage systems per family to min of available and requested
+   NofMyPhase=count(material_phase==phase)
-           source_damage_anisoDuctile_Nslip(1:lattice_maxNslipFamily,instance))
+   instance = source_damage_anisoDuctile_instance(phase)
-         source_damage_anisoDuctile_totalNslip(instance) = sum(source_damage_anisoDuctile_Nslip(:,instance))
+   sourceOffset = source_damage_anisoDuctile_offset(phase)
-     if (source_damage_anisoDuctile_aTol(instance) < 0.0_pReal) &
+
-       source_damage_anisoDuctile_aTol(instance) = 1.0e-3_pReal                                     ! default absolute tolerance 1e-3
+   call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1_pInt,1_pInt,0_pInt)
-     if (source_damage_anisoDuctile_sdot_0(instance) <= 0.0_pReal) &
+   sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_anisoDuctile_sizePostResult(:,instance))
-       call IO_error(211_pInt,el=instance,ext_msg='sdot_0 ('//SOURCE_damage_anisoDuctile_LABEL//')')
+   sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol
-     if (any(source_damage_anisoDuctile_critPlasticStrain(:,instance) < 0.0_pReal)) &
+   
-       call IO_error(211_pInt,el=instance,ext_msg='criticaPlasticStrain ('//SOURCE_damage_anisoDuctile_LABEL//')')
+   source_damage_anisoDuctile_Nslip(1:size(param(instance)%Nslip),instance) = param(instance)%Nslip
-     if (source_damage_anisoDuctile_N(instance) <= 0.0_pReal) &
+   
-       call IO_error(211_pInt,el=instance,ext_msg='rate_sensitivity ('//SOURCE_damage_anisoDuctile_LABEL//')')
+ enddo
-   endif myPhase
- enddo sanityChecks
   
-
- initializeInstances: do phase = 1_pInt, material_Nphase
-   if (any(phase_source(:,phase) == SOURCE_damage_anisoDuctile_ID)) then
-     NofMyPhase=count(material_phase==phase)
-     instance = source_damage_anisoDuctile_instance(phase)
-     sourceOffset = source_damage_anisoDuctile_offset(phase)
-
-!--------------------------------------------------------------------------------------------------
-!  Determine size of postResults array
-     outputsLoop: do o = 1_pInt,source_damage_anisoDuctile_Noutput(instance)
-       select case(source_damage_anisoDuctile_outputID(o,instance))
-         case(damage_drivingforce_ID)
-           mySize = 1_pInt
-       end select
- 
-       if (mySize > 0_pInt) then  ! any meaningful output found
-          source_damage_anisoDuctile_sizePostResult(o,instance) = mySize
-          source_damage_anisoDuctile_sizePostResults(instance)  = source_damage_anisoDuctile_sizePostResults(instance) + mySize
-       endif
-     enddo outputsLoop
-
-!--------------------------------------------------------------------------------------------------
-! Determine size of state array
-     sizeDotState              = 1_pInt
-     sizeDeltaState            = 0_pInt
-     sizeState                 = 1_pInt
-     sourceState(phase)%p(sourceOffset)%sizeState      = sizeState
-     sourceState(phase)%p(sourceOffset)%sizeDotState   = sizeDotState
-     sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState
-     sourceState(phase)%p(sourceOffset)%sizePostResults = source_damage_anisoDuctile_sizePostResults(instance)
-     allocate(sourceState(phase)%p(sourceOffset)%aTolState           (sizeState),                &
-              source=source_damage_anisoDuctile_aTol(instance))
-     allocate(sourceState(phase)%p(sourceOffset)%state0              (sizeState,NofMyPhase),     source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%partionedState0     (sizeState,NofMyPhase),     source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%subState0           (sizeState,NofMyPhase),     source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%state               (sizeState,NofMyPhase),     source=0.0_pReal)
-
-     allocate(sourceState(phase)%p(sourceOffset)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
-     if (any(numerics_integrator == 1_pInt)) then
-       allocate(sourceState(phase)%p(sourceOffset)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)
-       allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     endif
-     if (any(numerics_integrator == 4_pInt)) &
-       allocate(sourceState(phase)%p(sourceOffset)%RK4dotState       (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     if (any(numerics_integrator == 5_pInt)) &
-       allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState    (6,sizeDotState,NofMyPhase),source=0.0_pReal)
-
-   endif
- 
- enddo initializeInstances
 end subroutine source_damage_anisoDuctile_init
 
 !--------------------------------------------------------------------------------------------------
@@ -326,8 +239,7 @@ subroutine source_damage_anisoDuctile_dotState(ipc, ip, el)
      sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = &
        sourceState(phase)%p(sourceOffset)%dotState(1,constituent) + &
        plasticState(phase)%slipRate(index,constituent)/ &
-       ((damage(homog)%p(damageOffset))**source_damage_anisoDuctile_N(instance))/ & 
+       ((damage(homog)%p(damageOffset))**param(instance)%N)/param(instance)%critPlasticStrain(index) 
-       source_damage_anisoDuctile_critPlasticStrain(f,instance) 
      
      index = index + 1_pInt
    enddo
@@ -338,31 +250,26 @@ end subroutine source_damage_anisoDuctile_dotState
 !--------------------------------------------------------------------------------------------------
 !> @brief returns local part of nonlocal damage driving force
 !--------------------------------------------------------------------------------------------------
-subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ipc, ip,  el)
+subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
  use material, only: &
-   phaseAt, phasememberAt, &
    sourceState
 
  implicit none
  integer(pInt), intent(in) :: &
-   ipc, &                                                                                           !< component-ID of integration point
+   phase, &
-   ip, &                                                                                            !< integration point
+   constituent
-   el                                                                                               !< element
  real(pReal),  intent(in) :: &
    phi
  real(pReal),  intent(out) :: &
    localphiDot, &
    dLocalphiDot_dPhi
  integer(pInt) :: &
-   phase, constituent, sourceOffset
+   sourceOffset
 
- phase = phaseAt(ipc,ip,el)
- constituent = phasememberAt(ipc,ip,el)
  sourceOffset = source_damage_anisoDuctile_offset(phase)
  
- localphiDot = 1.0_pReal - &
+ localphiDot = 1.0_pReal &
-               sourceState(phase)%p(sourceOffset)%state(1,constituent)* &
+             - sourceState(phase)%p(sourceOffset)%state(1,constituent) * phi
-               phi
  
  dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent)
  
@@ -371,33 +278,28 @@ end subroutine source_damage_anisoDuctile_getRateAndItsTangent
 !--------------------------------------------------------------------------------------------------
 !> @brief return array of local damage results
 !--------------------------------------------------------------------------------------------------
-function source_damage_anisoDuctile_postResults(ipc,ip,el)
+function source_damage_anisoDuctile_postResults(phase, constituent)
  use material, only: &
-   phaseAt, phasememberAt, &
    sourceState
 
  implicit none
- integer(pInt),              intent(in) :: &
+ integer(pInt), intent(in) :: &
-   ipc, &                                                                                           !< component-ID of integration point
+   phase, &
-   ip, &                                                                                            !< integration point
+   constituent
-   el                                                                                               !< element
+ real(pReal), dimension(sum(source_damage_anisoDuctile_sizePostResult(:, &
- real(pReal), dimension(source_damage_anisoDuctile_sizePostResults( &
+                          source_damage_anisoDuctile_instance(phase)))) :: &
-                          source_damage_anisoDuctile_instance(phaseAt(ipc,ip,el)))) :: &
    source_damage_anisoDuctile_postResults
 
  integer(pInt) :: &
-   instance, phase, constituent, sourceOffset, o, c
+   instance, sourceOffset, o, c
    
- phase = phaseAt(ipc,ip,el)
- constituent = phasememberAt(ipc,ip,el)
  instance = source_damage_anisoDuctile_instance(phase)
  sourceOffset = source_damage_anisoDuctile_offset(phase)
 
  c = 0_pInt
- source_damage_anisoDuctile_postResults = 0.0_pReal
 
- do o = 1_pInt,source_damage_anisoDuctile_Noutput(instance)
+ do o = 1_pInt,size(param(instance)%outputID)
-    select case(source_damage_anisoDuctile_outputID(o,instance))
+    select case(param(instance)%outputID(o))
       case (damage_drivingforce_ID)
         source_damage_anisoDuctile_postResults(c+1_pInt) = &
           sourceState(phase)%p(sourceOffset)%state(1,constituent)

src/source_damage_isoBrittle.f90

@@ -12,7 +12,6 @@ module source_damage_isoBrittle
  implicit none
  private
  integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
-   source_damage_isoBrittle_sizePostResults, &                                                        !< cumulative size of post results
    source_damage_isoBrittle_offset, &                                                                 !< which source is my current damage mechanism?
    source_damage_isoBrittle_instance                                                                  !< instance of damage source mechanism
 
@@ -21,22 +20,23 @@ module source_damage_isoBrittle
 
  character(len=64),                   dimension(:,:),         allocatable, target, public :: &
    source_damage_isoBrittle_output                                                                    !< name of each post result output
-   
- integer(pInt),                       dimension(:),           allocatable, target, public :: &
-   source_damage_isoBrittle_Noutput                                                                   !< number of outputs per instance of this damage 
-
- real(pReal),                         dimension(:),     allocatable,         private :: &
-   source_damage_isoBrittle_aTol, &
-   source_damage_isoBrittle_N, &
-   source_damage_isoBrittle_critStrainEnergy
 
  enum, bind(c) 
    enumerator :: undefined_ID, &
                  damage_drivingforce_ID
  end enum                                                 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11 ToDo
- 
+
- integer(kind(undefined_ID)),         dimension(:,:),         allocatable,          private :: & 
+
-   source_damage_isoBrittle_outputID                                                                  !< ID of each post result output
+ type, private :: tParameters                                                                       !< container type for internal constitutive parameters
+   real(pReal) :: &
+     critStrainEnergy, &
+     N, &
+     aTol
+   integer(kind(undefined_ID)), allocatable, dimension(:) :: &
+     outputID
+ end type tParameters
+
+ type(tParameters), dimension(:), allocatable, private :: param                                     !< containers of constitutive parameters (len Ninstance)
 
 
  public :: &
@@ -52,30 +52,22 @@ contains
 !> @brief module initialization
 !> @details reads in material parameters, allocates arrays, and does sanity checks
 !--------------------------------------------------------------------------------------------------
-subroutine source_damage_isoBrittle_init(fileUnit)
+subroutine source_damage_isoBrittle_init
 #if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
  use, intrinsic :: iso_fortran_env, only: &
    compiler_version, &
    compiler_options
 #endif
+ use prec, only: &
+   pStringLen
  use debug, only: &
    debug_level,&
    debug_constitutive,&
    debug_levelBasic
  use IO, only: &
-   IO_read, &
+   IO_error
-   IO_lc, &
-   IO_getTag, &
-   IO_isBlank, &
-   IO_stringPos, &
-   IO_stringValue, &
-   IO_floatValue, &
-   IO_intValue, &
-   IO_warning, &
-   IO_error, &
-   IO_timeStamp, &
-   IO_EOF
  use material, only: &
+   material_allocateSourceState, &
    phase_source, &
    phase_Nsources, &
    phase_Noutput, &
@@ -84,31 +76,31 @@ subroutine source_damage_isoBrittle_init(fileUnit)
    material_phase, &  
    sourceState
  use config, only: &
+   config_phase, &
    material_Nphase, &
    MATERIAL_partPhase
- use numerics,only: &
-   numerics_integrator
 
  implicit none
- integer(pInt), intent(in) :: fileUnit
 
- integer(pInt), allocatable, dimension(:) :: chunkPos
+ integer(pInt) :: Ninstance,phase,instance,source,sourceOffset
- integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,source,sourceOffset,o
+ integer(pInt) :: NofMyPhase,p,i   
- integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
+ character(len=65536),   dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
- integer(pInt) :: NofMyPhase   
+ integer(kind(undefined_ID)) :: &
- character(len=65536) :: &
+   outputID
-   tag  = '', &
-   line = ''
 
- write(6,'(/,a)')   ' <<<+-  source_'//SOURCE_damage_isoBrittle_label//' init  -+>>>'
+ character(len=pStringLen) :: &
- write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
+   extmsg = ''
+ character(len=65536), dimension(:), allocatable :: &
+   outputs
+ 
+ write(6,'(/,a)')   ' <<<+-  source_'//SOURCE_DAMAGE_ISOBRITTLE_LABEL//' init  -+>>>'
 #include "compilation_info.f90"
 
- maxNinstance = int(count(phase_source == SOURCE_damage_isoBrittle_ID),pInt)
+ Ninstance = int(count(phase_source == SOURCE_damage_isoBrittle_ID),pInt)
- if (maxNinstance == 0_pInt) return
+ if (Ninstance == 0_pInt) return
  
  if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
-   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
+   write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
  
  allocate(source_damage_isoBrittle_offset(material_Nphase), source=0_pInt)
  allocate(source_damage_isoBrittle_instance(material_Nphase), source=0_pInt)
@@ -120,121 +112,64 @@ subroutine source_damage_isoBrittle_init(fileUnit)
    enddo    
  enddo
    
- allocate(source_damage_isoBrittle_sizePostResults(maxNinstance),                     source=0_pInt)
+ allocate(source_damage_isoBrittle_sizePostResult(maxval(phase_Noutput),Ninstance),source=0_pInt)
- allocate(source_damage_isoBrittle_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
+ allocate(source_damage_isoBrittle_output(maxval(phase_Noutput),Ninstance))
- allocate(source_damage_isoBrittle_output(maxval(phase_Noutput),maxNinstance))
           source_damage_isoBrittle_output = ''
- allocate(source_damage_isoBrittle_outputID(maxval(phase_Noutput),maxNinstance),      source=undefined_ID)
- allocate(source_damage_isoBrittle_Noutput(maxNinstance),                             source=0_pInt) 
- allocate(source_damage_isoBrittle_critStrainEnergy(maxNinstance),                    source=0.0_pReal) 
- allocate(source_damage_isoBrittle_N(maxNinstance),                                   source=1.0_pReal) 
- allocate(source_damage_isoBrittle_aTol(maxNinstance),                                source=0.0_pReal) 
 
- rewind(fileUnit)
+ allocate(param(Ninstance))
- phase = 0_pInt
+ 
- do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase)         ! wind forward to <phase>
+ do p=1, size(config_phase)
-   line = IO_read(fileUnit)
+   if (all(phase_source(:,p) /= SOURCE_DAMAGE_ISOBRITTLE_ID)) cycle
+   associate(prm => param(source_damage_isoBrittle_instance(p)), &
+             config => config_phase(p))
+             
+   prm%aTol             = config%getFloat('isobrittle_atol',defaultVal = 1.0e-3_pReal)
+   
+   prm%N                = config%getFloat('isobrittle_n')
+   prm%critStrainEnergy = config%getFloat('isobrittle_criticalstrainenergy')
+   
+   ! sanity checks
+   if (prm%aTol                < 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_atol'
+   
+   if (prm%N                  <= 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_n'
+   if (prm%critStrainEnergy   <= 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_criticalstrainenergy'
+   
+!--------------------------------------------------------------------------------------------------
+!  exit if any parameter is out of range
+   if (extmsg /= '') &
+     call IO_error(211_pInt,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ISOBRITTLE_LABEL//')')
+
+!--------------------------------------------------------------------------------------------------
+!  output pararameters
+   outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
+   allocate(prm%outputID(0))
+   do i=1_pInt, size(outputs)
+     outputID = undefined_ID
+     select case(outputs(i))
+     
+       case ('isobrittle_drivingforce')
+         source_damage_isoBrittle_sizePostResult(i,source_damage_isoBrittle_instance(p)) = 1_pInt
+         source_damage_isoBrittle_output(i,source_damage_isoBrittle_instance(p)) = outputs(i)
+         prm%outputID = [prm%outputID, damage_drivingforce_ID]
+           
+     end select
+
+   enddo
+
+   end associate
+   
+  phase = p
+  
+  NofMyPhase=count(material_phase==phase)
+  instance = source_damage_isoBrittle_instance(phase)
+  sourceOffset = source_damage_isoBrittle_offset(phase)
+     
+  call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1_pInt,1_pInt,1_pInt)
+  sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_isoBrittle_sizePostResult(:,instance))
+  sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol
+ 
  enddo
  
- parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
-   line = IO_read(fileUnit)
-   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
-   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
-     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
-     exit                                                                                           
-   endif   
-   if (IO_getTag(line,'[',']') /= '') then                                                          ! next phase section
-     phase = phase + 1_pInt                                                                         ! advance phase section counter
-     cycle                                                                                          ! skip to next line
-   endif
-   if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_damage_isoBrittle_ID)) then   ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
-     instance = source_damage_isoBrittle_instance(phase)                                            ! which instance of my damage is present phase
-     chunkPos = IO_stringPos(line)
-     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
-     select case(tag)
-       case ('(output)')
-         select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
-           case ('isobrittle_drivingforce')
-             source_damage_isoBrittle_Noutput(instance) = source_damage_isoBrittle_Noutput(instance) + 1_pInt
-             source_damage_isoBrittle_outputID(source_damage_isoBrittle_Noutput(instance),instance) = damage_drivingforce_ID
-             source_damage_isoBrittle_output(source_damage_isoBrittle_Noutput(instance),instance) = &
-                                                       IO_lc(IO_stringValue(line,chunkPos,2_pInt))
-          end select
-
-       case ('isobrittle_criticalstrainenergy')
-         source_damage_isoBrittle_critStrainEnergy(instance) = IO_floatValue(line,chunkPos,2_pInt)
-
-       case ('isobrittle_n')
-         source_damage_isoBrittle_N(instance) = IO_floatValue(line,chunkPos,2_pInt)
-
-       case ('isobrittle_atol')
-         source_damage_isoBrittle_aTol(instance) = IO_floatValue(line,chunkPos,2_pInt)
-
-     end select
-   endif; endif
- enddo parsingFile
-
-
-!--------------------------------------------------------------------------------------------------
-!  sanity checks
- sanityChecks: do phase = 1_pInt, material_Nphase 
-   myPhase: if (any(phase_source(:,phase) == SOURCE_damage_isoBrittle_ID)) then
-     instance = source_damage_isoBrittle_instance(phase)
-     if (source_damage_isoBrittle_aTol(instance) < 0.0_pReal) &
-       source_damage_isoBrittle_aTol(instance) = 1.0e-3_pReal                                              ! default absolute tolerance 1e-3
-     if (source_damage_isoBrittle_critStrainEnergy(instance) <= 0.0_pReal) &
-       call IO_error(211_pInt,el=instance,ext_msg='criticalStrainEnergy ('//SOURCE_damage_isoBrittle_LABEL//')')
-   endif myPhase
- enddo sanityChecks
-
- initializeInstances: do phase = 1_pInt, material_Nphase
-   if (any(phase_source(:,phase) == SOURCE_damage_isoBrittle_ID)) then
-     NofMyPhase=count(material_phase==phase)
-     instance = source_damage_isoBrittle_instance(phase)
-     sourceOffset = source_damage_isoBrittle_offset(phase)
-!--------------------------------------------------------------------------------------------------
-!  Determine size of postResults array
-     outputsLoop: do o = 1_pInt,source_damage_isoBrittle_Noutput(instance)
-       select case(source_damage_isoBrittle_outputID(o,instance))
-         case(damage_drivingforce_ID)
-           mySize = 1_pInt
-       end select
- 
-       if (mySize > 0_pInt) then  ! any meaningful output found
-          source_damage_isoBrittle_sizePostResult(o,instance) = mySize
-          source_damage_isoBrittle_sizePostResults(instance)  = source_damage_isoBrittle_sizePostResults(instance) + mySize
-       endif
-     enddo outputsLoop
-! Determine size of state array
-     sizeDotState              =   1_pInt
-     sizeDeltaState            =   1_pInt
-     sizeState                 =   1_pInt
-                
-     sourceState(phase)%p(sourceOffset)%sizeState      = sizeState
-     sourceState(phase)%p(sourceOffset)%sizeDotState   = sizeDotState
-     sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState
-     sourceState(phase)%p(sourceOffset)%sizePostResults = source_damage_isoBrittle_sizePostResults(instance)
-     allocate(sourceState(phase)%p(sourceOffset)%aTolState           (sizeState),                &
-              source=source_damage_isoBrittle_aTol(instance))
-     allocate(sourceState(phase)%p(sourceOffset)%state0              (sizeState,NofMyPhase),     source=.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%partionedState0     (sizeState,NofMyPhase),     source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%subState0           (sizeState,NofMyPhase),     source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%state               (sizeState,NofMyPhase),     source=0.0_pReal)
-
-     allocate(sourceState(phase)%p(sourceOffset)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
-     if (any(numerics_integrator == 1_pInt)) then
-       allocate(sourceState(phase)%p(sourceOffset)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)
-       allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     endif
-     if (any(numerics_integrator == 4_pInt)) &
-       allocate(sourceState(phase)%p(sourceOffset)%RK4dotState       (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     if (any(numerics_integrator == 5_pInt)) &
-       allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState    (6,sizeDotState,NofMyPhase),source=0.0_pReal)
-
-   endif
- 
- enddo initializeInstances
 end subroutine source_damage_isoBrittle_init
 
 !--------------------------------------------------------------------------------------------------
@@ -243,15 +178,11 @@ end subroutine source_damage_isoBrittle_init
 subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
  use material, only: &
    phaseAt, phasememberAt, &
-   sourceState, &
+   sourceState
-   material_homog, &
-   phase_NstiffnessDegradations, &
-   phase_stiffnessDegradation
  use math, only : &
+   math_sym33to6, &
    math_mul33x33, &
    math_mul66x6, &
-   math_Mandel33to6, &
-   math_transpose33, &
    math_I3
 
  implicit none
@@ -264,10 +195,9 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
  real(pReal),  intent(in), dimension(6,6) :: &
    C
  integer(pInt) :: &
-   phase, constituent, instance, sourceOffset, mech
+   phase, constituent, instance, sourceOffset
  real(pReal) :: &
    strain(6), &
-   stiffness(6,6), &
    strainenergy
 
  phase = phaseAt(ipc,ip,el)                                                            !< phase ID at ipc,ip,el
@@ -276,11 +206,11 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
  instance = source_damage_isoBrittle_instance(phase)                                                 !< instance of damage_isoBrittle source
  sourceOffset = source_damage_isoBrittle_offset(phase)
 
- stiffness = C                                            
+                                   
- strain = 0.5_pReal*math_Mandel33to6(math_mul33x33(math_transpose33(Fe),Fe)-math_I3)
+ strain = 0.5_pReal*math_sym33to6(math_mul33x33(transpose(Fe),Fe)-math_I3)
 
- strainenergy = 2.0_pReal*sum(strain*math_mul66x6(stiffness,strain))/ &
+ strainenergy = 2.0_pReal*sum(strain*math_mul66x6(C,strain))/param(instance)%critStrainEnergy
-                source_damage_isoBrittle_critStrainEnergy(instance) 
+ 
  if (strainenergy > sourceState(phase)%p(sourceOffset)%subState0(1,constituent)) then
    sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = &
      strainenergy - sourceState(phase)%p(sourceOffset)%state(1,constituent)
@@ -295,33 +225,29 @@ end subroutine source_damage_isoBrittle_deltaState
 !--------------------------------------------------------------------------------------------------
 !> @brief returns local part of nonlocal damage driving force
 !--------------------------------------------------------------------------------------------------
-subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ipc, ip, el)
+subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
  use material, only: &
-   phaseAt, phasememberAt, &
    sourceState
 
  implicit none
  integer(pInt), intent(in) :: &
-   ipc, &                                                                                           !< component-ID of integration point
+   phase, &
-   ip, &                                                                                            !< integration point
+   constituent
-   el                                                                                               !< element
  real(pReal),  intent(in) :: &
    phi
  real(pReal),  intent(out) :: &
    localphiDot, &
    dLocalphiDot_dPhi
  integer(pInt) :: &
-   phase, constituent, instance, sourceOffset
+   instance, sourceOffset
 
- phase = phaseAt(ipc,ip,el)
- constituent = phasememberAt(ipc,ip,el)
  instance = source_damage_isoBrittle_instance(phase)
  sourceOffset = source_damage_isoBrittle_offset(phase)
  
- localphiDot = (1.0_pReal - phi)**(source_damage_isoBrittle_N(instance) - 1.0_pReal) - &
+ localphiDot = (1.0_pReal - phi)**(param(instance)%N - 1.0_pReal) - &
                phi*sourceState(phase)%p(sourceOffset)%state(1,constituent)
- dLocalphiDot_dPhi = - (source_damage_isoBrittle_N(instance) - 1.0_pReal)* &
+ dLocalphiDot_dPhi = - (param(instance)%N - 1.0_pReal)* &
-                       (1.0_pReal - phi)**max(0.0_pReal,source_damage_isoBrittle_N(instance) - 2.0_pReal) &
+                       (1.0_pReal - phi)**max(0.0_pReal,param(instance)%N - 2.0_pReal) &
                      - sourceState(phase)%p(sourceOffset)%state(1,constituent)
  
 end subroutine source_damage_isoBrittle_getRateAndItsTangent
@@ -329,33 +255,28 @@ end subroutine source_damage_isoBrittle_getRateAndItsTangent
 !--------------------------------------------------------------------------------------------------
 !> @brief return array of local damage results
 !--------------------------------------------------------------------------------------------------
-function source_damage_isoBrittle_postResults(ipc,ip,el)
+function source_damage_isoBrittle_postResults(phase, constituent)
  use material, only: &
-   phaseAt, phasememberAt, &
    sourceState
 
  implicit none
- integer(pInt),              intent(in) :: &
+ integer(pInt), intent(in) :: &
-   ipc, &                                                                                           !< component-ID of integration point
+   phase, &
-   ip, &                                                                                            !< integration point
+   constituent
-   el                                                                                               !< element
+ real(pReal), dimension(sum(source_damage_isoBrittle_sizePostResult(:, &
- real(pReal), dimension(source_damage_isoBrittle_sizePostResults( &
+                         source_damage_isoBrittle_instance(phase)))) :: &
-                         source_damage_isoBrittle_instance(phaseAt(ipc,ip,el)))) :: &
    source_damage_isoBrittle_postResults
 
  integer(pInt) :: &
-   instance, phase, constituent, sourceOffset, o, c
+   instance, sourceOffset, o, c
    
- phase = phaseAt(ipc,ip,el)
- constituent = phasememberAt(ipc,ip,el)
  instance = source_damage_isoBrittle_instance(phase)
  sourceOffset = source_damage_isoBrittle_offset(phase)
 
  c = 0_pInt
- source_damage_isoBrittle_postResults = 0.0_pReal
 
- do o = 1_pInt,source_damage_isoBrittle_Noutput(instance)
+ do o = 1_pInt,size(param(instance)%outputID)
-    select case(source_damage_isoBrittle_outputID(o,instance))
+    select case(param(instance)%outputID(o))
       case (damage_drivingforce_ID)
         source_damage_isoBrittle_postResults(c+1_pInt) = sourceState(phase)%p(sourceOffset)%state(1,constituent)
         c = c + 1

src/source_damage_isoDuctile.f90

@@ -12,7 +12,6 @@ module source_damage_isoDuctile
  implicit none
  private
  integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
-   source_damage_isoDuctile_sizePostResults, &                                                        !< cumulative size of post results
    source_damage_isoDuctile_offset, &                                                                 !< which source is my current damage mechanism?
    source_damage_isoDuctile_instance                                                                  !< instance of damage source mechanism
 
@@ -21,22 +20,23 @@ module source_damage_isoDuctile
 
  character(len=64),                   dimension(:,:),         allocatable, target, public :: &
    source_damage_isoDuctile_output                                                                    !< name of each post result output
-   
- integer(pInt),                       dimension(:),           allocatable, target, public :: &
-   source_damage_isoDuctile_Noutput                                                                   !< number of outputs per instance of this damage 
 
- real(pReal),                         dimension(:),     allocatable,         private :: &
-   source_damage_isoDuctile_aTol, &
-   source_damage_isoDuctile_critPlasticStrain, &
-   source_damage_isoDuctile_N
 
  enum, bind(c) 
    enumerator :: undefined_ID, &
                  damage_drivingforce_ID
  end enum                                                 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11 ToDo
- 
+
- integer(kind(undefined_ID)),         dimension(:,:),         allocatable,          private :: & 
+ type, private :: tParameters                                                                       !< container type for internal constitutive parameters
-   source_damage_isoDuctile_outputID                                                                  !< ID of each post result output
+   real(pReal) :: &
+     critPlasticStrain, &
+     N, &
+     aTol
+   integer(kind(undefined_ID)), allocatable, dimension(:) :: &
+     outputID
+ end type tParameters
+
+ type(tParameters), dimension(:), allocatable, private :: param                                     !< containers of constitutive parameters (len Ninstance)
 
 
  public :: &
@@ -52,30 +52,23 @@ contains
 !> @brief module initialization
 !> @details reads in material parameters, allocates arrays, and does sanity checks
 !--------------------------------------------------------------------------------------------------
-subroutine source_damage_isoDuctile_init(fileUnit)
+subroutine source_damage_isoDuctile_init
 #if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
  use, intrinsic :: iso_fortran_env, only: &
    compiler_version, &
    compiler_options
 #endif
+ use prec, only: &
+   pStringLen
  use debug, only: &
    debug_level,&
    debug_constitutive,&
    debug_levelBasic
  use IO, only: &
-   IO_read, &
-   IO_lc, &
-   IO_getTag, &
-   IO_isBlank, &
-   IO_stringPos, &
-   IO_stringValue, &
-   IO_floatValue, &
-   IO_intValue, &
    IO_warning, &
-   IO_error, &
+   IO_error
-   IO_timeStamp, &
-   IO_EOF
  use material, only: &
+   material_allocateSourceState, &
    phase_source, &
    phase_Nsources, &
    phase_Noutput, &
@@ -84,32 +77,31 @@ subroutine source_damage_isoDuctile_init(fileUnit)
    material_phase, &  
    sourceState
  use config, only: &
+   config_phase, &
    material_Nphase, &
    MATERIAL_partPhase
- 
- use numerics,only: &
-   numerics_integrator
 
  implicit none
- integer(pInt), intent(in) :: fileUnit
 
- integer(pInt), allocatable, dimension(:) :: chunkPos
+ integer(pInt) :: Ninstance,phase,instance,source,sourceOffset
- integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,source,sourceOffset,o
+ integer(pInt) :: NofMyPhase,p,i
- integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
+ character(len=65536),   dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
- integer(pInt) :: NofMyPhase   
+ integer(kind(undefined_ID)) :: &
- character(len=65536) :: &
+   outputID
-   tag  = '', &
-   line = ''
 
- write(6,'(/,a)')   ' <<<+-  source_'//SOURCE_damage_isoDuctile_label//' init  -+>>>'
+ character(len=pStringLen) :: &
- write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
+   extmsg = ''
+ character(len=65536), dimension(:), allocatable :: &
+   outputs
+
+ write(6,'(/,a)')   ' <<<+-  source_'//SOURCE_DAMAGE_ISODUCTILE_LABEL//' init  -+>>>'
 #include "compilation_info.f90"
 
- maxNinstance = int(count(phase_source == SOURCE_damage_isoDuctile_ID),pInt)
+ Ninstance = int(count(phase_source == SOURCE_damage_isoDuctile_ID),pInt)
- if (maxNinstance == 0_pInt) return
+ if (Ninstance == 0_pInt) return
  
  if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
-   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
+   write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
  
  allocate(source_damage_isoDuctile_offset(material_Nphase), source=0_pInt)
  allocate(source_damage_isoDuctile_instance(material_Nphase), source=0_pInt)
@@ -121,121 +113,64 @@ subroutine source_damage_isoDuctile_init(fileUnit)
    enddo    
  enddo
    
- allocate(source_damage_isoDuctile_sizePostResults(maxNinstance),                     source=0_pInt)
+ allocate(source_damage_isoDuctile_sizePostResult(maxval(phase_Noutput),Ninstance),source=0_pInt)
- allocate(source_damage_isoDuctile_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
+ allocate(source_damage_isoDuctile_output(maxval(phase_Noutput),Ninstance))
- allocate(source_damage_isoDuctile_output(maxval(phase_Noutput),maxNinstance))
           source_damage_isoDuctile_output = ''
- allocate(source_damage_isoDuctile_outputID(maxval(phase_Noutput),maxNinstance),      source=undefined_ID)
- allocate(source_damage_isoDuctile_Noutput(maxNinstance),                             source=0_pInt) 
- allocate(source_damage_isoDuctile_critPlasticStrain(maxNinstance),                     source=0.0_pReal) 
- allocate(source_damage_isoDuctile_N(maxNinstance),                                     source=0.0_pReal) 
- allocate(source_damage_isoDuctile_aTol(maxNinstance),                                source=0.0_pReal) 
 
- rewind(fileUnit)
+ allocate(param(Ninstance))
- phase = 0_pInt
- do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase)         ! wind forward to <phase>
-   line = IO_read(fileUnit)
- enddo
  
- parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
+ do p=1, size(config_phase)
-   line = IO_read(fileUnit)
+   if (all(phase_source(:,p) /= SOURCE_DAMAGE_ISODUCTILE_ID)) cycle
-   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
+   associate(prm => param(source_damage_isoDuctile_instance(p)), &
-   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
+             config => config_phase(p))
-     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
+             
-     exit                                                                                           
+   prm%aTol              = config%getFloat('isoductile_atol',defaultVal = 1.0e-3_pReal)
-   endif   
-   if (IO_getTag(line,'[',']') /= '') then                                                          ! next phase section
-     phase = phase + 1_pInt                                                                         ! advance phase section counter
-     cycle                                                                                          ! skip to next line
-   endif
-   if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_damage_isoDuctile_ID)) then   ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
-     instance = source_damage_isoDuctile_instance(phase)                                            ! which instance of my damage is present phase
-     chunkPos = IO_stringPos(line)
-     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
-     select case(tag)
-       case ('(output)')
-         select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
-           case ('isoductile_drivingforce')
-             source_damage_isoDuctile_Noutput(instance) = source_damage_isoDuctile_Noutput(instance) + 1_pInt
-             source_damage_isoDuctile_outputID(source_damage_isoDuctile_Noutput(instance),instance) = damage_drivingforce_ID
-             source_damage_isoDuctile_output(source_damage_isoDuctile_Noutput(instance),instance) = &
-                                                       IO_lc(IO_stringValue(line,chunkPos,2_pInt))
-          end select
 
-       case ('isoductile_criticalplasticstrain')
+   prm%N                 = config%getFloat('isoductile_ratesensitivity')
-         source_damage_isoDuctile_critPlasticStrain(instance) = IO_floatValue(line,chunkPos,2_pInt)
+   prm%critPlasticStrain = config%getFloat('isoductile_criticalplasticstrain')
+   
+   ! sanity checks
+   if (prm%aTol                 < 0.0_pReal) extmsg = trim(extmsg)//' isoductile_atol'
+   
+   if (prm%N                   <= 0.0_pReal) extmsg = trim(extmsg)//' isoductile_ratesensitivity'
+   if (prm%critPlasticStrain   <= 0.0_pReal) extmsg = trim(extmsg)//' isoductile_criticalplasticstrain'
+   
+!--------------------------------------------------------------------------------------------------
+!  exit if any parameter is out of range
+   if (extmsg /= '') &
+     call IO_error(211_pInt,ext_msg=trim(extmsg)//'('//SOURCE_DAMAGE_ISODUCTILE_LABEL//')')
 
-       case ('isoductile_ratesensitivity')
+!--------------------------------------------------------------------------------------------------
-         source_damage_isoDuctile_N(instance) = IO_floatValue(line,chunkPos,2_pInt)
+!  output pararameters
-
+   outputs = config%getStrings('(output)',defaultVal=emptyStringArray)
-       case ('isoductile_atol')
+   allocate(prm%outputID(0))
-         source_damage_isoDuctile_aTol(instance) = IO_floatValue(line,chunkPos,2_pInt)
+   do i=1_pInt, size(outputs)
+     outputID = undefined_ID
+     select case(outputs(i))
+     
+       case ('isoductile_drivingforce')
+         source_damage_isoDuctile_sizePostResult(i,source_damage_isoDuctile_instance(p)) = 1_pInt
+         source_damage_isoDuctile_output(i,source_damage_isoDuctile_instance(p)) = outputs(i)
+         prm%outputID = [prm%outputID, damage_drivingforce_ID]
 
      end select
-   endif; endif
- enddo parsingFile
 
+   enddo
 
-!--------------------------------------------------------------------------------------------------
+   end associate
-!  sanity checks
+   
- sanityChecks: do phase = 1_pInt, material_Nphase 
+   phase = p
-   myPhase: if (any(phase_source(:,phase) == SOURCE_damage_isoDuctile_ID)) then
+   NofMyPhase=count(material_phase==phase)
-     instance = source_damage_isoDuctile_instance(phase)
+   instance = source_damage_isoDuctile_instance(phase)
-     if (source_damage_isoDuctile_aTol(instance) < 0.0_pReal) &
+   sourceOffset = source_damage_isoDuctile_offset(phase)
-       source_damage_isoDuctile_aTol(instance) = 1.0e-3_pReal                                              ! default absolute tolerance 1e-3
-     if (source_damage_isoDuctile_critPlasticStrain(instance) <= 0.0_pReal) &
-       call IO_error(211_pInt,el=instance,ext_msg='critical plastic strain ('//SOURCE_damage_isoDuctile_LABEL//')')
-   endif myPhase
- enddo sanityChecks
 
- initializeInstances: do phase = 1_pInt, material_Nphase
+   call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1_pInt,1_pInt,0_pInt)
-   if (any(phase_source(:,phase) == SOURCE_damage_isoDuctile_ID)) then
+   sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_isoDuctile_sizePostResult(:,instance))
-     NofMyPhase=count(material_phase==phase)
+   sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol
-     instance = source_damage_isoDuctile_instance(phase)
+              
-     sourceOffset = source_damage_isoDuctile_offset(phase)
-!--------------------------------------------------------------------------------------------------
-!  Determine size of postResults array
-     outputsLoop: do o = 1_pInt,source_damage_isoDuctile_Noutput(instance)
-       select case(source_damage_isoDuctile_outputID(o,instance))
-         case(damage_drivingforce_ID)
-           mySize = 1_pInt
-       end select
  
-       if (mySize > 0_pInt) then  ! any meaningful output found
+ enddo
-          source_damage_isoDuctile_sizePostResult(o,instance) = mySize
-          source_damage_isoDuctile_sizePostResults(instance)  = source_damage_isoDuctile_sizePostResults(instance) + mySize
-       endif
-     enddo outputsLoop
-! Determine size of state array
-     sizeDotState              =   1_pInt
-     sizeDeltaState            =   0_pInt
-     sizeState                 =   1_pInt
-                
-     sourceState(phase)%p(sourceOffset)%sizeState      = sizeState
-     sourceState(phase)%p(sourceOffset)%sizeDotState   = sizeDotState
-     sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState
-     sourceState(phase)%p(sourceOffset)%sizePostResults = source_damage_isoDuctile_sizePostResults(instance)
-     allocate(sourceState(phase)%p(sourceOffset)%aTolState           (sizeState),                &
-              source=source_damage_isoDuctile_aTol(instance))
-     allocate(sourceState(phase)%p(sourceOffset)%state0              (sizeState,NofMyPhase),     source=.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%partionedState0     (sizeState,NofMyPhase),     source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%subState0           (sizeState,NofMyPhase),     source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%state               (sizeState,NofMyPhase),     source=0.0_pReal)
-
-     allocate(sourceState(phase)%p(sourceOffset)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     allocate(sourceState(phase)%p(sourceOffset)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
-     if (any(numerics_integrator == 1_pInt)) then
-       allocate(sourceState(phase)%p(sourceOffset)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)
-       allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     endif
-     if (any(numerics_integrator == 4_pInt)) &
-       allocate(sourceState(phase)%p(sourceOffset)%RK4dotState       (sizeDotState,NofMyPhase),  source=0.0_pReal)
-     if (any(numerics_integrator == 5_pInt)) &
-       allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState    (6,sizeDotState,NofMyPhase),source=0.0_pReal)
-
-   endif
  
- enddo initializeInstances
 end subroutine source_damage_isoDuctile_init
 
 !--------------------------------------------------------------------------------------------------
@@ -267,39 +202,34 @@ subroutine source_damage_isoDuctile_dotState(ipc, ip, el)
 
  sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = &
    sum(plasticState(phase)%slipRate(:,constituent))/ &
-   ((damage(homog)%p(damageOffset))**source_damage_isoDuctile_N(instance))/ & 
+   ((damage(homog)%p(damageOffset))**param(instance)%N)/ & 
-   source_damage_isoDuctile_critPlasticStrain(instance) 
+   param(instance)%critPlasticStrain 
 
 end subroutine source_damage_isoDuctile_dotState
  
 !--------------------------------------------------------------------------------------------------
 !> @brief returns local part of nonlocal damage driving force
 !--------------------------------------------------------------------------------------------------
-subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ipc, ip, el)
+subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
  use material, only: &
-   phaseAt, phasememberAt, &
    sourceState
 
  implicit none
  integer(pInt), intent(in) :: &
-   ipc, &                                                                                           !< component-ID of integration point
+   phase, &
-   ip, &                                                                                            !< integration point
+   constituent
-   el                                                                                               !< element
  real(pReal),  intent(in) :: &
    phi
  real(pReal),  intent(out) :: &
    localphiDot, &
    dLocalphiDot_dPhi
  integer(pInt) :: &
-   phase, constituent, sourceOffset
+   sourceOffset
 
- phase = phaseAt(ipc,ip,el)
- constituent = phasememberAt(ipc,ip,el)
  sourceOffset = source_damage_isoDuctile_offset(phase)
  
- localphiDot = 1.0_pReal - &
+ localphiDot = 1.0_pReal &
-               sourceState(phase)%p(sourceOffset)%state(1,constituent)* &
+             - sourceState(phase)%p(sourceOffset)%state(1,constituent) * phi
-               phi
  
  dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent)
  
@@ -308,33 +238,28 @@ end subroutine source_damage_isoDuctile_getRateAndItsTangent
 !--------------------------------------------------------------------------------------------------
 !> @brief return array of local damage results
 !--------------------------------------------------------------------------------------------------
-function source_damage_isoDuctile_postResults(ipc,ip,el)
+function source_damage_isoDuctile_postResults(phase, constituent)
  use material, only: &
-   phaseAt, phasememberAt, &
    sourceState
 
  implicit none
- integer(pInt),              intent(in) :: &
+ integer(pInt), intent(in) :: &
-   ipc, &                                                                                           !< component-ID of integration point
+   phase, &
-   ip, &                                                                                            !< integration point
+   constituent
-   el                                                                                               !< element
+ real(pReal), dimension(sum(source_damage_isoDuctile_sizePostResult(:, &
- real(pReal), dimension(source_damage_isoDuctile_sizePostResults( &
+                          source_damage_isoDuctile_instance(phase)))) :: &
-                          source_damage_isoDuctile_instance(phaseAt(ipc,ip,el)))) :: &
    source_damage_isoDuctile_postResults
 
  integer(pInt) :: &
-   instance, phase, constituent, sourceOffset, o, c
+   instance, sourceOffset, o, c
    
- phase = phaseAt(ipc,ip,el)
- constituent = phasememberAt(ipc,ip,el)
  instance = source_damage_isoDuctile_instance(phase)
  sourceOffset = source_damage_isoDuctile_offset(phase)
 
  c = 0_pInt
- source_damage_isoDuctile_postResults = 0.0_pReal
 
- do o = 1_pInt,source_damage_isoDuctile_Noutput(instance)
+ do o = 1_pInt,size(param(instance)%outputID)
-    select case(source_damage_isoDuctile_outputID(o,instance))
+    select case(param(instance)%outputID(o))
       case (damage_drivingforce_ID)
         source_damage_isoDuctile_postResults(c+1_pInt) = sourceState(phase)%p(sourceOffset)%state(1,constituent)
         c = c + 1