Merge remote-tracking branch 'origin/development' into rotate-Voigt

PRIVATE

@@ -1 +1 @@
-Subproject commit 76bb51348de75207d483d369628670e5ae51dca9
+Subproject commit bc6de828cc4ee9c941b37113ca49fcf51abd3512

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

@@ -6,12 +6,12 @@ references:
 
 output: [xi_sl, xi_tw]
 
-N_sl: [3, 3, 0, 6, 0, 6]   # basal, 1. prism, -, 1. pyr<a>, -, 2. pyr<c+a>
+N_sl: [3, 3, 0, 6, 0, 6]   # basal, prism, -, 1. pyr<a>, -, 2. pyr<c+a>
-N_tw: [6, 0, 6]   # tension, -, compression
+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_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., 60.e+6]
+xi_0_tw:   [40.e+6, 0., 60.e+6]
 
 a_sl: 2.25
 dot_gamma_0_sl: 0.001
@@ -21,9 +21,18 @@ 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-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_sl-sl: [+1.0,  1.0,  1.0,  1.0,  1.0,  1.0, -1.0, -1.0, -1.0, -1.0,
-h_tw-tw: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
+          -1.0, -1.0,  1.0,  1.0, -1.0,  1.0,  1.0, -1.0,  1.0,  1.0,
-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]
+          -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
-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]
+          +1.0,  1.0, -1.0,  1.0, -1.0,  1.0,  1.0, -1.0,  1.0, -1.0,
+          +1.0,  1.0]  # unused entries are indicated by -1.0
+h_tw-tw: [+1.0,  1.0, -1.0, -1.0, -1.0, -1.0,  1.0, -1.0,  1.0,  1.0,
+          -1.0,  1.0]  # unused entries are indicated by -1.0
+h_tw-sl: [+1.0, -1.0,  1.0, -1.0,  1.0, -1.0,  1.0, -1.0, -1.0, -1.0,
+          -1.0, -1.0,  1.0, -1.0,  1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
+          +1.0, -1.0,  1.0, -1.0]  # unused entries are indicated by -1.0
+h_sl-tw: [+1.0, -1.0,  1.0, -1.0,  1.0, -1.0,  1.0, -1.0, -1.0, -1.0,
+          -1.0, -1.0,  1.0, -1.0,  1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
+          +1.0, -1.0,  1.0]  # unused entries are indicated by -1.0

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

@@ -8,7 +8,7 @@ references:
     https://doi.org/10.1016/j.actamat.2017.05.015
 output: [gamma_sl]
 
-N_sl: [3, 3, 0, 0, 12]   # basal, 1. prism, -, -, 2. pyr<c+a>
+N_sl: [3, 3, 0, 0, 12]   # basal, 1. prism, -, -, 1. pyr<c+a>
 n_sl: 20
 a_sl: 2.0
 dot_gamma_0_sl: 0.001
@@ -20,4 +20,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]
+h_sl-sl: [+1.0,  1.0,  1.0,  1.0,  1.0,  1.0, -1.0, -1.0, -1.0, -1.0,
+          -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0,
+          +1.0,  1.0, -1.0, -1.0,  1.0,  1.0, -1.0, -1.0,  1.0,  1.0] # unused entries are indicated by -1.0

python/damask/VERSION

@@ -1 +1 @@
-v3.0.0-alpha5-135-ga3f74994b
+v3.0.0-alpha5-191-gf32a78861

python/damask/_colormap.py

@@ -2,6 +2,8 @@ import os
 import json
 import functools
 import colorsys
+from pathlib import Path
+from typing import Sequence, Union, TextIO
 
 import numpy as np
 import matplotlib as mpl
@@ -14,9 +16,9 @@ from PIL import Image
 from . import util
 from . import Table
 
-_eps   = 216./24389.
+_EPS   = 216./24389.
-_kappa = 24389./27.
+_KAPPA = 24389./27.
-_ref_white = np.array([.95047, 1.00000, 1.08883])                                                   # Observer = 2, Illuminant = D65
+_REF_WHITE = np.array([.95047, 1.00000, 1.08883])                                                   # Observer = 2, Illuminant = D65
 
 # ToDo (if needed)
 # - support alpha channel (paraview/ASCII/input)
@@ -39,20 +41,20 @@ class Colormap(mpl.colors.ListedColormap):
 
     """
 
-    def __add__(self,other):
+    def __add__(self, other: "Colormap") -> "Colormap":
         """Concatenate."""
         return Colormap(np.vstack((self.colors,other.colors)),
                         f'{self.name}+{other.name}')
 
-    def __iadd__(self,other):
+    def __iadd__(self, other: "Colormap") -> "Colormap":
         """Concatenate (in-place)."""
         return self.__add__(other)
 
-    def __invert__(self):
+    def __invert__(self) -> "Colormap":
         """Reverse."""
         return self.reversed()
 
-    def __repr__(self):
+    def __repr__(self) -> str:
         """Show as matplotlib figure."""
         fig = plt.figure(self.name,figsize=(5,.5))
         ax1 = fig.add_axes([0, 0, 1, 1])
@@ -64,7 +66,11 @@ class Colormap(mpl.colors.ListedColormap):
 
 
     @staticmethod
-    def from_range(low,high,name='DAMASK colormap',N=256,model='rgb'):
+    def from_range(low: Sequence[float],
+                   high: Sequence[float],
+                   name: str = 'DAMASK colormap',
+                   N: int = 256,
+                   model: str = 'rgb') -> "Colormap":
         """
         Create a perceptually uniform colormap between given (inclusive) bounds.
 
@@ -145,7 +151,7 @@ class Colormap(mpl.colors.ListedColormap):
 
 
     @staticmethod
-    def from_predefined(name,N=256):
+    def from_predefined(name: str, N: int = 256) -> "Colormap":
         """
         Select from a set of predefined colormaps.
 
@@ -185,7 +191,10 @@ class Colormap(mpl.colors.ListedColormap):
             return Colormap.from_range(definition['low'],definition['high'],name,N)
 
 
-    def shade(self,field,bounds=None,gap=None):
+    def shade(self,
+              field: np.ndarray,
+              bounds: Sequence[float] = None,
+              gap: float = None) -> Image:
         """
         Generate PIL image of 2D field using colormap.
 
@@ -226,7 +235,7 @@ class Colormap(mpl.colors.ListedColormap):
             mode='RGBA')
 
 
-    def reversed(self,name=None):
+    def reversed(self, name: str = None) -> "Colormap":
         """
         Reverse.
 
@@ -251,7 +260,7 @@ class Colormap(mpl.colors.ListedColormap):
         return Colormap(np.array(rev.colors),rev.name[:-4] if rev.name.endswith('_r_r') else rev.name)
 
 
-    def _get_file_handle(self,fname,extension):
+    def _get_file_handle(self, fname: Union[TextIO, str, Path, None], suffix: str) -> TextIO:
         """
         Provide file handle.
 
@@ -259,8 +268,7 @@ class Colormap(mpl.colors.ListedColormap):
         ----------
         fname : file, str, pathlib.Path, or None
             Filename or filehandle, will be name of the colormap+extension if None.
-
+        suffix: str
-        extension: str
             Extension of the filename.
 
         Returns
@@ -270,17 +278,14 @@ class Colormap(mpl.colors.ListedColormap):
 
         """
         if fname is None:
-            fhandle = open(self.name.replace(' ','_')+'.'+extension,'w',newline='\n')
+            return open(self.name.replace(' ','_')+suffix, 'w', newline='\n')
+        elif isinstance(fname, (str, Path)):
+            return open(fname, 'w', newline='\n')
         else:
-            try:
+            return fname
-                fhandle = open(fname,'w',newline='\n')
-            except TypeError:
-                fhandle = fname
-
-        return fhandle
 
 
-    def save_paraview(self,fname=None):
+    def save_paraview(self, fname: Union[TextIO, str, Path] = None):
         """
         Save as JSON file for use in Paraview.
 
@@ -303,10 +308,10 @@ class Colormap(mpl.colors.ListedColormap):
                 'RGBPoints':colors
                }]
 
-        json.dump(out,self._get_file_handle(fname,'json'),indent=4)
+        json.dump(out,self._get_file_handle(fname,'.json'),indent=4)
 
 
-    def save_ASCII(self,fname=None):
+    def save_ASCII(self, fname: Union[TextIO, str, Path] = None):
         """
         Save as ASCII file.
 
@@ -319,10 +324,10 @@ class Colormap(mpl.colors.ListedColormap):
         """
         labels = {'RGBA':4} if self.colors.shape[1] == 4 else {'RGB': 3}
         t = Table(self.colors,labels,f'Creator: {util.execution_stamp("Colormap")}')
-        t.save(self._get_file_handle(fname,'txt'))
+        t.save(self._get_file_handle(fname,'.txt'))
 
 
-    def save_GOM(self,fname=None):
+    def save_GOM(self, fname: Union[TextIO, str, Path] = None):
         """
         Save as ASCII file for use in GOM Aramis.
 
@@ -340,10 +345,10 @@ class Colormap(mpl.colors.ListedColormap):
                 + ' '.join([f' 0 {c[0]} {c[1]} {c[2]} 255 1' for c in reversed((self.colors*255).astype(int))]) \
                 + '\n'
 
-        self._get_file_handle(fname,'legend').write(GOM_str)
+        self._get_file_handle(fname,'.legend').write(GOM_str)
 
 
-    def save_gmsh(self,fname=None):
+    def save_gmsh(self, fname: Union[TextIO, str, Path] = None):
         """
         Save as ASCII file for use in gmsh.
 
@@ -358,11 +363,13 @@ class Colormap(mpl.colors.ListedColormap):
         gmsh_str = 'View.ColorTable = {\n' \
                  +'\n'.join([f'{c[0]},{c[1]},{c[2]},' for c in self.colors[:,:3]*255]) \
                  +'\n}\n'
-        self._get_file_handle(fname,'msh').write(gmsh_str)
+        self._get_file_handle(fname,'.msh').write(gmsh_str)
 
 
     @staticmethod
-    def _interpolate_msh(frac,low,high):
+    def _interpolate_msh(frac,
+                         low: np.ndarray,
+                         high: np.ndarray) -> np.ndarray:
         """
         Interpolate in Msh color space.
 
@@ -439,31 +446,31 @@ class Colormap(mpl.colors.ListedColormap):
 
 
     @staticmethod
-    def _hsv2rgb(hsv):
+    def _hsv2rgb(hsv: np.ndarray) -> np.ndarray:
         """H(ue) S(aturation) V(alue) to R(red) G(reen) B(lue)."""
         return np.array(colorsys.hsv_to_rgb(hsv[0]/360.,hsv[1],hsv[2]))
 
     @staticmethod
-    def _rgb2hsv(rgb):
+    def _rgb2hsv(rgb: np.ndarray) -> np.ndarray:
         """R(ed) G(reen) B(lue) to H(ue) S(aturation) V(alue)."""
         h,s,v = colorsys.rgb_to_hsv(rgb[0],rgb[1],rgb[2])
         return np.array([h*360,s,v])
 
 
     @staticmethod
-    def _hsl2rgb(hsl):
+    def _hsl2rgb(hsl: np.ndarray) -> np.ndarray:
         """H(ue) S(aturation) L(uminance) to R(red) G(reen) B(lue)."""
         return np.array(colorsys.hls_to_rgb(hsl[0]/360.,hsl[2],hsl[1]))
 
     @staticmethod
-    def _rgb2hsl(rgb):
+    def _rgb2hsl(rgb: np.ndarray) -> np.ndarray:
         """R(ed) G(reen) B(lue) to H(ue) S(aturation) L(uminance)."""
         h,l,s = colorsys.rgb_to_hls(rgb[0],rgb[1],rgb[2])
         return np.array([h*360,s,l])
 
 
     @staticmethod
-    def _xyz2rgb(xyz):
+    def _xyz2rgb(xyz: np.ndarray) -> np.ndarray:
         """
         CIE Xyz to R(ed) G(reen) B(lue).
 
@@ -483,7 +490,7 @@ class Colormap(mpl.colors.ListedColormap):
         return np.clip(rgb,0.,1.)
 
     @staticmethod
-    def _rgb2xyz(rgb):
+    def _rgb2xyz(rgb: np.ndarray) -> np.ndarray:
         """
         R(ed) G(reen) B(lue) to CIE Xyz.
 
@@ -501,7 +508,7 @@ class Colormap(mpl.colors.ListedColormap):
 
 
     @staticmethod
-    def _lab2xyz(lab,ref_white=None):
+    def _lab2xyz(lab: np.ndarray, ref_white: np.ndarray = None) -> np.ndarray:
         """
         CIE Lab to CIE Xyz.
 
@@ -514,13 +521,13 @@ class Colormap(mpl.colors.ListedColormap):
         f_z = (lab[0]+16.)/116. - lab[2]/200.
 
         return np.array([
-                         f_x**3.                if f_x**3. > _eps     else (116.*f_x-16.)/_kappa,
+                         f_x**3.                if f_x**3. > _EPS     else (116.*f_x-16.)/_KAPPA,
-                         ((lab[0]+16.)/116.)**3 if lab[0]>_kappa*_eps else lab[0]/_kappa,
+                         ((lab[0]+16.)/116.)**3 if lab[0]>_KAPPA*_EPS else lab[0]/_KAPPA,
-                         f_z**3.                if f_z**3. > _eps     else (116.*f_z-16.)/_kappa
+                         f_z**3.                if f_z**3. > _EPS     else (116.*f_z-16.)/_KAPPA
-                        ])*(ref_white if ref_white is not None else _ref_white)
+                        ])*(ref_white if ref_white is not None else _REF_WHITE)
 
     @staticmethod
-    def _xyz2lab(xyz,ref_white=None):
+    def _xyz2lab(xyz: np.ndarray, ref_white: np.ndarray = None) -> np.ndarray:
         """
         CIE Xyz to CIE Lab.
 
@@ -529,8 +536,8 @@ class Colormap(mpl.colors.ListedColormap):
         http://www.brucelindbloom.com/index.html?Eqn_Lab_to_XYZ.html
 
         """
-        ref_white = ref_white if ref_white is not None else _ref_white
+        ref_white = ref_white if ref_white is not None else _REF_WHITE
-        f = np.where(xyz/ref_white > _eps,(xyz/ref_white)**(1./3.),(_kappa*xyz/ref_white+16.)/116.)
+        f = np.where(xyz/ref_white > _EPS,(xyz/ref_white)**(1./3.),(_KAPPA*xyz/ref_white+16.)/116.)
 
         return np.array([
                          116.0 *  f[1] - 16.0,
@@ -540,7 +547,7 @@ class Colormap(mpl.colors.ListedColormap):
 
 
     @staticmethod
-    def _lab2msh(lab):
+    def _lab2msh(lab: np.ndarray) -> np.ndarray:
         """
         CIE Lab to Msh.
 
@@ -558,7 +565,7 @@ class Colormap(mpl.colors.ListedColormap):
                         ])
 
     @staticmethod
-    def _msh2lab(msh):
+    def _msh2lab(msh: np.ndarray) -> np.ndarray:
         """
         Msh to CIE Lab.
 
@@ -575,29 +582,29 @@ class Colormap(mpl.colors.ListedColormap):
                         ])
 
     @staticmethod
-    def _lab2rgb(lab):
+    def _lab2rgb(lab: np.ndarray) -> np.ndarray:
         return Colormap._xyz2rgb(Colormap._lab2xyz(lab))
 
     @staticmethod
-    def _rgb2lab(rgb):
+    def _rgb2lab(rgb: np.ndarray) -> np.ndarray:
         return Colormap._xyz2lab(Colormap._rgb2xyz(rgb))
 
     @staticmethod
-    def _msh2rgb(msh):
+    def _msh2rgb(msh: np.ndarray) -> np.ndarray:
         return Colormap._lab2rgb(Colormap._msh2lab(msh))
 
     @staticmethod
-    def _rgb2msh(rgb):
+    def _rgb2msh(rgb: np.ndarray) -> np.ndarray:
         return Colormap._lab2msh(Colormap._rgb2lab(rgb))
 
     @staticmethod
-    def _hsv2msh(hsv):
+    def _hsv2msh(hsv: np.ndarray) -> np.ndarray:
         return Colormap._rgb2msh(Colormap._hsv2rgb(hsv))
 
     @staticmethod
-    def _hsl2msh(hsl):
+    def _hsl2msh(hsl: np.ndarray) -> np.ndarray:
         return Colormap._rgb2msh(Colormap._hsl2rgb(hsl))
 
     @staticmethod
-    def _xyz2msh(xyz):
+    def _xyz2msh(xyz: np.ndarray) -> np.ndarray:
         return Colormap._lab2msh(Colormap._xyz2lab(xyz))

python/tests/test_Colormap.py

@@ -77,6 +77,12 @@ class TestColormap:
             # xyz2msh
             assert np.allclose(Colormap._xyz2msh(xyz),msh,atol=1.e-6,rtol=0)
 
+    @pytest.mark.parametrize('low,high',[((0,0,0),(1,1,1)),
+                                         ([0,0,0],[1,1,1])])
+    def test_from_range_types(self,low,high):
+        a = Colormap.from_range(low,high)
+        b = Colormap.from_range(np.array(low),np.array(high))
+        assert np.all(a.colors == b.colors) 
 
     @pytest.mark.parametrize('format',['ASCII','paraview','GOM','gmsh'])
     @pytest.mark.parametrize('model',['rgb','hsv','hsl','xyz','lab','msh'])

src/commercialFEM_fileList.f90

@@ -4,6 +4,7 @@
 !> @details List of files needed by MSC.Marc
 !--------------------------------------------------------------------------------------------------
 #include "parallelization.f90"
+#include "constants.f90"
 #include "IO.f90"
 #include "YAML_types.f90"
 #include "YAML_parse.f90"

src/constants.f90

@@ -0,0 +1,15 @@
+!--------------------------------------------------------------------------------------------------
+!> @author Martin Diehl, KU Leuven
+!> @brief  physical constants
+!--------------------------------------------------------------------------------------------------
+module constants
+  use prec
+
+  implicit none
+  public
+
+  real(pReal), parameter :: &
+    T_ROOM = 300.0_pReal, &                                                                         !< Room temperature in K
+    K_B = 1.38e-23_pReal                                                                             !< Boltzmann constant in J/Kelvin
+
+end module constants

src/lattice.f90

@@ -219,18 +219,18 @@ module lattice
        2, -1, -1,  0,     0,  1, -1,  0, &
       -1,  2, -1,  0,    -1,  0,  1,  0, &
       -1, -1,  2,  0,     1, -1,  0,  0, &
-    ! <-11.0>{11.0}/2nd order prismatic compound systems (plane normal independent of c/a-ratio)
+    ! <-11.0>{11.0}/2. order prismatic compound systems (plane normal independent of c/a-ratio)
       -1,  1,  0,  0,     1,  1, -2,  0, &
        0, -1,  1,  0,    -2,  1,  1,  0, &
        1,  0, -1,  0,     1, -2,  1,  0, &
-    ! <-1-1.0>{-11.1}/1st order pyramidal <a> systems (direction independent of c/a-ratio)
+    ! <-1-1.0>{-11.1}/1. order pyramidal <a> systems (direction independent of c/a-ratio)
       -1,  2, -1,  0,     1,  0, -1,  1, &
       -2,  1,  1,  0,     0,  1, -1,  1, &
       -1, -1,  2,  0,    -1,  1,  0,  1, &
        1, -2,  1,  0,    -1,  0,  1,  1, &
        2, -1, -1,  0,     0, -1,  1,  1, &
        1,  1, -2,  0,     1, -1,  0,  1, &
-    ! <11.3>{-10.1}/1st order pyramidal <c+a> systems (direction independent of c/a-ratio)
+    ! <11.3>{-10.1}/1. order pyramidal <c+a> systems (direction independent of c/a-ratio)
       -2,  1,  1,  3,     1,  0, -1,  1, &
       -1, -1,  2,  3,     1,  0, -1,  1, &
       -1, -1,  2,  3,     0,  1, -1,  1, &
@@ -243,7 +243,7 @@ module lattice
       -1,  2, -1,  3,     0, -1,  1,  1, &
       -1,  2, -1,  3,     1, -1,  0,  1, &
       -2,  1,  1,  3,     1, -1,  0,  1, &
-    ! <11.3>{-1-1.2}/2nd order pyramidal <c+a> systems
+    ! <11.3>{-1-1.2}/2. order pyramidal <c+a> systems
       -1, -1,  2,  3,     1,  1, -2,  2, &
        1, -2,  1,  3,    -1,  2, -1,  2, &
        2, -1, -1,  3,    -2,  1,  1,  2, &
@@ -491,7 +491,7 @@ end function lattice_characteristicShear_Twin
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Rotated elasticity matrices for twinning in 66-vector notation
+!> @brief Rotated elasticity matrices for twinning in 6x6-matrix notation
 !--------------------------------------------------------------------------------------------------
 function lattice_C66_twin(Ntwin,C66,lattice,CoverA)
 
@@ -529,7 +529,7 @@ end function lattice_C66_twin
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Rotated elasticity matrices for transformation in 66-vector notation
+!> @brief Rotated elasticity matrices for transformation in 6x6-matrix notation
 !--------------------------------------------------------------------------------------------------
 function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
                            cOverA_trans,a_bcc,a_fcc)
@@ -753,22 +753,23 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(
 
   integer, dimension(HEX_NSLIP,HEX_NSLIP), parameter :: &
     HEX_INTERACTIONSLIPSLIP = reshape( [&
+    ! basal      prism      2. prism   1. pyr<a>           1. pyr<c+a>                           2. pyr<c+a>
        1, 2, 2,   3, 3, 3,   7, 7, 7,  13,13,13,13,13,13,  21,21,21,21,21,21,21,21,21,21,21,21,  31,31,31,31,31,31, & ! -----> acting (forest)
-       2, 1, 2,   3, 3, 3,   7, 7, 7,  13,13,13,13,13,13,  21,21,21,21,21,21,21,21,21,21,21,21,  31,31,31,31,31,31, & ! |
+       2, 1, 2,   3, 3, 3,   7, 7, 7,  13,13,13,13,13,13,  21,21,21,21,21,21,21,21,21,21,21,21,  31,31,31,31,31,31, & ! | basal
        2, 2, 1,   3, 3, 3,   7, 7, 7,  13,13,13,13,13,13,  21,21,21,21,21,21,21,21,21,21,21,21,  31,31,31,31,31,31, & ! |
                                                                                                                       ! v
        6, 6, 6,   4, 5, 5,   8, 8, 8,  14,14,14,14,14,14,  22,22,22,22,22,22,22,22,22,22,22,22,  32,32,32,32,32,32, & ! reacting (primary)
-       6, 6, 6,   5, 4, 5,   8, 8, 8,  14,14,14,14,14,14,  22,22,22,22,22,22,22,22,22,22,22,22,  32,32,32,32,32,32, &
+       6, 6, 6,   5, 4, 5,   8, 8, 8,  14,14,14,14,14,14,  22,22,22,22,22,22,22,22,22,22,22,22,  32,32,32,32,32,32, & ! prism
        6, 6, 6,   5, 5, 4,   8, 8, 8,  14,14,14,14,14,14,  22,22,22,22,22,22,22,22,22,22,22,22,  32,32,32,32,32,32, &
 
       12,12,12,  11,11,11,   9,10,10,  15,15,15,15,15,15,  23,23,23,23,23,23,23,23,23,23,23,23,  33,33,33,33,33,33, &
-      12,12,12,  11,11,11,  10, 9,10,  15,15,15,15,15,15,  23,23,23,23,23,23,23,23,23,23,23,23,  33,33,33,33,33,33, &
+      12,12,12,  11,11,11,  10, 9,10,  15,15,15,15,15,15,  23,23,23,23,23,23,23,23,23,23,23,23,  33,33,33,33,33,33, & ! 2. prism
       12,12,12,  11,11,11,  10,10, 9,  15,15,15,15,15,15,  23,23,23,23,23,23,23,23,23,23,23,23,  33,33,33,33,33,33, &
 
       20,20,20,  19,19,19,  18,18,18,  16,17,17,17,17,17,  24,24,24,24,24,24,24,24,24,24,24,24,  34,34,34,34,34,34, &
       20,20,20,  19,19,19,  18,18,18,  17,16,17,17,17,17,  24,24,24,24,24,24,24,24,24,24,24,24,  34,34,34,34,34,34, &
       20,20,20,  19,19,19,  18,18,18,  17,17,16,17,17,17,  24,24,24,24,24,24,24,24,24,24,24,24,  34,34,34,34,34,34, &
-      20,20,20,  19,19,19,  18,18,18,  17,17,17,16,17,17,  24,24,24,24,24,24,24,24,24,24,24,24,  34,34,34,34,34,34, &
+      20,20,20,  19,19,19,  18,18,18,  17,17,17,16,17,17,  24,24,24,24,24,24,24,24,24,24,24,24,  34,34,34,34,34,34, & ! 1. pyr<a>
       20,20,20,  19,19,19,  18,18,18,  17,17,17,17,16,17,  24,24,24,24,24,24,24,24,24,24,24,24,  34,34,34,34,34,34, &
       20,20,20,  19,19,19,  18,18,18,  17,17,17,17,17,16,  24,24,24,24,24,24,24,24,24,24,24,24,  34,34,34,34,34,34, &
 
@@ -778,7 +779,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(
       30,30,30,  29,29,29,  28,28,28,  27,27,27,27,27,27,  26,26,26,25,26,26,26,26,26,26,26,26,  35,35,35,35,35,35, &
       30,30,30,  29,29,29,  28,28,28,  27,27,27,27,27,27,  26,26,26,26,25,26,26,26,26,26,26,26,  35,35,35,35,35,35, &
       30,30,30,  29,29,29,  28,28,28,  27,27,27,27,27,27,  26,26,26,26,26,25,26,26,26,26,26,26,  35,35,35,35,35,35, &
-      30,30,30,  29,29,29,  28,28,28,  27,27,27,27,27,27,  26,26,26,26,26,26,25,26,26,26,26,26,  35,35,35,35,35,35, &
+      30,30,30,  29,29,29,  28,28,28,  27,27,27,27,27,27,  26,26,26,26,26,26,25,26,26,26,26,26,  35,35,35,35,35,35, & ! 1. pyr<c+a>
       30,30,30,  29,29,29,  28,28,28,  27,27,27,27,27,27,  26,26,26,26,26,26,26,25,26,26,26,26,  35,35,35,35,35,35, &
       30,30,30,  29,29,29,  28,28,28,  27,27,27,27,27,27,  26,26,26,26,26,26,26,26,25,26,26,26,  35,35,35,35,35,35, &
       30,30,30,  29,29,29,  28,28,28,  27,27,27,27,27,27,  26,26,26,26,26,26,26,26,26,25,26,26,  35,35,35,35,35,35, &
@@ -788,7 +789,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(
       42,42,42,  41,41,41,  40,40,40,  39,39,39,39,39,39,  38,38,38,38,38,38,38,38,38,38,38,38,  36,37,37,37,37,37, &
       42,42,42,  41,41,41,  40,40,40,  39,39,39,39,39,39,  38,38,38,38,38,38,38,38,38,38,38,38,  37,36,37,37,37,37, &
       42,42,42,  41,41,41,  40,40,40,  39,39,39,39,39,39,  38,38,38,38,38,38,38,38,38,38,38,38,  37,37,36,37,37,37, &
-      42,42,42,  41,41,41,  40,40,40,  39,39,39,39,39,39,  38,38,38,38,38,38,38,38,38,38,38,38,  37,37,37,36,37,37, &
+      42,42,42,  41,41,41,  40,40,40,  39,39,39,39,39,39,  38,38,38,38,38,38,38,38,38,38,38,38,  37,37,37,36,37,37, & ! 2. pyr<c+a>
       42,42,42,  41,41,41,  40,40,40,  39,39,39,39,39,39,  38,38,38,38,38,38,38,38,38,38,38,38,  37,37,37,37,36,37, &
       42,42,42,  41,41,41,  40,40,40,  39,39,39,39,39,39,  38,38,38,38,38,38,38,38,38,38,38,38,  37,37,37,37,37,36  &
       ],shape(HEX_INTERACTIONSLIPSLIP))                                                             !< Slip-slip interaction types for hex (onion peel naming scheme)
@@ -934,31 +935,32 @@ function lattice_interaction_TwinByTwin(Ntwin,interactionValues,lattice) result(
                                                                                                     !< 3: other interaction
   integer, dimension(HEX_NTWIN,HEX_NTWIN), parameter :: &
     HEX_INTERACTIONTWINTWIN = reshape( [&
+    ! <-10.1>{10.2}       <11.6>{-1-1.1}      <10.-2>{10.1}       <11.-3>{11.2}
        1, 2, 2, 2, 2, 2,   3, 3, 3, 3, 3, 3,   7, 7, 7, 7, 7, 7,  13,13,13,13,13,13, & ! -----> acting
        2, 1, 2, 2, 2, 2,   3, 3, 3, 3, 3, 3,   7, 7, 7, 7, 7, 7,  13,13,13,13,13,13, & ! |
        2, 2, 1, 2, 2, 2,   3, 3, 3, 3, 3, 3,   7, 7, 7, 7, 7, 7,  13,13,13,13,13,13, & ! |
-       2, 2, 2, 1, 2, 2,   3, 3, 3, 3, 3, 3,   7, 7, 7, 7, 7, 7,  13,13,13,13,13,13, & ! v
+       2, 2, 2, 1, 2, 2,   3, 3, 3, 3, 3, 3,   7, 7, 7, 7, 7, 7,  13,13,13,13,13,13, & ! v <-10.1>{10.2}
        2, 2, 2, 2, 1, 2,   3, 3, 3, 3, 3, 3,   7, 7, 7, 7, 7, 7,  13,13,13,13,13,13, & ! reacting
        2, 2, 2, 2, 2, 1,   3, 3, 3, 3, 3, 3,   7, 7, 7, 7, 7, 7,  13,13,13,13,13,13, &
 
        6, 6, 6, 6, 6, 6,   4, 5, 5, 5, 5, 5,   8, 8, 8, 8, 8, 8,  14,14,14,14,14,14, &
        6, 6, 6, 6, 6, 6,   5, 4, 5, 5, 5, 5,   8, 8, 8, 8, 8, 8,  14,14,14,14,14,14, &
        6, 6, 6, 6, 6, 6,   5, 5, 4, 5, 5, 5,   8, 8, 8, 8, 8, 8,  14,14,14,14,14,14, &
-       6, 6, 6, 6, 6, 6,   5, 5, 5, 4, 5, 5,   8, 8, 8, 8, 8, 8,  14,14,14,14,14,14, &
+       6, 6, 6, 6, 6, 6,   5, 5, 5, 4, 5, 5,   8, 8, 8, 8, 8, 8,  14,14,14,14,14,14, & ! <11.6>{-1-1.1}
        6, 6, 6, 6, 6, 6,   5, 5, 5, 5, 4, 5,   8, 8, 8, 8, 8, 8,  14,14,14,14,14,14, &
        6, 6, 6, 6, 6, 6,   5, 5, 5, 5, 5, 4,   8, 8, 8, 8, 8, 8,  14,14,14,14,14,14, &
 
       12,12,12,12,12,12,  11,11,11,11,11,11,   9,10,10,10,10,10,  15,15,15,15,15,15, &
       12,12,12,12,12,12,  11,11,11,11,11,11,  10, 9,10,10,10,10,  15,15,15,15,15,15, &
       12,12,12,12,12,12,  11,11,11,11,11,11,  10,10, 9,10,10,10,  15,15,15,15,15,15, &
-      12,12,12,12,12,12,  11,11,11,11,11,11,  10,10,10, 9,10,10,  15,15,15,15,15,15, &
+      12,12,12,12,12,12,  11,11,11,11,11,11,  10,10,10, 9,10,10,  15,15,15,15,15,15, & ! <10.-2>{10.1}
       12,12,12,12,12,12,  11,11,11,11,11,11,  10,10,10,10, 9,10,  15,15,15,15,15,15, &
       12,12,12,12,12,12,  11,11,11,11,11,11,  10,10,10,10,10, 9,  15,15,15,15,15,15, &
 
       20,20,20,20,20,20,  19,19,19,19,19,19,  18,18,18,18,18,18,  16,17,17,17,17,17, &
       20,20,20,20,20,20,  19,19,19,19,19,19,  18,18,18,18,18,18,  17,16,17,17,17,17, &
       20,20,20,20,20,20,  19,19,19,19,19,19,  18,18,18,18,18,18,  17,17,16,17,17,17, &
-      20,20,20,20,20,20,  19,19,19,19,19,19,  18,18,18,18,18,18,  17,17,17,16,17,17, &
+      20,20,20,20,20,20,  19,19,19,19,19,19,  18,18,18,18,18,18,  17,17,17,16,17,17, & ! <11.-3>{11.2}
       20,20,20,20,20,20,  19,19,19,19,19,19,  18,18,18,18,18,18,  17,17,17,17,16,17, &
       20,20,20,20,20,20,  19,19,19,19,19,19,  18,18,18,18,18,18,  17,17,17,17,17,16  &
       ],shape(HEX_INTERACTIONTWINTWIN))                                                             !< Twin-twin interaction types for hex
@@ -1012,7 +1014,7 @@ function lattice_interaction_TransByTrans(Ntrans,interactionValues,lattice) resu
       2,2,2,2,2,2,2,2,2,1,1,1  &
       ],shape(FCC_INTERACTIONTRANSTRANS))                                                           !< Trans-trans interaction types for fcc
 
-  if(lattice == 'cF') then
+  if (lattice == 'cF') then
     interactionTypes = FCC_INTERACTIONTRANSTRANS
     NtransMax        = FCC_NTRANSSYSTEM
   else
@@ -1125,22 +1127,23 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) r
 
   integer, dimension(HEX_NTWIN,HEX_NSLIP), parameter :: &
     HEX_INTERACTIONSLIPTWIN = reshape( [&
+    ! <-10.1>{10.2}       <11.6>{-1-1.1}      <10.-2>{10.1}       <11.-3>{11.2}
        1, 1, 1, 1, 1, 1,   2, 2, 2, 2, 2, 2,   3, 3, 3, 3, 3, 3,   4, 4, 4, 4, 4, 4, & ! ----> twin (acting)
-       1, 1, 1, 1, 1, 1,   2, 2, 2, 2, 2, 2,   3, 3, 3, 3, 3, 3,   4, 4, 4, 4, 4, 4, & ! |
+       1, 1, 1, 1, 1, 1,   2, 2, 2, 2, 2, 2,   3, 3, 3, 3, 3, 3,   4, 4, 4, 4, 4, 4, & ! | basal
        1, 1, 1, 1, 1, 1,   2, 2, 2, 2, 2, 2,   3, 3, 3, 3, 3, 3,   4, 4, 4, 4, 4, 4, & ! |
                                                                                        ! v
        5, 5, 5, 5, 5, 5,   6, 6, 6, 6, 6, 6,   7, 7, 7, 7, 7, 7,   8, 8, 8, 8, 8, 8, & ! slip (reacting)
-       5, 5, 5, 5, 5, 5,   6, 6, 6, 6, 6, 6,   7, 7, 7, 7, 7, 7,   8, 8, 8, 8, 8, 8, &
+       5, 5, 5, 5, 5, 5,   6, 6, 6, 6, 6, 6,   7, 7, 7, 7, 7, 7,   8, 8, 8, 8, 8, 8, & ! prism
        5, 5, 5, 5, 5, 5,   6, 6, 6, 6, 6, 6,   7, 7, 7, 7, 7, 7,   8, 8, 8, 8, 8, 8, &
 
        9, 9, 9, 9, 9, 9,  10,10,10,10,10,10,  11,11,11,11,11,11,  12,12,12,12,12,12, &
-       9, 9, 9, 9, 9, 9,  10,10,10,10,10,10,  11,11,11,11,11,11,  12,12,12,12,12,12, &
+       9, 9, 9, 9, 9, 9,  10,10,10,10,10,10,  11,11,11,11,11,11,  12,12,12,12,12,12, & ! 2.prism
        9, 9, 9, 9, 9, 9,  10,10,10,10,10,10,  11,11,11,11,11,11,  12,12,12,12,12,12, &
 
       13,13,13,13,13,13,  14,14,14,14,14,14,  15,15,15,15,15,15,  16,16,16,16,16,16, &
       13,13,13,13,13,13,  14,14,14,14,14,14,  15,15,15,15,15,15,  16,16,16,16,16,16, &
       13,13,13,13,13,13,  14,14,14,14,14,14,  15,15,15,15,15,15,  16,16,16,16,16,16, &
-      13,13,13,13,13,13,  14,14,14,14,14,14,  15,15,15,15,15,15,  16,16,16,16,16,16, &
+      13,13,13,13,13,13,  14,14,14,14,14,14,  15,15,15,15,15,15,  16,16,16,16,16,16, & ! 1. pyr<a>
       13,13,13,13,13,13,  14,14,14,14,14,14,  15,15,15,15,15,15,  16,16,16,16,16,16, &
       13,13,13,13,13,13,  14,14,14,14,14,14,  15,15,15,15,15,15,  16,16,16,16,16,16, &
 
@@ -1150,7 +1153,7 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) r
       17,17,17,17,17,17,  18,18,18,18,18,18,  19,19,19,19,19,19,  20,20,20,20,20,20, &
       17,17,17,17,17,17,  18,18,18,18,18,18,  19,19,19,19,19,19,  20,20,20,20,20,20, &
       17,17,17,17,17,17,  18,18,18,18,18,18,  19,19,19,19,19,19,  20,20,20,20,20,20, &
-      17,17,17,17,17,17,  18,18,18,18,18,18,  19,19,19,19,19,19,  20,20,20,20,20,20, &
+      17,17,17,17,17,17,  18,18,18,18,18,18,  19,19,19,19,19,19,  20,20,20,20,20,20, & ! 1. pyr<c+a>
       17,17,17,17,17,17,  18,18,18,18,18,18,  19,19,19,19,19,19,  20,20,20,20,20,20, &
       17,17,17,17,17,17,  18,18,18,18,18,18,  19,19,19,19,19,19,  20,20,20,20,20,20, &
       17,17,17,17,17,17,  18,18,18,18,18,18,  19,19,19,19,19,19,  20,20,20,20,20,20, &
@@ -1160,7 +1163,7 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) r
       21,21,21,21,21,21,  22,22,22,22,22,22,  23,23,23,23,23,23,  24,24,24,24,24,24, &
       21,21,21,21,21,21,  22,22,22,22,22,22,  23,23,23,23,23,23,  24,24,24,24,24,24, &
       21,21,21,21,21,21,  22,22,22,22,22,22,  23,23,23,23,23,23,  24,24,24,24,24,24, &
-      21,21,21,21,21,21,  22,22,22,22,22,22,  23,23,23,23,23,23,  24,24,24,24,24,24, &
+      21,21,21,21,21,21,  22,22,22,22,22,22,  23,23,23,23,23,23,  24,24,24,24,24,24, & ! 2. pyr<c+a>
       21,21,21,21,21,21,  22,22,22,22,22,22,  23,23,23,23,23,23,  24,24,24,24,24,24, &
       21,21,21,21,21,21,  22,22,22,22,22,22,  23,23,23,23,23,23,  24,24,24,24,24,24  &
       ],shape(HEX_INTERACTIONSLIPTWIN))                                                             !< Slip-twin interaction types for hex
@@ -1264,31 +1267,32 @@ function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,lattice) r
 
   integer, dimension(HEX_NSLIP,HEX_NTWIN), parameter :: &
     HEX_INTERACTIONTWINSLIP = reshape( [&
+    ! basal      prism     2. prism   1. pyr<a>           1. pyr<c+a>                           2. pyr<c+a>
       1, 1, 1,   5, 5, 5,   9, 9, 9,  13,13,13,13,13,13,  17,17,17,17,17,17,17,17,17,17,17,17,  21,21,21,21,21,21, & ! ----> slip (acting)
       1, 1, 1,   5, 5, 5,   9, 9, 9,  13,13,13,13,13,13,  17,17,17,17,17,17,17,17,17,17,17,17,  21,21,21,21,21,21, & ! |
       1, 1, 1,   5, 5, 5,   9, 9, 9,  13,13,13,13,13,13,  17,17,17,17,17,17,17,17,17,17,17,17,  21,21,21,21,21,21, & ! |
-      1, 1, 1,   5, 5, 5,   9, 9, 9,  13,13,13,13,13,13,  17,17,17,17,17,17,17,17,17,17,17,17,  21,21,21,21,21,21, & ! v
+      1, 1, 1,   5, 5, 5,   9, 9, 9,  13,13,13,13,13,13,  17,17,17,17,17,17,17,17,17,17,17,17,  21,21,21,21,21,21, & ! v <-10.1>{10.2}
       1, 1, 1,   5, 5, 5,   9, 9, 9,  13,13,13,13,13,13,  17,17,17,17,17,17,17,17,17,17,17,17,  21,21,21,21,21,21, & ! twin (reacting)
       1, 1, 1,   5, 5, 5,   9, 9, 9,  13,13,13,13,13,13,  17,17,17,17,17,17,17,17,17,17,17,17,  21,21,21,21,21,21, &
 
       2, 2, 2,   6, 6, 6,  10,10,10,  14,14,14,14,14,14,  18,18,18,18,18,18,18,18,18,18,18,18,  22,22,22,22,22,22, &
       2, 2, 2,   6, 6, 6,  10,10,10,  14,14,14,14,14,14,  18,18,18,18,18,18,18,18,18,18,18,18,  22,22,22,22,22,22, &
       2, 2, 2,   6, 6, 6,  10,10,10,  14,14,14,14,14,14,  18,18,18,18,18,18,18,18,18,18,18,18,  22,22,22,22,22,22, &
-      2, 2, 2,   6, 6, 6,  10,10,10,  14,14,14,14,14,14,  18,18,18,18,18,18,18,18,18,18,18,18,  22,22,22,22,22,22, &
+      2, 2, 2,   6, 6, 6,  10,10,10,  14,14,14,14,14,14,  18,18,18,18,18,18,18,18,18,18,18,18,  22,22,22,22,22,22, & ! <11.6>{-1-1.1}
       2, 2, 2,   6, 6, 6,  10,10,10,  14,14,14,14,14,14,  18,18,18,18,18,18,18,18,18,18,18,18,  22,22,22,22,22,22, &
       2, 2, 2,   6, 6, 6,  10,10,10,  14,14,14,14,14,14,  18,18,18,18,18,18,18,18,18,18,18,18,  22,22,22,22,22,22, &
 
       3, 3, 3,   7, 7, 7,  11,11,11,  15,15,15,15,15,15,  19,19,19,19,19,19,19,19,19,19,19,19,  23,23,23,23,23,23, &
       3, 3, 3,   7, 7, 7,  11,11,11,  15,15,15,15,15,15,  19,19,19,19,19,19,19,19,19,19,19,19,  23,23,23,23,23,23, &
       3, 3, 3,   7, 7, 7,  11,11,11,  15,15,15,15,15,15,  19,19,19,19,19,19,19,19,19,19,19,19,  23,23,23,23,23,23, &
-      3, 3, 3,   7, 7, 7,  11,11,11,  15,15,15,15,15,15,  19,19,19,19,19,19,19,19,19,19,19,19,  23,23,23,23,23,23, &
+      3, 3, 3,   7, 7, 7,  11,11,11,  15,15,15,15,15,15,  19,19,19,19,19,19,19,19,19,19,19,19,  23,23,23,23,23,23, & ! <10.-2>{10.1}
       3, 3, 3,   7, 7, 7,  11,11,11,  15,15,15,15,15,15,  19,19,19,19,19,19,19,19,19,19,19,19,  23,23,23,23,23,23, &
       3, 3, 3,   7, 7, 7,  11,11,11,  15,15,15,15,15,15,  19,19,19,19,19,19,19,19,19,19,19,19,  23,23,23,23,23,23, &
 
       4, 4, 4,   8, 8, 8,  12,12,12,  16,16,16,16,16,16,  20,20,20,20,20,20,20,20,20,20,20,20,  24,24,24,24,24,24, &
       4, 4, 4,   8, 8, 8,  12,12,12,  16,16,16,16,16,16,  20,20,20,20,20,20,20,20,20,20,20,20,  24,24,24,24,24,24, &
       4, 4, 4,   8, 8, 8,  12,12,12,  16,16,16,16,16,16,  20,20,20,20,20,20,20,20,20,20,20,20,  24,24,24,24,24,24, &
-      4, 4, 4,   8, 8, 8,  12,12,12,  16,16,16,16,16,16,  20,20,20,20,20,20,20,20,20,20,20,20,  24,24,24,24,24,24, &
+      4, 4, 4,   8, 8, 8,  12,12,12,  16,16,16,16,16,16,  20,20,20,20,20,20,20,20,20,20,20,20,  24,24,24,24,24,24, & ! <11.-3>{11.2}
       4, 4, 4,   8, 8, 8,  12,12,12,  16,16,16,16,16,16,  20,20,20,20,20,20,20,20,20,20,20,20,  24,24,24,24,24,24, &
       4, 4, 4,   8, 8, 8,  12,12,12,  16,16,16,16,16,16,  20,20,20,20,20,20,20,20,20,20,20,20,  24,24,24,24,24,24  &
       ],shape(HEX_INTERACTIONTWINSLIP))                                                             !< Twin-slip interaction types for hex
@@ -1413,7 +1417,7 @@ end function lattice_SchmidMatrix_twin
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Schmid matrix for twinning
+!> @brief Schmid matrix for transformation
 !> details only active twin systems are considered
 !--------------------------------------------------------------------------------------------------
 function lattice_SchmidMatrix_trans(Ntrans,lattice_target,cOverA,a_bcc,a_fcc) result(SchmidMatrix)
@@ -1539,7 +1543,7 @@ end function lattice_slip_transverse
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Labels for slip systems
+!> @brief Labels of slip systems
 !> details only active slip systems are considered
 !--------------------------------------------------------------------------------------------------
 function lattice_labels_slip(Nslip,lattice) result(labels)
@@ -2089,12 +2093,12 @@ function lattice_equivalent_nu(C,assumption) result(nu)
   real(pReal), dimension(6,6) :: S
 
 
-  if    (IO_lc(assumption) == 'voigt') then
+  if     (IO_lc(assumption) == 'voigt') then
     K = (C(1,1)+C(2,2)+C(3,3) +2.0_pReal*(C(1,2)+C(2,3)+C(1,3))) &
       / 9.0_pReal
-  elseif(IO_lc(assumption) == 'reuss') then
+  elseif (IO_lc(assumption) == 'reuss') then
-    call math_invert(S,error,C)                                                                     ! ToDo: correct for Voigt?
+    call math_invert(S,error,C)
-    if(error) error stop 'matrix inversion failed'
+    if (error) error stop 'matrix inversion failed'
     K = 1.0_pReal &
       / (S(1,1)+S(2,2)+S(3,3) +2.0_pReal*(S(1,2)+S(2,3)+S(1,3)))
   else
@@ -2102,7 +2106,7 @@ function lattice_equivalent_nu(C,assumption) result(nu)
   endif
 
   mu = lattice_equivalent_mu(C,assumption)
-  nu = (1.5_pReal*K -mu)/(3.0_pReal*K+mu)
+  nu = (1.5_pReal*K-mu)/(3.0_pReal*K+mu)
 
 end function lattice_equivalent_nu
 
@@ -2121,12 +2125,12 @@ function lattice_equivalent_mu(C,assumption) result(mu)
   real(pReal), dimension(6,6) :: S
 
 
-  if    (IO_lc(assumption) == 'voigt') then
+  if     (IO_lc(assumption) == 'voigt') then
     mu = (1.0_pReal*(C(1,1)+C(2,2)+C(3,3)) -1.0_pReal*(C(1,2)+C(2,3)+C(1,3)) +3.0_pReal*(C(4,4)+C(5,5)+C(6,6))) &
        / 15.0_pReal
-  elseif(IO_lc(assumption) == 'reuss') then
+  elseif (IO_lc(assumption) == 'reuss') then
-    call math_invert(S,error,C)                                                                     ! ToDo: correct for Voigt?
+    call math_invert(S,error,C)
-    if(error) error stop 'matrix inversion failed'
+    if (error) error stop 'matrix inversion failed'
     mu = 15.0_pReal &
        / (4.0_pReal*(S(1,1)+S(2,2)+S(3,3)) -4.0_pReal*(S(1,2)+S(2,3)+S(1,3)) +3.0_pReal*(S(4,4)+S(5,5)+S(6,6)))
   else
@@ -2155,7 +2159,7 @@ subroutine selfTest
 
   system = reshape([1.0_pReal+r(1),0.0_pReal,0.0_pReal, 0.0_pReal,1.0_pReal+r(2),0.0_pReal],[6,1])
   CoSy   = buildCoordinateSystem([1],[1],system,'cF',0.0_pReal)
-  if(any(dNeq(CoSy(1:3,1:3,1),math_I3))) error stop 'buildCoordinateSystem'
+  if (any(dNeq(CoSy(1:3,1:3,1),math_I3))) error stop 'buildCoordinateSystem'
 
   do i = 1, 10
     call random_number(C)
@@ -2201,13 +2205,13 @@ subroutine selfTest
   C(1,1) = C(1,1) + C(1,2) + 0.1_pReal
   C(4,4) = 0.5_pReal * (C(1,1) - C(1,2))
   C = lattice_symmetrize_C66(C,'cI')
-  if(dNeq(C(4,4),lattice_equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/voigt'
+  if (dNeq(C(4,4),lattice_equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/voigt'
-  if(dNeq(C(4,4),lattice_equivalent_mu(C,'reuss'),1.0e-12_pReal)) error stop 'equivalent_mu/reuss'
+  if (dNeq(C(4,4),lattice_equivalent_mu(C,'reuss'),1.0e-12_pReal)) error stop 'equivalent_mu/reuss'
 
   lambda = C(1,2)
-  if(dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'voigt')), &
+  if (dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'voigt')), &
           lattice_equivalent_nu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_nu/voigt'
-  if(dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'reuss')), &
+  if (dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'reuss')), &
           lattice_equivalent_nu(C,'reuss'),1.0e-12_pReal)) error stop 'equivalent_nu/reuss'
 
 end subroutine selfTest

src/math.f90

@@ -15,15 +15,15 @@ module math
   implicit none
   public
 #if __INTEL_COMPILER >= 1900
-  ! do not make use associated entities available to other modules
+  ! do not make use of associated entities available to other modules
   private :: &
     IO, &
     config
 #endif
 
   real(pReal),    parameter :: PI = acos(-1.0_pReal)                                                !< ratio of a circle's circumference to its diameter
-  real(pReal),    parameter :: INDEG = 180.0_pReal/PI                                               !< conversion from radian into degree
+  real(pReal),    parameter :: INDEG = 180.0_pReal/PI                                               !< conversion from radian to degree
-  real(pReal),    parameter :: INRAD = PI/180.0_pReal                                               !< conversion from degree into radian
+  real(pReal),    parameter :: INRAD = PI/180.0_pReal                                               !< conversion from degree to radian
   complex(pReal), parameter :: TWOPIIMG = cmplx(0.0_pReal,2.0_pReal*PI)                             !< Re(0.0), Im(2xPi)
 
   real(pReal), dimension(3,3), parameter :: &
@@ -132,19 +132,19 @@ pure recursive subroutine math_sort(a, istart, iend, sortDim)
   integer, intent(in),optional :: istart,iend, sortDim
   integer :: ipivot,s,e,d
 
-  if(present(istart)) then
+  if (present(istart)) then
     s = istart
   else
     s = lbound(a,2)
   endif
 
-  if(present(iend)) then
+  if (present(iend)) then
     e = iend
   else
     e = ubound(a,2)
   endif
 
-  if(present(sortDim)) then
+  if (present(sortDim)) then
     d = sortDim
   else
     d = 1
@@ -467,7 +467,7 @@ pure function math_inv33(A)
   logical     :: error
 
   call math_invert33(math_inv33,DetA,error,A)
-  if(error) math_inv33 = 0.0_pReal
+  if (error) math_inv33 = 0.0_pReal
 
 end function math_inv33
 
@@ -698,7 +698,7 @@ pure function math_sym33to6(m33,weighted)
   integer :: i
 
 
-  if(present(weighted)) then
+  if (present(weighted)) then
     w = merge(NRMMANDEL,1.0_pReal,weighted)
   else
     w = NRMMANDEL
@@ -725,7 +725,7 @@ pure function math_6toSym33(v6,weighted)
   integer :: i
 
 
-  if(present(weighted)) then
+  if (present(weighted)) then
     w = merge(INVNRMMANDEL,1.0_pReal,weighted)
   else
     w = INVNRMMANDEL
@@ -802,11 +802,11 @@ pure function math_sym3333to66(m3333,weighted)
   integer :: i,j
 
 
-  if(present(weighted)) then
+  if (present(weighted)) then
     w = merge(NRMMANDEL,1.0_pReal,weighted)
   else
     w = NRMMANDEL
-  endif
+  end if
 
 #ifndef __INTEL_COMPILER
   do concurrent(i=1:6, j=1:6)
@@ -822,7 +822,7 @@ end function math_sym3333to66
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief convert 66 matrix into symmetric 3x3x3x3 matrix
+!> @brief convert 6x6 matrix into symmetric 3x3x3x3 matrix
 !> @details Weighted conversion (default) rearranges according to Nye and weights shear
 ! components according to Mandel. Advisable for matrix operations.
 ! Unweighted conversion only rearranges order according to Nye
@@ -837,7 +837,7 @@ pure function math_66toSym3333(m66,weighted)
   integer :: i,j
 
 
-  if(present(weighted)) then
+  if (present(weighted)) then
     w = merge(INVNRMMANDEL,1.0_pReal,weighted)
   else
     w = INVNRMMANDEL
@@ -1047,7 +1047,7 @@ subroutine math_eigh33(w,v,m)
                v(2,2) + m(2, 3) * w(1), &
               (m(1,1) - w(1)) * (m(2,2) - w(1)) - v(3,2)]
   norm = norm2(v(1:3, 1))
-  fallback1: if(norm < threshold) then
+  fallback1: if (norm < threshold) then
     call math_eigh(w,v,error,m)
   else fallback1
     v(1:3,1) = v(1:3, 1) / norm
@@ -1055,7 +1055,7 @@ subroutine math_eigh33(w,v,m)
                  v(2,2) + m(2, 3) * w(2), &
                 (m(1,1) - w(2)) * (m(2,2) - w(2)) - v(3,2)]
     norm = norm2(v(1:3, 2))
-    fallback2: if(norm < threshold) then
+    fallback2: if (norm < threshold) then
       call math_eigh(w,v,error,m)
     else fallback2
       v(1:3,2) = v(1:3, 2) / norm
@@ -1092,7 +1092,7 @@ pure function math_rotationalPart(F) result(R)
   I_F = [math_trace33(F), 0.5*(math_trace33(F)**2 - math_trace33(matmul(F,F)))]
 
   x = math_clip(I_C(1)**2 -3.0_pReal*I_C(2),0.0_pReal)**(3.0_pReal/2.0_pReal)
-  if(dNeq0(x)) then
+  if (dNeq0(x)) then
     Phi = acos(math_clip((I_C(1)**3 -4.5_pReal*I_C(1)*I_C(2) +13.5_pReal*I_C(3))/x,-1.0_pReal,1.0_pReal))
     lambda = I_C(1) +(2.0_pReal * sqrt(math_clip(I_C(1)**2-3.0_pReal*I_C(2),0.0_pReal))) &
                     *cos((Phi-2.0_pReal * PI*[1.0_pReal,2.0_pReal,3.0_pReal])/3.0_pReal)
@@ -1154,7 +1154,7 @@ function math_eigvalsh33(m)
     - 2.0_pReal/27.0_pReal*I(1)**3.0_pReal &
     - I(3)                                                                                          ! different from http://arxiv.org/abs/physics/0610206 (this formulation was in DAMASK)
 
-  if(all(abs([P,Q]) < TOL)) then
+  if (all(abs([P,Q]) < TOL)) then
     math_eigvalsh33 = math_eigvalsh(m)
   else
     rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal
@@ -1277,7 +1277,7 @@ real(pReal) pure elemental function math_clip(a, left, right)
   if (present(left))  math_clip = max(left,math_clip)
   if (present(right)) math_clip = min(right,math_clip)
   if (present(left) .and. present(right)) then
-    if(left>right) error stop 'left > right'
+    if (left>right) error stop 'left > right'
   endif
 
 end function math_clip
@@ -1323,38 +1323,38 @@ subroutine selfTest
     error stop 'math_expand [1,2] by [1,2,3] => [1,2,2,1,1,1]'
 
   call math_sort(sort_in_,1,3,2)
-  if(any(sort_in_ /= sort_out_)) &
+  if (any(sort_in_ /= sort_out_)) &
     error stop 'math_sort'
 
-  if(any(math_range(5) /= range_out_)) &
+  if (any(math_range(5) /= range_out_)) &
     error stop 'math_range'
 
-  if(any(dNeq(math_exp33(math_I3,0),math_I3))) &
+  if (any(dNeq(math_exp33(math_I3,0),math_I3))) &
     error stop 'math_exp33(math_I3,1)'
-  if(any(dNeq(math_exp33(math_I3,128),exp(1.0_pReal)*math_I3))) &
+  if (any(dNeq(math_exp33(math_I3,128),exp(1.0_pReal)*math_I3))) &
     error stop 'math_exp33(math_I3,128)'
 
   call random_number(v9)
-  if(any(dNeq(math_33to9(math_9to33(v9)),v9))) &
+  if (any(dNeq(math_33to9(math_9to33(v9)),v9))) &
     error stop 'math_33to9/math_9to33'
 
   call random_number(t99)
-  if(any(dNeq(math_3333to99(math_99to3333(t99)),t99))) &
+  if (any(dNeq(math_3333to99(math_99to3333(t99)),t99))) &
     error stop 'math_3333to99/math_99to3333'
 
   call random_number(v6)
-  if(any(dNeq(math_sym33to6(math_6toSym33(v6)),v6))) &
+  if (any(dNeq(math_sym33to6(math_6toSym33(v6)),v6))) &
     error stop 'math_sym33to6/math_6toSym33'
 
   call random_number(t66)
-  if(any(dNeq(math_sym3333to66(math_66toSym3333(t66)),t66,1.0e-15_pReal))) &
+  if (any(dNeq(math_sym3333to66(math_66toSym3333(t66)),t66,1.0e-15_pReal))) &
     error stop 'math_sym3333to66/math_66toSym3333'
 
-  if(any(dNeq(math_3333toVoigt66(math_Voigt66to3333(t66)),t66,1.0e-15_pReal))) &
+  if (any(dNeq(math_3333toVoigt66(math_Voigt66to3333(t66)),t66,1.0e-15_pReal))) &
     error stop 'math_3333toVoigt66/math_Voigt66to3333'
 
   call random_number(v6)
-  if(any(dNeq0(math_6toSym33(v6) - math_symmetric33(math_6toSym33(v6))))) &
+  if (any(dNeq0(math_6toSym33(v6) - math_symmetric33(math_6toSym33(v6))))) &
     error stop 'math_symmetric33'
 
   call random_number(v3_1)
@@ -1362,34 +1362,34 @@ subroutine selfTest
   call random_number(v3_3)
   call random_number(v3_4)
 
-  if(dNeq(abs(dot_product(math_cross(v3_1-v3_4,v3_2-v3_4),v3_3-v3_4))/6.0, &
+  if (dNeq(abs(dot_product(math_cross(v3_1-v3_4,v3_2-v3_4),v3_3-v3_4))/6.0, &
           math_volTetrahedron(v3_1,v3_2,v3_3,v3_4),tol=1.0e-12_pReal)) &
     error stop 'math_volTetrahedron'
 
   call random_number(t33)
-  if(dNeq(math_det33(math_symmetric33(t33)),math_detSym33(math_symmetric33(t33)),tol=1.0e-12_pReal)) &
+  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))))) &
+  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)))) &
+  if (any(dNeq0(math_eye(3),math_inv33(math_I3)))) &
     error stop 'math_inv33(math_I3)'
 
   do while(abs(math_det33(t33))<1.0e-9_pReal)
     call random_number(t33)
   enddo
-  if(any(dNeq0(matmul(t33,math_inv33(t33)) - math_eye(3),tol=1.0e-9_pReal))) &
+  if (any(dNeq0(matmul(t33,math_inv33(t33)) - math_eye(3),tol=1.0e-9_pReal))) &
     error stop 'math_inv33'
 
   call math_invert33(t33_2,det,e,t33)
-  if(any(dNeq0(matmul(t33,t33_2) - math_eye(3),tol=1.0e-9_pReal)) .or. e) &
+  if (any(dNeq0(matmul(t33,t33_2) - math_eye(3),tol=1.0e-9_pReal)) .or. e) &
     error stop 'math_invert33: T:T^-1 != I'
-  if(dNeq(det,math_det33(t33),tol=1.0e-12_pReal)) &
+  if (dNeq(det,math_det33(t33),tol=1.0e-12_pReal)) &
     error stop 'math_invert33 (determinant)'
 
   call math_invert(t33_2,e,t33)
-  if(any(dNeq0(matmul(t33,t33_2) - math_eye(3),tol=1.0e-9_pReal)) .or. e) &
+  if (any(dNeq0(matmul(t33,t33_2) - math_eye(3),tol=1.0e-9_pReal)) .or. e) &
     error stop 'math_invert t33'
 
   do while(math_det33(t33)<1.0e-2_pReal)                                                            ! O(det(F)) = 1
@@ -1397,7 +1397,7 @@ subroutine selfTest
   enddo
   t33_2 = math_rotationalPart(transpose(t33))
   t33   = math_rotationalPart(t33)
-  if(any(dNeq0(matmul(t33_2,t33) - math_I3,tol=1.0e-10_pReal))) &
+  if (any(dNeq0(matmul(t33_2,t33) - math_I3,tol=1.0e-10_pReal))) &
     error stop 'math_rotationalPart'
 
   call random_number(r)
@@ -1405,33 +1405,33 @@ subroutine selfTest
   txx = math_eye(d)
   allocate(txx_2(d,d))
   call math_invert(txx_2,e,txx)
-  if(any(dNeq0(txx_2,txx) .or. e)) &
+  if (any(dNeq0(txx_2,txx) .or. e)) &
     error stop 'math_invert(txx)/math_eye'
 
   call math_invert(t99_2,e,t99) ! not sure how likely it is that we get a singular matrix
-  if(any(dNeq0(matmul(t99_2,t99)-math_eye(9),tol=1.0e-9_pReal)) .or. e) &
+  if (any(dNeq0(matmul(t99_2,t99)-math_eye(9),tol=1.0e-9_pReal)) .or. e) &
     error stop 'math_invert(t99)'
 
-  if(any(dNeq(math_clip([4.0_pReal,9.0_pReal],5.0_pReal,6.5_pReal),[5.0_pReal,6.5_pReal]))) &
+  if (any(dNeq(math_clip([4.0_pReal,9.0_pReal],5.0_pReal,6.5_pReal),[5.0_pReal,6.5_pReal]))) &
     error stop 'math_clip'
 
-  if(math_factorial(10) /= 3628800) &
+  if (math_factorial(10) /= 3628800) &
     error stop 'math_factorial'
 
-  if(math_binomial(49,6) /= 13983816) &
+  if (math_binomial(49,6) /= 13983816) &
     error stop 'math_binomial'
 
-  if(math_multinomial([1,2,3,4]) /= 12600) &
+  if (math_multinomial([1,2,3,4]) /= 12600) &
     error stop 'math_multinomial'
 
   ijk = cshift([1,2,3],int(r*1.0e2_pReal))
-  if(dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),+1.0_pReal)) &
+  if (dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),+1.0_pReal)) &
     error stop 'math_LeviCivita(even)'
   ijk = cshift([3,2,1],int(r*2.0e2_pReal))
-  if(dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),-1.0_pReal)) &
+  if (dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),-1.0_pReal)) &
     error stop 'math_LeviCivita(odd)'
   ijk = cshift([2,2,1],int(r*2.0e2_pReal))
-  if(dNeq0(math_LeviCivita(ijk(1),ijk(2),ijk(3)))) &
+  if (dNeq0(math_LeviCivita(ijk(1),ijk(2),ijk(3)))) &
     error stop 'math_LeviCivita'
 
 end subroutine selfTest

src/phase.f90

@@ -5,6 +5,7 @@
 !--------------------------------------------------------------------------------------------------
 module phase
   use prec
+  use constants
   use math
   use rotations
   use IO

src/phase_mechanical_eigen_thermalexpansion.f90

@@ -58,14 +58,14 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics)
         associate(prm  => param(kinematics_thermal_expansion_instance(p)))
           kinematic_type => kinematics%get(k)
 
-          prm%T_ref = kinematic_type%get_asFloat('T_ref', defaultVal=0.0_pReal)
+          prm%T_ref = kinematic_type%get_asFloat('T_ref', defaultVal=T_ROOM)
 
           prm%A(1,1,1) = kinematic_type%get_asFloat('A_11')
-          prm%A(1,1,2) = kinematic_type%get_asFloat('A_11,T',defaultVal=0.0_pReal)
+          prm%A(1,1,2) = kinematic_type%get_asFloat('A_11,T',  defaultVal=0.0_pReal)
           prm%A(1,1,3) = kinematic_type%get_asFloat('A_11,T^2',defaultVal=0.0_pReal)
           if (any(phase_lattice(p) == ['hP','tI'])) then
             prm%A(3,3,1) = kinematic_type%get_asFloat('A_33')
-            prm%A(3,3,2) = kinematic_type%get_asFloat('A_33,T',defaultVal=0.0_pReal)
+            prm%A(3,3,2) = kinematic_type%get_asFloat('A_33,T',  defaultVal=0.0_pReal)
             prm%A(3,3,3) = kinematic_type%get_asFloat('A_33,T^2',defaultVal=0.0_pReal)
           end if
           do i=1, size(prm%A,3)
@@ -98,14 +98,14 @@ module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me)
 
   associate(prm => param(kinematics_thermal_expansion_instance(ph)))
     Li = dot_T * ( &
-                  prm%A(1:3,1:3,1)*(T - prm%T_ref)**0 &                                             ! constant  coefficient
+                  prm%A(1:3,1:3,1) &                                                                ! constant  coefficient
-                + prm%A(1:3,1:3,2)*(T - prm%T_ref)**1 &                                             ! linear    coefficient
+                + prm%A(1:3,1:3,2)*(T - prm%T_ref)**1         &                                     ! linear    coefficient
-                + prm%A(1:3,1:3,3)*(T - prm%T_ref)**2 &                                             ! quadratic coefficient
+                + prm%A(1:3,1:3,3)*(T - prm%T_ref)**2         &                                     ! quadratic coefficient
                 ) / &
          (1.0_pReal &
-               + prm%A(1:3,1:3,1)*(T - prm%T_ref)**1 / 1. &
+               + prm%A(1:3,1:3,1)*(T - prm%T_ref)**1 / 1.0_pReal &
-               + prm%A(1:3,1:3,2)*(T - prm%T_ref)**2 / 2. &
+               + prm%A(1:3,1:3,2)*(T - prm%T_ref)**2 / 2.0_pReal &
-               + prm%A(1:3,1:3,3)*(T - prm%T_ref)**3 / 3. &
+               + prm%A(1:3,1:3,3)*(T - prm%T_ref)**3 / 3.0_pReal &
          )
   end associate
   dLi_dTstar = 0.0_pReal

src/phase_mechanical_elastic.f90

@@ -1,13 +1,15 @@
 submodule(phase:mechanical) elastic
 
   type :: tParameters
-    real(pReal) :: &
+    real(pReal),dimension(3) :: &
       C_11 = 0.0_pReal, &
       C_12 = 0.0_pReal, &
       C_13 = 0.0_pReal, &
       C_33 = 0.0_pReal, &
       C_44 = 0.0_pReal, &
       C_66 = 0.0_pReal
+    real(pReal) :: &
+      T_ref
   end type tParameters
 
   type(tParameters), allocatable, dimension(:) :: param
@@ -15,7 +17,7 @@ submodule(phase:mechanical) elastic
 contains
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Initialize elasticity.
+!> @brief initialize elasticity
 !--------------------------------------------------------------------------------------------------
 module subroutine elastic_init(phases)
 
@@ -28,7 +30,7 @@ module subroutine elastic_init(phases)
     phase, &
     mech, &
     elastic
-
+  logical :: thermal_active
 
   print'(/,1x,a)', '<<<+-  phase:mechanical:elastic init  -+>>>'
   print'(/,1x,a)', '<<<+-  phase:mechanical:elastic:Hooke init  -+>>>'
@@ -45,15 +47,35 @@ module subroutine elastic_init(phases)
 
     associate(prm => param(ph))
 
-      prm%C_11 = elastic%get_asFloat('C_11')
+      prm%T_ref = elastic%get_asFloat('T_ref', defaultVal=T_ROOM)
-      prm%C_12 = elastic%get_asFloat('C_12')
-      prm%C_44 = elastic%get_asFloat('C_44')
 
+      prm%C_11(1) = elastic%get_asFloat('C_11')
+      prm%C_11(2) = elastic%get_asFloat('C_11,T',  defaultVal=0.0_pReal)
+      prm%C_11(3) = elastic%get_asFloat('C_11,T^2',defaultVal=0.0_pReal)
+ 
+      prm%C_12(1) = elastic%get_asFloat('C_12')
+      prm%C_12(2) = elastic%get_asFloat('C_12,T',  defaultVal=0.0_pReal)
+      prm%C_12(3) = elastic%get_asFloat('C_12,T^2',defaultVal=0.0_pReal)
+ 
+      prm%C_44(1) = elastic%get_asFloat('C_44')
+      prm%C_44(2) = elastic%get_asFloat('C_44,T',  defaultVal=0.0_pReal)
+      prm%C_44(3) = elastic%get_asFloat('C_44,T^2',defaultVal=0.0_pReal)
+ 
       if (any(phase_lattice(ph) == ['hP','tI'])) then
-        prm%C_13 = elastic%get_asFloat('C_13')
+        prm%C_13(1) = elastic%get_asFloat('C_13')
-        prm%C_33 = elastic%get_asFloat('C_33')
+        prm%C_13(2) = elastic%get_asFloat('C_13,T',  defaultVal=0.0_pReal)
+        prm%C_13(3) = elastic%get_asFloat('C_13,T^2',defaultVal=0.0_pReal)
+ 
+        prm%C_33(1) = elastic%get_asFloat('C_33')
+        prm%C_33(2) = elastic%get_asFloat('C_33,T',  defaultVal=0.0_pReal)
+        prm%C_33(3) = elastic%get_asFloat('C_33,T^2',defaultVal=0.0_pReal)
+      end if
+
+      if (phase_lattice(ph) == 'tI') then
+        prm%C_66(1) = elastic%get_asFloat('C_66')
+        prm%C_66(2) = elastic%get_asFloat('C_66,T',  defaultVal=0.0_pReal)
+        prm%C_66(3) = elastic%get_asFloat('C_66,T^2',defaultVal=0.0_pReal)
       end if
-      if (phase_lattice(ph) == 'tI') prm%C_66 = elastic%get_asFloat('C_66')
 
     end associate
   end do
@@ -62,28 +84,51 @@ end subroutine elastic_init
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Return 6x6 elasticity tensor.
+!> @brief return 6x6 elasticity tensor
 !--------------------------------------------------------------------------------------------------
 module function elastic_C66(ph,en) result(C66)
 
   integer, intent(in) :: &
     ph, &
     en
+
   real(pReal), dimension(6,6) :: C66
+  real(pReal) :: T
 
 
   associate(prm => param(ph))
     C66 = 0.0_pReal
-    C66(1,1) = prm%C_11
+    T = thermal_T(ph,en)
-    C66(1,2) = prm%C_12
+
-    C66(4,4) = prm%C_44
+    C66(1,1) = prm%C_11(1) &
+             + prm%C_11(2)*(T - prm%T_ref)**1 &
+             + prm%C_11(3)*(T - prm%T_ref)**2
+
+    C66(1,2) = prm%C_12(1) &
+             + prm%C_12(2)*(T - prm%T_ref)**1 &
+             + prm%C_12(3)*(T - prm%T_ref)**2
+
+    C66(4,4) = prm%C_44(1) &
+             + prm%C_44(2)*(T - prm%T_ref)**1 &
+             + prm%C_44(3)*(T - prm%T_ref)**2
+
 
     if (any(phase_lattice(ph) == ['hP','tI'])) then
-      C66(1,3) = prm%C_13
+      C66(1,3) = prm%C_13(1) &
-      C66(3,3) = prm%C_33
+               + prm%C_13(2)*(T - prm%T_ref)**1 &
+               + prm%C_13(3)*(T - prm%T_ref)**2
+
+      C66(3,3) = prm%C_33(1) &
+               + prm%C_33(2)*(T - prm%T_ref)**1 &
+               + prm%C_33(3)*(T - prm%T_ref)**2
+
     end if
 
-    if (phase_lattice(ph) == 'tI') C66(6,6) = prm%C_66
+    if (phase_lattice(ph) == 'tI') then
+      C66(6,6) = prm%C_66(1) &
+               + prm%C_66(2)*(T - prm%T_ref)**1 &
+               + prm%C_66(3)*(T - prm%T_ref)**2
+    end if
 
     C66 = lattice_symmetrize_C66(C66,phase_lattice(ph))
 
@@ -93,7 +138,7 @@ end function elastic_C66
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Return shear modulus.
+!> @brief return shear modulus
 !--------------------------------------------------------------------------------------------------
 module function elastic_mu(ph,en) result(mu)
 
@@ -110,7 +155,7 @@ end function elastic_mu
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Return Poisson ratio.
+!> @brief return Poisson ratio
 !--------------------------------------------------------------------------------------------------
 module function elastic_nu(ph,en) result(nu)
 
@@ -128,7 +173,7 @@ end function elastic_nu
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
+!> @brief return the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
 !> the elastic and intermediate deformation gradients using Hooke's law
 ! ToDo: Use Voigt matrix directly
 !--------------------------------------------------------------------------------------------------

src/phase_mechanical_plastic_dislotungsten.f90

@@ -7,9 +7,6 @@
 !--------------------------------------------------------------------------------------------------
 submodule(phase:plastic) dislotungsten
 
-  real(pReal), parameter :: &
-    kB = 1.38e-23_pReal                                                                             !< Boltzmann constant in J/Kelvin
-
   type :: tParameters
     real(pReal) :: &
       D    = 1.0_pReal, &                                                                           !< grain size
@@ -344,7 +341,7 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
       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*mu*prm%D_0*exp(-prm%Q_cl/(kB*T))*prm%f_at/(2.0_pReal*PI*kB*T)) &
+      v_cl = (3.0_pReal*mu*prm%D_0*exp(-prm%Q_cl/(K_B*T))*prm%f_at/(2.0_pReal*PI*K_B*T)) &
            * (1.0_pReal/(d_hat+prm%d_caron))
       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
@@ -475,7 +472,7 @@ pure subroutine kinetics(Mp,T,ph,en, &
     if (present(tau_pos_out)) tau_pos_out = tau_pos
     if (present(tau_neg_out)) tau_neg_out = tau_neg
 
-    associate(BoltzmannRatio  => prm%Q_s/(kB*T), &
+    associate(BoltzmannRatio  => prm%Q_s/(K_B*T), &
               b_rho_half      => stt%rho_mob(:,en) * prm%b_sl * 0.5_pReal, &
               effectiveLength => dst%Lambda_sl(:,en) - prm%w)
 

src/phase_mechanical_plastic_dislotwin.f90

@@ -9,9 +9,6 @@
 !--------------------------------------------------------------------------------------------------
 submodule(phase:plastic) dislotwin
 
-  real(pReal), parameter :: &
-    kB = 1.38e-23_pReal                                                                             !< Boltzmann constant in J/Kelvin
-
   type :: tParameters
     real(pReal) :: &
       Q_cl                = 1.0_pReal, &                                                            !< activation energy for dislocation climb
@@ -31,7 +28,7 @@ submodule(phase:plastic) dislotwin
       delta_G             = 1.0_pReal, &                                                            !< Free energy difference between austensite and martensite
       i_tr                = 1.0_pReal, &                                                            !< adjustment parameter to calculate MFP for transformation
       h                   = 1.0_pReal, &                                                            !< Stack height of hex nucleus
-      T_ref               = 0.0_pReal, &
+      T_ref               = T_ROOM, &
       a_cI                = 1.0_pReal, &
       a_cF                = 1.0_pReal
     real(pReal),                            dimension(2) :: &
@@ -150,7 +147,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
 
 
   myPlasticity = plastic_active('dislotwin')
-  if(count(myPlasticity) == 0) return
+  if (count(myPlasticity) == 0) return
 
   print'(/,1x,a)', '<<<+-  phase:mechanical:plastic:dislotwin init  -+>>>'
   print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
@@ -173,7 +170,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
 
 
   do ph = 1, phases%length
-    if(.not. myPlasticity(ph)) cycle
+    if (.not. myPlasticity(ph)) cycle
 
     associate(prm => param(ph), dot => dotState(ph), stt => state(ph), dst => dependentState(ph))
 
@@ -368,7 +365,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
 
 !--------------------------------------------------------------------------------------------------
 ! parameters required for several mechanisms and their interactions
-    if(prm%sum_N_sl + prm%sum_N_tw + prm%sum_N_tw > 0) &
+    if (prm%sum_N_sl + prm%sum_N_tw + prm%sum_N_tw > 0) &
       prm%D = pl%get_asFloat('D')
 
     if (prm%sum_N_tw + prm%sum_N_tr > 0) &
@@ -425,7 +422,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
     stt%gamma_sl=>plasticState(ph)%state(startIndex:endIndex,:)
     dot%gamma_sl=>plasticState(ph)%dotState(startIndex:endIndex,:)
     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'
+    if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
 
     startIndex = endIndex + 1
     endIndex   = endIndex + prm%sum_N_tw
@@ -473,9 +470,7 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC)
   integer,     intent(in) :: &
     ph, en
   real(pReal), dimension(6,6) :: &
-    homogenizedC
+    homogenizedC, &
-
-  real(pReal), dimension(6,6) :: &
     C
   real(pReal), dimension(:,:,:), allocatable :: &
     C66_tw, &
@@ -500,8 +495,8 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC)
         homogenizedC = homogenizedC &
                      + stt%f_tw(i,en)*C66_tw(1:6,1:6,i)
       end do
-     end if twinActive 
+     end if twinActive
-    
+
     transActive: if (prm%sum_N_tr > 0) then
       C66_tr    = lattice_C66_trans(prm%N_tr,C,prm%lattice_tr,0.0_pReal,prm%a_cI,prm%a_cF)
       do i=1,prm%sum_N_tr
@@ -597,9 +592,9 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
   Lp      = Lp      * f_unrotated
   dLp_dMp = dLp_dMp * f_unrotated
 
-  shearBandingContribution: if(dNeq0(prm%v_sb)) then
+  shearBandingContribution: if (dNeq0(prm%v_sb)) then
 
-    E_kB_T = prm%E_sb/(kB*T)
+    E_kB_T = prm%E_sb/(K_B*T)
     call math_eigh33(eigValues,eigVectors,Mp)                                                       ! is Mp symmetric by design?
 
     do i = 1,6
@@ -696,8 +691,8 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
                         * (prm%Gamma_sf(1) + prm%Gamma_sf(2) * T) / (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)) &
+          v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(K_B*T)) &
-               * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal)
+               * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(K_B*T)) - 1.0_pReal)
 
           dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,en) &
                                / (d_hat-prm%d_caron(i))
@@ -848,7 +843,7 @@ module subroutine plastic_dislotwin_results(ph,group)
                                        'threshold stress for twinning','Pa',prm%systems_tw)
 
         case('f_tr')
-          if(prm%sum_N_tr>0) call results_writeDataset(stt%f_tr,group,trim(prm%output(ou)), &
+          if (prm%sum_N_tr>0) call results_writeDataset(stt%f_tr,group,trim(prm%output(ou)), &
                                                        'martensite volume fraction','m³/m³')
 
       end select
@@ -909,7 +904,7 @@ pure subroutine kinetics_sl(Mp,T,ph,en, &
     significantStress: where(tau_eff > tol_math_check)
       stressRatio    = tau_eff/prm%tau_0
       StressRatio_p  = stressRatio** prm%p
-      Q_kB_T = prm%Q_sl/(kB*T)
+      Q_kB_T = prm%Q_sl/(K_B*T)
       v_wait_inverse = exp(Q_kB_T*(1.0_pReal-StressRatio_p)** prm%q) &
                      / prm%v_0
       v_run_inverse  = prm%B/(tau_eff*prm%b_sl)
@@ -931,8 +926,8 @@ pure subroutine kinetics_sl(Mp,T,ph,en, &
 
   end associate
 
-  if(present(ddot_gamma_dtau_sl)) ddot_gamma_dtau_sl = ddot_gamma_dtau
+  if (present(ddot_gamma_dtau_sl)) ddot_gamma_dtau_sl = ddot_gamma_dtau
-  if(present(tau_sl))             tau_sl             = tau
+  if (present(tau_sl))             tau_sl             = tau
 
 end subroutine kinetics_sl
 
@@ -982,7 +977,7 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,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)))/&
                   (prm%L_tw*prm%b_sl(i))*&
-                  (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tw(i,en)-tau(i))))
+                  (1.0_pReal-exp(-prm%V_cs/(K_B*T)*(dst%tau_r_tw(i,en)-tau(i))))
         else
           Ndot0=0.0_pReal
         end if
@@ -1002,7 +997,7 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
 
   end associate
 
-  if(present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw = ddot_gamma_dtau
+  if (present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw = ddot_gamma_dtau
 
 end subroutine kinetics_tw
 
@@ -1051,7 +1046,7 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,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)))/&
                   (prm%L_tr*prm%b_sl(i))*&
-                  (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tr(i,en)-tau(i))))
+                  (1.0_pReal-exp(-prm%V_cs/(K_B*T)*(dst%tau_r_tr(i,en)-tau(i))))
         else
           Ndot0=0.0_pReal
         end if
@@ -1071,7 +1066,7 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
 
   end associate
 
-  if(present(ddot_gamma_dtau_tr)) ddot_gamma_dtau_tr = ddot_gamma_dtau
+  if (present(ddot_gamma_dtau_tr)) ddot_gamma_dtau_tr = ddot_gamma_dtau
 
 end subroutine kinetics_tr
 

src/phase_mechanical_plastic_nonlocal.f90

@@ -19,9 +19,6 @@ submodule(phase:plastic) nonlocal
 
   type(tGeometry), dimension(:), allocatable :: geom
 
-  real(pReal), parameter :: &
-    kB = 1.38e-23_pReal                                                                             !< Boltzmann constant in J/Kelvin
-
   ! storage order of dislocation types
   integer, dimension(*), parameter :: &
     sgl = [1,2,3,4,5,6,7,8]                                                                         !< signed (single)
@@ -1094,9 +1091,9 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
 
   ! thermally activated annihilation of edge dipoles by climb
   rhoDotThermalAnnihilation = 0.0_pReal
-  D_SD = prm%D_0 * exp(-prm%Q_cl / (kB * Temperature))                                              ! eq. 3.53
+  D_SD = prm%D_0 * exp(-prm%Q_cl / (K_B * Temperature))                                              ! eq. 3.53
   v_climb = D_SD * mu * prm%V_at &
-          / (PI * (1.0_pReal-nu) * (dUpper(:,1) + dLower(:,1)) * kB * Temperature)              ! eq. 3.54
+          / (PI * (1.0_pReal-nu) * (dUpper(:,1) + dLower(:,1)) * K_B * Temperature)              ! eq. 3.54
   forall (s = 1:prm%sum_N_sl, dUpper(s,1) > dLower(s,1)) &
     rhoDotThermalAnnihilation(s,9) = max(- 4.0_pReal * rhoDip(s,1) * v_climb(s) / (dUpper(s,1) - dLower(s,1)), &
                                          - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) &
@@ -1671,9 +1668,9 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, T,
       activationEnergy_P = criticalStress_P * activationVolume_P
       tauRel_P = min(1.0_pReal, tauEff / criticalStress_P)
       tPeierls = 1.0_pReal / prm%nu_a &
-               * exp(activationEnergy_P / (kB * T) &
+               * exp(activationEnergy_P / (K_B * T) &
                      * (1.0_pReal - tauRel_P**prm%p)**prm%q)
-      dtPeierls_dtau = merge(tPeierls * prm%p * prm%q * activationVolume_P / (kB * T) &
+      dtPeierls_dtau = merge(tPeierls * prm%p * prm%q * activationVolume_P / (K_B * T) &
                              * (1.0_pReal - tauRel_P**prm%p)**(prm%q-1.0_pReal) * tauRel_P**(prm%p-1.0_pReal), &
                              0.0_pReal, &
                              tauEff < criticalStress_P)
@@ -1685,8 +1682,8 @@ pure subroutine kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, tauThreshold, c, T,
       criticalStress_S = prm%Q_sol / activationVolume_S
       tauRel_S = min(1.0_pReal, tauEff / criticalStress_S)
       tSolidSolution = 1.0_pReal /  prm%nu_a &
-                     * exp(prm%Q_sol / (kB * T)* (1.0_pReal - tauRel_S**prm%p)**prm%q)
+                     * exp(prm%Q_sol / (K_B * T)* (1.0_pReal - tauRel_S**prm%p)**prm%q)
-      dtSolidSolution_dtau = merge(tSolidSolution * prm%p * prm%q * activationVolume_S / (kB * T) &
+      dtSolidSolution_dtau = merge(tSolidSolution * prm%p * prm%q * activationVolume_S / (K_B * T) &
                                    * (1.0_pReal - tauRel_S**prm%p)**(prm%q-1.0_pReal)* tauRel_S**(prm%p-1.0_pReal), &
                                    0.0_pReal, &
                                    tauEff < criticalStress_S)