Merge branch 'development' into pytest-curl-div-grad-2

PRIVATE

@@ -1 +1 @@
-Subproject commit 232a094c715bcbbd1c6652c4dc4a4a50d402b82f
+Subproject commit c595994cd8880acadf50b5dedb79156d04d35b91

VERSION

@@ -1 +1 @@
-v2.0.3-2375-gc6d58954
+v2.0.3-2402-gb88f5ec0

installation/symlink_Processing.py

@@ -16,10 +16,6 @@ if not os.path.isdir(binDir):
 #define ToDo list
 processing_subDirs    = ['pre',
                          'post',
-                         'misc',
-                        ]
-processing_extensions = ['.py',
-                         '.sh',
                         ]
 
 sys.stdout.write('\nsymbolic linking...\n')
@@ -31,7 +27,7 @@ for subDir in processing_subDirs:
 
   for theFile in os.listdir(theDir):
     theName,theExt = os.path.splitext(theFile)
-    if theExt in processing_extensions:                                                             # only consider files with proper extensions
+    if theExt in ['.py']:
 
       src      = os.path.abspath(os.path.join(theDir,theFile))
       sym_link = os.path.abspath(os.path.join(binDir,theName))

processing/misc/ang_toTable.py

@@ -1,30 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import sys
-from io import StringIO
-from optparse import OptionParser
-
-import damask
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-#--------------------------------------------------------------------------------------------------
-#                                MAIN
-#--------------------------------------------------------------------------------------------------
-
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog [angfile[s]]', description = """
-Convert TSL/EDAX *.ang file to ASCIItable
-
-""", version = scriptID)
-
-(options, filenames) = parser.parse_args()
-if filenames == []: filenames = [None]
-
-for name in filenames:
-    damask.util.report(scriptName,name)
-
-    table = damask.Table.from_ang(StringIO(''.join(sys.stdin.read())) if name is None else name)
-    table.to_ASCII(sys.stdout if name is None else os.path.splitext(name)[0]+'.txt')

processing/post/addCauchy.py

@@ -1,49 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import sys
-from io import StringIO
-from optparse import OptionParser
-
-import damask
-
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-# --------------------------------------------------------------------
-#                                MAIN
-# --------------------------------------------------------------------
-
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
-Add column containing Cauchy stress based on deformation gradient and first Piola--Kirchhoff stress.
-
-""", version = scriptID)
-
-parser.add_option('-f','--defgrad',
-                  dest = 'defgrad',
-                  type = 'string', metavar = 'string',
-                  help = 'heading of columns containing deformation gradient [%default]')
-parser.add_option('-p','--stress',
-                  dest = 'stress',
-                  type = 'string', metavar = 'string',
-                  help = 'heading of columns containing first Piola--Kirchhoff stress [%default]')
-
-parser.set_defaults(defgrad = 'f',
-                    stress  = 'p',
-                   )
-
-(options,filenames) = parser.parse_args()
-if filenames == []: filenames = [None]
-
-for name in filenames:
-    damask.util.report(scriptName,name)
-
-    table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
-    table.add('Cauchy',
-              damask.mechanics.Cauchy(table.get(options.stress ).reshape(-1,3,3),
-                                      table.get(options.defgrad).reshape(-1,3,3)).reshape(-1,9),
-              scriptID+' '+' '.join(sys.argv[1:]))
-
-    table.to_ASCII(sys.stdout if name is None else name)

processing/post/addDeterminant.py

@@ -1,45 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import sys
-from io import StringIO
-from optparse import OptionParser
-
-import numpy as np
-
-import damask
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-# --------------------------------------------------------------------
-#                                MAIN
-# --------------------------------------------------------------------
-
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
-Add column(s) containing determinant of requested tensor column(s).
-
-""", version = scriptID)
-
-parser.add_option('-t','--tensor',
-                  dest = 'tensor',
-                  action = 'extend', metavar = '<string LIST>',
-                  help = 'heading of columns containing tensor field values')
-
-(options,filenames) = parser.parse_args()
-if filenames == []: filenames = [None]
-
-if options.tensor is None:
-  parser.error('no data column specified.')
-
-for name in filenames:
-    damask.util.report(scriptName,name)
-
-    table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
-    for tensor in options.tensor:
-         table.add('det({})'.format(tensor),
-                   np.linalg.det(table.get(tensor).reshape(-1,3,3)),
-                   scriptID+' '+' '.join(sys.argv[1:]))
-
-    table.to_ASCII(sys.stdout if name is None else name)

processing/post/addDeviator.py

@@ -1,51 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import sys
-from io import StringIO
-from optparse import OptionParser
-
-import damask
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-# --------------------------------------------------------------------
-#                                MAIN
-# --------------------------------------------------------------------
-
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(2)]', description = """
-Add column(s) containing deviator of requested tensor column(s).
-
-""", version = scriptID)
-
-parser.add_option('-t','--tensor',
-                  dest = 'tensor',
-                  action = 'extend', metavar='<string LIST>',
-                  help = 'heading of columns containing tensor field values')
-parser.add_option('-s','--spherical',
-                  dest = 'spherical',
-                  action = 'store_true',
-                  help = 'report spherical part of tensor (hydrostatic component, pressure)')
-
-(options,filenames) = parser.parse_args()
-if filenames == []: filenames = [None]
-
-if options.tensor is None:
-    parser.error('no data column specified...')
-
-for name in filenames:
-    damask.util.report(scriptName,name)
-    
-    table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
-    for tensor in options.tensor:
-         table.add('dev({})'.format(tensor),
-                   damask.mechanics.deviatoric_part(table.get(tensor).reshape(-1,3,3)).reshape(-1,9),
-                   scriptID+' '+' '.join(sys.argv[1:]))
-         if options.spherical:
-             table.add('sph({})'.format(tensor),
-                       damask.mechanics.spherical_part(table.get(tensor).reshape(-1,3,3)),
-                       scriptID+' '+' '.join(sys.argv[1:]))
-
-    table.to_ASCII(sys.stdout if name is None else name)

processing/post/addInfo.py

@@ -1,41 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import sys
-from io import StringIO
-from optparse import OptionParser
-
-import damask
-
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-# --------------------------------------------------------------------
-#                                MAIN
-# --------------------------------------------------------------------
-
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
-Add info lines to ASCIItable header.
-
-""", version = scriptID)
-
-parser.add_option('-i',
-                  '--info',
-                  dest = 'info', action = 'extend', metavar = '<string LIST>',
-                  help    = 'items to add')
-
-(options,filenames) = parser.parse_args()
-if filenames == []: filenames = [None]
-
-if options.info is None:
-  parser.error('no info specified.')
-
-for name in filenames:
-    damask.util.report(scriptName,name)
-
-    table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
-    table.comments += options.info
-
-    table.to_ASCII(sys.stdout if name is None else name)

processing/post/addMises.py

@@ -1,56 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import sys
-from io import StringIO
-from optparse import OptionParser
-
-import damask
-
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-# --------------------------------------------------------------------
-#                                MAIN
-# --------------------------------------------------------------------
-
-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)
-
-parser.add_option('-e','--strain',
-                  dest = 'strain',
-                  action = 'extend', metavar = '<string LIST>',
-                  help = 'heading(s) of columns containing strain tensors')
-parser.add_option('-s','--stress',
-                  dest = '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 options.stress is [] and options.strain is []:
-    parser.error('no data column specified...')
-
-if filenames == []: filenames = [None]
-
-for name in filenames:
-    damask.util.report(scriptName,name)
-
-    table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
-    for strain in options.strain:
-        table.add('Mises({})'.format(strain),
-                  damask.mechanics.Mises_strain(damask.mechanics.symmetric(table.get(strain).reshape(-1,3,3))),
-                  scriptID+' '+' '.join(sys.argv[1:]))
-    for stress in options.stress:
-        table.add('Mises({})'.format(stress),
-                  damask.mechanics.Mises_stress(damask.mechanics.symmetric(table.get(stress).reshape(-1,3,3))),
-                  scriptID+' '+' '.join(sys.argv[1:]))
-    
-    table.to_ASCII(sys.stdout if name is None else name)

processing/post/addOrientations.py

@@ -172,7 +172,7 @@ for name in filenames:
       
     elif inputtype == 'matrix':
       d = representations['matrix'][1]
-      o = damask.Rotation.fromMatrix(list(map(float,table.data[column:column+d])))
+      o = damask.Rotation.fromMatrix(np.array(list(map(float,table.data[column:column+d]))).reshape(3,3))
 
     elif inputtype == 'frame':
       M = np.array(list(map(float,table.data[column[0]:column[0]+3] + \

processing/post/addSchmidfactors.py

@@ -214,7 +214,7 @@ for name in filenames:
 
   outputAlive = True
   while outputAlive and table.data_read():                                                          # read next data line of ASCII table
-    o = damask.Rotation(list(map(float,table.data[column:column+4])))
+    o = damask.Rotation(np.array(list(map(float,table.data[column:column+4]))))
 
     table.data_append(  np.abs(  np.sum(slip_direction * (o * force) ,axis=1) \
                                * np.sum(slip_normal    * (o * normal),axis=1)))

processing/post/addTable.py

@@ -1,44 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import sys
-from io import StringIO
-from optparse import OptionParser
-
-import damask
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-# --------------------------------------------------------------------
-#                                MAIN
-# --------------------------------------------------------------------
-
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
-Append data of ASCIItable(s) column-wise.
-
-""", version = scriptID)
-
-parser.add_option('-a', '--add','--table',
-                  dest = 'table',
-                  action = 'extend', metavar = '<string LIST>',
-                  help = 'tables to add')
-
-(options,filenames) = parser.parse_args()
-if filenames == []: filenames = [None]
-
-if options.table is None:
-      parser.error('no table specified.')
-
-for name in filenames:
-    damask.util.report(scriptName,name)
-
-    table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
-
-    for addTable in options.table:
-        table2 = damask.Table.from_ASCII(addTable)
-        table2.data = table2.data[:table.data.shape[0]]
-        table.join(table2)
-
-    table.to_ASCII(sys.stdout if name is None else name)

processing/post/growTable.py

@@ -1,43 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import sys
-from io import StringIO
-from optparse import OptionParser
-
-import damask
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-# --------------------------------------------------------------------
-#                                MAIN
-# --------------------------------------------------------------------
-
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
-Append data of ASCIItable(s) row-wise.
-
-""", version = scriptID)
-
-parser.add_option('-a', '--add','--table',
-                  dest = 'table',
-                  action = 'extend', metavar = '<string LIST>',
-                  help = 'tables to add')
-
-(options,filenames) = parser.parse_args()
-if filenames == []: filenames = [None]
-
-if options.table is None:
-    parser.error('no table specified.')
-
-for name in filenames:
-    damask.util.report(scriptName,name)
-
-    table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
-
-    for growTable in options.table:
-        table2 = damask.Table.from_ASCII(growTable)
-        table.append(table2)
-
-    table.to_ASCII(sys.stdout if name is None else name)

processing/post/permuteData.py

@@ -2,6 +2,7 @@
 
 import os
 import sys
+from io import StringIO
 from optparse import OptionParser
 
 import numpy as np
@@ -41,73 +42,20 @@ parser.set_defaults(label = [],
                    )
 
 (options,filenames) = parser.parse_args()
-
-if len(options.label) == 0:
-  parser.error('no labels specified.')
-
-# --- loop over input files -------------------------------------------------------------------------
-
 if filenames == []: filenames = [None]
 
 for name in filenames:
-  try:
+    damask.util.report(scriptName,name)
-    table = damask.ASCIItable(name = name)
-  except IOError:
-    continue
-  damask.util.report(scriptName,name)
 
-# ------------------------------------------ read header ------------------------------------------
+    table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
 
-  table.head_read()
+    randomSeed = int(os.urandom(4).hex(), 16) if options.randomSeed is None else options.randomSeed # random seed per file
+    rng = np.random.default_rng(randomSeed)
 
-# ------------------------------------------ process labels ---------------------------------------
+    for label in options.label:
+        data = table.get(label)
+        uniques,inverse  = np.unique(data,return_inverse=True,axis=0) if options.unique else (data,np.arange(len(data)))
+        rng.shuffle(uniques)
+        table.set(label,uniques[inverse], scriptID+' '+' '.join(sys.argv[1:]))
 
-  errors  = []
+    table.to_ASCII(sys.stdout if name is None else name)
-  remarks = []
-  columns = []
-  dims    = []
-
-  indices    = table.label_index    (options.label)
-  dimensions = table.label_dimension(options.label)
-  for i,index in enumerate(indices):
-    if index == -1: remarks.append('label "{}" not present...'.format(options.label[i]))
-    else:
-      columns.append(index)
-      dims.append(dimensions[i])
-
-  if remarks != []: damask.util.croak(remarks)
-  if errors  != []:
-    damask.util.croak(errors)
-    table.close(dismiss = True)
-    continue
-
-# ------------------------------------------ assemble header ---------------------------------------
-
-  randomSeed = int(os.urandom(4).hex(), 16) if options.randomSeed is None else options.randomSeed   # random seed per file
-  np.random.seed(randomSeed)
-
-  table.info_append([scriptID + '\t' + ' '.join(sys.argv[1:]),
-                     'random seed {}'.format(randomSeed),
-                    ])
-  table.head_write()
-
-# ------------------------------------------ process data ------------------------------------------
-
-  table.data_readArray()                                                                            # read all data at once
-  for col,dim in zip(columns,dims):
-    if options.unique:
-      s = set(map(tuple,table.data[:,col:col+dim]))                                                 # generate set of (unique) values
-      uniques = np.array(list(map(np.array,s)))                                                     # translate set to np.array
-      shuffler = dict(zip(s,np.random.permutation(len(s))))                                         # random permutation
-      table.data[:,col:col+dim] = uniques[np.array(list(map(lambda x: shuffler[tuple(x)],
-                                                            table.data[:,col:col+dim])))]           # fill table with mapped uniques
-    else:
-      np.random.shuffle(table.data[:,col:col+dim])                                                  # independently shuffle every row
-
-# ------------------------------------------ output result -----------------------------------------
-
-  table.data_writeArray()
-
-# ------------------------------------------ output finalization -----------------------------------
-
-  table.close()                                                                                     # close ASCII tables

processing/post/reLabel.py

@@ -1,50 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import sys
-from io import StringIO
-from optparse import OptionParser
-
-import damask
-
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-# --------------------------------------------------------------------
-#                                MAIN
-# --------------------------------------------------------------------
-
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
-Rename scalar, vectorial, and/or tensorial data header labels.
-
-""", version = scriptID)
-
-parser.add_option('-l','--label',
-                  dest = 'label',
-                  action = 'extend', metavar='<string LIST>',
-                  help = 'column(s) to rename')
-parser.add_option('-s','--substitute',
-                  dest = 'substitute',
-                  action = 'extend', metavar='<string LIST>',
-                  help = 'new column label(s)')
-
-parser.set_defaults(label = [],
-                    substitute = [],
-                   )
-
-(options,filenames) = parser.parse_args()
-if filenames == []: filenames = [None]
-
-if len(options.label) != len(options.substitute):
-    parser.error('number of column labels and substitutes do not match.')
-
-for name in filenames:
-    damask.util.report(scriptName,name)
-
-    table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
-    for label,substitute in zip(options.label,options.substitute):
-        table.rename(label,substitute,scriptID+' '+' '.join(sys.argv[1:]))
-
-    table.to_ASCII(sys.stdout if name is None else name)

processing/post/scaleData.py

@@ -1,49 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import sys
-from io import StringIO
-from optparse import OptionParser
-
-import damask
-
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-# --------------------------------------------------------------------
-#                                MAIN
-# --------------------------------------------------------------------
-
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
-Uniformly scale column values by given factor.
-
-""", version = scriptID)
-
-parser.add_option('-l','--label',
-                  dest = 'labels',
-                  action = 'extend', metavar = '<string LIST>',
-                  help  ='column(s) to scale')
-parser.add_option('-f','--factor',
-                  dest = 'factor',
-                  action = 'extend', metavar='<float LIST>',
-                  help = 'factor(s) per column')
-
-parser.set_defaults(label  = [],
-                    factor = [])
-
-(options,filenames) = parser.parse_args()
-if filenames == []: filenames = [None]
-
-if len(options.labels) != len(options.factor):
-    parser.error('number of column labels and factors do not match.')
-
-for name in filenames:
-    damask.util.report(scriptName,name)
-
-    table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
-    for label,factor in zip(options.labels,options.factor):
-        table.set(label,table.get(label)*float(factor),scriptID+' '+' '.join(sys.argv[1:]))
-
-    table.to_ASCII(sys.stdout if name is None else name)

processing/post/shiftData.py

@@ -1,49 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import sys
-from io import StringIO
-from optparse import OptionParser
-
-import damask
-
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-# --------------------------------------------------------------------
-#                                MAIN
-# --------------------------------------------------------------------
-
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
-Uniformly shift column values by given offset.
-
-""", version = scriptID)
-
-parser.add_option('-l','--label',
-                  dest = 'labels',
-                  action = 'extend', metavar = '<string LIST>',
-                  help  ='column(s) to shift')
-parser.add_option('-o','--offset',
-                  dest = 'offset',
-                  action = 'extend', metavar='<float LIST>',
-                  help = 'offset(s) per column')
-
-parser.set_defaults(label  = [],
-                    offset = [])
-
-(options,filenames) = parser.parse_args()
-if filenames == []: filenames = [None]
-
-if len(options.labels) != len(options.offset):
-    parser.error('number of column labels and offsets do not match.')
-
-for name in filenames:
-    damask.util.report(scriptName,name)
-
-    table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
-    for label,offset in zip(options.labels,options.offset):
-        table.set(label,table.get(label)+float(offset),scriptID+' '+' '.join(sys.argv[1:]))
-
-    table.to_ASCII(sys.stdout if name is None else name)

processing/post/sortTable.py

@@ -1,50 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import sys
-from io import StringIO
-from optparse import OptionParser
-
-import damask
-
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-# --------------------------------------------------------------------
-#                                MAIN
-# --------------------------------------------------------------------
-
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
-Sort rows by given (or all) column label(s).
-
-Examples:
-With coordinates in columns "x", "y", and "z"; sorting with x slowest and z fastest varying index: --label x,y,z.
-""", version = scriptID)
-
-
-parser.add_option('-l','--label',
-                  dest   = 'labels',
-                  action = 'extend', metavar = '<string LIST>',
-                  help   = 'list of column labels (a,b,c,...)')
-parser.add_option('-r','--reverse',
-                  dest   = 'reverse',
-                  action = 'store_true',
-                  help   = 'sort in reverse')
-
-parser.set_defaults(reverse = False,
-                   )
-
-(options,filenames) = parser.parse_args()
-if filenames == []: filenames = [None]
-
-if options.labels is None:
-    parser.error('no labels specified.')
-for name in filenames:
-    damask.util.report(scriptName,name)
-
-    table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
-    table.sort_by(options.labels,not options.reverse)
-
-    table.to_ASCII(sys.stdout if name is None else name)

processing/pre/geom_clean.py

@@ -1,42 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import sys
-from io import StringIO
-from optparse import OptionParser
-
-import damask
-
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-#--------------------------------------------------------------------------------------------------
-#                                MAIN
-#--------------------------------------------------------------------------------------------------
-
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [geomfile(s)]', description = """
-Smooth microstructure by selecting most frequent index within given stencil at each location.
-
-""", version=scriptID)
-
-parser.add_option('-s','--stencil',
-                  dest = 'stencil',
-                  type = 'int', metavar = 'int',
-                  help = 'size of smoothing stencil [%default]')
-
-parser.set_defaults(stencil = 3)
-
-(options, filenames) = parser.parse_args()
-
-
-if filenames == []: filenames = [None]
-
-for name in filenames:
-  damask.util.report(scriptName,name)
-
-  geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name)
-  damask.util.croak(geom.clean(options.stencil))
-  geom.add_comments(scriptID + ' ' + ' '.join(sys.argv[1:]))
-  geom.to_file(sys.stdout if name is None else name,pack=False)

processing/pre/geom_fromScratch.py

@@ -1,66 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import sys
-from optparse import OptionParser
-
-import numpy as np
-
-import damask
-
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-# --------------------------------------------------------------------
-#                                MAIN
-# --------------------------------------------------------------------
-
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [geomfile(s)]', description = """
-Generate homogeneous geometry.
-
-""", version = scriptID)
-
-parser.add_option('-g','--grid',
-                  dest = 'grid',
-                  type = 'int', nargs = 3, metavar = ' '.join(['int']*3),
-                  help = 'a,b,c grid of hexahedral box %default')
-parser.add_option('-s', '--size',
-                  dest = 'size',
-                  type = 'float', nargs = 3, metavar = ' '.join(['float']*3),
-                  help = 'x,y,z of geometry size')
-parser.add_option('-o','--origin',
-                  dest = 'origin',
-                  type = 'float', nargs = 3, metavar = ' '.join(['float']*3),
-                  help = 'x,y,z of geometry origin %default')
-parser.add_option('--homogenization',
-                  dest = 'homogenization',
-                  type = 'int', metavar = 'int',
-                  help = 'homogenization index [%default]')
-parser.add_option('-f','--fill',
-                  dest = 'fill',
-                  type = 'float', metavar = 'int',
-                  help = 'microstructure index [%default]')
-
-parser.set_defaults(grid = (16,16,16),
-                    origin = (0.,0.,0.),
-                    homogenization = 1,
-                    fill = 1,
-                   )
-
-(options, filename) = parser.parse_args()
-
-
-name = None if filename == [] else filename[0]
-damask.util.report(scriptName,name)
-
-dtype = float if np.isnan(options.fill) or int(options.fill) != options.fill else int
-geom = damask.Geom(microstructure=np.full(options.grid,options.fill,dtype=dtype),
-                   size=options.size,
-                   origin=options.origin,
-                   homogenization=options.homogenization,
-                   comments=scriptID + ' ' + ' '.join(sys.argv[1:]))
-damask.util.croak(geom)
-
-geom.to_file(sys.stdout if name is None else name,pack=False)

processing/pre/geom_mirror.py

@@ -1,47 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import sys
-from io import StringIO
-from optparse import OptionParser
-
-import damask
-
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-#--------------------------------------------------------------------------------------------------
-#                                MAIN
-#--------------------------------------------------------------------------------------------------
-
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [geomfile(s)]', description = """
-Mirror along given directions.
-
-""", version=scriptID)
-
-validDirections = ['x','y','z']
-
-parser.add_option('-d','--direction',
-                  dest = 'directions',
-                  action = 'extend', metavar = '<string LIST>',
-                  help = "directions in which to mirror {{{}}}".format(','.join(validDirections)))
-parser.add_option(    '--reflect',
-                  dest = 'reflect',
-                  action = 'store_true',
-                  help = 'reflect (include) outermost layers')
-
-parser.set_defaults(reflect = False)
-
-(options, filenames) = parser.parse_args()
-
-if filenames == []: filenames = [None]
-
-for name in filenames:
-  damask.util.report(scriptName,name)
-
-  geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name)
-  damask.util.croak(geom.mirror(options.directions,options.reflect))
-  geom.add_comments(scriptID + ' ' + ' '.join(sys.argv[1:]))
-  geom.to_file(sys.stdout if name is None else name,pack=False)

processing/pre/geom_pack.py

@@ -1,37 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import sys
-from io import StringIO
-from optparse import OptionParser
-
-import damask
-
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-#--------------------------------------------------------------------------------------------------
-#                                MAIN
-#--------------------------------------------------------------------------------------------------
-
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog [geomfile(s)]', description = """
-Pack ranges to "a to b" and/or multiples to "n of x".
-
-""", version = scriptID)
-
-(options, filenames) = parser.parse_args()
-
-
-if filenames == []: filenames = [None]
-
-for name in filenames:
-  damask.util.report(scriptName,name)
-
-  geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name)
-
-  damask.util.croak(geom)
-  geom.add_comments(scriptID + ' ' + ' '.join(sys.argv[1:]))
- 
-  geom.to_file(sys.stdout if name is None else name,pack=True)

processing/pre/geom_rescale.py

@@ -1,60 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import sys
-from io import StringIO
-from optparse import OptionParser
-
-import numpy as np
-
-import damask
-
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-#--------------------------------------------------------------------------------------------------
-#                                MAIN
-#--------------------------------------------------------------------------------------------------
-
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [geomfile(s)]', description = """
-Scales independently in x, y, and z direction in terms of grid and/or size.
-Either absolute values or relative factors (like "0.25x") can be used.
-
-""", version = scriptID)
-
-parser.add_option('-g', '--grid',
-                  dest = 'grid',
-                  type = 'string', nargs = 3, metavar = 'string string string',
-                  help = 'a,b,c grid of hexahedral box')
-parser.add_option('-s', '--size',
-                  dest = 'size',
-                  type = 'string', nargs = 3, metavar = 'string string string',
-                  help = 'x,y,z size of hexahedral box')
-
-(options, filenames) = parser.parse_args()
-
-
-if filenames == []: filenames = [None]
-
-for name in filenames:
-  damask.util.report(scriptName,name)
-
-  geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name)
-
-  grid = geom.get_grid()
-  size = geom.get_size()
-
-  new_grid = grid if options.grid is None else \
-             np.array([int(o*float(n.lower().replace('x',''))) if n.lower().endswith('x') \
-                  else int(n) for o,n in zip(grid,options.grid)],dtype=int)
-
-  new_size = size if options.size is None else \
-             np.array([o*float(n.lower().replace('x','')) if n.lower().endswith('x') \
-                  else float(n) for o,n in zip(size,options.size)],dtype=float)
-
-  geom.scale(new_grid)
-  damask.util.croak(geom.update(microstructure = None,size = new_size))
-  geom.add_comments(scriptID + ' ' + ' '.join(sys.argv[1:]))
-  geom.to_file(sys.stdout if name is None else name,pack=False)

processing/pre/geom_toTable.py

@@ -1,45 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import sys
-from io import StringIO
-from optparse import OptionParser
-
-import damask
-
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-#--------------------------------------------------------------------------------------------------
-#                                MAIN
-#--------------------------------------------------------------------------------------------------
-
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog [geomfile(s)]', description = """
-Translate geom description into ASCIItable containing position and microstructure.
-
-""", version = scriptID)
-
-(options, filenames) = parser.parse_args()
-if filenames == []: filenames = [None]
-
-for name in filenames:
-    damask.util.report(scriptName,name)
-
-    geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name)
-    damask.util.croak(geom)
-
-    coord0 = damask.grid_filters.cell_coord0(geom.grid,geom.size,geom.origin).reshape(-1,3)
-
-    comments = geom.comments \
-             + [scriptID + ' ' + ' '.join(sys.argv[1:]),
-                'grid\ta {}\tb {}\tc {}'.format(*geom.grid),
-                'size\tx {}\ty {}\tz {}'.format(*geom.size),
-                'origin\tx {}\ty {}\tz {}'.format(*geom.origin),
-                'homogenization\t{}'.format(geom.homogenization)]
-
-    table = damask.Table(coord0,{'pos':(3,)},comments)
-    table.add('microstructure',geom.microstructure.reshape((-1,1),order='F'))
-
-    table.to_ASCII(sys.stdout if name is None else os.path.splitext(name)[0]+'.txt')

processing/pre/geom_unpack.py

@@ -1,37 +0,0 @@
-#!/usr/bin/env python3
-
-import os
-import sys
-from optparse import OptionParser
-from io import StringIO
-
-import damask
-
-
-scriptName = os.path.splitext(os.path.basename(__file__))[0]
-scriptID   = ' '.join([scriptName,damask.version])
-
-
-#--------------------------------------------------------------------------------------------------
-#                                MAIN
-#--------------------------------------------------------------------------------------------------
-
-parser = OptionParser(option_class=damask.extendableOption, usage='%prog [geomfile(s)]', description = """
-Unpack ranges "a to b" and/or "n of x" multiples (exclusively in one line).
-
-""", version = scriptID)
-
-(options, filenames) = parser.parse_args()
-
-
-if filenames == []: filenames = [None]
-
-for name in filenames:
-  damask.util.report(scriptName,name)
-
-  geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name)
-
-  damask.util.croak(geom)
-  geom.add_comments(scriptID + ' ' + ' '.join(sys.argv[1:]))
-
-  geom.to_file(sys.stdout if name is None else name,pack=False)

python/.coveragerc

@@ -1,2 +1,5 @@
 [run]
 omit = tests/*
+       damask/_asciitable.py
+       damask/_test.py
+       damask/config/*

python/damask/_orientation.py

@@ -38,6 +38,9 @@ class Orientation:
         else:
             self.rotation = Rotation.fromQuaternion(rotation)                                       # assume quaternion
 
+        if self.rotation.quaternion.shape != (4,):
+            raise NotImplementedError('Support for multiple rotations missing')
+
     def disorientation(self,
                        other,
                        SST = True,

python/damask/_result.py

@@ -465,6 +465,14 @@ class Result:
             with h5py.File(self.fname,'r') as f:
                 return f['geometry/x_c'][()]
 
+    def node_coordinates(self):
+        """Return initial coordinates of the cell centers."""
+        if self.structured:
+            return grid_filters.node_coord0(self.grid,self.size,self.origin).reshape(-1,3,order='F')
+        else:
+            with h5py.File(self.fname,'r') as f:
+                return f['geometry/x_n'][()]
+
 
     @staticmethod
     def _add_absolute(x):

python/damask/_rotation.py

@@ -1,6 +1,7 @@
 import numpy as np
 
 from ._Lambert import ball_to_cube, cube_to_ball
+from . import mechanics
 
 _P = -1
 
@@ -61,6 +62,8 @@ class Rotation:
 
     def __repr__(self):
         """Orientation displayed as unit quaternion, rotation matrix, and Bunge-Euler angles."""
+        if self.quaternion.shape != (4,):
+            raise NotImplementedError('Support for multiple rotations missing')
         return '\n'.join([
                'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)),
                'Matrix:\n{}'.format(self.asMatrix()),
@@ -83,6 +86,8 @@ class Rotation:
         considere rotation of (3,3,3,3)-matrix
 
         """
+        if self.quaternion.shape != (4,):
+            raise NotImplementedError('Support for multiple rotations missing')
         if isinstance(other, Rotation):                                                             # rotate a rotation
             self_q  = self.quaternion[0]
             self_p  = self.quaternion[1:]
@@ -107,7 +112,7 @@ class Rotation:
             elif other.shape == (3,3,):                                                             # rotate a single (3x3)-matrix
                 return np.dot(self.asMatrix(),np.dot(other,self.asMatrix().T))
             elif other.shape == (3,3,3,3,):
-                raise NotImplementedError
+                raise NotImplementedError('Support for rotation of 4th order tensors missing')
             else:
                 return NotImplemented
         else:
@@ -116,7 +121,7 @@ class Rotation:
 
     def inverse(self):
         """In-place inverse rotation/backward rotation."""
-        self.quaternion[1:] *= -1
+        self.quaternion[...,1:] *= -1
         return self
 
     def inversed(self):
@@ -125,12 +130,12 @@ class Rotation:
 
 
     def standardize(self):
-        """In-place quaternion representation with positive q."""
+        """In-place quaternion representation with positive real part."""
-        if self.quaternion[0] < 0.0: self.quaternion*=-1
+        self.quaternion[self.quaternion[...,0] < 0.0] *= -1
         return self
 
     def standardized(self):
-        """Quaternion representation with positive q."""
+        """Quaternion representation with positive real part."""
         return self.copy().standardize()
 
 
@@ -157,15 +162,17 @@ class Rotation:
             Rotation from which the average is rotated.
 
         """
+        if self.quaternion.shape != (4,) or other.quaternion.shape != (4,):
+            raise NotImplementedError('Support for multiple rotations missing')
         return Rotation.fromAverage([self,other])
 
 
     ################################################################################################
     # convert to different orientation representations (numpy arrays)
 
-    def asQuaternion(self):
+    def as_quaternion(self):
         """
-        Unit quaternion [q, p_1, p_2, p_3] unless quaternion == True: damask.quaternion object.
+        Unit quaternion [q, p_1, p_2, p_3].
 
         Parameters
         ----------
@@ -175,8 +182,8 @@ class Rotation:
         """
         return self.quaternion
 
-    def asEulers(self,
+    def as_Eulers(self,
-                 degrees = False):
+                  degrees = False):
         """
         Bunge-Euler angles: (φ_1, ϕ, φ_2).
 
@@ -190,9 +197,9 @@ class Rotation:
         if degrees: eu = np.degrees(eu)
         return eu
 
-    def asAxisAngle(self,
+    def as_axis_angle(self,
-                    degrees = False,
+                      degrees = False,
-                    pair = False):
+                      pair = False):
         """
         Axis angle representation [n_1, n_2, n_3, ω] unless pair == True: ([n_1, n_2, n_3], ω).
 
@@ -205,15 +212,15 @@ class Rotation:
 
         """
         ax = Rotation.qu2ax(self.quaternion)
-        if degrees: ax[3] = np.degrees(ax[3])
+        if degrees: ax[...,3] = np.degrees(ax[...,3])
-        return (ax[:3],ax[3]) if pair else ax
+        return (ax[...,:3],ax[...,3]) if pair else ax
 
-    def asMatrix(self):
+    def as_matrix(self):
         """Rotation matrix."""
         return Rotation.qu2om(self.quaternion)
 
-    def asRodrigues(self,
+    def as_Rodrigues(self,
-                    vector = False):
+                     vector = False):
         """
         Rodrigues-Frank vector representation [n_1, n_2, n_3, tan(ω/2)] unless vector == True: [n_1, n_2, n_3] * tan(ω/2).
 
@@ -224,9 +231,9 @@ class Rotation:
 
         """
         ro = Rotation.qu2ro(self.quaternion)
-        return ro[:3]*ro[3] if vector else ro
+        return ro[...,:3]*ro[...,3] if vector else ro
 
-    def asHomochoric(self):
+    def as_homochoric(self):
         """Homochoric vector: (h_1, h_2, h_3)."""
         return Rotation.qu2ho(self.quaternion)
 
@@ -234,7 +241,7 @@ class Rotation:
         """Cubochoric vector: (c_1, c_2, c_3)."""
         return Rotation.qu2cu(self.quaternion)
 
-    def asM(self):
+    def M(self): # ToDo not sure about the name: as_M or M? we do not have a from_M
         """
         Intermediate representation supporting quaternion averaging.
 
@@ -244,114 +251,133 @@ class Rotation:
         https://doi.org/10.2514/1.28949
 
         """
-        return np.outer(self.quaternion,self.quaternion)
+        return np.einsum('...i,...j',self.quaternion,self.quaternion)
 
+    # for compatibility (old names do not follow convention)
+    asM          = M
+    asQuaternion = as_quaternion
+    asEulers     = as_Eulers
+    asAxisAngle  = as_axis_angle
+    asMatrix     = as_matrix
+    asRodrigues  = as_Rodrigues
+    asHomochoric = as_homochoric
 
     ################################################################################################
-    # static constructors. The input data needs to follow the convention, options allow to
+    # Static constructors. The input data needs to follow the conventions, options allow to
-    # relax these convections
+    # relax the conventions.
     @staticmethod
-    def fromQuaternion(quaternion,
+    def from_quaternion(quaternion,
-                       acceptHomomorph = False,
+                        acceptHomomorph = False,
-                       P = -1):
+                        P = -1):
 
-        qu =   quaternion if isinstance(quaternion,np.ndarray) and quaternion.dtype == np.dtype(float) \
+        qu = np.array(quaternion,dtype=float)
-                          else np.array(quaternion,dtype=float)
+        if qu.shape[:-2:-1] != (4,):
-        if P > 0: qu[1:4] *= -1                                                                     # convert from P=1 to P=-1
+            raise ValueError('Invalid shape.')
-        if qu[0] < 0.0:
+
-            if acceptHomomorph:
+        if P > 0: qu[...,1:4] *= -1                                                                 # convert from P=1 to P=-1
-                qu *= -1.
+        if acceptHomomorph:
-            else:
+            qu[qu[...,0] < 0.0] *= -1
-                raise ValueError('Quaternion has negative first component: {}.'.format(qu[0]))
+        else:
-        if not np.isclose(np.linalg.norm(qu), 1.0):
+            if np.any(qu[...,0] < 0.0):
-            raise ValueError('Quaternion is not of unit length: {} {} {} {}.'.format(*qu))
+                raise ValueError('Quaternion with negative first (real) component.')
+        if not np.all(np.isclose(np.linalg.norm(qu,axis=-1), 1.0)):
+            raise ValueError('Quaternion is not of unit length.')
 
         return Rotation(qu)
 
     @staticmethod
-    def fromEulers(eulers,
+    def from_Eulers(eulers,
-                   degrees = False):
+                    degrees = False):
+
+        eu = np.array(eulers,dtype=float)
+        if eu.shape[:-2:-1] != (3,):
+            raise ValueError('Invalid shape.')
 
-        eu = eulers if isinstance(eulers, np.ndarray) and eulers.dtype == np.dtype(float) \
-                    else np.array(eulers,dtype=float)
         eu = np.radians(eu) if degrees else eu
-        if np.any(eu < 0.0) or np.any(eu > 2.0*np.pi) or eu[1] > np.pi:
+        if np.any(eu < 0.0) or np.any(eu > 2.0*np.pi) or np.any(eu[...,1] > np.pi):                 # ToDo: No separate check for PHI
-            raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π]: {} {} {}.'.format(*eu))
+            raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π].')
 
         return Rotation(Rotation.eu2qu(eu))
 
     @staticmethod
-    def fromAxisAngle(angleAxis,
+    def from_axis_angle(axis_angle,
-                      degrees = False,
+                        degrees = False,
-                      normalise = False,
+                        normalise = False,
-                      P = -1):
+                        P = -1):
 
-        ax = angleAxis if isinstance(angleAxis, np.ndarray) and angleAxis.dtype == np.dtype(float) \
+        ax = np.array(axis_angle,dtype=float)
-                       else np.array(angleAxis,dtype=float)
+        if ax.shape[:-2:-1] != (4,):
-        if P > 0:     ax[0:3] *= -1                                                                 # convert from P=1 to P=-1
+            raise ValueError('Invalid shape.')
-        if degrees:   ax[  3]  = np.radians(ax[3])
+
-        if normalise: ax[0:3] /= np.linalg.norm(ax[0:3])
+        if P > 0:     ax[...,0:3] *= -1                                                             # convert from P=1 to P=-1
-        if ax[3] < 0.0 or ax[3] > np.pi:
+        if degrees:   ax[...,  3]  = np.radians(ax[...,3])
-            raise ValueError('Axis angle rotation angle outside of [0..π]: {}.'.format(ax[3]))
+        if normalise: ax[...,0:3] /= np.linalg.norm(ax[...,0:3],axis=-1)
-        if not np.isclose(np.linalg.norm(ax[0:3]), 1.0):
+        if np.any(ax[...,3] < 0.0) or np.any(ax[...,3] > np.pi):
-            raise ValueError('Axis angle rotation axis is not of unit length: {} {} {}.'.format(*ax[0:3]))
+            raise ValueError('Axis angle rotation angle outside of [0..π].')
+        if not np.all(np.isclose(np.linalg.norm(ax[...,0:3],axis=-1), 1.0)):
+            raise ValueError('Axis angle rotation axis is not of unit length.')
 
         return Rotation(Rotation.ax2qu(ax))
 
     @staticmethod
-    def fromBasis(basis,
+    def from_basis(basis,
-                  orthonormal = True,
+                   orthonormal = True,
-                  reciprocal = False,
+                   reciprocal = False):
-                 ):
+
+        om = np.array(basis,dtype=float)
+        if om.shape[:-3:-1] != (3,3):
+            raise ValueError('Invalid shape.')
 
-        om = basis if isinstance(basis, np.ndarray) else np.array(basis).reshape(3,3)
         if reciprocal:
-            om = np.linalg.inv(om.T/np.pi)                                                          # transform reciprocal basis set
+            om = np.linalg.inv(mechanics.transpose(om)/np.pi)                                       # transform reciprocal basis set
             orthonormal = False                                                                     # contains stretch
         if not orthonormal:
             (U,S,Vh) = np.linalg.svd(om)                                                            # singular value decomposition
-            om = np.dot(U,Vh)
+            om = np.einsum('...ij,...jl->...il',U,Vh)
-        if not np.isclose(np.linalg.det(om),1.0):
+        if not np.all(np.isclose(np.linalg.det(om),1.0)):
-            raise ValueError('matrix is not a proper rotation: {}.'.format(om))
+            raise ValueError('Orientation matrix has determinant ≠ 1.')
-        if    not np.isclose(np.dot(om[0],om[1]), 0.0) \
+        if    not np.all(np.isclose(np.einsum('...i,...i',om[...,0],om[...,1]), 0.0)) \
-           or not np.isclose(np.dot(om[1],om[2]), 0.0) \
+           or not np.all(np.isclose(np.einsum('...i,...i',om[...,1],om[...,2]), 0.0)) \
-           or not np.isclose(np.dot(om[2],om[0]), 0.0):
+           or not np.all(np.isclose(np.einsum('...i,...i',om[...,2],om[...,0]), 0.0)):
-            raise ValueError('matrix is not orthogonal: {}.'.format(om))
+            raise ValueError('Orientation matrix is not orthogonal.')
 
         return Rotation(Rotation.om2qu(om))
 
     @staticmethod
-    def fromMatrix(om,
+    def from_matrix(om):
-                  ):
 
-        return Rotation.fromBasis(om)
+        return Rotation.from_basis(om)
 
     @staticmethod
-    def fromRodrigues(rodrigues,
+    def from_Rodrigues(rodrigues,
-                      normalise = False,
+                       normalise = False,
-                      P = -1):
+                       P = -1):
 
-        ro = rodrigues if isinstance(rodrigues, np.ndarray) and rodrigues.dtype == np.dtype(float) \
+        ro = np.array(rodrigues,dtype=float)
-                       else np.array(rodrigues,dtype=float)
+        if ro.shape[:-2:-1] != (4,):
-        if P > 0:     ro[0:3] *= -1                                                                 # convert from P=1 to P=-1
+            raise ValueError('Invalid shape.')
-        if normalise: ro[0:3] /= np.linalg.norm(ro[0:3])
+
-        if not np.isclose(np.linalg.norm(ro[0:3]), 1.0):
+        if P > 0:     ro[...,0:3] *= -1                                                             # convert from P=1 to P=-1
-            raise ValueError('Rodrigues rotation axis is not of unit length: {} {} {}.'.format(*ro[0:3]))
+        if normalise: ro[...,0:3] /= np.linalg.norm(ro[...,0:3],axis=-1)
-        if ro[3] < 0.0:
+        if np.any(ro[...,3] < 0.0):
-            raise ValueError('Rodrigues rotation angle not positive: {}.'.format(ro[3]))
+            raise ValueError('Rodrigues vector rotation angle not positive.')
+        if not np.all(np.isclose(np.linalg.norm(ro[...,0:3],axis=-1), 1.0)):
+            raise ValueError('Rodrigues vector rotation axis is not of unit length.')
 
         return Rotation(Rotation.ro2qu(ro))
 
     @staticmethod
-    def fromHomochoric(homochoric,
+    def from_homochoric(homochoric,
-                       P = -1):
+                        P = -1):
+
+        ho = np.array(homochoric,dtype=float)
+        if ho.shape[:-2:-1] != (3,):
+            raise ValueError('Invalid shape.')
 
-        ho = homochoric if isinstance(homochoric, np.ndarray) and homochoric.dtype == np.dtype(float) \
-                        else np.array(homochoric,dtype=float)
         if P > 0: ho *= -1                                                                          # convert from P=1 to P=-1
 
-        if np.linalg.norm(ho) > (3.*np.pi/4.)**(1./3.)+1e-9:
+        if np.any(np.linalg.norm(ho,axis=-1) > (3.*np.pi/4.)**(1./3.)+1e-9):
-            raise ValueError('Coordinate outside of the sphere: {} {} {}.'.format(ho))
+            raise ValueError('Homochoric coordinate outside of the sphere.')
 
         return Rotation(Rotation.ho2qu(ho))
 
@@ -359,11 +385,12 @@ class Rotation:
     def fromCubochoric(cubochoric,
                        P = -1):
 
-        cu = cubochoric if isinstance(cubochoric, np.ndarray) and cubochoric.dtype == np.dtype(float) \
+        cu = np.array(cubochoric,dtype=float)
-                        else np.array(cubochoric,dtype=float)
+        if cu.shape[:-2:-1] != (3,):
+            raise ValueError('Invalid shape.')
 
         if np.abs(np.max(cu))>np.pi**(2./3.) * 0.5+1e-9:
-            raise ValueError('Coordinate outside of the cube: {} {} {}.'.format(*cu))
+            raise ValueError('Cubochoric coordinate outside of the cube: {} {} {}.'.format(*cu))
 
         ho = Rotation.cu2ho(cu)
         if P > 0: ho *= -1                                                                          # convert from P=1 to P=-1
@@ -403,17 +430,34 @@ class Rotation:
 
         return Rotation.fromQuaternion(np.real(vec.T[eig.argmax()]),acceptHomomorph = True)
 
-
     @staticmethod
-    def fromRandom():
+    def from_random(shape=None):
-        r = np.random.random(3)
+        if shape is None:
-        A = np.sqrt(r[2])
+            r = np.random.random(3)
-        B = np.sqrt(1.0-r[2])
+        elif hasattr(shape, '__iter__'):
-        return Rotation(np.array([np.cos(2.0*np.pi*r[0])*A,
+            r = np.random.random(tuple(shape)+(3,))
-                                  np.sin(2.0*np.pi*r[1])*B,
+        else:
-                                  np.cos(2.0*np.pi*r[1])*B,
+            r = np.random.random((shape,3))
-                                  np.sin(2.0*np.pi*r[0])*A])).standardize()
 
+        A = np.sqrt(r[...,2])
+        B = np.sqrt(1.0-r[...,2])
+        q = np.stack([np.cos(2.0*np.pi*r[...,0])*A,
+                      np.sin(2.0*np.pi*r[...,1])*B,
+                      np.cos(2.0*np.pi*r[...,1])*B,
+                      np.sin(2.0*np.pi*r[...,0])*A],axis=-1)
+
+        return Rotation(q.reshape(r.shape[:-1]+(4,)) if shape is not None else q).standardize()
+
+
+    # for compatibility (old names do not follow convention)
+    fromQuaternion = from_quaternion
+    fromEulers     = from_Eulers
+    fromAxisAngle  = from_axis_angle
+    fromBasis      = from_basis
+    fromMatrix     = from_matrix
+    fromRodrigues  = from_Rodrigues
+    fromHomochoric = from_homochoric
+    fromRandom     = from_random
 
 ####################################################################################################
 # Code below available according to the following conditions on https://github.com/MarDiehl/3Drotations
@@ -808,12 +852,11 @@ class Rotation:
                 c = np.cos(ax[3]*0.5)
                 s = np.sin(ax[3]*0.5)
                 qu = np.array([ c, ax[0]*s, ax[1]*s, ax[2]*s ])
-            return qu
         else:
             c = np.cos(ax[...,3:4]*.5)
             s = np.sin(ax[...,3:4]*.5)
             qu = np.where(np.abs(ax[...,3:4])<1.e-6,[1.0, 0.0, 0.0, 0.0],np.block([c, ax[...,:3]*s]))
-            return qu
+        return qu
 
     @staticmethod
     def ax2om(ax):
@@ -859,7 +902,7 @@ class Rotation:
                 # 180 degree case
                 ro += [np.inf] if np.isclose(ax[3],np.pi,atol=1.0e-15,rtol=0.0) else \
                       [np.tan(ax[3]*0.5)]
-            return np.array(ro)
+            ro = np.array(ro)
         else:
             ro = np.block([ax[...,:3],
                            np.where(np.isclose(ax[...,3:4],np.pi,atol=1.e-15,rtol=.0),
@@ -867,7 +910,7 @@ class Rotation:
                                     np.tan(ax[...,3:4]*0.5))
                           ])
             ro[np.abs(ax[...,3])<1.e-6] = [.0,.0,_P,.0]
-            return ro
+        return ro
 
     @staticmethod
     def ax2ho(ax):
@@ -875,11 +918,10 @@ class Rotation:
         if len(ax.shape) == 1:
             f = (0.75 * ( ax[3] - np.sin(ax[3]) ))**(1.0/3.0)
             ho = ax[0:3] * f
-            return ho
         else:
             f = (0.75 * ( ax[...,3:4] - np.sin(ax[...,3:4]) ))**(1.0/3.0)
             ho = ax[...,:3] * f
-            return ho
+        return ho
 
     @staticmethod
     def ax2cu(ax):
@@ -936,7 +978,6 @@ class Rotation:
             f = np.where(np.isfinite(ro[...,3:4]),2.0*np.arctan(ro[...,3:4]) -np.sin(2.0*np.arctan(ro[...,3:4])),np.pi)
             ho = np.where(np.broadcast_to(np.sum(ro[...,0:3]**2.0,axis=-1,keepdims=True) < 1.e-6,ro[...,0:3].shape),
                           np.zeros(3), ro[...,0:3]* (0.75*f)**(1.0/3.0))
-
         return ho
 
     @staticmethod
@@ -1010,7 +1051,7 @@ class Rotation:
         if len(ho.shape) == 1:
             return ball_to_cube(ho)
         else:
-            raise NotImplementedError
+            raise NotImplementedError('Support for multiple rotations missing')
 
 
     #---------- Cubochoric ----------
@@ -1045,4 +1086,4 @@ class Rotation:
         if len(cu.shape) == 1:
             return cube_to_ball(cu)
         else:
-            raise NotImplementedError
+            raise NotImplementedError('Support for multiple rotations missing')

python/damask/mechanics.py

@@ -135,16 +135,16 @@ def PK2(P,F):
 
     Parameters
     ----------
-    P : numpy.ndarray of shape (:,3,3) or (3,3)
+    P : numpy.ndarray of shape (...,3,3) or (3,3)
         First Piola-Kirchhoff stress.
-    F : numpy.ndarray of shape (:,3,3) or (3,3)
+    F : numpy.ndarray of shape (...,3,3) or (3,3)
         Deformation gradient.
 
     """
     if _np.shape(F) == _np.shape(P) == (3,3):
         S = _np.dot(_np.linalg.inv(F),P)
     else:
-        S = _np.einsum('ijk,ikl->ijl',_np.linalg.inv(F),P)
+        S = _np.einsum('...jk,...kl->...jl',_np.linalg.inv(F),P)
     return symmetric(S)
 
 
@@ -241,7 +241,7 @@ def symmetric(T):
 
     Parameters
     ----------
-    T : numpy.ndarray of shape (:,3,3) or (3,3)
+    T : numpy.ndarray of shape (...,3,3) or (3,3)
         Tensor of which the symmetrized values are computed.
 
     """
@@ -254,12 +254,12 @@ def transpose(T):
 
     Parameters
     ----------
-    T : numpy.ndarray of shape (:,3,3) or (3,3)
+    T : numpy.ndarray of shape (...,3,3) or (3,3)
         Tensor of which the transpose is computed.
 
     """
     return T.T if _np.shape(T) == (3,3) else \
-           _np.transpose(T,(0,2,1))
+           _np.swapaxes(T,axis2=-2,axis1=-1)
 
 
 def _polar_decomposition(T,requested):

python/tests/test_Rotation.py

@@ -157,6 +157,30 @@ class TestRotation:
             print(m,o,rot.asQuaternion())
             assert ok and o.max() < np.pi**(2./3.)*0.5+1.e-9
 
+    @pytest.mark.parametrize('function',[Rotation.from_quaternion,
+                                         Rotation.from_Eulers,
+                                         Rotation.from_axis_angle,
+                                         Rotation.from_matrix,
+                                         Rotation.from_Rodrigues,
+                                         Rotation.from_homochoric])
+    def test_invalid_shape(self,function):
+        invalid_shape = np.random.random(np.random.randint(8,32,(3)))
+        with pytest.raises(ValueError):
+            function(invalid_shape)
+
+    @pytest.mark.parametrize('function,invalid',[(Rotation.from_quaternion, np.array([-1,0,0,0])),
+                                                 (Rotation.from_quaternion, np.array([1,1,1,0])),
+                                                 (Rotation.from_Eulers,     np.array([1,4,0])),
+                                                 (Rotation.from_axis_angle, np.array([1,0,0,4])),
+                                                 (Rotation.from_axis_angle, np.array([1,1,0,1])),
+                                                 (Rotation.from_matrix,     np.random.rand(3,3)),
+                                                 (Rotation.from_Rodrigues,  np.array([1,0,0,-1])),
+                                                 (Rotation.from_Rodrigues,  np.array([1,1,0,1])),
+                                                 (Rotation.from_homochoric, np.array([2,2,2]))  ])
+    def test_invalid(self,function,invalid):
+        with pytest.raises(ValueError):
+            function(invalid)
+
     @pytest.mark.parametrize('conversion',[Rotation.qu2om,
                                            Rotation.qu2eu,
                                            Rotation.qu2ax,

python/tests/test_mechanics.py

@@ -1,4 +1,6 @@
+import pytest
 import numpy as np
+
 from damask import mechanics
 
 class TestMechanics:
@@ -7,127 +9,77 @@ class TestMechanics:
     c = np.random.randint(n)
 
 
-    def test_vectorize_Cauchy(self):
+    @pytest.mark.parametrize('function',[mechanics.deviatoric_part,
-        P = np.random.random((self.n,3,3))
+                                         mechanics.eigenvalues,
-        F = np.random.random((self.n,3,3))
+                                         mechanics.eigenvectors,
-        assert np.allclose(mechanics.Cauchy(P,F)[self.c],
+                                         mechanics.left_stretch,
-                           mechanics.Cauchy(P[self.c],F[self.c]))
+                                         mechanics.maximum_shear,
+                                         mechanics.Mises_strain,
+                                         mechanics.Mises_stress,
+                                         mechanics.right_stretch,
+                                         mechanics.rotational_part,
+                                         mechanics.spherical_part,
+                                         mechanics.symmetric,
+                                         mechanics.transpose,
+                                        ])
+    def test_vectorize_1_arg(self,function):
+        epsilon = np.random.rand(self.n,3,3)
+        assert np.allclose(function(epsilon)[self.c],function(epsilon[self.c]))
 
-    def test_vectorize_deviatoric_part(self):
+    @pytest.mark.parametrize('function',[mechanics.Cauchy,
-        x = np.random.random((self.n,3,3))
+                                         mechanics.PK2,
-        assert np.allclose(mechanics.deviatoric_part(x)[self.c],
+                                        ])
-                           mechanics.deviatoric_part(x[self.c]))
+    def test_vectorize_2_arg(self,function):
-
+        P = np.random.rand(self.n,3,3)
-    def test_vectorize_eigenvalues(self):
+        F = np.random.rand(self.n,3,3)
-        x = np.random.random((self.n,3,3))
+        assert np.allclose(function(P,F)[self.c],function(P[self.c],F[self.c]))
-        assert np.allclose(mechanics.eigenvalues(x)[self.c],
-                           mechanics.eigenvalues(x[self.c]))
-
-    def test_vectorize_eigenvectors(self):
-        x = np.random.random((self.n,3,3))
-        assert np.allclose(mechanics.eigenvectors(x)[self.c],
-                           mechanics.eigenvectors(x[self.c]))
-
-    def test_vectorize_left_stretch(self):
-        x = np.random.random((self.n,3,3))
-        assert np.allclose(mechanics.left_stretch(x)[self.c],
-                           mechanics.left_stretch(x[self.c]))
-
-    def test_vectorize_maximum_shear(self):
-        x = np.random.random((self.n,3,3))
-        assert np.allclose(mechanics.maximum_shear(x)[self.c],
-                           mechanics.maximum_shear(x[self.c]))
-
-    def test_vectorize_Mises_strain(self):
-        epsilon = np.random.random((self.n,3,3))
-        assert np.allclose(mechanics.Mises_strain(epsilon)[self.c],
-                           mechanics.Mises_strain(epsilon[self.c]))
-
-    def test_vectorize_Mises_stress(self):
-        sigma = np.random.random((self.n,3,3))
-        assert np.allclose(mechanics.Mises_stress(sigma)[self.c],
-                           mechanics.Mises_stress(sigma[self.c]))
-
-    def test_vectorize_PK2(self):
-        F = np.random.random((self.n,3,3))
-        P = np.random.random((self.n,3,3))
-        assert np.allclose(mechanics.PK2(P,F)[self.c],
-                           mechanics.PK2(P[self.c],F[self.c]))
-
-    def test_vectorize_right_stretch(self):
-        x = np.random.random((self.n,3,3))
-        assert np.allclose(mechanics.right_stretch(x)[self.c],
-                           mechanics.right_stretch(x[self.c]))
-
-    def test_vectorize_rotational_part(self):
-        x = np.random.random((self.n,3,3))
-        assert np.allclose(mechanics.rotational_part(x)[self.c],
-                           mechanics.rotational_part(x[self.c]))
-
-    def test_vectorize_spherical_part(self):
-        x = np.random.random((self.n,3,3))
-        assert np.allclose(mechanics.spherical_part(x,True)[self.c],
-                           mechanics.spherical_part(x[self.c],True))
 
     def test_vectorize_strain_tensor(self):
-        F = np.random.random((self.n,3,3))
+        F = np.random.rand(self.n,3,3)
         t = ['V','U'][np.random.randint(0,2)]
         m = np.random.random()*10. -5.0
         assert np.allclose(mechanics.strain_tensor(F,t,m)[self.c],
                            mechanics.strain_tensor(F[self.c],t,m))
 
-    def test_vectorize_symmetric(self):
+    @pytest.mark.parametrize('function',[mechanics.Cauchy,
-        x = np.random.random((self.n,3,3))
+                                         mechanics.PK2,
-        assert np.allclose(mechanics.symmetric(x)[self.c],
+                                        ])
-                           mechanics.symmetric(x[self.c]))
+    def test_stress_measures(self,function):
-
+        """Ensure that all stress measures are equivalent for no deformation."""
-    def test_vectorize_transpose(self):
+        P = np.random.rand(self.n,3,3)
-        x = np.random.random((self.n,3,3))
+        assert np.allclose(function(P,np.broadcast_to(np.eye(3),(self.n,3,3))),mechanics.symmetric(P))
-        assert np.allclose(mechanics.transpose(x)[self.c],
-                           mechanics.transpose(x[self.c]))
-
-
-    def test_Cauchy(self):
-        """Ensure Cauchy stress is symmetrized 1. Piola-Kirchhoff stress for no deformation."""
-        P = np.random.random((self.n,3,3))
-        assert np.allclose(mechanics.Cauchy(P,np.broadcast_to(np.eye(3),(self.n,3,3))),
-                           mechanics.symmetric(P))
 
+    def test_deviatoric_part(self):
+        I_n = np.broadcast_to(np.eye(3),(self.n,3,3))
+        r   = np.logical_not(I_n)*np.random.rand(self.n,3,3)
+        assert np.allclose(mechanics.deviatoric_part(I_n+r),r)
 
     def test_polar_decomposition(self):
         """F = RU = VR."""
-        F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.random((self.n,3,3))
+        F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.rand(self.n,3,3)
         R = mechanics.rotational_part(F)
         V = mechanics.left_stretch(F)
         U = mechanics.right_stretch(F)
         assert np.allclose(np.matmul(R,U),
                            np.matmul(V,R))
 
-
-    def test_PK2(self):
-        """Ensure 2. Piola-Kirchhoff stress is symmetrized 1. Piola-Kirchhoff stress for no deformation."""
-        P = np.random.random((self.n,3,3))
-        assert np.allclose(mechanics.PK2(P,np.broadcast_to(np.eye(3),(self.n,3,3))),
-                           mechanics.symmetric(P))
-
-
     def test_strain_tensor_no_rotation(self):
         """Ensure that left and right stretch give same results for no rotation."""
-        F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.random((self.n,3,3))
+        F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.rand(self.n,3,3)
         m = np.random.random()*20.0-10.0
         assert np.allclose(mechanics.strain_tensor(F,'U',m),
                            mechanics.strain_tensor(F,'V',m))
 
     def test_strain_tensor_rotation_equivalence(self):
         """Ensure that left and right strain differ only by a rotation."""
-        F = np.broadcast_to(np.eye(3),[self.n,3,3]) + (np.random.random((self.n,3,3))*0.5 - 0.25)
+        F = np.broadcast_to(np.eye(3),[self.n,3,3]) + (np.random.rand(self.n,3,3)*0.5 - 0.25)
         m = np.random.random()*5.0-2.5
         assert np.allclose(np.linalg.det(mechanics.strain_tensor(F,'U',m)),
                            np.linalg.det(mechanics.strain_tensor(F,'V',m)))
 
     def test_strain_tensor_rotation(self):
         """Ensure that pure rotation results in no strain."""
-        F = mechanics.rotational_part(np.random.random((self.n,3,3)))
+        F = mechanics.rotational_part(np.random.rand(self.n,3,3))
         t = ['V','U'][np.random.randint(0,2)]
         m = np.random.random()*2.0 - 1.0
         assert np.allclose(mechanics.strain_tensor(F,t,m),
@@ -139,21 +91,20 @@ class TestMechanics:
 
         Should be +1, but random F might contain a reflection.
         """
-        x = np.random.random((self.n,3,3))
+        x = np.random.rand(self.n,3,3)
         assert np.allclose(np.abs(np.linalg.det(mechanics.rotational_part(x))),
                            1.0)
 
-
     def test_spherical_deviatoric_part(self):
         """Ensure that full tensor is sum of spherical and deviatoric part."""
-        x = np.random.random((self.n,3,3))
+        x = np.random.rand(self.n,3,3)
         sph = mechanics.spherical_part(x,True)
         assert np.allclose(sph + mechanics.deviatoric_part(x),
                            x)
 
     def test_deviatoric_Mises(self):
         """Ensure that Mises equivalent stress depends only on deviatoric part."""
-        x = np.random.random((self.n,3,3))
+        x = np.random.rand(self.n,3,3)
         full = mechanics.Mises_stress(x)
         dev  = mechanics.Mises_stress(mechanics.deviatoric_part(x))
         assert np.allclose(full,
@@ -161,7 +112,7 @@ class TestMechanics:
 
     def test_spherical_mapping(self):
         """Ensure that mapping to tensor is correct."""
-        x = np.random.random((self.n,3,3))
+        x = np.random.rand(self.n,3,3)
         tensor = mechanics.spherical_part(x,True)
         scalar = mechanics.spherical_part(x)
         assert np.allclose(np.linalg.det(tensor),
@@ -169,35 +120,32 @@ class TestMechanics:
 
     def test_spherical_Mises(self):
         """Ensure that Mises equivalent strrain of spherical strain is 0."""
-        x = np.random.random((self.n,3,3))
+        x = np.random.rand(self.n,3,3)
         sph = mechanics.spherical_part(x,True)
         assert np.allclose(mechanics.Mises_strain(sph),
                            0.0)
 
     def test_symmetric(self):
         """Ensure that a symmetric tensor is half of the sum of a tensor and its transpose."""
-        x = np.random.random((self.n,3,3))
+        x = np.random.rand(self.n,3,3)
         assert np.allclose(mechanics.symmetric(x)*2.0,
                            mechanics.transpose(x)+x)
 
-
     def test_transpose(self):
         """Ensure that a symmetric tensor equals its transpose."""
-        x = mechanics.symmetric(np.random.random((self.n,3,3)))
+        x = mechanics.symmetric(np.random.rand(self.n,3,3))
         assert np.allclose(mechanics.transpose(x),
                            x)
 
-
     def test_Mises(self):
         """Ensure that equivalent stress is 3/2 of equivalent strain."""
-        x = np.random.random((self.n,3,3))
+        x = np.random.rand(self.n,3,3)
         assert np.allclose(mechanics.Mises_stress(x)/mechanics.Mises_strain(x),
                            1.5)
 
-
     def test_eigenvalues(self):
         """Ensure that the characteristic polynomial can be solved."""
-        A = mechanics.symmetric(np.random.random((self.n,3,3)))
+        A = mechanics.symmetric(np.random.rand(self.n,3,3))
         lambd = mechanics.eigenvalues(A)
         s = np.random.randint(self.n)
         for i in range(3):
@@ -205,7 +153,7 @@ class TestMechanics:
 
     def test_eigenvalues_and_vectors(self):
         """Ensure that eigenvalues and -vectors are the solution to the characteristic polynomial."""
-        A = mechanics.symmetric(np.random.random((self.n,3,3)))
+        A = mechanics.symmetric(np.random.rand(self.n,3,3))
         lambd = mechanics.eigenvalues(A)
         x     = mechanics.eigenvectors(A)
         s = np.random.randint(self.n)
@@ -214,12 +162,12 @@ class TestMechanics:
 
     def test_eigenvectors_RHS(self):
         """Ensure that RHS coordinate system does only change sign of determinant."""
-        A = mechanics.symmetric(np.random.random((self.n,3,3)))
+        A = mechanics.symmetric(np.random.rand(self.n,3,3))
         LRHS = np.linalg.det(mechanics.eigenvectors(A,RHS=False))
         RHS  = np.linalg.det(mechanics.eigenvectors(A,RHS=True))
         assert np.allclose(np.abs(LRHS),RHS)
 
     def test_spherical_no_shear(self):
         """Ensure that sherical stress has max shear of 0.0."""
-        A = mechanics.spherical_part(mechanics.symmetric(np.random.random((self.n,3,3))),True)
+        A = mechanics.spherical_part(mechanics.symmetric(np.random.rand(self.n,3,3)),True)
         assert np.allclose(mechanics.maximum_shear(A),0.0)

src/CPFEM.f90

@@ -73,28 +73,24 @@ subroutine CPFEM_initAll(el,ip)
   integer(pInt), intent(in) :: el, &                                                                !< FE el number
                                ip                                                                   !< FE integration point number
 
-  !$OMP CRITICAL(init)
+  CPFEM_init_done = .true.
-  if (.not. CPFEM_init_done) then
+  call DAMASK_interface_init
-    call DAMASK_interface_init
+  call prec_init
-    call prec_init
+  call IO_init
-    call IO_init
+  call numerics_init
-    call numerics_init
+  call debug_init
-    call debug_init
+  call config_init
-    call config_init
+  call math_init
-    call math_init
+  call rotations_init
-    call rotations_init
+  call HDF5_utilities_init
-    call HDF5_utilities_init
+  call results_init
-    call results_init
+  call discretization_marc_init(ip, el)
-    call discretization_marc_init(ip, el)
+  call lattice_init
-    call lattice_init
+  call material_init
-    call material_init
+  call constitutive_init
-    call constitutive_init
+  call crystallite_init
-    call crystallite_init
+  call homogenization_init
-    call homogenization_init
+  call CPFEM_init
-    call CPFEM_init
-    CPFEM_init_done = .true.
-  endif
-  !$OMP END CRITICAL(init)
 
 end subroutine CPFEM_initAll
 

src/DAMASK_marc.f90

@@ -261,11 +261,10 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
  endif
 
  !$ defaultNumThreadsInt = omp_get_num_threads()                                                    ! remember number of threads set by Marc
+ !$ call omp_set_num_threads(DAMASK_NumThreadsInt)                                                  ! set number of threads for parallel execution set by DAMASK_NUM_THREADS
 
  if (.not. CPFEM_init_done) call CPFEM_initAll(m(1),nn)
 
- !$ call omp_set_num_threads(DAMASK_NumThreadsInt)                                                  ! set number of threads for parallel execution set by DAMASK_NUM_THREADS
-
  computationMode = 0                                                                                ! save initialization value, since it does not result in any calculation
  if (lovl == 4 ) then                                                                               ! jacobian requested by marc
    if (timinc < theDelta .and. theInc == inc .and. lastLovl /= lovl) &                              ! first after cutback

src/crystallite.f90

@@ -303,6 +303,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
     s
   logical, dimension(homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: todo     !ToDo: need to set some values to false for different Ngrains
 
+    todo = .false.
 #ifdef DEBUG
   if (iand(debug_level(debug_crystallite),debug_levelSelective) /= 0 &
       .and. FEsolving_execElem(1) <= debug_e &