Merge remote-tracking branch 'origin/development' into simple-int-formatting

VERSION

@@ -1 +1 @@
-v2.0.3-1177-ga41871e2
+v2.0.3-1228-g3e269f04

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

@@ -66,7 +66,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 = []
@@ -132,7 +132,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

@@ -52,6 +52,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 = []
@@ -111,7 +112,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()

python/damask/dadf5.py

@@ -30,7 +30,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 +47,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 +173,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 +191,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 +209,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')
 
 
@@ -422,6 +439,76 @@ class DADF5():
         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 +540,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 +694,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 +753,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 +787,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.

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'),

src/grid/DAMASK_grid.f90

@@ -28,7 +28,6 @@ program DAMASK_spectral
  use grid_damage_spectral
  use grid_thermal_spectral
  use results
- use rotations
 
  implicit none
 
@@ -78,7 +77,6 @@ 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
@@ -189,6 +187,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 +243,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
 
@@ -295,14 +293,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
@@ -462,7 +458,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)
@@ -481,7 +477,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)

src/grid/grid_mech_FEM.f90

@@ -206,8 +206,7 @@ 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)
@@ -240,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
 !--------------------------------------------------------------------------------------------------
@@ -295,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(:,:,:,:) :: &
@@ -480,7 +480,7 @@ subroutine formResidual(da_local,x_local, &
             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)
@@ -496,7 +496,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

@@ -172,8 +172,7 @@ 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
@@ -210,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
 !--------------------------------------------------------------------------------------------------
@@ -267,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
@@ -299,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])
@@ -309,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
@@ -444,7 +444,7 @@ subroutine formResidual(in, F, &
             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)
@@ -469,7 +469,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
@@ -185,8 +186,7 @@ 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
@@ -225,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
 !--------------------------------------------------------------------------------------------------
@@ -286,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
@@ -322,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])
     
@@ -336,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), &
@@ -347,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
@@ -512,7 +513,7 @@ subroutine formResidual(in, FandF_tau, &
             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)
@@ -531,7 +532,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
 
 !--------------------------------------------------------------------------------------------------
@@ -549,7 +550,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)

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/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/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))))

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