Merge remote-tracking branch 'origin/development' into 134-output_none

.gitlab-ci.yml

@@ -3,20 +3,18 @@ stages:
   - prepare
   - python
   - compile
-  - setup
   - fortran
   - performance
   - deploy
-  - update_master
+  - finalize
 
 
 ###################################################################################################
 default:
   before_script:
     - ${LOCAL_HOME}/bin/queue ${CI_JOB_ID}
-    - source $DAMASKROOT/env/DAMASK.sh
+    - source env/DAMASK.sh
-    - export PATH=${TESTROOT}/bin:$PATH
+    - export PATH=${TESTROOT}/bin:${PATH}
-    - cd $DAMASKROOT/PRIVATE/testing
     - echo Job start:" $(date)"
   after_script:
     - echo Job end:" $(date)"
@@ -33,20 +31,19 @@ variables:
   # Shortcut names
   # ===============================================================================================
   TESTROOT:   "$LOCAL_HOME/GitLabCI_Pipeline_$CI_PIPELINE_ID"
-  DAMASKROOT: "$LOCAL_HOME/GitLabCI_Pipeline_$CI_PIPELINE_ID/DAMASK"
 
   # ===============================================================================================
   # Names of module files to load
   # ===============================================================================================
   # ++++++++++++ Compiler +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-  IntelCompiler:  "Compiler/Intel/19.1.2 Libraries/IMKL/2020"
+  COMPILER_INTEL: "Compiler/Intel/19.1.2 Libraries/IMKL/2020"
-  GNUCompiler:    "Compiler/GNU/10"
+  COMPILER_GNU:   "Compiler/GNU/10"
   # ++++++++++++ MPI ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-  MPI_Intel:      "MPI/Intel/19.1.2/IntelMPI/2019"
+  MPI_INTEL:      "MPI/Intel/19.1.2/IntelMPI/2019"
   MPI_GNU:        "MPI/GNU/10/OpenMPI/4.1.1"
   # ++++++++++++ PETSc ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-  PETSc_Intel:    "Libraries/PETSc/3.16.1/Intel-19.1.2-IntelMPI-2019"
+  PETSC_INTEL:    "Libraries/PETSc/3.16.1/Intel-19.1.2-IntelMPI-2019"
-  PETSc_GNU:      "Libraries/PETSc/3.16.1/GNU-10-OpenMPI-4.1.1"
+  PETSC_GNU:      "Libraries/PETSc/3.16.1/GNU-10-OpenMPI-4.1.1"
   # ++++++++++++ MSC Marc +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   MSC:            "FEM/MSC/2021.2"
   IntelMarc:      "Compiler/Intel/19.1.2 Libraries/IMKL/2020"
@@ -54,97 +51,107 @@ variables:
 
 
 ###################################################################################################
-checkout:
+create_testroot:
   stage: prepare
   before_script:
     - ${LOCAL_HOME}/bin/queue ${CI_JOB_ID}
     - echo Job start:" $(date)"
   script:
-    - mkdir -p ${DAMASKROOT}
+    - mkdir -p ${TESTROOT}
-    - cd ${DAMASKROOT}
-    - git clone -q git@git.damask.mpie.de:damask/DAMASK.git .
-    - git checkout ${CI_COMMIT_SHA}
-    - git submodule update --init
 
 
 ###################################################################################################
 pytest:
   stage: python
   script:
-    - PYTHONPATH=${CI_PROJECT_DIR}/python
+    - cd python
-    - cd ${CI_PROJECT_DIR}/python
     - pytest --basetemp ${TESTROOT}/python -v --cov --cov-report=term
     - coverage report --fail-under=90
 
 mypy:
   stage: python
   script:
-    - cd ${CI_PROJECT_DIR}/python
+    - cd python
     - mypy damask
 
 
 ###################################################################################################
-compile_grid_Intel:
+test_grid_Intel:
   stage: compile
   script:
-    - module load $IntelCompiler $MPI_Intel $PETSc_Intel
+    - module load ${COMPILER_INTEL} ${MPI_INTEL} ${PETSC_INTEL}
-    - cd pytest
+    - cd PRIVATE/testing/pytest
     - pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_Intel
 
-compile_mesh_Intel:
+test_mesh_Intel:
   stage: compile
   script:
-    - module load $IntelCompiler $MPI_Intel $PETSc_Intel
+    - module load ${COMPILER_INTEL} ${MPI_INTEL} ${PETSC_INTEL}
-    - cd pytest
+    - cd PRIVATE/testing/pytest
     - pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_Intel
 
-compile_grid_GNU:
+test_grid_GNU:
   stage: compile
   script:
-    - module load $GNUCompiler $MPI_GNU $PETSc_GNU
+    - module load ${COMPILER_GNU} ${MPI_GNU} ${PETSC_GNU}
-    - cd pytest
+    - cd PRIVATE/testing/pytest
     - pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_GNU
 
-compile_mesh_GNU:
+test_mesh_GNU:
   stage: compile
   script:
-    - module load $GNUCompiler $MPI_GNU $PETSc_GNU
+    - module load ${COMPILER_GNU} ${MPI_GNU} ${PETSC_GNU}
-    - cd pytest
+    - cd PRIVATE/testing/pytest
     - pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_GNU
 
-compile_Marc:
+test_Marc:
   stage: compile
   script:
     - module load $IntelMarc $HDF5Marc $MSC
-    - cd pytest
+    - cd PRIVATE/testing/pytest
     - pytest -k 'compile and Marc' --basetemp ${TESTROOT}/compile_Marc
 
 
-###################################################################################################
 setup_grid:
-  stage: setup
+  stage: compile
   script:
-    - module load $IntelCompiler $MPI_Intel $PETSc_Intel
+    - module load ${COMPILER_INTEL} ${MPI_INTEL} ${PETSC_INTEL}
     - cd $(mktemp -d)
     - cmake -DDAMASK_SOLVER=GRID -DCMAKE_INSTALL_PREFIX=${TESTROOT} ${CI_PROJECT_DIR}
     - make -j2 all install
 
 setup_mesh:
-  stage: setup
+  stage: compile
   script:
-    - module load $IntelCompiler $MPI_Intel $PETSc_Intel
+    - module load ${COMPILER_INTEL} ${MPI_INTEL} ${PETSC_INTEL}
     - cd $(mktemp -d)
     - cmake -DDAMASK_SOLVER=MESH -DCMAKE_INSTALL_PREFIX=${TESTROOT} ${CI_PROJECT_DIR}
     - make -j2 all install
 
+compile_Marc:
+  stage: compile
+  script:
+    - module load $IntelMarc $HDF5Marc $MSC
+    - cd $(mktemp -d)
+    - cp ${CI_PROJECT_DIR}/examples/Marc/* .
+    - python3 -c "import damask;damask.solver.Marc().submit_job('r-value','texture',True,'h')"
+    - mkdir ${TESTROOT}/src
+    - mv ${CI_PROJECT_DIR}/src/DAMASK_Marc.marc ${TESTROOT}/src
+
 
 ###################################################################################################
-core:
+open-source:
   stage: fortran
   script:
-    - module load $IntelCompiler $MPI_Intel $PETSc_Intel
+    - module load ${COMPILER_INTEL} ${MPI_INTEL} ${PETSC_INTEL}
-    - cd pytest
+    - cd PRIVATE/testing/pytest
-    - pytest -k 'not compile' --basetemp ${TESTROOT}/fortran -v
+    - pytest -k 'not compile and not Marc' --basetemp ${TESTROOT}/open-source -v
+
+Marc:
+  stage: fortran
+  script:
+    - cd PRIVATE/testing/pytest
+    - pytest -k 'not compile and Marc' --damask-root=${TESTROOT} --basetemp ${TESTROOT}/Marc -v
 
 # Needs closer look
 # Phenopowerlaw_singleSlip:
@@ -156,20 +163,25 @@ core:
 grid_runtime:
   stage: performance
   script:
-    - module load $IntelCompiler $MPI_Intel $PETSc_Intel
+    - module load ${COMPILER_INTEL} ${MPI_INTEL} ${PETSC_INTEL}
     - cd $(mktemp -d)
-    - cmake -DOPTIMIZATION=AGGRESSIVE -DDAMASK_SOLVER=GRID -DCMAKE_INSTALL_PREFIX=${TESTROOT} ${CI_PROJECT_DIR}
+    - cmake -DOPTIMIZATION=AGGRESSIVE -DDAMASK_SOLVER=GRID -DCMAKE_INSTALL_PREFIX=./ ${CI_PROJECT_DIR}
     - make -j2 all install
+    - export PATH=${PWD}/bin:${PATH}
     - cd $(mktemp -d)
     - git clone -q git@git.damask.mpie.de:damask/performance.git .
-    - ${DAMASKROOT}/PRIVATE/testing/runtime.py --input_dir $DAMASKROOT/examples/grid --output_dir . --tag ${CI_COMMIT_SHA}
+    - >
+        ${CI_PROJECT_DIR}/PRIVATE/testing/runtime.py
+        --input_dir ${CI_PROJECT_DIR}/examples/grid
+        --output_dir ./
+        --tag ${CI_COMMIT_SHA}
     - if [ ${CI_COMMIT_BRANCH} == development ]; then git commit -am ${CI_PIPELINE_ID}_${CI_COMMIT_SHA}; git push; fi
   before_script:
     - ${LOCAL_HOME}/bin/queue ${CI_JOB_ID} --blocking
-    - source $DAMASKROOT/env/DAMASK.sh
+    - source env/DAMASK.sh
-    - export PATH=${TESTROOT}/bin:$PATH
     - echo Job start:" $(date)"
 
+
 ###################################################################################################
 source_distribution:
   stage: deploy
@@ -179,16 +191,19 @@ source_distribution:
 
 
 ###################################################################################################
-merge_into_master:
+update_revision:
-  stage: update_master
+  stage: finalize
+  before_script:
+    - ${LOCAL_HOME}/bin/queue ${CI_JOB_ID}
+    - echo Job start:" $(date)"
   script:
-    - cd ${DAMASKROOT}
+    - cd $(mktemp -d)
-    - export TESTEDREV=$(git describe)                                                              # might be detached from development branch
+    - git clone -q git@git.damask.mpie.de:damask/DAMASK.git .
-    - echo ${TESTEDREV} > python/damask/VERSION
+    - git checkout ${CI_COMMIT_SHA}
+    - export VERSION=$(git describe)
+    - echo ${VERSION} > python/damask/VERSION
     - git add python/damask/VERSION
-    - >
+    - git commit -m "[skip ci] updated version information after successful test of $VERSION"
-        git diff-index --quiet HEAD ||
-        git commit -m "[skip ci] updated version information after successful test of $TESTEDREV"
     - export UPDATEDREV=$(git describe)                                                             # tested state + 1 commit
     - git checkout master
     - git pull

PRIVATE

@@ -1 +1 @@
-Subproject commit f730bcb8ddd7224011e70c57d0a5c03068532d2d
+Subproject commit f177ee4d1e8619c41756d8e06df4d7054d2526da

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

@@ -6,7 +6,7 @@ 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
 
 xi_0_sl:   [10.e+6,  55.e+6, 0.,  60.e+6, 0.,  60.e+6]
@@ -23,7 +23,16 @@ 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_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

install/MarcMentat/apply_DAMASK_modifications.py

@@ -8,7 +8,7 @@ from pathlib import Path
 
 import damask
 
-def copy_and_patch(patch,orig,msc_root,editor):
+def copy_and_patch(patch,orig,marc_root,editor):
     try:
         shutil.copyfile(orig,orig.parent/patch.stem)
     except shutil.SameFileError:
@@ -17,31 +17,31 @@ def copy_and_patch(patch,orig,msc_root,editor):
     with open(orig.parent/patch.stem) as f_in:
         content = f_in.read()
     with open(orig.parent/patch.stem,'w') as f_out:
-        f_out.write(content.replace('%INSTALLDIR%',msc_root).replace('%EDITOR%',editor))
+        f_out.write(content.replace('%INSTALLDIR%',marc_root).replace('%EDITOR%',editor))
 
 
 parser = argparse.ArgumentParser(
-                  description='Apply DAMASK modification to MSC.Marc/Mentat',
+                  description='Apply DAMASK modification to MSC Marc/Mentat',
                   prog = Path(__file__).name,
                   formatter_class=argparse.ArgumentDefaultsHelpFormatter)
 
 parser.add_argument('--editor', dest='editor', metavar='string', default='vi',
-                    help='Name of the editor for MSC.Mentat (executable)')
+                    help='Name of the editor for Marc Mentat (executable)')
-parser.add_argument('--msc-root', dest='msc_root', metavar='string',
+parser.add_argument('--marc-root', dest='marc_root', metavar='string',
-                    default=damask.solver._marc._msc_root,
+                    default=damask.solver._marc._marc_root,
-                    help='MSC.Marc/Mentat root directory')
+                    help='Marc root directory')
-parser.add_argument('--msc-version', dest='msc_version', type=float, metavar='float',
+parser.add_argument('--marc-version', dest='marc_version', type=float, metavar='float',
-                    default=damask.solver._marc._msc_version,
+                    default=damask.solver._marc._marc_version,
-                    help='MSC.Marc/Mentat version')
+                    help='Marc version')
 parser.add_argument('--damask-root', dest='damask_root', metavar = 'string',
                     default=damask.solver._marc._damask_root,
                     help='DAMASK root directory')
 
 args = parser.parse_args()
-msc_root    = Path(args.msc_root).expanduser()
+marc_root = Path(args.marc_root).expanduser()
 damask_root = Path(args.damask_root).expanduser()
-msc_version = int(args.msc_version) if str(args.msc_version).split('.')[1] == '0' else \
+marc_version = int(args.marc_version) if str(args.marc_version).split('.')[1] == '0' else \
-              args.msc_version
+               args.marc_version
 
 matches = {'Marc_tools':  [['comp_user','comp_damask_*mp'],
                            ['run_marc','run_damask_*mp'],
@@ -54,15 +54,15 @@ matches = {'Marc_tools':  [['comp_user','comp_damask_*mp'],
 
 print('patching files...\n')
 
-for directory in glob.glob(str(damask_root/f'install/MarcMentat/{msc_version}/*')):
+for directory in glob.glob(str(damask_root/f'install/MarcMentat/{marc_version}/*')):
     for orig, mods in matches[Path(directory).name]:
-        product,subfolder = (msc_root/Path(directory)).name.split('_')
+        product,subfolder = (marc_root/Path(directory)).name.split('_')
-        orig = msc_root/f'{product.lower()}{msc_version}/{subfolder}/{orig}'
+        orig = marc_root/f'{product.lower()}{marc_version}/{subfolder}/{orig}'
         for patch in glob.glob(f'{directory}/{mods}.patch'):
-            copy_and_patch(Path(patch),orig,msc_root,args.editor)
+            copy_and_patch(Path(patch),orig,marc_root,args.editor)
 
 print('compiling Mentat menu binaries...')
 
-executable = msc_root/f'mentat{msc_version}/bin/mentat'
+executable = marc_root/f'mentat{marc_version}/bin/mentat'
-menu_file  = msc_root/f'mentat{msc_version}/menus/linux64/main.msb'
+menu_file  = marc_root/f'mentat{marc_version}/menus/linux64/main.msb'
 os.system(f'xvfb-run -a {executable} -compile {menu_file}')

python/damask/VERSION

@@ -1 +1 @@
-v3.0.0-alpha5-90-gbb6655045
+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/damask/_configmaterial.py

@@ -440,9 +440,10 @@ class ConfigMaterial(Config):
         """
         N,n,shaped = 1,1,{}
 
+        map_dim = {'O':-1,'F_i':-2}
         for k,v in kwargs.items():
             shaped[k] = np.array(v)
-            s = shaped[k].shape[:-1] if k=='O' else shaped[k].shape
+            s = shaped[k].shape[:map_dim.get(k,None)]
             N = max(N,s[0]) if len(s)>0 else N
             n = max(n,s[1]) if len(s)>1 else n
 
@@ -451,11 +452,12 @@ class ConfigMaterial(Config):
         if 'v' not in kwargs:
             shaped['v'] = np.broadcast_to(1/n,(N,n))
 
+        map_shape = {'O':(N,n,4),'F_i':(N,n,3,3)}
         for k,v in shaped.items():
-            target = (N,n,4) if k=='O' else (N,n)
+            target = map_shape.get(k,(N,n))
             obj = np.broadcast_to(v.reshape(util.shapeshifter(v.shape,target,mode='right')),target)
             for i in range(N):
-                if k in ['phase','O','v']:
+                if k in ['phase','O','v','F_i']:
                     for j in range(n):
                         mat[i]['constituents'][j][k] = obj[i,j].item() if isinstance(obj[i,j],np.generic) else obj[i,j]
                 else:

python/damask/_grid.py

@@ -54,9 +54,9 @@ class Grid:
         mat_max = np.nanmax(self.material)
         mat_N   = self.N_materials
         return util.srepr([
-               f'cells  a b c: {util.srepr(self.cells, " x ")}',
+               f'cells : {util.srepr(self.cells, " x ")}',
-               f'size   x y z: {util.srepr(self.size,  " x ")}',
+               f'size  : {util.srepr(self.size,  " x ")} / m³',
-               f'origin x y z: {util.srepr(self.origin,"   ")}',
+               f'origin: {util.srepr(self.origin,"   ")} / m',
                f'# materials: {mat_N}' + ('' if mat_min == 0 and mat_max+1 == mat_N else
                                           f' (min: {mat_min}, max: {mat_max})')
               ])

python/damask/solver/_marc.py

@@ -3,14 +3,14 @@ import shlex
 import re
 from pathlib import Path
 
-_msc_version = 2021.2
+_marc_version = 2021.2
-_msc_root = '/opt/msc'
+_marc_root = '/opt/msc'
 _damask_root = str(Path(__file__).parents[3])
 
 class Marc:
-    """Wrapper to run DAMASK with MSC.Marc."""
+    """Wrapper to run DAMASK with MSC Marc."""
 
-    def __init__(self,msc_version=_msc_version,msc_root=_msc_root,damask_root=_damask_root):
+    def __init__(self,marc_version=_marc_version,marc_root=_marc_root,damask_root=_damask_root):
         """
         Create a Marc solver object.
 
@@ -20,14 +20,14 @@ class Marc:
             Marc version
 
         """
-        self.msc_version = msc_version
+        self.marc_version = marc_version
-        self.msc_root    = Path(msc_root)
+        self.marc_root    = Path(marc_root)
         self.damask_root = Path(damask_root)
 
     @property
     def library_path(self):
 
-        path_lib = self.msc_root/f'mentat{self.msc_version}/shlib/linux64'
+        path_lib = self.marc_root/f'mentat{self.marc_version}/shlib/linux64'
         if not path_lib.is_dir():
             raise FileNotFoundError(f'library path "{path_lib}" not found')
 
@@ -37,7 +37,7 @@ class Marc:
     @property
     def tools_path(self):
 
-        path_tools = self.msc_root/f'marc{self.msc_version}/tools'
+        path_tools = self.marc_root/f'marc{self.marc_version}/tools'
         if not path_tools.is_dir():
             raise FileNotFoundError(f'tools path "{path_tools}" not found')
 

python/tests/conftest.py

@@ -49,9 +49,8 @@ def patch_plt_show(monkeypatch):
 
 
 def pytest_addoption(parser):
-    parser.addoption("--update",
+    parser.addoption('--update', action='store_true', default=False,
-                     action="store_true",
+                     help='Update reference results.')
-                     default=False)
 
 
 @pytest.fixture

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

python/tests/test_ConfigMaterial.py

@@ -99,12 +99,15 @@ class TestConfigMaterial:
     @pytest.mark.parametrize('N,n,kw',[
                                         (1,1,{'phase':'Gold',
                                               'O':[1,0,0,0],
+                                              'F_i':np.eye(3),
                                               'homogenization':'SX'}),
                                         (3,1,{'phase':'Gold',
                                               'O':Rotation.from_random(3),
+                                              'F_i':np.broadcast_to(np.eye(3),(3,3,3)),
                                               'homogenization':'SX'}),
                                         (2,3,{'phase':np.broadcast_to(['a','b','c'],(2,3)),
                                               'O':Rotation.from_random((2,3)),
+                                              'F_i':np.broadcast_to(np.eye(3),(2,3,3,3)),
                                               'homogenization':['SX','PX']}),
                                         ])
     def test_material_add(self,kw,N,n):

python/tests/test_Rotation.py

@@ -989,11 +989,19 @@ class TestRotation:
         with pytest.raises(TypeError):
             R@data
 
-    def test_misorientation(self):
+    def test_misorientation_invariant(self):
         R = Rotation.from_random()
         assert np.allclose(R.misorientation(R).as_matrix(),np.eye(3))
 
-    def test_misorientation360(self):
+    def test_misorientation_average(self):
+        """2 times the average is the misorientation."""
+        r = Rotation.from_random(2)
+        a = r[0].misorientation(r[1]).as_axis_angle()
+        b = r.average().misorientation(r[1]).as_axis_angle()
+        b[3] = (b[3]*2)%np.pi
+        assert np.allclose(a,b)
+
+    def test_misorientation_360deg(self):
         R_1 = Rotation()
         R_2 = Rotation.from_Euler_angles([360,0,0],degrees=True)
         assert np.allclose(R_1.misorientation(R_2).as_matrix(),np.eye(3))

src/CPFEM.f90

@@ -103,7 +103,7 @@ subroutine CPFEM_init
   class(tNode), pointer :: &
     debug_CPFEM
 
-  print'(/,a)', ' <<<+-  CPFEM init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  CPFEM init  -+>>>'; flush(IO_STDOUT)
 
   allocate(CPFEM_cs(               6,discretization_nIPs,discretization_Nelems), source= 0.0_pReal)
   allocate(CPFEM_dcsdE(          6,6,discretization_nIPs,discretization_Nelems), source= 0.0_pReal)

src/CPFEM2.f90

@@ -81,7 +81,7 @@ subroutine CPFEM_init
   integer(HID_T) :: fileHandle
 
 
-  print'(/,a)', ' <<<+-  CPFEM init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  CPFEM init  -+>>>'; flush(IO_STDOUT)
 
 
   if (interface_restartInc > 0) then

src/DAMASK_Marc.f90

@@ -46,7 +46,7 @@ subroutine DAMASK_interface_init
   integer                 :: ierr
   character(len=pPathLen) :: wd
 
-  print'(/,a)', ' <<<+-  DAMASK_marc init -+>>>'
+  print'(/,1x,a)', '<<<+-  DAMASK_marc init -+>>>'
 
   print*, 'Roters et al., Computational Materials Science 158:420–478, 2019'
   print*, 'https://doi.org/10.1016/j.commatsci.2018.04.030'

src/DAMASK_interface.f90

@@ -70,7 +70,7 @@ subroutine DAMASK_interface_init
   external :: &
     quit
 
-  print'(/,a)', ' <<<+-  DAMASK_interface init  -+>>>'
+  print'(/,1x,a)', '<<<+-  DAMASK_interface init  -+>>>'
 
   if(worldrank == 0) open(OUTPUT_UNIT, encoding='UTF-8')                                            ! for special characters in output
 

src/HDF5_utilities.f90

@@ -109,7 +109,7 @@ subroutine HDF5_utilities_init
   integer(SIZE_T) :: typeSize
 
 
-  print'(/,a)', ' <<<+-  HDF5_Utilities init  -+>>>'
+  print'(/,1x,a)', '<<<+-  HDF5_Utilities init  -+>>>'
 
 
   call h5open_f(hdferr)

src/IO.f90

@@ -56,7 +56,7 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine IO_init
 
-  print'(/,a)', ' <<<+-  IO init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  IO init  -+>>>'; flush(IO_STDOUT)
 
   call selfTest
 

src/YAML_types.f90

@@ -187,7 +187,7 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine YAML_types_init
 
-  print'(/,a)', ' <<<+-  YAML_types init  -+>>>'
+  print'(/,1x,a)', '<<<+-  YAML_types init  -+>>>'
 
   call selfTest
 

src/base64.f90

@@ -27,7 +27,7 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine base64_init
 
-  print'(/,a)', ' <<<+-  base64 init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  base64 init  -+>>>'; flush(IO_STDOUT)
 
   call selfTest
 

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/config.f90

@@ -30,8 +30,7 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine config_init
 
-  print'(/,a)', ' <<<+-  config init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  config init  -+>>>'; flush(IO_STDOUT)
-
 
   call parse_material
   call parse_numerics
@@ -53,7 +52,7 @@ subroutine parse_material()
   if (.not. fileExists) call IO_error(100,ext_msg='material.yaml')
 
   if (worldrank == 0) then
-    print*, 'reading material.yaml'; flush(IO_STDOUT)
+    print'(/,1x,a)', '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')
@@ -81,7 +80,7 @@ subroutine parse_numerics()
   if (fileExists) then
 
     if (worldrank == 0) then
-      print*, 'reading numerics.yaml'; flush(IO_STDOUT)
+      print'(1x,a)', 'reading numerics.yaml'; flush(IO_STDOUT)
       fileContent = IO_read('numerics.yaml')
       if (len(fileContent) > 0) then
         call results_openJobFile(parallel=.false.)
@@ -113,7 +112,7 @@ subroutine parse_debug()
   if (fileExists) then
 
     if (worldrank == 0) then
-      print*, 'reading debug.yaml'; flush(IO_STDOUT)
+      print'(1x,a)', 'reading debug.yaml'; flush(IO_STDOUT)
       fileContent = IO_read('debug.yaml')
       if (len(fileContent) > 0) then
         call results_openJobFile(parallel=.false.)

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/discretization.f90

@@ -49,7 +49,7 @@ subroutine discretization_init(materialAt,&
   integer, optional,           intent(in) :: &
     sharedNodesBegin                                                                                !< index of first node shared among different processes (MPI)
 
-  print'(/,a)', ' <<<+-  discretization init  -+>>>'; flush(6)
+  print'(/,1x,a)', '<<<+-  discretization init  -+>>>'; flush(6)
 
   discretization_Nelems = size(materialAt,1)
   discretization_nIPs   = size(IPcoords0,2)/discretization_Nelems

src/element.f90

@@ -923,7 +923,7 @@ subroutine tElement_init(self,elemType)
 
   self%nIPneighbors = size(self%IPneighbor,1)
 
-  print'(/,a)', ' <<<+-  element_init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  element_init  -+>>>'; flush(IO_STDOUT)
 
   print*, 'element type:      ',self%elemType
   print*, '  geom type:       ',self%geomType

src/grid/DAMASK_grid.f90

@@ -113,10 +113,10 @@ program DAMASK_grid
 ! init DAMASK (all modules)
 
   call CPFEM_initAll
-  print'(/,a)',   ' <<<+-  DAMASK_grid init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  DAMASK_grid init  -+>>>'; flush(IO_STDOUT)
 
-  print*, 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019'
+  print'(/,1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019'
-  print*, 'https://doi.org/10.1007/978-981-10-6855-3_80'
+  print'(  1x,a)', 'https://doi.org/10.1007/978-981-10-6855-3_80'
 
 
 !-------------------------------------------------------------------------------------------------
@@ -223,7 +223,7 @@ program DAMASK_grid
     loadCases(l)%r         = step_discretization%get_asFloat('r',         defaultVal= 1.0_pReal)
 
     loadCases(l)%f_restart = load_step%get_asInt('f_restart', defaultVal=huge(0))
-    if (load_step%get_asString('f_out',defaultVal='DAMASK') == 'none') then
+    if (load_step%get_asString('f_out',defaultVal='n/a') == 'none') then
       loadCases(l)%f_out = huge(0)
     else
       loadCases(l)%f_out = load_step%get_asInt('f_out', defaultVal=1)
@@ -231,12 +231,12 @@ program DAMASK_grid
     loadCases(l)%estimate_rate = (load_step%get_asBool('estimate_rate',defaultVal=.true.) .and. l>1)
 
     reportAndCheck: if (worldrank == 0) then
-      print'(/,a,i0)', ' load case: ', l
+      print'(/,1x,a,1x,i0)', 'load case:', l
-      print*, ' estimate_rate:', loadCases(l)%estimate_rate
+      print'(2x,a,1x,l1)', 'estimate_rate:', loadCases(l)%estimate_rate
       if (loadCases(l)%deformation%myType == 'F') then
-        print*, ' F:'
+        print'(2x,a)', 'F:'
       else
-        print*, ' '//loadCases(l)%deformation%myType//' / 1/s:'
+        print'(2x,a)', loadCases(l)%deformation%myType//' / 1/s:'
       end if
       do i = 1, 3; do j = 1, 3
         if (loadCases(l)%deformation%mask(i,j)) then
@@ -250,8 +250,8 @@ program DAMASK_grid
       if (any(.not.(loadCases(l)%stress%mask .or. transpose(loadCases(l)%stress%mask)) .and. (math_I3<1))) &
         errorID = 838                                                                               ! no rotation is allowed by stress BC
 
-      if (loadCases(l)%stress%myType == 'P')     print*, ' P / MPa:'
+      if (loadCases(l)%stress%myType == 'P')     print'(2x,a)', 'P / MPa:'
-      if (loadCases(l)%stress%myType == 'dot_P') print*, ' dot_P / MPa/s:'
+      if (loadCases(l)%stress%myType == 'dot_P') print'(2x,a)', 'dot_P / MPa/s:'
 
       if (loadCases(l)%stress%myType /= '') then
         do i = 1, 3; do j = 1, 3
@@ -274,16 +274,16 @@ program DAMASK_grid
       if (loadCases(l)%f_restart < 1)  errorID = 839
 
       if (dEq(loadCases(l)%r,1.0_pReal,1.e-9_pReal)) then
-        print'(a)', '  r: 1 (constant step width)'
+        print'(2x,a)', 'r: 1 (constant step width)'
       else
-        print'(a,f0.3)', '  r: ', loadCases(l)%r
+        print'(2x,a,1x,f0.3)', 'r:', loadCases(l)%r
       endif
-      print'(a,f0.3)',   '  t: ', loadCases(l)%t
+      print'(2x,a,1x,f0.3)',   't:', loadCases(l)%t
-      print'(a,i0)',     '  N: ', loadCases(l)%N
+      print'(2x,a,1x,i0)',     'N:', loadCases(l)%N
       if (loadCases(l)%f_out < huge(0)) &
-        print'(a,i0)',   '  f_out: ', loadCases(l)%f_out
+        print'(2x,a,1x,i0)',   'f_out:', loadCases(l)%f_out
       if (loadCases(l)%f_restart < huge(0)) &
-        print'(a,i0)',   '  f_restart: ', loadCases(l)%f_restart
+        print'(2x,a,1x,i0)',   'f_restart:', loadCases(l)%f_restart
 
       if (errorID > 0) call IO_error(error_ID = errorID, el = l)
 
@@ -322,7 +322,7 @@ program DAMASK_grid
   endif
 
   writeUndeformed: if (interface_restartInc < 1) then
-    print'(/,a)', ' ... writing initial configuration to file ........................'
+    print'(/,1x,a)', '... writing initial configuration to file .................................'
     flush(IO_STDOUT)
     call CPFEM_results(0,0.0_pReal)
   endif writeUndeformed
@@ -358,8 +358,8 @@ program DAMASK_grid
 
 !--------------------------------------------------------------------------------------------------
 ! report begin of new step
-          print'(/,a)', ' ###########################################################################'
+          print'(/,1x,a)', '###########################################################################'
-          print'(1x,a,es12.5,6(a,i0))', &
+          print'(1x,a,1x,es12.5,6(a,i0))', &
                   'Time', t, &
                   's: Increment ', inc,'/',loadCases(l)%N,&
                   '-', stepFraction,'/',subStepFactor**cutBackLevel,&
@@ -430,7 +430,7 @@ program DAMASK_grid
             cutBackLevel = cutBackLevel + 1
             t = t - Delta_t
             Delta_t = Delta_t/real(subStepFactor,pReal)                                             ! cut timestep
-            print'(/,a)', ' cutting back '
+            print'(/,1x,a)', 'cutting back '
           else                                                                                      ! no more options to continue
             if (worldrank == 0) close(statUnit)
             call IO_error(950)
@@ -441,15 +441,15 @@ program DAMASK_grid
         cutBackLevel = max(0, cutBackLevel - 1)                                                     ! try half number of subincs next inc
 
         if (all(solres(:)%converged)) then
-          print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' converged'
+          print'(/,1x,a,i0,a)', 'increment ', totalIncsCounter, ' converged'
         else
-          print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' NOT converged'
+          print'(/,1x,a,i0,a)', 'increment ', totalIncsCounter, ' NOT converged'
         endif; flush(IO_STDOUT)
 
         call MPI_Allreduce(interface_SIGUSR1,signal,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr)
         if (ierr /= 0) error stop 'MPI error'
         if (mod(inc,loadCases(l)%f_out) == 0 .or. signal) then
-          print'(1/,a)', ' ... writing results to file ......................................'
+          print'(/,1x,a)', '... writing results to file ...............................................'
           flush(IO_STDOUT)
           call CPFEM_results(totalIncsCounter,t)
         endif
@@ -473,7 +473,7 @@ program DAMASK_grid
 
 !--------------------------------------------------------------------------------------------------
 ! report summary of whole calculation
-  print'(/,a)', ' ###########################################################################'
+  print'(/,1x,a)', '###########################################################################'
   if (worldrank == 0) close(statUnit)
 
   call quit(0)                                                                                      ! no complains ;)

src/grid/discretization_grid.f90

@@ -72,7 +72,7 @@ subroutine discretization_grid_init(restart)
     fileContent, fname
 
 
-  print'(/,a)', ' <<<+-  discretization_grid init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  discretization_grid init  -+>>>'; flush(IO_STDOUT)
 
 
   if (worldrank == 0) then
@@ -96,9 +96,9 @@ subroutine discretization_grid_init(restart)
   call MPI_Bcast(origin,3,MPI_DOUBLE,0,MPI_COMM_WORLD, ierr)
   if (ierr /= 0) error stop 'MPI error'
 
-  print'(/,a,3(i12  ))',  ' cells  a b c: ', grid
+  print'(/,1x,a,3(i12,1x))',    'cells  a b c: ', grid
-  print'(a,3(es12.5))',   ' size   x y z: ', geomSize
+  print  '(1x,a,3(es12.5,1x))', 'size   x y z: ', geomSize
-  print'(a,3(es12.5))',   ' origin x y z: ', origin
+  print  '(1x,a,3(es12.5,1x))', 'origin x y z: ', origin
 
   if (worldsize>grid(3)) call IO_error(894, ext_msg='number of processes exceeds grid(3)')
 

src/grid/grid_damage_spectral.f90

@@ -76,10 +76,10 @@ subroutine grid_damage_spectral_init()
   character(len=pStringLen) :: &
     snes_type
 
-  print'(/,a)', ' <<<+-  grid_spectral_damage init  -+>>>'
+  print'(/,1x,a)', '<<<+-  grid_spectral_damage init  -+>>>'
 
-  print*, 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019'
+  print'(/,1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019'
-  print*, 'https://doi.org/10.1007/978-981-10-6855-3_80'
+  print'(  1x,a)', 'https://doi.org/10.1007/978-981-10-6855-3_80'
 
 !-------------------------------------------------------------------------------------------------
 ! read numerical parameters and do sanity checks
@@ -206,9 +206,9 @@ function grid_damage_spectral_solution(Delta_t) result(solution)
   call VecMin(solution_vec,devNull,phi_min,ierr); CHKERRQ(ierr)
   call VecMax(solution_vec,devNull,phi_max,ierr); CHKERRQ(ierr)
   if (solution%converged) &
-    print'(/,a)', ' ... nonlocal damage converged .....................................'
+    print'(/,1x,a)', '... nonlocal damage converged .....................................'
-  print'(/,a,f8.6,2x,f8.6,2x,e11.4)', ' Minimum|Maximum|Delta Damage      = ', phi_min, phi_max, stagNorm
+  print'(/,1x,a,f8.6,2x,f8.6,2x,e11.4)', 'Minimum|Maximum|Delta Damage      = ', phi_min, phi_max, stagNorm
-  print'(/,a)', ' ==========================================================================='
+  print'(/,1x,a)', '==========================================================================='
   flush(IO_STDOUT)
 
 end function grid_damage_spectral_solution

src/grid/grid_mech_FEM.f90

@@ -116,7 +116,7 @@ subroutine grid_mechanical_FEM_init
     num_grid, &
     debug_grid
 
-  print'(/,a)', ' <<<+-  grid_mechanical_FEM init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  grid_mechanical_FEM init  -+>>>'; flush(IO_STDOUT)
 
 !-------------------------------------------------------------------------------------------------
 ! debugging options
@@ -234,7 +234,7 @@ subroutine grid_mechanical_FEM_init
 !--------------------------------------------------------------------------------------------------
 ! init fields
   restartRead: if (interface_restartInc > 0) then
-    print'(/,a,i0,a)', ' reading restart data of increment ', interface_restartInc, ' from file'
+    print'(/,1x,a,i0,a)', 'reading restart data of increment ', interface_restartInc, ' from file'
 
     fileHandle  = HDF5_openFile(getSolverJobName()//'_restart.hdf5','r')
     groupHandle = HDF5_openGroup(fileHandle,'solver')
@@ -272,7 +272,7 @@ subroutine grid_mechanical_FEM_init
   CHKERRQ(ierr)
 
   restartRead2: if (interface_restartInc > 0) then
-    print'(a,i0,a)', ' reading more restart data of increment ', interface_restartInc, ' from file'
+    print'(1x,a,i0,a)', 'reading more restart data of increment ', interface_restartInc, ' from file'
     call HDF5_read(C_volAvg,groupHandle,'C_volAvg',.false.)
     call MPI_Bcast(C_volAvg,81,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
     if(ierr /=0) error stop 'MPI error'
@@ -442,7 +442,7 @@ subroutine grid_mechanical_FEM_restartWrite
   call DMDAVecGetArrayF90(mechanical_grid,solution_lastInc,u_lastInc,ierr)
   CHKERRQ(ierr)
 
-  print*, 'writing solver data required for restart to file'; flush(IO_STDOUT)
+  print'(1x,a)', 'writing solver data required for restart to file'; flush(IO_STDOUT)
 
   fileHandle  = HDF5_openFile(getSolverJobName()//'_restart.hdf5','w')
   groupHandle = HDF5_addGroup(fileHandle,'solver')
@@ -506,12 +506,12 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,i
     reason = 0
   endif
 
-  print'(1/,a)', ' ... reporting .............................................................'
+  print'(/,1x,a)', '... reporting .............................................................'
-  print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', &
+  print'(/,1x,a,f12.2,a,es8.2,a,es9.2,a)', 'error divergence = ', &
           err_div/divTol,  ' (',err_div,' / m, tol = ',divTol,')'
-  print'(a,f12.2,a,es8.2,a,es9.2,a)',    ' error stress BC  = ', &
+  print'(1x,a,f12.2,a,es8.2,a,es9.2,a)',    'error stress BC  = ', &
           err_BC/BCTol,    ' (',err_BC, ' Pa,  tol = ',BCTol,')'
-  print'(/,a)', ' ==========================================================================='
+  print'(/,1x,a)', '==========================================================================='
   flush(IO_STDOUT)
 
 end subroutine converged
@@ -547,9 +547,9 @@ subroutine formResidual(da_local,x_local, &
   newIteration: if (totalIter <= PETScIter) then
     totalIter = totalIter + 1
     print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter+1, '≤', num%itmax
-    if (debugRotation) print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
+    if (debugRotation) print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
       'deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
-    print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
+    print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
       'deformation gradient aim       =', transpose(F_aim)
     flush(IO_STDOUT)
   endif newIteration

src/grid/grid_mech_spectral_basic.f90

@@ -111,13 +111,13 @@ subroutine grid_mechanical_spectral_basic_init
     num_grid, &
     debug_grid
 
-  print'(/,a)', ' <<<+-  grid_mechanical_spectral_basic init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  grid_mechanical_spectral_basic init  -+>>>'; flush(IO_STDOUT)
 
-  print*, 'P. Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013'
+  print'(/,1x,a)', 'P. Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013'
-  print*, 'https://doi.org/10.1016/j.ijplas.2012.09.012'//IO_EOL
+  print'(  1x,a)', 'https://doi.org/10.1016/j.ijplas.2012.09.012'//IO_EOL
 
-  print*, 'P. Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
+  print'(  1x,a)', 'P. Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
-  print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006'
+  print'(  1x,a)', 'https://doi.org/10.1016/j.ijplas.2014.02.006'
 
 !-------------------------------------------------------------------------------------------------
 ! debugging options
@@ -186,7 +186,7 @@ subroutine grid_mechanical_spectral_basic_init
   call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)                                   ! places pointer on PETSc data
 
   restartRead: if (interface_restartInc > 0) then
-    print'(/,a,i0,a)', ' reading restart data of increment ', interface_restartInc, ' from file'
+    print'(/,1x,a,i0,a)', 'reading restart data of increment ', interface_restartInc, ' from file'
 
     fileHandle  = HDF5_openFile(getSolverJobName()//'_restart.hdf5','r')
     groupHandle = HDF5_openGroup(fileHandle,'solver')
@@ -219,7 +219,7 @@ subroutine grid_mechanical_spectral_basic_init
   call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)                                ! deassociate pointer
 
   restartRead2: if (interface_restartInc > 0) then
-    print'(a,i0,a)', ' reading more restart data of increment ', interface_restartInc, ' from file'
+    print'(1x,a,i0,a)', 'reading more restart data of increment ', interface_restartInc, ' from file'
     call HDF5_read(C_volAvg,groupHandle,'C_volAvg',.false.)
     call MPI_Bcast(C_volAvg,81,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
     if (ierr /=0) error stop 'MPI error'
@@ -385,7 +385,7 @@ subroutine grid_mechanical_spectral_basic_restartWrite
 
   call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
 
-  print*, 'writing solver data required for restart to file'; flush(IO_STDOUT)
+  print'(1x,a)', 'writing solver data required for restart to file'; flush(IO_STDOUT)
 
   fileHandle  = HDF5_openFile(getSolverJobName()//'_restart.hdf5','w')
   groupHandle = HDF5_addGroup(fileHandle,'solver')
@@ -445,12 +445,12 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
     reason = 0
   end if
 
-  print'(1/,a)', ' ... reporting .............................................................'
+  print'(/,1x,a)', '... reporting .............................................................'
-  print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', &
+  print'(/,1x,a,f12.2,a,es8.2,a,es9.2,a)', 'error divergence = ', &
           err_div/divTol,  ' (',err_div,' / m, tol = ',divTol,')'
-  print'(a,f12.2,a,es8.2,a,es9.2,a)',    ' error stress BC  = ', &
+  print'(1x,a,f12.2,a,es8.2,a,es9.2,a)',    'error stress BC  = ', &
           err_BC/BCTol,    ' (',err_BC, ' Pa,  tol = ',BCTol,')'
-  print'(/,a)', ' ==========================================================================='
+  print'(/,1x,a)', '==========================================================================='
   flush(IO_STDOUT)
 
 end subroutine converged
@@ -485,9 +485,9 @@ subroutine formResidual(in, F, &
   newIteration: if (totalIter <= PETScIter) then
     totalIter = totalIter + 1
     print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax
-    if (debugRotation) print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
+    if (debugRotation) print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
       'deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
-    print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
+    print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
       'deformation gradient aim       =', transpose(F_aim)
     flush(IO_STDOUT)
   end if newIteration

src/grid/grid_mech_spectral_polarisation.f90

@@ -124,10 +124,10 @@ subroutine grid_mechanical_spectral_polarisation_init
     num_grid, &
     debug_grid
 
-  print'(/,a)', ' <<<+-  grid_mechanical_spectral_polarization init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  grid_mechanical_spectral_polarization init  -+>>>'; flush(IO_STDOUT)
 
-  print*, 'P. Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
+  print'(/,1x,a)', 'P. Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
-  print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006'
+  print'(  1x,a)', 'https://doi.org/10.1016/j.ijplas.2014.02.006'
 
 !-------------------------------------------------------------------------------------------------
 ! debugging options
@@ -208,7 +208,7 @@ subroutine grid_mechanical_spectral_polarisation_init
   F_tau => FandF_tau(9:17,:,:,:)
 
   restartRead: if (interface_restartInc > 0) then
-    print'(/,a,i0,a)', ' reading restart data of increment ', interface_restartInc, ' from file'
+    print'(/,1x,a,i0,a)', 'reading restart data of increment ', interface_restartInc, ' from file'
 
     fileHandle  = HDF5_openFile(getSolverJobName()//'_restart.hdf5','r')
     groupHandle = HDF5_openGroup(fileHandle,'solver')
@@ -245,7 +245,7 @@ subroutine grid_mechanical_spectral_polarisation_init
   call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)                        ! deassociate pointer
 
   restartRead2: if (interface_restartInc > 0) then
-    print'(a,i0,a)', ' reading more restart data of increment ', interface_restartInc, ' from file'
+    print'(1x,a,i0,a)', 'reading more restart data of increment ', interface_restartInc, ' from file'
     call HDF5_read(C_volAvg,groupHandle,'C_volAvg',.false.)
     call MPI_Bcast(C_volAvg,81,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
     if (ierr /=0) error stop 'MPI error'
@@ -441,7 +441,7 @@ subroutine grid_mechanical_spectral_polarisation_restartWrite
   F     => FandF_tau(0: 8,:,:,:)
   F_tau => FandF_tau(9:17,:,:,:)
 
-  print*, 'writing solver data required for restart to file'; flush(IO_STDOUT)
+  print'(1x,a)', 'writing solver data required for restart to file'; flush(IO_STDOUT)
 
   fileHandle  = HDF5_openFile(getSolverJobName()//'_restart.hdf5','w')
   groupHandle = HDF5_addGroup(fileHandle,'solver')
@@ -504,14 +504,14 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
     reason = 0
   end if
 
-  print'(1/,a)', ' ... reporting .............................................................'
+  print'(/,1x,a)', '... reporting .............................................................'
-  print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', &
+  print'(/,1x,a,f12.2,a,es8.2,a,es9.2,a)', 'error divergence = ', &
             err_div/divTol,  ' (',err_div, ' / m, tol = ',divTol,')'
-  print  '(a,f12.2,a,es8.2,a,es9.2,a)', ' error curl       = ', &
+  print  '(1x,a,f12.2,a,es8.2,a,es9.2,a)', 'error curl       = ', &
             err_curl/curlTol,' (',err_curl,' -,   tol = ',curlTol,')'
-  print  '(a,f12.2,a,es8.2,a,es9.2,a)', ' error stress BC  = ', &
+  print  '(1x,a,f12.2,a,es8.2,a,es9.2,a)', 'error stress BC  = ', &
             err_BC/BCTol,    ' (',err_BC,  ' Pa,  tol = ',BCTol,')'
-  print'(/,a)', ' ==========================================================================='
+  print'(/,1x,a)', '==========================================================================='
   flush(IO_STDOUT)
 
 end subroutine converged
@@ -565,9 +565,9 @@ subroutine formResidual(in, FandF_tau, &
   newIteration: if (totalIter <= PETScIter) then
     totalIter = totalIter + 1
     print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax
-    if (debugRotation) print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
+    if (debugRotation) print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
       'deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
-    print'(/,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
+    print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
       'deformation gradient aim       =', transpose(F_aim)
     flush(IO_STDOUT)
   end if newIteration

src/grid/grid_thermal_spectral.f90

@@ -75,10 +75,10 @@ subroutine grid_thermal_spectral_init(T_0)
   class(tNode), pointer :: &
     num_grid
 
-  print'(/,a)', ' <<<+-  grid_thermal_spectral init  -+>>>'
+  print'(/,1x,a)', '<<<+-  grid_thermal_spectral init  -+>>>'
 
-  print*, 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019'
+  print'(/,1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019'
-  print*, 'https://doi.org/10.1007/978-981-10-6855-3_80'
+  print'(  1x,a)', 'https://doi.org/10.1007/978-981-10-6855-3_80'
 
 !-------------------------------------------------------------------------------------------------
 ! read numerical parameters and do sanity checks
@@ -201,9 +201,9 @@ function grid_thermal_spectral_solution(Delta_t) result(solution)
   call VecMin(solution_vec,devNull,T_min,ierr); CHKERRQ(ierr)
   call VecMax(solution_vec,devNull,T_max,ierr); CHKERRQ(ierr)
   if (solution%converged) &
-    print'(/,a)', ' ... thermal conduction converged ..................................'
+    print'(/,1x,a)', '... thermal conduction converged ..................................'
-  print'(/,a,f8.4,2x,f8.4,2x,f8.4)', ' Minimum|Maximum|Delta Temperature / K = ', T_min, T_max, stagNorm
+  print'(/,1x,a,f8.4,2x,f8.4,2x,f8.4)', 'Minimum|Maximum|Delta Temperature / K = ', T_min, T_max, stagNorm
-  print'(/,a)', ' ==========================================================================='
+  print'(/,1x,a)', '==========================================================================='
   flush(IO_STDOUT)
 
 end function grid_thermal_spectral_solution

src/grid/spectral_utilities.f90

@@ -177,19 +177,19 @@ subroutine spectral_utilities_init
     num_grid, &
     debug_grid                                                                                      ! pointer to grid  debug options
 
-  print'(/,a)', ' <<<+-  spectral_utilities init  -+>>>'
+  print'(/,1x,a)', '<<<+-  spectral_utilities init  -+>>>'
 
-  print*, 'M. Diehl, Diploma Thesis TU München, 2010'
+  print'(/,1x,a)', 'M. Diehl, Diploma Thesis TU München, 2010'
-  print*, 'https://doi.org/10.13140/2.1.3234.3840'//IO_EOL
+  print'(  1x,a)', 'https://doi.org/10.13140/2.1.3234.3840'//IO_EOL
 
-  print*, 'P. Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013'
+  print'(  1x,a)', 'P. Eisenlohr et al., International Journal of Plasticity 46:37–53, 2013'
-  print*, 'https://doi.org/10.1016/j.ijplas.2012.09.012'//IO_EOL
+  print'(  1x,a)', 'https://doi.org/10.1016/j.ijplas.2012.09.012'//IO_EOL
 
-  print*, 'P. Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
+  print'(  1x,a)', 'P. Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
-  print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006'//IO_EOL
+  print'(  1x,a)', 'https://doi.org/10.1016/j.ijplas.2014.02.006'//IO_EOL
 
-  print*, 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019'
+  print'(  1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019'
-  print*, 'https://doi.org/10.1007/978-981-10-6855-3_80'
+  print'(  1x,a)', 'https://doi.org/10.1007/978-981-10-6855-3_80'
 
 !--------------------------------------------------------------------------------------------------
 ! set debugging parameters
@@ -200,7 +200,7 @@ subroutine spectral_utilities_init
   debugRotation   =  debug_grid%contains('rotation')
   debugPETSc      =  debug_grid%contains('PETSc')
 
-  if(debugPETSc) print'(3(/,a),/)', &
+  if (debugPETSc) print'(3(/,1x,a),/)', &
                  'Initializing PETSc with debug options: ', &
                  trim(PETScDebug), &
                  'add more using the "PETSc_options" keyword in numerics.yaml'
@@ -271,7 +271,7 @@ subroutine spectral_utilities_init
   if (pReal /= C_DOUBLE .or. kind(1) /= C_INT) error stop 'C and Fortran datatypes do not match'
   call fftw_set_timelimit(num_grid%get_asFloat('fftw_timelimit',defaultVal=-1.0_pReal))
 
-  print*, 'FFTW initialized'; flush(IO_STDOUT)
+  print'(/,1x,a)', 'FFTW initialized'; flush(IO_STDOUT)
 
 !--------------------------------------------------------------------------------------------------
 ! MPI allocation
@@ -497,7 +497,7 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
   logical :: err
 
 
-  print'(/,a)', ' ... doing gamma convolution ...............................................'
+  print'(/,1x,a)', '... doing gamma convolution ...............................................'
   flush(IO_STDOUT)
 
 !--------------------------------------------------------------------------------------------------
@@ -567,7 +567,7 @@ real(pReal) function utilities_divergenceRMS()
   integer :: i, j, k, ierr
   complex(pReal), dimension(3)   :: rescaledGeom
 
-  print'(/,a)', ' ... calculating divergence ................................................'
+  print'(/,1x,a)', '... calculating divergence ................................................'
   flush(IO_STDOUT)
 
   rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal)
@@ -610,7 +610,7 @@ real(pReal) function utilities_curlRMS()
   complex(pReal), dimension(3,3) :: curl_fourier
   complex(pReal), dimension(3)   :: rescaledGeom
 
-  print'(/,a)', ' ... calculating curl ......................................................'
+  print'(/,1x,a)', '... calculating curl ......................................................'
   flush(IO_STDOUT)
 
   rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal)
@@ -690,8 +690,8 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
     temp99_real = math_3333to99(rot_BC%rotate(C))
 
     if (debugGeneral) then
-      print'(/,a)', ' ... updating masked compliance ............................................'
+      print'(/,1x,a)', '... updating masked compliance ............................................'
-      print'(/,a,/,8(9(2x,f12.7,1x)/),9(2x,f12.7,1x))', &
+      print'(/,1x,a,/,8(9(2x,f12.7,1x)/),9(2x,f12.7,1x))', &
         'Stiffness C (load) / GPa =', transpose(temp99_Real)*1.0e-9_pReal
       flush(IO_STDOUT)
     endif
@@ -711,7 +711,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
     errmatinv = errmatinv .or. any(dNeq(sTimesC,math_eye(size_reduced),1.0e-12_pReal))
     if (debugGeneral .or. errmatinv) then
       write(formatString, '(i2)') size_reduced
-      formatString = '(/,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))'
+      formatString = '(/,1x,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))'
       print trim(formatString), 'C * S (load) ', transpose(matmul(c_reduced,s_reduced))
       print trim(formatString), 'S (load) ', transpose(s_reduced)
       if (errmatinv) error stop 'matrix inversion error'
@@ -724,7 +724,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
   utilities_maskedCompliance = math_99to3333(temp99_Real)
 
   if (debugGeneral) then
-    print'(/,a,/,9(9(2x,f10.5,1x)/),9(2x,f10.5,1x))', &
+    print'(/,1x,a,/,9(9(2x,f10.5,1x)/),9(2x,f10.5,1x))', &
       'Masked Compliance (load) * GPa =', transpose(temp99_Real)*1.0e9_pReal
     flush(IO_STDOUT)
   endif
@@ -810,7 +810,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
   real(pReal)                     :: dPdF_norm_max, dPdF_norm_min
   real(pReal), dimension(2) :: valueAndRank                                                         !< pair of min/max norm of dPdF to synchronize min/max of dPdF
 
-  print'(/,a)', ' ... evaluating constitutive response ......................................'
+  print'(/,1x,a)', '... evaluating constitutive response ......................................'
   flush(IO_STDOUT)
 
   homogenization_F  = reshape(F,[3,3,product(grid(1:2))*grid3])                                     ! set materialpoint target F to estimated field
@@ -824,10 +824,10 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
   P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3])
   P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt
   call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr)
-  if (debugRotation) print'(/,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &
+  if (debugRotation) print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &
     'Piola--Kirchhoff stress (lab) / MPa =', transpose(P_av)*1.e-6_pReal
   if (present(rotation_BC)) P_av = rotation_BC%rotate(P_av)
-  print'(/,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &
+  print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &
     'Piola--Kirchhoff stress       / MPa =', transpose(P_av)*1.e-6_pReal
   flush(IO_STDOUT)
 
@@ -1094,7 +1094,7 @@ subroutine utilities_saveReferenceStiffness
     fileUnit,ierr
 
   if (worldrank == 0) then
-    print'(a)', ' writing reference stiffness data required for restart to file'; flush(IO_STDOUT)
+    print'(/,1x,a)', '... writing reference stiffness data required for restart to file .........'; flush(IO_STDOUT)
     open(newunit=fileUnit, file=getSolverJobName()//'.C_ref',&
          status='replace',access='stream',action='write',iostat=ierr)
     if (ierr /=0) call IO_error(100,ext_msg='could not open file '//getSolverJobName()//'.C_ref')

src/homogenization.f90

@@ -199,7 +199,7 @@ subroutine homogenization_init()
     num_homog, &
     num_homogGeneric
 
-  print'(/,a)', ' <<<+-  homogenization init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  homogenization init  -+>>>'; flush(IO_STDOUT)
 
 
   allocate(homogState      (size(material_name_homogenization)))

src/homogenization_damage.f90

@@ -41,7 +41,7 @@ module subroutine damage_init()
   integer :: ho,Nmembers
 
 
-  print'(/,a)', ' <<<+-  homogenization:damage init  -+>>>'
+  print'(/,1x,a)', '<<<+-  homogenization:damage init  -+>>>'
 
 
   configHomogenizations => config_material%get('homogenization')

src/homogenization_damage_pass.f90

@@ -8,7 +8,7 @@ contains
 
 module subroutine pass_init()
   
-  print'(/,a)', ' <<<+-  homogenization:damage:pass init  -+>>>'
+  print'(/,1x,a)', '<<<+-  homogenization:damage:pass init  -+>>>'
 
 end subroutine pass_init
 

src/homogenization_mechanical.f90

@@ -68,7 +68,7 @@ module subroutine mechanical_init(num_homog)
   class(tNode), pointer :: &
     num_homogMech
 
-  print'(/,a)', ' <<<+-  homogenization:mechanical init  -+>>>'
+  print'(/,1x,a)', '<<<+-  homogenization:mechanical init  -+>>>'
 
   call material_parseHomogenization2()
 

src/homogenization_mechanical_RGC.f90

@@ -87,16 +87,16 @@ module subroutine RGC_init(num_homogMech)
     homog, &
     homogMech
 
-  print'(/,a)', ' <<<+-  homogenization:mechanical:RGC init  -+>>>'
+  print'(/,1x,a)', '<<<+-  homogenization:mechanical:RGC init  -+>>>'
 
-  print'(a,i0)', ' # homogenizations: ',count(homogenization_type == HOMOGENIZATION_RGC_ID)
+  print'(/,a,i0)', ' # homogenizations: ',count(homogenization_type == HOMOGENIZATION_RGC_ID)
   flush(IO_STDOUT)
 
-  print*, 'D.D. Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009'
+  print'(/,1x,a)', 'D.D. Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009'
-  print*, 'https://doi.org/10.1007/s12289-009-0619-1'//IO_EOL
+  print'(  1x,a)', 'https://doi.org/10.1007/s12289-009-0619-1'//IO_EOL
 
-  print*, 'D.D. Tjahjanto et al., Modelling and Simulation in Materials Science and Engineering 18:015006, 2010'
+  print'(/,1x,a)', 'D.D. Tjahjanto et al., Modelling and Simulation in Materials Science and Engineering 18:015006, 2010'
-  print*, 'https://doi.org/10.1088/0965-0393/18/1/015006'//IO_EOL
+  print'(  1x,a)', 'https://doi.org/10.1088/0965-0393/18/1/015006'//IO_EOL
 
 
   material_homogenization => config_material%get('homogenization')
@@ -643,16 +643,16 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
   !-------------------------------------------------------------------------------------------------
   !> @brief compute the equivalent shear and bulk moduli from the elasticity tensor
   !-------------------------------------------------------------------------------------------------
-  real(pReal) function equivalentMu(grainID,ce)
+  real(pReal) function equivalentMu(co,ce)
 
     integer, intent(in)    :: &
-      grainID,&
+      co,&
       ce
 
     real(pReal), dimension(6,6) :: C
 
 
-    C = phase_homogenizedC(material_phaseID(grainID,ce),material_phaseEntry(grainID,ce))
+    C = phase_homogenizedC66(material_phaseID(co,ce),material_phaseEntry(co,ce))                    ! damage not included!
     equivalentMu = lattice_equivalent_mu(C,'voigt')
 
   end function equivalentMu

src/homogenization_mechanical_isostrain.f90

@@ -17,9 +17,9 @@ module subroutine isostrain_init
     ho, &
     Nmembers
 
-  print'(/,a)', ' <<<+-  homogenization:mechanical:isostrain init  -+>>>'
+  print'(/,1x,a)', '<<<+-  homogenization:mechanical:isostrain init  -+>>>'
 
-  print'(a,i0)', ' # homogenizations: ',count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)
+  print'(/,a,i0)', ' # homogenizations: ',count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)
   flush(IO_STDOUT)
 
   do ho = 1, size(homogenization_type)

src/homogenization_mechanical_pass.f90

@@ -17,9 +17,9 @@ module subroutine pass_init
     ho, &
     Nmembers
 
-  print'(/,a)', ' <<<+-  homogenization:mechanical:pass init  -+>>>'
+  print'(/,1x,a)', '<<<+-  homogenization:mechanical:pass init  -+>>>'
 
-  print'(a,i0)', ' # homogenizations: ',count(homogenization_type == HOMOGENIZATION_NONE_ID)
+  print'(/,a,i0)', ' # homogenizations: ',count(homogenization_type == HOMOGENIZATION_NONE_ID)
   flush(IO_STDOUT)
 
   do ho = 1, size(homogenization_type)

src/homogenization_thermal.f90

@@ -44,7 +44,7 @@ module subroutine thermal_init()
   integer :: ho
 
 
-  print'(/,a)', ' <<<+-  homogenization:thermal init  -+>>>'
+  print'(/,1x,a)', '<<<+-  homogenization:thermal init  -+>>>'
 
 
   configHomogenizations => config_material%get('homogenization')

src/homogenization_thermal_pass.f90

@@ -8,7 +8,7 @@ contains
 
 module subroutine pass_init()
   
-  print'(/,a)', ' <<<+-  homogenization:thermal:pass init  -+>>>'
+  print'(/,1x,a)', '<<<+-  homogenization:thermal:pass init  -+>>>'
 
 end subroutine pass_init
 

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, &
@@ -405,11 +405,11 @@ module lattice
 contains
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Module initialization
+!> @brief module initialization
 !--------------------------------------------------------------------------------------------------
 subroutine lattice_init
 
-  print'(/,a)', ' <<<+-  lattice init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  lattice init  -+>>>'; flush(IO_STDOUT)
 
   call selfTest
 
@@ -417,7 +417,7 @@ end subroutine lattice_init
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Characteristic shear for twinning
+!> @brief characteristic shear for twinning
 !--------------------------------------------------------------------------------------------------
 function lattice_characteristicShear_Twin(Ntwin,lattice,CoverA) result(characteristicShear)
 
@@ -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)
 
@@ -505,6 +505,7 @@ function lattice_C66_twin(Ntwin,C66,lattice,CoverA)
   type(rotation)                        :: R
   integer                               :: i
 
+
   select case(lattice)
     case('cF')
       coordinateSystem = buildCoordinateSystem(Ntwin,FCC_NSLIPSYSTEM,FCC_SYSTEMTWIN,&
@@ -521,14 +522,14 @@ function lattice_C66_twin(Ntwin,C66,lattice,CoverA)
 
   do i = 1, sum(Ntwin)
     call R%fromAxisAngle([coordinateSystem(1:3,2,i),PI],P=1)                                        ! ToDo: Why always 180 deg?
-    lattice_C66_twin(1:6,1:6,i) = R%rotTensor4sym(C66)
+    lattice_C66_twin(1:6,1:6,i) = math_3333toVoigt66(R%rotTensor4(math_Voigt66to3333(C66)))
   enddo
 
 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)
@@ -580,14 +581,14 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
 
   do i = 1,sum(Ntrans)
     call R%fromMatrix(Q(1:3,1:3,i))
-    lattice_C66_trans(1:6,1:6,i) = R%rotTensor4sym(C_target_unrotated66)
+    lattice_C66_trans(1:6,1:6,i) = math_3333toVoigt66(R%rotTensor4(math_Voigt66to3333(C_target_unrotated66)))
   enddo
 
  end function lattice_C66_trans
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Non-schmid projections for bcc with up to 6 coefficients
+!> @brief non-Schmid projections for bcc with up to 6 coefficients
 ! Koester et al. 2012, Acta Materialia 60 (2012) 3894–3901, eq. (17)
 ! Gröger et al. 2008, Acta Materialia 56 (2008) 5412–5425, table 1
 !--------------------------------------------------------------------------------------------------
@@ -634,7 +635,7 @@ end function lattice_nonSchmidMatrix
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Slip-slip interaction matrix
+!> @brief slip-slip interaction matrix
 !> details only active slip systems are considered
 !--------------------------------------------------------------------------------------------------
 function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(interactionMatrix)
@@ -750,22 +751,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, &
 
@@ -775,7 +777,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, &
@@ -785,7 +787,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)
@@ -882,7 +884,7 @@ end function lattice_interaction_SlipBySlip
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Twin-twin interaction matrix
+!> @brief twin-twin interaction matrix
 !> details only active twin systems are considered
 !--------------------------------------------------------------------------------------------------
 function lattice_interaction_TwinByTwin(Ntwin,interactionValues,lattice) result(interactionMatrix)
@@ -931,31 +933,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
@@ -980,7 +983,7 @@ end function lattice_interaction_TwinByTwin
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Trans-trans interaction matrix
+!> @brief trans-trans interaction matrix
 !> details only active trans systems are considered
 !--------------------------------------------------------------------------------------------------
 function lattice_interaction_TransByTrans(Ntrans,interactionValues,lattice) result(interactionMatrix)
@@ -1022,7 +1025,7 @@ end function lattice_interaction_TransByTrans
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Slip-twin interaction matrix
+!> @brief slip-twin interaction matrix
 !> details only active slip and twin systems are considered
 !--------------------------------------------------------------------------------------------------
 function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) result(interactionMatrix)
@@ -1122,22 +1125,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, &
 
@@ -1147,7 +1151,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, &
@@ -1157,7 +1161,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
@@ -1185,7 +1189,7 @@ end function lattice_interaction_SlipByTwin
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Slip-trans interaction matrix
+!> @brief slip-trans interaction matrix
 !> details only active slip and trans systems are considered
 !--------------------------------------------------------------------------------------------------
 function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,lattice) result(interactionMatrix)
@@ -1238,7 +1242,7 @@ function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,lattice)
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Twin-slip interaction matrix
+!> @brief twin-slip interaction matrix
 !> details only active twin and slip systems are considered
 !--------------------------------------------------------------------------------------------------
 function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,lattice) result(interactionMatrix)
@@ -1261,31 +1265,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
@@ -1410,7 +1415,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)
@@ -1482,7 +1487,7 @@ end function lattice_SchmidMatrix_cleavage
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Slip direction of slip systems (|| b)
+!> @brief slip direction of slip systems (|| b)
 !--------------------------------------------------------------------------------------------------
 function lattice_slip_direction(Nslip,lattice,cOverA) result(d)
 
@@ -1500,7 +1505,7 @@ end function lattice_slip_direction
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Normal direction of slip systems (|| n)
+!> @brief normal direction of slip systems (|| n)
 !--------------------------------------------------------------------------------------------------
 function lattice_slip_normal(Nslip,lattice,cOverA) result(n)
 
@@ -1518,7 +1523,7 @@ end function lattice_slip_normal
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Transverse direction of slip systems (|| t = b x n)
+!> @brief transverse direction of slip systems (|| t = b x n)
 !--------------------------------------------------------------------------------------------------
 function lattice_slip_transverse(Nslip,lattice,cOverA) result(t)
 
@@ -1536,7 +1541,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)
@@ -1577,7 +1582,7 @@ end function lattice_labels_slip
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Return 3x3 tensor with symmetry according to given Bravais lattice
+!> @brief return 3x3 tensor with symmetry according to given Bravais lattice
 !--------------------------------------------------------------------------------------------------
 pure function lattice_symmetrize_33(T,lattice) result(T_sym)
 
@@ -1604,7 +1609,7 @@ end function lattice_symmetrize_33
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Return stiffness matrix in 6x6 notation with symmetry according to given Bravais lattice
+!> @brief return stiffness matrix in 6x6 notation with symmetry according to given Bravais lattice
 !> @details J. A. Rayne and B. S. Chandrasekhar Phys. Rev. 120, 1658 Erratum Phys. Rev. 122, 1962
 !--------------------------------------------------------------------------------------------------
 pure function lattice_symmetrize_C66(C66,lattice) result(C66_sym)
@@ -1650,7 +1655,7 @@ end function lattice_symmetrize_C66
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Labels for twin systems
+!> @brief labels of twin systems
 !> details only active twin systems are considered
 !--------------------------------------------------------------------------------------------------
 function lattice_labels_twin(Ntwin,lattice) result(labels)
@@ -1688,7 +1693,7 @@ end function lattice_labels_twin
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Projection of the transverse direction onto the slip plane
+!> @brief projection of the transverse direction onto the slip plane
 !> @details: This projection is used to calculate forest hardening for edge dislocations
 !--------------------------------------------------------------------------------------------------
 function slipProjection_transverse(Nslip,lattice,cOverA) result(projection)
@@ -1712,7 +1717,7 @@ end function slipProjection_transverse
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Projection of the slip direction onto the slip plane
+!> @brief projection of the slip direction onto the slip plane
 !> @details: This projection is used to calculate forest hardening for screw dislocations
 !--------------------------------------------------------------------------------------------------
 function slipProjection_direction(Nslip,lattice,cOverA) result(projection)
@@ -1778,7 +1783,7 @@ end function coordinateSystem_slip
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Populate reduced interaction matrix
+!> @brief populate reduced interaction matrix
 !--------------------------------------------------------------------------------------------------
 function buildInteraction(reacting_used,acting_used,reacting_max,acting_max,values,matrix)
 
@@ -1821,7 +1826,7 @@ end function buildInteraction
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Build a local coordinate system on slip, twin, trans, cleavage systems
+!> @brief build a local coordinate system on slip, twin, trans, cleavage systems
 !> @details Order: Direction, plane (normal), and common perpendicular
 !--------------------------------------------------------------------------------------------------
 function buildCoordinateSystem(active,potential,system,lattice,cOverA)
@@ -1888,7 +1893,7 @@ end function buildCoordinateSystem
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Helper function to define transformation systems
+!> @brief helper function to define transformation systems
 ! Needed to calculate Schmid matrix and rotated stiffness matrices.
 ! @details: set c/a   = 0.0 for fcc -> bcc transformation
 !           set a_Xcc = 0.0 for fcc -> hex transformation
@@ -2072,7 +2077,7 @@ end function getlabels
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Equivalent Poisson's ratio (ν)
+!> @brief equivalent Poisson's ratio (ν)
 !> @details https://doi.org/10.1143/JPSJ.20.635
 !--------------------------------------------------------------------------------------------------
 function lattice_equivalent_nu(C,assumption) result(nu)
@@ -2105,7 +2110,7 @@ end function lattice_equivalent_nu
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Equivalent shear modulus (μ)
+!> @brief equivalent shear modulus (μ)
 !> @details https://doi.org/10.1143/JPSJ.20.635
 !--------------------------------------------------------------------------------------------------
 function lattice_equivalent_mu(C,assumption) result(mu)
@@ -2134,7 +2139,7 @@ end function lattice_equivalent_mu
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Check correctness of some lattice functions.
+!> @brief check correctness of some lattice functions
 !--------------------------------------------------------------------------------------------------
 subroutine selfTest
 

src/material.f90

@@ -26,6 +26,7 @@ module material
 
 
   type(tRotationContainer), dimension(:), allocatable :: material_O_0
+  type(tTensorContainer),   dimension(:), allocatable :: material_F_i_0
 
   integer, dimension(:), allocatable, public, protected :: &
     homogenization_Nconstituents                                                                    !< number of grains in each homogenization
@@ -48,6 +49,7 @@ module material
   public :: &
     tTensorContainer, &
     tRotationContainer, &
+    material_F_i_0, &
     material_O_0, &
     material_init
 
@@ -61,11 +63,11 @@ subroutine material_init(restart)
   logical, intent(in) :: restart
 
 
-  print'(/,a)', ' <<<+-  material init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  material init  -+>>>'; flush(IO_STDOUT)
 
 
   call parse
-  print*, 'parsed material.yaml'
+  print'(/,1x,a)', 'parsed material.yaml'
 
 
   if (.not. restart) then
@@ -159,14 +161,17 @@ subroutine parse()
   end do
 
   allocate(material_O_0(materials%length))
+  allocate(material_F_i_0(materials%length))
 
   do ma = 1, materials%length
     material     => materials%get(ma)
     constituents => material%get('constituents')
     allocate(material_O_0(ma)%data(constituents%length))
+    allocate(material_F_i_0(ma)%data(1:3,1:3,constituents%length))
     do co = 1, constituents%length
       constituent => constituents%get(co)
       call material_O_0(ma)%data(co)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4))
+      material_F_i_0(ma)%data(1:3,1:3,co) = constituent%get_as2dFloat('F_i',defaultVal=math_I3,requiredShape=[3,3])
     enddo
  enddo
 

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 :: &
@@ -91,7 +91,7 @@ subroutine math_init
   class(tNode), pointer :: &
     num_generic
 
-  print'(/,a)', ' <<<+-  math init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  math init  -+>>>'; flush(IO_STDOUT)
 
   num_generic => config_numerics%get('generic',defaultVal=emptyDict)
   randomSeed  = num_generic%get_asInt('random_seed', defaultVal = 0)
@@ -109,7 +109,7 @@ subroutine math_init
   call random_seed(put = randInit)
   call random_number(randTest)
 
-  print'(a,i2)',                ' size  of random seed:     ', randSize
+  print'(/,a,i2)',                ' size  of random seed:     ', randSize
   print'(  a,i0)',                ' value of random seed:     ', randInit(1)
   print'(  a,4(/,26x,f17.14),/)', ' start of random sequence: ', randTest
 
@@ -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
@@ -854,12 +854,13 @@ end function math_66toSym3333
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief convert 66 Voigt matrix into symmetric 3x3x3x3 matrix
+!> @brief convert 6x6 Voigt matrix into symmetric 3x3x3x3 matrix
 !--------------------------------------------------------------------------------------------------
 pure function math_Voigt66to3333(m66)
 
   real(pReal), dimension(3,3,3,3) :: math_Voigt66to3333
   real(pReal), dimension(6,6), intent(in) :: m66                                                    !< 6x6 matrix
+
   integer :: i,j
 
 
@@ -873,6 +874,31 @@ pure function math_Voigt66to3333(m66)
 end function math_Voigt66to3333
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief convert symmetric 3x3x3x3 matrix into 6x6 Voigt matrix
+!--------------------------------------------------------------------------------------------------
+pure function math_3333toVoigt66(m3333)
+
+  real(pReal), dimension(6,6)                 :: math_3333toVoigt66
+  real(pReal), dimension(3,3,3,3), intent(in) :: m3333                                              !< symmetric 3x3x3x3 matrix (no internal check)
+
+  integer :: i,j
+
+
+#ifndef __INTEL_COMPILER
+  do concurrent(i=1:6, j=1:6)
+    math_3333toVoigt66(i,j) = m3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j))
+  end do
+#else
+  do i=1,6; do j=1,6
+    math_3333toVoigt66(i,j) = m3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j))
+  end do; end do
+#endif
+
+end function math_3333toVoigt66
+
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief draw a random sample from Gauss variable
 !--------------------------------------------------------------------------------------------------
@@ -1261,6 +1287,9 @@ subroutine selfTest
   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))) &
+    error stop 'math_3333toVoigt66/math_Voigt66to3333'
+
   call random_number(v6)
   if (any(dNeq0(math_6toSym33(v6) - math_symmetric33(math_6toSym33(v6))))) &
     error stop 'math_symmetric33'

src/mesh/DAMASK_mesh.f90

@@ -86,7 +86,7 @@ program DAMASK_mesh
 !--------------------------------------------------------------------------------------------------
 ! init DAMASK (all modules)
   call CPFEM_initAll
-  print'(/,a)', ' <<<+-  DAMASK_mesh init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  DAMASK_mesh init  -+>>>'; flush(IO_STDOUT)
 
 !---------------------------------------------------------------------
 ! reading field information from numerics file and do sanity checks
@@ -204,25 +204,25 @@ program DAMASK_mesh
   errorID = 0
   checkLoadcases: do currentLoadCase = 1, size(loadCases)
     write (loadcase_string, '(i0)' ) currentLoadCase
-    print'(a,i0)', ' load case: ', currentLoadCase
+    print'(/,1x,a,i0)', 'load case: ', currentLoadCase
     if (.not. loadCases(currentLoadCase)%followFormerTrajectory) &
-      print'(a)', '  drop guessing along trajectory'
+      print'(2x,a)', 'drop guessing along trajectory'
-    print'(a)', '  Field '//trim(FIELD_MECH_label)
+    print'(2x,a)', 'Field '//trim(FIELD_MECH_label)
 
     do faceSet = 1, mesh_Nboundaries
        do component = 1, loadCases(currentLoadCase)%fieldBC(1)%nComponents
          if (loadCases(currentLoadCase)%fieldBC(1)%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, &
-                                       ' Value ', loadCases(currentLoadCase)%fieldBC(1)% &
+           ' Value ', loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Value(faceSet)
-                                                          componentBC(component)%Value(faceSet)
       end do
     end do
-    print'(a,f12.6)', '  time:       ', loadCases(currentLoadCase)%time
+    print'(2x,a,f12.6)', 'time:       ', loadCases(currentLoadCase)%time
     if (loadCases(currentLoadCase)%incs < 1)             errorID = 835                            ! non-positive incs count
-    print'(a,i5)',    '  increments: ', loadCases(currentLoadCase)%incs
+    print'(2x,a,i5)',    'increments: ', loadCases(currentLoadCase)%incs
     if (loadCases(currentLoadCase)%outputfrequency < 1)  errorID = 836                            ! non-positive result frequency
-    print'(a,i5)',    '  output  frequency:  ', &
+    print'(2x,a,i5)',    'output frequency: ', &
                loadCases(currentLoadCase)%outputfrequency
     if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string)                 ! exit with error message
   end do checkLoadcases
@@ -237,7 +237,7 @@ program DAMASK_mesh
     write(statUnit,'(a)') 'Increment Time CutbackLevel Converged IterationsNeeded'                  ! statistics file
   end if
 
-  print'(/,a)', ' ... writing initial configuration to file ........................'
+  print'(/,1x,a)', '... writing initial configuration to file .................................'
   flush(IO_STDOUT)
   call CPFEM_results(0,0.0_pReal)
 
@@ -262,7 +262,7 @@ program DAMASK_mesh
 
 !--------------------------------------------------------------------------------------------------
 ! report begin of new step
-        print'(/,a)', ' ###########################################################################'
+        print'(/,1x,a)', '###########################################################################'
         print'(1x,a,es12.5,6(a,i0))',&
                 'Time', time, &
                 's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,&
@@ -298,7 +298,7 @@ program DAMASK_mesh
             cutBackLevel = cutBackLevel + 1
             time    = time - timeinc                                                                ! rewind time
             timeinc = timeinc/2.0_pReal
-            print'(/,a)', ' cutting back'
+            print'(/,1x,a)', 'cutting back'
           else                                                                                      ! default behavior, exit if spectral solver does not converge
             if (worldrank == 0) close(statUnit)
             call IO_error(950)
@@ -315,13 +315,13 @@ program DAMASK_mesh
       cutBackLevel = max(0, cutBackLevel - 1)                                                       ! try half number of subincs next inc
 
       if (all(solres(:)%converged)) then
-        print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' converged'
+        print'(/,1x,a,i0,a)', 'increment ', totalIncsCounter, ' converged'
       else
-        print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' NOT converged'
+        print'(/,1x,a,i0,a)', 'increment ', totalIncsCounter, ' NOT converged'
       end if; flush(IO_STDOUT)
 
       if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then                            ! at output frequency
-        print'(/,a)', ' ... writing results to file ......................................'
+        print'(/,1x,a)', '... writing results to file ...............................................'
         call FEM_mechanical_updateCoords
         call CPFEM_results(totalIncsCounter,time)
       end if
@@ -334,7 +334,7 @@ program DAMASK_mesh
 
 !--------------------------------------------------------------------------------------------------
 ! report summary of whole calculation
-  print'(/,a)', ' ###########################################################################'
+  print'(/,1x,a)', '###########################################################################'
   if (worldrank == 0) close(statUnit)
 
   call quit(0)                                                                                      ! no complains ;)

src/mesh/FEM_quadrature.f90

@@ -41,10 +41,10 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine FEM_quadrature_init()
 
-  print'(/,a)', ' <<<+-  FEM_quadrature init  -+>>>'; flush(6)
+  print'(/,1x,a)', '<<<+-  FEM_quadrature init  -+>>>'; flush(6)
 
-  print*, 'L. Zhang et al., Journal of Computational Mathematics 27(1):89-96, 2009'
+  print'(/,1x,a)', 'L. Zhang et al., Journal of Computational Mathematics 27(1):89-96, 2009'
-  print*, 'https://www.jstor.org/stable/43693493'
+  print'(  1x,a)', 'https://www.jstor.org/stable/43693493'
 
 !--------------------------------------------------------------------------------------------------
 ! 2D linear

src/mesh/FEM_utilities.f90

@@ -96,7 +96,7 @@ subroutine FEM_utilities_init
   logical :: debugPETSc                                                                             !< use some in debug defined options for more verbose PETSc solution
 
 
-  print'(/,a)',   ' <<<+-  FEM_utilities init  -+>>>'
+  print'(/,1x,a)',   '<<<+-  FEM_utilities init  -+>>>'
 
   num_mesh    => config_numerics%get('mesh',defaultVal=emptyDict)
 
@@ -111,7 +111,7 @@ subroutine FEM_utilities_init
   debug_mesh  => config_debug%get('mesh',defaultVal=emptyList)
   debugPETSc  =  debug_mesh%contains('PETSc')
 
-  if(debugPETSc) print'(3(/,a),/)', &
+  if(debugPETSc) print'(3(/,1x,a),/)', &
                  'Initializing PETSc with debug options: ', &
                  trim(PETScDebug), &
                  'add more using the "PETSc_options" keyword in numerics.yaml'
@@ -149,8 +149,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
 
   PetscErrorCode :: ierr
 
-
+  print'(/,1x,a)', '... evaluating constitutive response ......................................'
-  print'(/,a)', ' ... evaluating constitutive response ......................................'
 
   call homogenization_mechanical_response(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems])                ! calculate P field
   if (.not. terminallyIll) &

src/mesh/discretization_mesh.f90

@@ -88,7 +88,7 @@ subroutine discretization_mesh_init(restart)
   integer :: p_i                                                                                    !< integration order (quadrature rule)
   type(tvec) :: coords_node0
 
-  print'(/,a)',   ' <<<+-  discretization_mesh init  -+>>>'
+  print'(/,1x,a)',   '<<<+-  discretization_mesh init  -+>>>'
 
 !--------------------------------------------------------------------------------
 ! read numerics parameter
@@ -106,6 +106,7 @@ subroutine discretization_mesh_init(restart)
   CHKERRQ(ierr)
   call DMGetStratumSize(globalMesh,'depth',dimPlex,mesh_NcpElemsGlobal,ierr)
   CHKERRQ(ierr)
+  print'()'
   call DMView(globalMesh, PETSC_VIEWER_STDOUT_WORLD,ierr)
   CHKERRQ(ierr)
 

src/mesh/mesh_mech_FEM.f90

@@ -113,7 +113,7 @@ subroutine FEM_mechanical_init(fieldBC)
   class(tNode), pointer :: &
     num_mesh
 
-  print'(/,a)', ' <<<+-  FEM_mech init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  FEM_mech init  -+>>>'; flush(IO_STDOUT)
 
 !-----------------------------------------------------------------------------
 ! read numerical parametes and do sanity checks
@@ -302,7 +302,7 @@ type(tSolutionState) function FEM_mechanical_solution( &
     CHKERRQ(ierr)
   endif
 
-  print'(/,a)', ' ==========================================================================='
+  print'(/,1x,a)', '==========================================================================='
   flush(IO_STDOUT)
 
 end function FEM_mechanical_solution
@@ -663,7 +663,7 @@ subroutine FEM_mechanical_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reaso
   print'(/,1x,a,a,i0,a,i0,f0.3)', trim(incInfo), &
                   ' @ Iteration ',PETScIter,' mechanical residual norm = ', &
                                                   int(fnorm/divTol),fnorm/divTol-int(fnorm/divTol)
-  print'(/,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &
+  print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &
     'Piola--Kirchhoff stress / MPa =',transpose(P_av)*1.e-6_pReal
   flush(IO_STDOUT)
 

src/parallelization.f90

@@ -76,11 +76,11 @@ subroutine parallelization_init
   call MPI_Comm_rank(MPI_COMM_WORLD,worldrank,err)
   if (err /= 0)                              error stop 'Could not determine worldrank'
 
-  if (worldrank == 0) print'(/,a)',  ' <<<+-  parallelization init  -+>>>'
+  if (worldrank == 0) print'(/,1x,a)', '<<<+-  parallelization init  -+>>>'
 
   call MPI_Comm_size(MPI_COMM_WORLD,worldsize,err)
   if (err /= 0)                              error stop 'Could not determine worldsize'
-  if (worldrank == 0) print'(a,i3)', ' MPI processes: ',worldsize
+  if (worldrank == 0) print'(/,1x,a,i3)', 'MPI processes: ',worldsize
 
   call MPI_Type_size(MPI_INTEGER,typeSize,err)
   if (err /= 0)                              error stop 'Could not determine MPI integer size'
@@ -97,16 +97,16 @@ subroutine parallelization_init
 
 !$ call get_environment_variable(name='OMP_NUM_THREADS',value=NumThreadsString,STATUS=got_env)
 !$ if(got_env /= 0) then
-!$   print*, 'Could not get $OMP_NUM_THREADS, using default'
+!$   print'(1x,a)', 'Could not get $OMP_NUM_THREADS, using default'
 !$   OMP_NUM_THREADS = 4_pI32
 !$ else
 !$   read(NumThreadsString,'(i6)') OMP_NUM_THREADS
 !$   if (OMP_NUM_THREADS < 1_pI32) then
-!$     print*, 'Invalid OMP_NUM_THREADS: "'//trim(NumThreadsString)//'", using default'
+!$     print'(1x,a)', 'Invalid OMP_NUM_THREADS: "'//trim(NumThreadsString)//'", using default'
 !$     OMP_NUM_THREADS = 4_pI32
 !$   endif
 !$ endif
-!$ print'(a,i2)',   ' OMP_NUM_THREADS: ',OMP_NUM_THREADS
+!$ print'(1x,a,1x,i2)',   'OMP_NUM_THREADS:',OMP_NUM_THREADS
 !$ call omp_set_num_threads(OMP_NUM_THREADS)
 
 end subroutine parallelization_init

src/phase.f90

@@ -5,6 +5,7 @@
 !--------------------------------------------------------------------------------------------------
 module phase
   use prec
+  use constants
   use math
   use rotations
   use IO
@@ -230,15 +231,15 @@ module phase
     end function phase_mechanical_constitutive
 
     !ToDo: Merge all the stiffness functions
-    module function phase_homogenizedC(ph,en) result(C)
+    module function phase_homogenizedC66(ph,en) result(C)
       integer, intent(in) :: ph, en
       real(pReal), dimension(6,6) :: C
-    end function phase_homogenizedC
+    end function phase_homogenizedC66
-    module function phase_damage_C(C_homogenized,ph,en) result(C)
+    module function phase_damage_C66(C66,ph,en) result(C66_degraded)
-      real(pReal), dimension(3,3,3,3), intent(in)  :: C_homogenized
+      real(pReal), dimension(6,6), intent(in)  :: C66
       integer,                     intent(in)  :: ph,en
-      real(pReal), dimension(3,3,3,3) :: C
+      real(pReal), dimension(6,6) :: C66_degraded
-    end function phase_damage_C
+    end function phase_damage_C66
 
     module function phase_f_phi(phi,co,ce) result(f)
       integer, intent(in) :: ce,co
@@ -299,7 +300,7 @@ module phase
 
   public :: &
     phase_init, &
-    phase_homogenizedC, &
+    phase_homogenizedC66, &
     phase_f_phi, &
     phase_f_T, &
     phase_K_phi, &
@@ -343,7 +344,7 @@ subroutine phase_init
     phase
 
 
-  print'(/,a)', ' <<<+-  phase init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  phase init  -+>>>'; flush(IO_STDOUT)
 
   debug_constitutive => config_debug%get('phase', defaultVal=emptyList)
   debugConstitutive%basic     = debug_constitutive%contains('basic')
@@ -495,7 +496,7 @@ subroutine crystallite_init()
     phases
 
 
-  print'(/,a)', ' <<<+-  crystallite init  -+>>>'
+  print'(/,1x,a)', '<<<+-  crystallite init  -+>>>'
 
   cMax = homogenization_maxNconstituents
   iMax = discretization_nIPs
@@ -534,7 +535,7 @@ subroutine crystallite_init()
 
   phases => config_material%get('phase')
 
-  print'(a42,1x,i10)', '    # of elements:                       ', eMax
+  print'(/,a42,1x,i10)', '    # of elements:                       ', eMax
   print'(  a42,1x,i10)', '    # of integration points/element:     ', iMax
   print'(  a42,1x,i10)', 'max # of constituents/integration point: ', cMax
   flush(IO_STDOUT)

src/phase_damage.f90

@@ -83,7 +83,7 @@ module subroutine damage_init
    source
   logical:: damage_active
 
-  print'(/,a)', ' <<<+-  phase:damage init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:damage init  -+>>>'
 
   phases => config_material%get('phase')
 
@@ -139,6 +139,7 @@ module function phase_damage_constitutive(Delta_t,co,ip,el) result(converged_)
   integer :: &
     ph, en
 
+
   ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
   en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
 
@@ -150,20 +151,21 @@ end function phase_damage_constitutive
 !--------------------------------------------------------------------------------------------------
 !> @brief returns the degraded/modified elasticity matrix
 !--------------------------------------------------------------------------------------------------
-module function phase_damage_C(C_homogenized,ph,en) result(C)
+module function phase_damage_C66(C66,ph,en) result(C66_degraded)
 
-  real(pReal), dimension(3,3,3,3), intent(in)  :: C_homogenized
+  real(pReal), dimension(6,6), intent(in)  :: C66
   integer,                     intent(in)  :: ph,en
-  real(pReal), dimension(3,3,3,3) :: C
+  real(pReal), dimension(6,6) :: C66_degraded
+
 
   damageType: select case (phase_damage(ph))
     case (DAMAGE_ISOBRITTLE_ID) damageType
-     C = C_homogenized * damage_phi(ph,en)**2
+     C66_degraded = C66 * damage_phi(ph,en)**2
     case default damageType
-     C = C_homogenized
+     C66_degraded = C66
   end select damageType
 
-end function phase_damage_C
+end function phase_damage_C66
 
 
 !--------------------------------------------------------------------------------------------------
@@ -417,7 +419,7 @@ function phase_damage_deltaState(Fe, ph, en) result(broken)
    sourceType: select case (phase_damage(ph))
 
     case (DAMAGE_ISOBRITTLE_ID) sourceType
-      call isobrittle_deltaState(phase_homogenizedC(ph,en), Fe, ph,en)
+      call isobrittle_deltaState(phase_homogenizedC66(ph,en), Fe, ph,en)
       broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,en)))
       if (.not. broken) then
         myOffset = damageState(ph)%offsetDeltaState

src/phase_damage_anisobrittle.f90

@@ -48,8 +48,8 @@ module function anisobrittle_init() result(mySources)
   mySources = source_active('anisobrittle')
   if (count(mySources) == 0) return
 
-  print'(/,a)', ' <<<+-  phase:damage:anisobrittle init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:damage:anisobrittle init  -+>>>'
-  print'(a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT)
+  print'(/,a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT)
 
 
   phases => config_material%get('phase')

src/phase_damage_isobrittle.f90

@@ -46,8 +46,8 @@ module function isobrittle_init() result(mySources)
   mySources = source_active('isobrittle')
   if (count(mySources) == 0) return
 
-  print'(/,a)', ' <<<+-  phase:damage:isobrittle init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:damage:isobrittle init  -+>>>'
-  print'(a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT)
+  print'(/,a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT)
 
 
   phases => config_material%get('phase')
@@ -96,6 +96,7 @@ end function isobrittle_init
 
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates derived quantities from state
+! ToDo: Use Voigt directly
 !--------------------------------------------------------------------------------------------------
 module subroutine isobrittle_deltaState(C, Fe, ph,en)
 
@@ -109,13 +110,16 @@ module subroutine isobrittle_deltaState(C, Fe, ph,en)
     epsilon
   real(pReal) :: &
     r_W
+  real(pReal), dimension(6,6) :: &
+    C_sym
 
 
+  C_sym = math_sym3333to66(math_Voigt66to3333(C))
   epsilon = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3)
 
   associate(prm => param(ph), stt => state(ph), dlt => deltaState(ph))
 
-    r_W = 2.0_pReal*dot_product(epsilon,matmul(C,epsilon))/prm%W_crit
+    r_W = 2.0_pReal*dot_product(epsilon,matmul(C_sym,epsilon))/prm%W_crit
     dlt%r_W(en) = merge(r_W - stt%r_W(en), 0.0_pReal, r_W > stt%r_W(en))
 
   end associate

src/phase_mechanical.f90

@@ -168,19 +168,19 @@ submodule(phase) mechanical
       integer,     intent(in) :: ph,en
     end function plastic_dislotwin_homogenizedC
 
-    module function elastic_C66(ph) result(C66)
+    module function elastic_C66(ph,en) result(C66)
       real(pReal), dimension(6,6) :: C66
-      integer,     intent(in) :: ph
+      integer,     intent(in) :: ph, en
     end function elastic_C66
 
-    module function elastic_mu(ph) result(mu)
+    module function elastic_mu(ph,en) result(mu)
       real(pReal) :: mu
-      integer, intent(in) :: ph
+      integer, intent(in) :: ph, en
     end function elastic_mu
 
-    module function elastic_nu(ph) result(nu)
+    module function elastic_nu(ph,en) result(nu)
       real(pReal) :: nu
-      integer, intent(in) :: ph
+      integer, intent(in) :: ph, en
     end function elastic_nu
 
   end interface
@@ -205,6 +205,9 @@ module subroutine mechanical_init(phases)
     phases
 
   integer :: &
+    ce, &
+    co, &
+    ma, &
     ph, &
     en, &
     Nmembers
@@ -214,7 +217,7 @@ module subroutine mechanical_init(phases)
     phase, &
     mech
 
-  print'(/,a)', ' <<<+-  phase:mechanical init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:mechanical init  -+>>>'
 
 !-------------------------------------------------------------------------------------------------
   allocate(output_constituent(phases%length))
@@ -261,15 +264,21 @@ module subroutine mechanical_init(phases)
 #endif
   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_mechanical_Fi0(ph)%data(1:3,1:3,material_phaseEntry(co,ce)) = material_F_i_0(ma)%data(1:3,1:3,co)
+    enddo
+  enddo
+
   do ph = 1, phases%length
     do en = 1, count(material_phaseID == ph)
 
       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
       phase_mechanical_F0(ph)%data(1:3,1:3,en) = math_I3
-
       phase_mechanical_Fe(ph)%data(1:3,1:3,en) = math_inv33(matmul(phase_mechanical_Fi0(ph)%data(1:3,1:3,en), &
                                                                    phase_mechanical_Fp0(ph)%data(1:3,1:3,en)))  ! assuming that euler angles are given in internal strain free configuration
     enddo

src/phase_mechanical_eigen.f90

@@ -44,7 +44,7 @@ module subroutine eigen_init(phases)
     kinematics, &
     mechanics
 
-  print'(/,a)', ' <<<+-  phase:mechanical:eigen init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:mechanical:eigen init  -+>>>'
 
 !--------------------------------------------------------------------------------------------------
 ! explicit eigen mechanisms

src/phase_mechanical_eigen_cleavageopening.f90

@@ -21,8 +21,8 @@ module function damage_anisobrittle_init() result(myKinematics)
   myKinematics = kinematics_active2('anisobrittle')
   if(count(myKinematics) == 0) return
 
-  print'(/,a)', ' <<<+-  phase:mechanical:eigen:cleavageopening init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:mechanical:eigen:cleavageopening init  -+>>>'
-  print'(a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT)
+  print'(/,a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT)
 
 end function damage_anisobrittle_init
 

src/phase_mechanical_eigen_thermalexpansion.f90

@@ -36,11 +36,11 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics)
     kinematics, &
     kinematic_type
 
-  print'(/,a)', ' <<<+-  phase:mechanical:eigen:thermalexpansion init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:mechanical:eigen:thermalexpansion init  -+>>>'
 
   myKinematics = kinematics_active('thermalexpansion',kinematics_length)
   Ninstances = count(myKinematics)
-  print'(a,i2)', ' # phases: ',Ninstances; flush(IO_STDOUT)
+  print'(/,a,i2)', ' # phases: ',Ninstances; flush(IO_STDOUT)
   if (Ninstances == 0) return
 
   phases => config_material%get('phase')
@@ -58,7 +58,7 @@ 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)
@@ -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,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,18 +1,24 @@
 submodule(phase:mechanical) elastic
 
   type :: tParameters
-    real(pReal), dimension(6,6) :: &
+    real(pReal),dimension(3) :: &
-      C66 = 0.0_pReal                                                                               !< Elastic constants in Voigt notation
+      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) :: &
-      mu, &
+      T_ref
-      nu
   end type tParameters
 
   type(tParameters), allocatable, dimension(:) :: param
 
 contains
 
-
+!--------------------------------------------------------------------------------------------------
+!> @brief initialize elasticity
+!--------------------------------------------------------------------------------------------------
 module subroutine elastic_init(phases)
 
   class(tNode), pointer :: &
@@ -24,12 +30,12 @@ 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  -+>>>'
 
-  print'(/,a)', ' <<<+-  phase:mechanical:elastic init  -+>>>'
+  print'(/,a,i0)', ' # phases: ',phases%length; flush(IO_STDOUT)
-  print'(/,a)', ' <<<+-  phase:mechanical:elastic:Hooke init  -+>>>'
-
-  print'(a,i0)', ' # phases: ',phases%length; flush(IO_STDOUT)
 
   allocate(param(phases%length))
 
@@ -41,22 +47,35 @@ module subroutine elastic_init(phases)
 
     associate(prm => param(ph))
 
-      prm%C66(1,1) = elastic%get_asFloat('C_11')
+      prm%T_ref = elastic%get_asFloat('T_ref', defaultVal=T_ROOM)
-      prm%C66(1,2) = elastic%get_asFloat('C_12')
+
-      prm%C66(4,4) = 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%C66(1,3) = elastic%get_asFloat('C_13')
+        prm%C_13(1) = elastic%get_asFloat('C_13')
-        prm%C66(3,3) = 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') prm%C66(6,6) = elastic%get_asFloat('C_66')
 
-      prm%C66 = lattice_symmetrize_C66(prm%C66,phase_lattice(ph))
+      if (phase_lattice(ph) == 'tI') then
-
+        prm%C_66(1) = elastic%get_asFloat('C_66')
-      prm%nu = lattice_equivalent_nu(prm%C66,'voigt')
+        prm%C_66(2) = elastic%get_asFloat('C_66,T',  defaultVal=0.0_pReal)
-      prm%mu = lattice_equivalent_mu(prm%C66,'voigt')
+        prm%C_66(3) = elastic%get_asFloat('C_66,T^2',defaultVal=0.0_pReal)
-      
+      end if
-      prm%C66 = math_sym3333to66(math_Voigt66to3333(prm%C66))                                       ! Literature data is in Voigt notation
 
     end associate
   end do
@@ -65,8 +84,98 @@ end subroutine elastic_init
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
+!> @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
+    T = thermal_T(ph,en)
+
+    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(1) &
+               + 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') 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))
+
+  end associate
+
+end function elastic_C66
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief return shear modulus
+!--------------------------------------------------------------------------------------------------
+module function elastic_mu(ph,en) result(mu)
+
+  integer, intent(in) :: &
+    ph, &
+    en
+  real(pReal) :: &
+    mu
+
+
+  mu = lattice_equivalent_mu(elastic_C66(ph,en),'voigt')
+
+end function elastic_mu
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief return Poisson ratio
+!--------------------------------------------------------------------------------------------------
+module function elastic_nu(ph,en) result(nu)
+
+  integer, intent(in) :: &
+    ph, &
+    en
+  real(pReal) :: &
+    nu
+
+
+  nu = lattice_equivalent_nu(elastic_C66(ph,en),'voigt')
+
+end function elastic_nu
+
+
+
+!--------------------------------------------------------------------------------------------------
+!> @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
 !--------------------------------------------------------------------------------------------------
 module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
                                               Fe, Fi, ph, en)
@@ -89,8 +198,7 @@ module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
     i, j
 
 
-  C = math_66toSym3333(phase_homogenizedC(ph,en))
+  C = math_Voigt66to3333(phase_damage_C66(phase_homogenizedC66(ph,en),ph,en))
-  C = phase_damage_C(C,ph,en)
 
   E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3)                                                  !< Green-Lagrange strain in unloaded configuration
   S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi))                                        !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration
@@ -104,50 +212,22 @@ end subroutine phase_hooke_SandItsTangents
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief returns the homogenized elasticity matrix
+!> @brief Return the homogenized elasticity matrix.
-!> ToDo: homogenizedC66 would be more consistent
 !--------------------------------------------------------------------------------------------------
-module function phase_homogenizedC(ph,en) result(C)
+module function phase_homogenizedC66(ph,en) result(C)
 
   real(pReal), dimension(6,6) :: C
   integer,      intent(in)    :: ph, en
 
+
   plasticType: select case (phase_plasticity(ph))
     case (PLASTICITY_DISLOTWIN_ID) plasticType
      C = plastic_dislotwin_homogenizedC(ph,en)
     case default plasticType
-     C = param(ph)%C66
+     C = elastic_C66(ph,en)
   end select plasticType
 
-end function phase_homogenizedC
+end function phase_homogenizedC66
 
-module function elastic_C66(ph) result(C66)
-  real(pReal), dimension(6,6) :: C66
-  integer,     intent(in) :: ph
-
-
-  C66 = param(ph)%C66
-
-end function elastic_C66
-
-module function elastic_mu(ph) result(mu)
-
-  real(pReal) :: mu
-  integer, intent(in) :: ph
-
-
-  mu = param(ph)%mu
-
-end function elastic_mu
-
-module function elastic_nu(ph) result(nu)
-
-  real(pReal) :: nu
-  integer, intent(in) :: ph
-
-
-  nu = param(ph)%nu
-
-end function elastic_nu
 
 end submodule elastic

src/phase_mechanical_plastic.f90

@@ -222,7 +222,7 @@ contains
 module subroutine plastic_init
 
 
-  print'(/,a)', ' <<<+-  phase:mechanical:plastic init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:mechanical:plastic init  -+>>>'
 
   where(plastic_none_init())              phase_plasticity = PLASTICITY_NONE_ID
   where(plastic_isotropic_init())         phase_plasticity = PLASTICITY_ISOTROPIC_ID
@@ -262,15 +262,16 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
     i, j
 
 
+  if (phase_plasticity(ph) == PLASTICITY_NONE_ID) then
+    Lp      = 0.0_pReal
+    dLp_dFi = 0.0_pReal
+    dLp_dS  = 0.0_pReal
+  else
+
     Mp = matmul(matmul(transpose(Fi),Fi),S)
 
-
     plasticType: select case (phase_plasticity(ph))
 
-    case (PLASTICITY_NONE_ID) plasticType
-      Lp = 0.0_pReal
-      dLp_dMp = 0.0_pReal
-
       case (PLASTICITY_ISOTROPIC_ID) plasticType
         call isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
 
@@ -297,6 +298,8 @@ module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
       dLp_dS(i,j,1:3,1:3)  = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3))                     ! ToDo: @PS: why not:   dLp_dMp:(FiT Fi)
     end do; end do
 
+  end if
+
 end subroutine plastic_LpAndItsTangents
 
 
@@ -318,6 +321,7 @@ module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken)
   logical :: broken
 
 
+  if (phase_plasticity(ph) /= PLASTICITY_NONE_ID) then
     Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,en)),&
                        phase_mechanical_Fi(ph)%data(1:3,1:3,en)),phase_mechanical_S(ph)%data(1:3,1:3,en))
 
@@ -341,8 +345,9 @@ module function plastic_dotState(subdt,co,ip,el,ph,en) result(broken)
       case (PLASTICITY_NONLOCAL_ID) plasticType
         call nonlocal_dotState(Mp,thermal_T(ph,en),subdt,ph,en,ip,el)
     end select plasticType
-  broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,en)))
+  end if
 
+  broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,en)))
 
 end function plastic_dotState
 
@@ -390,8 +395,7 @@ module function plastic_deltaState(ph, en) result(broken)
   integer, intent(in) :: &
     ph, &
     en
-  logical :: &
+  logical :: broken
-    broken
 
   real(pReal),               dimension(3,3) :: &
     Mp
@@ -399,36 +403,35 @@ module function plastic_deltaState(ph, en) result(broken)
     myOffset, &
     mySize
 
+  broken = .false.
+
+  select case (phase_plasticity(ph))
+    case (PLASTICITY_NONLOCAL_ID,PLASTICITY_KINEHARDENING_ID)
 
       Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,en)),&
-                     phase_mechanical_Fi(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_S(ph)%data(1:3,1:3,en))
 
       plasticType: select case (phase_plasticity(ph))
 
         case (PLASTICITY_KINEHARDENING_ID) plasticType
           call plastic_kinehardening_deltaState(Mp,ph,en)
-      broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,en)))
 
         case (PLASTICITY_NONLOCAL_ID) plasticType
           call plastic_nonlocal_deltaState(Mp,ph,en)
-      broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,en)))
-
-    case default
-      broken = .false.
 
       end select plasticType
 
+      broken = any(IEEE_is_NaN(plasticState(ph)%deltaState(:,en)))
       if (.not. broken) then
-    select case(phase_plasticity(ph))
-      case (PLASTICITY_NONLOCAL_ID,PLASTICITY_KINEHARDENING_ID)
-
         myOffset = plasticState(ph)%offsetDeltaState
         mySize   = plasticState(ph)%sizeDeltaState
         plasticState(ph)%state(myOffset + 1:myOffset + mySize,en) = &
         plasticState(ph)%state(myOffset + 1:myOffset + mySize,en) + plasticState(ph)%deltaState(1:mySize,en)
-    end select
       end if
 
+  end select
+
 end function plastic_deltaState
 
 
@@ -453,7 +456,7 @@ function plastic_active(plastic_label)  result(active_plastic)
     phase => phases%get(ph)
     mech  => phase%get('mechanical')
     pl    => mech%get('plastic',defaultVal = emptyDict)
-    if(pl%get_asString('type',defaultVal='none') == plastic_label) active_plastic(ph) = .true.
+    active_plastic(ph) = pl%get_asString('type',defaultVal='none') == plastic_label
   enddo
 
 end function plastic_active

src/phase_mechanical_plastic_dislotungsten.f90

@@ -7,13 +7,9 @@
 !--------------------------------------------------------------------------------------------------
 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
-      mu   = 1.0_pReal, &                                                                           !< equivalent shear modulus
       D_0  = 1.0_pReal, &                                                                           !< prefactor for self-diffusion coefficient
       Q_cl = 1.0_pReal                                                                              !< activation energy for dislocation climb
     real(pReal),               allocatable, dimension(:) :: &
@@ -101,11 +97,11 @@ module function plastic_dislotungsten_init() result(myPlasticity)
   myPlasticity = plastic_active('dislotungsten')
   if (count(myPlasticity) == 0) return
 
-  print'(/,a)', ' <<<+-  phase:mechanical:plastic:dislotungsten init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:mechanical:plastic:dislotungsten init  -+>>>'
-  print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
+  print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
 
-  print*, 'D. Cereceda et al., International Journal of Plasticity 78:242–256, 2016'
+  print'(/,1x,a)', 'D. Cereceda et al., International Journal of Plasticity 78:242–256, 2016'
-  print*, 'https://doi.org/10.1016/j.ijplas.2015.09.002'
+  print'(  1x,a)', 'https://doi.org/10.1016/j.ijplas.2015.09.002'
 
 
   phases => config_material%get('phase')
@@ -130,8 +126,6 @@ module function plastic_dislotungsten_init() result(myPlasticity)
     prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray)
 #endif
 
-    prm%mu = elastic_mu(ph)
-
 !--------------------------------------------------------------------------------------------------
 ! slip related parameters
     N_sl         = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
@@ -324,10 +318,13 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
     dot_rho_dip_formation, &
     dot_rho_dip_climb, &
     d_hat
-
+  real(pReal) :: &
+    mu
 
   associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph))
 
+    mu = elastic_mu(ph,en)
+
     call kinetics(Mp,T,ph,en,&
                   dot_gamma_pos,dot_gamma_neg, &
                   tau_pos_out = tau_pos,tau_neg_out = tau_neg)
@@ -338,13 +335,13 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
       dot_rho_dip_formation = 0.0_pReal
       dot_rho_dip_climb     = 0.0_pReal
     else where
-      d_hat = math_clip(3.0_pReal*prm%mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos+tau_neg)*0.5_pReal), &
+      d_hat = math_clip(3.0_pReal*mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos+tau_neg)*0.5_pReal), &
                         prm%d_caron, &                                                              ! lower limit
                         dst%Lambda_sl(:,en))                                                        ! upper limit
       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*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
@@ -376,7 +373,7 @@ module subroutine dislotungsten_dependentState(ph,en)
 
   associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
 
-    dst%tau_pass(:,en) = prm%mu*prm%b_sl &
+    dst%tau_pass(:,en) = elastic_mu(ph,en)*prm%b_sl &
                        * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
 
     Lambda_sl_inv = 1.0_pReal/prm%D &
@@ -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,13 +9,8 @@
 !--------------------------------------------------------------------------------------------------
 submodule(phase:plastic) dislotwin
 
-  real(pReal), parameter :: &
-    kB = 1.38e-23_pReal                                                                             !< Boltzmann constant in J/Kelvin
-
   type :: tParameters
     real(pReal) :: &
-      mu                  = 1.0_pReal, &                                                            !< equivalent shear modulus
-      nu                  = 1.0_pReal, &                                                            !< equivalent shear Poisson's ratio
       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
@@ -33,7 +28,9 @@ 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) :: &
       Gamma_sf = 0.0_pReal
     real(pReal),               allocatable, dimension(:) :: &
@@ -62,20 +59,22 @@ submodule(phase:plastic) dislotwin
       h_sl_tr, &                                                                                    !< components of slip-trans interaction matrix
       h_tr_tr, &                                                                                    !< components of trans-trans interaction matrix
       n0_sl, &                                                                                      !< slip system normal
-      forestProjection, &
+      forestProjection
-      C66
     real(pReal),               allocatable, dimension(:,:,:) :: &
       P_sl, &
       P_tw, &
-      P_tr, &
+      P_tr
-      C66_tw, &
-      C66_tr
     integer :: &
       sum_N_sl, &                                                                                   !< total number of active slip system
       sum_N_tw, &                                                                                   !< total number of active twin system
       sum_N_tr                                                                                      !< total number of active transformation system
+    integer,                   allocatable, dimension(:)   :: &
+      N_tw, &
+      N_tr
     integer,                   allocatable, dimension(:,:) :: &
       fcc_twinNucleationSlipPair                                                                    ! ToDo: Better name? Is also use for trans
+    character(len=:),          allocatable                 :: &
+      lattice_tr
     character(len=pStringLen), allocatable, dimension(:) :: &
       output
     logical :: &
@@ -134,7 +133,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
     sizeState, sizeDotState, &
     startIndex, endIndex
   integer,     dimension(:), allocatable :: &
-    N_sl, N_tw, N_tr
+    N_sl
   real(pReal), allocatable, dimension(:) :: &
     rho_mob_0, &                                                                                    !< initial unipolar dislocation density per slip system
     rho_dip_0                                                                                       !< initial dipole dislocation density per slip system
@@ -150,17 +149,17 @@ module function plastic_dislotwin_init() result(myPlasticity)
   myPlasticity = plastic_active('dislotwin')
   if (count(myPlasticity) == 0) return
 
-  print'(/,a)', ' <<<+-  phase:mechanical:plastic:dislotwin init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:mechanical:plastic:dislotwin init  -+>>>'
-  print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
+  print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
 
-  print*, 'A. Ma and F. Roters, Acta Materialia 52(12):3603–3612, 2004'
+  print'(/,1x,a)', 'A. Ma and F. Roters, Acta Materialia 52(12):3603–3612, 2004'
-  print*, 'https://doi.org/10.1016/j.actamat.2004.04.012'//IO_EOL
+  print'(  1x,a)', 'https://doi.org/10.1016/j.actamat.2004.04.012'//IO_EOL
 
-  print*, 'F. Roters et al., Computational Materials Science 39:91–95, 2007'
+  print'(/,1x,a)', 'F. Roters et al., Computational Materials Science 39:91–95, 2007'
-  print*, 'https://doi.org/10.1016/j.commatsci.2006.04.014'//IO_EOL
+  print'(  1x,a)', 'https://doi.org/10.1016/j.commatsci.2006.04.014'//IO_EOL
 
-  print*, 'S.L. Wong et al., Acta Materialia 118:140–151, 2016'
+  print'(/,1x,a)', 'S.L. Wong et al., Acta Materialia 118:140–151, 2016'
-  print*, 'https://doi.org/10.1016/j.actamat.2016.07.032'
+  print'(  1x,a)', 'https://doi.org/10.1016/j.actamat.2016.07.032'
 
 
   phases => config_material%get('phase')
@@ -185,10 +184,6 @@ module function plastic_dislotwin_init() result(myPlasticity)
     prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray)
 #endif
 
-    ! This data is read in already in lattice
-    prm%mu  = elastic_mu(ph)
-    prm%nu  = elastic_nu(ph)
-    prm%C66 = elastic_C66(ph)
 
 !--------------------------------------------------------------------------------------------------
 ! slip related parameters
@@ -260,35 +255,33 @@ module function plastic_dislotwin_init() result(myPlasticity)
 
 !--------------------------------------------------------------------------------------------------
 ! twin related parameters
-    N_tw         = pl%get_as1dInt('N_tw', defaultVal=emptyIntArray)
+    prm%N_tw     = pl%get_as1dInt('N_tw', defaultVal=emptyIntArray)
-    prm%sum_N_tw = sum(abs(N_tw))
+    prm%sum_N_tw = sum(abs(prm%N_tw))
     twinActive: if (prm%sum_N_tw > 0) then
-      prm%systems_tw = lattice_labels_twin(N_tw,phase_lattice(ph))
+      prm%systems_tw = lattice_labels_twin(prm%N_tw,phase_lattice(ph))
-      prm%P_tw  = lattice_SchmidMatrix_twin(N_tw,phase_lattice(ph),phase_cOverA(ph))
+      prm%P_tw  = lattice_SchmidMatrix_twin(prm%N_tw,phase_lattice(ph),phase_cOverA(ph))
-      prm%h_tw_tw   = lattice_interaction_TwinByTwin(N_tw,pl%get_as1dFloat('h_tw-tw'), &
+      prm%h_tw_tw   = lattice_interaction_TwinByTwin(prm%N_tw,pl%get_as1dFloat('h_tw-tw'), &
                                                      phase_lattice(ph))
 
-      prm%b_tw      = pl%get_as1dFloat('b_tw',     requiredSize=size(N_tw))
+      prm%b_tw      = pl%get_as1dFloat('b_tw',     requiredSize=size(prm%N_tw))
-      prm%t_tw      = pl%get_as1dFloat('t_tw',     requiredSize=size(N_tw))
+      prm%t_tw      = pl%get_as1dFloat('t_tw',     requiredSize=size(prm%N_tw))
-      prm%r         = pl%get_as1dFloat('p_tw',     requiredSize=size(N_tw))
+      prm%r         = pl%get_as1dFloat('p_tw',     requiredSize=size(prm%N_tw))
 
       prm%x_c_tw    = pl%get_asFloat('x_c_tw')
       prm%L_tw      = pl%get_asFloat('L_tw')
       prm%i_tw      = pl%get_asFloat('i_tw')
 
-      prm%gamma_char= lattice_characteristicShear_Twin(N_tw,phase_lattice(ph),phase_cOverA(ph))
+      prm%gamma_char= lattice_characteristicShear_Twin(prm%N_tw,phase_lattice(ph),phase_cOverA(ph))
-
-      prm%C66_tw    = lattice_C66_twin(N_tw,prm%C66,phase_lattice(ph),phase_cOverA(ph))
 
       if (.not. prm%fccTwinTransNucleation) then
         prm%dot_N_0_tw = pl%get_as1dFloat('dot_N_0_tw')
-        prm%dot_N_0_tw = math_expand(prm%dot_N_0_tw,N_tw)
+        prm%dot_N_0_tw = math_expand(prm%dot_N_0_tw,prm%N_tw)
       endif
 
       ! expand: family => system
-      prm%b_tw = math_expand(prm%b_tw,N_tw)
+      prm%b_tw = math_expand(prm%b_tw,prm%N_tw)
-      prm%t_tw = math_expand(prm%t_tw,N_tw)
+      prm%t_tw = math_expand(prm%t_tw,prm%N_tw)
-      prm%r    = math_expand(prm%r,N_tw)
+      prm%r    = math_expand(prm%r,prm%N_tw)
 
       ! sanity checks
       if (    prm%x_c_tw        < 0.0_pReal)  extmsg = trim(extmsg)//' x_c_tw'
@@ -307,39 +300,38 @@ module function plastic_dislotwin_init() result(myPlasticity)
 
 !--------------------------------------------------------------------------------------------------
 ! transformation related parameters
-    N_tr         = pl%get_as1dInt('N_tr', defaultVal=emptyIntArray)
+    prm%N_tr         = pl%get_as1dInt('N_tr', defaultVal=emptyIntArray)
-    prm%sum_N_tr = sum(abs(N_tr))
+    prm%sum_N_tr = sum(abs(prm%N_tr))
     transActive: if (prm%sum_N_tr > 0) then
       prm%b_tr = pl%get_as1dFloat('b_tr')
-      prm%b_tr = math_expand(prm%b_tr,N_tr)
+      prm%b_tr = math_expand(prm%b_tr,prm%N_tr)
 
       prm%h             = pl%get_asFloat('h',       defaultVal=0.0_pReal) ! ToDo: How to handle that???
       prm%i_tr          = pl%get_asFloat('i_tr',    defaultVal=0.0_pReal) ! ToDo: How to handle that???
       prm%delta_G       = pl%get_asFloat('delta_G')
       prm%x_c_tr        = pl%get_asFloat('x_c_tr',  defaultVal=0.0_pReal) ! ToDo: How to handle that???
       prm%L_tr          = pl%get_asFloat('L_tr')
+      prm%a_cI          = pl%get_asFloat('a_cI',    defaultVal=0.0_pReal)
+      prm%a_cF          = pl%get_asFloat('a_cF',    defaultVal=0.0_pReal)
 
-      prm%h_tr_tr = lattice_interaction_TransByTrans(N_tr,pl%get_as1dFloat('h_tr-tr'),&
+      prm%lattice_tr    = pl%get_asString('lattice_tr')
+
+      prm%h_tr_tr = lattice_interaction_TransByTrans(prm%N_tr,pl%get_as1dFloat('h_tr-tr'),&
                                                      phase_lattice(ph))
 
-      prm%C66_tr  = lattice_C66_trans(N_tr,prm%C66,pl%get_asString('lattice_tr'), &
+      prm%P_tr    = lattice_SchmidMatrix_trans(prm%N_tr,prm%lattice_tr, &
                                                0.0_pReal, &
-                                      pl%get_asFloat('a_cI', defaultVal=0.0_pReal), &
+                                               prm%a_cI, &
-                                      pl%get_asFloat('a_cF', defaultVal=0.0_pReal))
+                                               prm%a_cF)
-
-      prm%P_tr    = lattice_SchmidMatrix_trans(N_tr,pl%get_asString('lattice_tr'), &
-                                               0.0_pReal, &
-                                               pl%get_asFloat('a_cI', defaultVal=0.0_pReal), &
-                                               pl%get_asFloat('a_cF', defaultVal=0.0_pReal))
 
       if (phase_lattice(ph) /= 'cF') then
         prm%dot_N_0_tr = pl%get_as1dFloat('dot_N_0_tr')
-        prm%dot_N_0_tr = math_expand(prm%dot_N_0_tr,N_tr)
+        prm%dot_N_0_tr = math_expand(prm%dot_N_0_tr,prm%N_tr)
       endif
       prm%t_tr = pl%get_as1dFloat('t_tr')
-      prm%t_tr = math_expand(prm%t_tr,N_tr)
+      prm%t_tr = math_expand(prm%t_tr,prm%N_tr)
       prm%s    = pl%get_as1dFloat('p_tr',defaultVal=[0.0_pReal])
-      prm%s    = math_expand(prm%s,N_tr)
+      prm%s    = math_expand(prm%s,prm%N_tr)
 
       ! sanity checks
       if (    prm%x_c_tr        < 0.0_pReal)  extmsg = trim(extmsg)//' x_c_tr'
@@ -386,15 +378,15 @@ module function plastic_dislotwin_init() result(myPlasticity)
     end if
 
     slipAndTwinActive: if (prm%sum_N_sl * prm%sum_N_tw > 0) then
-      prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,N_tw,pl%get_as1dFloat('h_sl-tw'), &
+      prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,prm%N_tw,pl%get_as1dFloat('h_sl-tw'), &
                                                    phase_lattice(ph))
-      if (prm%fccTwinTransNucleation .and. size(N_tw) /= 1) extmsg = trim(extmsg)//' N_tw: nucleation'
+      if (prm%fccTwinTransNucleation .and. size(prm%N_tw) /= 1) extmsg = trim(extmsg)//' N_tw: nucleation'
     endif slipAndTwinActive
 
     slipAndTransActive: if (prm%sum_N_sl * prm%sum_N_tr > 0) then
-      prm%h_sl_tr = lattice_interaction_SlipByTrans(N_sl,N_tr,pl%get_as1dFloat('h_sl-tr'), &
+      prm%h_sl_tr = lattice_interaction_SlipByTrans(N_sl,prm%N_tr,pl%get_as1dFloat('h_sl-tr'), &
                                                     phase_lattice(ph))
-      if (prm%fccTwinTransNucleation .and. size(N_tr) /= 1) extmsg = trim(extmsg)//' N_tr: nucleation'
+      if (prm%fccTwinTransNucleation .and. size(prm%N_tr) /= 1) extmsg = trim(extmsg)//' N_tr: nucleation'
     endif slipAndTransActive
 
 !--------------------------------------------------------------------------------------------------
@@ -478,27 +470,40 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC)
   integer,     intent(in) :: &
     ph, en
   real(pReal), dimension(6,6) :: &
-    homogenizedC
+    homogenizedC, &
-
+    C
+  real(pReal), dimension(:,:,:), allocatable :: &
+    C66_tw, &
+    C66_tr
   integer :: i
   real(pReal) :: f_unrotated
 
 
+  C = elastic_C66(ph,en)
+
   associate(prm => param(ph), stt => state(ph))
 
     f_unrotated = 1.0_pReal &
                 - sum(stt%f_tw(1:prm%sum_N_tw,en)) &
                 - sum(stt%f_tr(1:prm%sum_N_tr,en))
 
-    homogenizedC = f_unrotated * prm%C66
+    homogenizedC = f_unrotated * C
+
+    twinActive: if (prm%sum_N_tw > 0) then
+      C66_tw    = lattice_C66_twin(prm%N_tw,C,phase_lattice(ph),phase_cOverA(ph))
       do i=1,prm%sum_N_tw
         homogenizedC = homogenizedC &
-                   + stt%f_tw(i,en)*prm%C66_tw(1:6,1:6,i)
+                     + stt%f_tw(i,en)*C66_tw(1:6,1:6,i)
       end do
+     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
         homogenizedC = homogenizedC &
-                   + stt%f_tr(i,en)*prm%C66_tr(1:6,1:6,i)
+                     + stt%f_tr(i,en)*C66_tr(1:6,1:6,i)
       end do
+    end if transActive
 
   end associate
 
@@ -589,7 +594,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
 
   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
@@ -647,10 +652,15 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
     dot_gamma_tw
   real(pReal), dimension(param(ph)%sum_N_tr) :: &
     dot_gamma_tr
-
+  real(pReal) :: &
+    mu, &
+    nu
 
   associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph))
 
+    mu = elastic_mu(ph,en)
+    nu = elastic_nu(ph,en)
+
     f_unrotated = 1.0_pReal &
                 - sum(stt%f_tw(1:prm%sum_N_tw,en)) &
                 - sum(stt%f_tr(1:prm%sum_N_tr,en))
@@ -665,7 +675,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
         dot_rho_dip_formation(i) = 0.0_pReal
         dot_rho_dip_climb(i) = 0.0_pReal
       else significantSlipStress
-        d_hat = 3.0_pReal*prm%mu*prm%b_sl(i)/(16.0_pReal*PI*abs(tau))
+        d_hat = 3.0_pReal*mu*prm%b_sl(i)/(16.0_pReal*PI*abs(tau))
         d_hat = math_clip(d_hat, right = dst%Lambda_sl(i,en))
         d_hat = math_clip(d_hat, left  = prm%d_caron(i))
 
@@ -677,12 +687,12 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
         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) &
+          b_d = merge(24.0_pReal*PI*(1.0_pReal - nu)/(2.0_pReal + nu) &
-                        * (prm%Gamma_sf(1) + prm%Gamma_sf(2) * T) / (prm%mu*prm%b_sl(i)), &
+                        * (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))
@@ -732,10 +742,15 @@ module subroutine dislotwin_dependentState(T,ph,en)
     f_over_t_tr
   real(pReal), dimension(:), allocatable :: &
     x0
-
+  real(pReal) :: &
+    mu, &
+    nu
 
   associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
 
+    mu = elastic_mu(ph,en)
+    nu = elastic_nu(ph,en)
+
     sumf_tw = sum(stt%f_tw(1:prm%sum_N_tw,en))
     sumf_tr = sum(stt%f_tr(1:prm%sum_N_tr,en))
 
@@ -760,24 +775,24 @@ module subroutine dislotwin_dependentState(T,ph,en)
     dst%Lambda_tr(:,en) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr)
 
     !* threshold stress for dislocation motion
-    dst%tau_pass(:,en) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
+    dst%tau_pass(:,en) = mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
 
     !* threshold stress for growing twin/martensite
     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*mu/(prm%L_tw*prm%b_tw)
     dst%tau_hat_tr(:,en) = Gamma/(3.0_pReal*prm%b_tr) &
-                         + 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_tr) &
+                         + 3.0_pReal*prm%b_tr*mu/(prm%L_tr*prm%b_tr) &
                          + 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
+    x0 = mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+nu)/(1.0_pReal-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)
+    dst%tau_r_tw(:,en) = 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
+    x0 = mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+nu)/(1.0_pReal-nu)                  ! ToDo: In the paper, this is the Burgers vector for slip
-    dst%tau_r_tr(:,en) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0)
+    dst%tau_r_tr(:,en) = mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0)
 
   end associate
 
@@ -889,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)
@@ -962,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
@@ -1031,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

src/phase_mechanical_plastic_isotropic.f90

@@ -68,11 +68,11 @@ module function plastic_isotropic_init() result(myPlasticity)
   myPlasticity = plastic_active('isotropic')
   if(count(myPlasticity) == 0) return
 
-  print'(/,a)', ' <<<+-  phase:mechanical:plastic:isotropic init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:mechanical:plastic:isotropic init  -+>>>'
-  print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
+  print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
 
-  print*, 'T. Maiti and P. Eisenlohr, Scripta Materialia 145:37–40, 2018'
+  print'(/,a)', 'T. Maiti and P. Eisenlohr, Scripta Materialia 145:37–40, 2018'
-  print*, 'https://doi.org/10.1016/j.scriptamat.2017.09.047'
+  print'(/,a)', 'https://doi.org/10.1016/j.scriptamat.2017.09.047'
 
   phases => config_material%get('phase')
   allocate(param(phases%length))

src/phase_mechanical_plastic_kinehardening.f90

@@ -83,8 +83,8 @@ module function plastic_kinehardening_init() result(myPlasticity)
   myPlasticity = plastic_active('kinehardening')
   if(count(myPlasticity) == 0) return
 
-  print'(/,a)', ' <<<+-  phase:mechanical:plastic:kinehardening init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:mechanical:plastic:kinehardening init  -+>>>'
-  print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
+  print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
 
 
   phases => config_material%get('phase')

src/phase_mechanical_plastic_none.f90

@@ -24,8 +24,8 @@ module function plastic_none_init() result(myPlasticity)
   myPlasticity = plastic_active('none')
   if (count(myPlasticity) == 0) return
 
-  print'(/,a)', ' <<<+-  phase:mechanical:plastic:none init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:mechanical:plastic:none init  -+>>>'
-  print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
+  print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
 
   phases => config_material%get('phase')
   do ph = 1, phases%length

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)
@@ -202,14 +199,14 @@ module function plastic_nonlocal_init() result(myPlasticity)
     return
   end if
 
-  print'(/,a)', ' <<<+-  phase:mechanical:plastic:nonlocal init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:mechanical:plastic:nonlocal init  -+>>>'
-  print'(a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT)
+  print'(/,a,i0)', ' # phases: ',Ninstances; flush(IO_STDOUT)
 
-  print*, 'C. Reuber et al., Acta Materialia 71:333–348, 2014'
+  print'(/,1x,a)', 'C. Reuber et al., Acta Materialia 71:333–348, 2014'
-  print*, 'https://doi.org/10.1016/j.actamat.2014.03.012'//IO_EOL
+  print'(  1x,a)', 'https://doi.org/10.1016/j.actamat.2014.03.012'//IO_EOL
 
-  print*, 'C. Kords, Dissertation RWTH Aachen, 2014'
+  print'(/,1x,a)', 'C. Kords, Dissertation RWTH Aachen, 2014'
-  print*, 'http://publications.rwth-aachen.de/record/229993'
+  print'(  1x,a)', 'http://publications.rwth-aachen.de/record/229993'
 
 
   phases => config_material%get('phase')
@@ -242,9 +239,6 @@ module function plastic_nonlocal_init() result(myPlasticity)
 
     prm%atol_rho = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
 
-    prm%mu = elastic_mu(ph)
-    prm%nu = elastic_nu(ph)
-
     ini%N_sl     = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
     prm%sum_N_sl = sum(abs(ini%N_sl))
     slipActive: if (prm%sum_N_sl > 0) then
@@ -570,7 +564,9 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
     n
   real(pReal) :: &
     FVsize, &
-    nRealNeighbors                                                                                  ! number of really existing neighbors
+    nRealNeighbors, &                                                                               ! number of really existing neighbors
+    mu, &
+    nu
   integer, dimension(2) :: &
     neighbors
   real(pReal), dimension(2) :: &
@@ -609,6 +605,8 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
 
   associate(prm => param(ph),dst => dependentState(ph), stt => state(ph))
 
+  mu = elastic_mu(ph,en)
+  nu = elastic_nu(ph,en)
   rho = getRho(ph,en)
 
   stt%rho_forest(:,en) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) &
@@ -627,7 +625,7 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
     myInteractionMatrix = prm%h_sl_sl
   end if
 
-  dst%tau_pass(:,en) = prm%mu * prm%b_sl &
+  dst%tau_pass(:,en) = mu * prm%b_sl &
                      * sqrt(matmul(myInteractionMatrix,sum(abs(rho),2)))
 
 !*** calculate the dislocation stress of the neighboring excess dislocation densities
@@ -728,8 +726,8 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
       where(rhoTotal > 0.0_pReal) rhoExcessGradient_over_rho = rhoExcessGradient / rhoTotal
 
         ! ... gives the local stress correction when multiplied with a factor
-      dst%tau_back(s,en) = - prm%mu * prm%b_sl(s) / (2.0_pReal * PI) &
+      dst%tau_back(s,en) = - mu * prm%b_sl(s) / (2.0_pReal * PI) &
-                         * (  rhoExcessGradient_over_rho(1) / (1.0_pReal - prm%nu) &
+                         * (  rhoExcessGradient_over_rho(1) / (1.0_pReal - nu) &
                             + rhoExcessGradient_over_rho(2))
     end do
   end if
@@ -855,6 +853,9 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
     c, &                                                                                            ! character of dislocation
     t, &                                                                                            ! type of dislocation
     s                                                                                               ! index of my current slip system
+  real(pReal) :: &
+    mu, &
+    nu
   real(pReal), dimension(param(ph)%sum_N_sl,10) :: &
     deltaRhoRemobilization, &                                                                       ! density increment by remobilization
     deltaRhoDipole2SingleStress                                                                     ! density increment by dipole dissociation (by stress change)
@@ -872,6 +873,9 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
    
   associate(prm => param(ph),dst => dependentState(ph),del => deltaState(ph))
 
+  mu = elastic_mu(ph,en)
+  nu = elastic_nu(ph,en)
+
   !*** shortcut to state variables
   forall (s = 1:prm%sum_N_sl, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en)
   forall (s = 1:prm%sum_N_sl, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,ph),en)
@@ -901,8 +905,8 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
     if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal
   end do
 
-  dUpper(:,1) = prm%mu * prm%b_sl/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau))
+  dUpper(:,1) = mu * prm%b_sl/(8.0_pReal * PI * (1.0_pReal - nu) * abs(tau))
-  dUpper(:,2) = prm%mu * prm%b_sl/(4.0_pReal * PI * abs(tau))
+  dUpper(:,2) = mu * prm%b_sl/(4.0_pReal * PI * abs(tau))
 
   where(dNeq0(sqrt(sum(abs(rho(:,edg)),2)))) &
     dUpper(:,1) = min(1.0_pReal/sqrt(sum(abs(rho(:,edg)),2)),dUpper(:,1))
@@ -975,7 +979,9 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
     dLower, &                                                                                       !< minimum stable dipole distance for edges and screws
     dUpper                                                                                          !< current maximum stable dipole distance for edges and screws
   real(pReal) :: &
-    D_SD
+    D_SD, &
+    mu, &
+    nu
 
   if (timestep <= 0.0_pReal) then
     plasticState(ph)%dotState = 0.0_pReal
@@ -984,6 +990,9 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
 
   associate(prm => param(ph), dst => dependentState(ph), dot => dotState(ph), stt => state(ph))
 
+  mu = elastic_mu(ph,en)
+  nu = elastic_nu(ph,en)
+
   tau = 0.0_pReal
   dot_gamma = 0.0_pReal
 
@@ -1005,8 +1014,8 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
   end do
 
   dLower = prm%minDipoleHeight
-  dUpper(:,1) = prm%mu * prm%b_sl/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau))
+  dUpper(:,1) = mu * prm%b_sl/(8.0_pReal * PI * (1.0_pReal - nu) * abs(tau))
-  dUpper(:,2) = prm%mu * prm%b_sl/(4.0_pReal * PI * abs(tau))
+  dUpper(:,2) = mu * prm%b_sl/(4.0_pReal * PI * abs(tau))
 
   where(dNeq0(sqrt(sum(abs(rho(:,edg)),2)))) &
     dUpper(:,1) = min(1.0_pReal/sqrt(sum(abs(rho(:,edg)),2)),dUpper(:,1))
@@ -1082,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 * prm%mu * prm%V_at &
+  v_climb = D_SD * mu * prm%V_at &
-          / (PI * (1.0_pReal-prm%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) &
@@ -1659,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)
@@ -1673,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)

src/phase_mechanical_plastic_phenopowerlaw.f90

@@ -95,8 +95,8 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
   myPlasticity = plastic_active('phenopowerlaw')
   if(count(myPlasticity) == 0) return
 
-  print'(/,a)', ' <<<+-  phase:mechanical:plastic:phenopowerlaw init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:mechanical:plastic:phenopowerlaw init  -+>>>'
-  print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
+  print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
 
 
   phases => config_material%get('phase')

src/phase_thermal.f90

@@ -86,7 +86,7 @@ module subroutine thermal_init(phases)
     Nmembers
 
 
-  print'(/,a)', ' <<<+-  phase:thermal init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:thermal init  -+>>>'
 
   allocate(current(phases%length))
 

src/phase_thermal_dissipation.f90

@@ -36,8 +36,8 @@ module function dissipation_init(source_length) result(mySources)
 
   mySources = thermal_active('dissipation',source_length)
   if(count(mySources) == 0) return
-  print'(/,a)', ' <<<+-  phase:thermal:dissipation init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:thermal:dissipation init  -+>>>'
-  print'(a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT)
+  print'(/,a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT)
 
 
   phases => config_material%get('phase')

src/phase_thermal_externalheat.f90

@@ -43,8 +43,8 @@ module function externalheat_init(source_length) result(mySources)
 
   mySources = thermal_active('externalheat',source_length)
   if(count(mySources) == 0) return
-  print'(/,a)', ' <<<+-  phase:thermal:externalheat init  -+>>>'
+  print'(/,1x,a)', '<<<+-  phase:thermal:externalheat init  -+>>>'
-  print'(a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT)
+  print'(/,a,i2)', ' # phases: ',count(mySources); flush(IO_STDOUT)
 
 
   phases => config_material%get('phase')

src/prec.f90

@@ -69,15 +69,15 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine prec_init
 
-  print'(/,a)', ' <<<+-  prec init  -+>>>'
+  print'(/,1x,a)', '<<<+-  prec init  -+>>>'
 
-  print'(a,i3)',    ' Size of integer in bit: ',bit_size(0)
+  print'(/,a,i3)',    ' integer size / bit:   ',bit_size(0)
-  print'(a,i19)',   '   Maximum value:        ',huge(0)
+  print'(  a,i19)',   '   maximum value:      ',huge(0)
-  print'(/,a,i3)',  ' Size of float in bit:   ',storage_size(0.0_pReal)
+  print'(/,a,i3)',    ' float size / bit:     ',storage_size(0.0_pReal)
-  print'(a,e10.3)', '   Maximum value:        ',huge(0.0_pReal)
+  print'(  a,e10.3)', '   maximum value:      ',huge(0.0_pReal)
-  print'(a,e10.3)', '   Minimum value:        ',PREAL_MIN
+  print'(  a,e10.3)', '   minimum value:      ',PREAL_MIN
-  print'(a,e10.3)', '   Epsilon value:        ',PREAL_EPSILON
+  print'(  a,e10.3)', '   epsilon value:      ',PREAL_EPSILON
-  print'(a,i3)',    '   Decimal precision:    ',precision(0.0_pReal)
+  print'(  a,i3)',    '   decimal precision:  ',precision(0.0_pReal)
 
   call selfTest
 

src/results.f90

@@ -70,10 +70,10 @@ subroutine results_init(restart)
   character(len=:), allocatable :: date
 
 
-  print'(/,a)', ' <<<+-  results init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  results init  -+>>>'; flush(IO_STDOUT)
 
-  print*, 'M. Diehl et al., Integrating Materials and Manufacturing Innovation 6(1):83–91, 2017'
+  print'(/,1x,a)', '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
+  print'(  1x,a)', 'https://doi.org/10.1007/s40192-017-0084-5'
 
   if (.not. restart) then
     resultsFile = HDF5_openFile(getSolverJobName()//'.hdf5','w')

src/rotations.f90

@@ -75,7 +75,6 @@ module rotations
       procedure, public  :: rotVector
       procedure, public  :: rotTensor2
       procedure, public  :: rotTensor4
-      procedure, public  :: rotTensor4sym
       procedure, public  :: misorientation
       procedure, public  :: standardize
   end type rotation
@@ -103,10 +102,10 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine rotations_init
 
-  print'(/,a)', ' <<<+-  rotations init  -+>>>'; flush(IO_STDOUT)
+  print'(/,1x,a)', '<<<+-  rotations init  -+>>>'; flush(IO_STDOUT)
 
-  print*, 'D. Rowenhorst et al., Modelling and Simulation in Materials Science and Engineering 23:083501, 2015'
+  print'(/,1x,a)', 'D. Rowenhorst et al., Modelling and Simulation in Materials Science and Engineering 23:083501, 2015'
-  print*, 'https://doi.org/10.1088/0965-0393/23/8/083501'
+  print'(  1x,a)', 'https://doi.org/10.1088/0965-0393/23/8/083501'
 
   call selfTest
 
@@ -371,27 +370,6 @@ pure function rotTensor4(self,T,active) result(tRot)
 end function rotTensor4
 
 
-!---------------------------------------------------------------------------------------------------
-!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief rotate a symmetric rank-4 tensor stored as (6,6) passively (default) or actively
-!! ToDo: Need to check active/passive !!!
-!---------------------------------------------------------------------------------------------------
-pure function rotTensor4sym(self,T,active) result(tRot)
-
-  real(pReal),                 dimension(6,6) :: tRot
-  class(rotation), intent(in)                 :: self
-  real(pReal),     intent(in), dimension(6,6) :: T
-  logical,         intent(in), optional       :: active
-
-  if (present(active)) then
-    tRot = math_sym3333to66(rotTensor4(self,math_66toSym3333(T),active))
-  else
-    tRot = math_sym3333to66(rotTensor4(self,math_66toSym3333(T)))
-  endif
-
-end function rotTensor4sym
-
-
 !---------------------------------------------------------------------------------------------------
 !> @brief misorientation
 !---------------------------------------------------------------------------------------------------
@@ -400,6 +378,7 @@ pure elemental function misorientation(self,other)
   type(rotation)              :: misorientation
   class(rotation), intent(in) :: self, other
 
+
   misorientation%q = multiply_quaternion(other%q, conjugate_quaternion(self%q))
 
 end function misorientation