Merge remote-tracking branch 'origin/development' into restructure-Orientation-2

PRIVATE

@@ -1 +1 @@
-Subproject commit 174ecac2d3ab7596bdb60184d6bb9e1a52cb7378
+Subproject commit 4ce625b4ac0da9d490620f8cf1694d0a057cfa47

examples/config/Phase_DisloUCLA_Tungsten.config

@@ -1,42 +0,0 @@
-[Tungsten]
-elasticity          hooke
-plasticity	    disloucla 
-
-(output)           edge_density
-(output)           dipole_density
-(output)           shear_rate_slip
-(output)           accumulated_shear_slip
-(output)           resolved_stress_slip
-(output)           threshold_stress_slip
-
-grainsize             2.7e-5        # Average grain size [m] 2.0e-5
-SolidSolutionStrength 0.0           # Strength due to elements in solid solution
-
-### Dislocation glide parameters ###
-#per family
-Nslip               12
-slipburgers         2.72e-10                # Burgers vector of slip system [m]
-rhoedge0            1.0e12                  # Initial edge dislocation density [m/m**3]
-rhoedgedip0         1.0                     # Initial edged dipole dislocation density [m/m**3]
-Qedge               2.61154e-19             # Activation energy for dislocation glide [J], 1.63 eV
-v0                  1                       # Initial glide velocity [m/s]
-p_slip              0.86                    # p-exponent in glide velocity
-q_slip              1.69                    # q-exponent in glide velocity
-tau_peierls         2.03e9                  # peierls stress [Pa]
-
-#mobility law
-kink_height	    2.567e-10                   # kink height sqrt(6)/3*lattice_parameter [m]
-omega		        9.1e11	                    # attemp frequency (from kMC paper) [s^(-1)]
-kink_width	    29.95e-10                   # kink pair width ~ 11 b (kMC paper) [m]
-dislolength	    78e-10	                    # dislocation length (ideally lambda) [m] initial value 11b
-friction_coeff  8.3e-5                      # friction coeff. B [Pa*s]
-
-#hardening
-dipoleformationfactor 0                     # to have hardening due to dipole formation off
-CLambdaSlip         10.0                    # Adj. parameter controlling dislocation mean free path
-D0                  4.0e-5                  # Vacancy diffusion prefactor [m**2/s]
-Qsd                 4.5e-19                 # Activation energy for climb [J]
-Catomicvolume       1.0                     # Adj. parameter controlling the atomic volume [in b]
-Cedgedipmindistance 1.0                     # Adj. parameter controlling the minimum dipole distance [in b]
-interaction_slipslip     0.009 0.72 0.009 0.05 0.05 0.06 0.09
-nonschmid_coefficients   0.938 0.71 4.43 0.0 0.0 0.0

examples/config/Phase_Dislotungsten_W.yaml

@@ -0,0 +1,26 @@
+type: dislotungsten
+
+N_sl: [12]
+
+rho_mob_0: [1.0e+9]
+rho_dip_0: [1.0]
+
+nu_a: [9.1e+11]
+b_sl: [2.72e-10]
+Delta_H_kp,0: [2.61154e-19]  # 1.63 eV, Delta_H0
+
+tau_Peierls: [2.03e+9]
+p_sl: [0.86]
+q_sl: [1.69]
+h: [2.566e-10]
+w: [2.992e-09]
+B: [8.3e-5]
+D_a: 1.0 # d_edge
+
+# climb (disabled)
+D_0: 0.0
+Q_cl: 0.0
+V_cl: [0.0]
+
+h_sl-sl: [0.009, 0.72, 0.009, 0.05, 0.05, 0.06, 0.09]
+a_nonSchmid: [0.938, 0.71, 4.43]

examples/config/Phase_Dislotwin_TWIP-Steel-FeMnC.yaml

@@ -6,7 +6,7 @@ b_sl: [2.56e-10]
 rho_mob_0: [1.0e+12]
 rho_dip_0: [1.0]
 v_0: [1.0e+4]
-Q_s: [3.7e-19]
+Q_sl: [3.7e-19]
 p_sl: [1.0]
 q_sl: [1.0]
 tau_0: [1.5e+8]

examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml

@@ -11,7 +11,7 @@ b_sl: [2.49e-10, 2.49e-10]
 rho_mob_0: [2.81e12, 2.8e+12]
 rho_dip_0: [1.0, 1.0] # not given
 v_0: [1.4e+3, 1.4e+3]
-Q_s: [1.57e-19, 1.57e-19] # Delta_F
+Q_sl: [1.57e-19, 1.57e-19] # Delta_F
 tau_0: [454.e+6, 454.e+6]
 p_sl: [0.325, 0.325]
 q_sl: [1.55, 1.55]
@@ -19,6 +19,5 @@ i_sl: [23.3, 23.3]
 D_a: 7.4 # C_anni
 B: [0.001, 0.001]
 h_sl-sl: [0.1, 0.72, 0.1, 0.053, 0.053, 0.073, 0.137, 0.72, 0.72, 0.053, 0.053, 0.053, 0.053, 0.073, 0.073, 0.073, 0.073, 0.073, 0.073, 0.137, 0.073, 0.073, 0.137, 0.073]
-D_0: 4.0e-05
 Q_cl: 5.4e-19 # no recovery!
 D: 40.e-6 # estimated

examples/config/phase/mechanical/plastic/phenopowerlaw_Mg.yaml

@@ -1,20 +1,29 @@
-N_sl: [3, 3, 0, 6, 0, 6]
-N_tw: [6, 0, 0, 6]
-h_0_tw-tw: 50.0e+6
-h_0_sl-sl: 500.0e+6
-h_0_tw-sl: 150.0e+6
-h_sl-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
-h_tw-tw: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
-h_sl-tw: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
-h_tw-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
-output: [xi_sl, xi_tw]
 type: phenopowerlaw
-xi_0_sl: [10.e+6, 55.e+6, 0., 60.e+6, 0., 60.e+6]
+references:
+  - F. Wang et al.,
+    Acta Materialia 80:77-93, 2014,
+    https://doi.org/10.1016/j.actamat.2014.07.048
+
+output: [xi_sl, xi_tw]
+
+N_sl: [3, 3, 0, 6, 0, 6]   # basal, 1. prism, -, 1. pyr<a>, -, 2. pyr<c+a>
+N_tw: [6, 0, 6]   # tension, -, compression
+
+xi_0_sl:   [10.e+6,  55.e+6, 0.,  60.e+6, 0., 60.e+6]
 xi_inf_sl: [40.e+6, 135.e+6, 0., 150.e+6, 0., 150.e+6]
-xi_0_tw: [40.e+6, 0., 0., 60.e+6]
+xi_0_tw: [40.e+6, 0., 60.e+6]
+
 a_sl: 2.25
 dot_gamma_0_sl: 0.001
 dot_gamma_0_tw: 0.001
 n_sl: 20
 n_tw: 20
 f_sat_sl-tw: 10.0
+
+h_0_sl-sl: 500.0e+6
+h_0_tw-tw: 50.0e+6
+h_0_tw-sl: 150.0e+6
+h_sl-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
+h_tw-tw: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
+h_tw-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
+h_sl-tw: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml

@@ -7,14 +7,17 @@ references:
     Acta Materialia 132:598-610, 2017,
     https://doi.org/10.1016/j.actamat.2017.05.015
 output: [gamma_sl]
-N_sl: [3, 3, 0, 0, 12]
+
+N_sl: [3, 3, 0, 0, 12]   # basal, 1. prism, -, -, 2. pyr<c+a>
 n_sl: 20
 a_sl: 2.0
 dot_gamma_0_sl: 0.001
 h_0_sl-sl: 200.e+6
+
 # C. Zambaldi et al.:
 xi_0_sl:   [349.e+6, 150.e+6, 0.0, 0.0, 1107.e+6]
 xi_inf_sl: [568.e+6, 150.e+7, 0.0, 0.0, 3420.e+6]
 # L. Wang et al. :
 # xi_0_sl:   [127.e+6, 96.e+6, 0.0, 0.0, 240.e+6]
+
 h_sl-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

examples/mesh/tensionZ.load

@@ -1,11 +1,11 @@
 #initial elastic step
-$Loadcase 1 time 0.0005 incs 1 frequency 5
+$Loadcase 1 t 0.0005 N 1 f_out 5
  Face 1 X 0.01
  Face 2 X 0.0
  Face 2 Y 0.0
  Face 2 Z 0.0
 $EndLoadcase
-$Loadcase 2 time 10.0 incs 200 frequency 5
+$Loadcase 2 t 10.0 N 200 f_out 5
  Face 1 X 0.01
  Face 2 X 0.0
  Face 2 Y 0.0

python/damask/VERSION

@@ -1 +1 @@
-v3.0.0-alpha4-137-gb69b85754
+v3.0.0-alpha4-235-g2635bb012

python/damask/_result.py

@@ -99,8 +99,10 @@ class Result:
             self.version_major = f.attrs['DADF5_version_major']
             self.version_minor = f.attrs['DADF5_version_minor']
 
-            if self.version_major != 0 or not 12 <= self.version_minor <= 13:
+            if self.version_major != 0 or not 12 <= self.version_minor <= 14:
                 raise TypeError(f'Unsupported DADF5 version {self.version_major}.{self.version_minor}')
+            if self.version_major == 0 and self.version_minor < 14:
+                self.export_setup = None
 
             self.structured = 'cells' in f['geometry'].attrs.keys()
 
@@ -1383,7 +1385,7 @@ class Result:
 
     def export_XDMF(self,output='*'):
         """
-        Write XDMF file to directly visualize data in DADF5 file.
+        Write XDMF file to directly visualize data from DADF5 file.
 
         The XDMF format is only supported for structured grids
         with single phase and single constituent.
@@ -1736,3 +1738,32 @@ class Result:
         if flatten: r = util.dict_flatten(r)
 
         return None if (type(r) == dict and r == {}) else r
+
+
+    def export_setup(self,output='*',overwrite=False):
+        """
+        Export configuration files.
+
+        Parameters
+        ----------
+        output : (list of) str, optional
+            Names of the datasets to export to the file.
+            Defaults to '*', in which case all datasets are exported.
+        overwrite : boolean, optional
+            Overwrite existing configuration files.
+            Defaults to False.
+
+        """
+        def export(name,obj,output,overwrite):
+            if type(obj) == h5py.Dataset and  _match(output,[name]):
+                d = obj.attrs['description'] if h5py3 else obj.attrs['description'].decode()
+                if not Path(name).exists() or overwrite:
+                    with open(name,'w') as f_out: f_out.write(obj[()].decode())
+                    print(f"Exported {d} to '{name}'.")
+                else:
+                    print(f"'{name}' exists, {d} not exported.")
+            elif type(obj) == h5py.Group:
+                os.makedirs(name, exist_ok=True)
+
+        with h5py.File(self.fname,'r') as f_in:
+            f_in['setup'].visititems(partial(export,output=output,overwrite=overwrite))

src/CPFEM.f90

@@ -4,22 +4,23 @@
 !> @brief CPFEM engine
 !--------------------------------------------------------------------------------------------------
 module CPFEM
+  use DAMASK_interface
   use prec
-  use math
+  use IO
-  use rotations
   use YAML_types
   use YAML_parse
-  use discretization_marc
-  use material
-  use config
-  use homogenization
-  use IO
-  use discretization
-  use DAMASK_interface
   use HDF5_utilities
   use results
+  use config
+  use math
+  use rotations
   use lattice
+  use material
   use phase
+  use homogenization
+
+  use discretization
+  use discretization_marc
 
   implicit none
   private
@@ -68,7 +69,7 @@ contains
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief call all module initializations
+!> @brief Initialize all modules.
 !--------------------------------------------------------------------------------------------------
 subroutine CPFEM_initAll
 
@@ -77,13 +78,13 @@ subroutine CPFEM_initAll
   call IO_init
   call YAML_types_init
   call YAML_parse_init
+  call HDF5_utilities_init
+  call results_init(.false.)
   call config_init
   call math_init
   call rotations_init
-  call HDF5_utilities_init
-  call results_init(.false.)
-  call discretization_marc_init
   call lattice_init
+  call discretization_marc_init
   call material_init(.false.)
   call phase_init
   call homogenization_init

src/CPFEM2.f90

@@ -4,28 +4,29 @@
 !> @brief needs a good name and description
 !--------------------------------------------------------------------------------------------------
 module CPFEM2
-  use prec
   use parallelization
-  use config
+  use DAMASK_interface
-  use math
+  use prec
-  use rotations
+  use IO
   use YAML_types
   use YAML_parse
-  use material
-  use lattice
-  use IO
-  use base64
-  use DAMASK_interface
-  use discretization
   use HDF5
   use HDF5_utilities
   use results
-  use homogenization
+  use config
+  use math
+  use rotations
+  use lattice
+  use material
   use phase
+  use homogenization
+
+  use discretization
 #if   defined(MESH)
   use FEM_quadrature
   use discretization_mesh
 #elif defined(GRID)
+  use base64
   use discretization_grid
 #endif
 
@@ -36,7 +37,7 @@ contains
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief call all module initializations
+!> @brief Initialize all modules.
 !--------------------------------------------------------------------------------------------------
 subroutine CPFEM_initAll
 
@@ -44,18 +45,19 @@ subroutine CPFEM_initAll
   call DAMASK_interface_init                                                                        ! Spectral and FEM interface to commandline
   call prec_init
   call IO_init
-  call base64_init
+#if   defined(MESH)
-#ifdef MESH
   call FEM_quadrature_init
+#elif defined(GRID)
+   call base64_init
 #endif
   call YAML_types_init
   call YAML_parse_init
+  call HDF5_utilities_init
+  call results_init(restart=interface_restartInc>0)
   call config_init
   call math_init
   call rotations_init
   call lattice_init
-  call HDF5_utilities_init
-  call results_init(restart=interface_restartInc>0)
 #if   defined(MESH)
   call discretization_mesh_init(restart=interface_restartInc>0)
 #elif defined(GRID)

src/DAMASK_interface.f90

@@ -139,10 +139,10 @@ subroutine DAMASK_interface_init
         print'(a)',  ' Optional arguments:'
         print'(/,a)','   --workingdirectory PathToWorkingDirectory'
         print'(a)',  '        Specifies the working directory and overwrites the default ./'
-        print'(a)',  '        Make sure the file "material.config" exists in the working'
+        print'(a)',  '        Make sure the file "material.yaml" exists in the working'
         print'(a)',  '            directory.'
-        print'(a)',  '        For further configuration place "numerics.config"'
+        print'(a)',  '        For further configuration place "numerics.yaml"'
-        print'(a)','            and "debug.config" in that directory.'
+        print'(a)','            and "debug.yaml" in that directory.'
         print'(/,a)','   --restart N'
         print'(a)',  '        Reads in increment N and continues with calculating'
         print'(a)',  '            increment N+1 based on this.'

src/HDF5_utilities.f90

@@ -71,6 +71,7 @@ module HDF5_utilities
     module procedure HDF5_addAttribute_str
     module procedure HDF5_addAttribute_int
     module procedure HDF5_addAttribute_real
+    module procedure HDF5_addAttribute_str_array
     module procedure HDF5_addAttribute_int_array
     module procedure HDF5_addAttribute_real_array
   end interface HDF5_addAttribute
@@ -84,6 +85,7 @@ module HDF5_utilities
     HDF5_utilities_init, &
     HDF5_read, &
     HDF5_write, &
+    HDF5_write_str, &
     HDF5_addAttribute, &
     HDF5_addGroup, &
     HDF5_openGroup, &
@@ -127,10 +129,11 @@ end subroutine HDF5_utilities_init
 !--------------------------------------------------------------------------------------------------
 !> @brief open and initializes HDF5 output file
 !--------------------------------------------------------------------------------------------------
-integer(HID_T) function HDF5_openFile(fileName,mode)
+integer(HID_T) function HDF5_openFile(fileName,mode,parallel)
 
   character(len=*), intent(in)           :: fileName
   character,        intent(in), optional :: mode
+  logical,          intent(in), optional :: parallel
 
   character                              :: m
   integer(HID_T)                         :: plist_id
@@ -147,7 +150,11 @@ integer(HID_T) function HDF5_openFile(fileName,mode)
   if(hdferr < 0) error stop 'HDF5 error'
 
 #ifdef PETSC
-  call h5pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL, hdferr)
+  if (present(parallel)) then
+    if (parallel) call h5pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL, hdferr)
+  else
+    call h5pset_fapl_mpio_f(plist_id, PETSC_COMM_WORLD, MPI_INFO_NULL, hdferr)
+  endif
   if(hdferr < 0) error stop 'HDF5 error'
 #endif
 
@@ -270,7 +277,7 @@ end subroutine HDF5_closeGroup
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief check whether a group or a dataset exists
+!> @brief Check whether a group or a dataset exists.
 !--------------------------------------------------------------------------------------------------
 logical function HDF5_objectExists(loc_id,path)
 
@@ -280,6 +287,7 @@ logical function HDF5_objectExists(loc_id,path)
   integer :: hdferr
   character(len=:), allocatable :: p
 
+
   if (present(path)) then
     p = trim(path)
   else
@@ -298,7 +306,7 @@ end function HDF5_objectExists
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief adds a string attribute to the path given relative to the location
+!> @brief Add a string attribute to the path given relative to the location.
 !--------------------------------------------------------------------------------------------------
 subroutine HDF5_addAttribute_str(loc_id,attrLabel,attrValue,path)
 
@@ -313,6 +321,7 @@ subroutine HDF5_addAttribute_str(loc_id,attrLabel,attrValue,path)
   character(len=:,kind=C_CHAR), allocatable,target :: attrValue_
   type(c_ptr), target, dimension(1) :: ptr
 
+
   if (present(path)) then
     p = trim(path)
   else
@@ -326,6 +335,7 @@ subroutine HDF5_addAttribute_str(loc_id,attrLabel,attrValue,path)
   if(hdferr < 0) error stop 'HDF5 error'
   call h5tcopy_f(H5T_STRING, type_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
+
   call h5aexists_by_name_f(loc_id,trim(p),attrLabel,attrExists,hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
   if (attrExists) then
@@ -336,6 +346,7 @@ subroutine HDF5_addAttribute_str(loc_id,attrLabel,attrValue,path)
   if(hdferr < 0) error stop 'HDF5 error'
   call h5awrite_f(attr_id, type_id, c_loc(ptr(1)), hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
+
   call h5aclose_f(attr_id,hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
   call h5tclose_f(type_id,hdferr)
@@ -347,7 +358,7 @@ end subroutine HDF5_addAttribute_str
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief adds a integer attribute to the path given relative to the location
+!> @brief Add an integer attribute to the path given relative to the location.
 !--------------------------------------------------------------------------------------------------
 subroutine HDF5_addAttribute_int(loc_id,attrLabel,attrValue,path)
 
@@ -361,6 +372,7 @@ subroutine HDF5_addAttribute_int(loc_id,attrLabel,attrValue,path)
   logical        :: attrExists
   character(len=:), allocatable :: p
 
+
   if (present(path)) then
     p = trim(path)
   else
@@ -369,6 +381,7 @@ subroutine HDF5_addAttribute_int(loc_id,attrLabel,attrValue,path)
 
   call h5screate_f(H5S_SCALAR_F,space_id,hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
+
   call h5aexists_by_name_f(loc_id,trim(p),attrLabel,attrExists,hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
   if (attrExists) then
@@ -379,6 +392,7 @@ subroutine HDF5_addAttribute_int(loc_id,attrLabel,attrValue,path)
   if(hdferr < 0) error stop 'HDF5 error'
   call h5awrite_f(attr_id, H5T_NATIVE_INTEGER, attrValue, int([1],HSIZE_T), hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
+
   call h5aclose_f(attr_id,hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
   call h5sclose_f(space_id,hdferr)
@@ -388,7 +402,7 @@ end subroutine HDF5_addAttribute_int
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief adds a integer attribute to the path given relative to the location
+!> @brief Add a real attribute to the path given relative to the location.
 !--------------------------------------------------------------------------------------------------
 subroutine HDF5_addAttribute_real(loc_id,attrLabel,attrValue,path)
 
@@ -402,6 +416,7 @@ subroutine HDF5_addAttribute_real(loc_id,attrLabel,attrValue,path)
   logical        :: attrExists
   character(len=:), allocatable :: p
 
+
   if (present(path)) then
     p = trim(path)
   else
@@ -410,6 +425,7 @@ subroutine HDF5_addAttribute_real(loc_id,attrLabel,attrValue,path)
 
   call h5screate_f(H5S_SCALAR_F,space_id,hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
+
   call h5aexists_by_name_f(loc_id,trim(p),attrLabel,attrExists,hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
   if (attrExists) then
@@ -420,6 +436,7 @@ subroutine HDF5_addAttribute_real(loc_id,attrLabel,attrValue,path)
   if(hdferr < 0) error stop 'HDF5 error'
   call h5awrite_f(attr_id, H5T_NATIVE_DOUBLE, attrValue, int([1],HSIZE_T), hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
+
   call h5aclose_f(attr_id,hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
   call h5sclose_f(space_id,hdferr)
@@ -429,7 +446,67 @@ end subroutine HDF5_addAttribute_real
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief adds a integer attribute to the path given relative to the location
+!> @brief Add a string array attribute to the path given relative to the location.
+!--------------------------------------------------------------------------------------------------
+subroutine HDF5_addAttribute_str_array(loc_id,attrLabel,attrValue,path)
+
+  integer(HID_T),   intent(in)               :: loc_id
+  character(len=*), intent(in)               :: attrLabel
+  character(len=*), intent(in), dimension(:) :: attrValue
+  character(len=*), intent(in), optional     :: path
+
+  integer(HID_T)                :: attr_id, space_id, filetype_id, memtype_id
+  integer                       :: hdferr
+  logical                       :: attrExists
+  character(len=:), allocatable :: p
+  type(C_PTR) :: f_ptr
+  character(len=:), allocatable, dimension(:), target :: attrValue_
+
+
+  if (present(path)) then
+    p = trim(path)
+  else
+    p = '.'
+  endif
+
+  attrValue_ = attrValue
+
+  call h5tcopy_f(H5T_C_S1,filetype_id,hdferr)
+  if(hdferr < 0) error stop 'HDF5 error'
+  call h5tset_size_f(filetype_id, int(len(attrValue_)+1,C_SIZE_T),hdferr)
+  if(hdferr < 0) error stop 'HDF5 error'
+  call h5tcopy_f(H5T_FORTRAN_S1, memtype_id, hdferr)
+  if(hdferr < 0) error stop 'HDF5 error'
+  call h5tset_size_f(memtype_id, int(len(attrValue_),C_SIZE_T), hdferr)
+  if(hdferr < 0) error stop 'HDF5 error'
+  call h5screate_simple_f(1,shape(attrValue_,kind=HSIZE_T),space_id, hdferr)
+  if(hdferr < 0) error stop 'HDF5 error'
+
+  call h5aexists_by_name_f(loc_id,trim(p),attrLabel,attrExists,hdferr)
+  if (attrExists) then
+    call h5adelete_by_name_f(loc_id, trim(p), attrLabel, hdferr)
+    if(hdferr < 0) error stop 'HDF5 error'
+  endif
+  call h5acreate_by_name_f(loc_id,trim(p),trim(attrLabel),filetype_id,space_id,attr_id,hdferr)
+  if(hdferr < 0) error stop 'HDF5 error'
+  f_ptr = c_loc(attrValue_)
+  call h5awrite_f(attr_id, memtype_id, f_ptr, hdferr)
+  if(hdferr < 0) error stop 'HDF5 error'
+
+  call h5tclose_f(memtype_id,hdferr)
+  if(hdferr < 0) error stop 'HDF5 error'
+  call h5tclose_f(filetype_id,hdferr)
+  if(hdferr < 0) error stop 'HDF5 error'
+  call h5aclose_f(attr_id,hdferr)
+  if(hdferr < 0) error stop 'HDF5 error'
+  call h5sclose_f(space_id,hdferr)
+  if(hdferr < 0) error stop 'HDF5 error'
+
+end subroutine HDF5_addAttribute_str_array
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Add an integer array attribute to the path given relative to the location.
 !--------------------------------------------------------------------------------------------------
 subroutine HDF5_addAttribute_int_array(loc_id,attrLabel,attrValue,path)
 
@@ -444,6 +521,7 @@ subroutine HDF5_addAttribute_int_array(loc_id,attrLabel,attrValue,path)
   logical                       :: attrExists
   character(len=:), allocatable :: p
 
+
   if (present(path)) then
     p = trim(path)
   else
@@ -454,6 +532,7 @@ subroutine HDF5_addAttribute_int_array(loc_id,attrLabel,attrValue,path)
 
   call h5screate_simple_f(1, array_size, space_id, hdferr, array_size)
   if(hdferr < 0) error stop 'HDF5 error'
+
   call h5aexists_by_name_f(loc_id,trim(p),attrLabel,attrExists,hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
   if (attrExists) then
@@ -464,6 +543,7 @@ subroutine HDF5_addAttribute_int_array(loc_id,attrLabel,attrValue,path)
   if(hdferr < 0) error stop 'HDF5 error'
   call h5awrite_f(attr_id, H5T_NATIVE_INTEGER, attrValue, array_size, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
+
   call h5aclose_f(attr_id,hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
   call h5sclose_f(space_id,hdferr)
@@ -473,7 +553,7 @@ end subroutine HDF5_addAttribute_int_array
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief adds a real attribute to the path given relative to the location
+!> @brief Add a real array attribute to the path given relative to the location.
 !--------------------------------------------------------------------------------------------------
 subroutine HDF5_addAttribute_real_array(loc_id,attrLabel,attrValue,path)
 
@@ -488,6 +568,7 @@ subroutine HDF5_addAttribute_real_array(loc_id,attrLabel,attrValue,path)
   logical                       :: attrExists
   character(len=:), allocatable :: p
 
+
   if (present(path)) then
     p = trim(path)
   else
@@ -498,6 +579,7 @@ subroutine HDF5_addAttribute_real_array(loc_id,attrLabel,attrValue,path)
 
   call h5screate_simple_f(1, array_size, space_id, hdferr, array_size)
   if(hdferr < 0) error stop 'HDF5 error'
+
   call h5aexists_by_name_f(loc_id,trim(p),attrLabel,attrExists,hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
   if (attrExists) then
@@ -508,6 +590,7 @@ subroutine HDF5_addAttribute_real_array(loc_id,attrLabel,attrValue,path)
   if(hdferr < 0) error stop 'HDF5 error'
   call h5awrite_f(attr_id, H5T_NATIVE_DOUBLE, attrValue, array_size, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
+
   call h5aclose_f(attr_id,hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
   call h5sclose_f(space_id,hdferr)
@@ -517,7 +600,7 @@ end subroutine HDF5_addAttribute_real_array
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief set link to object in results file
+!> @brief Set link to object in results file.
 !--------------------------------------------------------------------------------------------------
 subroutine HDF5_setLink(loc_id,target_name,link_name)
 
@@ -549,7 +632,7 @@ subroutine HDF5_read_real1(dataset,loc_id,datasetName,parallel)
   logical, intent(in), optional :: parallel                                                         !< dataset is distributed over multiple processes
 
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id, aplist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &                                            ! ToDo: Fortran 2018 size(shape(A)) = rank(A)
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -589,7 +672,7 @@ subroutine HDF5_read_real2(dataset,loc_id,datasetName,parallel)
   logical, intent(in), optional :: parallel                                                         !< dataset is distributed over multiple processes
 
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id, aplist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -629,7 +712,7 @@ subroutine HDF5_read_real3(dataset,loc_id,datasetName,parallel)
   logical, intent(in), optional :: parallel                                                         !< dataset is distributed over multiple processes
 
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id, aplist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -669,7 +752,7 @@ subroutine HDF5_read_real4(dataset,loc_id,datasetName,parallel)
   logical, intent(in), optional :: parallel                                                         !< dataset is distributed over multiple processes
 
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id, aplist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -709,7 +792,7 @@ subroutine HDF5_read_real5(dataset,loc_id,datasetName,parallel)
   logical, intent(in), optional :: parallel                                                         !< dataset is distributed over multiple processes
 
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id, aplist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -749,7 +832,7 @@ subroutine HDF5_read_real6(dataset,loc_id,datasetName,parallel)
   logical, intent(in), optional :: parallel                                                         !< dataset is distributed over multiple processes
 
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id, aplist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -789,7 +872,7 @@ subroutine HDF5_read_real7(dataset,loc_id,datasetName,parallel)
   logical, intent(in), optional :: parallel                                                         !< dataset is distributed over multiple processes
 
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id, aplist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -831,7 +914,7 @@ subroutine HDF5_read_int1(dataset,loc_id,datasetName,parallel)
 
 
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id, aplist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -871,7 +954,7 @@ subroutine HDF5_read_int2(dataset,loc_id,datasetName,parallel)
   logical, intent(in), optional :: parallel                                                         !< dataset is distributed over multiple processes
 
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id, aplist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -911,7 +994,7 @@ subroutine HDF5_read_int3(dataset,loc_id,datasetName,parallel)
   logical, intent(in), optional :: parallel                                                         !< dataset is distributed over multiple processes
 
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id, aplist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -951,7 +1034,7 @@ subroutine HDF5_read_int4(dataset,loc_id,datasetName,parallel)
   logical, intent(in), optional :: parallel                                                         !< dataset is distributed over multiple processes
 
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id, aplist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -991,7 +1074,7 @@ subroutine HDF5_read_int5(dataset,loc_id,datasetName,parallel)
   logical, intent(in), optional :: parallel                                                         !< dataset is distributed over multiple processes
 
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id, aplist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -1031,7 +1114,7 @@ subroutine HDF5_read_int6(dataset,loc_id,datasetName,parallel)
   logical, intent(in), optional :: parallel                                                         !< dataset is distributed over multiple processes
 
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id, aplist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -1071,7 +1154,7 @@ subroutine HDF5_read_int7(dataset,loc_id,datasetName,parallel)
   logical, intent(in), optional :: parallel                                                         !< dataset is distributed over multiple processes
 
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id, aplist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -1114,7 +1197,7 @@ subroutine HDF5_write_real1(dataset,loc_id,datasetName,parallel)
 
   integer :: hdferr
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -1155,7 +1238,7 @@ subroutine HDF5_write_real2(dataset,loc_id,datasetName,parallel)
 
   integer :: hdferr
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -1196,7 +1279,7 @@ subroutine HDF5_write_real3(dataset,loc_id,datasetName,parallel)
 
   integer :: hdferr
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -1237,7 +1320,7 @@ subroutine HDF5_write_real4(dataset,loc_id,datasetName,parallel)
 
   integer :: hdferr
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -1279,7 +1362,7 @@ subroutine HDF5_write_real5(dataset,loc_id,datasetName,parallel)
 
   integer :: hdferr
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -1320,7 +1403,7 @@ subroutine HDF5_write_real6(dataset,loc_id,datasetName,parallel)
 
   integer :: hdferr
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -1361,7 +1444,7 @@ subroutine HDF5_write_real7(dataset,loc_id,datasetName,parallel)
 
   integer :: hdferr
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -1390,6 +1473,48 @@ subroutine HDF5_write_real7(dataset,loc_id,datasetName,parallel)
 end subroutine HDF5_write_real7
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief Write dataset of type string (scalar).
+!> @details Not collective, must be called by one process at at time.
+!--------------------------------------------------------------------------------------------------
+subroutine HDF5_write_str(dataset,loc_id,datasetName)
+
+  character(len=*), intent(in) :: dataset
+  integer(HID_T),   intent(in) :: loc_id
+  character(len=*), intent(in) :: datasetName                                                      !< name of the dataset in the file
+
+  integer(HID_T)  :: filetype_id, space_id, dataset_id
+  integer :: hdferr
+  character(len=len_trim(dataset)+1,kind=C_CHAR), dimension(1), target :: dataset_
+  type(C_PTR), target, dimension(1) :: ptr
+
+
+  dataset_(1) = trim(dataset)//C_NULL_CHAR
+  ptr(1) = c_loc(dataset_(1))
+
+  call h5tcopy_f(H5T_STRING, filetype_id, hdferr)
+  if(hdferr < 0) error stop 'HDF5 error'
+  call h5tset_size_f(filetype_id, int(len(dataset_),HSIZE_T), hdferr)
+  if(hdferr < 0) error stop 'HDF5 error'
+
+  call h5screate_f(H5S_SCALAR_F, space_id, hdferr)
+  if(hdferr < 0) error stop 'HDF5 error'
+  call h5dcreate_f(loc_id, datasetName, H5T_STRING, space_id, dataset_id, hdferr)
+  if(hdferr < 0) error stop 'HDF5 error'
+
+  call h5dwrite_f(dataset_id, H5T_STRING, c_loc(ptr), hdferr)
+  if(hdferr < 0) error stop 'HDF5 error'
+
+  call h5dclose_f(dataset_id, hdferr)
+  if(hdferr < 0) error stop 'HDF5 error'
+  call h5sclose_f(space_id, hdferr)
+  if(hdferr < 0) error stop 'HDF5 error'
+  call h5tclose_f(filetype_id, hdferr)
+  if(hdferr < 0) error stop 'HDF5 error'
+
+end subroutine HDF5_write_str
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief write dataset of type integer with 1 dimension
 !--------------------------------------------------------------------------------------------------
@@ -1403,7 +1528,7 @@ subroutine HDF5_write_int1(dataset,loc_id,datasetName,parallel)
 
   integer :: hdferr
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -1444,7 +1569,7 @@ subroutine HDF5_write_int2(dataset,loc_id,datasetName,parallel)
 
   integer :: hdferr
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -1485,7 +1610,7 @@ subroutine HDF5_write_int3(dataset,loc_id,datasetName,parallel)
 
   integer :: hdferr
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -1526,7 +1651,7 @@ subroutine HDF5_write_int4(dataset,loc_id,datasetName,parallel)
 
   integer :: hdferr
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -1567,7 +1692,7 @@ subroutine HDF5_write_int5(dataset,loc_id,datasetName,parallel)
 
   integer :: hdferr
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -1608,7 +1733,7 @@ subroutine HDF5_write_int6(dataset,loc_id,datasetName,parallel)
 
   integer :: hdferr
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -1649,7 +1774,7 @@ subroutine HDF5_write_int7(dataset,loc_id,datasetName,parallel)
 
   integer :: hdferr
   integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id
-  integer(HSIZE_T), dimension(size(shape(dataset))) :: &
+  integer(HSIZE_T), dimension(rank(dataset)) :: &
     myStart, &
     myShape, &                                                                                      !< shape of the dataset (this process)
     totalShape                                                                                      !< shape of the dataset (all processes)
@@ -1795,7 +1920,7 @@ subroutine initialize_write(dset_id, filespace_id, memspace_id, plist_id, &
   integer(HSIZE_T), parameter :: chunkSize = 1024_HSIZE_T**2/8_HSIZE_T
 
 !-------------------------------------------------------------------------------------------------
-! creating a property list for transfer properties (is collective when reading in parallel)
+! creating a property list for transfer properties (is collective when writing in parallel)
   call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 #ifdef PETSC

src/IO.f90

@@ -16,7 +16,8 @@ module IO
   private
 
   character(len=*), parameter, public :: &
-    IO_WHITESPACE = achar(44)//achar(32)//achar(9)//achar(10)//achar(13)                            !< whitespace characters
+    IO_WHITESPACE = achar(44)//achar(32)//achar(9)//achar(10)//achar(13), &                         !< whitespace characters
+    IO_QUOTES  = "'"//'"'
   character, parameter, public :: &
     IO_EOL = new_line('DAMASK'), &                                                                  !< end of line character
     IO_COMMENT = '#'
@@ -119,27 +120,28 @@ function IO_read(fileName) result(fileContent)
   character(len=:), allocatable :: fileContent
 
   integer ::  &
-    fileLength, &
     fileUnit, &
     myStat
+  integer(pI64) ::  &
+    fileLength
 
 
   inquire(file = fileName, size=fileLength)
   open(newunit=fileUnit, file=fileName, access='stream',&
        status='old', position='rewind', action='read',iostat=myStat)
-  if(myStat /= 0) call IO_error(100,ext_msg=trim(fileName))
+  if (myStat /= 0) call IO_error(100,ext_msg=trim(fileName))
   allocate(character(len=fileLength)::fileContent)
-  if(fileLength==0) then
+  if (fileLength==0) then
     close(fileUnit)
     return
   endif
 
   read(fileUnit,iostat=myStat) fileContent
-  if(myStat /= 0) call IO_error(102,ext_msg=trim(fileName))
+  if (myStat /= 0) call IO_error(102,ext_msg=trim(fileName))
   close(fileUnit)
 
   if (scan(fileContent(:index(fileContent,LF)),CR//LF) /= 0) fileContent = CRLF2LF(fileContent)
-  if(fileContent(fileLength:fileLength) /= IO_EOL)           fileContent = fileContent//IO_EOL      ! ensure EOL@EOF
+  if (fileContent(fileLength:fileLength) /= IO_EOL)          fileContent = fileContent//IO_EOL      ! ensure EOL@EOF
 
 end function IO_read
 
@@ -494,6 +496,8 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
       msg = '--- expected after YAML file header'
     case (709)
       msg = 'Length mismatch'
+    case (710)
+      msg = 'Closing quotation mark missing in string'
 
 !-------------------------------------------------------------------------------------------------
 ! errors related to the grid solver

src/Marc/discretization_Marc.f90

@@ -216,7 +216,13 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,materialAt)
     mapElemSet                                                                                      !< list of elements in elementSet
 
 
-  inputFile = IO_readlines(trim(getSolverJobName())//trim(InputFileExtension))
+  call results_openJobFile
+  call results_writeDataset_str(IO_read(trim(getSolverJobName())//InputFileExtension), 'setup', &
+                                        trim(getSolverJobName())//InputFileExtension, &
+                                        'MSC.Marc input deck')
+  call results_closeJobFile
+
+  inputFile = IO_readlines(trim(getSolverJobName())//InputFileExtension)
   call inputRead_fileFormat(fileFormatVersion, &
                             inputFile)
   call inputRead_tableStyles(initialcondTableStyle,hypoelasticTableStyle, &

src/YAML_parse.f90

@@ -14,7 +14,7 @@ module YAML_parse
 
   public :: &
     YAML_parse_init, &
-    YAML_parse_file
+    YAML_parse_str
 
 contains
 
@@ -29,16 +29,16 @@ end subroutine YAML_parse_init
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Parse a YAML file into a a structure of nodes.
+!> @brief Parse a YAML string into a a structure of nodes.
 !--------------------------------------------------------------------------------------------------
-function YAML_parse_file(fname) result(node)
+function YAML_parse_str(str) result(node)
 
-  character(len=*), intent(in) :: fname
+  character(len=*), intent(in) :: str
   class (tNode), pointer :: node
 
-  node => parse_flow(to_flow(IO_read(fname)))
+  node => parse_flow(to_flow(str))
 
-end function YAML_parse_file
+end function YAML_parse_str
 
 
 !--------------------------------------------------------------------------------------------------
@@ -71,8 +71,8 @@ recursive function parse_flow(YAML_flow) result(node)
       s = e
       d = s + scan(flow_string(s+1:),':')
       e = d + find_end(flow_string(d+1:),'}')
-
       key = trim(adjustl(flow_string(s+1:d-1)))
+      if(quotedString(key)) key = key(2:len(key)-1)
       myVal => parse_flow(flow_string(d+1:e-1))                                                     ! parse items (recursively)
 
       select type (node)
@@ -97,7 +97,11 @@ recursive function parse_flow(YAML_flow) result(node)
     allocate(tScalar::node)
       select type (node)
         class is (tScalar)
-          node = trim(adjustl(flow_string))
+          if(quotedString(flow_string)) then
+            node = trim(adjustl(flow_string(2:len(flow_string)-1)))
+          else
+            node = trim(adjustl(flow_string))
+          endif
       end select
   endif
 
@@ -119,18 +123,38 @@ integer function find_end(str,e_char)
 
   N_sq = 0
   N_cu = 0
-  do i = 1, len_trim(str)
+  i = 1
+  do while(i<=len_trim(str))
+    if (scan(str(i:i),IO_QUOTES) == 1)  i = i + scan(str(i+1:),str(i:i))
     if (N_sq==0 .and. N_cu==0 .and. scan(str(i:i),e_char//',') == 1) exit
     N_sq = N_sq + merge(1,0,str(i:i) == '[')
     N_cu = N_cu + merge(1,0,str(i:i) == '{')
     N_sq = N_sq - merge(1,0,str(i:i) == ']')
     N_cu = N_cu - merge(1,0,str(i:i) == '}')
+    i = i + 1
   enddo
   find_end = i
 
 end function find_end
 
 
+!--------------------------------------------------------------------------------------------------
+! @brief check whether a string is enclosed with single or double quotes
+!--------------------------------------------------------------------------------------------------
+logical function quotedString(line)
+
+  character(len=*), intent(in) :: line
+
+  quotedString = .false.
+
+  if (scan(line(:1),IO_QUOTES) == 1) then
+    quotedString = .true.
+    if(line(len(line):len(line)) /= line(:1)) call IO_error(710,ext_msg=line)
+  endif
+
+end function quotedString
+
+
 !--------------------------------------------------------------------------------------------------
 ! @brief Returns Indentation.
 ! @details It determines the indentation level for a given block/line.
@@ -333,6 +357,36 @@ subroutine remove_line_break(blck,s_blck,e_char,flow_line)
 end subroutine remove_line_break
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief return the scalar list item without line break
+!--------------------------------------------------------------------------------------------------
+subroutine list_item_inline(blck,s_blck,inline)    !ToDo: SR: merge with remove_line_break eventually
+
+  character(len=*), intent(in)               :: blck                                                !< YAML in mixed style
+  integer,          intent(inout)            :: s_blck
+  character(len=:), allocatable, intent(out) :: inline
+
+  character(len=:), allocatable :: line
+  integer :: indent,indent_next
+ 
+  line   = IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))
+  indent = indentDepth(blck(s_blck:))
+  inline = line(indent+3:)
+  s_blck = s_blck + index(blck(s_blck:),IO_EOL)
+
+  indent_next = indentDepth(blck(s_blck:))
+
+  do while(indent_next > indent)
+    inline = inline//' '//trim(adjustl(IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))))
+    s_blck = s_blck + index(blck(s_blck:),IO_EOL)
+    indent_next = indentDepth(blck(s_blck:))
+  enddo
+
+  if(scan(inline,",") > 0) inline = '"'//inline//'"' 
+
+end subroutine list_item_inline
+
+
 !--------------------------------------------------------------------------------------------------
 ! @brief reads a line of YAML block which is already in flow style
 ! @details Dicts should be enlcosed within '{}' for it to be consistent with DAMASK YAML parser
@@ -463,7 +517,7 @@ recursive subroutine lst(blck,flow,s_blck,s_flow,offset)
   integer,          intent(inout) :: s_blck, &                                                      !< start position in blck
                                      s_flow, &                                                      !< start position in flow
                                      offset                                                         !< stores leading '- ' in nested lists
-  character(len=:), allocatable :: line,flow_line
+  character(len=:), allocatable :: line,flow_line,inline
   integer :: e_blck,indent
 
   indent = indentDepth(blck(s_blck:),offset)
@@ -509,8 +563,8 @@ recursive subroutine lst(blck,flow,s_blck,s_flow,offset)
       else                                                                                          ! list item in the same line
         line = line(indentDepth(line)+3:)
         if(isScalar(line)) then
-          call line_toFlow(flow,s_flow,line)
+          call list_item_inline(blck,s_blck,inline)
-          s_blck = e_blck +2
+          call line_toFlow(flow,s_flow,inline)
           offset = 0
         elseif(isFlow(line)) then
           s_blck = s_blck + index(blck(s_blck:),'-')
@@ -723,6 +777,8 @@ subroutine selfTest
   if (indentDepth('a')  /= 0)     error stop 'indentDepth'
   if (indentDepth('x ') /= 0)     error stop 'indentDepth'
 
+  if (.not. quotedString("'a'"))  error stop 'quotedString'
+
   if (      isFlow(' a'))         error stop 'isFLow'
   if (.not. isFlow('{'))          error stop 'isFlow'
   if (.not. isFlow(' ['))         error stop 'isFlow'
@@ -809,14 +865,14 @@ subroutine selfTest
 
   multi_line_flow1: block
   character(len=*), parameter :: flow_multi = &
-    "%YAML 1.1"//IO_EOL//&
+    '%YAML 1.1'//IO_EOL//&
-    "---"//IO_EOL//&
+    '---'//IO_EOL//&
-    "a: [b,"//IO_EOL//&
+    'a:     ["b",'//IO_EOL//&
-    "c: "//IO_EOL//&
+    'c: '//IO_EOL//&
-    "d, e]"//IO_EOL
+    '"d",                               "e"]'//IO_EOL
 
   character(len=*), parameter :: flow = &
-    "{a: [b, {c: d}, e]}"
+    '{a: ["b", {c: "d"}, "e"]}'
   
   if( .not. to_flow(flow_multi)        == flow) error stop 'to_flow'
   end block multi_line_flow1
@@ -848,14 +904,15 @@ subroutine selfTest
     " "//IO_EOL//&
     " "//IO_EOL//&
     "                 param_1: [a:                   b, c, {d: {e: [f: g, h]}}]"//IO_EOL//&
-    " - c: d"//IO_EOL//&
+    " - c:d"//IO_EOL//&
+    "  e.f,"//IO_EOL//&
     " bb:"//IO_EOL//&
     " "//IO_EOL//&
     "  - "//IO_EOL//&
     "   {param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}"//IO_EOL//&
     "..."//IO_EOL
   character(len=*), parameter :: mixed_flow = &
-    "{aa: [{param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}, {c: d}], bb: [{param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}]}"
+    '{aa: [{param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}, "c:d e.f,"], bb: [{param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}]}'
 
   if(.not. to_flow(block_flow) == mixed_flow)    error stop 'to_flow'
   end block basic_mixed

src/YAML_types.f90

@@ -849,7 +849,7 @@ function tNode_get_byKey_as1dFloat(self,k,defaultVal,requiredSize) result(nodeAs
 
   if (self%contains(k)) then
     node => self%get(k)
-    select type(self)
+    select type(node)
       class is(tList)
         list => node%asList()
         nodeAs1dFloat = list%as1dFloat()
@@ -872,11 +872,12 @@ end function tNode_get_byKey_as1dFloat
 !--------------------------------------------------------------------------------------------------
 !> @brief Access by key and convert to float array (2D)
 !--------------------------------------------------------------------------------------------------
-function tNode_get_byKey_as2dFloat(self,k,defaultVal) result(nodeAs2dFloat)
+function tNode_get_byKey_as2dFloat(self,k,defaultVal,requiredShape) result(nodeAs2dFloat)
 
   class(tNode),     intent(in), target                   :: self
   character(len=*), intent(in)                           :: k
   real(pReal),      intent(in), dimension(:,:), optional :: defaultVal
+  integer,          intent(in), dimension(2),   optional :: requiredShape
 
   real(pReal), dimension(:,:), allocatable :: nodeAs2dFloat
 
@@ -898,6 +899,10 @@ function tNode_get_byKey_as2dFloat(self,k,defaultVal) result(nodeAs2dFloat)
     call IO_error(143,ext_msg=k)
   endif
 
+  if (present(requiredShape)) then
+    if (any(requiredShape /= shape(nodeAs2dFloat))) call IO_error(146,ext_msg=k)
+  endif
+
 end function tNode_get_byKey_as2dFloat
 
 

src/commercialFEM_fileList.f90

@@ -7,18 +7,18 @@
 #include "IO.f90"
 #include "YAML_types.f90"
 #include "YAML_parse.f90"
+#include "HDF5_utilities.f90"
+#include "results.f90"
 #include "config.f90"
 #include "LAPACK_interface.f90"
 #include "math.f90"
 #include "rotations.f90"
+#include "lattice.f90"
 #include "element.f90"
-#include "HDF5_utilities.f90"
-#include "results.f90"
 #include "geometry_plastic_nonlocal.f90"
 #include "discretization.f90"
 #include "Marc/discretization_Marc.f90"
 #include "material.f90"
-#include "lattice.f90"
 #include "phase.f90"
 #include "phase_mechanical.f90"
 #include "phase_mechanical_elastic.f90"

src/config.f90

@@ -8,7 +8,8 @@ module config
   use IO
   use YAML_parse
   use YAML_types
-
+  use results
+  use parallelization
 
   implicit none
   private
@@ -31,6 +32,7 @@ subroutine config_init
 
   print'(/,a)', ' <<<+-  config init  -+>>>'; flush(IO_STDOUT)
 
+
   call parse_material
   call parse_numerics
   call parse_debug
@@ -41,15 +43,25 @@ end subroutine config_init
 !--------------------------------------------------------------------------------------------------
 !> @brief Read material.yaml or <jobname>.yaml.
 !--------------------------------------------------------------------------------------------------
-subroutine parse_material
+subroutine parse_material()
 
   logical :: fileExists
+  character(len=:), allocatable :: fileContent
 
 
   inquire(file='material.yaml',exist=fileExists)
   if(.not. fileExists) call IO_error(100,ext_msg='material.yaml')
-  print*, 'reading material.yaml'; flush(IO_STDOUT)
+
-  config_material => YAML_parse_file('material.yaml')
+  if (worldrank == 0) then
+    print*, 'reading material.yaml'; flush(IO_STDOUT)
+    fileContent = IO_read('material.yaml')
+    call results_openJobFile(parallel=.false.)
+    call results_writeDataset_str(fileContent,'setup','material.yaml','main configuration')
+    call results_closeJobFile
+  endif
+  call parallelization_bcast_str(fileContent)
+
+  config_material => YAML_parse_str(fileContent)
 
 end subroutine parse_material
 
@@ -57,16 +69,28 @@ end subroutine parse_material
 !--------------------------------------------------------------------------------------------------
 !> @brief Read numerics.yaml.
 !--------------------------------------------------------------------------------------------------
-subroutine parse_numerics
+subroutine parse_numerics()
 
-  logical :: fexist
+  logical :: fileExists
+  character(len=:), allocatable :: fileContent
 
 
   config_numerics => emptyDict
-  inquire(file='numerics.yaml', exist=fexist)
+
-  if (fexist) then
+  inquire(file='numerics.yaml', exist=fileExists)
-    print*, 'reading numerics.yaml'; flush(IO_STDOUT)
+  if (fileExists) then
-    config_numerics => YAML_parse_file('numerics.yaml')
+
+    if (worldrank == 0) then
+      print*, 'reading numerics.yaml'; flush(IO_STDOUT)
+      fileContent = IO_read('numerics.yaml')
+      call results_openJobFile(parallel=.false.)
+      call results_writeDataset_str(fileContent,'setup','numerics.yaml','numerics configuration')
+      call results_closeJobFile
+    endif
+    call parallelization_bcast_str(fileContent)
+
+    config_numerics => YAML_parse_str(fileContent)
+
   endif
 
 end subroutine parse_numerics
@@ -75,17 +99,29 @@ end subroutine parse_numerics
 !--------------------------------------------------------------------------------------------------
 !> @brief Read debug.yaml.
 !--------------------------------------------------------------------------------------------------
-subroutine parse_debug
+subroutine parse_debug()
 
-  logical :: fexist
+  logical :: fileExists
+  character(len=:), allocatable :: fileContent
 
 
   config_debug => emptyDict
-  inquire(file='debug.yaml', exist=fexist)
+
-  fileExists: if (fexist) then
+  inquire(file='debug.yaml', exist=fileExists)
-    print*, 'reading debug.yaml'; flush(IO_STDOUT)
+  if (fileExists) then
-    config_debug => YAML_parse_file('debug.yaml')
+
-  endif fileExists
+    if (worldrank == 0) then
+      print*, 'reading debug.yaml'; flush(IO_STDOUT)
+      fileContent = IO_read('debug.yaml')
+      call results_openJobFile(parallel=.false.)
+      call results_writeDataset_str(fileContent,'setup','debug.yaml','debug configuration')
+      call results_closeJobFile
+    endif
+    call parallelization_bcast_str(fileContent)
+
+    config_debug => YAML_parse_str(fileContent)
+
+  endif
 
 end subroutine parse_debug
 

src/element.f90

@@ -151,6 +151,7 @@ module element
   integer, dimension(NIPNEIGHBOR(CELLTYPE(1)),NIP(1)),   parameter :: IPNEIGHBOR1 = &
     reshape([&
       -2,-3,-1  &
+! Note: This fix is for gfortran 9 only. gfortran 8 supports neither, gfortran > 9 both variants
 #if !defined(__GFORTRAN__)
     ],shape(IPNEIGHBOR1))
 #else

src/grid/DAMASK_grid.f90

@@ -107,6 +107,8 @@ program DAMASK_grid
     step_bc, &
     step_mech, &
     step_discretization
+  character(len=:), allocatable :: &
+    fileContent, fname
 
 !--------------------------------------------------------------------------------------------------
 ! init DAMASK (all modules)
@@ -127,7 +129,17 @@ program DAMASK_grid
   if (stagItMax < 0)    call IO_error(301,ext_msg='maxStaggeredIter')
   if (maxCutBack < 0)   call IO_error(301,ext_msg='maxCutBack')
 
-  config_load => YAML_parse_file(trim(interface_loadFile))
+  if (worldrank == 0) then
+    fileContent = IO_read(interface_loadFile)
+    fname = interface_loadFile
+    if (scan(fname,'/') /= 0) fname = fname(scan(fname,'/',.true.)+1:)
+    call results_openJobFile(parallel=.false.)
+    call results_writeDataset_str(fileContent,'setup',fname,'load case definition (grid solver)')
+    call results_closeJobFile
+  endif
+
+  call parallelization_bcast_str(fileContent)
+  config_load => YAML_parse_str(fileContent)
   solver => config_load%get('solver')
 
 !--------------------------------------------------------------------------------------------------

src/grid/discretization_grid.f90

@@ -68,11 +68,21 @@ subroutine discretization_grid_init(restart)
     devNull, z, z_offset
   integer, dimension(worldsize) :: &
     displs, sendcounts
+  character(len=:), allocatable :: &
+    fileContent, fname
+
 
   print'(/,a)', ' <<<+-  discretization_grid init  -+>>>'; flush(IO_STDOUT)
 
+
   if(worldrank == 0) then
-    call readVTI(grid,geomSize,origin,materialAt_global)
+    fileContent = IO_read(interface_geomFile)
+    call readVTI(grid,geomSize,origin,materialAt_global,fileContent)
+    fname = interface_geomFile
+    if (scan(fname,'/') /= 0) fname = fname(scan(fname,'/',.true.)+1:)
+    call results_openJobFile(parallel=.false.)
+    call results_writeDataset_str(fileContent,'setup',fname,'geometry definition (grid solver)')
+    call results_closeJobFile
   else
     allocate(materialAt_global(0))                                                                  ! needed for IntelMPI
   endif
@@ -157,7 +167,8 @@ end subroutine discretization_grid_init
 !> @brief Parse vtk image data (.vti)
 !> @details https://vtk.org/Wiki/VTK_XML_Formats
 !--------------------------------------------------------------------------------------------------
-subroutine readVTI(grid,geomSize,origin,material)
+subroutine readVTI(grid,geomSize,origin,material, &
+                   fileContent)
 
   integer,     dimension(3), intent(out) :: &
     grid                                                                                            ! grid   (across all processes!)
@@ -166,28 +177,19 @@ subroutine readVTI(grid,geomSize,origin,material)
     origin                                                                                          ! origin (across all processes!)
   integer,     dimension(:), intent(out), allocatable :: &
     material
+  character(len=*),          intent(in) :: &
+    fileContent
 
-  character(len=:), allocatable :: fileContent, dataType, headerType
+  character(len=:), allocatable :: dataType, headerType
   logical :: inFile,inImage,gotCellData,compressed
-  integer :: fileUnit, myStat
   integer(pI64) :: &
-    fileLength, &                                                                                   !< length of the geom file (in characters)
     startPos, endPos, &
     s
 
+
   grid = -1
   geomSize = -1.0_pReal
 
-!--------------------------------------------------------------------------------------------------
-! read raw data as stream
-  inquire(file = trim(interface_geomFile), size=fileLength)
-  open(newunit=fileUnit, file=trim(interface_geomFile), access='stream',&
-       status='old', position='rewind', action='read',iostat=myStat)
-  if(myStat /= 0) call IO_error(100,ext_msg=trim(interface_geomFile))
-  allocate(character(len=fileLength)::fileContent)
-  read(fileUnit) fileContent
-  close(fileUnit)
-
   inFile         = .false.
   inImage        = .false.
   gotCelldata    = .false.

src/material.f90

@@ -8,6 +8,7 @@ module material
   use prec
   use config
   use results
+  use math
   use IO
   use rotations
   use discretization
@@ -19,8 +20,12 @@ module material
   type :: tRotationContainer
     type(Rotation), dimension(:),  allocatable :: data
   end type
+  type :: tTensorContainer
+    real(pReal), dimension(:,:,:), allocatable :: data
+  end type
 
-  type(tRotationContainer), dimension(:), allocatable :: material_orientation0
+
+  type(tRotationContainer), dimension(:), allocatable :: material_O_0
 
   integer, dimension(:), allocatable, public, protected :: &
     homogenization_Nconstituents                                                                    !< number of grains in each homogenization
@@ -41,8 +46,9 @@ module material
   integer, dimension(:,:,:), allocatable, public, protected :: &                                    ! (constituent,IP,elem)
     material_phaseMemberAt                                                                          !TODO: remove
   public :: &
+    tTensorContainer, &
     tRotationContainer, &
-    material_orientation0, &
+    material_O_0, &
     material_init
 
 contains
@@ -152,15 +158,15 @@ subroutine parse()
 
   enddo
 
-  allocate(material_orientation0(materials%length))
+  allocate(material_O_0(materials%length))
 
   do ma = 1, materials%length
     material     => materials%get(ma)
     constituents => material%get('constituents')
-    allocate(material_orientation0(ma)%data(constituents%length))
+    allocate(material_O_0(ma)%data(constituents%length))
     do co = 1, constituents%length
       constituent => constituents%get(co)
-      call material_orientation0(ma)%data(co)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4))
+      call material_O_0(ma)%data(co)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4))
     enddo
  enddo
 

src/math.f90

@@ -233,7 +233,7 @@ end function math_range
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief second rank identity tensor of specified dimension
+!> @brief Rank two identity tensor of specified dimension.
 !--------------------------------------------------------------------------------------------------
 pure function math_eye(d)
 
@@ -250,21 +250,25 @@ end function math_eye
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief symmetric fourth rank identity tensor of specified dimension
+!> @brief Symmetric rank four identity tensor.
 !  from http://en.wikipedia.org/wiki/Tensor_derivative_(continuum_mechanics)#Derivative_of_a_second-order_tensor_with_respect_to_itself
 !--------------------------------------------------------------------------------------------------
-pure function math_identity4th(d)
+pure function math_identity4th()
+
+  real(pReal), dimension(3,3,3,3) :: math_identity4th
 
-  integer, intent(in) :: d                                                                          !< tensor dimension
   integer :: i,j,k,l
-  real(pReal), dimension(d,d,d,d) :: math_identity4th
-  real(pReal), dimension(d,d)     :: identity2nd
 
-  identity2nd = math_eye(d)
+
-  do i=1,d; do j=1,d; do k=1,d; do l=1,d
+#ifndef __INTEL_COMPILER
-    math_identity4th(i,j,k,l) = 0.5_pReal &
+  do concurrent(i=1:3, j=1:3, k=1:3, l=1:3)
-                              *(identity2nd(i,k)*identity2nd(j,l)+identity2nd(i,l)*identity2nd(j,k))
+    math_identity4th(i,j,k,l) = 0.5_pReal*(math_I3(i,k)*math_I3(j,l)+math_I3(i,l)*math_I3(j,k))
+  enddo
+#else
+  do i=1,3; do j=1,3; do k=1,3; do l=1,3
+    math_identity4th(i,j,k,l) = 0.5_pReal*(math_I3(i,k)*math_I3(j,l)+math_I3(i,l)*math_I3(j,k))
   enddo; enddo; enddo; enddo
+#endif
 
 end function math_identity4th
 
@@ -331,9 +335,16 @@ pure function math_outer(A,B)
   real(pReal), dimension(size(A,1),size(B,1)) ::  math_outer
   integer :: i,j
 
+
+#ifndef __INTEL_COMPILER
+  do concurrent(i=1:size(A,1), j=1:size(B,1))
+    math_outer(i,j) = A(i)*B(j)
+  enddo
+#else
   do i=1,size(A,1); do j=1,size(B,1)
     math_outer(i,j) = A(i)*B(j)
   enddo; enddo
+#endif
 
 end function math_outer
 
@@ -373,9 +384,16 @@ pure function math_mul3333xx33(A,B)
   real(pReal), dimension(3,3) :: math_mul3333xx33
   integer :: i,j
 
+
+#ifndef __INTEL_COMPILER
+  do concurrent(i=1:3, j=1:3)
+    math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
+  enddo
+#else
   do i=1,3; do j=1,3
     math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
   enddo; enddo
+#endif
 
 end function math_mul3333xx33
 
@@ -390,9 +408,16 @@ pure function math_mul3333xx3333(A,B)
   real(pReal), dimension(3,3,3,3), intent(in) :: B
   real(pReal), dimension(3,3,3,3) :: math_mul3333xx3333
 
+
+#ifndef __INTEL_COMPILER
+  do concurrent(i=1:3, j=1:3, k=1:3, l=1:3)
+    math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l))
+  enddo
+#else
   do i=1,3; do j=1,3; do k=1,3; do l=1,3
     math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l))
   enddo; enddo; enddo; enddo
+#endif
 
 end function math_mul3333xx3333
 
@@ -545,19 +570,6 @@ pure function math_symmetric33(m)
 end function math_symmetric33
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief symmetrize a 6x6 matrix
-!--------------------------------------------------------------------------------------------------
-pure function math_symmetric66(m)
-
-  real(pReal), dimension(6,6) :: math_symmetric66
-  real(pReal), dimension(6,6), intent(in) :: m
-
-  math_symmetric66 = 0.5_pReal * (m + transpose(m))
-
-end function math_symmetric66
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief skew part of a 3x3 matrix
 !--------------------------------------------------------------------------------------------------
@@ -712,6 +724,7 @@ pure function math_6toSym33(v6,weighted)
   real(pReal), dimension(6) :: w
   integer :: i
 
+
   if(present(weighted)) then
     w = merge(INVNRMMANDEL,1.0_pReal,weighted)
   else
@@ -736,9 +749,16 @@ pure function math_3333to99(m3333)
 
   integer :: i,j
 
+
+#ifndef __INTEL_COMPILER
+  do concurrent(i=1:9, j=1:9)
+    math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j))
+  enddo
+#else
   do i=1,9; do j=1,9
     math_3333to99(i,j) = m3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j))
   enddo; enddo
+#endif
 
 end function math_3333to99
 
@@ -753,9 +773,15 @@ pure function math_99to3333(m99)
 
   integer :: i,j
 
+#ifndef __INTEL_COMPILER
+  do concurrent(i=1:9, j=1:9)
+    math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j)
+  enddo
+#else
   do i=1,9; do j=1,9
     math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j)
   enddo; enddo
+#endif
 
 end function math_99to3333
 
@@ -775,15 +801,22 @@ pure function math_sym3333to66(m3333,weighted)
   real(pReal), dimension(6) :: w
   integer :: i,j
 
+
   if(present(weighted)) then
     w = merge(NRMMANDEL,1.0_pReal,weighted)
   else
     w = NRMMANDEL
   endif
 
+#ifndef __INTEL_COMPILER
+  do concurrent(i=1:6, j=1:6)
+    math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j))
+  enddo
+#else
   do i=1,6; do j=1,6
     math_sym3333to66(i,j) = w(i)*w(j)*m3333(MAPNYE(1,i),MAPNYE(2,i),MAPNYE(1,j),MAPNYE(2,j))
   enddo; enddo
+#endif
 
 end function math_sym3333to66
 
@@ -803,6 +836,7 @@ pure function math_66toSym3333(m66,weighted)
   real(pReal), dimension(6) :: w
   integer :: i,j
 
+
   if(present(weighted)) then
     w = merge(INVNRMMANDEL,1.0_pReal,weighted)
   else
@@ -828,6 +862,7 @@ pure function math_Voigt66to3333(m66)
   real(pReal), dimension(6,6), intent(in) :: m66                                                    !< 6x6 matrix
   integer :: i,j
 
+
   do i=1,6; do j=1, 6
     math_Voigt66to3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = m66(i,j)
     math_Voigt66to3333(MAPVOIGT(2,i),MAPVOIGT(1,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) = m66(i,j)
@@ -885,6 +920,7 @@ subroutine math_eigh(w,v,error,m)
   integer :: ierr
   real(pReal), dimension(size(m,1)**2) :: work
 
+
   v = m                                                                                             ! copy matrix to input (doubles as output) array
   call dsyev('V','U',size(m,1),v,size(m,1),w,work,size(work,1),ierr)
   error = (ierr /= 0)
@@ -1222,7 +1258,7 @@ subroutine selfTest
     error stop 'math_sym33to6/math_6toSym33'
 
   call random_number(t66)
-  if(any(dNeq(math_sym3333to66(math_66toSym3333(t66)),t66))) &
+  if(any(dNeq(math_sym3333to66(math_66toSym3333(t66)),t66,1.0e-15_pReal))) &
     error stop 'math_sym3333to66/math_66toSym3333'
 
   call random_number(v6)
@@ -1242,6 +1278,9 @@ subroutine selfTest
   if(dNeq(math_det33(math_symmetric33(t33)),math_detSym33(math_symmetric33(t33)),tol=1.0e-12_pReal)) &
     error stop 'math_det33/math_detSym33'
 
+  if(any(dNeq(t33+transpose(t33),math_mul3333xx33(math_identity4th(),t33+transpose(t33))))) &
+    error stop 'math_mul3333xx33/math_identity4th'
+
   if(any(dNeq0(math_eye(3),math_inv33(math_I3)))) &
     error stop 'math_inv33(math_I3)'
 

src/mesh/DAMASK_mesh.f90

@@ -22,6 +22,15 @@ program DAMASK_mesh
 
   implicit none
 
+  type :: tLoadCase
+    real(pReal)  :: time                   = 0.0_pReal                                              !< length of increment
+    integer      :: incs                   = 0, &                                                   !< number of increments
+                    outputfrequency        = 1                                                      !< frequency of result writes
+    logical      :: followFormerTrajectory = .true.                                                 !< follow trajectory of former loadcase
+    integer,        allocatable, dimension(:) :: faceID
+    type(tFieldBC), allocatable, dimension(:) :: fieldBC
+  end type tLoadCase
+
 !--------------------------------------------------------------------------------------------------
 ! variables related to information from load case and geom file
   integer, allocatable, dimension(:) :: chunkPos                                                    ! this is longer than needed for geometry parsing
@@ -68,7 +77,7 @@ program DAMASK_mesh
 
   type(tLoadCase), allocatable, dimension(:) :: loadCases                                           !< array of all load cases
   type(tSolutionState), allocatable, dimension(:) :: solres
-  PetscInt :: faceSet, currentFaceSet, field, dimPlex
+  PetscInt :: faceSet, currentFaceSet, dimPlex
   PetscErrorCode :: ierr
   integer(kind(COMPONENT_UNDEFINED_ID)) :: ID
   external :: &
@@ -91,8 +100,7 @@ program DAMASK_mesh
 ! reading basic information from load case file and allocate data structure containing load cases
   call DMGetDimension(geomMesh,dimPlex,ierr)                                                        !< dimension of mesh (2D or 3D)
   CHKERRA(ierr)
-  nActiveFields = 1
+  allocate(solres(1))
-  allocate(solres(nActiveFields))
 
 !--------------------------------------------------------------------------------------------------
 ! reading basic information from load case file and allocate data structure containing load cases
@@ -103,8 +111,8 @@ program DAMASK_mesh
 
     chunkPos = IO_stringPos(line)
     do i = 1, chunkPos(1)                                                                           ! reading compulsory parameters for loadcase
-      select case (IO_lc(IO_stringValue(line,chunkPos,i)))
+      select case (IO_stringValue(line,chunkPos,i))
-        case('$loadcase')
+        case('$Loadcase')
           N_def = N_def + 1
       end select
     enddo                                                                                           ! count all identifiers to allocate memory and do sanity check
@@ -114,32 +122,26 @@ program DAMASK_mesh
   allocate(loadCases(N_def))
 
   do i = 1, size(loadCases)
-    allocate(loadCases(i)%fieldBC(nActiveFields))
+    allocate(loadCases(i)%fieldBC(1))
-    field = 1
+    loadCases(i)%fieldBC(1)%ID = FIELD_MECH_ID
-    loadCases(i)%fieldBC(field)%ID = FIELD_MECH_ID
   enddo
 
   do i = 1, size(loadCases)
-    do field = 1, nActiveFields
+    loadCases(i)%fieldBC(1)%nComponents = dimPlex                                                   !< X, Y (, Z) displacements
-      select case (loadCases(i)%fieldBC(field)%ID)
+    allocate(loadCases(i)%fieldBC(1)%componentBC(loadCases(i)%fieldBC(1)%nComponents))
-        case(FIELD_MECH_ID)
+    do component = 1, loadCases(i)%fieldBC(1)%nComponents
-          loadCases(i)%fieldBC(field)%nComponents = dimPlex                                         !< X, Y (, Z) displacements
+      select case (component)
-          allocate(loadCases(i)%fieldBC(field)%componentBC(loadCases(i)%fieldBC(field)%nComponents))
+        case (1)
-          do component = 1, loadCases(i)%fieldBC(field)%nComponents
+          loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_X_ID
-            select case (component)
+        case (2)
-              case (1)
+          loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_Y_ID
-                loadCases(i)%fieldBC(field)%componentBC(component)%ID = COMPONENT_MECH_X_ID
+        case (3)
-              case (2)
+          loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_Z_ID
-                loadCases(i)%fieldBC(field)%componentBC(component)%ID = COMPONENT_MECH_Y_ID
-              case (3)
-                loadCases(i)%fieldBC(field)%componentBC(component)%ID = COMPONENT_MECH_Z_ID
-            end select
-          enddo
       end select
-      do component = 1, loadCases(i)%fieldBC(field)%nComponents
+    enddo
-        allocate(loadCases(i)%fieldBC(field)%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pReal)
+    do component = 1, loadCases(i)%fieldBC(1)%nComponents
-        allocate(loadCases(i)%fieldBC(field)%componentBC(component)%Mask (mesh_Nboundaries), source = .false.)
+      allocate(loadCases(i)%fieldBC(1)%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pReal)
-      enddo
+      allocate(loadCases(i)%fieldBC(1)%componentBC(component)%Mask (mesh_Nboundaries), source = .false.)
     enddo
   enddo
 
@@ -151,52 +153,45 @@ program DAMASK_mesh
 
     chunkPos = IO_stringPos(line)
     do i = 1, chunkPos(1)
-      select case (IO_lc(IO_stringValue(line,chunkPos,i)))
+      select case (IO_stringValue(line,chunkPos,i))
 !--------------------------------------------------------------------------------------------------
 ! loadcase information
-        case('$loadcase')
+        case('$Loadcase')
           currentLoadCase = IO_intValue(line,chunkPos,i+1)
-        case('face')
+        case('Face')
           currentFace = IO_intValue(line,chunkPos,i+1)
           currentFaceSet = -1
           do faceSet = 1, mesh_Nboundaries
             if (mesh_boundaries(faceSet) == currentFace) currentFaceSet = faceSet
           enddo
           if (currentFaceSet < 0) call IO_error(error_ID = 837, ext_msg = 'invalid BC')
-        case('t','time','delta')                                                                    ! increment time
+        case('t')
           loadCases(currentLoadCase)%time = IO_floatValue(line,chunkPos,i+1)
-        case('n','incs','increments','steps')                                                       ! number of increments
+        case('N')
           loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1)
-        case('logincs','logincrements','logsteps')                                                  ! number of increments (switch to log time scaling)
+        case('f_out')
-          loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1)
-          loadCases(currentLoadCase)%logscale = 1
-        case('freq','frequency','outputfreq')                                                       ! frequency of result writings
           loadCases(currentLoadCase)%outputfrequency = IO_intValue(line,chunkPos,i+1)
-        case('guessreset','dropguessing')
+        case('estimate_rate')
           loadCases(currentLoadCase)%followFormerTrajectory = .false.                                ! do not continue to predict deformation along former trajectory
 
 !--------------------------------------------------------------------------------------------------
 ! boundary condition information
-        case('x','y','z')
+        case('X','Y','Z')
-          select case(IO_lc(IO_stringValue(line,chunkPos,i)))
+          select case(IO_stringValue(line,chunkPos,i))
-            case('x')
+            case('X')
               ID = COMPONENT_MECH_X_ID
-            case('y')
+            case('Y')
               ID = COMPONENT_MECH_Y_ID
-            case('z')
+            case('Z')
               ID = COMPONENT_MECH_Z_ID
            end select
 
-           do field = 1, nActiveFields
+           do component = 1, loadcases(currentLoadCase)%fieldBC(1)%nComponents
-             if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then
+             if (loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%ID == ID) then
-               do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
+               loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Mask (currentFaceSet) = &
-                 if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == ID) then
+                   .true.
-                   loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
+               loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Value(currentFaceSet) = &
-                       .true.
+                   IO_floatValue(line,chunkPos,i+1)
-                   loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
-                       IO_floatValue(line,chunkPos,i+1)
-                 endif
-               enddo
              endif
            enddo
       end select
@@ -212,21 +207,16 @@ program DAMASK_mesh
     print'(a,i0)', ' load case: ', currentLoadCase
     if (.not. loadCases(currentLoadCase)%followFormerTrajectory) &
       print'(a)', '  drop guessing along trajectory'
-    do field = 1, nActiveFields
+    print'(a)', '  Field '//trim(FIELD_MECH_label)
-      select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
-        case(FIELD_MECH_ID)
-          print'(a)', '  Field '//trim(FIELD_MECH_label)
 
-      end select
+    do faceSet = 1, mesh_Nboundaries
-      do faceSet = 1, mesh_Nboundaries
+       do component = 1, loadCases(currentLoadCase)%fieldBC(1)%nComponents
-         do component = 1, loadCases(currentLoadCase)%fieldBC(field)%nComponents
+         if (loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Mask(faceSet)) &
-           if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask(faceSet)) &
+           print'(a,i2,a,i2,a,f12.7)', '    Face  ', mesh_boundaries(faceSet), &
-             print'(a,i2,a,i2,a,f12.7)', '    Face  ', mesh_boundaries(faceSet), &
+                                       ' Component ', component, &
-                                               ' Component ', component, &
+                                       ' Value ', loadCases(currentLoadCase)%fieldBC(1)% &
-                                               ' Value ', loadCases(currentLoadCase)%fieldBC(field)% &
+                                                          componentBC(component)%Value(faceSet)
-                                                            componentBC(component)%Value(faceSet)
+      enddo
-         enddo
-       enddo
     enddo
     print'(a,f12.6)', '  time:       ', loadCases(currentLoadCase)%time
     if (loadCases(currentLoadCase)%incs < 1)             errorID = 835                            ! non-positive incs count
@@ -240,12 +230,7 @@ program DAMASK_mesh
 !--------------------------------------------------------------------------------------------------
 ! doing initialization depending on active solvers
   call FEM_Utilities_init
-  do field = 1, nActiveFields
+  call FEM_mechanical_init(loadCases(1)%fieldBC(1))
-    select case (loadCases(1)%fieldBC(field)%ID)
-      case(FIELD_MECH_ID)
-        call FEM_mechanical_init(loadCases(1)%fieldBC(field))
-    end select
-  enddo
 
   if (worldrank == 0) then
     open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE')
@@ -266,32 +251,14 @@ program DAMASK_mesh
 !--------------------------------------------------------------------------------------------------
 ! forwarding time
       timeIncOld = timeinc                                                                          ! last timeinc that brought former inc to an end
-      if (loadCases(currentLoadCase)%logscale == 0) then                                            ! linear scale
+      timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal)
-        timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal)
-      else
-        if (currentLoadCase == 1) then                                                              ! 1st load case of logarithmic scale
-          if (inc == 1) then                                                                        ! 1st inc of 1st load case of logarithmic scale
-            timeinc = loadCases(1)%time*(2.0_pReal**real(    1-loadCases(1)%incs ,pReal))           ! assume 1st inc is equal to 2nd
-          else                                                                                      ! not-1st inc of 1st load case of logarithmic scale
-            timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1-loadCases(1)%incs ,pReal))
-          endif
-        else                                                                                        ! not-1st load case of logarithmic scale
-          timeinc = time0 * &
-               ( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real(          inc,pReal)/&
-                                                     real(loadCases(currentLoadCase)%incs ,pReal))&
-                -(1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc-1  ,pReal)/&
-                                                      real(loadCases(currentLoadCase)%incs ,pReal)))
-        endif
-      endif
       timeinc = timeinc * real(subStepFactor,pReal)**real(-cutBackLevel,pReal)                      ! depending on cut back level, decrease time step
-
+      stepFraction = 0                                                                              ! fraction scaled by stepFactor**cutLevel
-
-      stepFraction = 0                                                                            ! fraction scaled by stepFactor**cutLevel
 
       subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel)
         remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time
-        time = time + timeinc                                                                     ! forward target time
+        time = time + timeinc                                                                       ! forward target time
-        stepFraction = stepFraction + 1                                                           ! count step
+        stepFraction = stepFraction + 1                                                             ! count step
 
 !--------------------------------------------------------------------------------------------------
 ! report begin of new step
@@ -306,33 +273,16 @@ program DAMASK_mesh
                '-',stepFraction, '/', subStepFactor**cutBackLevel
         flush(IO_STDOUT)
 
-!--------------------------------------------------------------------------------------------------
+        call FEM_mechanical_forward(guess,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(1))
-! forward fields
-        do field = 1, nActiveFields
-          select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
-            case(FIELD_MECH_ID)
-              call FEM_mechanical_forward (&
-                  guess,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field))
-
-         end select
-        enddo
 
 !--------------------------------------------------------------------------------------------------
 ! solve fields
         stagIter = 0
         stagIterate = .true.
         do while (stagIterate)
-          do field = 1, nActiveFields
+          solres(1) = FEM_mechanical_solution(incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(1))
-            select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
+          if(.not. solres(1)%converged) exit
-              case(FIELD_MECH_ID)
-                solres(field) = FEM_mechanical_solution (&
-                      incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field))
 
-            end select
-
-            if(.not. solres(field)%converged) exit                                                  ! no solution found
-
-          enddo
           stagIter = stagIter + 1
           stagIterate =            stagIter < stagItMax &
                        .and.       all(solres(:)%converged) &

src/mesh/FEM_quadrature.f90

@@ -9,7 +9,7 @@ module FEM_quadrature
   private
 
   integer, parameter :: &
-    maxOrder = 5                                                                                    !< current max interpolation set at cubic (intended to be arbitrary)
+    maxOrder = 5                                                                                    !< maximum integration order
   real(pReal),   dimension(2,3),       parameter :: &
     triangle    = reshape([-1.0_pReal, -1.0_pReal, &
                             1.0_pReal, -1.0_pReal, &
@@ -20,8 +20,12 @@ module FEM_quadrature
                            -1.0_pReal,  1.0_pReal, -1.0_pReal, &
                            -1.0_pReal, -1.0_pReal,  1.0_pReal], shape=[3,4])
 
+  type :: group_float                                                                               !< variable length datatype
+    real(pReal), dimension(:), allocatable :: p
+  end type group_float
+
   integer,             dimension(2:3,maxOrder), public, protected :: &
-    FEM_nQuadrature                                                                                 !< number of quadrature points for a given spatial dimension(2-3) and interpolation order(1-maxOrder)
+    FEM_nQuadrature                                                                                 !< number of quadrature points for spatial dimension(2-3) and interpolation order (1-maxOrder)
   type(group_float),   dimension(2:3,maxOrder), public, protected :: &
     FEM_quadrature_weights, &                                                                       !< quadrature weights for each quadrature rule
     FEM_quadrature_points                                                                           !< quadrature point coordinates (in simplical system) for each quadrature rule
@@ -35,145 +39,146 @@ contains
 !--------------------------------------------------------------------------------------------------
 !> @brief initializes FEM interpolation data
 !--------------------------------------------------------------------------------------------------
-subroutine FEM_quadrature_init
+subroutine FEM_quadrature_init()
 
   print'(/,a)', ' <<<+-  FEM_quadrature init  -+>>>'; flush(6)
 
+  print*, 'L. Zhang et al., Journal of Computational Mathematics 27(1):89-96, 2009'
+  print*, 'https://www.jstor.org/stable/43693493'
+
 !--------------------------------------------------------------------------------------------------
 ! 2D linear
   FEM_nQuadrature(2,1) = 1
 
-  allocate(FEM_quadrature_weights(2,1)%p(1))
+  allocate(FEM_quadrature_weights(2,1)%p(FEM_nQuadrature(2,1)))
-  FEM_quadrature_weights(2,1)%p(1)   = 1.0_pReal
+  FEM_quadrature_weights(2,1)%p(1) = 1._pReal
 
-  allocate(FEM_quadrature_points (2,1)%p(2))
+  FEM_quadrature_points (2,1)%p = permutationStar3([1._pReal/3._pReal])
-  FEM_quadrature_points (2,1)%p(1:2) = permutationStar3([1.0_pReal/3.0_pReal])
 
 !--------------------------------------------------------------------------------------------------
 ! 2D quadratic
   FEM_nQuadrature(2,2) = 3
 
-  allocate(FEM_quadrature_weights(2,2)%p(3))
+  allocate(FEM_quadrature_weights(2,2)%p(FEM_nQuadrature(2,2)))
-  FEM_quadrature_weights(2,2)%p(1:3) = 1.0_pReal/3.0_pReal
+  FEM_quadrature_weights(2,2)%p(1:3) = 1._pReal/3._pReal
 
-  allocate(FEM_quadrature_points (2,2)%p(6))
+  FEM_quadrature_points (2,2)%p = permutationStar21([1._pReal/6._pReal])
-  FEM_quadrature_points (2,2)%p(1:6) = permutationStar21([1.0_pReal/6.0_pReal])
 
 !--------------------------------------------------------------------------------------------------
 ! 2D cubic
   FEM_nQuadrature(2,3) = 6
 
-  allocate(FEM_quadrature_weights(2,3)%p(6))
+  allocate(FEM_quadrature_weights(2,3)%p(FEM_nQuadrature(2,3)))
-  FEM_quadrature_weights(2,3)%p(1:3) = 0.22338158967801146570_pReal
+  FEM_quadrature_weights(2,3)%p(1:3) = 2.2338158967801147e-1_pReal
-  FEM_quadrature_weights(2,3)%p(4:6) = 0.10995174365532186764_pReal
+  FEM_quadrature_weights(2,3)%p(4:6) = 1.0995174365532187e-1_pReal
 
-  allocate(FEM_quadrature_points (2,3)%p(12))
+  FEM_quadrature_points (2,3)%p = [ &
-  FEM_quadrature_points (2,3)%p(1:6) = permutationStar21([0.44594849091596488632_pReal])
+    permutationStar21([4.4594849091596489e-1_pReal]), &
-  FEM_quadrature_points (2,3)%p(7:12)= permutationStar21([0.091576213509770743460_pReal])
+    permutationStar21([9.157621350977074e-2_pReal]) ]
 
 !--------------------------------------------------------------------------------------------------
 ! 2D quartic
   FEM_nQuadrature(2,4) = 12
 
-  allocate(FEM_quadrature_weights(2,4)%p(12))
+  allocate(FEM_quadrature_weights(2,4)%p(FEM_nQuadrature(2,4)))
-  FEM_quadrature_weights(2,4)%p(1:3)  = 0.11678627572638_pReal
+  FEM_quadrature_weights(2,4)%p(1:3)  = 1.1678627572637937e-1_pReal
-  FEM_quadrature_weights(2,4)%p(4:6)  = 0.05084490637021_pReal
+  FEM_quadrature_weights(2,4)%p(4:6)  = 5.0844906370206817e-2_pReal
-  FEM_quadrature_weights(2,4)%p(7:12) = 0.08285107561837_pReal
+  FEM_quadrature_weights(2,4)%p(7:12) = 8.285107561837358e-2_pReal
 
-  allocate(FEM_quadrature_points (2,4)%p(24))
+  FEM_quadrature_points (2,4)%p = [ &
-  FEM_quadrature_points (2,4)%p(1:6)  = permutationStar21([0.24928674517091_pReal])
+    permutationStar21([2.4928674517091042e-1_pReal]), &
-  FEM_quadrature_points (2,4)%p(7:12) = permutationStar21([0.06308901449150_pReal])
+    permutationStar21([6.308901449150223e-2_pReal]), &
-  FEM_quadrature_points (2,4)%p(13:24)= permutationStar111([0.31035245103378_pReal, 0.63650249912140_pReal])
+    permutationStar111([3.1035245103378440e-1_pReal, 5.3145049844816947e-2_pReal]) ]
 
 !--------------------------------------------------------------------------------------------------
 ! 2D quintic
   FEM_nQuadrature(2,5) = 16
 
-  allocate(FEM_quadrature_weights(2,5)%p(16))
+  allocate(FEM_quadrature_weights(2,5)%p(FEM_nQuadrature(2,5)))
-  FEM_quadrature_weights(2,5)%p(1  )  = 0.14431560767779_pReal
+  FEM_quadrature_weights(2,5)%p(1:1)   = 1.4431560767778717e-1_pReal
-  FEM_quadrature_weights(2,5)%p(2:4)  = 0.09509163426728_pReal
+  FEM_quadrature_weights(2,5)%p(2:4)   = 9.509163426728463e-2_pReal
-  FEM_quadrature_weights(2,5)%p(5:7)  = 0.10321737053472_pReal
+  FEM_quadrature_weights(2,5)%p(5:7)   = 1.0321737053471825e-1_pReal
-  FEM_quadrature_weights(2,5)%p(8:10) = 0.03245849762320_pReal
+  FEM_quadrature_weights(2,5)%p(8:10)  = 3.2458497623198080e-2_pReal
-  FEM_quadrature_weights(2,5)%p(11:16)= 0.02723031417443_pReal
+  FEM_quadrature_weights(2,5)%p(11:16) = 2.7230314174434994e-2_pReal
 
-  allocate(FEM_quadrature_points (2,5)%p(32))
+  FEM_quadrature_points (2,5)%p = [ &
-  FEM_quadrature_points (2,5)%p(1:2)  = permutationStar3([0.33333333333333_pReal])
+    permutationStar3([1._pReal/3._pReal]), &
-  FEM_quadrature_points (2,5)%p(3:8)  = permutationStar21([0.45929258829272_pReal])
+    permutationStar21([4.5929258829272316e-1_pReal]), &
-  FEM_quadrature_points (2,5)%p(9:14) = permutationStar21([0.17056930775176_pReal])
+    permutationStar21([1.705693077517602e-1_pReal]), &
-  FEM_quadrature_points (2,5)%p(15:20)= permutationStar21([0.05054722831703_pReal])
+    permutationStar21([5.0547228317030975e-2_pReal]), &
-  FEM_quadrature_points (2,5)%p(21:32)= permutationStar111([0.26311282963464_pReal, 0.72849239295540_pReal])
+    permutationStar111([2.631128296346381e-1_pReal, 8.3947774099576053e-2_pReal]) ]
 
 !--------------------------------------------------------------------------------------------------
 ! 3D linear
   FEM_nQuadrature(3,1) = 1
 
-  allocate(FEM_quadrature_weights(3,1)%p(1))
+  allocate(FEM_quadrature_weights(3,1)%p(FEM_nQuadrature(3,1)))
-  FEM_quadrature_weights(3,1)%p(1)  = 1.0_pReal
+  FEM_quadrature_weights(3,1)%p(1) = 1.0_pReal
 
-  allocate(FEM_quadrature_points (3,1)%p(3))
+  FEM_quadrature_points (3,1)%p = permutationStar4([0.25_pReal])
-  FEM_quadrature_points (3,1)%p(1:3)= permutationStar4([0.25_pReal])
 
 !--------------------------------------------------------------------------------------------------
 ! 3D quadratic
   FEM_nQuadrature(3,2) = 4
 
-  allocate(FEM_quadrature_weights(3,2)%p(4))
+  allocate(FEM_quadrature_weights(3,2)%p(FEM_nQuadrature(3,2)))
   FEM_quadrature_weights(3,2)%p(1:4) = 0.25_pReal
 
-  allocate(FEM_quadrature_points (3,2)%p(12))
+  FEM_quadrature_points (3,2)%p = permutationStar31([1.3819660112501052e-1_pReal])
-  FEM_quadrature_points (3,2)%p(1:12)= permutationStar31([0.13819660112501051518_pReal])
 
 !--------------------------------------------------------------------------------------------------
 ! 3D cubic
   FEM_nQuadrature(3,3) = 14
 
-  allocate(FEM_quadrature_weights(3,3)%p(14))
+  allocate(FEM_quadrature_weights(3,3)%p(FEM_nQuadrature(3,3)))
-  FEM_quadrature_weights(3,3)%p(5:8)  = 0.11268792571801585080_pReal
+  FEM_quadrature_weights(3,3)%p(1:4)  = 7.3493043116361949e-2_pReal
-  FEM_quadrature_weights(3,3)%p(1:4)  = 0.073493043116361949544_pReal
+  FEM_quadrature_weights(3,3)%p(5:8)  = 1.1268792571801585e-1_pReal
-  FEM_quadrature_weights(3,3)%p(9:14) = 0.042546020777081466438_pReal
+  FEM_quadrature_weights(3,3)%p(9:14) = 4.2546020777081467e-2_pReal
 
-  allocate(FEM_quadrature_points (3,3)%p(42))
+  FEM_quadrature_points (3,3)%p = [ &
-  FEM_quadrature_points (3,3)%p(1:12) = permutationStar31([0.092735250310891226402_pReal])
+    permutationStar31([9.273525031089123e-2_pReal]), &
-  FEM_quadrature_points (3,3)%p(13:24)= permutationStar31([0.31088591926330060980_pReal])
+    permutationStar31([3.108859192633006e-1_pReal]), &
-  FEM_quadrature_points (3,3)%p(25:42)= permutationStar22([0.045503704125649649492_pReal])
+    permutationStar22([4.5503704125649649e-2_pReal]) ]
 
 !--------------------------------------------------------------------------------------------------
-! 3D quartic
+! 3D quartic (lower precision/unknown source)
   FEM_nQuadrature(3,4) = 35
 
-  allocate(FEM_quadrature_weights(3,4)%p(35))
+  allocate(FEM_quadrature_weights(3,4)%p(FEM_nQuadrature(3,4)))
-  FEM_quadrature_weights(3,4)%p(1:4)    = 0.0021900463965388_pReal
+  FEM_quadrature_weights(3,4)%p(1:4)   = 0.0021900463965388_pReal
-  FEM_quadrature_weights(3,4)%p(5:16)   = 0.0143395670177665_pReal
+  FEM_quadrature_weights(3,4)%p(5:16)  = 0.0143395670177665_pReal
-  FEM_quadrature_weights(3,4)%p(17:22)  = 0.0250305395686746_pReal
+  FEM_quadrature_weights(3,4)%p(17:22) = 0.0250305395686746_pReal
-  FEM_quadrature_weights(3,4)%p(23:34)  = 0.0479839333057554_pReal
+  FEM_quadrature_weights(3,4)%p(23:34) = 0.0479839333057554_pReal
-  FEM_quadrature_weights(3,4)%p(35)     = 0.0931745731195340_pReal
+  FEM_quadrature_weights(3,4)%p(35)    = 0.0931745731195340_pReal
 
-  allocate(FEM_quadrature_points (3,4)%p(105))
+  FEM_quadrature_points (3,4)%p = [ &
-  FEM_quadrature_points (3,4)%p(1:12)   = permutationStar31([0.0267367755543735_pReal])
+    permutationStar31([0.0267367755543735_pReal]), &
-  FEM_quadrature_points (3,4)%p(13:48)  = permutationStar211([0.0391022406356488_pReal, 0.7477598884818090_pReal])
+    permutationStar211([0.0391022406356488_pReal, 0.7477598884818090_pReal]), &
-  FEM_quadrature_points (3,4)%p(49:66)  = permutationStar22([0.4547545999844830_pReal])
+    permutationStar22([0.4547545999844830_pReal]), &
-  FEM_quadrature_points (3,4)%p(67:102) = permutationStar211([0.2232010379623150_pReal, 0.0504792790607720_pReal])
+    permutationStar211([0.2232010379623150_pReal, 0.0504792790607720_pReal]), &
-  FEM_quadrature_points (3,4)%p(103:105)= permutationStar4([0.25_pReal])
+    permutationStar4([0.25_pReal]) ]
 
 !--------------------------------------------------------------------------------------------------
-! 3D quintic
+! 3D quintic (lower precision/unknown source)
   FEM_nQuadrature(3,5) = 56
 
-  allocate(FEM_quadrature_weights(3,5)%p(56))
+  allocate(FEM_quadrature_weights(3,5)%p(FEM_nQuadrature(3,5)))
-  FEM_quadrature_weights(3,5)%p(1:4)    = 0.0010373112336140_pReal
+  FEM_quadrature_weights(3,5)%p(1:4)   = 0.0010373112336140_pReal
-  FEM_quadrature_weights(3,5)%p(5:16)   = 0.0096016645399480_pReal
+  FEM_quadrature_weights(3,5)%p(5:16)  = 0.0096016645399480_pReal
-  FEM_quadrature_weights(3,5)%p(17:28)  = 0.0164493976798232_pReal
+  FEM_quadrature_weights(3,5)%p(17:28) = 0.0164493976798232_pReal
-  FEM_quadrature_weights(3,5)%p(29:40)  = 0.0153747766513310_pReal
+  FEM_quadrature_weights(3,5)%p(29:40) = 0.0153747766513310_pReal
-  FEM_quadrature_weights(3,5)%p(41:52)  = 0.0293520118375230_pReal
+  FEM_quadrature_weights(3,5)%p(41:52) = 0.0293520118375230_pReal
-  FEM_quadrature_weights(3,5)%p(53:56)  = 0.0366291366405108_pReal
+  FEM_quadrature_weights(3,5)%p(53:56) = 0.0366291366405108_pReal
 
-  allocate(FEM_quadrature_points (3,5)%p(168))
+  FEM_quadrature_points (3,5)%p = [ &
-  FEM_quadrature_points (3,5)%p(1:12)   = permutationStar31([0.0149520651530592_pReal])
+    permutationStar31([0.0149520651530592_pReal]), &
-  FEM_quadrature_points (3,5)%p(13:48)  = permutationStar211([0.0340960211962615_pReal, 0.1518319491659370_pReal])
+    permutationStar211([0.0340960211962615_pReal, 0.1518319491659370_pReal]), &
-  FEM_quadrature_points (3,5)%p(49:84)  = permutationStar211([0.0462051504150017_pReal, 0.3549340560639790_pReal])
+    permutationStar211([0.0462051504150017_pReal, 0.3549340560639790_pReal]), &
-  FEM_quadrature_points (3,5)%p(85:120) = permutationStar211([0.2281904610687610_pReal, 0.0055147549744775_pReal])
+    permutationStar211([0.2281904610687610_pReal, 0.0055147549744775_pReal]), &
-  FEM_quadrature_points (3,5)%p(121:156)= permutationStar211([0.3523052600879940_pReal, 0.0992057202494530_pReal])
+    permutationStar211([0.3523052600879940_pReal, 0.0992057202494530_pReal]), &
-  FEM_quadrature_points (3,5)%p(157:168)= permutationStar31([0.1344783347929940_pReal])
+    permutationStar31([0.1344783347929940_pReal]) ]
+
+  call selfTest
 
 end subroutine FEM_quadrature_init
 
@@ -186,11 +191,9 @@ pure function permutationStar3(point) result(qPt)
   real(pReal), dimension(2)             :: qPt
   real(pReal), dimension(1), intent(in) :: point
 
-  real(pReal), dimension(3,1) :: temp
 
-  temp(:,1) = [point(1), point(1), point(1)]
+  qPt = pack(matmul(triangle,reshape([ &
-
+    point(1), point(1), point(1)],[3,1])),.true.)
-  qPt = reshape(matmul(triangle, temp),[2])
 
 end function permutationStar3
 
@@ -203,13 +206,11 @@ pure function permutationStar21(point) result(qPt)
   real(pReal), dimension(6)             :: qPt
   real(pReal), dimension(1), intent(in) :: point
 
-  real(pReal), dimension(3,3) :: temp
 
-  temp(:,1) = [point(1), point(1), 1.0_pReal - 2.0_pReal*point(1)]
+  qPt = pack(matmul(triangle,reshape([ &
-  temp(:,2) = [point(1), 1.0_pReal - 2.0_pReal*point(1), point(1)]
+    point(1), point(1), 1.0_pReal - 2.0_pReal*point(1), &
-  temp(:,3) = [1.0_pReal - 2.0_pReal*point(1), point(1), point(1)]
+    point(1), 1.0_pReal - 2.0_pReal*point(1), point(1), &
-
+    1.0_pReal - 2.0_pReal*point(1), point(1), point(1)],[3,3])),.true.)
- qPt = reshape(matmul(triangle, temp),[6])
 
 end function permutationStar21
 
@@ -222,16 +223,14 @@ pure function permutationStar111(point) result(qPt)
   real(pReal), dimension(12)            :: qPt
   real(pReal), dimension(2), intent(in) :: point
 
-  real(pReal), dimension(3,6) :: temp
 
-  temp(:,1) = [point(1), point(2), 1.0_pReal - point(1) - point(2)]
+  qPt = pack(matmul(triangle,reshape([ &
-  temp(:,2) = [point(1), 1.0_pReal - point(1) - point(2), point(2)]
+    point(1), point(2), 1.0_pReal - point(1) - point(2), &
-  temp(:,3) = [point(2), point(1), 1.0_pReal - point(1) - point(2)]
+    point(1), 1.0_pReal - point(1) - point(2), point(2), &
-  temp(:,4) = [point(2), 1.0_pReal - point(1) - point(2), point(1)]
+    point(2), point(1), 1.0_pReal - point(1) - point(2), &
-  temp(:,5) = [1.0_pReal - point(1) - point(2), point(2), point(1)]
+    point(2), 1.0_pReal - point(1) - point(2), point(1), &
-  temp(:,6) = [1.0_pReal - point(1) - point(2), point(1), point(2)]
+    1.0_pReal - point(1) - point(2), point(2), point(1), &
-
+    1.0_pReal - point(1) - point(2), point(1), point(2)],[3,6])),.true.)
-  qPt = reshape(matmul(triangle, temp),[12])
 
 end function permutationStar111
 
@@ -244,11 +243,9 @@ pure function permutationStar4(point) result(qPt)
   real(pReal), dimension(3)             :: qPt
   real(pReal), dimension(1), intent(in) :: point
 
-  real(pReal), dimension(4,1) :: temp
 
-  temp(:,1) = [point(1), point(1), point(1), point(1)]
+  qPt = pack(matmul(tetrahedron,reshape([ &
-
+    point(1), point(1), point(1), point(1)],[4,1])),.true.)
-  qPt = reshape(matmul(tetrahedron, temp),[3])
 
 end function permutationStar4
 
@@ -261,14 +258,12 @@ pure function permutationStar31(point) result(qPt)
   real(pReal), dimension(12)            :: qPt
   real(pReal), dimension(1), intent(in) :: point
 
-  real(pReal), dimension(4,4) :: temp
 
-  temp(:,1) = [point(1), point(1), point(1), 1.0_pReal - 3.0_pReal*point(1)]
+  qPt = pack(matmul(tetrahedron,reshape([ &
-  temp(:,2) = [point(1), point(1), 1.0_pReal - 3.0_pReal*point(1), point(1)]
+    point(1), point(1), point(1), 1.0_pReal - 3.0_pReal*point(1), &
-  temp(:,3) = [point(1), 1.0_pReal - 3.0_pReal*point(1), point(1), point(1)]
+    point(1), point(1), 1.0_pReal - 3.0_pReal*point(1), point(1), &
-  temp(:,4) = [1.0_pReal - 3.0_pReal*point(1), point(1), point(1), point(1)]
+    point(1), 1.0_pReal - 3.0_pReal*point(1), point(1), point(1), &
-
+    1.0_pReal - 3.0_pReal*point(1), point(1), point(1), point(1)],[4,4])),.true.)
-  qPt = reshape(matmul(tetrahedron, temp),[12])
 
 end function permutationStar31
 
@@ -276,21 +271,19 @@ end function permutationStar31
 !--------------------------------------------------------------------------------------------------
 !> @brief star 22 permutation of input
 !--------------------------------------------------------------------------------------------------
-pure function permutationStar22(point) result(qPt)
+function permutationStar22(point) result(qPt)
 
   real(pReal), dimension(18)            :: qPt
   real(pReal), dimension(1), intent(in) :: point
 
-  real(pReal), dimension(4,6) :: temp
 
-  temp(:,1) = [point(1), point(1), 0.5_pReal - point(1), 0.5_pReal - point(1)]
+  qPt = pack(matmul(tetrahedron,reshape([ &
-  temp(:,2) = [point(1), 0.5_pReal - point(1), point(1), 0.5_pReal - point(1)]
+    point(1), point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), &
-  temp(:,3) = [0.5_pReal - point(1), point(1), point(1), 0.5_pReal - point(1)]
+    point(1), 0.5_pReal - point(1), point(1), 0.5_pReal - point(1), &
-  temp(:,4) = [0.5_pReal - point(1), point(1), 0.5_pReal - point(1), point(1)]
+    0.5_pReal - point(1), point(1), point(1), 0.5_pReal - point(1), &
-  temp(:,5) = [0.5_pReal - point(1), 0.5_pReal - point(1), point(1), point(1)]
+    0.5_pReal - point(1), point(1), 0.5_pReal - point(1), point(1), &
-  temp(:,6) = [point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), point(1)]
+    0.5_pReal - point(1), 0.5_pReal - point(1), point(1), point(1), &
-
+    point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), point(1)],[4,6])),.true.)
-  qPt = reshape(matmul(tetrahedron, temp),[18])
 
 end function permutationStar22
 
@@ -303,22 +296,20 @@ pure function permutationStar211(point) result(qPt)
   real(pReal), dimension(36)            :: qPt
   real(pReal), dimension(2), intent(in) :: point
 
-  real(pReal), dimension(4,12) :: temp
 
-  temp(:,1 ) = [point(1), point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2)]
+  qPt = pack(matmul(tetrahedron,reshape([ &
-  temp(:,2 ) = [point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2)]
+    point(1), point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), &
-  temp(:,3 ) = [point(1), point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2)]
+    point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), &
-  temp(:,4 ) = [point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1)]
+    point(1), point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), &
-  temp(:,5 ) = [point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2)]
+    point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), &
-  temp(:,6 ) = [point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1)]
+    point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2), &
-  temp(:,7 ) = [point(2), point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2)]
+    point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1), &
-  temp(:,8 ) = [point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1)]
+    point(2), point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), &
-  temp(:,9 ) = [point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1)]
+    point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), &
-  temp(:,10) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1), point(2)]
+    point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1), &
-  temp(:,11) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2), point(1)]
+    1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1), point(2), &
-  temp(:,12) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1), point(1)]
+    1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2), point(1), &
-
+    1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1), point(1)],[4,12])),.true.)
-  qPt = reshape(matmul(tetrahedron, temp),[36])
 
 end function permutationStar211
 
@@ -331,35 +322,60 @@ pure function permutationStar1111(point) result(qPt)
   real(pReal), dimension(72)            :: qPt
   real(pReal), dimension(3), intent(in) :: point
 
-  real(pReal), dimension(4,24) :: temp
 
-  temp(:,1 ) = [point(1), point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3)]
+  qPt = pack(matmul(tetrahedron,reshape([ &
-  temp(:,2 ) = [point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3)]
+    point(1), point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3), &
-  temp(:,3 ) = [point(1), point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3)]
+    point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3), &
-  temp(:,4 ) = [point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2)]
+    point(1), point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3), &
-  temp(:,5 ) = [point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(3)]
+    point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2), &
-  temp(:,6 ) = [point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(2)]
+    point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(3), &
-  temp(:,7 ) = [point(2), point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3)]
+    point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(2), &
-  temp(:,8 ) = [point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3)]
+    point(2), point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3), &
-  temp(:,9 ) = [point(2), point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3)]
+    point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3), &
-  temp(:,10) = [point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1)]
+    point(2), point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3), &
-  temp(:,11) = [point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(3)]
+    point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1), &
-  temp(:,12) = [point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(1)]
+    point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(3), &
-  temp(:,13) = [point(3), point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3)]
+    point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(1), &
-  temp(:,14) = [point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2)]
+    point(3), point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3), &
-  temp(:,15) = [point(3), point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3)]
+    point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2), &
-  temp(:,16) = [point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1)]
+    point(3), point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3), &
-  temp(:,17) = [point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(2)]
+    point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1), &
-  temp(:,18) = [point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(1)]
+    point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(2), &
-  temp(:,19) = [1.0_pReal - point(1) - point(2)- point(3), point(1), point(2), point(3)]
+    point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(1), &
-  temp(:,20) = [1.0_pReal - point(1) - point(2)- point(3), point(1), point(3), point(2)]
+    1.0_pReal - point(1) - point(2)- point(3), point(1), point(2), point(3), &
-  temp(:,21) = [1.0_pReal - point(1) - point(2)- point(3), point(2), point(1), point(3)]
+    1.0_pReal - point(1) - point(2)- point(3), point(1), point(3), point(2), &
-  temp(:,22) = [1.0_pReal - point(1) - point(2)- point(3), point(2), point(3), point(1)]
+    1.0_pReal - point(1) - point(2)- point(3), point(2), point(1), point(3), &
-  temp(:,23) = [1.0_pReal - point(1) - point(2)- point(3), point(3), point(1), point(2)]
+    1.0_pReal - point(1) - point(2)- point(3), point(2), point(3), point(1), &
-  temp(:,24) = [1.0_pReal - point(1) - point(2)- point(3), point(3), point(2), point(1)]
+    1.0_pReal - point(1) - point(2)- point(3), point(3), point(1), point(2), &
-
+    1.0_pReal - point(1) - point(2)- point(3), point(3), point(2), point(1)],[4,24])),.true.)
-  qPt = reshape(matmul(tetrahedron, temp),[72])
 
 end function permutationStar1111
 
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Check correctness of quadrature weights and points.
+!--------------------------------------------------------------------------------------------------
+subroutine selfTest
+
+  integer :: o, d, n
+  real(pReal), dimension(2:3), parameter :: w = [3.0_pReal,2.0_pReal]
+
+
+  do d = lbound(FEM_quadrature_weights,1), ubound(FEM_quadrature_weights,1)
+    do o = lbound(FEM_quadrature_weights(d,:),1), ubound(FEM_quadrature_weights(d,:),1)
+      if (dNeq(sum(FEM_quadrature_weights(d,o)%p),1.0_pReal,5e-15_pReal)) &
+        error stop 'quadrature weights'
+    enddo
+  enddo
+
+  do d = lbound(FEM_quadrature_points,1), ubound(FEM_quadrature_points,1)
+    do o = lbound(FEM_quadrature_points(d,:),1), ubound(FEM_quadrature_points(d,:),1)
+      n = size(FEM_quadrature_points(d,o)%p,1)/d
+      if (any(dNeq(sum(reshape(FEM_quadrature_points(d,o)%p,[d,n]),2),-real(n,pReal)/w(d),1.e-14_pReal))) &
+        error stop 'quadrature points'
+    enddo
+  enddo
+
+end subroutine selfTest
+
 end module FEM_quadrature

src/mesh/FEM_utilities.f90

@@ -23,14 +23,8 @@ module FEM_utilities
   implicit none
   private
 
-!--------------------------------------------------------------------------------------------------
+  logical,     public             :: cutBack = .false.                                              !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
-  logical, public             :: cutBack = .false.                                                  !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
+  real(pReal), public, protected  :: wgt                                                            !< weighting factor 1/Nelems
-  integer, public, parameter  :: maxFields = 6
-  integer, public             :: nActiveFields = 0
-
-!--------------------------------------------------------------------------------------------------
-! grid related information information
-  real(pReal),   public       :: wgt                                                                !< weighting factor 1/Nelems
 
 
 !--------------------------------------------------------------------------------------------------
@@ -49,10 +43,6 @@ module FEM_utilities
     COMPONENT_MECH_Z_ID
   end enum
 
-!--------------------------------------------------------------------------------------------------
-! variables controlling debugging
- logical :: &
-   debugPETSc                                                                                       !< use some in debug defined options for more verbose PETSc solution
 
 !--------------------------------------------------------------------------------------------------
 ! derived types
@@ -63,27 +53,17 @@ module FEM_utilities
   end type tSolutionState
 
   type, public :: tComponentBC
-    integer(kind(COMPONENT_UNDEFINED_ID))          :: ID
+    integer(kind(COMPONENT_UNDEFINED_ID)) :: ID
-    real(pReal),                       allocatable, dimension(:) :: Value
+    real(pReal), allocatable, dimension(:) :: Value
-    logical,                           allocatable, dimension(:) :: Mask
+    logical,     allocatable, dimension(:) :: Mask
   end type tComponentBC
 
   type, public :: tFieldBC
     integer(kind(FIELD_UNDEFINED_ID))  :: ID
     integer                            :: nComponents = 0
-    type(tComponentBC),    allocatable :: componentBC(:)
+    type(tComponentBC), allocatable, dimension(:) :: componentBC
   end type tFieldBC
 
-  type, public :: tLoadCase
-    real(pReal)  :: time                   = 0.0_pReal                                              !< length of increment
-    integer      :: incs                   = 0, &                                                   !< number of increments
-                    outputfrequency        = 1, &                                                   !< frequency of result writes
-                    logscale               = 0                                                      !< linear/logarithmic time inc flag
-    logical      :: followFormerTrajectory = .true.                                                 !< follow trajectory of former loadcase
-    integer,        allocatable, dimension(:) :: faceID
-    type(tFieldBC), allocatable, dimension(:) :: fieldBC
-  end type tLoadCase
-
   public :: &
     FEM_utilities_init, &
     utilities_constitutiveResponse, &
@@ -109,8 +89,9 @@ subroutine FEM_utilities_init
   integer :: structOrder                                                                            !< order of displacement shape functions
   character(len=*), parameter :: &
     PETSCDEBUG = ' -snes_view -snes_monitor '
-
   PetscErrorCode            :: ierr
+  logical :: debugPETSc                                                                             !< use some in debug defined options for more verbose PETSc solution
+
 
   print'(/,a)',   ' <<<+-  FEM_utilities init  -+>>>'
 

src/mesh/discretization_mesh.f90

@@ -40,6 +40,11 @@ module discretization_mesh
     mesh_maxNips                                                                                    !< max number of IPs in any CP element
 !!!! BEGIN DEPRECATED !!!!!
 
+  DM, public :: geomMesh
+
+  PetscInt, dimension(:), allocatable, public, protected :: &
+    mesh_boundaries
+
   real(pReal), dimension(:,:), allocatable :: &
     mesh_ipVolume, &                                                                                !< volume associated with IP (initially!)
     mesh_node0                                                                                      !< node x,y,z coordinates (initially!)
@@ -50,11 +55,6 @@ module discretization_mesh
   real(pReal), dimension(:,:,:), allocatable :: &
     mesh_ipCoordinates                                                                              !< IP x,y,z coordinates (after deformation!)
 
-  DM, public :: geomMesh
-
-  PetscInt, dimension(:), allocatable, public, protected :: &
-    mesh_boundaries
-
   public :: &
     discretization_mesh_init, &
     mesh_FEM_build_ipVolumes, &
@@ -71,16 +71,14 @@ subroutine discretization_mesh_init(restart)
 
   logical, intent(in) :: restart
 
-  integer, allocatable, dimension(:) :: chunkPos
   integer :: dimPlex, &
     mesh_Nnodes, &                                                                                  !< total number of nodes in mesh
-    j, l, &
+    j, &
     debug_element, debug_ip
   PetscSF :: sf
   DM :: globalMesh
   PetscInt :: nFaceSets
   PetscInt, pointer, dimension(:) :: pFaceSets
-  character(len=pStringLen), dimension(:), allocatable :: fileContent
   IS :: faceSetIS
   PetscErrorCode :: ierr
   integer, dimension(:), allocatable :: &
@@ -88,7 +86,7 @@ subroutine discretization_mesh_init(restart)
   class(tNode), pointer :: &
     num_mesh
   integer :: integrationOrder                                                                       !< order of quadrature rule required
-  type(tvec)  :: coords_node0
+  type(tvec) :: coords_node0
 
   print'(/,a)',   ' <<<+-  discretization_mesh init  -+>>>'
 

src/mesh/mesh_mech_FEM.f90

@@ -109,7 +109,7 @@ subroutine FEM_mechanical_init(fieldBC)
 
   character(len=*), parameter            :: prefix = 'mechFE_'
   PetscErrorCode                         :: ierr
-
+  real(pReal), dimension(3,3) :: devNull
   class(tNode), pointer :: &
     num_mesh
    
@@ -258,6 +258,7 @@ subroutine FEM_mechanical_init(fieldBC)
     call DMPlexVecSetClosure(mechanical_mesh,section,solution_local,cell,px_scal,5,ierr)
     CHKERRQ(ierr)
   enddo
+  call utilities_constitutiveResponse(0.0_pReal,devNull,.true.)
 
 end subroutine FEM_mechanical_init
 
@@ -288,8 +289,8 @@ type(tSolutionState) function FEM_mechanical_solution( &
   params%timeinc = timeinc
   params%fieldBC = fieldBC
 
-  call SNESSolve(mechanical_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr)                             ! solve mechanical_snes based on solution guess (result in solution)
+  call SNESSolve(mechanical_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr)                       ! solve mechanical_snes based on solution guess (result in solution)
-  call SNESGetConvergedReason(mechanical_snes,reason,ierr); CHKERRQ(ierr)                                 ! solution converged?
+  call SNESGetConvergedReason(mechanical_snes,reason,ierr); CHKERRQ(ierr)                           ! solution converged?
   terminallyIll = .false.
 
   if (reason < 1) then                                                                              ! 0: still iterating (will not occur), negative -> convergence error
@@ -397,7 +398,7 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,ierr)
 
 !--------------------------------------------------------------------------------------------------
 ! evaluate constitutive response
-  call Utilities_constitutiveResponse(params%timeinc,P_av,ForwardData)
+  call utilities_constitutiveResponse(params%timeinc,P_av,ForwardData)
   call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr)
   ForwardData = .false.
 
@@ -670,7 +671,7 @@ end subroutine FEM_mechanical_converged
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Calculate current coordinates (FEM nodal coordinates only at the moment)
+!> @brief Calculate current coordinates (both nodal and ip coordinates)
 !--------------------------------------------------------------------------------------------------
 subroutine FEM_mechanical_updateCoords()
 
@@ -678,21 +679,35 @@ subroutine FEM_mechanical_updateCoords()
     nodeCoords_linear                                                                               !< nodal coordinates (dimPlex*Nnodes)
   real(pReal), pointer, dimension(:,:) :: &
     nodeCoords                                                                                      !< nodal coordinates (3,Nnodes)
+  real(pReal), pointer, dimension(:,:,:) :: &                                                       
+    ipCoords                                                                                        !< ip coordinates (3,nQuadrature,mesh_NcpElems)
+
+  integer :: &
+    qPt, &
+    comp, &
+    qOffset, &
+    nOffset
 
   DM  :: dm_local
   Vec :: x_local
   PetscErrorCode :: ierr
-  PetscInt :: dimPlex, pStart, pEnd, p, s, e
+  PetscInt :: pStart, pEnd, p, s, e, q, &
+              cellStart, cellEnd, c, n
   PetscSection :: section
+  PetscQuadrature :: mechQuad
+  PetscReal,      dimension(:),  pointer :: basisField, basisFieldDer
+  PetscScalar, dimension(:),   pointer :: x_scal
 
   call SNESGetDM(mechanical_snes,dm_local,ierr); CHKERRQ(ierr)
+  call DMGetDS(dm_local,mechQuad,ierr); CHKERRQ(ierr)
   call DMGetLocalSection(dm_local,section,ierr); CHKERRQ(ierr)
   call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
   call DMGetDimension(dm_local,dimPlex,ierr); CHKERRQ(ierr)
+
+  ! write cell vertex displacements
   call DMPlexGetDepthStratum(dm_local,0,pStart,pEnd,ierr); CHKERRQ(ierr)
   allocate(nodeCoords(3,pStart:pEnd-1),source=0.0_pReal)
   call VecGetArrayF90(x_local,nodeCoords_linear,ierr); CHKERRQ(ierr)
-
   do p=pStart, pEnd-1
     call DMPlexGetPointLocal(dm_local, p, s, e, ierr); CHKERRQ(ierr)
     nodeCoords(1:dimPlex,p)=nodeCoords_linear(s+1:e)
@@ -700,6 +715,31 @@ subroutine FEM_mechanical_updateCoords()
 
   call discretization_setNodeCoords(nodeCoords)
   call VecRestoreArrayF90(x_local,nodeCoords_linear,ierr); CHKERRQ(ierr)
+
+  ! write ip displacements
+  call DMPlexGetHeightStratum(dm_local,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
+  call PetscDSGetTabulation(mechQuad,0,basisField,basisFieldDer,ierr); CHKERRQ(ierr)
+  allocate(ipCoords(3,nQuadrature,mesh_NcpElems),source=0.0_pReal)
+  do c=cellStart,cellEnd-1                                                                                      
+    qOffset=0
+    call DMPlexVecGetClosure(dm_local,section,x_local,c,x_scal,ierr); CHKERRQ(ierr)                             !< get nodal coordinates of each element
+    do qPt=0,nQuadrature-1
+      qOffset= qPt * (size(basisField)/nQuadrature)
+      do comp=0,dimPlex-1                                                                                       !< loop over components
+        nOffset=0
+        q = comp
+        do n=0,nBasis-1
+          ipCoords(comp+1,qPt+1,c+1)=ipCoords(comp+1,qPt+1,c+1)+&
+                                     sum(basisField(qOffset+(q*dimPlex)+1:qOffset+(q*dimPlex)+dimPlex)*&
+                                     x_scal(nOffset+1:nOffset+dimPlex))
+          q = q+dimPlex
+          nOffset = nOffset+dimPlex
+        enddo
+      enddo
+    enddo
+    call DMPlexVecRestoreClosure(dm_local,section,x_local,c,x_scal,ierr); CHKERRQ(ierr)     
+  end do    
+  call discretization_setIPcoords(reshape(ipCoords,[3,mesh_NcpElems*nQuadrature]))
   call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
 
 end subroutine FEM_mechanical_updateCoords

src/parallelization.f90

@@ -24,9 +24,18 @@ module parallelization
     worldrank = 0, &                                                                                !< MPI worldrank (/=0 for MPI simulations only)
     worldsize = 1                                                                                   !< MPI worldsize (/=1 for MPI simulations only)
 
-#ifdef PETSC
+#ifndef PETSC
+public :: parallelization_bcast_str
+
+contains
+subroutine parallelization_bcast_str(string)
+  character(len=:), allocatable, intent(inout) :: string
+end subroutine parallelization_bcast_str
+
+#else
   public :: &
-    parallelization_init
+    parallelization_init, &
+    parallelization_bcast_str
 
 contains
 
@@ -101,6 +110,27 @@ subroutine parallelization_init
 !$ call omp_set_num_threads(OMP_NUM_THREADS)
 
 end subroutine parallelization_init
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Broadcast a string from process 0.
+!--------------------------------------------------------------------------------------------------
+subroutine parallelization_bcast_str(string)
+
+  character(len=:), allocatable, intent(inout) :: string
+
+  integer :: strlen, ierr                                                                           ! pI64 for strlen not supported by MPI
+
+
+  if (worldrank == 0) strlen = len(string)
+  call MPI_Bcast(strlen,1,MPI_INTEGER,0,MPI_COMM_WORLD, ierr)
+  if (worldrank /= 0) allocate(character(len=strlen)::string)
+
+  call MPI_Bcast(string,strlen,MPI_CHARACTER,0,MPI_COMM_WORLD, ierr)
+
+
+end subroutine parallelization_bcast_str
+
 #endif
 
 end module parallelization

src/phase.f90

@@ -26,12 +26,8 @@ module phase
   real(pReal),      allocatable, dimension(:) :: phase_rho
 
   type(tRotationContainer), dimension(:), allocatable :: &
-    phase_orientation0, &
+    phase_O_0, &
-    phase_orientation
+    phase_O
-
-  type :: tTensorContainer
-    real(pReal), dimension(:,:,:), allocatable :: data
-  end type
 
   type :: tNumerics
     integer :: &
@@ -364,7 +360,7 @@ subroutine phase_init
   allocate(phase_lattice(phases%length))
   allocate(phase_cOverA(phases%length),source=-1.0_pReal)
   allocate(phase_rho(phases%length))
-  allocate(phase_orientation0(phases%length))
+  allocate(phase_O_0(phases%length))
 
   do ph = 1,phases%length
     phase => phases%get(ph)
@@ -374,20 +370,20 @@ subroutine phase_init
     if (any(phase_lattice(ph) == ['hP','tI'])) &
       phase_cOverA(ph) = phase%get_asFloat('c/a')
     phase_rho(ph) = phase%get_asFloat('rho',defaultVal=0.0_pReal)
-    allocate(phase_orientation0(ph)%data(count(material_phaseID==ph)))
+    allocate(phase_O_0(ph)%data(count(material_phaseID==ph)))
   enddo
 
   do ce = 1, size(material_phaseID,2)
     ma = discretization_materialAt((ce-1)/discretization_nIPs+1)
     do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
       ph = material_phaseID(co,ce)
-      phase_orientation0(ph)%data(material_phaseEntry(co,ce)) = material_orientation0(ma)%data(co)
+      phase_O_0(ph)%data(material_phaseEntry(co,ce)) = material_O_0(ma)%data(co)
     enddo
   enddo
 
-  allocate(phase_orientation(phases%length))
+  allocate(phase_O(phases%length))
   do ph = 1,phases%length
-    phase_orientation(ph)%data = phase_orientation0(ph)%data
+    phase_O(ph)%data = phase_O_0(ph)%data
   enddo
 
   call mechanical_init(materials,phases)
@@ -577,10 +573,10 @@ subroutine crystallite_orientations(co,ip,el)
   ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
   en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
 
-  call phase_orientation(ph)%data(en)%fromMatrix(transpose(math_rotationalPart(mechanical_F_e(ph,en))))
+  call phase_O(ph)%data(en)%fromMatrix(transpose(math_rotationalPart(mechanical_F_e(ph,en))))
 
   if (plasticState(material_phaseAt(1,el))%nonlocal) &
-    call plastic_nonlocal_updateCompatibility(phase_orientation,material_phaseAt(1,el),ip,el)
+    call plastic_nonlocal_updateCompatibility(phase_O,material_phaseAt(1,el),ip,el)
 
 
 end subroutine crystallite_orientations
@@ -602,7 +598,7 @@ function crystallite_push33ToRef(co,ce, tensor33)
 
   ph = material_phaseID(co,ce)
   en = material_phaseEntry(co,ce)
-  T = matmul(phase_orientation0(ph)%data(en)%asMatrix(),transpose(math_inv33(phase_F(co,ce))))      ! ToDo: initial orientation correct?
+  T = matmul(phase_O_0(ph)%data(en)%asMatrix(),transpose(math_inv33(phase_F(co,ce))))               ! ToDo: initial orientation correct?
 
   crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T))
 

src/phase_mechanical.f90

@@ -229,8 +229,8 @@ module subroutine mechanical_init(materials,phases)
   allocate(phase_mechanical_F0(phases%length))
   allocate(phase_mechanical_Li(phases%length))
   allocate(phase_mechanical_Li0(phases%length))
-  allocate(phase_mechanical_Lp0(phases%length))
   allocate(phase_mechanical_Lp(phases%length))
+  allocate(phase_mechanical_Lp0(phases%length))
   allocate(phase_mechanical_S(phases%length))
   allocate(phase_mechanical_P(phases%length))
   allocate(phase_mechanical_S0(phases%length))
@@ -238,20 +238,20 @@ module subroutine mechanical_init(materials,phases)
   do ph = 1, phases%length
     Nmembers = count(material_phaseID == ph)
 
-    allocate(phase_mechanical_Fi(ph)%data(3,3,Nmembers))
     allocate(phase_mechanical_Fe(ph)%data(3,3,Nmembers))
+    allocate(phase_mechanical_Fi(ph)%data(3,3,Nmembers))
     allocate(phase_mechanical_Fi0(ph)%data(3,3,Nmembers))
     allocate(phase_mechanical_Fp(ph)%data(3,3,Nmembers))
     allocate(phase_mechanical_Fp0(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_Li(ph)%data(3,3,Nmembers))
+    allocate(phase_mechanical_F(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_Li0(ph)%data(3,3,Nmembers))
+    allocate(phase_mechanical_F0(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_Lp0(ph)%data(3,3,Nmembers))
+    allocate(phase_mechanical_Li(ph)%data(3,3,Nmembers),source=0.0_pReal)
-    allocate(phase_mechanical_Lp(ph)%data(3,3,Nmembers))
+    allocate(phase_mechanical_Li0(ph)%data(3,3,Nmembers),source=0.0_pReal)
+    allocate(phase_mechanical_Lp(ph)%data(3,3,Nmembers),source=0.0_pReal)
+    allocate(phase_mechanical_Lp0(ph)%data(3,3,Nmembers),source=0.0_pReal)
     allocate(phase_mechanical_S(ph)%data(3,3,Nmembers),source=0.0_pReal)
     allocate(phase_mechanical_P(ph)%data(3,3,Nmembers),source=0.0_pReal)
     allocate(phase_mechanical_S0(ph)%data(3,3,Nmembers),source=0.0_pReal)
-    allocate(phase_mechanical_F(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_F0(ph)%data(3,3,Nmembers))
 
     phase   => phases%get(ph)
     mech    => phase%get('mechanical')
@@ -265,7 +265,7 @@ module subroutine mechanical_init(materials,phases)
   do ph = 1, phases%length
     do en = 1, count(material_phaseID == ph)
 
-      phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = phase_orientation0(ph)%data(en)%asMatrix()        ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005)
+      phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = phase_O_0(ph)%data(en)%asMatrix()                 ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005)
       phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) &
                                                 / math_det33(phase_mechanical_Fp0(ph)%data(1:3,1:3,en))**(1.0_pReal/3.0_pReal)
       phase_mechanical_Fi0(ph)%data(1:3,1:3,en) = math_I3
@@ -508,7 +508,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
     enddo LpLoop
 
     call phase_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, &
-                                       S, Fi_new, ph,en)
+                                S, Fi_new, ph,en)
 
     !* update current residuum and check for convergence of loop
     atol_Li = max(num%rtol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), &        ! absolute tolerance from largest acceptable relative error
@@ -902,7 +902,6 @@ subroutine crystallite_results(group,ph)
   integer,          intent(in) :: ph
 
   integer :: ou
-  real(pReal), allocatable, dimension(:,:)   :: selected_rotations
 
 
   call results_closeGroup(results_addGroup(group//'/mechanical'))
@@ -935,8 +934,7 @@ subroutine crystallite_results(group,ph)
         call results_writeDataset(phase_mechanical_S(ph)%data,group//'/mechanical/','S', &
                                  'second Piola-Kirchhoff stress','Pa')
       case('O')
-        selected_rotations = select_rotations(phase_orientation(ph)%data)
+        call results_writeDataset(to_quaternion(phase_O(ph)%data),group//'/mechanical',output_constituent(ph)%label(ou),&
-        call results_writeDataset(selected_rotations,group//'/mechanical',output_constituent(ph)%label(ou),&
                                  'crystal orientation as quaternion','q_0 (q_1 q_2 q_3)')
         call results_addAttribute('lattice',phase_lattice(ph),group//'/mechanical/'//output_constituent(ph)%label(ou))
         if (any(phase_lattice(ph) == ['hP', 'tI'])) &
@@ -948,19 +946,21 @@ subroutine crystallite_results(group,ph)
   contains
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Convert orientation for output: ToDo: implement in HDF5/results
+!> @brief Convert orientation array to quaternion array
 !--------------------------------------------------------------------------------------------------
-  function select_rotations(dataset)
+  function to_quaternion(dataset)
 
     type(rotation), dimension(:), intent(in) :: dataset
-    real(pReal), dimension(4,size(dataset,1)) :: select_rotations
+    real(pReal), dimension(4,size(dataset,1)) :: to_quaternion
-    integer :: en
 
-    do en = 1, size(dataset,1)
+    integer :: i
-      select_rotations(:,en) = dataset(en)%asQuaternion()
-   enddo
 
- end function select_rotations
+
+    do i = 1, size(dataset,1)
+      to_quaternion(:,i) = dataset(i)%asQuaternion()
+    enddo
+
+ end function to_quaternion
 
 end subroutine crystallite_results
 
@@ -1021,7 +1021,7 @@ module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged
 
   subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,en)
   subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,en)
-  subState0 = plasticState(ph)%State0(:,en)
+  allocate(subState0,source=plasticState(ph)%State0(:,en))
   subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,en)
   subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,en)
   subF0  = phase_mechanical_F0(ph)%data(1:3,1:3,en)
@@ -1144,12 +1144,12 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
   en = material_phaseEntry(co,ce)
 
   call phase_hooke_SandItsTangents(devNull,dSdFe,dSdFi, &
-                                          phase_mechanical_Fe(ph)%data(1:3,1:3,en), &
+                                   phase_mechanical_Fe(ph)%data(1:3,1:3,en), &
-                                          phase_mechanical_Fi(ph)%data(1:3,1:3,en),ph,en)
+                                   phase_mechanical_Fi(ph)%data(1:3,1:3,en),ph,en)
   call phase_LiAndItsTangents(devNull,dLidS,dLidFi, &
-                                     phase_mechanical_S(ph)%data(1:3,1:3,en), &
+                              phase_mechanical_S(ph)%data(1:3,1:3,en), &
-                                     phase_mechanical_Fi(ph)%data(1:3,1:3,en), &
+                              phase_mechanical_Fi(ph)%data(1:3,1:3,en), &
-                                     ph,en)
+                              ph,en)
 
   invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,en))
   invFi = math_inv33(phase_mechanical_Fi(ph)%data(1:3,1:3,en))

src/phase_mechanical_plastic_dislotungsten.f90

@@ -18,7 +18,7 @@ submodule(phase:plastic) dislotungsten
       Q_cl = 1.0_pReal                                                                              !< activation energy for dislocation climb
     real(pReal),               allocatable, dimension(:) :: &
       b_sl, &                                                                                       !< magnitude of Burgers vector [m]
-      D_a, &
+      d_caron, &                                                                                    !< distance of spontaneous annhihilation
       i_sl, &                                                                                       !< Adj. parameter for distance between 2 forest dislocations
       f_at, &                                                                                       !< factor to calculate atomic volume
       tau_Peierls, &                                                                                !< Peierls stress
@@ -56,7 +56,7 @@ submodule(phase:plastic) dislotungsten
   type :: tDisloTungstendependentState
     real(pReal), dimension(:,:), allocatable :: &
       Lambda_sl, &
-      threshold_stress
+      tau_pass
   end type tDisloTungstendependentState
 
 !--------------------------------------------------------------------------------------------------
@@ -172,7 +172,6 @@ module function plastic_dislotungsten_init() result(myPlasticity)
       prm%D_0  = pl%get_asFloat('D_0')
       prm%Q_cl = pl%get_asFloat('Q_cl')
       prm%f_at = pl%get_asFloat('f_at') * prm%b_sl**3.0_pReal
-      prm%D_a  = pl%get_asFloat('D_a') * prm%b_sl
 
       prm%dipoleformation = .not. pl%get_asBool('no_dipole_formation', defaultVal = .false.)
 
@@ -191,7 +190,7 @@ module function plastic_dislotungsten_init() result(myPlasticity)
       prm%B              = math_expand(prm%B,           N_sl)
       prm%i_sl           = math_expand(prm%i_sl,        N_sl)
       prm%f_at           = math_expand(prm%f_at,        N_sl)
-      prm%D_a            = math_expand(prm%D_a,         N_sl)
+      prm%d_caron        = pl%get_asFloat('D_a') * prm%b_sl
 
       ! sanity checks
       if (    prm%D_0          <= 0.0_pReal)  extmsg = trim(extmsg)//' D_0'
@@ -202,12 +201,13 @@ module function plastic_dislotungsten_init() result(myPlasticity)
       if (any(prm%b_sl         <= 0.0_pReal)) extmsg = trim(extmsg)//' b_sl'
       if (any(prm%Q_s          <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_s'
       if (any(prm%tau_Peierls  <  0.0_pReal)) extmsg = trim(extmsg)//' tau_Peierls'
-      if (any(prm%D_a          <= 0.0_pReal)) extmsg = trim(extmsg)//' D_a or b_sl'
+      if (any(prm%B            <  0.0_pReal)) extmsg = trim(extmsg)//' B'
+      if (any(prm%d_caron      <  0.0_pReal)) extmsg = trim(extmsg)//' d_caron(D_a,b_sl)'
       if (any(prm%f_at         <= 0.0_pReal)) extmsg = trim(extmsg)//' f_at or b_sl'
 
     else slipActive
       rho_mob_0= emptyRealArray; rho_dip_0 = emptyRealArray
-      allocate(prm%b_sl,prm%D_a,prm%i_sl,prm%f_at,prm%tau_Peierls, &
+      allocate(prm%b_sl,prm%d_caron,prm%i_sl,prm%f_at,prm%tau_Peierls, &
                prm%Q_s,prm%v_0,prm%p,prm%q,prm%B,prm%h,prm%w,prm%omega, &
                source = emptyRealArray)
       allocate(prm%forestProjection(0,0))
@@ -246,8 +246,8 @@ module function plastic_dislotungsten_init() result(myPlasticity)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
     if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
 
-    allocate(dst%Lambda_sl(prm%sum_N_sl,Nmembers),         source=0.0_pReal)
+    allocate(dst%Lambda_sl(prm%sum_N_sl,Nmembers), source=0.0_pReal)
-    allocate(dst%threshold_stress(prm%sum_N_sl,Nmembers),  source=0.0_pReal)
+    allocate(dst%tau_pass(prm%sum_N_sl,Nmembers),  source=0.0_pReal)
 
     end associate
 
@@ -316,8 +316,6 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
     ph, &
     en
 
-  real(pReal) :: &
-    VacancyDiffusion
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
     dot_gamma_pos, dot_gamma_neg,&
     tau_pos,&
@@ -325,38 +323,36 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
     v_cl, &
     dot_rho_dip_formation, &
     dot_rho_dip_climb, &
-    dip_distance
+    d_hat
 
-  associate(prm => param(ph), stt => state(ph),&
+  associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph))
-            dot => dotState(ph), dst => dependentState(ph))
 
   call kinetics(Mp,T,ph,en,&
                 dot_gamma_pos,dot_gamma_neg, &
                 tau_pos_out = tau_pos,tau_neg_out = tau_neg)
 
-  dot%gamma_sl(:,en) = (dot_gamma_pos+dot_gamma_neg)                                                          ! ToDo: needs to be abs
+  dot%gamma_sl(:,en) = abs(dot_gamma_pos+dot_gamma_neg)
-  VacancyDiffusion = prm%D_0*exp(-prm%Q_cl/(kB*T))
 
-  where(dEq0(tau_pos))                                                                              ! ToDo: use avg of pos and neg
+  where(dEq0(tau_pos))                                                                              ! ToDo: use avg of +/-
     dot_rho_dip_formation = 0.0_pReal
     dot_rho_dip_climb     = 0.0_pReal
   else where
-    dip_distance = math_clip(3.0_pReal*prm%mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos)), &
+    d_hat = math_clip(3.0_pReal*prm%mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos)), &                     ! ToDo: use avg of +/-
-                             prm%D_a, &                                                             ! lower limit
+                      prm%d_caron, &                                                                ! lower limit
-                             dst%Lambda_sl(:,en))                                                   ! upper limit
+                      dst%Lambda_sl(:,en))                                                          ! upper limit
-    dot_rho_dip_formation = merge(2.0_pReal*dip_distance* stt%rho_mob(:,en)*abs(dot%gamma_sl(:,en))/prm%b_sl, & ! ToDo: ignore region of spontaneous annihilation
+    dot_rho_dip_formation = merge(2.0_pReal*(d_hat-prm%d_caron)*stt%rho_mob(:,en)*dot%gamma_sl(:,en)/prm%b_sl, &
                                   0.0_pReal, &
                                   prm%dipoleformation)
-    v_cl = (3.0_pReal*prm%mu*VacancyDiffusion*prm%f_at/(2.0_pReal*pi*kB*T)) &
+    v_cl = (3.0_pReal*prm%mu*prm%D_0*exp(-prm%Q_cl/(kB*T))*prm%f_at/(2.0_pReal*PI*kB*T)) &
-                  * (1.0_pReal/(dip_distance+prm%D_a))
+         * (1.0_pReal/(d_hat+prm%d_caron))
-    dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,en))/(dip_distance-prm%D_a)                   ! ToDo: Discuss with Franz: Stress dependency?
+    dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,en))/(d_hat-prm%d_caron)                      ! ToDo: Discuss with Franz: Stress dependency?
   end where
 
-  dot%rho_mob(:,en) = abs(dot%gamma_sl(:,en))/(prm%b_sl*dst%Lambda_sl(:,en)) &                      ! multiplication
+  dot%rho_mob(:,en) = dot%gamma_sl(:,en)/(prm%b_sl*dst%Lambda_sl(:,en)) &                           ! multiplication
                     - dot_rho_dip_formation &
-                    - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_mob(:,en)*abs(dot%gamma_sl(:,en))        ! Spontaneous annihilation of 2 single edge dislocations
+                    - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_mob(:,en)*dot%gamma_sl(:,en)         ! Spontaneous annihilation of 2 edges
   dot%rho_dip(:,en) = dot_rho_dip_formation &
-                    - (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_dip(:,en)*abs(dot%gamma_sl(:,en)) &      ! Spontaneous annihilation of a single edge dislocation with a dipole constituent
+                    - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_dip(:,en)*dot%gamma_sl(:,en) &       ! Spontaneous annihilation of an edge with a dipole
                     - dot_rho_dip_climb
 
   end associate
@@ -377,11 +373,11 @@ module subroutine dislotungsten_dependentState(ph,en)
     dislocationSpacing
 
 
-  associate(prm => param(ph), stt => state(ph),dst => dependentState(ph))
+  associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
 
     dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
-    dst%threshold_stress(:,en) = prm%mu*prm%b_sl &
+    dst%tau_pass(:,en) = prm%mu*prm%b_sl &
-                               * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
+                       * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
 
     dst%Lambda_sl(:,en) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl)
 
@@ -416,7 +412,7 @@ module subroutine plastic_dislotungsten_results(ph,group)
         if(prm%sum_N_sl>0) call results_writeDataset(dst%Lambda_sl,group,trim(prm%output(o)), &
                                                      'mean free path for slip','m')
       case('tau_pass')
-        if(prm%sum_N_sl>0) call results_writeDataset(dst%threshold_stress,group,trim(prm%output(o)), &
+        if(prm%sum_N_sl>0) call results_writeDataset(dst%tau_pass,group,trim(prm%output(o)), &
                                                      'threshold stress for slip','Pa')
     end select
   enddo outputsLoop
@@ -456,8 +452,7 @@ pure subroutine kinetics(Mp,T,ph,en, &
     StressRatio_p,StressRatio_pminus1, &
     dvel, vel, &
     tau_pos,tau_neg, &
-    t_n, t_k, dtk,dtn, &
+    t_n, t_k, dtk,dtn
-    needsGoodName                                                                                   ! ToDo: @Karo: any idea?
   integer :: j
 
   associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
@@ -475,13 +470,12 @@ pure subroutine kinetics(Mp,T,ph,en, &
             dot_gamma_0     => stt%rho_mob(:,en)*prm%b_sl*prm%v_0, &
             effectiveLength => dst%Lambda_sl(:,en) - prm%w)
 
-  significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,en) > tol_math_check)
+  significantPositiveTau: where(abs(tau_pos)-dst%tau_pass(:,en) > tol_math_check)
-    StressRatio = (abs(tau_pos)-dst%threshold_stress(:,en))/prm%tau_Peierls
+    StressRatio = (abs(tau_pos)-dst%tau_pass(:,en))/prm%tau_Peierls
     StressRatio_p       = StressRatio** prm%p
     StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal)
-    needsGoodName       = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)
 
-    t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength)
+    t_n = prm%b_sl/(exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)*prm%omega*effectiveLength)
     t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos)
 
     vel = prm%h/(t_n + t_k)
@@ -492,7 +486,7 @@ pure subroutine kinetics(Mp,T,ph,en, &
   end where significantPositiveTau
 
   if (present(ddot_gamma_dtau_pos)) then
-    significantPositiveTau2: where(abs(tau_pos)-dst%threshold_stress(:,en) > tol_math_check)
+    significantPositiveTau2: where(abs(tau_pos)-dst%tau_pass(:,en) > tol_math_check)
       dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) &
           * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls
       dtk = -1.0_pReal * t_k / tau_pos
@@ -505,13 +499,12 @@ pure subroutine kinetics(Mp,T,ph,en, &
     end where significantPositiveTau2
   endif
 
-  significantNegativeTau: where(abs(tau_neg)-dst%threshold_stress(:,en) > tol_math_check)
+  significantNegativeTau: where(abs(tau_neg)-dst%tau_pass(:,en) > tol_math_check)
-    StressRatio = (abs(tau_neg)-dst%threshold_stress(:,en))/prm%tau_Peierls
+    StressRatio = (abs(tau_neg)-dst%tau_pass(:,en))/prm%tau_Peierls
     StressRatio_p       = StressRatio** prm%p
     StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal)
-    needsGoodName       = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)
 
-    t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength)
+    t_n = prm%b_sl/(exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)*prm%omega*effectiveLength)
     t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos)
 
     vel = prm%h/(t_n + t_k)
@@ -522,7 +515,7 @@ pure subroutine kinetics(Mp,T,ph,en, &
   end where significantNegativeTau
 
   if (present(ddot_gamma_dtau_neg)) then
-    significantNegativeTau2: where(abs(tau_neg)-dst%threshold_stress(:,en) > tol_math_check)
+    significantNegativeTau2: where(abs(tau_neg)-dst%tau_pass(:,en) > tol_math_check)
       dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) &
           * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls
       dtk = -1.0_pReal * t_k / tau_neg

src/phase_mechanical_plastic_dislotwin.f90

@@ -16,13 +16,11 @@ submodule(phase:plastic) dislotwin
     real(pReal) :: &
       mu                  = 1.0_pReal, &                                                            !< equivalent shear modulus
       nu                  = 1.0_pReal, &                                                            !< equivalent shear Poisson's ratio
-      D_0                 = 1.0_pReal, &                                                            !< prefactor for self-diffusion coefficient
       Q_cl                = 1.0_pReal, &                                                            !< activation energy for dislocation climb
       omega               = 1.0_pReal, &                                                            !< frequency factor for dislocation climb
       D                   = 1.0_pReal, &                                                            !< grain size
       p_sb                = 1.0_pReal, &                                                            !< p-exponent in shear band velocity
       q_sb                = 1.0_pReal, &                                                            !< q-exponent in shear band velocity
-      D_a                 = 1.0_pReal, &                                                            !< adjustment parameter to calculate minimum dipole distance
       i_tw                = 1.0_pReal, &                                                            !< adjustment parameter to calculate MFP for twinning
       L_tw                = 1.0_pReal, &                                                            !< Length of twin nuclei in Burgers vectors
       L_tr                = 1.0_pReal, &                                                            !< Length of trans nuclei in Burgers vectors
@@ -42,7 +40,7 @@ submodule(phase:plastic) dislotwin
       b_sl, &                                                                                       !< absolute length of Burgers vector [m] for each slip system
       b_tw, &                                                                                       !< absolute length of Burgers vector [m] for each twin system
       b_tr, &                                                                                       !< absolute length of Burgers vector [m] for each transformation system
-      Q_s,&                                                                                         !< activation energy for glide [J] for each slip system
+      Q_sl,&                                                                                        !< activation energy for glide [J] for each slip system
       v_0, &                                                                                        !< dislocation velocity prefactor [m/s] for each slip system
       dot_N_0_tw, &                                                                                 !< twin nucleation rate [1/m³s] for each twin system
       dot_N_0_tr, &                                                                                 !< trans nucleation rate [1/m³s] for each trans system
@@ -55,7 +53,8 @@ submodule(phase:plastic) dislotwin
       s, &                                                                                          !< s-exponent in trans nucleation rate
       tau_0, &                                                                                      !< strength due to elements in solid solution
       gamma_char, &                                                                                 !< characteristic shear for twins
-      B                                                                                             !< drag coefficient
+      B, &                                                                                          !< drag coefficient
+      d_caron                                                                                       !< distance of spontaneous annhihilation
     real(pReal),               allocatable, dimension(:,:) :: &
       h_sl_sl, &                                                                                    !< components of slip-slip interaction matrix
       h_sl_tw, &                                                                                    !< components of slip-twin interaction matrix
@@ -206,7 +205,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
       rho_dip_0                = pl%get_as1dFloat('rho_dip_0',   requiredSize=size(N_sl))
       prm%v_0                  = pl%get_as1dFloat('v_0',         requiredSize=size(N_sl))
       prm%b_sl                 = pl%get_as1dFloat('b_sl',        requiredSize=size(N_sl))
-      prm%Q_s                  = pl%get_as1dFloat('Q_s',         requiredSize=size(N_sl))
+      prm%Q_sl                 = pl%get_as1dFloat('Q_sl',        requiredSize=size(N_sl))
       prm%i_sl                 = pl%get_as1dFloat('i_sl',        requiredSize=size(N_sl))
       prm%p                    = pl%get_as1dFloat('p_sl',        requiredSize=size(N_sl))
       prm%q                    = pl%get_as1dFloat('q_sl',        requiredSize=size(N_sl))
@@ -214,9 +213,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
       prm%B                    = pl%get_as1dFloat('B',           requiredSize=size(N_sl), &
                                                   defaultVal=[(0.0_pReal, i=1,size(N_sl))])
 
-      prm%D_a                  = pl%get_asFloat('D_a')
+      prm%Q_cl = pl%get_asFloat('Q_cl')
-      prm%D_0                  = pl%get_asFloat('D_0')
-      prm%Q_cl                 = pl%get_asFloat('Q_cl')
 
       prm%ExtendedDislocations = pl%get_asBool('extend_dislocations',defaultVal = .false.)
       prm%omitDipoles          = pl%get_asBool('omit_dipoles',defaultVal = .false.)
@@ -231,28 +228,29 @@ module function plastic_dislotwin_init() result(myPlasticity)
       rho_dip_0        = math_expand(rho_dip_0,       N_sl)
       prm%v_0          = math_expand(prm%v_0,         N_sl)
       prm%b_sl         = math_expand(prm%b_sl,        N_sl)
-      prm%Q_s          = math_expand(prm%Q_s,         N_sl)
+      prm%Q_sl         = math_expand(prm%Q_sl,        N_sl)
       prm%i_sl         = math_expand(prm%i_sl,        N_sl)
       prm%p            = math_expand(prm%p,           N_sl)
       prm%q            = math_expand(prm%q,           N_sl)
       prm%tau_0        = math_expand(prm%tau_0,       N_sl)
       prm%B            = math_expand(prm%B,           N_sl)
+      prm%d_caron      = pl%get_asFloat('D_a') * prm%b_sl
 
       ! sanity checks
-      if (    prm%D_0           <= 0.0_pReal)          extmsg = trim(extmsg)//' D_0'
       if (    prm%Q_cl          <= 0.0_pReal)          extmsg = trim(extmsg)//' Q_cl'
       if (any(rho_mob_0         <  0.0_pReal))         extmsg = trim(extmsg)//' rho_mob_0'
       if (any(rho_dip_0         <  0.0_pReal))         extmsg = trim(extmsg)//' rho_dip_0'
       if (any(prm%v_0           <  0.0_pReal))         extmsg = trim(extmsg)//' v_0'
       if (any(prm%b_sl          <= 0.0_pReal))         extmsg = trim(extmsg)//' b_sl'
-      if (any(prm%Q_s           <= 0.0_pReal))         extmsg = trim(extmsg)//' Q_s'
+      if (any(prm%Q_sl           <= 0.0_pReal))        extmsg = trim(extmsg)//' Q_sl'
       if (any(prm%i_sl          <= 0.0_pReal))         extmsg = trim(extmsg)//' i_sl'
       if (any(prm%B             <  0.0_pReal))         extmsg = trim(extmsg)//' B'
+      if (any(prm%d_caron       <  0.0_pReal))         extmsg = trim(extmsg)//' d_caron(D_a,b_sl)'
       if (any(prm%p<=0.0_pReal .or. prm%p>1.0_pReal))  extmsg = trim(extmsg)//' p_sl'
       if (any(prm%q< 1.0_pReal .or. prm%q>2.0_pReal))  extmsg = trim(extmsg)//' q_sl'
     else slipActive
       rho_mob_0 = emptyRealArray; rho_dip_0 = emptyRealArray
-      allocate(prm%b_sl,prm%Q_s,prm%v_0,prm%i_sl,prm%p,prm%q,prm%B,source=emptyRealArray)
+      allocate(prm%b_sl,prm%Q_sl,prm%v_0,prm%i_sl,prm%p,prm%q,prm%B,source=emptyRealArray)
       allocate(prm%forestProjection(0,0),prm%h_sl_sl(0,0))
     endif slipActive
 
@@ -516,7 +514,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
   integer :: i,k,l,m,n
   real(pReal) :: &
      f_unrotated,StressRatio_p,&
-     BoltzmannRatio, &
+     E_kB_T, &
      ddot_gamma_dtau, &
      tau
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
@@ -586,7 +584,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
 
   shearBandingContribution: if(dNeq0(prm%v_sb)) then
 
-    BoltzmannRatio = prm%E_sb/(kB*T)
+    E_kB_T = prm%E_sb/(kB*T)
     call math_eigh33(eigValues,eigVectors,Mp)                                                       ! is Mp symmetric by design?
 
     do i = 1,6
@@ -596,8 +594,8 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
 
       significantShearBandStress: if (abs(tau) > tol_math_check) then
         StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb
-        dot_gamma_sb = sign(prm%v_sb*exp(-BoltzmannRatio*(1-StressRatio_p)**prm%q_sb), tau)
+        dot_gamma_sb = sign(prm%v_sb*exp(-E_kB_T*(1-StressRatio_p)**prm%q_sb), tau)
-        ddot_gamma_dtau = abs(dot_gamma_sb)*BoltzmannRatio* prm%p_sb*prm%q_sb/ prm%xi_sb &
+        ddot_gamma_dtau = abs(dot_gamma_sb)*E_kB_T*prm%p_sb*prm%q_sb/prm%xi_sb &
                         * (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) &
                         * (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal)
 
@@ -631,7 +629,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
   integer :: i
   real(pReal) :: &
     f_unrotated, &
-    rho_dip_distance, &
+    d_hat, &
     v_cl, &                                                                                         !< climb velocity
     tau, &
     sigma_cl, &                                                                                     !< climb stress
@@ -639,70 +637,67 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
     dot_rho_dip_formation, &
     dot_rho_dip_climb, &
-    rho_dip_distance_min, &
     dot_gamma_sl
   real(pReal), dimension(param(ph)%sum_N_tw) :: &
     dot_gamma_tw
   real(pReal), dimension(param(ph)%sum_N_tr) :: &
     dot_gamma_tr
 
-  associate(prm => param(ph),    stt => state(ph), &
-            dot => dotState(ph), dst => dependentState(ph))
 
-  f_unrotated = 1.0_pReal &
+  associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph))
-              - sum(stt%f_tw(1:prm%sum_N_tw,en)) &
-              - sum(stt%f_tr(1:prm%sum_N_tr,en))
 
-  call kinetics_sl(Mp,T,ph,en,dot_gamma_sl)
+    f_unrotated = 1.0_pReal &
-  dot%gamma_sl(:,en) = abs(dot_gamma_sl)
+                - sum(stt%f_tw(1:prm%sum_N_tw,en)) &
+                - sum(stt%f_tr(1:prm%sum_N_tr,en))
 
-  rho_dip_distance_min = prm%D_a*prm%b_sl
+    call kinetics_sl(Mp,T,ph,en,dot_gamma_sl)
+    dot%gamma_sl(:,en) = abs(dot_gamma_sl)
 
-  slipState: do i = 1, prm%sum_N_sl
+    slipState: do i = 1, prm%sum_N_sl
-    tau = math_tensordot(Mp,prm%P_sl(1:3,1:3,i))
+      tau = math_tensordot(Mp,prm%P_sl(1:3,1:3,i))
 
-    significantSlipStress: if (dEq0(tau) .or. prm%omitDipoles) then
+      significantSlipStress: if (dEq0(tau) .or. prm%omitDipoles) then
-      dot_rho_dip_formation(i) = 0.0_pReal
+        dot_rho_dip_formation(i) = 0.0_pReal
-      dot_rho_dip_climb(i) = 0.0_pReal
-    else significantSlipStress
-      rho_dip_distance = 3.0_pReal*prm%mu*prm%b_sl(i)/(16.0_pReal*PI*abs(tau))
-      rho_dip_distance = math_clip(rho_dip_distance, right = dst%Lambda_sl(i,en))
-      rho_dip_distance = math_clip(rho_dip_distance, left  = rho_dip_distance_min(i))
-
-      dot_rho_dip_formation(i) = 2.0_pReal*(rho_dip_distance-rho_dip_distance_min(i))/prm%b_sl(i) &
-                               * stt%rho_mob(i,en)*abs(dot_gamma_sl(i))
-
-      if (dEq(rho_dip_distance,rho_dip_distance_min(i))) then
         dot_rho_dip_climb(i) = 0.0_pReal
-      else
+      else significantSlipStress
-        ! Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981
+        d_hat = 3.0_pReal*prm%mu*prm%b_sl(i)/(16.0_pReal*PI*abs(tau))
-        sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i)))
+        d_hat = math_clip(d_hat, right = dst%Lambda_sl(i,en))
-        b_d = merge(24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu) &
+        d_hat = math_clip(d_hat, left  = prm%d_caron(i))
-                      * (prm%Gamma_sf(1) + prm%Gamma_sf(2) * T) / (prm%mu*prm%b_sl(i)), &
-                    1.0_pReal, &
-                    prm%ExtendedDislocations)
-        v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(kB*T)) &
-             * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal)
 
-        dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,en) &
+        dot_rho_dip_formation(i) = 2.0_pReal*(d_hat-prm%d_caron(i))/prm%b_sl(i) &
-                             / (rho_dip_distance-rho_dip_distance_min(i))
+                                 * stt%rho_mob(i,en)*abs(dot_gamma_sl(i))
-      endif
-    endif significantSlipStress
-  enddo slipState
 
-  dot%rho_mob(:,en) = abs(dot_gamma_sl)/(prm%b_sl*dst%Lambda_sl(:,en)) &
+        if (dEq(d_hat,prm%d_caron(i))) then
-                    - dot_rho_dip_formation &
+          dot_rho_dip_climb(i) = 0.0_pReal
-                    - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_mob(:,en)*abs(dot_gamma_sl)
+        else
+          ! Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981
+          sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i)))
+          b_d = merge(24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu) &
+                        * (prm%Gamma_sf(1) + prm%Gamma_sf(2) * T) / (prm%mu*prm%b_sl(i)), &
+                      1.0_pReal, &
+                      prm%ExtendedDislocations)
+          v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(kB*T)) &
+               * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal)
 
-  dot%rho_dip(:,en) = dot_rho_dip_formation &
+          dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,en) &
-                    - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,en)*abs(dot_gamma_sl) &
+                               / (d_hat-prm%d_caron(i))
-                    - dot_rho_dip_climb
+        endif
+      endif significantSlipStress
+    enddo slipState
 
-  call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw)
+    dot%rho_mob(:,en) = abs(dot_gamma_sl)/(prm%b_sl*dst%Lambda_sl(:,en)) &
-  dot%f_tw(:,en) = f_unrotated*dot_gamma_tw/prm%gamma_char
+                      - dot_rho_dip_formation &
+                      - 2.0_pReal*prm%d_caron/prm%b_sl * stt%rho_mob(:,en)*abs(dot_gamma_sl)
 
-  call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr)
+    dot%rho_dip(:,en) = dot_rho_dip_formation &
-  dot%f_tr(:,en) = f_unrotated*dot_gamma_tr
+                      - 2.0_pReal*prm%d_caron/prm%b_sl * stt%rho_dip(:,en)*abs(dot_gamma_sl) &
+                      - dot_rho_dip_climb
+
+    call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw)
+    dot%f_tw(:,en) = f_unrotated*dot_gamma_tw/prm%gamma_char
+
+    call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr)
+    dot%f_tr(:,en) = f_unrotated*dot_gamma_tr
 
   end associate
 
@@ -763,19 +758,17 @@ module subroutine dislotwin_dependentState(T,ph,en)
     dst%tau_pass(:,en) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
 
     !* threshold stress for growing twin/martensite
-    if(prm%sum_N_tw == prm%sum_N_sl) &
+    dst%tau_hat_tw(:,en) = Gamma/(3.0_pReal*prm%b_tw) &
-      dst%tau_hat_tw(:,en) = Gamma/(3.0_pReal*prm%b_tw) &
+                         + 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_tw)
-                           + 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_sl) ! slip Burgers here correct?
+    dst%tau_hat_tr(:,en) = Gamma/(3.0_pReal*prm%b_tr) &
-    if(prm%sum_N_tr == prm%sum_N_sl) &
+                         + 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_tr) &
-      dst%tau_hat_tr(:,en) = Gamma/(3.0_pReal*prm%b_tr) &
+                         + prm%h*prm%delta_G/(3.0_pReal*prm%b_tr)
-                           + 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip Burgers here correct?
-                           + prm%h*prm%delta_G/ (3.0_pReal*prm%b_tr)
 
     dst%V_tw(:,en) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,en)**2.0_pReal
     dst%V_tr(:,en) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,en)**2.0_pReal
 
 
-    x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu)      ! ToDo: In the paper, this is the Burgers vector for slip and is the same for twin and trans
+    x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu)      ! ToDo: In the paper, this is the Burgers vector for slip
     dst%tau_r_tw(:,en) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0)
 
     x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu)      ! ToDo: In the paper, this is the Burgers vector for slip
@@ -867,7 +860,7 @@ pure subroutine kinetics_sl(Mp,T,ph,en, &
     tau, &
     stressRatio, &
     StressRatio_p, &
-    BoltzmannRatio, &
+    Q_kB_T, &
     v_wait_inverse, &                                                                               !< inverse of the effective velocity of a dislocation waiting at obstacles (unsigned)
     v_run_inverse, &                                                                                !< inverse of the velocity of a free moving dislocation (unsigned)
     dV_wait_inverse_dTau, &
@@ -876,33 +869,34 @@ pure subroutine kinetics_sl(Mp,T,ph,en, &
     tau_eff                                                                                         !< effective resolved stress
   integer :: i
 
+
   associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
 
-  tau = [(math_tensordot(Mp,prm%P_sl(1:3,1:3,i)),i = 1, prm%sum_N_sl)]
+    tau = [(math_tensordot(Mp,prm%P_sl(1:3,1:3,i)),i = 1, prm%sum_N_sl)]
 
-  tau_eff = abs(tau)-dst%tau_pass(:,en)
+    tau_eff = abs(tau)-dst%tau_pass(:,en)
 
-  significantStress: where(tau_eff > tol_math_check)
+    significantStress: where(tau_eff > tol_math_check)
-    stressRatio    = tau_eff/prm%tau_0
+      stressRatio    = tau_eff/prm%tau_0
-    StressRatio_p  = stressRatio** prm%p
+      StressRatio_p  = stressRatio** prm%p
-    BoltzmannRatio = prm%Q_s/(kB*T)
+      Q_kB_T = prm%Q_sl/(kB*T)
-    v_wait_inverse = prm%v_0**(-1.0_pReal) * exp(BoltzmannRatio*(1.0_pReal-StressRatio_p)** prm%q)
+      v_wait_inverse = prm%v_0**(-1.0_pReal) * exp(Q_kB_T*(1.0_pReal-StressRatio_p)** prm%q)
-    v_run_inverse  = prm%B/(tau_eff*prm%b_sl)
+      v_run_inverse  = prm%B/(tau_eff*prm%b_sl)
 
-    dot_gamma_sl = sign(stt%rho_mob(:,en)*prm%b_sl/(v_wait_inverse+v_run_inverse),tau)
+      dot_gamma_sl = sign(stt%rho_mob(:,en)*prm%b_sl/(v_wait_inverse+v_run_inverse),tau)
 
-    dV_wait_inverse_dTau = -1.0_pReal * v_wait_inverse * prm%p * prm%q * BoltzmannRatio &
+      dV_wait_inverse_dTau = -1.0_pReal * v_wait_inverse * prm%p * prm%q * Q_kB_T &
-                         * (stressRatio**(prm%p-1.0_pReal)) &
+                           * (stressRatio**(prm%p-1.0_pReal)) &
-                         * (1.0_pReal-StressRatio_p)**(prm%q-1.0_pReal) &
+                           * (1.0_pReal-StressRatio_p)**(prm%q-1.0_pReal) &
-                         / prm%tau_0
+                           / prm%tau_0
-    dV_run_inverse_dTau  = -1.0_pReal * v_run_inverse/tau_eff
+      dV_run_inverse_dTau  = -1.0_pReal * v_run_inverse/tau_eff
-    dV_dTau              = -1.0_pReal * (dV_wait_inverse_dTau+dV_run_inverse_dTau) &
+      dV_dTau              = -1.0_pReal * (dV_wait_inverse_dTau+dV_run_inverse_dTau) &
-                         / (v_wait_inverse+v_run_inverse)**2.0_pReal
+                           / (v_wait_inverse+v_run_inverse)**2.0_pReal
-    ddot_gamma_dtau = dV_dTau*stt%rho_mob(:,en)*prm%b_sl
+      ddot_gamma_dtau = dV_dTau*stt%rho_mob(:,en)*prm%b_sl
-  else where significantStress
+    else where significantStress
-    dot_gamma_sl    = 0.0_pReal
+      dot_gamma_sl    = 0.0_pReal
-    ddot_gamma_dtau = 0.0_pReal
+      ddot_gamma_dtau = 0.0_pReal
-  end where significantStress
+    end where significantStress
 
   end associate
 
@@ -945,34 +939,35 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
 
   integer :: i,s1,s2
 
+
   associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
 
-  do i = 1, prm%sum_N_tw
+    do i = 1, prm%sum_N_tw
-    tau(i) = math_tensordot(Mp,prm%P_tw(1:3,1:3,i))
+      tau(i) = math_tensordot(Mp,prm%P_tw(1:3,1:3,i))
-    isFCC: if (prm%fccTwinTransNucleation) then
+      isFCC: if (prm%fccTwinTransNucleation) then
-      s1=prm%fcc_twinNucleationSlipPair(1,i)
+        s1=prm%fcc_twinNucleationSlipPair(1,i)
-      s2=prm%fcc_twinNucleationSlipPair(2,i)
+        s2=prm%fcc_twinNucleationSlipPair(2,i)
-      if (tau(i) < dst%tau_r_tw(i,en)) then                                                         ! ToDo: correct?
+        if (tau(i) < dst%tau_r_tw(i,en)) then                                                       ! ToDo: correct?
-        Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+&
+          Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+&
-               abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/&                     ! ToDo: MD: it would be more consistent to use shearrates from state
+                 abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/&
-                (prm%L_tw*prm%b_sl(i))*&
+                  (prm%L_tw*prm%b_sl(i))*&
-                (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tw(i,en)-tau(i))))                       ! P_ncs
+                  (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tw(i,en)-tau(i))))
-      else
+        else
-        Ndot0=0.0_pReal
+          Ndot0=0.0_pReal
-      end if
+        end if
-    else isFCC
+      else isFCC
-      Ndot0=prm%dot_N_0_tw(i)
+        Ndot0=prm%dot_N_0_tw(i)
-    endif isFCC
+      endif isFCC
-  enddo
+    enddo
 
-  significantStress: where(tau > tol_math_check)
+    significantStress: where(tau > tol_math_check)
-    StressRatio_r   = (dst%tau_hat_tw(:,en)/tau)**prm%r
+      StressRatio_r   = (dst%tau_hat_tw(:,en)/tau)**prm%r
-    dot_gamma_tw    = prm%gamma_char * dst%V_tw(:,en) * Ndot0*exp(-StressRatio_r)
+      dot_gamma_tw    = prm%gamma_char * dst%V_tw(:,en) * Ndot0*exp(-StressRatio_r)
-    ddot_gamma_dtau = (dot_gamma_tw*prm%r/tau)*StressRatio_r
+      ddot_gamma_dtau = (dot_gamma_tw*prm%r/tau)*StressRatio_r
-  else where significantStress
+    else where significantStress
-    dot_gamma_tw    = 0.0_pReal
+      dot_gamma_tw    = 0.0_pReal
-    ddot_gamma_dtau = 0.0_pReal
+      ddot_gamma_dtau = 0.0_pReal
-  end where significantStress
+    end where significantStress
 
   end associate
 
@@ -1011,36 +1006,37 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
     Ndot0, &
     stressRatio_s, &
     ddot_gamma_dtau
-
   integer :: i,s1,s2
+
+
   associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
 
-  do i = 1, prm%sum_N_tr
+    do i = 1, prm%sum_N_tr
-    tau(i) = math_tensordot(Mp,prm%P_tr(1:3,1:3,i))
+      tau(i) = math_tensordot(Mp,prm%P_tr(1:3,1:3,i))
-    isFCC: if (prm%fccTwinTransNucleation) then
+      isFCC: if (prm%fccTwinTransNucleation) then
-      s1=prm%fcc_twinNucleationSlipPair(1,i)
+        s1=prm%fcc_twinNucleationSlipPair(1,i)
-      s2=prm%fcc_twinNucleationSlipPair(2,i)
+        s2=prm%fcc_twinNucleationSlipPair(2,i)
-      if (tau(i) < dst%tau_r_tr(i,en)) then                                                         ! ToDo: correct?
+        if (tau(i) < dst%tau_r_tr(i,en)) then                                                       ! ToDo: correct?
-        Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+&
+          Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+&
-               abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/&                     ! ToDo: MD: it would be more consistent to use shearrates from state
+                 abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/&
-                (prm%L_tr*prm%b_sl(i))*&
+                  (prm%L_tr*prm%b_sl(i))*&
-                (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tr(i,en)-tau(i))))                       ! P_ncs
+                  (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tr(i,en)-tau(i))))
-      else
+        else
-        Ndot0=0.0_pReal
+          Ndot0=0.0_pReal
-      end if
+        end if
-    else isFCC
+      else isFCC
-      Ndot0=prm%dot_N_0_tr(i)
+        Ndot0=prm%dot_N_0_tr(i)
-    endif isFCC
+      endif isFCC
-  enddo
+    enddo
 
-  significantStress: where(tau > tol_math_check)
+    significantStress: where(tau > tol_math_check)
-    StressRatio_s   = (dst%tau_hat_tr(:,en)/tau)**prm%s
+      StressRatio_s   = (dst%tau_hat_tr(:,en)/tau)**prm%s
-    dot_gamma_tr    = dst%V_tr(:,en) * Ndot0*exp(-StressRatio_s)
+      dot_gamma_tr    = dst%V_tr(:,en) * Ndot0*exp(-StressRatio_s)
-    ddot_gamma_dtau = (dot_gamma_tr*prm%s/tau)*StressRatio_s
+      ddot_gamma_dtau = (dot_gamma_tr*prm%s/tau)*StressRatio_s
-  else where significantStress
+    else where significantStress
-    dot_gamma_tr  = 0.0_pReal
+      dot_gamma_tr  = 0.0_pReal
-    ddot_gamma_dtau = 0.0_pReal
+      ddot_gamma_dtau = 0.0_pReal
-  end where significantStress
+    end where significantStress
 
   end associate
 

src/phase_mechanical_plastic_kinehardening.f90

@@ -359,11 +359,11 @@ module subroutine plastic_kinehardening_results(ph,group)
      case ('xi')
        if(prm%sum_N_sl>0) call results_writeDataset(stt%xi,group,trim(prm%output(o)), &
                                                     'resistance against plastic slip','Pa')
-     case ('tau_b') !ToDo: chi
+     case ('chi')
        if(prm%sum_N_sl>0) call results_writeDataset(stt%chi,group,trim(prm%output(o)), &
                                                     'back stress','Pa')
      case ('sgn(gamma)')
-       if(prm%sum_N_sl>0) call results_writeDataset(stt%sgn_gamma,group,trim(prm%output(o)), & ! ToDo: could be int
+       if(prm%sum_N_sl>0) call results_writeDataset(int(stt%sgn_gamma),group,trim(prm%output(o)), &
                                                     'sense of shear','1')
      case ('chi_0')
        if(prm%sum_N_sl>0) call results_writeDataset(stt%chi_0,group,trim(prm%output(o)), &

src/phase_mechanical_plastic_nonlocal.f90

@@ -1402,8 +1402,8 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
         forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = 0.0_pReal
     elseif (prm%chi_GB >= 0.0_pReal) then
       !* GRAIN BOUNDARY
-      if (any(dNeq(phase_orientation0(ph)%data(en)%asQuaternion(), &
+      if (any(dNeq(phase_O_0(ph)%data(en)%asQuaternion(), &
-                   phase_orientation0(neighbor_phase)%data(neighbor_me)%asQuaternion())) .and. &
+                   phase_O_0(neighbor_phase)%data(neighbor_me)%asQuaternion())) .and. &
           plasticState(neighbor_phase)%nonlocal) &
         forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB)
     else

src/prec.f90

@@ -27,11 +27,6 @@ module prec
 
   real(pReal), parameter :: tol_math_check = 1.0e-8_pReal                                           !< tolerance for internal math self-checks (rotation)
 
-
-  type :: group_float                                                                               !< variable length datatype used for storage of state
-    real(pReal), dimension(:), pointer :: p
-  end type group_float
-
   type :: tState
     integer :: &
       sizeState        = 0, &                                                                       !< size of state
@@ -94,6 +89,7 @@ end subroutine prec_init
 ! replaces "==" but for certain (relative) tolerance. Counterpart to dNeq
 ! https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
 ! AlmostEqualRelative
+! ToDo: Use 'spacing': https://gcc.gnu.org/onlinedocs/gfortran/SPACING.html#SPACING
 !--------------------------------------------------------------------------------------------------
 logical elemental pure function dEq(a,b,tol)
 

src/quit.f90

@@ -6,7 +6,10 @@
 !--------------------------------------------------------------------------------------------------
 subroutine quit(stop_id)
 #include <petsc/finclude/petscsys.h>
-  use PetscSys
+  use PETScSys
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  use MPI_f08
+#endif
   use HDF5
 
   implicit none

src/results.f90

@@ -52,6 +52,7 @@ module results
     results_openGroup, &
     results_closeGroup, &
     results_writeDataset, &
+    results_writeDataset_str, &
     results_setLink, &
     results_addAttribute, &
     results_removeLink, &
@@ -64,16 +65,20 @@ subroutine results_init(restart)
   logical, intent(in) :: restart
 
   character(len=pPathLen) :: commandLine
+  integer :: hdferr
+  integer(HID_T) :: group_id
+  character(len=:), allocatable :: date
+
 
   print'(/,a)', ' <<<+-  results init  -+>>>'; flush(IO_STDOUT)
 
   print*, 'M. Diehl et al., Integrating Materials and Manufacturing Innovation 6(1):83–91, 2017'
   print*, 'https://doi.org/10.1007/s40192-017-0084-5'//IO_EOL
 
-  if(.not. restart) then
+  if (.not. restart) then
     resultsFile = HDF5_openFile(getSolverJobName()//'.hdf5','w')
     call results_addAttribute('DADF5_version_major',0)
-    call results_addAttribute('DADF5_version_minor',13)
+    call results_addAttribute('DADF5_version_minor',14)
     call get_command_argument(0,commandLine)
     call results_addAttribute('creator',trim(commandLine)//' '//DAMASKVERSION)
     call results_addAttribute('created',now())
@@ -81,18 +86,34 @@ subroutine results_init(restart)
     call results_addAttribute('call',trim(commandLine))
     call results_closeGroup(results_addGroup('cell_to'))
     call results_addAttribute('description','mappings to place data in space','cell_to')
-    call results_closeJobFile
+    call results_closeGroup(results_addGroup('setup'))
+    call results_addAttribute('description','input data used to run the simulation','setup')
+  else
+    date = now()
+    call results_openJobFile
+    call get_command(commandLine)
+    call results_addAttribute('call (restart at '//date//')',trim(commandLine))
+    call h5gmove_f(resultsFile,'setup','tmp',hdferr)
+    call results_addAttribute('description','input data used to run the simulation up to restart at '//date,'tmp')
+    call results_closeGroup(results_addGroup('setup'))
+    call results_addAttribute('description','input data used to run the simulation','setup')
+    call h5gmove_f(resultsFile,'tmp','setup/previous',hdferr)
   endif
 
+  call results_closeJobFile
+
 end subroutine results_init
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief opens the results file to append data
 !--------------------------------------------------------------------------------------------------
-subroutine results_openJobFile
+subroutine results_openJobFile(parallel)
 
-  resultsFile = HDF5_openFile(getSolverJobName()//'.hdf5','a')
+  logical, intent(in), optional :: parallel
+
+
+  resultsFile = HDF5_openFile(getSolverJobName()//'.hdf5','a',parallel)
 
 end subroutine results_openJobFile
 
@@ -297,6 +318,25 @@ subroutine results_removeLink(link)
 end subroutine results_removeLink
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief Store string dataset.
+!> @details Not collective, must be called by one process at at time.
+!--------------------------------------------------------------------------------------------------
+subroutine results_writeDataset_str(dataset,group,label,description)
+
+  character(len=*), intent(in) :: label,group,description,dataset
+
+  integer(HID_T) :: groupHandle
+
+
+  groupHandle = results_openGroup(group)
+  call HDF5_write_str(dataset,groupHandle,label)
+  call executionStamp(group//'/'//label,description)
+  call HDF5_closeGroup(groupHandle)
+
+end subroutine results_writeDataset_str
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief Store real scalar dataset with associated metadata.
 !--------------------------------------------------------------------------------------------------