Merge branch 'development' into grid-filters

.gitattributes

@@ -3,8 +3,8 @@
 # always use LF, even if the files are edited on windows, they need to be compiled/used on unix
 * text eol=lf
 
-# Denote all files that are truly binary and should not be modified.
+# Denote all files that are binary and should not be modified.
 *.png binary
 *.jpg binary
-*.cae binary
 *.hdf5 binary
+*.pdf binary

.gitlab-ci.yml

@@ -115,13 +115,6 @@ Pytest:
     - release
 
 ###################################################################################################
-OrientationRelationship:
-  stage: preprocessing
-  script: OrientationRelationship/test.py
-  except:
-    - master
-    - release
-
 Pre_SeedGeneration:
   stage: preprocessing
   script: PreProcessing_SeedGeneration/test.py
@@ -506,18 +499,6 @@ GridSolver:
     - master
     - release
 
-Processing:
-  stage: createDocumentation
-  script:
-    - cd $DAMASKROOT/processing/pre
-    - $DAMASKROOT/PRIVATE/documenting/scriptHelpToWiki.py --debug *.py
-    - cd $DAMASKROOT/processing/post
-    - rm vtk2ang.py DAD*.py
-    - $DAMASKROOT/PRIVATE/documenting/scriptHelpToWiki.py --debug *.py
-  except:
-    - master
-    - release
-
 ##################################################################################################
 backupData:
   stage: saveDocumentation
@@ -528,7 +509,6 @@ backupData:
     - mv $TESTROOT/performance/time.png $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 b580fd4e992c2ed457f16d65bcba2e6d099dc29d
+Subproject commit 036faecca39b46fd2328597ca858cbb04e37f79a

VERSION

@@ -1 +1 @@
-v2.0.3-1177-ga41871e2
+v2.0.3-1294-g034367fa

examples/ConfigFiles/Homogenization_Thermal_Conduction.config

@@ -1,3 +1,3 @@
 thermal        conduction
-initialT       300.0
+t0             270.0
 (output) temperature

examples/SpectralMethod/EshelbyInclusion/material.config

@@ -6,7 +6,7 @@
 mech	none                               # isostrain 1 grain
 
 thermal                 adiabatic        # thermal strain (stress) induced mass transport
-initialT                300.0
+t0                      330.0
 (output)                temperature
 
 #-------------------#

processing/post/DADF5_postResults.py

@@ -47,6 +47,8 @@ for filename in options.filenames:
 
   coords = np.concatenate((z[:,:,:,None],y[:,:,:,None],x[:,:,:,None]),axis = 3) 
   
+  N_digits = int(np.floor(np.log10(int(results.increments[-1][3:]))))+1
+  N_digits = 5 # hack to keep test intact
   for i,inc in enumerate(results.iter_visible('increments')):
     print('Output step {}/{}'.format(i+1,len(results.increments)))
 
@@ -92,5 +94,6 @@ for filename in options.filenames:
     dirname  = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir))
     if not os.path.isdir(dirname):
       os.mkdir(dirname,0o755)
-    file_out = '{}_{}.txt'.format(os.path.splitext(os.path.split(filename)[-1])[0],inc)
+    file_out = '{}_inc{}.txt'.format(os.path.splitext(os.path.split(filename)[-1])[0],
+                                     inc[3:].zfill(N_digits))
     np.savetxt(os.path.join(dirname,file_out),data,header=header,comments='')

processing/post/DADF5_vtk_cells.py

@@ -2,6 +2,7 @@
 
 import os
 import argparse
+import re
 
 import h5py
 import numpy as np
@@ -66,7 +67,7 @@ for filename in options.filenames:
       for i in f['/geometry/T_c']:
         grid.InsertNextCell(vtk.VTK_HEXAHEDRON,8,i-1)
 
-
+  N_digits = int(np.floor(np.log10(int(results.increments[-1][3:]))))+1
   for i,inc in enumerate(results.iter_visible('increments')):
     print('Output step {}/{}'.format(i+1,len(results.increments)))
     vtk_data = []
@@ -89,10 +90,12 @@ for filename in options.filenames:
           x = results.get_dataset_location(label)
           if len(x) == 0:
             continue
+          ph_name = re.compile(r'(?<=(constituent\/))(.*?)(?=(generic))') #looking for phase name in dataset name
           array = results.read_dataset(x,0)
           shape = [array.shape[0],np.product(array.shape[1:])]
           vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE))
-          vtk_data[-1].SetName('1_'+x[0].split('/',1)[1])
+          dset_name = '1_' + re.sub(ph_name,r'',x[0].split('/',1)[1])           #removing phase name from generic dataset
+          vtk_data[-1].SetName(dset_name)
           grid.GetCellData().AddArray(vtk_data[-1])
     
     results.set_visible('constituents',  False)
@@ -132,7 +135,9 @@ for filename in options.filenames:
     dirname  = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir))
     if not os.path.isdir(dirname):
       os.mkdir(dirname,0o755)
-    file_out = '{}_{}.{}'.format(os.path.splitext(os.path.split(filename)[-1])[0],inc,writer.GetDefaultFileExtension())
+    file_out = '{}_inc{}.{}'.format(os.path.splitext(os.path.split(filename)[-1])[0],
+                                    inc[3:].zfill(N_digits),
+                                    writer.GetDefaultFileExtension())
     
     writer.SetCompressorTypeToZLib()
     writer.SetDataModeToBinary()

processing/post/DADF5_vtk_points.py

@@ -2,6 +2,7 @@
 
 import os
 import argparse
+import re
 
 import numpy as np
 import vtk
@@ -52,6 +53,7 @@ for filename in options.filenames:
   Polydata.SetVerts(Vertices)
   Polydata.Modified()
   
+  N_digits = int(np.floor(np.log10(int(results.increments[-1][3:]))))+1
   for i,inc in enumerate(results.iter_visible('increments')):
     print('Output step {}/{}'.format(i+1,len(results.increments)))
     vtk_data = []
@@ -75,10 +77,12 @@ for filename in options.filenames:
           x = results.get_dataset_location(label)
           if len(x) == 0:
             continue
+          ph_name = re.compile(r'(?<=(constituent\/))(.*?)(?=(generic))') #looking for phase name in dataset name
           array = results.read_dataset(x,0)
           shape = [array.shape[0],np.product(array.shape[1:])]
           vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True,array_type= vtk.VTK_DOUBLE))
-          vtk_data[-1].SetName('1_'+x[0].split('/',1)[1])
+          dset_name = '1_' + re.sub(ph_name,r'',x[0].split('/',1)[1])           #removing phase name from generic dataset
+          vtk_data[-1].SetName(dset_name)
           Polydata.GetCellData().AddArray(vtk_data[-1])
     
     results.set_visible('constituents',  False)
@@ -111,7 +115,9 @@ for filename in options.filenames:
     dirname  = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir))
     if not os.path.isdir(dirname):
       os.mkdir(dirname,0o755)
-    file_out = '{}_{}.{}'.format(os.path.splitext(os.path.split(filename)[-1])[0],inc,writer.GetDefaultFileExtension())
+    file_out = '{}_inc{}.{}'.format(os.path.splitext(os.path.split(filename)[-1])[0],
+                                   inc[3:].zfill(N_digits),
+                                   writer.GetDefaultFileExtension())
     
     writer.SetCompressorTypeToZLib()
     writer.SetDataModeToBinary()

processing/pre/geom_fromDREAM3D.py

@@ -145,7 +145,6 @@ for name in filenames:
   config_header += ['<microstructure>']
   for i in range(np.nanmax(microstructure)):
     config_header += ['[{}{}]'.format(label,i+1),
-                      'crystallite 1',
                       '(constituent)\tphase {}\ttexture {}\tfraction 1.0'.format(phase[i],i+1),
                       ]
 

processing/pre/geom_fromOsteonGeometry.py

@@ -126,15 +126,12 @@ for i in range(3,np.max(microstructure)):
 
 config_header = ['<microstructure>',
                  '[canal]',
-                 'crystallite 1',
                  '(constituent)\tphase 1\ttexture 1\tfraction 1.0',
                  '[interstitial]',
-                 'crystallite 1',
                  '(constituent)\tphase 2\ttexture 2\tfraction 1.0'
                 ]
 for i in range(3,np.max(microstructure)):
   config_header += ['[Point{}]'.format(i-2),
-                    'crystallite 1',
                     '(constituent)\tphase 3\ttexture {}\tfraction 1.0'.format(i)
                    ]
 

processing/pre/geom_fromVoronoiTessellation.py

@@ -290,7 +290,6 @@ for name in filenames:
     config_header += ['<microstructure>']
     for ID in grainIDs:
       config_header += ['[Grain{}]'.format(ID),
-                        'crystallite 1',
                         '(constituent)\tphase {}\ttexture {}\tfraction 1.0'.format(options.phase,ID)
                        ]
 

processing/pre/hybridIA_linODFsampling.py

@@ -19,7 +19,7 @@ def integerFactorization(i):
   return j
 
 def binAsBins(bin,intervals):
-  """Explode compound bin into 3D bins list"""
+  """Explode compound bin into 3D bins list."""
   bins = [0]*3
   bins[0] = (bin//(intervals[1] * intervals[2])) % intervals[0]
   bins[1] = (bin//intervals[2]) % intervals[1]
@@ -27,17 +27,17 @@ def binAsBins(bin,intervals):
   return bins
   
 def binsAsBin(bins,intervals):
-  """Implode 3D bins into compound bin"""
+  """Implode 3D bins into compound bin."""
   return (bins[0]*intervals[1] + bins[1])*intervals[2] + bins[2]
 
 def EulersAsBins(Eulers,intervals,deltas,center):
-  """Return list of Eulers translated into 3D bins list"""
+  """Return list of Eulers translated into 3D bins list."""
   return [int((euler+(0.5-center)*delta)//delta)%interval \
                   for euler,delta,interval in zip(Eulers,deltas,intervals) \
          ]
 
 def binAsEulers(bin,intervals,deltas,center):
-  """Compound bin number translated into list of Eulers"""
+  """Compound bin number translated into list of Eulers."""
   Eulers = [0.0]*3
   Eulers[2] = (bin%intervals[2] + center)*deltas[2]
   Eulers[1] = (bin//intervals[2]%intervals[1] + center)*deltas[1]
@@ -45,7 +45,7 @@ def binAsEulers(bin,intervals,deltas,center):
   return Eulers
 
 def directInvRepetitions(probability,scale):
-  """Calculate number of samples drawn by direct inversion"""
+  """Calculate number of samples drawn by direct inversion."""
   nDirectInv = 0
   for bin in range(len(probability)):                                                               # loop over bins
     nDirectInv += int(round(probability[bin]*scale))                                                # calc repetition
@@ -56,7 +56,7 @@ def directInvRepetitions(probability,scale):
 
 # ----- efficient algorithm ---------
 def directInversion (ODF,nSamples):
-  """ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians)"""
+  """ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians)."""
   nOptSamples = max(ODF['nNonZero'],nSamples)                                                       # random subsampling if too little samples requested
 
   nInvSamples = 0
@@ -118,7 +118,7 @@ def directInversion (ODF,nSamples):
 # ----- trial and error algorithms ---------
 
 def MonteCarloEulers (ODF,nSamples):
-  """ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians)"""
+  """ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians)."""
   countMC = 0
   maxdV_V = max(ODF['dV_V'])
   orientations     = np.zeros((nSamples,3),'f')
@@ -141,7 +141,7 @@ def MonteCarloEulers (ODF,nSamples):
 
 
 def MonteCarloBins (ODF,nSamples):
-  """ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians)"""
+  """ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians)."""
   countMC = 0
   maxdV_V = max(ODF['dV_V'])
   orientations     = np.zeros((nSamples,3),'f')
@@ -163,7 +163,7 @@ def MonteCarloBins (ODF,nSamples):
 
 
 def TothVanHoutteSTAT (ODF,nSamples):
-  """ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians)"""
+  """ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians)."""
   orientations     = np.zeros((nSamples,3),'f')
   reconstructedODF = np.zeros(ODF['nBins'],'f')
   unitInc = 1.0/nSamples
@@ -211,10 +211,6 @@ parser.add_option('-p','--phase',
                   dest = 'phase',
                   type = 'int', metavar = 'int',
                   help = 'phase index to be used [%default]')
-parser.add_option('--crystallite',
-                  dest = 'crystallite',
-                  type = 'int', metavar = 'int',
-                  help = 'crystallite index to be used [%default]')
 parser.add_option('-r', '--rnd',
                   dest = 'randomSeed',
                   type = 'int', metavar = 'int', \
@@ -223,7 +219,6 @@ parser.set_defaults(randomSeed = None,
                     number      = 500,
                     algorithm   = 'IA',
                     phase       = 1,
-                    crystallite = 1,
                     ang  = True,
                    )
 
@@ -240,7 +235,7 @@ if filenames == []: filenames = [None]
 for name in filenames:
   try:
     table = damask.ASCIItable(name = name, buffered = False, readonly=True)
-  except:
+  except IOError:
     continue
   damask.util.report(scriptName,name)
 
@@ -351,7 +346,6 @@ for name in filenames:
   
   for i,ID in enumerate(range(nSamples)):
     materialConfig += ['[Grain%s]'%(str(ID+1).zfill(formatwidth)),
-                      'crystallite %i'%options.crystallite,
                       '(constituent)   phase %i   texture %s   fraction 1.0'%(options.phase,str(ID+1).rjust(formatwidth)),
                      ]
   

processing/pre/patchFromReconstructedBoundaries.py

@@ -78,13 +78,11 @@ def rcbOrientationParser(content,idcolumn):
         damask.util.croak('You might not have chosen the correct column for the grain IDs! '+
                           'Please check the "--id" option.')
         raise
-      except:
-        raise
 
   return grains
 
 def rcbParser(content,M,size,tolerance,idcolumn,segmentcolumn):
-  """Parser for TSL-OIM reconstructed boundary files"""
+  """Parser for TSL-OIM reconstructed boundary files."""
 # find bounding box
   boxX = [1.*sys.maxint,-1.*sys.maxint]
   boxY = [1.*sys.maxint,-1.*sys.maxint]
@@ -99,8 +97,6 @@ def rcbParser(content,M,size,tolerance,idcolumn,segmentcolumn):
       damask.util.croak('You might not have chosen the correct column for the segment end points! '+
                         'Please check the "--segment" option.')
       raise
-    except:
-      raise
     (x[0],y[0]) = (M[0]*x[0]+M[1]*y[0],M[2]*x[0]+M[3]*y[0])   # apply transformation to coordinates
     (x[1],y[1]) = (M[0]*x[1]+M[1]*y[1],M[2]*x[1]+M[3]*y[1])   # to get rcb --> Euler system
     boxX[0] = min(boxX[0],x[0],x[1])
@@ -728,7 +724,7 @@ def image(name,imgsize,marginX,marginY,rcData):
 
 # -------------------------
 def inside(x,y,points):
-  """Tests whether point(x,y) is within polygon described by points"""
+  """Tests whether point(x,y) is within polygon described by points."""
   inside = False
   npoints=len(points)
   (x1,y1) = points[npoints-1]                                                # start with last point of points
@@ -750,7 +746,7 @@ def inside(x,y,points):
   
 # -------------------------
 def fftbuild(rcData,height,xframe,yframe,grid,extrusion):
-  """Build array of grain numbers"""
+  """Build array of grain numbers."""
   maxX = -1.*sys.maxint
   maxY = -1.*sys.maxint
   for line in rcData['point']:                                               # find data range
@@ -883,7 +879,7 @@ try:
   boundaryFile = open(args[0])
   boundarySegments = boundaryFile.readlines()
   boundaryFile.close()
-except:
+except IOError:
   damask.util.croak('unable to read boundary file "{}".'.format(args[0]))
   raise
 
@@ -941,19 +937,15 @@ if any(output in options.output for output in ['spectral','mentat']):
   
   for i,grain in enumerate(rcData['grainMapping']):
     config+=['[grain{}]'.format(grain),
-             'crystallite\t1',
              '(constituent)\tphase 1\ttexture {}\tfraction 1.0'.format(i+1)]
   if (options.xmargin > 0.0):
     config+=['[x-margin]',
-             'crystallite\t1',
              '(constituent)\tphase 2\ttexture {}\tfraction 1.0\n'.format(len(rcData['grainMapping'])+1)]
   if (options.ymargin > 0.0):
     config+=['[y-margin]',
-             'crystallite\t1',
              '(constituent)\tphase 2\ttexture {}\tfraction 1.0\n'.format(len(rcData['grainMapping'])+1)]
   if (options.xmargin > 0.0 and options.ymargin > 0.0):
     config+=['[xy-margin]',
-             'crystallite\t1',
              '(constituent)\tphase 2\ttexture {}\tfraction 1.0\n'.format(len(rcData['grainMapping'])+1)]
   
   if (options.xmargin > 0.0 or options.ymargin > 0.0):

python/damask/dadf5.py

@@ -1,7 +1,10 @@
 from queue import Queue
 import re
 import glob
+import os
 
+import vtk
+from vtk.util import numpy_support
 import h5py
 import numpy as np
 
@@ -30,7 +33,14 @@ class DADF5():
     """
     with h5py.File(fname,'r') as f:
       
-      if f.attrs['DADF5-major'] != 0 or not 2 <= f.attrs['DADF5-minor'] <= 3:
+      try:
+        self.version_major = f.attrs['DADF5_version_major']
+        self.version_minor = f.attrs['DADF5_version_minor']
+      except KeyError:
+        self.version_major = f.attrs['DADF5-major']
+        self.version_minor = f.attrs['DADF5-minor']
+
+      if self.version_major != 0 or not 2 <= self.version_minor <= 4:
         raise TypeError('Unsupported DADF5 version {} '.format(f.attrs['DADF5-version']))
     
       self.structured = 'grid' in f['geometry'].attrs.keys()
@@ -40,7 +50,8 @@ class DADF5():
         self.size = f['geometry'].attrs['size']
         
       r=re.compile('inc[0-9]+')
-      self.increments = [i for i in f.keys() if r.match(i)]
+      increments_unsorted = {int(i[3:]):i for i in f.keys() if r.match(i)}
+      self.increments     = [increments_unsorted[i] for i in sorted(increments_unsorted)]
       self.times          = [round(f[i].attrs['time/s'],12) for i in self.increments]
 
       self.Nmaterialpoints, self.Nconstituents =   np.shape(f['mapping/cellResults/constituent'])
@@ -165,6 +176,9 @@ class DADF5():
       end increment (included)
 
     """
+    if self.version_minor >= 4:
+      self.__manage_visible([    'inc{}'.format(i) for i in range(start,end+1)],'increments','set')
+    else:
       self.__manage_visible(['inc{:05d}'.format(i) for i in range(start,end+1)],'increments','set')
 
 
@@ -180,6 +194,9 @@ class DADF5():
       end increment (included)
 
     """
+    if self.version_minor >= 4:
+      self.__manage_visible([    'inc{}'.format(i) for i in range(start,end+1)],'increments','add')
+    else:
       self.__manage_visible(['inc{:05d}'.format(i) for i in range(start,end+1)],'increments','add')
 
 
@@ -195,6 +212,9 @@ class DADF5():
       end increment (included)
 
     """
+    if self.version_minor >= 4:
+      self.__manage_visible([    'inc{}'.format(i) for i in range(start,end+1)],'increments','del')
+    else:
       self.__manage_visible(['inc{:05d}'.format(i) for i in range(start,end+1)],'increments','del')
 
 
@@ -343,7 +363,7 @@ class DADF5():
           f[k]
           path.append(k)
         except KeyError as e:
-          print('unable to locate geometry dataset: {}'.format(str(e)))
+          pass
         for o,p in zip(['constituents','materialpoints'],['con_physics','mat_physics']):
           for oo in self.iter_visible(o):
             for pp in self.iter_visible(p):
@@ -352,7 +372,7 @@ class DADF5():
                 f[k]
                 path.append(k)
               except KeyError as e:
-                print('unable to locate {} dataset: {}'.format(o,str(e)))
+                pass
     return path
     
     
@@ -416,12 +436,82 @@ class DADF5():
                             np.linspace(delta[1],self.size[1]-delta[1],self.grid[1]),
                             np.linspace(delta[0],self.size[0]-delta[0],self.grid[0]),
                            )
-      return np.concatenate((x[:,:,:,None],y[:,:,:,None],y[:,:,:,None]),axis = 3).reshape([np.product(self.grid),3])
+      return np.concatenate((x[:,:,:,None],y[:,:,:,None],z[:,:,:,None]),axis = 3).reshape([np.product(self.grid),3])
     else:
       with h5py.File(self.fname,'r') as f:
         return f['geometry/x_c'][()]
     
     
+  def add_absolute(self,x):
+    """
+    Add absolute value.
+
+    Parameters
+    ----------
+    x : str
+      Label of the dataset containing a scalar, vector, or tensor.
+
+    """
+    def __add_absolute(x):
+
+      return {
+              'data':  np.abs(x['data']),
+              'label': '|{}|'.format(x['label']),
+              'meta':  {
+                        'Unit':        x['meta']['Unit'],
+                        'Description': 'Absolute value of {} ({})'.format(x['label'],x['meta']['Description']),
+                        'Creator':     'dadf5.py:add_abs v{}'.format(version)
+                        }
+               }
+
+    requested = [{'label':x,'arg':'x'}]
+    
+    self.__add_generic_pointwise(__add_absolute,requested)
+    
+
+  def add_calculation(self,formula,label,unit='n/a',description=None,vectorized=True):
+    """
+    Add result of a general formula.
+    
+    Parameters
+    ----------
+    formula : str
+      Formula, refer to datasets by ‘#Label#‘.
+    label : str
+      Label of the dataset containing the result of the calculation.
+    unit : str, optional
+      Physical unit of the result.
+    description : str, optional
+      Human readable description of the result.
+    vectorized : bool, optional
+      Indicate whether the formula is written in vectorized form. Default is ‘True’.
+
+    """
+    if vectorized is not True:
+      raise NotImplementedError
+    
+    def __add_calculation(**kwargs):
+      
+      formula = kwargs['formula']
+      for d in re.findall(r'#(.*?)#',formula):
+        formula = formula.replace('#{}#'.format(d),"kwargs['{}']['data']".format(d))
+        
+      return {
+              'data':  eval(formula), 
+              'label': kwargs['label'],
+              'meta':  {
+                        'Unit':        kwargs['unit'],
+                        'Description': '{} (formula: {})'.format(kwargs['description'],kwargs['formula']),
+                        'Creator':     'dadf5.py:add_calculation v{}'.format(version)
+                        }
+               }
+    
+    requested    = [{'label':d,'arg':d} for d in set(re.findall(r'#(.*?)#',formula))]               # datasets used in the formula
+    pass_through = {'formula':formula,'label':label,'unit':unit,'description':description} 
+    
+    self.__add_generic_pointwise(__add_calculation,requested,pass_through)
+
+
   def add_Cauchy(self,P='P',F='F'):
     """
     Add Cauchy stress calculated from 1. Piola-Kirchhoff stress and deformation gradient.
@@ -453,6 +543,90 @@ class DADF5():
     self.__add_generic_pointwise(__add_Cauchy,requested)
 
 
+  def add_determinant(self,x):
+    """
+    Add the determinant of a tensor.
+
+    Parameters
+    ----------
+    x : str
+      Label of the dataset containing a tensor.
+
+    """
+    def __add_determinant(x):
+      
+      return {
+              'data':  np.linalg.det(x['data']),
+              'label': 'det({})'.format(x['label']),
+              'meta':  {
+                        'Unit':        x['meta']['Unit'],
+                        'Description': 'Determinant of tensor {} ({})'.format(x['label'],x['meta']['Description']),
+                        'Creator':     'dadf5.py:add_determinant v{}'.format(version)
+                        }
+              }
+      
+    requested = [{'label':x,'arg':'x'}]
+    
+    self.__add_generic_pointwise(__add_determinant,requested)
+
+
+  def add_deviator(self,x):
+    """
+    Add the deviatoric part of a tensor.
+
+    Parameters
+    ----------
+    x : str
+      Label of the dataset containing a tensor.
+
+    """
+    def __add_deviator(x):
+      
+      if not np.all(np.array(x['data'].shape[1:]) == np.array([3,3])):
+        raise ValueError
+      
+      return {
+              'data':  mechanics.deviatoric_part(x['data']),
+              'label': 's_{}'.format(x['label']),
+              'meta':  {
+                        'Unit':        x['meta']['Unit'],
+                        'Description': 'Deviator of tensor {} ({})'.format(x['label'],x['meta']['Description']),
+                        'Creator':     'dadf5.py:add_deviator v{}'.format(version)
+                        }
+               }
+      
+    requested = [{'label':x,'arg':'x'}]
+    
+    self.__add_generic_pointwise(__add_deviator,requested)
+    
+    
+  def add_maximum_shear(self,x):
+    """
+    Add maximum shear components of symmetric tensor.
+    
+    Parameters
+    ----------
+    x : str
+      Label of the dataset containing a symmetric tensor.
+
+    """
+    def __add_maximum_shear(x):
+
+      return {
+              'data':  mechanics.maximum_shear(x['data']),
+              'label': 'max_shear({})'.format(x['label']),
+              'meta':  {
+                        'Unit':        x['meta']['Unit'],
+                        'Description': 'Maximum shear component of of {} ({})'.format(x['label'],x['meta']['Description']),
+                        'Creator':     'dadf5.py:add_maximum_shear v{}'.format(version)
+                        }
+               }
+
+    requested = [{'label':x,'arg':'x'}]
+
+    self.__add_generic_pointwise(__add_maximum_shear,requested)
+
+
   def add_Mises(self,x):
     """
     Add the equivalent Mises stress or strain of a symmetric tensor.
@@ -523,58 +697,33 @@ class DADF5():
     self.__add_generic_pointwise(__add_norm,requested,{'ord':ord})
     
     
-  def add_absolute(self,x):
+  def add_principal_components(self,x):
     """
-    Add absolute value.
+    Add principal components of symmetric tensor.
+    
+    The principal components are sorted in descending order, each repeated according to its multiplicity.
     
     Parameters
     ----------
     x : str
-      Label of the dataset containing a scalar, vector, or tensor.
+      Label of the dataset containing a symmetric tensor.
 
     """
-    def __add_absolute(x):
+    def __add_principal_components(x):
 
       return {
-              'data':  np.abs(x['data']),
+              'data':  mechanics.principal_components(x['data']),
-              'label': '|{}|'.format(x['label']),
+              'label': 'lambda_{}'.format(x['label']),
               'meta':  {
                         'Unit':        x['meta']['Unit'],
-                        'Description': 'Absolute value of {} ({})'.format(x['label'],x['meta']['Description']),
+                        'Description': 'Pricipal components of {} ({})'.format(x['label'],x['meta']['Description']),
-                        'Creator':     'dadf5.py:add_abs v{}'.format(version)
+                        'Creator':     'dadf5.py:add_principal_components v{}'.format(version)
                         }
                }
 
     requested = [{'label':x,'arg':'x'}]
 
-    self.__add_generic_pointwise(__add_absolute,requested)
+    self.__add_generic_pointwise(__add_principal_components,requested)
-
-
-  def add_determinant(self,x):
-    """
-    Add the determinant of a tensor.
-
-    Parameters
-    ----------
-    x : str
-      Label of the dataset containing a tensor.
-
-    """
-    def __add_determinant(x):
-      
-      return {
-              'data':  np.linalg.det(x['data']),
-              'label': 'det({})'.format(x['label']),
-              'meta':  {
-                        'Unit':        x['meta']['Unit'],
-                        'Description': 'Determinant of tensor {} ({})'.format(x['label'],x['meta']['Description']),
-                        'Creator':     'dadf5.py:add_determinant v{}'.format(version)
-                        }
-              }
-      
-    requested = [{'label':x,'arg':'x'}]
-    
-    self.__add_generic_pointwise(__add_determinant,requested)
 
 
   def add_spherical(self,x):
@@ -607,79 +756,6 @@ class DADF5():
     self.__add_generic_pointwise(__add_spherical,requested)
 
 
-  def add_deviator(self,x):
-    """
-    Add the deviatoric part of a tensor.
-
-    Parameters
-    ----------
-    x : str
-      Label of the dataset containing a tensor.
-
-    """
-    def __add_deviator(x):
-      
-      if not np.all(np.array(x['data'].shape[1:]) == np.array([3,3])):
-        raise ValueError
-      
-      return {
-              'data':  mechanics.deviatoric_part(x['data']),
-              'label': 's_{}'.format(x['label']),
-              'meta':  {
-                        'Unit':        x['meta']['Unit'],
-                        'Description': 'Deviator of tensor {} ({})'.format(x['label'],x['meta']['Description']),
-                        'Creator':     'dadf5.py:add_deviator v{}'.format(version)
-                        }
-               }
-      
-    requested = [{'label':x,'arg':'x'}]
-    
-    self.__add_generic_pointwise(__add_deviator,requested)
-    
-
-  def add_calculation(self,formula,label,unit='n/a',description=None,vectorized=True):
-    """
-    Add result of a general formula.
-    
-    Parameters
-    ----------
-    formula : str
-      Formula, refer to datasets by ‘#Label#‘.
-    label : str
-      Label of the dataset containing the result of the calculation.
-    unit : str, optional
-      Physical unit of the result.
-    description : str, optional
-      Human readable description of the result.
-    vectorized : bool, optional
-      Indicate whether the formula is written in vectorized form. Default is ‘True’.
-
-    """
-    if vectorized is not True:
-      raise NotImplementedError
-    
-    def __add_calculation(**kwargs):
-      
-      formula = kwargs['formula']
-      for d in re.findall(r'#(.*?)#',formula):
-        formula = formula.replace('#{}#'.format(d),"kwargs['{}']['data']".format(d))
-        
-      return {
-              'data':  eval(formula), 
-              'label': kwargs['label'],
-              'meta':  {
-                        'Unit':        kwargs['unit'],
-                        'Description': '{} (formula: {})'.format(kwargs['description'],kwargs['formula']),
-                        'Creator':     'dadf5.py:add_calculation v{}'.format(version)
-                        }
-               }
-    
-    requested    = [{'label':d,'arg':d} for d in set(re.findall(r'#(.*?)#',formula))]               # datasets used in the formula
-    pass_through = {'formula':formula,'label':label,'unit':unit,'description':description} 
-    
-    self.__add_generic_pointwise(__add_calculation,requested,pass_through)
-
-
   def add_strain_tensor(self,F='F',t='U',m=0):
     """
     Add strain tensor calculated from a deformation gradient.
@@ -714,62 +790,6 @@ class DADF5():
     self.__add_generic_pointwise(__add_strain_tensor,requested,{'t':t,'m':m})
 
 
-  def add_principal_components(self,x):
-    """
-    Add principal components of symmetric tensor.
-    
-    The principal components are sorted in descending order, each repeated according to its multiplicity.
-    
-    Parameters
-    ----------
-    x : str
-      Label of the dataset containing a symmetric tensor.
-
-    """
-    def __add_principal_components(x):
-
-      return {
-              'data':  mechanics.principal_components(x['data']),
-              'label': 'lambda_{}'.format(x['label']),
-              'meta':  {
-                        'Unit':        x['meta']['Unit'],
-                        'Description': 'Pricipal components of {} ({})'.format(x['label'],x['meta']['Description']),
-                        'Creator':     'dadf5.py:add_principal_components v{}'.format(version)
-                        }
-               }
-
-    requested = [{'label':x,'arg':'x'}]
-
-    self.__add_generic_pointwise(__add_principal_components,requested)
-    
-    
-  def add_maximum_shear(self,x):
-    """
-    Add maximum shear components of symmetric tensor.
-    
-    Parameters
-    ----------
-    x : str
-      Label of the dataset containing a symmetric tensor.
-
-    """
-    def __add_maximum_shear(x):
-
-      return {
-              'data':  mechanics.maximum_shear(x['data']),
-              'label': 'max_shear({})'.format(x['label']),
-              'meta':  {
-                        'Unit':        x['meta']['Unit'],
-                        'Description': 'Maximum shear component of of {} ({})'.format(x['label'],x['meta']['Description']),
-                        'Creator':     'dadf5.py:add_maximum_shear v{}'.format(version)
-                        }
-               }
-
-    requested = [{'label':x,'arg':'x'}]
-
-    self.__add_generic_pointwise(__add_maximum_shear,requested)
-
-
   def __add_generic_pointwise(self,func,datasets_requested,extra_args={}):
     """
     General function to add pointwise data.
@@ -824,3 +844,142 @@ class DADF5():
         N_added +=1
 
     pool.wait_completion()
+
+
+  def to_vtk(self,labels,mode='Cell'):
+    """
+    Export to vtk cell/point data.
+    
+    Parameters
+    ----------
+    labels : list of str
+      Labels of the datasets to be exported.
+    mode : str, either 'Cell' or 'Point'
+      Export in cell format or point format.
+      Default value is 'Cell'.
+
+    """
+    if mode=='Cell':
+
+      if self.structured:
+  
+        coordArray = [vtk.vtkDoubleArray(),vtk.vtkDoubleArray(),vtk.vtkDoubleArray()]
+        for dim in [0,1,2]:
+          for c in np.linspace(0,self.size[dim],1+self.grid[dim]):
+            coordArray[dim].InsertNextValue(c)
+  
+        vtk_geom = vtk.vtkRectilinearGrid()
+        vtk_geom.SetDimensions(*(self.grid+1))
+        vtk_geom.SetXCoordinates(coordArray[0])
+        vtk_geom.SetYCoordinates(coordArray[1])
+        vtk_geom.SetZCoordinates(coordArray[2])
+  
+      else:
+        
+        nodes = vtk.vtkPoints()
+        with h5py.File(self.fname) as f:
+          nodes.SetData(numpy_support.numpy_to_vtk(f['/geometry/x_n'][()],deep=True))
+          
+          vtk_geom = vtk.vtkUnstructuredGrid()
+          vtk_geom.SetPoints(nodes)
+          vtk_geom.Allocate(f['/geometry/T_c'].shape[0])
+          for i in f['/geometry/T_c']:
+            vtk_geom.InsertNextCell(vtk.VTK_HEXAHEDRON,8,i-1) # not for all elements!
+    elif mode == 'Point':
+      Points   = vtk.vtkPoints()
+      Vertices = vtk.vtkCellArray()  
+      for c in self.cell_coordinates():
+        pointID = Points.InsertNextPoint(c)
+        Vertices.InsertNextCell(1)
+        Vertices.InsertCellPoint(pointID)
+      
+      vtk_geom = vtk.vtkPolyData()
+      vtk_geom.SetPoints(Points)
+      vtk_geom.SetVerts(Vertices)
+      vtk_geom.Modified()
+  
+    N_digits = int(np.floor(np.log10(int(self.increments[-1][3:]))))+1
+    
+    for i,inc in enumerate(self.iter_visible('increments')):
+      vtk_data = []
+      
+      materialpoints_backup = self.visible['materialpoints'].copy()
+      self.set_visible('materialpoints',False)
+      for label in labels:
+        for p in self.iter_visible('con_physics'):
+          if p != 'generic':
+            for c in self.iter_visible('constituents'):
+              x = self.get_dataset_location(label)
+              if len(x) == 0:
+                continue
+              array = self.read_dataset(x,0)
+              shape = [array.shape[0],np.product(array.shape[1:])]
+              vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),
+                              deep=True,array_type= vtk.VTK_DOUBLE))
+              vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) #ToDo: hard coded 1!
+              vtk_geom.GetCellData().AddArray(vtk_data[-1])
+
+          else:
+            x = self.get_dataset_location(label)
+            if len(x) == 0:
+              continue
+            array = self.read_dataset(x,0)
+            shape = [array.shape[0],np.product(array.shape[1:])]
+            vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),
+                            deep=True,array_type= vtk.VTK_DOUBLE))
+            ph_name = re.compile(r'(?<=(constituent\/))(.*?)(?=(generic))')                         # identify  phase name
+            dset_name = '1_' + re.sub(ph_name,r'',x[0].split('/',1)[1])                             # removing phase name
+            vtk_data[-1].SetName(dset_name)
+            vtk_geom.GetCellData().AddArray(vtk_data[-1])
+
+      self.set_visible('materialpoints',materialpoints_backup)
+  
+      constituents_backup = self.visible['constituents'].copy()
+      self.set_visible('constituents',False)
+      for label in labels:
+        for p in self.iter_visible('mat_physics'):
+          if p != 'generic':
+            for m in self.iter_visible('materialpoints'):
+              x = self.get_dataset_location(label)
+              if len(x) == 0:
+                continue
+              array = self.read_dataset(x,0)
+              shape = [array.shape[0],np.product(array.shape[1:])]
+              vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),
+                              deep=True,array_type= vtk.VTK_DOUBLE))
+              vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) #ToDo: why 1_?
+              vtk_geom.GetCellData().AddArray(vtk_data[-1])
+          else:
+            x = self.get_dataset_location(label)
+            if len(x) == 0:
+              continue
+            array = self.read_dataset(x,0)
+            shape = [array.shape[0],np.product(array.shape[1:])]
+            vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),
+                            deep=True,array_type= vtk.VTK_DOUBLE))
+            vtk_data[-1].SetName('1_'+x[0].split('/',1)[1])
+            vtk_geom.GetCellData().AddArray(vtk_data[-1])
+      self.set_visible('constituents',constituents_backup)
+     
+      if mode=='Cell': 
+        writer = vtk.vtkXMLRectilinearGridWriter() if self.structured else \
+                 vtk.vtkXMLUnstructuredGridWriter()
+        x = self.get_dataset_location('u_n')
+        vtk_data.append(numpy_support.numpy_to_vtk(num_array=self.read_dataset(x,0),
+                        deep=True,array_type=vtk.VTK_DOUBLE))
+        vtk_data[-1].SetName('u')
+        vtk_geom.GetPointData().AddArray(vtk_data[-1])
+      elif mode == 'Point':
+        writer = vtk.vtkXMLPolyDataWriter()
+
+      
+      file_out = '{}_inc{}.{}'.format(os.path.splitext(os.path.basename(self.fname))[0],
+                                      inc[3:].zfill(N_digits),
+                                      writer.GetDefaultFileExtension())
+      
+      writer.SetCompressorTypeToZLib()
+      writer.SetDataModeToBinary()
+      writer.SetFileName(file_out)
+      writer.SetInputData(vtk_geom)
+
+      writer.Write()

python/damask/orientation.py

@@ -701,14 +701,14 @@ class Symmetry:
 
     v = np.array(vector,dtype=float)
     if proper:                                                                                      # check both improper ...
-      theComponents = np.dot(basis['improper'],v)
+      theComponents = np.around(np.dot(basis['improper'],v),12)
       inSST = np.all(theComponents >= 0.0)
       if not inSST:                                                                                 # ... and proper SST
-        theComponents = np.dot(basis['proper'],v)
+        theComponents = np.around(np.dot(basis['proper'],v),12)
         inSST = np.all(theComponents >= 0.0)
     else:      
       v[2] = abs(v[2])                                                                              # z component projects identical 
-      theComponents = np.dot(basis['improper'],v)                                                   # for positive and negative values
+      theComponents = np.around(np.dot(basis['improper'],v),12)                                     # for positive and negative values
       inSST = np.all(theComponents >= 0.0)
 
     if color:                                                                                       # have to return color array
@@ -875,7 +875,7 @@ class Lattice:
       [[ 17, 12,  5],[ 17,  7, 17]],
       [[  5, 17, 12],[ 17, 17,  7]],
       [[ 12, -5,-17],[  7,-17,-17]],
-      [[-17,-12,  5],[-17,  7, 17]]],dtype='float')}
+      [[-17,-12,  5],[-17,-7, 17]]],dtype='float')}
  
   # Greninger--Troiano' orientation relationship for fcc <-> bcc transformation
   # from Y. He et al., Journal of Applied Crystallography 39:72-81, 2006 
@@ -901,7 +901,7 @@ class Lattice:
       [[-17,-17,  7],[-17, -5, 12]],
       [[  7,-17,-17],[ 12,-17, -5]],
       [[ 17, -7,-17],[ 5, -12,-17]],
-      [[ 17,-17,  7],[ 17, -5,-12]],
+      [[ 17,-17, -7],[ 17, -5,-12]],
       [[ -7, 17,-17],[-12, 17, -5]],
       [[-17,  7,-17],[ -5, 12,-17]],
       [[-17, 17, -7],[-17,  5,-12]]],dtype='float'),
@@ -957,7 +957,7 @@ class Lattice:
       [[  2,  1, -1],[  0, -1,  1]],
       [[ -1, -2, -1],[  0, -1,  1]],
       [[ -1,  1,  2],[  0, -1,  1]],
-      [[ -1,  2,  1],[  0, -1,  1]],
+      [[  2, -1,  1],[  0, -1,  1]], #It is wrong in the paper, but matrix is correct
       [[ -1,  2,  1],[  0, -1,  1]],
       [[ -1, -1, -2],[  0, -1,  1]]],dtype='float')}
                
@@ -1025,7 +1025,7 @@ class Lattice:
     https://doi.org/10.1016/j.actamat.2004.11.021
 
     """
-    models={'KS':self.KS, 'GT':self.GT, "GT'":self.GTprime, 
+    models={'KS':self.KS, 'GT':self.GT, 'GT_prime':self.GTprime, 
             'NW':self.NW, 'Pitsch': self.Pitsch, 'Bain':self.Bain}
     try:
       relationship = models[model]
@@ -1046,13 +1046,13 @@ class Lattice:
     for miller in np.hstack((relationship['planes'],relationship['directions'])):
       myPlane     = miller[myPlane_id]/    np.linalg.norm(miller[myPlane_id])
       myDir       = miller[myDir_id]/      np.linalg.norm(miller[myDir_id])
-      myMatrix    = np.array([myDir,np.cross(myPlane,myDir),myPlane]).T
+      myMatrix    = np.array([myDir,np.cross(myPlane,myDir),myPlane])
 
       otherPlane  = miller[otherPlane_id]/ np.linalg.norm(miller[otherPlane_id])
       otherDir    = miller[otherDir_id]/   np.linalg.norm(miller[otherDir_id])
-      otherMatrix = np.array([otherDir,np.cross(otherPlane,otherDir),otherPlane]).T
+      otherMatrix = np.array([otherDir,np.cross(otherPlane,otherDir),otherPlane])
 
-      r['rotations'].append(Rotation.fromMatrix(np.dot(otherMatrix,myMatrix.T)))
+      r['rotations'].append(Rotation.fromMatrix(np.dot(otherMatrix.T,myMatrix)))
 
     return r
 
@@ -1126,11 +1126,10 @@ class Orientation:
     return (Orientation(r,self.lattice), i,j, k == 1) if symmetries else r                          # disorientation ...
                                                                                                     # ... own sym, other sym,
                                                                                                     # self-->other: True, self<--other: False
-
-
   def inFZ(self):
     return self.lattice.symmetry.inFZ(self.rotation.asRodrigues(vector=True))
 
+  
   def equivalentOrientations(self,members=[]):
     """List of orientations which are symmetrically equivalent."""
     try:
@@ -1144,7 +1143,8 @@ class Orientation:
   def relatedOrientations(self,model):
     """List of orientations related by the given orientation relationship."""
     r = self.lattice.relationOperations(model)
-    return [self.__class__(self.rotation*o,r['lattice']) for o in r['rotations']]
+    return [self.__class__(o*self.rotation,r['lattice']) for o in r['rotations']]
+
     
   def reduced(self):
     """Transform orientation to fall into fundamental zone according to symmetry."""
@@ -1153,6 +1153,7 @@ class Orientation:
 
     return self.__class__(me.rotation,self.lattice)
    
+ 
   def inversePole(self,
                   axis,
                   proper = False,

python/damask/table.py

@@ -202,7 +202,7 @@ class Table():
                                                  '' if info is None else ': {}'.format(info),
                                                  ))
 
-        self.shapes[label_new] = self.shapes.pop(label_old)
+        self.shapes = {(label if label is not label_old else label_new):self.shapes[label] for label in self.shapes}
 
 
     def sort_by(self,labels,ascending=True):
@@ -233,8 +233,9 @@ class Table():
             Filename or file for reading.
 
         """
+        seen = set()
         labels = []
-        for l in self.shapes:
+        for l in [x for x in self.data.columns if not (x in seen or seen.add(x))]:
             if(self.shapes[l] == (1,)):
                 labels.append('{}'.format(l))
             elif(len(self.shapes[l]) == 1):

python/damask/util.py

@@ -221,263 +221,6 @@ class return_message():
     return srepr(self.message)
 
 
-def leastsqBound(func, x0, args=(), bounds=None, Dfun=None, full_output=0,
-                col_deriv=0, ftol=1.49012e-8, xtol=1.49012e-8,
-                gtol=0.0, maxfev=0, epsfcn=None, factor=100, diag=None):
-  from scipy.optimize import _minpack
-  """
-  Non-linear least square fitting (Levenberg-Marquardt method) with 
-  bounded parameters.
-  the codes of transformation between int <-> ext refers to the work of
-  Jonathan J. Helmus: https://github.com/jjhelmus/leastsqbound-scipy
-  other codes refer to the source code of minpack.py:
-
-  An internal parameter list is used to enforce contraints on the fitting 
-  parameters. The transfomation is based on that of MINUIT package. 
-  please see: F. James and M. Winkler. MINUIT User's Guide, 2004.
-
-  bounds : list
-    (min, max) pairs for each parameter, use None for 'min' or 'max' 
-    when there is no bound in that direction. 
-    For example: if there are two parameters needed to be fitting, then
-    bounds is [(min1,max1), (min2,max2)]
-    
-  This function is based on 'leastsq' of minpack.py, the annotation of
-  other parameters can be found in 'least_squares.py'.
-  """
- 
-  def _check_func(checker, argname, thefunc, x0, args, numinputs,
-                output_shape=None):
-    from numpy import shape
-    """The same as that of minpack.py"""
-    res = np.atleast_1d(thefunc(*((x0[:numinputs],) + args)))
-    if (output_shape is not None) and (shape(res) != output_shape):
-        if (output_shape[0] != 1):
-            if len(output_shape) > 1:
-                if output_shape[1] == 1:
-                    return shape(res)
-            msg = "%s: there is a mismatch between the input and output " \
-                  "shape of the '%s' argument" % (checker, argname)
-            func_name = getattr(thefunc, '__name__', None)
-            if func_name:
-                msg += " '%s'." % func_name
-            else:
-                msg += "."
-            raise TypeError(msg)
-    if np.issubdtype(res.dtype, np.inexact):
-        dt = res.dtype
-    else:
-        dt = dtype(float)
-    return shape(res), dt
-  
-  def _int2extGrad(p_int, bounds):
-    """Calculate the gradients of transforming the internal (unconstrained) to external (constrained) parameter."""
-    grad = np.empty_like(p_int)
-    for i, (x, bound) in enumerate(zip(p_int, bounds)):
-        lower, upper = bound
-        if lower is None and upper is None:  # No constraints
-            grad[i] = 1.0
-        elif upper is None:                  # only lower bound
-            grad[i] =  x/np.sqrt(x*x + 1.0)
-        elif lower is None:                  # only upper bound
-            grad[i] = -x/np.sqrt(x*x + 1.0)
-        else:                                # lower and upper bounds
-            grad[i] = (upper - lower)*np.cos(x)/2.0
-    return grad
-  
-  def _int2extFunc(bounds):
-    """Transform internal parameters into external parameters."""
-    local = [_int2extLocal(b) for b in bounds]
-
-    def _transform_i2e(p_int):
-        p_ext = np.empty_like(p_int)
-        p_ext[:] = [i(j) for i, j in zip(local, p_int)]
-        return p_ext
-    return _transform_i2e
-  
-  def _ext2intFunc(bounds):
-    """Transform external parameters into internal parameters."""
-    local = [_ext2intLocal(b) for b in bounds]
-
-    def _transform_e2i(p_ext):
-        p_int = np.empty_like(p_ext)
-        p_int[:] = [i(j) for i, j in zip(local, p_ext)]
-        return p_int
-    return _transform_e2i
-
-  def _int2extLocal(bound):
-    """Transform a single internal parameter to an external parameter."""
-    lower, upper = bound
-    if lower is None and upper is None:      # no constraints
-        return lambda x: x
-    elif upper is None:                      # only lower bound
-        return lambda x: lower - 1.0 + np.sqrt(x*x + 1.0)
-    elif lower is None:                      # only upper bound
-        return lambda x: upper + 1.0 - np.sqrt(x*x + 1.0)
-    else:
-        return lambda x: lower + ((upper - lower)/2.0)*(np.sin(x) + 1.0)
-  
-  def _ext2intLocal(bound):
-    """Transform a single external parameter to an internal parameter."""
-    lower, upper = bound
-    if lower is None and upper is None:  # no constraints
-        return lambda x: x
-    elif upper is None:                  # only lower bound
-        return lambda x: np.sqrt((x - lower + 1.0)**2 - 1.0)
-    elif lower is None:                  # only upper bound
-        return lambda x: np.sqrt((x - upper - 1.0)**2 - 1.0)
-    else:
-        return lambda x: np.arcsin((2.0*(x - lower)/(upper - lower)) - 1.0)
-  
-  i2e = _int2extFunc(bounds)
-  e2i = _ext2intFunc(bounds)
-  
-  x0 = np.asarray(x0).flatten()
-  n = len(x0)
-
-  if len(bounds) != n:
-      raise ValueError('the length of bounds is inconsistent with the number of parameters ')
-  
-  if not isinstance(args, tuple):
-      args = (args,)
-  
-  shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
-  m = shape[0]
-
-  if n > m:
-      raise TypeError('Improper input: N=%s must not exceed M=%s' % (n, m))
-  if epsfcn is None:
-      epsfcn = np.finfo(dtype).eps
-
-  def funcWarp(x, *args):
-      return func(i2e(x), *args)
-
-  xi0 = e2i(x0)
-  
-  if Dfun is None:
-      if maxfev == 0:
-          maxfev = 200*(n + 1)
-      retval = _minpack._lmdif(funcWarp, xi0, args, full_output, ftol, xtol,
-                               gtol, maxfev, epsfcn, factor, diag)
-  else:
-      if col_deriv:
-          _check_func('leastsq', 'Dfun', Dfun, x0, args, n, (n, m))
-      else:
-          _check_func('leastsq', 'Dfun', Dfun, x0, args, n, (m, n))
-      if maxfev == 0:
-          maxfev = 100*(n + 1)
-
-      def DfunWarp(x, *args):
-          return Dfun(i2e(x), *args)
-
-      retval = _minpack._lmder(funcWarp, DfunWarp, xi0, args, full_output, col_deriv,
-                               ftol, xtol, gtol, maxfev, factor, diag)
-
-  errors = {0: ["Improper input parameters.", TypeError],
-            1: ["Both actual and predicted relative reductions "
-                "in the sum of squares\n  are at most %f" % ftol, None],
-            2: ["The relative error between two consecutive "
-                "iterates is at most %f" % xtol, None],
-            3: ["Both actual and predicted relative reductions in "
-                "the sum of squares\n  are at most %f and the "
-                "relative error between two consecutive "
-                "iterates is at \n  most %f" % (ftol, xtol), None],
-            4: ["The cosine of the angle between func(x) and any "
-                "column of the\n  Jacobian is at most %f in "
-                "absolute value" % gtol, None],
-            5: ["Number of calls to function has reached "
-                "maxfev = %d." % maxfev, ValueError],
-            6: ["ftol=%f is too small, no further reduction "
-                "in the sum of squares\n  is possible.""" % ftol,
-                ValueError],
-            7: ["xtol=%f is too small, no further improvement in "
-                "the approximate\n  solution is possible." % xtol,
-                ValueError],
-            8: ["gtol=%f is too small, func(x) is orthogonal to the "
-                "columns of\n  the Jacobian to machine "
-                "precision." % gtol, ValueError],
-            'unknown': ["Unknown error.", TypeError]}
-
-  info = retval[-1]    # The FORTRAN return value
-  
-  if info not in [1, 2, 3, 4] and not full_output:
-      if info in [5, 6, 7, 8]:
-          np.warnings.warn(errors[info][0], RuntimeWarning)
-      else:
-          try:
-              raise errors[info][1](errors[info][0])
-          except KeyError:
-              raise errors['unknown'][1](errors['unknown'][0])
-
-  mesg = errors[info][0]
-  x = i2e(retval[0])
-
-  if full_output:
-      grad = _int2extGrad(retval[0], bounds)
-      retval[1]['fjac'] = (retval[1]['fjac'].T / np.take(grad,
-                           retval[1]['ipvt'] - 1)).T
-      cov_x = None
-      if info in [1, 2, 3, 4]:
-          from numpy.dual import inv
-          from numpy.linalg import LinAlgError
-          perm = np.take(np.eye(n), retval[1]['ipvt'] - 1, 0)
-          r = np.triu(np.transpose(retval[1]['fjac'])[:n, :])
-          R = np.dot(r, perm)
-          try:
-              cov_x = inv(np.dot(np.transpose(R), R))
-          except LinAlgError as inverror:
-              print(inverror)
-              pass
-      return (x, cov_x) + retval[1:-1] + (mesg, info)
-  else:
-      return (x, info)
-
-def _general_function(params, ydata, xdata, function):
-    return  function(xdata, *params) - ydata
-def _weighted_general_function(params, ydata, xdata, function, weights):
-    return (function(xdata, *params) - ydata)*weights
-
-def curve_fit_bound(f, xdata, ydata, p0=None, sigma=None, bounds=None, **kw):
-    """Similar as 'curve_fit' in minpack.py."""
-    if p0 is None:
-        # determine number of parameters by inspecting the function
-        import inspect
-        args, varargs, varkw, defaults = inspect.getargspec(f)
-        if len(args) < 2:
-            msg = "Unable to determine number of fit parameters."
-            raise ValueError(msg)
-        if 'self' in args:
-            p0 = [1.0] * (len(args)-2)
-        else:
-            p0 = [1.0] * (len(args)-1)
-
-    if np.isscalar(p0):
-        p0 = np.array([p0])
-
-    args = (ydata, xdata, f)
-    if sigma is None:
-        func = _general_function
-    else:
-        func = _weighted_general_function
-        args += (1.0/np.asarray(sigma),)
-
-    return_full = kw.pop('full_output', False)
-    res = leastsqBound(func, p0, args=args, bounds = bounds, full_output=True, **kw)
-    (popt, pcov, infodict, errmsg, ier) = res
-
-    if ier not in [1, 2, 3, 4]:
-        msg = "Optimal parameters not found: " + errmsg
-        raise RuntimeError(msg)
-
-    if (len(ydata) > len(p0)) and pcov is not None:
-        s_sq = (func(popt, *args)**2).sum()/(len(ydata)-len(p0))
-        pcov = pcov * s_sq
-    else:
-        pcov = np.inf
-
-    return (popt, pcov, infodict, errmsg, ier) if return_full else (popt, pcov)
-    
-
 class ThreadPool:
   """Pool of threads consuming tasks from a queue."""
 

python/setup.py

@@ -15,6 +15,7 @@ setuptools.setup(
     packages=setuptools.find_packages(),
     include_package_data=True,
     install_requires = [
+        "pandas",
         "scipy",
         "h5py",
         "vtk"

python/tests/reference/Rotation/PoleFigures_OR.m

@@ -0,0 +1,38 @@
+% Start MTEX first in Matlab
+
+tmp = matlab.desktop.editor.getActive;
+cd(fileparts(tmp.Filename));
+
+symmetry = {crystalSymmetry('m-3m', [1 1 1], 'mineral', 'Iron', 'color', 'light blue')}
+
+% plotting convention
+setMTEXpref('xAxisDirection','north');
+setMTEXpref('zAxisDirection','outOfPlane');
+
+
+lattice_types = {'BCC','FCC'};
+models        = {'Bain','GT','GT_prime','KS','NW','Pitsch'};
+
+rotation = containers.Map;
+rotation('BCC') = 'Passive Rotation';
+rotation('FCC') = 'Active Rotation';
+
+for lattice = lattice_types
+  for p = 0:length(models)/3-1
+    EBSD_data = {loadEBSD(strcat(lattice,'_',models{p*3+1},'.txt'),symmetry,'interface','generic',...
+                          'ColumnNames', { 'phi1' 'Phi' 'phi2' 'x' 'y'}, 'Bunge', rotation(char(lattice))),
+                 loadEBSD(strcat(lattice,'_',models{p*3+2},'.txt'),symmetry,'interface','generic',...
+                          'ColumnNames', { 'phi1' 'Phi' 'phi2' 'x' 'y'}, 'Bunge', rotation(char(lattice))),
+                 loadEBSD(strcat(lattice,'_',models{p*3+3},'.txt'),symmetry,'interface','generic',...
+                          'ColumnNames', { 'phi1' 'Phi' 'phi2' 'x' 'y'}, 'Bunge', rotation(char(lattice)))}
+    h = [Miller(1,0,0,symmetry{1}),Miller(1,1,0,symmetry{1}),Miller(1,1,1,symmetry{1})];            % 3 pole figures
+    plotPDF(EBSD_data{1}.orientations,h,'MarkerSize',5,'MarkerColor','r','DisplayName',models{p*3+1})
+    hold on
+    plotPDF(EBSD_data{2}.orientations,h,'MarkerSize',4,'MarkerColor','b','DisplayName',models{p*3+2})
+    plotPDF(EBSD_data{3}.orientations,h,'MarkerSize',3,'MarkerColor','g','DisplayName',models{p*3+3})
+    legend('show','location','southoutside','Interpreter', 'none') 
+    orient('landscape')
+    print('-bestfit',strcat(int2str(p+1),'_',char(lattice),'.pdf'),'-dpdf')
+    close
+  end
+end

python/tests/reference/Rotation/bcc_Bain.txt

@@ -0,0 +1,5 @@
+1 header
+1_Eulers 2_Eulers 3_Eulers 1_pos 2_pos
+0.0 45.00000000000001 0.0 1 1
+90.0 45.00000000000001 270.0 1 2
+45.00000000000001 0.0 0.0 1 3

python/tests/reference/Rotation/bcc_GT.txt

@@ -0,0 +1,26 @@
+1 header
+1_Eulers 2_Eulers 3_Eulers 1_pos 2_pos
+283.60440567265294 9.976439066337804 33.24637065555936 1 1
+167.8261034151001 43.397849654402556 183.40022280897963 1 2
+262.1156357053931 43.82007387041961 104.07478363123654 1 3
+103.604405672653 9.976439066337804 213.24637065555936 1 4
+347.8261034151001 43.39784965440255 3.400222808979685 1 5
+82.11563570539313 43.82007387041961 284.0747836312365 1 6
+76.39559432734703 9.976439066337806 326.75362934444064 1 7
+192.17389658489986 43.397849654402556 176.59977719102034 1 8
+97.88436429460687 43.82007387041961 255.92521636876344 1 9
+256.395594327347 9.976439066337804 146.75362934444064 1 10
+12.173896584899929 43.39784965440254 356.59977719102034 1 11
+277.8843642946069 43.82007387041961 75.92521636876346 1 12
+102.17389658489992 43.39784965440254 266.59977719102034 1 13
+346.395594327347 9.976439066337804 56.75362934444064 1 14
+7.884364294606862 43.82007387041961 345.9252163687635 1 15
+282.17389658489986 43.39784965440254 86.59977719102032 1 16
+166.39559432734703 9.976439066337804 236.75362934444058 1 17
+187.88436429460683 43.82007387041961 165.92521636876344 1 18
+257.8261034151001 43.39784965440255 93.40022280897969 1 19
+13.604405672652977 9.976439066337804 303.24637065555936 1 20
+352.1156357053931 43.82007387041961 14.074783631236542 1 21
+77.82610341510008 43.397849654402556 273.4002228089796 1 22
+193.60440567265297 9.976439066337806 123.24637065555939 1 23
+172.11563570539317 43.82007387041961 194.07478363123653 1 24

python/tests/reference/Rotation/bcc_GT_prime.txt

@@ -0,0 +1,26 @@
+1 header
+1_Eulers 2_Eulers 3_Eulers 1_pos 2_pos
+303.24637065555936 9.976439066337804 13.604405672652977 1 1
+165.92521636876344 43.82007387041961 187.88436429460683 1 2
+266.59977719102034 43.39784965440254 102.17389658489992 1 3
+123.24637065555939 9.976439066337804 193.604405672653 1 4
+345.9252163687635 43.82007387041961 7.884364294606862 1 5
+86.59977719102032 43.39784965440254 282.17389658489986 1 6
+56.75362934444064 9.976439066337804 346.395594327347 1 7
+194.07478363123653 43.82007387041961 172.11563570539317 1 8
+93.40022280897969 43.39784965440255 257.8261034151001 1 9
+236.75362934444058 9.976439066337804 166.39559432734697 1 10
+14.074783631236542 43.82007387041961 352.1156357053931 1 11
+273.4002228089796 43.397849654402556 77.82610341510008 1 12
+104.07478363123654 43.82007387041961 262.1156357053931 1 13
+326.75362934444064 9.976439066337806 76.39559432734703 1 14
+3.400222808979685 43.39784965440255 347.8261034151001 1 15
+284.0747836312365 43.82007387041961 82.11563570539313 1 16
+146.75362934444064 9.976439066337804 256.395594327347 1 17
+183.40022280897963 43.397849654402556 167.8261034151001 1 18
+255.92521636876344 43.82007387041961 97.88436429460687 1 19
+33.24637065555936 9.976439066337804 283.60440567265294 1 20
+356.59977719102034 43.39784965440254 12.173896584899929 1 21
+75.92521636876346 43.82007387041961 277.8843642946069 1 22
+213.24637065555936 9.976439066337804 103.604405672653 1 23
+176.59977719102034 43.397849654402556 192.17389658489986 1 24

python/tests/reference/Rotation/bcc_KS.txt

@@ -0,0 +1,26 @@
+1 header
+1_Eulers 2_Eulers 3_Eulers 1_pos 2_pos
+335.7965716606702 10.528779365509317 65.79657166067024 1 1
+228.77270547567446 80.40593177313953 85.64260312151849 1 2
+131.22729452432552 80.40593177313954 4.357396878481506 1 3
+24.20342833932977 10.52877936550932 24.20342833932976 1 4
+221.95489158457983 85.70366403943002 80.37863910890589 1 5
+138.04510841542015 85.70366403943004 9.621360891094124 1 6
+131.22729452432552 80.40593177313953 94.35739687848151 1 7
+24.203428339329765 10.52877936550932 114.20342833932976 1 8
+221.95489158457983 85.70366403943004 170.37863910890587 1 9
+138.04510841542015 85.70366403943004 99.62136089109411 1 10
+335.7965716606702 10.52877936550932 155.79657166067025 1 11
+228.77270547567448 80.40593177313954 175.6426031215185 1 12
+335.7965716606702 10.52877936550932 335.7965716606702 1 13
+228.77270547567448 80.40593177313954 355.6426031215185 1 14
+131.2272945243255 80.40593177313954 274.35739687848144 1 15
+24.203428339329747 10.52877936550932 294.2034283393298 1 16
+221.95489158457985 85.70366403943004 350.3786391089059 1 17
+138.04510841542015 85.70366403943004 279.6213608910941 1 18
+41.95489158457986 94.29633596056998 9.621360891094133 1 19
+318.04510841542015 94.29633596056996 80.37863910890589 1 20
+155.79657166067025 169.4712206344907 24.203428339329754 1 21
+48.77270547567448 99.59406822686046 4.357396878481504 1 22
+311.2272945243255 99.59406822686046 85.64260312151852 1 23
+204.20342833932975 169.4712206344907 65.79657166067024 1 24

python/tests/reference/Rotation/bcc_NW.txt

@@ -0,0 +1,14 @@
+1 header
+1_Eulers 2_Eulers 3_Eulers 1_pos 2_pos
+225.41555594321144 83.13253115922213 83.08266205989301 1 1
+134.58444405678856 83.13253115922211 6.917337940107012 1 2
+4.702125169424418e-15 9.735610317245317 45.0 1 3
+134.58444405678856 83.13253115922213 276.91733794010696 1 4
+225.4155559432114 83.13253115922213 353.082662059893 1 5
+0.0 9.735610317245317 315.0 1 6
+134.58444405678858 83.13253115922213 96.91733794010702 1 7
+225.41555594321142 83.13253115922213 173.082662059893 1 8
+0.0 9.735610317245317 135.0 1 9
+99.59803029876785 45.81931182053557 166.36129272052355 1 10
+260.40196970123213 45.81931182053556 283.6387072794765 1 11
+180.0 99.73561031724535 225.0 1 12

python/tests/reference/Rotation/bcc_Pitsch.txt

@@ -0,0 +1,14 @@
+1 header
+1_Eulers 2_Eulers 3_Eulers 1_pos 2_pos
+6.9173379401070045 83.13253115922213 44.58444405678856 1 1
+45.0 89.99999999999999 279.7356103172453 1 2
+166.36129272052352 45.819311820535574 279.59803029876787 1 3
+83.08266205989301 83.13253115922213 225.41555594321144 1 4
+256.3612927205235 45.819311820535574 189.59803029876787 1 5
+315.0 90.0 9.735610317245369 1 6
+186.917337940107 83.13253115922213 224.58444405678856 1 7
+315.0 90.0 80.26438968275463 1 8
+13.638707279476478 45.81931182053557 260.40196970123213 1 9
+263.082662059893 83.13253115922213 45.415555943211444 1 10
+103.63870727947646 45.819311820535574 170.40196970123213 1 11
+224.99999999999997 90.0 170.26438968275465 1 12

python/tests/reference/Rotation/fcc_Bain.txt

@@ -0,0 +1,5 @@
+1 header
+1_Eulers 2_Eulers 3_Eulers 1_pos 2_pos
+180.0 45.00000000000001 180.0 1 1
+270.0 45.00000000000001 90.0 1 2
+315.0 0.0 0.0 1 3

python/tests/reference/Rotation/fcc_GT.txt

@@ -0,0 +1,26 @@
+1 header
+1_Eulers 2_Eulers 3_Eulers 1_pos 2_pos
+146.75362934444064 9.976439066337804 256.395594327347 1 1
+356.59977719102034 43.39784965440254 12.173896584899929 1 2
+75.92521636876346 43.82007387041961 277.8843642946069 1 3
+326.75362934444064 9.976439066337806 76.39559432734703 1 4
+176.59977719102034 43.397849654402556 192.17389658489986 1 5
+255.92521636876344 43.82007387041961 97.88436429460687 1 6
+213.24637065555936 9.976439066337804 103.604405672653 1 7
+3.400222808979685 43.39784965440255 347.8261034151001 1 8
+284.0747836312365 43.82007387041961 82.11563570539313 1 9
+33.24637065555936 9.976439066337804 283.60440567265294 1 10
+183.40022280897963 43.397849654402556 167.8261034151001 1 11
+104.07478363123654 43.82007387041961 262.1156357053931 1 12
+273.4002228089796 43.397849654402556 77.82610341510008 1 13
+123.24637065555939 9.976439066337806 193.60440567265297 1 14
+194.07478363123653 43.82007387041961 172.11563570539317 1 15
+93.40022280897969 43.39784965440255 257.8261034151001 1 16
+303.24637065555936 9.976439066337804 13.604405672652977 1 17
+14.074783631236542 43.82007387041961 352.1156357053931 1 18
+86.59977719102032 43.39784965440254 282.17389658489986 1 19
+236.75362934444058 9.976439066337804 166.39559432734703 1 20
+165.92521636876344 43.82007387041961 187.88436429460683 1 21
+266.59977719102034 43.39784965440254 102.17389658489992 1 22
+56.75362934444064 9.976439066337804 346.395594327347 1 23
+345.9252163687635 43.82007387041961 7.884364294606862 1 24

python/tests/reference/Rotation/fcc_GT_prime.txt

@@ -0,0 +1,26 @@
+1 header
+1_Eulers 2_Eulers 3_Eulers 1_pos 2_pos
+166.39559432734697 9.976439066337804 236.75362934444058 1 1
+352.1156357053931 43.82007387041961 14.074783631236542 1 2
+77.82610341510008 43.397849654402556 273.4002228089796 1 3
+346.395594327347 9.976439066337804 56.75362934444064 1 4
+172.11563570539317 43.82007387041961 194.07478363123653 1 5
+257.8261034151001 43.39784965440255 93.40022280897969 1 6
+193.604405672653 9.976439066337804 123.24637065555939 1 7
+7.884364294606862 43.82007387041961 345.9252163687635 1 8
+282.17389658489986 43.39784965440254 86.59977719102032 1 9
+13.604405672652977 9.976439066337804 303.24637065555936 1 10
+187.88436429460683 43.82007387041961 165.92521636876344 1 11
+102.17389658489992 43.39784965440254 266.59977719102034 1 12
+277.8843642946069 43.82007387041961 75.92521636876346 1 13
+103.604405672653 9.976439066337804 213.24637065555936 1 14
+192.17389658489986 43.397849654402556 176.59977719102034 1 15
+97.88436429460687 43.82007387041961 255.92521636876344 1 16
+283.60440567265294 9.976439066337804 33.24637065555936 1 17
+12.173896584899929 43.39784965440254 356.59977719102034 1 18
+82.11563570539313 43.82007387041961 284.0747836312365 1 19
+256.395594327347 9.976439066337804 146.75362934444064 1 20
+167.8261034151001 43.397849654402556 183.40022280897963 1 21
+262.1156357053931 43.82007387041961 104.07478363123654 1 22
+76.39559432734703 9.976439066337806 326.75362934444064 1 23
+347.8261034151001 43.39784965440255 3.400222808979685 1 24

python/tests/reference/Rotation/fcc_KS.txt

@@ -0,0 +1,26 @@
+1 header
+1_Eulers 2_Eulers 3_Eulers 1_pos 2_pos
+114.20342833932975 10.52877936550932 204.20342833932972 1 1
+94.3573968784815 80.40593177313954 311.22729452432543 1 2
+175.6426031215185 80.40593177313954 48.77270547567447 1 3
+155.79657166067025 10.52877936550932 155.79657166067025 1 4
+99.62136089109411 85.70366403943004 318.04510841542015 1 5
+170.37863910890587 85.70366403943002 41.954891584579855 1 6
+85.64260312151852 80.40593177313954 48.77270547567448 1 7
+65.79657166067024 10.52877936550932 155.79657166067025 1 8
+9.621360891094124 85.70366403943004 318.04510841542015 1 9
+80.37863910890587 85.70366403943004 41.95489158457987 1 10
+24.203428339329758 10.52877936550932 204.20342833932975 1 11
+4.357396878481486 80.40593177313954 311.2272945243255 1 12
+204.20342833932972 10.52877936550932 204.20342833932972 1 13
+184.35739687848147 80.40593177313954 311.2272945243255 1 14
+265.64260312151845 80.40593177313953 48.77270547567449 1 15
+245.79657166067025 10.528779365509317 155.79657166067025 1 16
+189.62136089109413 85.70366403943004 318.04510841542015 1 17
+260.3786391089059 85.70366403943002 41.954891584579855 1 18
+170.37863910890587 94.29633596056996 138.04510841542015 1 19
+99.62136089109411 94.29633596056998 221.95489158457983 1 20
+155.79657166067025 169.4712206344907 24.203428339329754 1 21
+175.64260312151848 99.59406822686046 131.22729452432552 1 22
+94.35739687848151 99.59406822686046 228.77270547567446 1 23
+114.20342833932975 169.4712206344907 335.7965716606702 1 24

python/tests/reference/Rotation/fcc_NW.txt

@@ -0,0 +1,14 @@
+1 header
+1_Eulers 2_Eulers 3_Eulers 1_pos 2_pos
+96.91733794010702 83.13253115922213 314.5844440567886 1 1
+173.082662059893 83.13253115922211 45.41555594321143 1 2
+135.0 9.735610317245317 180.0 1 3
+263.082662059893 83.13253115922213 45.415555943211444 1 4
+186.91733794010702 83.13253115922211 314.5844440567886 1 5
+224.99999999999997 9.735610317245317 180.0 1 6
+83.082662059893 83.13253115922213 45.415555943211444 1 7
+6.917337940106983 83.13253115922211 314.5844440567886 1 8
+45.0 9.73561031724532 180.0 1 9
+13.638707279476469  45.81931182053557 80.40196970123216 1 10
+256.36129272052347 45.81931182053556 279.59803029876775 1 11
+315.0 99.73561031724536 0.0 1 12

python/tests/reference/Rotation/fcc_Pitsch.txt

@@ -0,0 +1,14 @@
+1 header
+1_Eulers 2_Eulers 3_Eulers 1_pos 2_pos
+135.41555594321144 83.13253115922213 173.082662059893 1 1
+260.26438968275465 90.0 135.0 1 2
+260.40196970123213 45.81931182053557 13.638707279476478 1 3
+314.5844440567886 83.13253115922213 96.91733794010702 1 4
+350.40196970123213 45.81931182053557 283.6387072794765 1 5
+170.26438968275465 90.0 224.99999999999997 1 6
+315.4155559432114 83.13253115922213 353.08266205989304 1 7
+99.73561031724536 90.0 225.0 1 8
+279.59803029876787 45.819311820535574 166.36129272052352 1 9
+134.58444405678856 83.13253115922213 276.91733794010696 1 10
+9.598030298767851 45.819311820535574 76.36129272052355 1 11
+9.735610317245369 90.0 315.0 1 12

python/tests/test_DADF5.py

@@ -24,12 +24,19 @@ def reference_dir(reference_dir_base):
 
 class TestDADF5:
 
-    def test_add_deviator(self,default):
+    def test_time_increments(self,default):
-        default.add_deviator('P')
+        shape = default.read_dataset(default.get_dataset_location('F'),0).shape
-        loc = {'P'  :default.get_dataset_location('P'),
+        default.set_by_time(0.0,20.0)
-               's_P':default.get_dataset_location('s_P')}
+        for i in default.iter_visible('increments'):
-        in_memory = mechanics.deviatoric_part(default.read_dataset(loc['P'],0))
+           assert shape == default.read_dataset(default.get_dataset_location('F'),0).shape
-        in_file   = default.read_dataset(loc['s_P'],0)
+
+
+    def test_add_absolute(self,default):
+        default.add_absolute('Fe')
+        loc = {'Fe':   default.get_dataset_location('Fe'),
+               '|Fe|': default.get_dataset_location('|Fe|')}
+        in_memory = np.abs(default.read_dataset(loc['Fe'],0))
+        in_file   = default.read_dataset(loc['|Fe|'],0)
         assert np.allclose(in_memory,in_file)
 
     def test_add_Cauchy(self,default):
@@ -42,22 +49,30 @@ class TestDADF5:
         in_file   = default.read_dataset(loc['sigma'],0)
         assert np.allclose(in_memory,in_file)
 
-    def test_add_absolute(self,default):
-        default.add_absolute('Fe')
-        loc = {'Fe':   default.get_dataset_location('Fe'),
-               '|Fe|': default.get_dataset_location('|Fe|')}
-        in_memory = np.abs(default.read_dataset(loc['Fe'],0))
-        in_file   = default.read_dataset(loc['|Fe|'],0)
-        assert np.allclose(in_memory,in_file)
-
     def test_add_determinant(self,default):
         default.add_determinant('P')
         loc = {'P':     default.get_dataset_location('P'),
-               'det(P)': default.get_dataset_location('det(P)')}
+               'det(P)':default.get_dataset_location('det(P)')}
-        in_memory = np.linalg.det(default.read_dataset(loc['P'],0)).reshape(-1,1)
+        in_memory = np.linalg.det(default.read_dataset(loc['P'],0)).reshape((-1,1))
         in_file   = default.read_dataset(loc['det(P)'],0)
         assert np.allclose(in_memory,in_file)
 
+    def test_add_deviator(self,default):
+        default.add_deviator('P')
+        loc = {'P'  :default.get_dataset_location('P'),
+               's_P':default.get_dataset_location('s_P')}
+        in_memory = mechanics.deviatoric_part(default.read_dataset(loc['P'],0))
+        in_file   = default.read_dataset(loc['s_P'],0)
+        assert np.allclose(in_memory,in_file)
+
+    def test_add_norm(self,default):
+        default.add_norm('F',1)
+        loc = {'F':    default.get_dataset_location('F'),
+               '|F|_1':default.get_dataset_location('|F|_1')}
+        in_memory = np.linalg.norm(default.read_dataset(loc['F'],0),ord=1,axis=(1,2),keepdims=True)
+        in_file   = default.read_dataset(loc['|F|_1'],0)
+        assert np.allclose(in_memory,in_file)
+
     def test_add_spherical(self,default):
         default.add_spherical('P')
         loc = {'P':   default.get_dataset_location('P'),

python/tests/test_Rotation.py

@@ -1,7 +1,11 @@
+import os
+
 import pytest
 import numpy as np
 
+import damask
 from damask import Rotation
+from damask import Orientation
    
 n = 1000
 
@@ -10,6 +14,11 @@ def default():
     """A set of n random rotations."""
     return [Rotation.fromRandom() for r in range(n)]
 
+@pytest.fixture
+def reference_dir(reference_dir_base):
+    """Directory containing reference results."""
+    return os.path.join(reference_dir_base,'Rotation')
+
 
 class TestRotation:
 
@@ -18,38 +27,55 @@ class TestRotation:
             assert np.allclose(rot.asQuaternion(), 
                                Rotation.fromEulers(rot.asEulers()).asQuaternion())
 
-
     def test_AxisAngle(self,default):
         for rot in default:
             assert np.allclose(rot.asEulers(),
                                Rotation.fromAxisAngle(rot.asAxisAngle()).asEulers())
 
-
     def test_Matrix(self,default):
         for rot in default:
             assert np.allclose(rot.asAxisAngle(),
                                Rotation.fromMatrix(rot.asMatrix()).asAxisAngle())
 
-
     def test_Rodriques(self,default):
         for rot in default:
             assert np.allclose(rot.asMatrix(),
                                Rotation.fromRodrigues(rot.asRodrigues()).asMatrix())
 
-
     def test_Homochoric(self,default):
         for rot in default:
             assert np.allclose(rot.asRodrigues(),
                                Rotation.fromHomochoric(rot.asHomochoric()).asRodrigues())
 
-
     def test_Cubochoric(self,default):
         for rot in default:
             assert np.allclose(rot.asHomochoric(),
                                Rotation.fromCubochoric(rot.asCubochoric()).asHomochoric())
 
-
     def test_Quaternion(self,default):
         for rot in default:
             assert np.allclose(rot.asCubochoric(),
                                Rotation.fromQuaternion(rot.asQuaternion()).asCubochoric())
+
+
+    @pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])
+    @pytest.mark.parametrize('lattice',['fcc','bcc'])
+    def test_relationship_forward_backward(self,model,lattice):
+        ori = Orientation(Rotation.fromRandom(),lattice)
+        for i,r in enumerate(ori.relatedOrientations(model)):
+            ori2 = r.relatedOrientations(model)[i]
+            misorientation = ori.rotation.misorientation(ori2.rotation)
+            assert misorientation.asAxisAngle(degrees=True)[3]<1.0e-5
+
+    @pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])
+    @pytest.mark.parametrize('lattice',['fcc','bcc'])
+    def test_relationship_reference(self,update,reference_dir,model,lattice):
+        reference = os.path.join(reference_dir,'{}_{}.txt'.format(lattice,model))
+        ori = Orientation(Rotation(),lattice)
+        eu = np.array([o.rotation.asEulers(degrees=True) for o in ori.relatedOrientations(model)])
+        if update: 
+            coords = np.array([(1,i+1) for i,x in enumerate(eu)])
+            table = damask.Table(eu,{'Eulers':(3,)})
+            table.add('pos',coords)
+            table.to_ASCII(reference)
+        assert np.allclose(eu,damask.Table.from_ASCII(reference).get('Eulers'))

python/tests/test_Table.py

@@ -65,11 +65,13 @@ class TestTable:
         default.add('nine',d,'random data')
         assert np.allclose(d,default.get('nine'))
 
-    def test_rename_equivalent(self,default):
+    def test_rename_equivalent(self):
-        v = default.get('v')
+        x = np.random.random((5,13))
-        default.rename('v','u')
+        t = Table(x,{'F':(3,3),'v':(3,),'s':(1,)},['random test data'])
-        u = default.get('u')
+        s = t.get('s')
-        assert np.all(v == u)
+        t.rename('s','u')
+        u = t.get('u')
+        assert np.all(s == u)
 
     def test_rename_gone(self,default):
         default.rename('v','V')

src/CPFEM2.f90

@@ -42,7 +42,7 @@ contains
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief call (thread safe) all module initializations
+!> @brief call all module initializations
 !--------------------------------------------------------------------------------------------------
 subroutine CPFEM_initAll
 
@@ -65,33 +65,27 @@ subroutine CPFEM_initAll
   call constitutive_init
   call crystallite_init
   call homogenization_init
-  call materialpoint_postResults
   call CPFEM_init
 
 end subroutine CPFEM_initAll
 
+
 !--------------------------------------------------------------------------------------------------
 !> @brief allocate the arrays defined in module CPFEM and initialize them
 !--------------------------------------------------------------------------------------------------
 subroutine CPFEM_init
  
-  integer :: ph,homog
+  integer :: i
-  character(len=1024) :: rankStr, PlasticItem, HomogItem
+  integer(HID_T) :: fileHandle, groupHandle
-  integer(HID_T) :: fileHandle, groupPlasticID, groupHomogID
+  character(len=pStringLen) :: fileName, datasetName
 
-  write(6,'(/,a)')   ' <<<+-  CPFEM init  -+>>>'
+  write(6,'(/,a)')   ' <<<+-  CPFEM init  -+>>>'; flush(6)
-  flush(6)
 
-  ! *** restore the last converged values of each essential variable
   if (interface_restartInc > 0) then
-    if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0) then
+    write(6,'(/,a,i0,a)') ' reading restart information of increment ', interface_restartInc, ' from file'
-      write(6,'(a)') '<< CPFEM >> restored state variables of last converged step from hdf5 file'
-      flush(6)
-    endif
 
-    write(rankStr,'(a1,i0)')'_',worldrank
+    write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5'
-
+    fileHandle = HDF5_openFile(fileName)
-    fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5')
     
     call HDF5_read(fileHandle,crystallite_F0, 'F')
     call HDF5_read(fileHandle,crystallite_Fp0,'Fp')
@@ -100,19 +94,19 @@ subroutine CPFEM_init
     call HDF5_read(fileHandle,crystallite_Li0,'Li')
     call HDF5_read(fileHandle,crystallite_S0, 'S')
     
-    groupPlasticID = HDF5_openGroup(fileHandle,'constituent')
+    groupHandle = HDF5_openGroup(fileHandle,'constituent')
-    do ph = 1,size(phase_plasticity)
+    do i = 1,size(phase_plasticity)
-      write(PlasticItem,*) ph,'_'
+      write(datasetName,'(i0,a)') i,'_omega_plastic'
-      call HDF5_read(groupPlasticID,plasticState(ph)%state0,trim(PlasticItem)//'omega_plastic')
+      call HDF5_read(groupHandle,plasticState(i)%state0,datasetName)
     enddo
-    call HDF5_closeGroup(groupPlasticID)
+    call HDF5_closeGroup(groupHandle)
 
-    groupHomogID = HDF5_openGroup(fileHandle,'materialpoint')
+    groupHandle = HDF5_openGroup(fileHandle,'materialpoint')
-    do homog = 1, material_Nhomogenization
+    do i = 1, material_Nhomogenization
-      write(HomogItem,*) homog,'_'
+      write(datasetName,'(i0,a)') i,'_omega_homogenization'
-      call HDF5_read(groupHomogID,homogState(homog)%state0, trim(HomogItem)//'omega_homogenization')
+      call HDF5_read(groupHandle,homogState(i)%state0,datasetName)
     enddo
-    call HDF5_closeGroup(groupHomogID)
+    call HDF5_closeGroup(groupHandle)
 
     call HDF5_closeFile(fileHandle)
   endif
@@ -126,7 +120,7 @@ end subroutine CPFEM_init
 !--------------------------------------------------------------------------------------------------
 subroutine CPFEM_forward
 
-  integer :: i, homog, mySource
+  integer :: i, j
 
   if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0) &
     write(6,'(a)') '<< CPFEM >> aging states'
@@ -142,32 +136,31 @@ subroutine CPFEM_forward
     plasticState(i)%state0 = plasticState(i)%state
   enddo 
   do i = 1, size(sourceState)
-    do mySource = 1,phase_Nsources(i)
+    do j = 1,phase_Nsources(i)
-      sourceState(i)%p(mySource)%state0 = sourceState(i)%p(mySource)%state
+      sourceState(i)%p(j)%state0 = sourceState(i)%p(j)%state
   enddo; enddo
-  do homog = 1, material_Nhomogenization
+  do i = 1, material_Nhomogenization
-    homogState       (homog)%state0 =  homogState       (homog)%state
+    homogState  (i)%state0 = homogState  (i)%state
-    thermalState     (homog)%state0 =  thermalState     (homog)%state
+    thermalState(i)%state0 = thermalState(i)%state
-    damageState      (homog)%state0 =  damageState      (homog)%state
+    damageState (i)%state0 = damageState (i)%state
   enddo
 
 end subroutine CPFEM_forward
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Write current constitutive variables for restart to file.
+!> @brief Write current  restart information (Field and constitutive data) to file.
 !--------------------------------------------------------------------------------------------------
 subroutine CPFEM_restartWrite
 
-  integer           :: ph, homog
+  integer :: i
-  character(len=32) :: rankStr, PlasticItem, HomogItem
+  integer(HID_T) :: fileHandle, groupHandle
-  integer(HID_T)    :: fileHandle, groupPlastic, groupHomog
+  character(len=pStringLen) :: fileName, datasetName
 
+  write(6,'(a)') ' writing field and constitutive data required for restart to file';flush(6)
     
-  write(6,'(a)') ' writing constitutive data required for restart to file';flush(6)
+  write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5'
-    
+  fileHandle = HDF5_openFile(fileName,'a')
-  write(rankStr,'(a1,i0)')'_',worldrank
-  fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','a')
     
   call HDF5_write(fileHandle,crystallite_partionedF,'F')
   call HDF5_write(fileHandle,crystallite_Fp,        'Fp')
@@ -176,19 +169,19 @@ subroutine CPFEM_restartWrite
   call HDF5_write(fileHandle,crystallite_Li,        'Li')
   call HDF5_write(fileHandle,crystallite_S,         'S')
     
-  groupPlastic = HDF5_addGroup(fileHandle,'constituent')
+  groupHandle = HDF5_addGroup(fileHandle,'constituent')
-  do ph = 1,size(phase_plasticity)
+  do i = 1,size(phase_plasticity)
-    write(PlasticItem,*) ph,'_'
+    write(datasetName,'(i0,a)') i,'_omega_plastic'
-    call HDF5_write(groupPlastic,plasticState(ph)%state,trim(PlasticItem)//'omega_plastic')
+    call HDF5_write(groupHandle,plasticState(i)%state,datasetName)
   enddo
-  call HDF5_closeGroup(groupPlastic)
+  call HDF5_closeGroup(groupHandle)
 
-  groupHomog = HDF5_addGroup(fileHandle,'materialpoint')
+  groupHandle = HDF5_addGroup(fileHandle,'materialpoint')
-  do homog = 1, material_Nhomogenization
+  do i = 1, material_Nhomogenization
-    write(HomogItem,*) homog,'_'
+    write(datasetName,'(i0,a)') i,'_omega_homogenization'
-    call HDF5_write(groupHomog,homogState(homog)%state,trim(HomogItem)//'omega_homogenization')
+    call HDF5_write(groupHandle,homogState(i)%state,datasetName)
   enddo
-  call HDF5_closeGroup(groupHomog)
+  call HDF5_closeGroup(groupHandle)
     
   call HDF5_closeFile(fileHandle)
 
@@ -209,7 +202,7 @@ subroutine CPFEM_results(inc,time)
   call crystallite_results
   call homogenization_results
   call discretization_results
-  call results_removeLink('current') ! ToDo: put this into closeJobFile
+  call results_removeLink('current') ! ToDo: put this into closeJobFile?
   call results_closeJobFile
 
 end subroutine CPFEM_results

src/DAMASK_abaqus.f

@@ -143,9 +143,6 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
    outdatedByNewInc, &
    outdatedFFN1, &
    lastStep
- use homogenization, only: &
-   materialpoint_sizeResults, &
-   materialpoint_results
 
  implicit none
  integer(pInt),                       intent(in) :: &
@@ -332,7 +329,7 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
    ddsdde(6,:) = ddsdde_h(5,:)
  end if
 
- statev = materialpoint_results(1:min(nstatv,materialpoint_sizeResults),npt,mesh_FEasCP('elem', noel))
+ statev = 0
 
  if (terminallyIll) pnewdt = 0.5_pReal                                                              ! force cutback directly ?
 !$ call omp_set_num_threads(defaultNumThreadsInt)                                                   ! reset number of threads to stored default value

src/IO.f90

@@ -32,8 +32,7 @@ module IO
     IO_intValue, &
     IO_lc, &
     IO_error, &
-    IO_warning, &
+    IO_warning
-    IO_intOut
 #if defined(Marc4DAMASK) || defined(Abaqus)
   public :: &
     IO_open_inputFile, &
@@ -542,25 +541,6 @@ pure function IO_lc(string)
 end function IO_lc
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief returns format string for integer values without leading zeros
-!--------------------------------------------------------------------------------------------------
-pure function IO_intOut(intToPrint)
-
-  integer, intent(in) :: intToPrint
-  character(len=41)   :: IO_intOut
-  integer             :: N_digits
-  character(len=19)   :: width                                                                      ! maximum digits for 64 bit integer
-  character(len=20)   :: min_width                                                                  ! longer for negative values
-
-  N_digits =  1 + int(log10(real(max(abs(intToPrint),1))))
-  write(width, '(I19.19)') N_digits
-  write(min_width, '(I20.20)') N_digits + merge(1,0,intToPrint < 0)
-  IO_intOut = 'I'//trim(min_width)//'.'//trim(width)
-
-end function IO_intOut
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief write error statements to standard out and terminate the Marc/spectral run with exit #9xxx
 !> in ABAQUS either time step is reduced or execution terminated

src/config.f90

@@ -27,7 +27,7 @@ module config
     config_numerics, &
     config_debug
  
-  character(len=64), dimension(:), allocatable, public, protected :: &
+  character(len=pStringLen), dimension(:), allocatable, public, protected :: &
     config_name_phase, &                                                                            !< name of each phase
     config_name_homogenization, &                                                                   !< name of each homogenization
     config_name_crystallite, &                                                                      !< name of each crystallite setting
@@ -54,7 +54,7 @@ subroutine config_init
   character(len=pStringLen), dimension(:), allocatable :: fileContent
   logical :: fileExists
 
-  write(6,'(/,a)') ' <<<+-  config init  -+>>>'
+  write(6,'(/,a)') ' <<<+-  config init  -+>>>'; flush(6)
 
   verbose = iand(debug_level(debug_material),debug_levelBasic) /= 0
 
@@ -214,7 +214,7 @@ end function read_materialConfig
 subroutine parse_materialConfig(sectionNames,part,line, &
                                 fileContent)
 
-  character(len=64),            allocatable, dimension(:), intent(out)   :: sectionNames
+  character(len=pStringLen),    allocatable, dimension(:), intent(out)   :: sectionNames
   type(tPartitionedStringList), allocatable, dimension(:), intent(inout) :: part
   character(len=pStringLen),                               intent(inout) :: line
   character(len=pStringLen),                 dimension(:), intent(in)    :: fileContent
@@ -222,7 +222,7 @@ subroutine parse_materialConfig(sectionNames,part,line, &
   integer, allocatable, dimension(:) :: partPosition                                                !< position of [] tags + last line in section
   integer :: i, j
   logical :: echo
-  character(len=pStringLen) :: section_ID
+  character(len=pStringLen) :: sectionName
 
   echo = .false. 
 
@@ -246,8 +246,8 @@ subroutine parse_materialConfig(sectionNames,part,line, &
   partPosition = [partPosition, i]                                                                  ! needed when actually storing content
 
   do i = 1, size(partPosition) -1
-    write(section_ID,'('//IO_intOut(size(partPosition))//')') i
+    write(sectionName,'(i0,a,a)') i,'_',trim(IO_getTag(fileContent(partPosition(i)),'[',']'))
-    sectionNames(i) = trim(section_ID)//'_'//trim(adjustl(IO_getTag(fileContent(partPosition(i)),'[',']')))
+    sectionNames(i) = sectionName
     do j = partPosition(i) + 1,  partPosition(i+1) -1
       call part(i)%add(trim(adjustl(fileContent(j))))
     enddo

src/constitutive.f90

@@ -37,7 +37,6 @@ module constitutive
   
   integer, public, protected :: &
     constitutive_plasticity_maxSizeDotState, &
-    constitutive_source_maxSizePostResults, &
     constitutive_source_maxSizeDotState
  
   public :: &
@@ -50,7 +49,6 @@ module constitutive
     constitutive_SandItsTangents, &
     constitutive_collectDotState, &
     constitutive_collectDeltaState, &
-    constitutive_postResults, &
     constitutive_results
  
 contains
@@ -61,17 +59,9 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine constitutive_init
 
-  integer, parameter :: FILEUNIT = 204
   integer :: &
-    o, &                                                                                            !< counter in output loop
     ph, &                                                                                           !< counter in phase loop
-    s, &                                                                                            !< counter in source loop
+    s                                                                                               !< counter in source loop
-    ins                                                                                             !< instance of plasticity/source
- 
-  integer, dimension(:,:), pointer :: thisSize
-  character(len=64), dimension(:,:), pointer :: thisOutput
-  character(len=32) :: outputName                                                                   !< name of output, intermediate fix until HDF5 output is ready
-  logical :: knownSource
 
 !--------------------------------------------------------------------------------------------------
 ! initialized plasticity
@@ -101,58 +91,10 @@ subroutine constitutive_init
   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
  
-  write(6,'(/,a)')   ' <<<+-  constitutive init  -+>>>'
+  write(6,'(/,a)')   ' <<<+-  constitutive init  -+>>>'; flush(6)
- 
-  mainProcess: if (worldrank == 0) then
-!--------------------------------------------------------------------------------------------------
-! write description file for constitutive output
-    call IO_write_jobFile(FILEUNIT,'outputConstitutive')
-    PhaseLoop: do ph = 1,material_Nphase
-      activePhase: if (any(material_phaseAt == ph)) then
-        write(FILEUNIT,'(/,a,/)') '['//trim(config_name_phase(ph))//']'
-        
-        SourceLoop: do s = 1, phase_Nsources(ph)
-          knownSource = .true.                                                                      ! assume valid
-          sourceType: select case (phase_source(s,ph))
-            case (SOURCE_damage_isoBrittle_ID) sourceType
-              ins = source_damage_isoBrittle_instance(ph)
-              outputName = SOURCE_damage_isoBrittle_label
-              thisOutput => source_damage_isoBrittle_output
-              thisSize   => source_damage_isoBrittle_sizePostResult
-            case (SOURCE_damage_isoDuctile_ID) sourceType
-              ins = source_damage_isoDuctile_instance(ph)
-              outputName = SOURCE_damage_isoDuctile_label
-              thisOutput => source_damage_isoDuctile_output
-              thisSize   => source_damage_isoDuctile_sizePostResult
-            case (SOURCE_damage_anisoBrittle_ID) sourceType
-              ins = source_damage_anisoBrittle_instance(ph)
-              outputName = SOURCE_damage_anisoBrittle_label
-              thisOutput => source_damage_anisoBrittle_output
-              thisSize   => source_damage_anisoBrittle_sizePostResult
-            case (SOURCE_damage_anisoDuctile_ID) sourceType
-              ins = source_damage_anisoDuctile_instance(ph)
-              outputName = SOURCE_damage_anisoDuctile_label
-              thisOutput => source_damage_anisoDuctile_output
-              thisSize   => source_damage_anisoDuctile_sizePostResult
-            case default sourceType
-              knownSource = .false.
-          end select sourceType
-          if (knownSource) then
-            write(FILEUNIT,'(a)') '(source)'//char(9)//trim(outputName)
-            OutputSourceLoop: do o = 1,size(thisOutput(:,ins))
-              if(len_trim(thisOutput(o,ins)) > 0) &
-                write(FILEUNIT,'(a,i4)') trim(thisOutput(o,ins))//char(9),thisSize(o,ins)
-            enddo OutputSourceLoop
-          endif
-        enddo SourceLoop
-      endif activePhase
-    enddo PhaseLoop
-    close(FILEUNIT)
-  endif mainProcess
  
   constitutive_plasticity_maxSizeDotState = 0
   constitutive_source_maxSizeDotState = 0
-  constitutive_source_maxSizePostResults = 0
  
   PhaseLoop2:do ph = 1,material_Nphase
 !--------------------------------------------------------------------------------------------------
@@ -169,11 +111,8 @@ subroutine constitutive_init
                                                      plasticState(ph)%sizeDotState)
     constitutive_source_maxSizeDotState        = max(constitutive_source_maxSizeDotState, &
                                                      maxval(sourceState(ph)%p(:)%sizeDotState))
-    constitutive_source_maxSizePostResults     = max(constitutive_source_maxSizePostResults, &
-                                                     maxval(sourceState(ph)%p(:)%sizePostResults))
   enddo PhaseLoop2
 
-
 end subroutine constitutive_init
 
 
@@ -639,51 +578,6 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
 end subroutine constitutive_collectDeltaState
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief returns array of constitutive results
-!--------------------------------------------------------------------------------------------------
-function constitutive_postResults(S, Fi, ipc, ip, el)
-
-  integer, intent(in) :: &
-    ipc, &                                                                                          !< component-ID of integration point
-    ip, &                                                                                           !< integration point
-    el                                                                                              !< element
-  real(pReal), dimension(sum(sourceState(material_phaseAt(ipc,el))%p(:)%sizePostResults)) :: &
-    constitutive_postResults
-  real(pReal),  intent(in), dimension(3,3) :: &
-    Fi                                                                                              !< intermediate deformation gradient
-  real(pReal),  intent(in), dimension(3,3) :: &
-    S                                                                                               !< 2nd Piola Kirchhoff stress
-  integer :: &
-    startPos, endPos
-  integer :: &
-    i, of, instance                                                                                 !< counter in source loop
-
-  constitutive_postResults = 0.0_pReal
-
-
-  endPos   = 0
-
-  SourceLoop: do i = 1, phase_Nsources(material_phaseAt(ipc,el))
-    startPos = endPos + 1
-    endPos = endPos + sourceState(material_phaseAt(ipc,el))%p(i)%sizePostResults
-    of = material_phasememberAt(ipc,ip,el)
-    sourceType: select case (phase_source(i,material_phaseAt(ipc,el)))
-      case (SOURCE_damage_isoBrittle_ID) sourceType
-        constitutive_postResults(startPos:endPos) = source_damage_isoBrittle_postResults(material_phaseAt(ipc,el),of)
-      case (SOURCE_damage_isoDuctile_ID) sourceType
-        constitutive_postResults(startPos:endPos) = source_damage_isoDuctile_postResults(material_phaseAt(ipc,el),of)
-      case (SOURCE_damage_anisoBrittle_ID) sourceType
-        constitutive_postResults(startPos:endPos) = source_damage_anisoBrittle_postResults(material_phaseAt(ipc,el),of)
-      case (SOURCE_damage_anisoDuctile_ID) sourceType
-        constitutive_postResults(startPos:endPos) = source_damage_anisoDuctile_postResults(material_phaseAt(ipc,el),of)
-    end select sourceType
-
-  enddo SourceLoop
-
-end function constitutive_postResults
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief writes constitutive results to HDF5 output file
 !--------------------------------------------------------------------------------------------------

src/crystallite.f90

@@ -108,7 +108,6 @@ module crystallite
     crystallite_stressTangent, &
     crystallite_orientations, &
     crystallite_push33ToRef, &
-    crystallite_postResults, &
     crystallite_results
 
 contains
@@ -119,7 +118,6 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine crystallite_init
  
-  integer, parameter :: FILEUNIT=434
   logical, dimension(:,:), allocatable :: devNull
   integer :: &
     c, &                                                                                            !< counter in integration point component loop
@@ -233,13 +231,6 @@ subroutine crystallite_init
 #endif
   enddo
 
-!--------------------------------------------------------------------------------------------------
-! write description file for crystallite output
-  if (worldrank == 0) then
-    call IO_write_jobFile(FILEUNIT,'outputCrystallite')
-    write(FILEUNIT,'(/,a,/)') '[not supported anymore]'
-    close(FILEUNIT)
-  endif
   call config_deallocate('material.config/phase')
 
 !--------------------------------------------------------------------------------------------------
@@ -732,37 +723,6 @@ function crystallite_push33ToRef(ipc,ip,el, tensor33)
 end function crystallite_push33ToRef
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief return results of particular grain
-!--------------------------------------------------------------------------------------------------
-function crystallite_postResults(ipc, ip, el)
-
-  integer, intent(in):: &
-    el, &                         !< element index
-    ip, &                         !< integration point index
-    ipc                           !< grain index
-
-  real(pReal), dimension(1+ &
-                         1+sum(sourceState(material_phaseAt(ipc,el))%p(:)%sizePostResults)) :: &
-    crystallite_postResults
-  integer :: &
-    c
-
-
-  crystallite_postResults = 0.0_pReal
-  crystallite_postResults(1) = 0.0_pReal                      ! header-like information (length)
-  c = 1
-
-  crystallite_postResults(c+1) = real(sum(sourceState(material_phaseAt(ipc,el))%p(:)%sizePostResults),pReal)  ! size of constitutive results
-  c = c + 1
-  if (size(crystallite_postResults)-c > 0) &
-    crystallite_postResults(c+1:size(crystallite_postResults)) = &
-       constitutive_postResults(crystallite_S(1:3,1:3,ipc,ip,el), crystallite_Fi(1:3,1:3,ipc,ip,el), &
-                                ipc, ip, el)
-
-end function crystallite_postResults
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief writes crystallite results to HDF5 output file
 !--------------------------------------------------------------------------------------------------

src/damage_local.f90

@@ -11,12 +11,11 @@ module damage_local
   use source_damage_isoDuctile
   use source_damage_anisoBrittle
   use source_damage_anisoDuctile
+  use results
 
   implicit none
   private
   
-  integer,                       dimension(:,:),         allocatable, target, public :: &
-    damage_local_sizePostResult
   character(len=64),             dimension(:,:),         allocatable, target, public :: &
     damage_local_output
   integer,                       dimension(:),           allocatable, target, public :: &
@@ -42,7 +41,7 @@ module damage_local
   public :: &
     damage_local_init, &
     damage_local_updateState, &
-    damage_local_postResults
+    damage_local_Results
 
 contains
 
@@ -66,7 +65,6 @@ subroutine damage_local_init
   maxNinstance = count(damage_type == DAMAGE_local_ID)
   if (maxNinstance == 0) return
   
-  allocate(damage_local_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0)
   allocate(damage_local_output         (maxval(homogenization_Noutput),maxNinstance))
            damage_local_output = ''
   allocate(damage_local_outputID       (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID)
@@ -90,7 +88,6 @@ subroutine damage_local_init
             case ('damage')
             damage_local_output(i,damage_typeInstance(h)) = outputs(i)
               damage_local_Noutput(instance) = damage_local_Noutput(instance) + 1
-             damage_local_sizePostResult(i,damage_typeInstance(h)) = 1
         prm%outputID = [prm%outputID , damage_ID]
            end select
       
@@ -106,7 +103,6 @@ subroutine damage_local_init
 ! allocate state arrays
       sizeState = 1
       damageState(homog)%sizeState = sizeState
-      damageState(homog)%sizePostResults = sum(damage_local_sizePostResult(:,instance))
       allocate(damageState(homog)%state0   (sizeState,NofMyHomog), source=damage_initialPhi(homog))
       allocate(damageState(homog)%subState0(sizeState,NofMyHomog), source=damage_initialPhi(homog))
       allocate(damageState(homog)%state    (sizeState,NofMyHomog), source=damage_initialPhi(homog))
@@ -211,35 +207,30 @@ end subroutine damage_local_getSourceAndItsTangent
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief return array of damage results
+!> @brief writes results to HDF5 output file
 !--------------------------------------------------------------------------------------------------
-function damage_local_postResults(ip,el)
+subroutine damage_local_results(homog,group)
 
-  integer, intent(in) :: &
+  integer,          intent(in) :: homog
-    ip, &                                                                                           !< integration point
+  character(len=*), intent(in) :: group
-    el                                                                                              !< element
+#if defined(PETSc) || defined(DAMASK_HDF5)  
-  real(pReal), dimension(sum(damage_local_sizePostResult(:,damage_typeInstance(material_homogenizationAt(el))))) :: &
+  integer :: o, instance
-    damage_local_postResults
   
-  integer :: instance, homog, offset, o, c
-    
-  homog     = material_homogenizationAt(el)
-  offset    = damageMapping(homog)%p(ip,el)
   instance  = damage_typeInstance(homog)
   associate(prm => param(instance))
-  c = 0
 
   outputsLoop: do o = 1,size(prm%outputID)
     select case(prm%outputID(o))
     
       case (damage_ID)
-         damage_local_postResults(c+1) = damage(homog)%p(offset)
+        call results_writeDataset(group,damage(homog)%p,'phi',&
-         c = c + 1
+                                  'damage indicator','-')
     end select
   enddo outputsLoop
-  
   end associate
+#endif
+
+end subroutine damage_local_results
 
-end function damage_local_postResults
 
 end module damage_local

src/damage_nonlocal.f90

@@ -14,12 +14,11 @@ module damage_nonlocal
   use source_damage_isoDuctile
   use source_damage_anisoBrittle
   use source_damage_anisoDuctile
+  use results
 
   implicit none
   private
   
-  integer,                       dimension(:,:), allocatable, target, public :: &
-    damage_nonlocal_sizePostResult
   character(len=64),             dimension(:,:), allocatable, target, public :: &
     damage_nonlocal_output
   integer,                       dimension(:),   allocatable, target, public :: &
@@ -45,7 +44,7 @@ module damage_nonlocal
     damage_nonlocal_getDiffusion33, &
     damage_nonlocal_getMobility, &
     damage_nonlocal_putNonLocalDamage, &
-    damage_nonlocal_postResults
+    damage_nonlocal_Results
 
 contains
 
@@ -55,7 +54,7 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine damage_nonlocal_init
 
-  integer :: maxNinstance,homog,instance,o,i
+  integer :: maxNinstance,homog,instance,i
   integer :: sizeState
   integer :: NofMyHomog, h
   integer(kind(undefined_ID)) :: &
@@ -69,7 +68,6 @@ subroutine damage_nonlocal_init
   maxNinstance = count(damage_type == DAMAGE_nonlocal_ID)
   if (maxNinstance == 0) return
   
-  allocate(damage_nonlocal_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0)
   allocate(damage_nonlocal_output         (maxval(homogenization_Noutput),maxNinstance))
            damage_nonlocal_output = ''
   allocate(damage_nonlocal_Noutput        (maxNinstance),                               source=0) 
@@ -92,7 +90,6 @@ subroutine damage_nonlocal_init
             case ('damage')
             damage_nonlocal_output(i,damage_typeInstance(h)) = outputs(i)
               damage_nonlocal_Noutput(instance) = damage_nonlocal_Noutput(instance) + 1
-             damage_nonlocal_sizePostResult(i,damage_typeInstance(h)) = 1
         prm%outputID = [prm%outputID , damage_ID]
            end select
       
@@ -107,7 +104,6 @@ subroutine damage_nonlocal_init
 !  allocate state arrays
       sizeState = 1
       damageState(homog)%sizeState = sizeState
-      damageState(homog)%sizePostResults = sum(damage_nonlocal_sizePostResult(:,instance))
       allocate(damageState(homog)%state0   (sizeState,NofMyHomog), source=damage_initialPhi(homog))
       allocate(damageState(homog)%subState0(sizeState,NofMyHomog), source=damage_initialPhi(homog))
       allocate(damageState(homog)%state    (sizeState,NofMyHomog), source=damage_initialPhi(homog))
@@ -247,35 +243,29 @@ end subroutine damage_nonlocal_putNonLocalDamage
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief return array of damage results
+!> @brief writes results to HDF5 output file
 !--------------------------------------------------------------------------------------------------
-function damage_nonlocal_postResults(ip,el)
+subroutine damage_nonlocal_results(homog,group)
 
-  integer,              intent(in) :: &
+  integer,          intent(in) :: homog
-    ip, &                                                                                           !< integration point
+  character(len=*), intent(in) :: group
-    el                                                                                              !< element
+#if defined(PETSc) || defined(DAMASK_HDF5)  
-  real(pReal), dimension(sum(damage_nonlocal_sizePostResult(:,damage_typeInstance(material_homogenizationAt(el))))) :: &
+  integer :: o, instance
-    damage_nonlocal_postResults
   
-  integer :: &
-    instance, homog, offset, o, c
-    
-  homog     = material_homogenizationAt(el)
-  offset    = damageMapping(homog)%p(ip,el)
   instance  = damage_typeInstance(homog)
   associate(prm => param(instance))
-  c = 0
 
   outputsLoop: do o = 1,size(prm%outputID)
     select case(prm%outputID(o))
     
       case (damage_ID)
-         damage_nonlocal_postResults(c+1) = damage(homog)%p(offset)
+        call results_writeDataset(group,damage(homog)%p,'phi',&
-         c = c + 1
+                                  'damage indicator','-')
     end select
   enddo outputsLoop
-
   end associate
-end function damage_nonlocal_postResults
+#endif
+
+end subroutine damage_nonlocal_results
 
 end module damage_nonlocal

src/grid/DAMASK_grid.f90

@@ -15,11 +15,7 @@ program DAMASK_spectral
  use config
  use debug
  use math
- use mesh_grid
  use CPFEM2
- use FEsolving
- use numerics
- use homogenization
  use material
  use spectral_utilities
  use grid_mech_spectral_basic
@@ -28,7 +24,6 @@ program DAMASK_spectral
  use grid_damage_spectral
  use grid_thermal_spectral
  use results
- use rotations
 
  implicit none
 
@@ -78,16 +73,9 @@ program DAMASK_spectral
  character(len=6)  :: loadcase_string
  character(len=1024) :: &
    incInfo
- type(rotation) :: R
  type(tLoadCase), allocatable, dimension(:) :: loadCases                                            !< array of all load cases
  type(tLoadCase) :: newLoadCase
  type(tSolutionState), allocatable, dimension(:) :: solres
- integer(MPI_OFFSET_KIND) :: fileOffset
- integer(MPI_OFFSET_KIND), dimension(:), allocatable :: outputSize
- integer, parameter :: maxByteOut = 2147483647-4096                                                 !< limit of one file output write https://trac.mpich.org/projects/mpich/ticket/1742
- integer, parameter :: maxRealOut = maxByteOut/pReal
- integer(pLongInt), dimension(2) :: outputIndex
- PetscErrorCode :: ierr
  procedure(grid_mech_spectral_basic_init), pointer :: &
    mech_init
  procedure(grid_mech_spectral_basic_forward), pointer :: &
@@ -105,7 +93,7 @@ program DAMASK_spectral
 !--------------------------------------------------------------------------------------------------
 ! init DAMASK (all modules)
  call CPFEM_initAll
- write(6,'(/,a)')   ' <<<+-  DAMASK_spectral init  -+>>>'
+ write(6,'(/,a)')   ' <<<+-  DAMASK_spectral init  -+>>>'; flush(6)
 
  write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, 2019'
  write(6,'(a)')   ' https://doi.org/10.1007/978-981-10-6855-3_80'
@@ -189,6 +177,7 @@ program DAMASK_spectral
      newLoadCase%ID(field) = FIELD_DAMAGE_ID
    endif damageActive
 
+   call newLoadCase%rot%fromEulers(real([0.0,0.0,0.0],pReal))
    readIn: do i = 1, chunkPos(1)
      select case (IO_lc(IO_stringValue(line,chunkPos,i)))
        case('fdot','dotf','l','f')                                                                  ! assign values for the deformation BC matrix
@@ -244,14 +233,13 @@ program DAMASK_spectral
          do j = 1, 3
            temp_valueVector(j) = IO_floatValue(line,chunkPos,i+k+j)
          enddo
-         call R%fromEulers(temp_valueVector(1:3),degrees=(l==1))
+         call newLoadCase%rot%fromEulers(temp_valueVector(1:3),degrees=(l==1))
-         newLoadCase%rotation = R%asMatrix()
        case('rotation','rot')                                                                       ! assign values for the rotation  matrix
          temp_valueVector = 0.0_pReal
          do j = 1, 9
            temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j)
          enddo
-         newLoadCase%rotation = math_9to33(temp_valueVector)
+         call newLoadCase%rot%fromMatrix(math_9to33(temp_valueVector))
      end select
    enddo readIn
 
@@ -282,10 +270,8 @@ program DAMASK_spectral
      enddo
      if (any(newLoadCase%stress%maskLogical .eqv. &
              newLoadCase%deformation%maskLogical)) errorID = 831                                    ! exclusive or masking only
-     if (any(newLoadCase%stress%maskLogical .and. &
+     if (any(newLoadCase%stress%maskLogical .and. transpose(newLoadCase%stress%maskLogical) &
-             transpose(newLoadCase%stress%maskLogical) .and. &
+         .and. (math_I3<1))) errorID = 838                                                          ! no rotation is allowed by stress BC
-             reshape([ .false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false.],[ 3,3]))) &
-             errorID = 838                                                                          ! no rotation is allowed by stress BC
      write(6,'(2x,a)') 'stress / GPa:'
      do i = 1, 3; do j = 1, 3
        if(newLoadCase%stress%maskLogical(i,j)) then
@@ -295,14 +281,12 @@ program DAMASK_spectral
        endif
        enddo; write(6,'(/)',advance='no')
      enddo
-    if (any(abs(matmul(newLoadCase%rotation, &
+     if (any(abs(matmul(newLoadCase%rot%asMatrix(), &
-                transpose(newLoadCase%rotation))-math_I3) > &
+                        transpose(newLoadCase%rot%asMatrix()))-math_I3) > &
-                reshape(spread(tol_math_check,1,9),[ 3,3]))&
+                reshape(spread(tol_math_check,1,9),[ 3,3]))) errorID = 846                          ! given rotation matrix contains strain
-                .or. abs(math_det33(newLoadCase%rotation)) > &
+     if (any(dNeq(newLoadCase%rot%asMatrix(), math_I3))) &
-                1.0_pReal + tol_math_check) errorID = 846                                           ! given rotation matrix contains strain
-     if (any(dNeq(newLoadCase%rotation, math_I3))) &
        write(6,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'rotation of loadframe:',&
-                transpose(newLoadCase%rotation)
+                transpose(newLoadCase%rot%asMatrix())
      if (newLoadCase%time < 0.0_pReal) errorID = 834                                                ! negative time increment
      write(6,'(2x,a,f12.6)') 'time:       ', newLoadCase%time
      if (newLoadCase%incs < 1)    errorID = 835                                                     ! non-positive incs count
@@ -339,26 +323,10 @@ program DAMASK_spectral
 ! write header of output file
  if (worldrank == 0) then
    writeHeader: if (interface_restartInc < 1) then
-     open(newunit=fileUnit,file=trim(getSolverJobName())//&
-                                 '.spectralOut',form='UNFORMATTED',status='REPLACE')
-     write(fileUnit) 'load:',       trim(loadCaseFile)                                              ! ... and write header
-     write(fileUnit) 'workingdir:', 'n/a'
-     write(fileUnit) 'geometry:',   trim(geometryFile)
-     write(fileUnit) 'grid:',       grid
-     write(fileUnit) 'size:',       geomSize
-     write(fileUnit) 'materialpoint_sizeResults:', materialpoint_sizeResults
-     write(fileUnit) 'loadcases:',  size(loadCases)
-     write(fileUnit) 'frequencies:', loadCases%outputfrequency                                      ! one entry per LoadCase
-     write(fileUnit) 'times:',      loadCases%time                                                  ! one entry per LoadCase
-     write(fileUnit) 'logscales:',  loadCases%logscale
-     write(fileUnit) 'increments:', loadCases%incs                                                  ! one entry per LoadCase
-     write(fileUnit) 'startingIncrement:', interface_restartInc                                     ! start with writing out the previous inc
-     write(fileUnit) 'eoh'
-     close(fileUnit)                                                                                ! end of header
      open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE')
      write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded'                 ! statistics file
      if (iand(debug_level(debug_spectral),debug_levelBasic) /= 0) &
-       write(6,'(/,a)') ' header of result and statistics file written out'
+       write(6,'(/,a)') ' header of statistics file written out'
      flush(6)
    else writeHeader
      open(newunit=statUnit,file=trim(getSolverJobName())//&
@@ -366,40 +334,11 @@ program DAMASK_spectral
    endif writeHeader
  endif
 
-!--------------------------------------------------------------------------------------------------
-! prepare MPI parallel out (including opening of file)
- allocate(outputSize(worldsize), source = 0_MPI_OFFSET_KIND)
- outputSize(worldrank+1) = size(materialpoint_results,kind=MPI_OFFSET_KIND)*int(pReal,MPI_OFFSET_KIND)
- call MPI_allreduce(MPI_IN_PLACE,outputSize,worldsize,MPI_LONG,MPI_SUM,PETSC_COMM_WORLD,ierr)       ! get total output size over each process
- if (ierr /= 0) call IO_error(error_ID=894, ext_msg='MPI_allreduce')
- call MPI_file_open(PETSC_COMM_WORLD, trim(getSolverJobName())//'.spectralOut', &
-                    MPI_MODE_WRONLY + MPI_MODE_APPEND, &
-                    MPI_INFO_NULL, &
-                    fileUnit, &
-                    ierr)
- if (ierr /= 0) call IO_error(error_ID=894, ext_msg='MPI_file_open')
- call MPI_file_get_position(fileUnit,fileOffset,ierr)                                               ! get offset from header
- if (ierr /= 0) call IO_error(error_ID=894, ext_msg='MPI_file_get_position')
- fileOffset = fileOffset + sum(outputSize(1:worldrank))                                             ! offset of my process in file (header + processes before me)
- call MPI_file_seek (fileUnit,fileOffset,MPI_SEEK_SET,ierr)
- if (ierr /= 0) call IO_error(error_ID=894, ext_msg='MPI_file_seek')
-
  writeUndeformed: if (interface_restartInc < 1) then
    write(6,'(1/,a)') ' ... writing initial configuration to file ........................'
    call CPFEM_results(0,0.0_pReal)
-   do i = 1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1         ! slice the output of my process in chunks not exceeding the limit for one output
-     outputIndex = int([(i-1)*((maxRealOut)/materialpoint_sizeResults)+1, &
-                             min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt)
-     call MPI_file_write(fileUnit,reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)), &
-                                 [(outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)]), &
-                         int((outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)), &
-                         MPI_DOUBLE, MPI_STATUS_IGNORE, ierr)
-     if (ierr /= 0) call IO_error(error_ID=894, ext_msg='MPI_file_write')
-   enddo
-   fileOffset = fileOffset + sum(outputSize)                                                        ! forward to current file position
  endif writeUndeformed
 
-
  loadCaseLooping: do currentLoadCase = 1, size(loadCases)
    time0 = time                                                                                     ! load case start time
    guess = loadCases(currentLoadCase)%followFormerTrajectory                                        ! change of load case? homogeneous guess for the first inc
@@ -443,19 +382,12 @@ program DAMASK_spectral
 !--------------------------------------------------------------------------------------------------
 ! report begin of new step
          write(6,'(/,a)') ' ###########################################################################'
-         write(6,'(1x,a,es12.5'//&
+         write(6,'(1x,a,es12.5,6(a,i0))') &
-                 ',a,'//IO_intOut(inc)            //',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//&
-                 ',a,'//IO_intOut(stepFraction)   //',a,'//IO_intOut(subStepFactor**cutBackLevel)//&
-                 ',a,'//IO_intOut(currentLoadCase)//',a,'//IO_intOut(size(loadCases))//')') &
                  'Time', time, &
                  's: Increment ', inc,'/',loadCases(currentLoadCase)%incs,&
                  '-', stepFraction,'/',subStepFactor**cutBackLevel,&
                  ' of load case ', currentLoadCase,'/',size(loadCases)
-         write(incInfo,&
+         write(incInfo,'(4(a,i0))') &
-                 '(a,'//IO_intOut(totalIncsCounter)//&
-                 ',a,'//IO_intOut(sum(loadCases%incs))//&
-                 ',a,'//IO_intOut(stepFraction)//&
-                 ',a,'//IO_intOut(subStepFactor**cutBackLevel)//')') &
                  'Increment ',totalIncsCounter,'/',sum(loadCases%incs),&
                  '-', stepFraction,'/',subStepFactor**cutBackLevel
          flush(6)
@@ -469,7 +401,7 @@ program DAMASK_spectral
                        cutBack,guess,timeinc,timeIncOld,remainingLoadCaseTime, &
                        deformation_BC     = loadCases(currentLoadCase)%deformation, &
                        stress_BC          = loadCases(currentLoadCase)%stress, &
-                       rotation_BC        = loadCases(currentLoadCase)%rotation)
+                       rotation_BC        = loadCases(currentLoadCase)%rot)
 
              case(FIELD_THERMAL_ID); call grid_thermal_spectral_forward(cutBack)
              case(FIELD_DAMAGE_ID);  call grid_damage_spectral_forward(cutBack)
@@ -488,7 +420,7 @@ program DAMASK_spectral
                  solres(field) = mech_solution (&
                                         incInfo,timeinc,timeIncOld, &
                                         stress_BC   = loadCases(currentLoadCase)%stress, &
-                                        rotation_BC = loadCases(currentLoadCase)%rotation)
+                                        rotation_BC = loadCases(currentLoadCase)%rot)
 
                case(FIELD_THERMAL_ID)
                  solres(field) = grid_thermal_spectral_solution(timeinc,timeIncOld)
@@ -530,7 +462,6 @@ program DAMASK_spectral
            write(6,'(/,a)') ' cutting back '
          else                                                                                       ! no more options to continue
            call IO_warning(850)
-           call MPI_File_close(fileUnit,ierr)
            close(statUnit)
            call quit(0)                                                                             ! quit
          endif
@@ -540,29 +471,14 @@ program DAMASK_spectral
        cutBackLevel = max(0, cutBackLevel - 1)                                                      ! try half number of subincs next inc
 
        if (all(solres(:)%converged)) then
-         write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') &                                     ! report converged inc
+         write(6,'(/,a,i0,a)') ' increment ', totalIncsCounter, ' converged'
-                                   ' increment ', totalIncsCounter, ' converged'
        else
-         write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') &                                     ! report non-converged inc
+         write(6,'(/,a,i0,a)') ' increment ', totalIncsCounter, ' NOT converged'
-                                   ' increment ', totalIncsCounter, ' NOT converged'
        endif; flush(6)
 
        if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then                           ! at output frequency
          write(6,'(1/,a)') ' ... writing results to file ......................................'
          flush(6)
-         call materialpoint_postResults()
-         call MPI_File_seek (fileUnit,fileOffset,MPI_SEEK_SET,ierr)
-         if (ierr /= 0) call IO_error(894, ext_msg='MPI_file_seek')
-         do i=1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1     ! slice the output of my process in chunks not exceeding the limit for one output
-           outputIndex=int([(i-1)*((maxRealOut)/materialpoint_sizeResults)+1, &
-                      min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt)
-           call MPI_file_write(fileUnit,reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)),&
-                                       [(outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)]), &
-                               int((outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)),&
-                               MPI_DOUBLE, MPI_STATUS_IGNORE, ierr)
-           if(ierr /=0) call IO_error(894, ext_msg='MPI_file_write')
-         enddo
-         fileOffset = fileOffset + sum(outputSize)                                                  ! forward to current file position
          call CPFEM_results(totalIncsCounter,time)
        endif
        if (mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0) then
@@ -579,7 +495,6 @@ program DAMASK_spectral
 !--------------------------------------------------------------------------------------------------
 ! report summary of whole calculation
  write(6,'(/,a)') ' ###########################################################################'
- call MPI_file_close(fileUnit,ierr)
  close(statUnit)
 
  call quit(0)                                                                                       ! no complains ;)

src/grid/grid_mech_FEM.f90

@@ -98,12 +98,12 @@ subroutine grid_mech_FEM_init
   PetscErrorCode :: ierr
   integer        :: rank
   integer(HID_T) :: fileHandle, groupHandle
-  character(len=1024) :: rankStr
+  character(len=pStringLen) :: fileName
   real(pReal), dimension(3,3,3,3) :: devNull
   PetscScalar, pointer, dimension(:,:,:,:) :: &
   u_current,u_lastInc
  
-  write(6,'(/,a)') ' <<<+-  grid_mech_FEM init  -+>>>'
+  write(6,'(/,a)') ' <<<+-  grid_mech_FEM init  -+>>>'; flush(6)
   
 !--------------------------------------------------------------------------------------------------
 ! set default and user defined options for PETSc
@@ -184,11 +184,10 @@ subroutine grid_mech_FEM_init
 !--------------------------------------------------------------------------------------------------
 ! init fields
   restartRead: if (interface_restartInc > 0) then
-    write(6,'(/,a,'//IO_intOut(interface_restartInc)//',a)') &
+    write(6,'(/,a,i0,a)') ' reading restart data of increment ', interface_restartInc, ' from file'
-      'reading values of increment ', interface_restartInc, ' from file'
  
-    write(rankStr,'(a1,i0)')'_',worldrank
+    write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5'
-    fileHandle  = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5')
+    fileHandle  = HDF5_openFile(fileName)
     groupHandle = HDF5_openGroup(fileHandle,'solver')
     
     call HDF5_read(groupHandle,F_aim,        'F_aim')
@@ -207,16 +206,14 @@ subroutine grid_mech_FEM_init
   call utilities_updateCoords(F)
   call utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, &                     ! stress field, stress avg, global average of stiffness and (min+max)/2
                                       F, &                                                          ! target F
-                                      0.0_pReal, &                                                  ! time increment
+                                      0.0_pReal)                                                    ! time increment
-                                      math_I3)                                                      ! no rotation of boundary condition
   call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr)
   CHKERRQ(ierr)
   call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr)
   CHKERRQ(ierr)
  
   restartRead2: if (interface_restartInc > 0) then
-    write(6,'(/,a,'//IO_intOut(interface_restartInc)//',a)') &
+    write(6,'(/,a,i0,a)') ' reading more restart data of increment ', interface_restartInc, ' from file'
-      'reading more values of increment ', interface_restartInc, ' from file'
     call HDF5_read(groupHandle,C_volAvg,       'C_volAvg')
     call HDF5_read(groupHandle,C_volAvgLastInc,'C_volAvgLastInc')
     
@@ -242,7 +239,8 @@ function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation
     timeinc_old                                                                                     !< time increment of last successful increment
   type(tBoundaryCondition),    intent(in) :: &
     stress_BC
-  real(pReal), dimension(3,3), intent(in) :: rotation_BC
+  type(rotation),              intent(in) :: &
+    rotation_BC
   type(tSolutionState)                    :: &
     solution
 !--------------------------------------------------------------------------------------------------
@@ -297,7 +295,7 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,
   type(tBoundaryCondition),    intent(in) :: &
     stress_BC, &
     deformation_BC
-  real(pReal), dimension(3,3), intent(in) :: &
+  type(rotation),              intent(in) :: &
     rotation_BC
   PetscErrorCode :: ierr
   PetscScalar, pointer, dimension(:,:,:,:) :: &
@@ -355,7 +353,7 @@ end subroutine grid_mech_FEM_forward
 !--------------------------------------------------------------------------------------------------
 !> @brief Age
 !--------------------------------------------------------------------------------------------------
-subroutine grid_mech_FEM_updateCoords()
+subroutine grid_mech_FEM_updateCoords
 
   call utilities_updateCoords(F)
 
@@ -365,20 +363,20 @@ end subroutine grid_mech_FEM_updateCoords
 !--------------------------------------------------------------------------------------------------
 !> @brief Write current solver and constitutive data for restart to file
 !--------------------------------------------------------------------------------------------------
-subroutine grid_mech_FEM_restartWrite()
+subroutine grid_mech_FEM_restartWrite
 
   PetscErrorCode :: ierr
-  PetscScalar, dimension(:,:,:,:), pointer :: u_current,u_lastInc
   integer(HID_T) :: fileHandle, groupHandle
-  character(len=32) :: rankStr
+  PetscScalar, dimension(:,:,:,:), pointer :: u_current,u_lastInc
+  character(len=pStringLen) :: fileName
   
   call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr)
   call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr)
 
-  write(6,'(a)') ' writing solver data required for restart to file';flush(6)
+  write(6,'(a)') ' writing solver data required for restart to file'; flush(6)
   
-  write(rankStr,'(a1,i0)')'_',worldrank
+  write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5'
-  fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w')
+  fileHandle  = HDF5_openFile(fileName,'w')
   groupHandle = HDF5_addGroup(fileHandle,'solver')
 
   call HDF5_write(groupHandle,F_aim,        'F_aim')
@@ -478,11 +476,10 @@ subroutine formResidual(da_local,x_local, &
 ! begin of new iteration
   newIteration: if (totalIter <= PETScIter) then
     totalIter = totalIter + 1
-    write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') &
+    write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter+1, '≤', itmax
-            trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter+1, '≤', itmax
     if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
       write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
-              ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC))
+              ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotTensor2(F_aim,active=.true.))
     write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
               ' deformation gradient aim       =', transpose(F_aim)
     flush(6)
@@ -498,7 +495,7 @@ subroutine formResidual(da_local,x_local, &
       x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk)
     enddo; enddo; enddo
     ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1
-    F(1:3,1:3,ii,jj,kk) = math_rotate_backward33(F_aim,params%rotation_BC) + transpose(matmul(BMat,x_elem)) 
+    F(1:3,1:3,ii,jj,kk) = params%rotation_BC%rotTensor2(F_aim,active=.true.) + transpose(matmul(BMat,x_elem)) 
   enddo; enddo; enddo
   call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr)
 

src/grid/grid_mech_spectral_basic.f90

@@ -96,9 +96,9 @@ subroutine grid_mech_spectral_basic_init
   PetscInt, dimension(worldsize) :: localK  
   integer(HID_T) :: fileHandle, groupHandle
   integer        :: fileUnit
-  character(len=1024) :: rankStr
+  character(len=pStringLen) :: fileName
  
-  write(6,'(/,a)') ' <<<+-  grid_mech_spectral_basic init  -+>>>'
+  write(6,'(/,a)') ' <<<+-  grid_mech_spectral_basic init  -+>>>'; flush(6)
  
   write(6,'(/,a)') ' Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013'
   write(6,'(a)')   ' https://doi.org/10.1016/j.ijplas.2012.09.012'
@@ -151,11 +151,10 @@ subroutine grid_mech_spectral_basic_init
   call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)                                   ! places pointer on PETSc data
  
   restartRead: if (interface_restartInc > 0) then
-    write(6,'(/,a,'//IO_intOut(interface_restartInc)//',a)') &
+    write(6,'(/,a,i0,a)') ' reading restart data of increment ', interface_restartInc, ' from file'
-      ' reading values of increment ', interface_restartInc, ' from file'
 
-    write(rankStr,'(a1,i0)')'_',worldrank
+    write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5'
-    fileHandle  = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5')
+    fileHandle  = HDF5_openFile(fileName)
     groupHandle = HDF5_openGroup(fileHandle,'solver')
  
     call HDF5_read(groupHandle,F_aim,        'F_aim')
@@ -173,13 +172,11 @@ subroutine grid_mech_spectral_basic_init
   call Utilities_updateCoords(reshape(F,shape(F_lastInc)))
   call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, &                         ! stress field, stress avg, global average of stiffness and (min+max)/2
                                       reshape(F,shape(F_lastInc)), &                                ! target F
-                                      0.0_pReal, &                                                  ! time increment
+                                      0.0_pReal)                                                    ! time increment
-                                      math_I3)                                                      ! no rotation of boundary condition
   call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)                                ! deassociate pointer
  
   restartRead2: if (interface_restartInc > 0) then
-    write(6,'(/,a,'//IO_intOut(interface_restartInc)//',a)') &
+    write(6,'(/,a,i0,a)') ' reading more restart data of increment ', interface_restartInc, ' from file'
-      'reading more values of increment ', interface_restartInc, ' from file'
     call HDF5_read(groupHandle,C_volAvg,       'C_volAvg')
     call HDF5_read(groupHandle,C_volAvgLastInc,'C_volAvgLastInc')
     
@@ -212,7 +209,8 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_
     timeinc_old                                                                                     !< time increment of last successful increment
   type(tBoundaryCondition),    intent(in) :: &
     stress_BC
-  real(pReal), dimension(3,3), intent(in) :: rotation_BC
+  type(rotation),              intent(in) :: &
+    rotation_BC
   type(tSolutionState)                    :: &
     solution
 !--------------------------------------------------------------------------------------------------
@@ -269,7 +267,7 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo
   type(tBoundaryCondition),    intent(in) :: &
     stress_BC, &
     deformation_BC
-  real(pReal), dimension(3,3), intent(in) :: &
+  type(rotation),              intent(in) :: &
     rotation_BC
   PetscErrorCode :: ierr
   PetscScalar, dimension(:,:,:,:), pointer :: F
@@ -301,7 +299,7 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo
 
     Fdot = utilities_calculateRate(guess, &
                                    F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, &
-                                    math_rotate_backward33(F_aimDot,rotation_BC))
+                                   rotation_BC%rotTensor2(F_aimDot,active=.true.))
     F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3])
     
     materialpoint_F0 = reshape(F, [3,3,1,product(grid(1:2))*grid3])
@@ -311,7 +309,7 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo
 ! update average and local deformation gradients
   F_aim = F_aim_lastInc + F_aimDot * timeinc
   F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, &                                       ! estimate of F at end of time+timeinc that matches rotated F_aim on average
-              math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid3])
+              rotation_BC%rotTensor2(F_aim,active=.true.)),[9,grid(1),grid(2),grid3])
   call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
   
 end subroutine grid_mech_spectral_basic_forward
@@ -320,7 +318,7 @@ end subroutine grid_mech_spectral_basic_forward
 !--------------------------------------------------------------------------------------------------
 !> @brief Age
 !--------------------------------------------------------------------------------------------------
-subroutine grid_mech_spectral_basic_updateCoords()
+subroutine grid_mech_spectral_basic_updateCoords
 
   PetscErrorCode :: ierr
   PetscScalar, dimension(:,:,:,:), pointer :: F
@@ -335,19 +333,19 @@ end subroutine grid_mech_spectral_basic_updateCoords
 !--------------------------------------------------------------------------------------------------
 !> @brief Write current solver and constitutive data for restart to file
 !--------------------------------------------------------------------------------------------------
-subroutine grid_mech_spectral_basic_restartWrite()
+subroutine grid_mech_spectral_basic_restartWrite
 
   PetscErrorCode :: ierr
-  PetscScalar, dimension(:,:,:,:), pointer :: F
   integer(HID_T) :: fileHandle, groupHandle
-  character(len=32) :: rankStr
+  PetscScalar, dimension(:,:,:,:), pointer :: F
+  character(len=pStringLen) :: fileName
 
   call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
 
-  write(6,'(a)') ' writing solver data required for restart to file';flush(6)
+  write(6,'(a)') ' writing solver data required for restart to file'; flush(6)
   
-  write(rankStr,'(a1,i0)')'_',worldrank
+  write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5'
-  fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w')
+  fileHandle  = HDF5_openFile(fileName,'w')
   groupHandle = HDF5_addGroup(fileHandle,'solver')
   
   call HDF5_write(groupHandle,F_aim,        'F_aim')
@@ -442,11 +440,10 @@ subroutine formResidual(in, F, &
 ! begin of new iteration
   newIteration: if (totalIter <= PETScIter) then
     totalIter = totalIter + 1
-    write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') &
+    write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
-            trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
     if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
       write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
-              ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC))
+              ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotTensor2(F_aim,active=.true.))
     write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
               ' deformation gradient aim       =', transpose(F_aim)
     flush(6)
@@ -471,7 +468,7 @@ subroutine formResidual(in, F, &
   tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residuum                                  ! store fPK field for subsequent FFT forward transform
   call utilities_FFTtensorForward                                                                   ! FFT forward of global "tensorField_real"
   err_div = Utilities_divergenceRMS()                                                               ! divRMS of tensorField_fourier for later use
-  call utilities_fourierGammaConvolution(math_rotate_backward33(deltaF_aim,params%rotation_BC))     ! convolution of Gamma and tensorField_fourier, with arg 
+  call utilities_fourierGammaConvolution(params%rotation_BC%rotTensor2(deltaF_aim,active=.true.))   ! convolution of Gamma and tensorField_fourier
   call utilities_FFTtensorBackward                                                                  ! FFT backward of global tensorField_fourier
  
 !--------------------------------------------------------------------------------------------------

src/grid/grid_mech_spectral_polarisation.f90

@@ -14,6 +14,7 @@ module grid_mech_spectral_polarisation
   use DAMASK_interface
   use HDF5_utilities
   use math
+  use rotations
   use spectral_utilities
   use IO
   use FEsolving
@@ -104,9 +105,9 @@ subroutine grid_mech_spectral_polarisation_init
   PetscInt, dimension(worldsize) :: localK 
   integer(HID_T) :: fileHandle, groupHandle
   integer        :: fileUnit
-  character(len=1024) :: rankStr
+  character(len=pStringLen) :: fileName
   
-  write(6,'(/,a)') ' <<<+-  grid_mech_spectral_polarisation init  -+>>>'
+  write(6,'(/,a)') ' <<<+-  grid_mech_spectral_polarisation init  -+>>>'; flush(6)
  
   write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
   write(6,'(a)')   ' https://doi.org/10.1016/j.ijplas.2014.02.006'
@@ -160,11 +161,10 @@ subroutine grid_mech_spectral_polarisation_init
   F_tau => FandF_tau(9:17,:,:,:)
  
   restartRead: if (interface_restartInc > 0) then
-    write(6,'(/,a,'//IO_intOut(interface_restartInc)//',a)') &
+    write(6,'(/,a,i0,a)') ' reading restart data of increment ', interface_restartInc, ' from file'
-      ' reading values of increment ', interface_restartInc, ' from file'
 
-    write(rankStr,'(a1,i0)')'_',worldrank
+    write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5'
-    fileHandle  = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5')
+    fileHandle  = HDF5_openFile(fileName)
     groupHandle = HDF5_openGroup(fileHandle,'solver')
  
     call HDF5_read(groupHandle,F_aim,        'F_aim')
@@ -186,13 +186,11 @@ subroutine grid_mech_spectral_polarisation_init
   call Utilities_updateCoords(reshape(F,shape(F_lastInc)))
   call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, &                         ! stress field, stress avg, global average of stiffness and (min+max)/2
                                       reshape(F,shape(F_lastInc)), &                                ! target F
-                                      0.0_pReal, &                                                  ! time increment
+                                      0.0_pReal)                                                    ! time increment
-                                      math_I3)                                                      ! no rotation of boundary condition
   call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)                        ! deassociate pointer
  
   restartRead2: if (interface_restartInc > 0) then
-    write(6,'(/,a,'//IO_intOut(interface_restartInc)//',a)') &
+    write(6,'(/,a,i0,a)') ' reading more restart data of increment ', interface_restartInc, ' from file'
-      ' reading more values of increment ', interface_restartInc, ' from file'
     call HDF5_read(groupHandle,C_volAvg,       'C_volAvg')
     call HDF5_read(groupHandle,C_volAvgLastInc,'C_volAvgLastInc')
 
@@ -227,7 +225,8 @@ function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,
     timeinc_old                                                                                     !< time increment of last successful increment
   type(tBoundaryCondition),    intent(in) :: &
     stress_BC
-  real(pReal), dimension(3,3), intent(in) :: rotation_BC
+  type(rotation),              intent(in) :: &
+    rotation_BC
   type(tSolutionState)                    :: &
     solution
 !--------------------------------------------------------------------------------------------------
@@ -288,7 +287,7 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
   type(tBoundaryCondition),      intent(in) :: &
     stress_BC, &
     deformation_BC
-  real(pReal), dimension(3,3), intent(in) ::&
+  type(rotation), intent(in) :: &
     rotation_BC
   PetscErrorCode :: ierr
   PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau, F, F_tau
@@ -324,10 +323,10 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
 
     Fdot     = utilities_calculateRate(guess, &
                                        F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, &
-                                       math_rotate_backward33(F_aimDot,rotation_BC))
+                                       rotation_BC%rotTensor2(F_aimDot,active=.true.))
     F_tauDot = utilities_calculateRate(guess, &
                                        F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid3]), timeinc_old, &
-                                       math_rotate_backward33(F_aimDot,rotation_BC))
+                                       rotation_BC%rotTensor2(F_aimDot,active=.true.))
     F_lastInc     = reshape(F,    [3,3,grid(1),grid(2),grid3])
     F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3])
     
@@ -338,7 +337,7 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
 ! update average and local deformation gradients
   F_aim = F_aim_lastInc + F_aimDot * timeinc
   F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, &                                      ! estimate of F at end of time+timeinc that matches rotated F_aim on average
-                                     math_rotate_backward33(F_aim,rotation_BC)),&
+                                     rotation_BC%rotTensor2(F_aim,active=.true.)),&
               [9,grid(1),grid(2),grid3])
   if (guess) then
      F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), &
@@ -349,7 +348,7 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
        F_lambda33 = math_mul3333xx33(S_scale,matmul(F_lambda33, &
                                    math_mul3333xx33(C_scale,&
                                                     matmul(transpose(F_lambda33),&
-                                                                  F_lambda33)-math_I3))*0.5_pReal)&
+                                                           F_lambda33)-math_I3))*0.5_pReal) &
                   + math_I3
        F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k)
     enddo; enddo; enddo
@@ -363,7 +362,7 @@ end subroutine grid_mech_spectral_polarisation_forward
 !--------------------------------------------------------------------------------------------------
 !> @brief Age
 !--------------------------------------------------------------------------------------------------
-subroutine grid_mech_spectral_polarisation_updateCoords()
+subroutine grid_mech_spectral_polarisation_updateCoords
 
   PetscErrorCode :: ierr
   PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau
@@ -378,21 +377,21 @@ end subroutine grid_mech_spectral_polarisation_updateCoords
 !--------------------------------------------------------------------------------------------------
 !> @brief Write current solver and constitutive data for restart to file
 !--------------------------------------------------------------------------------------------------
-subroutine grid_mech_spectral_polarisation_restartWrite()
+subroutine grid_mech_spectral_polarisation_restartWrite
 
   PetscErrorCode :: ierr
-  PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau, F, F_tau
   integer(HID_T) :: fileHandle, groupHandle
-  character(len=32) :: rankStr
+  PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau, F, F_tau
+  character(len=pStringLen) :: fileName
 
   call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
   F     => FandF_tau(0: 8,:,:,:)
   F_tau => FandF_tau(9:17,:,:,:)
 
-  write(6,'(a)') ' writing solver data required for restart to file';flush(6)
+  write(6,'(a)') ' writing solver data required for restart to file'; flush(6)
 
-  write(rankStr,'(a1,i0)')'_',worldrank
+  write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5'
-  fileHandle  = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w')
+  fileHandle  = HDF5_openFile(fileName,'w')
   groupHandle = HDF5_addGroup(fileHandle,'solver')
   
   call HDF5_write(groupHandle,F_aim,        'F_aim')
@@ -510,11 +509,10 @@ subroutine formResidual(in, FandF_tau, &
 ! begin of new iteration
   newIteration: if (totalIter <= PETScIter) then
     totalIter = totalIter + 1
-    write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') &
+    write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
-            trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
     if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
       write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
-              ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC))
+              ' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotTensor2(F_aim,active=.true.))
     write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
               ' deformation gradient aim       =', transpose(F_aim)
     flush(6)
@@ -533,7 +531,7 @@ subroutine formResidual(in, FandF_tau, &
 !--------------------------------------------------------------------------------------------------
 ! doing convolution in Fourier space 
   call utilities_FFTtensorForward
-  call utilities_fourierGammaConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) 
+  call utilities_fourierGammaConvolution(params%rotation_BC%rotTensor2(polarBeta*F_aim,active=.true.)) 
   call utilities_FFTtensorBackward
 
 !--------------------------------------------------------------------------------------------------
@@ -551,7 +549,7 @@ subroutine formResidual(in, FandF_tau, &
 ! stress BC handling
   F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC)))                                  ! S = 0.0 for no bc
   err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim &
-                                                 -math_rotate_forward33(F_av,params%rotation_BC)) + &
+                                                 -params%rotation_BC%rotTensor2(F_av)) + &
                                  params%stress_mask  * (P_av-params%stress_BC)))                    ! mask = 0.0 for no bc
 ! calculate divergence
   tensorField_real = 0.0_pReal

src/grid/spectral_utilities.f90

@@ -10,6 +10,7 @@ module spectral_utilities
 
   use prec
   use math
+  use rotations
   use IO
   use mesh_grid
   use numerics
@@ -90,7 +91,7 @@ module spectral_utilities
   end type tBoundaryCondition
 
   type, public :: tLoadCase
-    real(pReal), dimension (3,3) :: rotation               = math_I3                                !< rotation of BC
+    type(rotation)               :: rot                                                             !< rotation of BC
     type(tBoundaryCondition) ::     stress, &                                                       !< stress BC
                                     deformation                                                     !< deformation BC (Fdot or L)
     real(pReal) ::                  time                   = 0.0_pReal                              !< length of increment
@@ -103,7 +104,8 @@ module spectral_utilities
   end type tLoadCase
 
   type, public :: tSolutionParams                                                                   !< @todo use here the type definition for a full loadcase
-    real(pReal), dimension(3,3) :: stress_mask, stress_BC, rotation_BC
+    real(pReal), dimension(3,3) :: stress_mask, stress_BC
+    type(rotation)              :: rotation_BC
     real(pReal) :: timeinc
     real(pReal) :: timeincOld
   end type tSolutionParams
@@ -685,9 +687,10 @@ end function utilities_curlRMS
 function utilities_maskedCompliance(rot_BC,mask_stress,C)
 
   real(pReal),                dimension(3,3,3,3) :: utilities_maskedCompliance                      !< masked compliance
-  real(pReal), intent(in) , dimension(3,3,3,3) :: C                                                 !< current average stiffness
+  real(pReal),    intent(in), dimension(3,3,3,3) :: C                                               !< current average stiffness
-  real(pReal), intent(in) , dimension(3,3)     :: rot_BC                                            !< rotation of load frame
+  type(rotation), intent(in)                     :: rot_BC                                          !< rotation of load frame
   logical,        intent(in), dimension(3,3)     :: mask_stress                                     !< mask of stress BC
+
   integer :: j, k, m, n
   logical, dimension(9) :: mask_stressVector
   real(pReal), dimension(9,9) :: temp99_Real
@@ -705,7 +708,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
     allocate (c_reduced(size_reduced,size_reduced), source =0.0_pReal)
     allocate (s_reduced(size_reduced,size_reduced), source =0.0_pReal)
     allocate (sTimesC(size_reduced,size_reduced),   source =0.0_pReal)
-    temp99_Real = math_3333to99(math_rotate_forward3333(C,rot_BC))
+    temp99_Real = math_3333to99(rot_BC%rotTensor4(C))
  
     if(debugGeneral) then
       write(6,'(/,a)') ' ... updating masked compliance ............................................'
@@ -834,12 +837,12 @@ end subroutine utilities_fourierTensorDivergence
 subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
                                           F,timeinc,rotation_BC)
  
-  real(pReal),intent(out), dimension(3,3,3,3)                     :: C_volAvg, C_minmaxAvg          !< average stiffness
+  real(pReal),    intent(out), dimension(3,3,3,3)                   :: C_volAvg, C_minmaxAvg        !< average stiffness
-  real(pReal),intent(out), dimension(3,3)                         :: P_av                           !< average PK stress
+  real(pReal),    intent(out), dimension(3,3)                       :: P_av                         !< average PK stress
-  real(pReal),intent(out), dimension(3,3,grid(1),grid(2),grid3)   :: P                              !< PK stress
+  real(pReal),    intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P                            !< PK stress
   real(pReal),    intent(in),  dimension(3,3,grid(1),grid(2),grid3) :: F                            !< deformation gradient target
   real(pReal),    intent(in)                                        :: timeinc                      !< loading time
-  real(pReal), intent(in), dimension(3,3)                         :: rotation_BC                    !< rotation of load frame
+  type(rotation), intent(in),  optional                             :: rotation_BC                  !< rotation of load frame
  
   
   integer :: &
@@ -861,7 +864,8 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
   if (debugRotation) &
   write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress (lab) / MPa =',&
                                                       transpose(P_av)*1.e-6_pReal
-  P_av = math_rotate_forward33(P_av,rotation_BC)
+  if(present(rotation_BC)) &
+    P_av = rotation_BC%rotTensor2(P_av)
   write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress       / MPa =',&
                                                       transpose(P_av)*1.e-6_pReal
   flush(6)
@@ -1119,7 +1123,7 @@ subroutine utilities_saveReferenceStiffness
     fileUnit
 
   if (worldrank == 0) then
-    write(6,'(a)') ' writing reference stiffness data required for restart to file';flush(6)
+    write(6,'(a)') ' writing reference stiffness data required for restart to file'; flush(6)
     fileUnit = IO_open_jobFile_binary('C_ref','w')
     write(fileUnit) C_ref
     close(fileUnit)

src/homogenization.f90

@@ -23,7 +23,6 @@ module homogenization
   use damage_local
   use damage_nonlocal
   use results
-  use HDF5_utilities
  
   implicit none
   private
@@ -36,12 +35,6 @@ module homogenization
     materialpoint_P                                                                                 !< first P--K stress of IP
   real(pReal),   dimension(:,:,:,:,:,:), allocatable, public ::  &
     materialpoint_dPdF                                                                              !< tangent of first P--K stress at IP
-  real(pReal),   dimension(:,:,:),       allocatable, public :: &
-    materialpoint_results                                                                           !< results array of material point
-  integer,                                            public, protected  :: &
-    materialpoint_sizeResults, &
-    thermal_maxSizePostResults, &
-    damage_maxSizePostResults
 
   real(pReal),   dimension(:,:,:,:),     allocatable :: &
     materialpoint_subF0, &                                                                          !< def grad of IP at beginning of homogenization increment
@@ -126,7 +119,6 @@ module homogenization
   public ::  &
     homogenization_init, &
     materialpoint_stressAndItsTangent, &
-    materialpoint_postResults, &
     homogenization_results
 
 contains
@@ -137,14 +129,6 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine homogenization_init
 
-  integer, parameter :: FILEUNIT = 200
-  integer :: e,i,p
-  integer, dimension(:,:), pointer :: thisSize
-  integer, dimension(:)  , pointer :: thisNoutput
-  character(len=64), dimension(:,:), pointer :: thisOutput
-  character(len=32) :: outputName                                                                   !< name of output, intermediate fix until HDF5 output is ready
-  logical :: valid
-
   if (any(homogenization_type == HOMOGENIZATION_NONE_ID))      call mech_none_init
   if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call mech_isostrain_init
   if (any(homogenization_type == HOMOGENIZATION_RGC_ID))       call mech_RGC_init
@@ -157,80 +141,6 @@ subroutine homogenization_init
   if (any(damage_type == DAMAGE_local_ID))     call damage_local_init
   if (any(damage_type == DAMAGE_nonlocal_ID))  call damage_nonlocal_init
 
-!--------------------------------------------------------------------------------------------------
-! write description file for homogenization output
-  mainProcess: if (worldrank == 0) then
-    call IO_write_jobFile(FILEUNIT,'outputHomogenization')
-    do p = 1,size(config_homogenization)
-      if (any(material_homogenizationAt == p)) then
-        write(FILEUNIT,'(/,a,/)')  '['//trim(config_name_homogenization(p))//']'
-        write(FILEUNIT,'(a)') '(type) n/a'
-        write(FILEUNIT,'(a,i4)') '(ngrains)'//char(9),homogenization_Ngrains(p)
-        
-        i = thermal_typeInstance(p)                                                                 ! which instance of this thermal type
-        valid = .true.                                                                              ! assume valid
-        select case(thermal_type(p))                                                                ! split per thermal type
-          case (THERMAL_isothermal_ID)
-            outputName = THERMAL_isothermal_label
-            thisNoutput => null()
-            thisOutput => null()
-            thisSize   => null()
-          case (THERMAL_adiabatic_ID)
-            outputName = THERMAL_adiabatic_label
-            thisNoutput => thermal_adiabatic_Noutput
-            thisOutput => thermal_adiabatic_output
-            thisSize   => thermal_adiabatic_sizePostResult
-          case (THERMAL_conduction_ID)
-            outputName = THERMAL_conduction_label
-            thisNoutput => thermal_conduction_Noutput
-            thisOutput => thermal_conduction_output
-            thisSize   => thermal_conduction_sizePostResult
-          case default
-            valid = .false.
-        end select
-        if (valid) then
-          write(FILEUNIT,'(a)') '(thermal)'//char(9)//trim(outputName)
-          if (thermal_type(p) /= THERMAL_isothermal_ID) then
-            do e = 1,thisNoutput(i)
-              write(FILEUNIT,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i)
-            enddo
-          endif
-        endif
-        
-        i = damage_typeInstance(p)                                                                  ! which instance of this damage type
-        valid = .true.                                                                              ! assume valid
-        select case(damage_type(p))                                                                 ! split per damage type
-          case (DAMAGE_none_ID)
-            outputName = DAMAGE_none_label
-            thisNoutput => null()
-            thisOutput => null()
-            thisSize   => null()
-          case (DAMAGE_local_ID)
-            outputName = DAMAGE_local_label
-            thisNoutput => damage_local_Noutput
-            thisOutput => damage_local_output
-            thisSize   => damage_local_sizePostResult
-          case (DAMAGE_nonlocal_ID)
-            outputName = DAMAGE_nonlocal_label
-            thisNoutput => damage_nonlocal_Noutput
-            thisOutput => damage_nonlocal_output
-            thisSize   => damage_nonlocal_sizePostResult
-          case default
-            valid = .false.
-        end select
-        if (valid) then
-          write(FILEUNIT,'(a)') '(damage)'//char(9)//trim(outputName)
-          if (damage_type(p) /= DAMAGE_none_ID) then
-            do e = 1,thisNoutput(i)
-              write(FILEUNIT,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i)
-            enddo
-          endif
-        endif
-      endif
-    enddo
-    close(FILEUNIT)
-  endif mainProcess
-
   call config_deallocate('material.config/homogenization')
 
 !--------------------------------------------------------------------------------------------------
@@ -250,23 +160,7 @@ subroutine homogenization_init
   allocate(materialpoint_converged(discretization_nIP,discretization_nElem),          source=.true.)
   allocate(materialpoint_doneAndHappy(2,discretization_nIP,discretization_nElem),     source=.true.)
 
-!--------------------------------------------------------------------------------------------------
+  write(6,'(/,a)')   ' <<<+-  homogenization init  -+>>>'; flush(6)
-! allocate and initialize global state and postresutls variables
-  thermal_maxSizePostResults        = 0
-  damage_maxSizePostResults         = 0
-  do p = 1,size(config_homogenization)
-    thermal_maxSizePostResults      = max(thermal_maxSizePostResults, thermalState(p)%sizePostResults)
-    damage_maxSizePostResults       = max(damage_maxSizePostResults,  damageState (p)%sizePostResults)
-  enddo
-
-  materialpoint_sizeResults = 1 &                                                                   ! grain count
-                            + 1 + thermal_maxSizePostResults        &
-                                + damage_maxSizePostResults         &
-                            + homogenization_maxNgrains * (  1 &     ! crystallite size
-                                                           + 1 + constitutive_source_maxSizePostResults)
-  allocate(materialpoint_results(materialpoint_sizeResults,discretization_nIP,discretization_nElem))
-
-  write(6,'(/,a)')   ' <<<+-  homogenization init  -+>>>'
 
   if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0) then
     write(6,'(a32,1x,7(i8,1x))')   'materialpoint_dPdF:             ', shape(materialpoint_dPdF)
@@ -582,52 +476,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
 end subroutine materialpoint_stressAndItsTangent
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief parallelized calculation of result array at material points
-!--------------------------------------------------------------------------------------------------
-subroutine materialpoint_postResults
-
-  integer :: &
-    thePos, &
-    theSize, &
-    myNgrains, &
-    g, &                                                                                            !< grain number
-    i, &                                                                                            !< integration point number
-    e                                                                                               !< element number
-
-  !$OMP PARALLEL DO PRIVATE(myNgrains,thePos,theSize)
-  elementLooping: do e = FEsolving_execElem(1),FEsolving_execElem(2)
-    myNgrains = homogenization_Ngrains(material_homogenizationAt(e))
-    IpLooping: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
-      thePos = 0
-
-      theSize = thermalState     (material_homogenizationAt(e))%sizePostResults &
-              + damageState      (material_homogenizationAt(e))%sizePostResults
-      materialpoint_results(thePos+1,i,e) = real(theSize,pReal)                                    ! tell size of homogenization results
-      thePos = thePos + 1
-
-      if (theSize > 0) then                                                                         ! any homogenization results to mention?
-        materialpoint_results(thePos+1:thePos+theSize,i,e) = postResults(i,e)
-        thePos = thePos + theSize
-      endif
-
-      materialpoint_results(thePos+1,i,e) = real(myNgrains,pReal)                                   ! tell number of grains at materialpoint
-      thePos = thePos + 1
-
-      grainLooping :do g = 1,myNgrains
-        theSize = 1 + &
-                  1 + &
-                  sum(sourceState(material_phaseAt(g,e))%p(:)%sizePostResults)
-        materialpoint_results(thePos+1:thePos+theSize,i,e) = crystallite_postResults(g,i,e)        ! tell crystallite results
-        thePos = thePos + theSize
-      enddo grainLooping
-    enddo IpLooping
-  enddo elementLooping
- !$OMP END PARALLEL DO
-
-end subroutine materialpoint_postResults
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief  partition material point def grad onto constituents
 !--------------------------------------------------------------------------------------------------
@@ -739,53 +587,6 @@ subroutine averageStressAndItsTangent(ip,el)
 end subroutine averageStressAndItsTangent
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief return array of homogenization results for post file inclusion. call only,
-!> if homogenization_sizePostResults(i,e) > 0 !!
-!--------------------------------------------------------------------------------------------------
-function postResults(ip,el)
-
- integer, intent(in) :: &
-   ip, &                                                                                            !< integration point
-   el                                                                                               !< element number
- real(pReal), dimension(  thermalState     (material_homogenizationAt(el))%sizePostResults &
-                        + damageState      (material_homogenizationAt(el))%sizePostResults) :: &
-   postResults
- integer :: &
-   startPos, endPos ,&
-   homog
-
-
- postResults = 0.0_pReal
- startPos = 1
- endPos   = thermalState(material_homogenizationAt(el))%sizePostResults
- chosenThermal: select case (thermal_type(material_homogenizationAt(el)))
-
-   case (THERMAL_adiabatic_ID) chosenThermal
-     homog = material_homogenizationAt(el)
-     postResults(startPos:endPos) = &
-       thermal_adiabatic_postResults(homog,thermal_typeInstance(homog),thermalMapping(homog)%p(ip,el))
-   case (THERMAL_conduction_ID) chosenThermal
-     homog = material_homogenizationAt(el)
-     postResults(startPos:endPos) = &
-       thermal_conduction_postResults(homog,thermal_typeInstance(homog),thermalMapping(homog)%p(ip,el))
-
- end select chosenThermal
-
- startPos = endPos + 1
- endPos   = endPos + damageState(material_homogenizationAt(el))%sizePostResults
- chosenDamage: select case (damage_type(material_homogenizationAt(el)))
-
-   case (DAMAGE_local_ID) chosenDamage
-     postResults(startPos:endPos) = damage_local_postResults(ip, el)
-   case (DAMAGE_nonlocal_ID) chosenDamage
-     postResults(startPos:endPos) = damage_nonlocal_postResults(ip, el)
-     
- end select chosenDamage
-
-end function postResults
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief writes homogenization results to HDF5 output file
 !--------------------------------------------------------------------------------------------------
@@ -795,25 +596,16 @@ subroutine homogenization_results
     material_homogenization_type => homogenization_type
     
   integer :: p
-  character(len=256) :: group
+  character(len=pStringLen) :: group_base,group
   
   !real(pReal), dimension(:,:,:), allocatable :: temp
                                              
   do p=1,size(config_name_homogenization)
-    group = trim('current/materialpoint')//'/'//trim(config_name_homogenization(p))
+    group_base = 'current/materialpoint/'//trim(config_name_homogenization(p))
-    call HDF5_closeGroup(results_addGroup(group))
+    call results_closeGroup(results_addGroup(group_base))
-    
-    group = trim(group)//'/mech'
-    
-    call HDF5_closeGroup(results_addGroup(group))  
-    select case(material_homogenization_type(p))
-      case(HOMOGENIZATION_rgc_ID)
-        call mech_RGC_results(homogenization_typeInstance(p),group)
-    end select
-    
-    group = trim('current/materialpoint')//'/'//trim(config_name_homogenization(p))//'/generic'
-    call HDF5_closeGroup(results_addGroup(group))
     
+    group = trim(group_base)//'/generic'
+    call results_closeGroup(results_addGroup(group))
     !temp = reshape(materialpoint_F,[3,3,discretization_nIP*discretization_nElem])
     !call results_writeDataset(group,temp,'F',&
     !                          'deformation gradient','1')  
@@ -821,7 +613,34 @@ subroutine homogenization_results
     !call results_writeDataset(group,temp,'P',&
     !                          '1st Piola-Kirchoff stress','Pa')  
     
+    group = trim(group_base)//'/mech'
+    call results_closeGroup(results_addGroup(group))
+    select case(material_homogenization_type(p))
+      case(HOMOGENIZATION_rgc_ID)
+        call mech_RGC_results(homogenization_typeInstance(p),group)
+    end select
+
+    group = trim(group_base)//'/damage'
+    call results_closeGroup(results_addGroup(group))
+    select case(damage_type(p))
+      case(DAMAGE_LOCAL_ID)
+        call damage_local_results(p,group)
+      case(DAMAGE_NONLOCAL_ID)
+        call damage_nonlocal_results(p,group)
+    end select
+    
+    group = trim(group_base)//'/thermal'
+    call results_closeGroup(results_addGroup(group))
+    select case(thermal_type(p))
+      case(THERMAL_ADIABATIC_ID)
+        call thermal_adiabatic_results(p,group)
+      case(THERMAL_CONDUCTION_ID)
+        call thermal_conduction_results(p,group)
+    end select
+    
+    
  enddo
+
 #endif
 end subroutine homogenization_results
 

src/math.f90

@@ -838,33 +838,6 @@ pure function math_Voigt66to3333(m66)
 end function math_Voigt66to3333
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief action of a quaternion on a vector (rotate vector v with Q)
-!> @details deprecated
-!--------------------------------------------------------------------------------------------------
-pure function math_qRot(Q,v)
-
-  real(pReal), dimension(4), intent(in) :: Q
-  real(pReal), dimension(3), intent(in) :: v
-  real(pReal), dimension(3) :: math_qRot
-  real(pReal), dimension(4,4) :: T
-  integer :: i, j
-
-  do i = 1,4
-    do j = 1,i
-      T(i,j) = Q(i) * Q(j)
-    enddo
-  enddo
-
-  math_qRot = [-v(1)*(T(3,3)+T(4,4)) + v(2)*(T(3,2)-T(4,1)) + v(3)*(T(4,2)+T(3,1)), &
-                v(1)*(T(3,2)+T(4,1)) - v(2)*(T(2,2)+T(4,4)) + v(3)*(T(4,3)-T(2,1)), &
-                v(1)*(T(4,2)-T(3,1)) + v(2)*(T(4,3)+T(2,1)) - v(3)*(T(2,2)+T(3,3))]
-
-  math_qRot = 2.0_pReal * math_qRot + v
-
-end function math_qRot
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief rotation matrix from Bunge-Euler (3-1-3) angles (in radians)
 !> @details deprecated
@@ -1328,55 +1301,6 @@ real(pReal) pure function math_areaTriangle(v1,v2,v3)
 end function math_areaTriangle
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief rotate 33 tensor forward
-!> @details deprecated
-!--------------------------------------------------------------------------------------------------
-pure function math_rotate_forward33(tensor,R)
-
-  real(pReal), dimension(3,3)             ::  math_rotate_forward33
-  real(pReal), dimension(3,3), intent(in) :: tensor, R
-
-  math_rotate_forward33 = matmul(R,matmul(tensor,transpose(R)))
-
-end function math_rotate_forward33
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief rotate 33 tensor backward
-!> @details deprecated
-!--------------------------------------------------------------------------------------------------
-pure function math_rotate_backward33(tensor,R)
-
-  real(pReal), dimension(3,3)             ::  math_rotate_backward33
-  real(pReal), dimension(3,3), intent(in) :: tensor, R
-
-  math_rotate_backward33 = matmul(transpose(R),matmul(tensor,R))
-
-end function math_rotate_backward33
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief rotate 3333 tensor C'_ijkl=g_im*g_jn*g_ko*g_lp*C_mnop
-!> @details deprecated
-!--------------------------------------------------------------------------------------------------
-pure function math_rotate_forward3333(tensor,R)
-
-  real(pReal), dimension(3,3,3,3)             ::  math_rotate_forward3333
-  real(pReal), dimension(3,3),     intent(in) :: R
-  real(pReal), dimension(3,3,3,3), intent(in) :: tensor
-  integer :: i,j,k,l,m,n,o,p
-
-  math_rotate_forward3333 = 0.0_pReal
-  do i = 1,3;do j = 1,3;do k = 1,3;do l = 1,3
-  do m = 1,3;do n = 1,3;do o = 1,3;do p = 1,3
-    math_rotate_forward3333(i,j,k,l) = math_rotate_forward3333(i,j,k,l) &
-                                     + R(i,m) * R(j,n) * R(k,o) * R(l,p) * tensor(m,n,o,p)
-  enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo
-
-end function math_rotate_forward3333
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief limits a scalar value to a certain range (either one or two sided)
 ! Will return NaN if left > right

src/mesh/DAMASK_FEM.f90

@@ -73,7 +73,7 @@ program DAMASK_FEM
 !--------------------------------------------------------------------------------------------------
 ! init DAMASK (all modules)
   call CPFEM_initAll
-  write(6,'(/,a)')   ' <<<+-  DAMASK_FEM init  -+>>>'
+  write(6,'(/,a)')   ' <<<+-  DAMASK_FEM init  -+>>>'; flush(6)
  
 ! reading basic information from load case file and allocate data structure containing load cases
   call DMGetDimension(geomMesh,dimPlex,ierr); CHKERRA(ierr)                                         !< dimension of mesh (2D or 3D)
@@ -296,19 +296,12 @@ program DAMASK_FEM
 !--------------------------------------------------------------------------------------------------
 ! report begin of new step
         write(6,'(/,a)') ' ###########################################################################'
-        write(6,'(1x,a,es12.5'//&
+        write(6,'(1x,a,es12.5,6(a,i0))')&
-                ',a,'//IO_intOut(inc)//',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//&
-                ',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//&
-                ',a,'//IO_intOut(currentLoadCase)//',a,'//IO_intOut(size(loadCases))//')') &
                 'Time', time, &
                 's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,&
                 '-', stepFraction, '/', subStepFactor**cutBackLevel,&
                 ' of load case ', currentLoadCase,'/',size(loadCases)
-        write(incInfo,&
+        write(incInfo,'(4(a,i0))') &
-               '(a,'//IO_intOut(totalIncsCounter)//&
-               ',a,'//IO_intOut(sum(loadCases%incs))//&
-               ',a,'//IO_intOut(stepFraction)//&
-               ',a,'//IO_intOut(subStepFactor**cutBackLevel)//')') &
                'Increment ',totalIncsCounter,'/',sum(loadCases%incs),&
                '-',stepFraction, '/', subStepFactor**cutBackLevel
         flush(6)
@@ -373,11 +366,9 @@ program DAMASK_FEM
       cutBackLevel = max(0, cutBackLevel - 1)                                                       ! try half number of subincs next inc
 
       if (all(solres(:)%converged)) then
-        write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') &                                      ! report converged inc
+        write(6,'(/,a,i0,a)') ' increment ', totalIncsCounter, ' converged'
-                                  ' increment ', totalIncsCounter, ' converged'
       else
-        write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') &                                      ! report non-converged inc
+        write(6,'(/,a,i0,a)') ' increment ', totalIncsCounter, ' NOT converged'
-                                  ' increment ', totalIncsCounter, ' NOT converged'
       endif; flush(6)
 
       if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then                            ! at output frequency

src/plastic_nonlocal.f90

@@ -24,10 +24,10 @@ module plastic_nonlocal
     
   implicit none
   private
-  real(pReal), parameter, private :: &
+  real(pReal), parameter :: &
     KB = 1.38e-23_pReal                                                                             !< Physical parameter, Boltzmann constant in J/Kelvin
   
-  character(len=64), dimension(:,:), allocatable, target, public :: &
+  character(len=64), dimension(:,:), allocatable :: &
     plastic_nonlocal_output                                                                         !< name of each post result output
 
   ! storage order of dislocation types
@@ -50,17 +50,17 @@ module plastic_nonlocal
     mob_scr_neg = 4                                                                                 !< mobile screw positive
 
   ! BEGIN DEPRECATES
-  integer, dimension(:,:,:), allocatable, private :: &
+  integer, dimension(:,:,:), allocatable :: &
     iRhoU, &                                                                                        !< state indices for unblocked density
     iRhoB, &                                                                                        !< state indices for blocked density
     iRhoD, &                                                                                        !< state indices for dipole density
     iV, &                                                                                           !< state indices for dislcation velocities
     iD                                                                                              !< state indices for stable dipole height
-  integer, dimension(:), allocatable, private, protected :: &
+  integer, dimension(:), allocatable :: &
     totalNslip                                                                                      !< total number of active slip systems for each instance
   !END DEPRECATED
   
-  real(pReal), dimension(:,:,:,:,:,:), allocatable, private :: &
+  real(pReal), dimension(:,:,:,:,:,:), allocatable :: &
     compatibility                                                                                   !< slip system compatibility between me and my neighbors
   
   enum, bind(c) 
@@ -89,7 +89,7 @@ module plastic_nonlocal
       gamma_ID
   end enum
   
-  type, private :: tParameters                                                                      !< container type for internal constitutive parameters
+  type :: tParameters                                                                               !< container type for internal constitutive parameters
     real(pReal) :: &
       atomicVolume, &                                                                               !< atomic volume
       Dsd0, &                                                                                       !< prefactor for self-diffusion coefficient
@@ -139,19 +139,19 @@ module plastic_nonlocal
       interactionSlipSlip ,&                                                                        !< coefficients for slip-slip interaction
       forestProjection_Edge, &                                                                      !< matrix of forest projections of edge dislocations
       forestProjection_Screw                                                                        !< matrix of forest projections of screw dislocations
-    real(pReal), dimension(:), allocatable, private :: &
+    real(pReal), dimension(:), allocatable :: &
       nonSchmidCoeff
-    real(pReal), dimension(:,:,:), allocatable, private :: &
+    real(pReal), dimension(:,:,:), allocatable :: &
       Schmid, &                                                                                     !< Schmid contribution
       nonSchmid_pos, &
       nonSchmid_neg                                                                                 !< combined projection of Schmid and non-Schmid contributions to the resolved shear stress (only for screws)
     integer :: &
       totalNslip
-    integer, dimension(:) ,allocatable , public:: &
+    integer, dimension(:) ,allocatable :: &
       Nslip,&
       colinearSystem                                                                                !< colinear system to the active slip system (only valid for fcc!)
  
-    logical, private :: &
+    logical :: &
       shortRangeStressCorrection, &                                                                 !< flag indicating the use of the short range stress correction by a excess density gradient term
       probabilisticMultiplication
     
@@ -160,13 +160,13 @@ module plastic_nonlocal
       
   end type tParameters
    
-  type, private :: tNonlocalMicrostructure
+  type :: tNonlocalMicrostructure
     real(pReal), allocatable, dimension(:,:) :: &
      tau_pass, & 
      tau_Back 
   end type tNonlocalMicrostructure
   
-  type, private :: tNonlocalState
+  type :: tNonlocalState
     real(pReal), pointer, dimension(:,:) :: &
       rho, &                                                                                        ! < all dislocations
         rhoSgl, &
@@ -192,16 +192,16 @@ module plastic_nonlocal
           v_scr_neg
   end type tNonlocalState
   
-  type(tNonlocalState), allocatable, dimension(:), private :: &
+  type(tNonlocalState), allocatable, dimension(:) :: &
     deltaState, &
     dotState, &
     state
  
-  type(tParameters), dimension(:), allocatable, private :: param                                    !< containers of constitutive parameters (len Ninstance)
+  type(tParameters), dimension(:), allocatable :: param                                             !< containers of constitutive parameters (len Ninstance)
  
-  type(tNonlocalMicrostructure), dimension(:), allocatable, private :: microstructure
+  type(tNonlocalMicrostructure), dimension(:), allocatable :: microstructure
 
-  integer(kind(undefined_ID)), dimension(:,:),   allocatable, private :: & 
+  integer(kind(undefined_ID)), dimension(:,:), allocatable :: & 
     plastic_nonlocal_outputID                                                                       !< ID of each post result output
   
   public :: &
@@ -1829,8 +1829,6 @@ subroutine plastic_nonlocal_updateCompatibility(orientation,i,e)
     ns, &                                                                                           ! number of active slip systems
     s1, &                                                                                           ! slip system index (me)
     s2                                                                                              ! slip system index (my neighbor)
-  real(pReal), dimension(4) :: &
-    absoluteMisorientation                                                                          ! absolute misorientation (without symmetry) between me and my neighbor
   real(pReal), dimension(2,totalNslip(phase_plasticityInstance(material_phaseAt(1,e))),&
                            totalNslip(phase_plasticityInstance(material_phaseAt(1,e))),&
                            nIPneighbors) :: &  
@@ -1841,7 +1839,7 @@ subroutine plastic_nonlocal_updateCompatibility(orientation,i,e)
     nThresholdValues
   logical, dimension(totalNslip(phase_plasticityInstance(material_phaseAt(1,e)))) :: & 
     belowThreshold
-  type(rotation) :: rot
+  type(rotation) :: mis
 
   Nneighbors = nIPneighbors
   ph = material_phaseAt(1,e)
@@ -1907,18 +1905,17 @@ subroutine plastic_nonlocal_updateCompatibility(orientation,i,e)
     !* Finally the smallest compatibility value is decreased until the sum is exactly equal to one. 
     !* All values below the threshold are set to zero.
     else
-      rot = orientation(1,i,e)%misorientation(orientation(1,neighbor_i,neighbor_e))
+      mis = orientation(1,i,e)%misorientation(orientation(1,neighbor_i,neighbor_e))
-      absoluteMisorientation = rot%asQuaternion()
       mySlipSystems: do s1 = 1,ns
         neighborSlipSystems: do s2 = 1,ns
           my_compatibility(1,s2,s1,n) =  math_inner(prm%slip_normal(1:3,s1), &
-          math_qRot(absoluteMisorientation, prm%slip_normal(1:3,s2))) &
+                                                       mis%rotate(prm%slip_normal(1:3,s2))) &
                                   * abs(math_inner(prm%slip_direction(1:3,s1), &
-                                  math_qRot(absoluteMisorientation, prm%slip_direction(1:3,s2))))
+                                                       mis%rotate(prm%slip_direction(1:3,s2))))
           my_compatibility(2,s2,s1,n) = abs(math_inner(prm%slip_normal(1:3,s1), &
-          math_qRot(absoluteMisorientation, prm%slip_normal(1:3,s2)))) &
+                                                       mis%rotate(prm%slip_normal(1:3,s2)))) &
                                   * abs(math_inner(prm%slip_direction(1:3,s1), &
-                                  math_qRot(absoluteMisorientation, prm%slip_direction(1:3,s2))))
+                                                       mis%rotate(prm%slip_direction(1:3,s2))))
         enddo neighborSlipSystems
         
         my_compatibilitySum = 0.0_pReal

src/prec.f90

@@ -42,8 +42,7 @@ module prec
       sizeState        = 0, &                                                                       !< size of state
       sizeDotState     = 0, &                                                                       !< size of dot state, i.e. state(1:sizeDot) follows time evolution by dotState rates
       offsetDeltaState = 0, &                                                                       !< index offset of delta state
-      sizeDeltaState   = 0, &                                                                       !< size of delta state, i.e. state(offset+1:offset+sizeDelta) follows time evolution by deltaState increments
+      sizeDeltaState   = 0                                                                          !< size of delta state, i.e. state(offset+1:offset+sizeDelta) follows time evolution by deltaState increments
-      sizePostResults  = 0                                                                          !< size of output data
     real(pReal), pointer,     dimension(:), contiguous :: &
       atolState
     real(pReal), pointer,     dimension(:,:), contiguous :: &                                       ! a pointer is needed here because we might point to state/doState. However, they will never point to something, but are rather allocated and, hence, contiguous

src/results.f90

@@ -69,10 +69,9 @@ subroutine results_init
   write(6,'(a)')   ' https://doi.org/10.1007/s40192-017-0084-5'
 
   resultsFile = HDF5_openFile(trim(getSolverJobName())//'.hdf5','w',.true.)
-  call HDF5_addAttribute(resultsFile,'DADF5-version',0.3_pReal)
+  call HDF5_addAttribute(resultsFile,'DADF5_version_major',0)
-  call HDF5_addAttribute(resultsFile,'DADF5-major',0)
+  call HDF5_addAttribute(resultsFile,'DADF5_version_minor',4)
-  call HDF5_addAttribute(resultsFile,'DADF5-minor',3)
+  call HDF5_addAttribute(resultsFile,'DAMASK_version',DAMASKVERSION)
-  call HDF5_addAttribute(resultsFile,'DAMASK',DAMASKVERSION)
   call get_command(commandLine)
   call HDF5_addAttribute(resultsFile,'call',trim(commandLine))
   call HDF5_closeGroup(results_addGroup('mapping'))
@@ -111,7 +110,7 @@ subroutine results_addIncrement(inc,time)
   real(pReal),   intent(in) :: time
   character(len=pStringLen) :: incChar
 
-  write(incChar,'(i5.5)') inc                                                                       ! allow up to 99999 increments
+  write(incChar,'(i10)') inc
   call HDF5_closeGroup(results_addGroup(trim('inc'//trim(adjustl(incChar)))))
   call results_setLink(trim('inc'//trim(adjustl(incChar))),'current')
   call HDF5_addAttribute(resultsFile,'time/s',time,trim('inc'//trim(adjustl(incChar))))
@@ -455,7 +454,7 @@ subroutine results_mapping_constituent(phaseAt,memberAt,label)
     
   integer,                   dimension(:,:),   intent(in) :: phaseAt                                !< phase section at (constituent,element)
   integer,                   dimension(:,:,:), intent(in) :: memberAt                               !< phase member at (constituent,IP,element)
-  character(len=64), dimension(:),     intent(in)  :: label                                         !< label of each phase section
+  character(len=pStringLen), dimension(:),     intent(in) :: label                                  !< label of each phase section
   
   integer, dimension(size(memberAt,1),size(memberAt,2),size(memberAt,3)) :: &
     phaseAt_perIP, &
@@ -590,7 +589,7 @@ subroutine results_mapping_materialpoint(homogenizationAt,memberAt,label)
     
   integer,                   dimension(:),   intent(in) :: homogenizationAt                         !< homogenization section at (element)
   integer,                   dimension(:,:), intent(in) :: memberAt                                 !< homogenization member at (IP,element)
-  character(len=64), dimension(:),   intent(in)  :: label                                           !< label of each homogenization section
+  character(len=pStringLen), dimension(:),   intent(in) :: label                                    !< label of each homogenization section
   
   integer, dimension(size(memberAt,1),size(memberAt,2)) :: &
     homogenizationAt_perIP, &

src/rotations.f90

@@ -64,6 +64,7 @@ module rotations
       procedure, public :: asRodrigues
       procedure, public :: asMatrix
       !------------------------------------------
+      procedure, public :: fromQuaternion
       procedure, public :: fromEulers
       procedure, public :: fromAxisAngle
       procedure, public :: fromMatrix
@@ -157,6 +158,18 @@ end function asHomochoric
 !---------------------------------------------------------------------------------------------------
 ! Initialize rotation from different representations
 !---------------------------------------------------------------------------------------------------
+subroutine fromQuaternion(self,qu)
+
+  class(rotation), intent(out)          :: self
+  real(pReal), dimension(4), intent(in) :: qu
+
+  if (dNeq(norm2(qu),1.0_pReal)) &
+    call IO_error(402,ext_msg='fromQuaternion')
+
+  self%q = qu
+
+end subroutine fromQuaternion
+!---------------------------------------------------------------------------------------------------
 subroutine fromEulers(self,eu,degrees)
 
   class(rotation), intent(out)          :: self

src/source_damage_anisoBrittle.f90

@@ -13,6 +13,7 @@ module source_damage_anisoBrittle
   use discretization
   use config
   use lattice
+  use results
 
   implicit none
   private
@@ -21,10 +22,7 @@ module source_damage_anisoBrittle
     source_damage_anisoBrittle_offset, &                                                            !< which source is my current source mechanism?
     source_damage_anisoBrittle_instance                                                             !< instance of source mechanism
 
-  integer,                       dimension(:,:),         allocatable, target, public  :: &
+  character(len=64),             dimension(:,:),         allocatable :: &
-    source_damage_anisoBrittle_sizePostResult                                                       !< size of each post result output
-
-  character(len=64),             dimension(:,:),         allocatable, target, public  :: &
     source_damage_anisoBrittle_output                                                               !< name of each post result output
     
   integer,                       dimension(:,:),         allocatable :: &
@@ -61,7 +59,7 @@ module source_damage_anisoBrittle
     source_damage_anisoBrittle_init, &
     source_damage_anisoBrittle_dotState, &
     source_damage_anisobrittle_getRateAndItsTangent, &
-    source_damage_anisoBrittle_postResults
+    source_damage_anisoBrittle_results
 
 contains
 
@@ -102,7 +100,6 @@ subroutine source_damage_anisoBrittle_init
     enddo    
   enddo
   
-  allocate(source_damage_anisoBrittle_sizePostResult(maxval(phase_Noutput),Ninstance), source=0)
   allocate(source_damage_anisoBrittle_output(maxval(phase_Noutput),Ninstance))
            source_damage_anisoBrittle_output = ''
 
@@ -154,7 +151,6 @@ subroutine source_damage_anisoBrittle_init
       select case(outputs(i))
       
         case ('anisobrittle_drivingforce')
-          source_damage_anisoBrittle_sizePostResult(i,source_damage_anisoBrittle_instance(p)) = 1
           source_damage_anisoBrittle_output(i,source_damage_anisoBrittle_instance(p)) = outputs(i)
           prm%outputID = [prm%outputID, damage_drivingforce_ID]
 
@@ -171,7 +167,6 @@ subroutine source_damage_anisoBrittle_init
 
 
     call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0)
-    sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_anisoBrittle_sizePostResult(:,instance))
     sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol
 
 
@@ -262,39 +257,32 @@ subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalph
   
   dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent)
  
-end subroutine source_damage_anisobrittle_getRateAndItsTangent
+end subroutine source_damage_anisoBrittle_getRateAndItsTangent
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief return array of local damage results
+!> @brief writes results to HDF5 output file
 !--------------------------------------------------------------------------------------------------
-function source_damage_anisoBrittle_postResults(phase, constituent)
+subroutine source_damage_anisoBrittle_results(phase,group)
 
-  integer, intent(in) :: &
+  integer, intent(in) :: phase
-    phase, &
+  character(len=*), intent(in) :: group
-    constituent
+#if defined(PETSc) || defined(DAMASK_HDF5)  
-
+  integer :: sourceOffset, o, instance
-  real(pReal), dimension(sum(source_damage_anisoBrittle_sizePostResult(:, &
-                           source_damage_anisoBrittle_instance(phase)))) :: &
-    source_damage_anisoBrittle_postResults
-
-  integer :: &
-    instance, sourceOffset, o, c
    
   instance     = source_damage_anisoBrittle_instance(phase)
   sourceOffset = source_damage_anisoBrittle_offset(phase)
 
-  c = 0
+   associate(prm => param(instance), stt => sourceState(phase)%p(sourceOffset)%state)
-
+   outputsLoop: do o = 1,size(prm%outputID)
-  do o = 1,size(param(instance)%outputID)
+     select case(prm%outputID(o))
-     select case(param(instance)%outputID(o))
        case (damage_drivingforce_ID)
-         source_damage_anisoBrittle_postResults(c+1) = &
+         call results_writeDataset(group,stt,'tbd','driving force','tbd')
-           sourceState(phase)%p(sourceOffset)%state(1,constituent)
-         c = c + 1
-
      end select
-  enddo
+   enddo outputsLoop
-end function source_damage_anisoBrittle_postResults
+   end associate
+#endif
+
+end subroutine source_damage_anisoBrittle_results
 
 end module source_damage_anisoBrittle

src/source_damage_anisoDuctile.f90

@@ -12,6 +12,7 @@ module source_damage_anisoDuctile
   use discretization
   use material
   use config
+  use results
  
   implicit none
   private
@@ -20,9 +21,6 @@ module source_damage_anisoDuctile
     source_damage_anisoDuctile_offset, &                                                            !< which source is my current damage mechanism?
     source_damage_anisoDuctile_instance                                                             !< instance of damage source mechanism
  
-  integer,                       dimension(:,:),         allocatable, target, public  :: &
-    source_damage_anisoDuctile_sizePostResult                                                       !< size of each post result output
- 
   character(len=64),             dimension(:,:),         allocatable, target, public  :: &
     source_damage_anisoDuctile_output                                                               !< name of each post result output
     
@@ -54,7 +52,7 @@ module source_damage_anisoDuctile
     source_damage_anisoDuctile_init, &
     source_damage_anisoDuctile_dotState, &
     source_damage_anisoDuctile_getRateAndItsTangent, &
-    source_damage_anisoDuctile_postResults
+    source_damage_anisoDuctile_results
 
 contains
 
@@ -96,7 +94,6 @@ subroutine source_damage_anisoDuctile_init
     enddo    
   enddo
     
-  allocate(source_damage_anisoDuctile_sizePostResult(maxval(phase_Noutput),Ninstance),source=0)
   allocate(source_damage_anisoDuctile_output(maxval(phase_Noutput),Ninstance))
            source_damage_anisoDuctile_output = ''
  
@@ -139,7 +136,6 @@ subroutine source_damage_anisoDuctile_init
       select case(outputs(i))
       
         case ('anisoductile_drivingforce')
-          source_damage_anisoDuctile_sizePostResult(i,source_damage_anisoDuctile_instance(p)) = 1
           source_damage_anisoDuctile_output(i,source_damage_anisoDuctile_instance(p)) = outputs(i)
           prm%outputID = [prm%outputID, damage_drivingforce_ID]
  
@@ -156,7 +152,6 @@ subroutine source_damage_anisoDuctile_init
     sourceOffset = source_damage_anisoDuctile_offset(phase)
  
     call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0)
-    sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_anisoDuctile_sizePostResult(:,instance))
       sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol
     
   enddo
@@ -226,35 +221,28 @@ end subroutine source_damage_anisoDuctile_getRateAndItsTangent
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief return array of local damage results
+!> @brief writes results to HDF5 output file
 !--------------------------------------------------------------------------------------------------
-function source_damage_anisoDuctile_postResults(phase, constituent)
+subroutine source_damage_anisoDuctile_results(phase,group)
 
-  integer, intent(in) :: &
+  integer, intent(in) :: phase
-    phase, &
+  character(len=*), intent(in) :: group
-    constituent
+#if defined(PETSc) || defined(DAMASK_HDF5)  
-  real(pReal), dimension(sum(source_damage_anisoDuctile_sizePostResult(:, &
+  integer :: sourceOffset, o, instance
-                           source_damage_anisoDuctile_instance(phase)))) :: &
-    source_damage_anisoDuctile_postResults
- 
-  integer :: &
-    instance, sourceOffset, o, c
    
   instance     = source_damage_anisoDuctile_instance(phase)
   sourceOffset = source_damage_anisoDuctile_offset(phase)
 
-  c = 0
+   associate(prm => param(instance), stt => sourceState(phase)%p(sourceOffset)%state)
- 
+   outputsLoop: do o = 1,size(prm%outputID)
-  do o = 1,size(param(instance)%outputID)
+     select case(prm%outputID(o))
-     select case(param(instance)%outputID(o))
        case (damage_drivingforce_ID)
-         source_damage_anisoDuctile_postResults(c+1) = &
+         call results_writeDataset(group,stt,'tbd','driving force','tbd')
-           sourceState(phase)%p(sourceOffset)%state(1,constituent)
-         c = c + 1
- 
      end select
-  enddo
+   enddo outputsLoop
+   end associate
+#endif
 
-end function source_damage_anisoDuctile_postResults
+end subroutine source_damage_anisoDuctile_results
 
 end module source_damage_anisoDuctile

src/source_damage_isoBrittle.f90

@@ -12,15 +12,14 @@ module source_damage_isoBrittle
   use discretization
   use material
   use config
+  use results
 
   implicit none
   private
   integer,                       dimension(:),           allocatable,         public, protected :: &
     source_damage_isoBrittle_offset, &
     source_damage_isoBrittle_instance
-  integer,                       dimension(:,:),         allocatable, target, public :: &
+  character(len=64),             dimension(:,:),         allocatable :: &
-    source_damage_isoBrittle_sizePostResult
-  character(len=64),             dimension(:,:),         allocatable, target, public :: &
     source_damage_isoBrittle_output
 
   enum, bind(c) 
@@ -46,7 +45,7 @@ module source_damage_isoBrittle
     source_damage_isoBrittle_init, &
     source_damage_isoBrittle_deltaState, &
     source_damage_isoBrittle_getRateAndItsTangent, &
-    source_damage_isoBrittle_postResults
+    source_damage_isoBrittle_Results
 
 contains
 
@@ -86,7 +85,6 @@ subroutine source_damage_isoBrittle_init
     enddo    
   enddo
 
-  allocate(source_damage_isoBrittle_sizePostResult(maxval(phase_Noutput),Ninstance),source=0)
   allocate(source_damage_isoBrittle_output(maxval(phase_Noutput),Ninstance))
            source_damage_isoBrittle_output = ''
  
@@ -122,7 +120,6 @@ subroutine source_damage_isoBrittle_init
       select case(outputs(i))
       
         case ('isobrittle_drivingforce')
-          source_damage_isoBrittle_sizePostResult(i,source_damage_isoBrittle_instance(p)) = 1
           source_damage_isoBrittle_output(i,source_damage_isoBrittle_instance(p)) = outputs(i)
           prm%outputID = [prm%outputID, damage_drivingforce_ID]
             
@@ -139,7 +136,6 @@ subroutine source_damage_isoBrittle_init
    sourceOffset = source_damage_isoBrittle_offset(phase)
       
    call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,1)
-   sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_isoBrittle_sizePostResult(:,instance))
    sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol
   
   enddo
@@ -215,34 +211,30 @@ subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiD
   
 end subroutine source_damage_isoBrittle_getRateAndItsTangent
 
-!--------------------------------------------------------------------------------------------------
-!> @brief return array of local damage results
-!--------------------------------------------------------------------------------------------------
-function source_damage_isoBrittle_postResults(phase, constituent)
 
-  integer, intent(in) :: &
+!--------------------------------------------------------------------------------------------------
-    phase, &
+!> @brief writes results to HDF5 output file
-    constituent
+!--------------------------------------------------------------------------------------------------
-  real(pReal), dimension(sum(source_damage_isoBrittle_sizePostResult(:, &
+subroutine source_damage_isoBrittle_results(phase,group)
-                          source_damage_isoBrittle_instance(phase)))) :: &
-    source_damage_isoBrittle_postResults
 
-  integer :: &
+  integer, intent(in) :: phase
-    instance, sourceOffset, o, c
+  character(len=*), intent(in) :: group
+#if defined(PETSc) || defined(DAMASK_HDF5)  
+  integer :: sourceOffset, o, instance
    
   instance     = source_damage_isoBrittle_instance(phase)
   sourceOffset = source_damage_isoBrittle_offset(phase)
 
-  c = 0
+   associate(prm => param(instance), stt => sourceState(phase)%p(sourceOffset)%state)
-
+   outputsLoop: do o = 1,size(prm%outputID)
-  do o = 1,size(param(instance)%outputID)
+     select case(prm%outputID(o))
-     select case(param(instance)%outputID(o))
        case (damage_drivingforce_ID)
-         source_damage_isoBrittle_postResults(c+1) = sourceState(phase)%p(sourceOffset)%state(1,constituent)
+         call results_writeDataset(group,stt,'tbd','driving force','tbd')
-         c = c + 1
-
      end select
-  enddo
+   enddo outputsLoop
-end function source_damage_isoBrittle_postResults
+   end associate
+#endif
+
+end subroutine source_damage_isoBrittle_results
 
 end module source_damage_isoBrittle

src/source_damage_isoDuctile.f90

@@ -11,6 +11,7 @@ module source_damage_isoDuctile
  use discretization
  use material
  use config
+ use results
 
  implicit none
  private
@@ -18,9 +19,6 @@ module source_damage_isoDuctile
    source_damage_isoDuctile_offset, &                                                                 !< which source is my current damage mechanism?
    source_damage_isoDuctile_instance                                                                  !< instance of damage source mechanism
 
- integer,                       dimension(:,:),         allocatable, target, public :: &
-   source_damage_isoDuctile_sizePostResult                                                            !< size of each post result output
-
  character(len=64),             dimension(:,:),         allocatable, target, public :: &
    source_damage_isoDuctile_output                                                                    !< name of each post result output
 
@@ -46,7 +44,7 @@ module source_damage_isoDuctile
    source_damage_isoDuctile_init, &
    source_damage_isoDuctile_dotState, &
    source_damage_isoDuctile_getRateAndItsTangent, &
-   source_damage_isoDuctile_postResults
+   source_damage_isoDuctile_Results
 
 contains
 
@@ -86,7 +84,6 @@ subroutine source_damage_isoDuctile_init
    enddo    
  enddo
    
- allocate(source_damage_isoDuctile_sizePostResult(maxval(phase_Noutput),Ninstance),source=0)
  allocate(source_damage_isoDuctile_output(maxval(phase_Noutput),Ninstance))
           source_damage_isoDuctile_output = ''
 
@@ -122,7 +119,6 @@ subroutine source_damage_isoDuctile_init
      select case(outputs(i))
      
        case ('isoductile_drivingforce')
-         source_damage_isoDuctile_sizePostResult(i,source_damage_isoDuctile_instance(p)) = 1
          source_damage_isoDuctile_output(i,source_damage_isoDuctile_instance(p)) = outputs(i)
          prm%outputID = [prm%outputID, damage_drivingforce_ID]
 
@@ -138,10 +134,8 @@ subroutine source_damage_isoDuctile_init
    sourceOffset = source_damage_isoDuctile_offset(phase)
 
    call material_allocateSourceState(phase,sourceOffset,NofMyPhase,1,1,0)
-   sourceState(phase)%p(sourceOffset)%sizePostResults = sum(source_damage_isoDuctile_sizePostResult(:,instance))
    sourceState(phase)%p(sourceOffset)%aTolState=param(instance)%aTol
  
- 
  enddo
  
 end subroutine source_damage_isoDuctile_init
@@ -197,34 +191,31 @@ subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiD
  
 end subroutine source_damage_isoDuctile_getRateAndItsTangent
 
-!--------------------------------------------------------------------------------------------------
-!> @brief return array of local damage results
-!--------------------------------------------------------------------------------------------------
-function source_damage_isoDuctile_postResults(phase, constituent)
 
- integer, intent(in) :: &
+!--------------------------------------------------------------------------------------------------
-   phase, &
+!> @brief writes results to HDF5 output file
-   constituent
+!--------------------------------------------------------------------------------------------------
- real(pReal), dimension(sum(source_damage_isoDuctile_sizePostResult(:, &
+subroutine source_damage_isoDuctile_results(phase,group)
-                          source_damage_isoDuctile_instance(phase)))) :: &
-   source_damage_isoDuctile_postResults
 
- integer :: &
+  integer, intent(in) :: phase
-   instance, sourceOffset, o, c
+  character(len=*), intent(in) :: group
+#if defined(PETSc) || defined(DAMASK_HDF5)  
+  integer :: sourceOffset, o, instance
    
   instance     = source_damage_isoDuctile_instance(phase)
   sourceOffset = source_damage_isoDuctile_offset(phase)
 
- c = 0
+   associate(prm => param(instance), stt => sourceState(phase)%p(sourceOffset)%state)
-
+   outputsLoop: do o = 1,size(prm%outputID)
- do o = 1,size(param(instance)%outputID)
+     select case(prm%outputID(o))
-    select case(param(instance)%outputID(o))
        case (damage_drivingforce_ID)
-        source_damage_isoDuctile_postResults(c+1) = sourceState(phase)%p(sourceOffset)%state(1,constituent)
+         call results_writeDataset(group,stt,'tbd','driving force','tbd')
-        c = c + 1
-
      end select
- enddo
+   enddo outputsLoop
-end function source_damage_isoDuctile_postResults
+   end associate
+#endif
+
+end subroutine source_damage_isoDuctile_results
+
 
 end module source_damage_isoDuctile

src/thermal_adiabatic.f90

@@ -7,6 +7,7 @@ module thermal_adiabatic
   use config
   use numerics
   use material
+  use results
   use source_thermal_dissipation
   use source_thermal_externalheat
   use crystallite
@@ -15,8 +16,6 @@ module thermal_adiabatic
   implicit none
   private
  
-  integer,                       dimension(:,:),  allocatable, target, public :: &
-    thermal_adiabatic_sizePostResult                                                                !< size of each post result output
   character(len=64),             dimension(:,:),  allocatable, target, public :: &
     thermal_adiabatic_output                                                                        !< name of each post result output
     
@@ -37,7 +36,7 @@ module thermal_adiabatic
     thermal_adiabatic_getSourceAndItsTangent, &
     thermal_adiabatic_getSpecificHeat, &
     thermal_adiabatic_getMassDensity, &
-    thermal_adiabatic_postResults
+    thermal_adiabatic_results
  
 contains
 
@@ -57,7 +56,6 @@ subroutine thermal_adiabatic_init
   maxNinstance = count(thermal_type == THERMAL_adiabatic_ID)
   if (maxNinstance == 0) return
   
-  allocate(thermal_adiabatic_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0)
   allocate(thermal_adiabatic_output         (maxval(homogenization_Noutput),maxNinstance))
            thermal_adiabatic_output = ''
   allocate(thermal_adiabatic_outputID       (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID)
@@ -75,14 +73,12 @@ subroutine thermal_adiabatic_init
               thermal_adiabatic_Noutput(instance) = thermal_adiabatic_Noutput(instance) + 1
               thermal_adiabatic_outputID(thermal_adiabatic_Noutput(instance),instance) = temperature_ID
               thermal_adiabatic_output(thermal_adiabatic_Noutput(instance),instance) = outputs(i)
-              thermal_adiabatic_sizePostResult(thermal_adiabatic_Noutput(instance),instance) = 1
       end select
     enddo
  
  ! allocate state arrays
     sizeState = 1
     thermalState(section)%sizeState = sizeState
-    thermalState(section)%sizePostResults = sum(thermal_adiabatic_sizePostResult(:,instance))
     allocate(thermalState(section)%state0   (sizeState,NofMyHomog), source=thermal_initialT(section))
     allocate(thermalState(section)%subState0(sizeState,NofMyHomog), source=thermal_initialT(section))
     allocate(thermalState(section)%state    (sizeState,NofMyHomog), source=thermal_initialT(section))
@@ -252,32 +248,27 @@ end function thermal_adiabatic_getMassDensity
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief return array of thermal results
+!> @brief writes results to HDF5 output file
 !--------------------------------------------------------------------------------------------------
-function thermal_adiabatic_postResults(homog,instance,of) result(postResults)
+subroutine thermal_adiabatic_results(homog,group)
 
-  integer, intent(in) :: &
+  integer,          intent(in) :: homog
-    homog, &
+  character(len=*), intent(in) :: group
-    instance, &
+#if defined(PETSc) || defined(DAMASK_HDF5)  
-    of
+  integer :: o, instance
   
-  real(pReal), dimension(sum(thermal_adiabatic_sizePostResult(:,instance))) :: &
+  instance  = thermal_typeInstance(homog)
-    postResults
 
-  integer :: &
+  outputsLoop: do o = 1,thermal_adiabatic_Noutput(instance)
-    o, c
- 
-  c = 0
- 
-  do o = 1,thermal_adiabatic_Noutput(instance)
      select case(thermal_adiabatic_outputID(o,instance))
     
       case (temperature_ID)
-         postResults(c+1) = temperature(homog)%p(of)
+        call results_writeDataset(group,temperature(homog)%p,'T',&
-         c = c + 1
+                                  'temperature','K')
     end select
-  enddo
+  enddo outputsLoop
+#endif
 
-end function thermal_adiabatic_postResults
+end subroutine thermal_adiabatic_results
 
 end module thermal_adiabatic

src/thermal_conduction.f90

@@ -7,6 +7,7 @@ module thermal_conduction
   use material
   use config
   use lattice
+  use results
   use crystallite
   use source_thermal_dissipation
   use source_thermal_externalheat
@@ -14,8 +15,6 @@ module thermal_conduction
   implicit none
   private
  
-  integer,                     dimension(:,:), allocatable, target, public :: &
-    thermal_conduction_sizePostResult                                                               !< size of each post result output
   character(len=64),           dimension(:,:), allocatable, target, public :: &
     thermal_conduction_output                                                                       !< name of each post result output
     
@@ -37,7 +36,7 @@ module thermal_conduction
     thermal_conduction_getSpecificHeat, &
     thermal_conduction_getMassDensity, &
     thermal_conduction_putTemperatureAndItsRate, &
-    thermal_conduction_postResults
+    thermal_conduction_results
 
 contains
 
@@ -60,7 +59,6 @@ subroutine thermal_conduction_init
   maxNinstance = count(thermal_type == THERMAL_conduction_ID)
   if (maxNinstance == 0) return
   
-  allocate(thermal_conduction_sizePostResult (maxval(homogenization_Noutput),maxNinstance),source=0)
   allocate(thermal_conduction_output         (maxval(homogenization_Noutput),maxNinstance))
            thermal_conduction_output = ''
   allocate(thermal_conduction_outputID       (maxval(homogenization_Noutput),maxNinstance),source=undefined_ID)
@@ -78,7 +76,6 @@ subroutine thermal_conduction_init
               thermal_conduction_Noutput(instance) = thermal_conduction_Noutput(instance) + 1
               thermal_conduction_outputID(thermal_conduction_Noutput(instance),instance) = temperature_ID
               thermal_conduction_output(thermal_conduction_Noutput(instance),instance) = outputs(i)
-              thermal_conduction_sizePostResult(thermal_conduction_Noutput(instance),instance) = 1
       end select
     enddo
  
@@ -86,7 +83,6 @@ subroutine thermal_conduction_init
  ! allocate state arrays
     sizeState = 0
     thermalState(section)%sizeState = sizeState
-    thermalState(section)%sizePostResults = sum(thermal_conduction_sizePostResult(:,instance))
     allocate(thermalState(section)%state0   (sizeState,NofMyHomog))
     allocate(thermalState(section)%subState0(sizeState,NofMyHomog))
     allocate(thermalState(section)%state    (sizeState,NofMyHomog))
@@ -265,31 +261,27 @@ end subroutine thermal_conduction_putTemperatureAndItsRate
  
 
 !--------------------------------------------------------------------------------------------------
-!> @brief return array of thermal results
+!> @brief writes results to HDF5 output file
 !--------------------------------------------------------------------------------------------------
-function thermal_conduction_postResults(homog,instance,of) result(postResults)
+subroutine thermal_conduction_results(homog,group)
 
-  integer,              intent(in) :: &
+  integer,          intent(in) :: homog
-    homog, &
+  character(len=*), intent(in) :: group
-    instance, &
+#if defined(PETSc) || defined(DAMASK_HDF5)  
-    of
+  integer :: o, instance
   
-  real(pReal), dimension(sum(thermal_conduction_sizePostResult(:,instance))) :: &
+  instance  = thermal_typeInstance(homog)
-    postResults
 
-  integer :: &
+  outputsLoop: do o = 1,thermal_conduction_Noutput(instance)
-    o, c
- 
-  c = 0
-  do o = 1,thermal_conduction_Noutput(instance)
      select case(thermal_conduction_outputID(o,instance))
     
       case (temperature_ID)
-         postResults(c+1) = temperature(homog)%p(of)
+        call results_writeDataset(group,temperature(homog)%p,'T',&
-         c = c + 1
+                                  'temperature','K')
     end select
-  enddo
+  enddo outputsLoop
+#endif
 
-end function thermal_conduction_postResults
+end subroutine thermal_conduction_results
 
 end module thermal_conduction