Merge remote-tracking branch 'origin/development' into mesh-clean

CMakeLists.txt

@@ -29,6 +29,8 @@ else()
 endif()
 add_definitions("-D${DAMASK_SOLVER}")
 
+set(CMAKE_Fortran_PREPROCESS "ON") # works only for CMake >= 3.18
+
 # EXPERIMENTAL: This might help to detect HDF5 and FFTW3 in the future if PETSc is not aware of them
 set(ENV{PKG_CONFIG_PATH} "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/externalpackages:$ENV{PKG_CONFIG_PATH}")
 pkg_check_modules(HDF5 hdf5)
@@ -91,24 +93,24 @@ if(CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY")
 endif()
 
 
-list(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake)
 if   (CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
-    include(Compiler-GNU)
     set(Fortran_COMPILER_VERSION_MIN 9.1)
 elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "Intel")
-    include(Compiler-Intel)
     set(Fortran_COMPILER_VERSION_MIN 19)
 elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "IntelLLVM")
-    include(Compiler-IntelLLVM)
+    set(Fortran_COMPILER_VERSION_MIN 19)
+elseif(CMAKE_Fortran_COMPILER_ID STREQUAL "LLVMFlang")
     set(Fortran_COMPILER_VERSION_MIN 19)
 else()
     message(FATAL_ERROR "Compiler '${CMAKE_Fortran_COMPILER_ID}' not supported")
 endif()
-
 if(CMAKE_Fortran_COMPILER_VERSION VERSION_LESS Fortran_COMPILER_VERSION_MIN)
     message(FATAL_ERROR "Version '${CMAKE_Fortran_COMPILER_VERSION}' of '${CMAKE_Fortran_COMPILER_ID}' is not supported")
 endif()
 
+list(APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake)
+include("Compiler-${CMAKE_Fortran_COMPILER_ID}")
+
 file(STRINGS "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/petsc/conf/petscvariables" PETSC_EXTERNAL_LIB REGEX "PETSC_EXTERNAL_LIB_BASIC = .*$?")
 string(REPLACE "PETSC_EXTERNAL_LIB_BASIC = " "" PETSC_EXTERNAL_LIB "${PETSC_EXTERNAL_LIB}")
 message("PETSC_EXTERNAL_LIB:\n${PETSC_EXTERNAL_LIB}\n")

PRIVATE

@@ -1 +1 @@
-Subproject commit 117b65d852659158c0f4ca3bf8ce8db51a7a1961
+Subproject commit 62df7f24f2a95fda255f7d20b130afcfeecb1b4a

VERSION

@@ -1 +1 @@
-3.0.0-alpha8-169-g85191cd02
+3.0.0-alpha8-211-gccf4867e0

cmake/Compiler-GNU.cmake

@@ -1,10 +1,6 @@
 ###################################################################################################
 # GNU Compiler
 ###################################################################################################
-if (CMAKE_Fortran_COMPILER_VERSION VERSION_LESS 9.0)
-  message (FATAL_ERROR "GCC Compiler version: ${CMAKE_Fortran_COMPILER_VERSION} not supported")
-endif ()
-
 if (OPENMP)
   set (OPENMP_FLAGS "-fopenmp")
 endif ()
@@ -23,8 +19,7 @@ set (STANDARD_CHECK "-std=f2018 -pedantic-errors" )
 
 #------------------------------------------------------------------------------------------------
 # Fine tuning compilation options
-set (COMPILE_FLAGS "${COMPILE_FLAGS} -cpp")
+set (COMPILE_FLAGS "${COMPILE_FLAGS} -cpp") # preprocessor, needed for CMake < 3.18
-# preprocessor
 
 set (COMPILE_FLAGS "${COMPILE_FLAGS} -fPIE")
 # position independent code

cmake/Compiler-Intel.cmake

@@ -1,10 +1,6 @@
 ###################################################################################################
 # Intel Compiler
 ###################################################################################################
-if (CMAKE_Fortran_COMPILER_VERSION VERSION_LESS 18.0)
-  message (FATAL_ERROR "Intel Compiler version: ${CMAKE_Fortran_COMPILER_VERSION} not supported")
-endif ()
-
 if (OPENMP)
   set (OPENMP_FLAGS "-qopenmp -parallel")
 endif ()
@@ -26,8 +22,7 @@ set (LINKER_FLAGS   "${LINKER_FLAGS} -shared-intel")
 
 #------------------------------------------------------------------------------------------------
 # Fine tuning compilation options
-set (COMPILE_FLAGS "${COMPILE_FLAGS} -fpp")
+set (COMPILE_FLAGS "${COMPILE_FLAGS} -fpp") # preprocessor, needed for CMake < 3.18
-# preprocessor
 
 set (COMPILE_FLAGS "${COMPILE_FLAGS} -no-ftz")
 # disable flush underflow to zero, will be set if -O[1,2,3]

cmake/Compiler-IntelLLVM.cmake

@@ -1,10 +1,6 @@
 ###################################################################################################
-# Intel Compiler
+# IntelLLVM Compiler
 ###################################################################################################
-if (CMAKE_Fortran_COMPILER_VERSION VERSION_LESS 18.0)
-  message (FATAL_ERROR "Intel Compiler version: ${CMAKE_Fortran_COMPILER_VERSION} not supported")
-endif ()
-
 if (OPENMP)
   set (OPENMP_FLAGS "-fiopenmp")
 endif ()
@@ -28,8 +24,7 @@ set (LINKER_FLAGS   "${LINKER_FLAGS} -shared-intel -fc=ifx")
 
 #------------------------------------------------------------------------------------------------
 # Fine tuning compilation options
-set (COMPILE_FLAGS "${COMPILE_FLAGS} -fpp")
+set (COMPILE_FLAGS "${COMPILE_FLAGS} -fpp") # preprocessor, needed for CMake < 3.18
-# preprocessor
 
 set (COMPILE_FLAGS "${COMPILE_FLAGS} -no-ftz")
 # disable flush underflow to zero, will be set if -O[1,2,3]

cmake/Compiler-LLVMFlang.cmake

@@ -0,0 +1,12 @@
+###################################################################################################
+# LLVM Compiler
+###################################################################################################
+if (OPENMP)
+  set (OPENMP_FLAGS "-fopenmp")
+endif ()
+
+set (STANDARD_CHECK "-std=f2018 -pedantic" )
+
+#------------------------------------------------------------------------------------------------
+# Fine tuning compilation options
+set (COMPILE_FLAGS "${COMPILE_FLAGS} -cpp") # preprocessor, needed for CMake < 3.18

examples/mesh/tensionY_mono.load

@@ -1,14 +0,0 @@
-# initial elastic step
-$Loadcase 1 t 0.0005 N 1 f_out 1
- Face 3 Y -0.025
- Face 4 X 0.0
- Face 4 Y 0.0
- Face 4 Z 0.0
-$EndLoadcase
-# plastic step
-$Loadcase 2 t 1.0 N 10 f_out 2
- Face 3 Y -0.025
- Face 4 X 0.0
- Face 4 Y 0.0
- Face 4 Z 0.0
-$EndLoadcase

examples/mesh/tensionY_mono.yaml

@@ -0,0 +1,22 @@
+---
+loadstep:
+  - boundary_conditions:
+      mechanical:
+        - dot_u: ['x', -0.025, 'x']
+          tag: 3
+        - dot_u: [0.0, 0.0, 0.0]
+          tag: 4
+    discretization:
+      t: 0.0005
+      N: 1
+    f_out: 1
+  - boundary_conditions:
+      mechanical:
+        - dot_u: ['x', -0.025, 'x']
+          tag: 3
+        - dot_u: [0.0, 0.0, 0.0]
+          tag: 4
+    discretization:
+      t: 1.0
+      N: 10
+    f_out: 2

examples/mesh/tensionZ_3g.load

@@ -1,14 +0,0 @@
-# initial elastic step
-$Loadcase 1 t 0.0005 N 1 f_out 1
- Face 1 Z 0.01
- Face 2 X 0.0
- Face 2 Y 0.0
- Face 2 Z 0.0
-$EndLoadcase
-# plastic step
-$Loadcase 2 t 1.0 N 10 f_out 2
- Face 1 Z 0.01
- Face 2 X 0.0
- Face 2 Y 0.0
- Face 2 Z 0.0
-$EndLoadcase

examples/mesh/tensionZ_3g.yaml

@@ -0,0 +1,22 @@
+---
+loadstep:
+  - boundary_conditions:
+      mechanical:
+        - dot_u: ['x', 'x', 0.01]
+          tag: 1
+        - dot_u: [0.0, 0.0, 0.0]
+          tag: 2
+    discretization:
+      t: 0.0005
+      N: 1
+    f_out: 1
+  - boundary_conditions:
+      mechanical:
+        - dot_u: ['x', 'x', 0.01]
+          tag: 1
+        - dot_u: [0.0, 0.0, 0.0]
+          tag: 2
+    discretization:
+      t: 1.0
+      N: 10
+    f_out: 2

examples/mesh/tensionZ_5g.load

@@ -1,18 +0,0 @@
-# initial elastic step
-$Loadcase 1 t 0.0005 N 1 f_out 1
- Face 1 X 0.0
- Face 1 Y 0.0
- Face 1 Z 0.0
- Face 2 X 0.0
- Face 2 Y 0.0
- Face 2 Z 0.0025
-$EndLoadcase
-# plastic step
-$Loadcase 2 t 1.0 N 10 f_out 2
- Face 1 X 0.0
- Face 1 Y 0.0
- Face 1 Z 0.0
- Face 2 X 0.0
- Face 2 Y 0.0
- Face 2 Z 0.0025
-$EndLoadcase

examples/mesh/tensionZ_5g.yaml

@@ -0,0 +1,22 @@
+---
+loadstep:
+  - boundary_conditions:
+      mechanical:
+        - dot_u: [0.0, 0.0, 0.0]
+          tag: 1
+        - dot_u: [0.0, 0.0, 0.0025]
+          tag: 2
+    discretization:
+      t: 0.0005
+      N: 1
+    f_out: 1
+  - boundary_conditions:
+      mechanical:
+        - dot_u: [0.0, 0.0, 0.0]
+          tag: 1
+        - dot_u: [0.0, 0.0, 0.0025]
+          tag: 2
+    discretization:
+      t: 1.0
+      N: 10
+    f_out: 2

python/damask/_colormap.py

@@ -2,7 +2,7 @@ import os
 import json
 import functools
 import colorsys
-from typing import Optional, Union, TextIO
+from typing import Optional, Union
 from itertools import chain
 
 import numpy as np
@@ -344,30 +344,6 @@ 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: Union[FileHandle, None],
-                         suffix: str = '') -> TextIO:
-        """
-        Provide file handle.
-
-        Parameters
-        ----------
-        fname : file, str, pathlib.Path, or None
-            Name or handle of file.
-            If None, colormap name + suffix.
-        suffix: str, optional
-            Extension to use for colormap file.
-            Defaults to empty.
-
-        Returns
-        -------
-        f : file object
-            File handle with write access.
-
-        """
-        return util.open_text(self.name.replace(' ','_')+suffix if fname is None else fname, 'w')
-
-
     def save_paraview(self,
                       fname: Optional[FileHandle] = None):
         """
@@ -387,9 +363,9 @@ class Colormap(mpl.colors.ListedColormap):
                 'RGBPoints':list(chain.from_iterable([(i,*c) for i,c in enumerate(self.colors.round(6))]))
                }]
 
-        fhandle = self._get_file_handle(fname,'.json')
+        with util.open_text(self.name.replace(' ','_')+'.json' if fname is None else fname, 'w') as fhandle:
-        json.dump(out,fhandle,indent=4)
+            json.dump(out,fhandle,indent=4)
-        fhandle.write('\n')
+            fhandle.write('\n')
 
 
     def save_ASCII(self,
@@ -405,7 +381,9 @@ class Colormap(mpl.colors.ListedColormap):
         """
         labels = {'RGBA':4} if self.colors.shape[1] == 4 else {'RGB': 3}
         t = Table(labels,self.colors,[f'Creator: {util.execution_stamp("Colormap")}'])
-        t.save(self._get_file_handle(fname,'.txt'))
+
+        with util.open_text(self.name.replace(' ','_')+'.txt' if fname is None else fname, 'w') as fhandle:
+            t.save(fhandle)
 
 
     def save_GOM(self, fname: Optional[FileHandle] = None):
@@ -425,7 +403,8 @@ class Colormap(mpl.colors.ListedColormap):
                 + ' '.join([f' 0 {c[0]} {c[1]} {c[2]} 255 1' for c in reversed((self.colors*255).astype(np.int64))]) \
                 + '\n'
 
-        self._get_file_handle(fname,'.legend').write(GOM_str)
+        with util.open_text(self.name.replace(' ','_')+'.legend' if fname is None else fname, 'w') as fhandle:
+            fhandle.write(GOM_str)
 
 
     def save_gmsh(self,
@@ -443,7 +422,9 @@ 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)
+
+        with util.open_text(self.name.replace(' ','_')+'.msh' if fname is None else fname, 'w') as fhandle:
+            fhandle.write(gmsh_str)
 
 
     @staticmethod

python/damask/_loadcasegrid.py

@@ -70,9 +70,9 @@ class LoadcaseGrid(YAML):
             if key not in kwargs:
                 kwargs[key] = default
 
-        fhandle = util.open_text(fname,'w')
+        with util.open_text(fname,'w') as fhandle:
-        try:
+            try:
-            fhandle.write(yaml.dump(self,Dumper=MaskedMatrixDumper,**kwargs))
+                fhandle.write(yaml.dump(self,Dumper=MaskedMatrixDumper,**kwargs))
-        except TypeError:                                                                           # compatibility with old pyyaml
+            except TypeError:                                                                       # compatibility with old pyyaml
-            del kwargs['sort_keys']
+                del kwargs['sort_keys']
-            fhandle.write(yaml.dump(self,Dumper=MaskedMatrixDumper,**kwargs))
+                fhandle.write(yaml.dump(self,Dumper=MaskedMatrixDumper,**kwargs))

python/damask/_result.py

@@ -115,8 +115,6 @@ class Result:
                 self.cells  = f['geometry'].attrs['cells']
                 self.size   = f['geometry'].attrs['size']
                 self.origin = f['geometry'].attrs['origin']
-            else:
-                self.add_curl = self.add_divergence = self.add_gradient = None                      # type: ignore
 
             r = re.compile(rf'{prefix_inc}([0-9]+)')
             self.increments = sorted([i for i in f.keys() if r.match(i)],key=util.natural_sort)
@@ -1313,8 +1311,8 @@ class Result:
 
         Notes
         -----
-        This function is only available for structured grids,
+        This function is implemented only for structured grids
-        i.e. fields resulting from the grid solver.
+        with one constituent and a single phase.
 
         """
         def curl(f: DADF5Dataset, size: np.ndarray) -> DADF5Dataset:
@@ -1342,8 +1340,8 @@ class Result:
 
         Notes
         -----
-        This function is only available for structured grids,
+        This function is implemented only for structured grids
-        i.e. fields resulting from the grid solver.
+        with one constituent and a single phase.
 
         """
         def divergence(f: DADF5Dataset, size: np.ndarray) -> DADF5Dataset:
@@ -1371,8 +1369,8 @@ class Result:
 
         Notes
         -----
-        This function is only available for structured grids,
+        This function is implemented only for structured grids
-        i.e. fields resulting from the grid solver.
+        with one constituent and a single phase.
 
         """
         def gradient(f: DADF5Dataset, size: np.ndarray) -> DADF5Dataset:
@@ -1410,13 +1408,13 @@ class Result:
             Arguments parsed to func.
 
         """
-        if len(datasets) != 1 or self.N_constituents != 1:
+        if self.N_constituents != 1 or len(datasets) != 1 or not self.structured:
-            raise NotImplementedError
+            raise NotImplementedError('not a structured grid with one constituent and a single phase')
 
         at_cell_ph,in_data_ph,at_cell_ho,in_data_ho = self._mappings()
 
         increments = self.place(list(datasets.values()),False)
-        if not increments: raise RuntimeError("received invalid dataset")
+        if not increments: raise RuntimeError('received invalid dataset')
         with h5py.File(self.fname, 'a') as f:
             for increment in increments.items():
                 for ty in increment[1].items():
@@ -1722,9 +1720,14 @@ class Result:
             Defaults to False, i.e. the XDMF file expects the
             DADF5 file at a stable relative path.
 
+        Notes
+        -----
+        This function is implemented only for structured grids with
+        one constituent and a single phase.
+
         """
         if self.N_constituents != 1 or len(self.phases) != 1 or not self.structured:
-            raise TypeError('XDMF output requires structured grid with single phase and single constituent.')
+            raise NotImplementedError('not a structured grid with one constituent and a single phase')
 
         attribute_type_map = defaultdict(lambda:'Matrix', ( ((),'Scalar'), ((3,),'Vector'), ((3,3),'Tensor')) )
 
@@ -1949,6 +1952,11 @@ class Result:
         target_dir : str or pathlib.Path, optional
             Directory to save DREAM3D files. Will be created if non-existent.
 
+        Notes
+        -----
+        This function is implemented only for structured grids with
+        one constituent.
+
         """
         def add_attribute(obj,name,data):
             """DREAM.3D requires fixed length string."""
@@ -1964,11 +1972,10 @@ class Result:
             return obj[name]
 
         if self.N_constituents != 1 or not self.structured:
-            raise TypeError('DREAM3D output requires structured grid with single constituent.')
+            raise NotImplementedError('not a structured grid with one constituent')
 
         N_digits = int(np.floor(np.log10(max(1,self.incs[-1]))))+1
 
-
         at_cell_ph,in_data_ph,_,_ = self._mappings()
 
         out_dir = Path.cwd() if target_dir is None else Path(target_dir)

python/damask/_table.py

@@ -277,28 +277,28 @@ class Table:
             Table data from file.
 
         """
-        f = util.open_text(fname)
+        with util.open_text(fname) as f:
-        f.seek(0)
+            f.seek(0)
 
-        comments = []
+            comments = []
-        while (line := f.readline().strip()).startswith('#'):
+            while (line := f.readline().strip()).startswith('#'):
-            comments.append(line.lstrip('#').strip())
+                comments.append(line.lstrip('#').strip())
-        labels = line.split()
+            labels = line.split()
 
-        shapes = {}
+            shapes = {}
-        for label in labels:
+            for label in labels:
-            tensor_column = re.search(r'[0-9,x]*?:[0-9]*?_',label)
+                tensor_column = re.search(r'[0-9,x]*?:[0-9]*?_',label)
-            if tensor_column:
+                if tensor_column:
-                my_shape = tensor_column.group().split(':',1)[0].split('x')
+                    my_shape = tensor_column.group().split(':',1)[0].split('x')
-                shapes[label.split('_',1)[1]] = tuple([int(d) for d in my_shape])
+                    shapes[label.split('_',1)[1]] = tuple([int(d) for d in my_shape])
-            else:
-                vector_column = re.match(r'[0-9]*?_',label)
-                if vector_column:
-                    shapes[label.split('_',1)[1]] = (int(label.split('_',1)[0]),)
                 else:
-                    shapes[label] = (1,)
+                    vector_column = re.match(r'[0-9]*?_',label)
+                    if vector_column:
+                        shapes[label.split('_',1)[1]] = (int(label.split('_',1)[0]),)
+                    else:
+                        shapes[label] = (1,)
 
-        data = pd.read_csv(f,names=list(range(len(labels))),sep=r'\s+')
+            data = pd.read_csv(f,names=list(range(len(labels))),sep=r'\s+')
 
         return Table(shapes,data,comments)
 
@@ -339,10 +339,9 @@ class Table:
             Table data from file.
 
         """
-        f = util.open_text(fname)
+        with util.open_text(fname) as f:
-        f.seek(0)
+            f.seek(0)
-
+            content = f.readlines()
-        content = f.readlines()
 
         comments = [util.execution_stamp('Table','from_ang')]
         for line in content:
@@ -605,10 +604,9 @@ class Table:
                     labels += [f'{util.srepr(self.shapes[l],"x")}:{i+1}_{l}' \
                               for i in range(np.prod(self.shapes[l]))]
 
-        f = util.open_text(fname,'w')
+        with util.open_text(fname,'w') as f:
-
+            f.write('\n'.join([f'# {c}' for c in self.comments] + [' '.join(labels)])+('\n' if labels else ''))
-        f.write('\n'.join([f'# {c}' for c in self.comments] + [' '.join(labels)])+('\n' if labels else ''))
+            try:                                                                                    # backward compatibility
-        try:                                                                                        # backward compatibility
+                self.data.to_csv(f,sep=' ',na_rep='nan',index=False,header=False,lineterminator='\n')
-            self.data.to_csv(f,sep=' ',na_rep='nan',index=False,header=False,lineterminator='\n')
+            except TypeError:
-        except TypeError:
+                self.data.to_csv(f,sep=' ',na_rep='nan',index=False,header=False,line_terminator='\n')
-            self.data.to_csv(f,sep=' ',na_rep='nan',index=False,header=False,line_terminator='\n')

python/damask/_yaml.py

@@ -197,7 +197,9 @@ class YAML(dict):
             YAML from file.
 
         """
-        return cls(yaml.load(util.open_text(fname), Loader=SafeLoader))
+        with util.open_text(fname) as fhandle:
+            return cls(yaml.load(fhandle, Loader=SafeLoader))
+
 
     def save(self,
              fname: FileHandle,
@@ -220,12 +222,12 @@ class YAML(dict):
         if 'sort_keys' not in kwargs:
             kwargs['sort_keys'] = False
 
-        fhandle = util.open_text(fname,'w')
+        with util.open_text(fname,'w') as fhandle:
-        try:
+            try:
-            fhandle.write(yaml.dump(self,Dumper=NiceDumper,**kwargs))
+                fhandle.write(yaml.dump(self,Dumper=NiceDumper,**kwargs))
-        except TypeError:                                                                           # compatibility with old pyyaml
+            except TypeError:                                                                       # compatibility with old pyyaml
-            del kwargs['sort_keys']
+                del kwargs['sort_keys']
-            fhandle.write(yaml.dump(self,Dumper=NiceDumper,**kwargs))
+                fhandle.write(yaml.dump(self,Dumper=NiceDumper,**kwargs))
 
 
     @property

python/damask/util.py

@@ -8,12 +8,13 @@ import shlex as _shlex
 import re as _re
 import signal as _signal
 import fractions as _fractions
+import contextlib as _contextlib
 from collections import abc as _abc, OrderedDict as _OrderedDict
 from functools import reduce as _reduce, partial as _partial, wraps as _wraps
 import inspect
 from typing import Optional as _Optional, Callable as _Callable, Union as _Union, Iterable as _Iterable, \
                    Dict as _Dict, List as _List, Tuple as _Tuple, Literal as _Literal, \
-                   Any as _Any, TextIO as _TextIO
+                   Any as _Any, TextIO as _TextIO, Generator as _Generator
 from pathlib import Path as _Path
 
 import numpy as _np
@@ -193,11 +194,15 @@ def run(cmd: str,
 
     return stdout, stderr
 
-
+@_contextlib.contextmanager
 def open_text(fname: _FileHandle,
-              mode: _Literal['r','w'] = 'r') -> _TextIO:                                            # noqa
+              mode: _Literal['r','w'] = 'r') -> _Generator[_TextIO, None, None]:                    # noqa
     """
-    Open a text file.
+    Open a text file with Unix line endings
+
+    If a path or string is given, a context manager ensures that
+    the file handle is closed.
+    If a file handle is given, it remains unmodified.
 
     Parameters
     ----------
@@ -211,8 +216,12 @@ def open_text(fname: _FileHandle,
     f : file handle
 
     """
-    return fname if not isinstance(fname, (str,_Path)) else \
+    if isinstance(fname, (str,_Path)):
-           open(_Path(fname).expanduser(),mode,newline=('\n' if mode == 'w' else None))
+        fhandle = open(_Path(fname).expanduser(),mode,newline=('\n' if mode == 'w' else None))
+        yield fhandle
+        fhandle.close()
+    else:
+        yield fname
 
 
 def execution_stamp(class_name: str,
@@ -618,7 +627,7 @@ def _docstringer(docstring: _Union[str, _Callable],
                                           adopted_,
                                           flags=_re.MULTILINE|_re.DOTALL).group('content')
             except AttributeError:
-                raise RuntimeError(f"Function docstring passed for docstring section '{key}' is invalid:\n{docstring}")
+                raise RuntimeError(f"function docstring passed for docstring section '{key}' is invalid:\n{docstring}")
 
             docstring_indent, adopted_indent = (min([len(line)-len(line.lstrip()) for line in section.split('\n') if line.strip()])
                                                 for section in [docstring_, adopted_])

python/tests/test_Result.py

@@ -63,7 +63,7 @@ def h5py_dataset_iterator():
     """Iterate over all datasets in an HDF5 file."""
     def _h5py_dataset_iterator(g, prefix=''):
         for key,item in g.items():
-            path = os.path.join(prefix, key)
+            path = '/'.join([prefix, key])
             if isinstance(item, h5py.Dataset): # test for dataset
                 yield (path, item)
             elif isinstance(item, h5py.Group): # test for group (go down)
@@ -472,12 +472,16 @@ class TestResult:
                     c = [_.decode() for _ in cur[path]]
                     r = ['Unknown Phase Type'] + result.phases
                     assert c == r
-                grp = os.path.split(path)[0]
+                grp = str(path).rpartition('/')[0]
                 for attr in ref[grp].attrs:
                     assert np.array_equal(ref[grp].attrs[attr],cur[grp].attrs[attr])
                 for attr in dset.attrs:
                     assert np.array_equal(dset.attrs[attr],cur[path].attrs[attr])
 
+    def test_export_DREAM3D_invalid(self,res_path):
+        with pytest.raises(NotImplementedError):
+            Result(res_path/'4grains2x4x3_compressionY.hdf5').export_DREAM3D()
+
 
     def test_XDMF_datatypes(self,tmp_path,single_phase,update,res_path):
         for what,shape in {'scalar':(),'vector':(3,),'tensor':(3,3),'matrix':(12,)}.items():
@@ -509,7 +513,7 @@ class TestResult:
         assert dim_vti == dim_xdmf and bounds_vti == bounds_xdmf
 
     def test_XDMF_invalid(self,default):
-        with pytest.raises(TypeError):
+        with pytest.raises(NotImplementedError):
             default.export_XDMF()
 
     def test_XDMF_custom_path(self,single_phase,tmp_path):

src/IO.f90

@@ -48,7 +48,8 @@ implicit none(type,external)
     IO_color, &
     IO_error, &
     IO_warning, &
-    IO_STDOUT
+    IO_STDOUT, &
+    tokenize
 
 contains
 
@@ -380,17 +381,11 @@ integer function IO_strAsInt(str)
 
   character(len=*), intent(in) :: str                                                               !< string for conversion to int value
 
-  integer                      :: readStatus
+  integer :: readStatus
-  character(len=*), parameter  :: VALIDCHARS = '0123456789+- '
 
 
-  valid: if (verify(str,VALIDCHARS) == 0) then
+  read(str,*,iostat=readStatus) IO_strAsInt
-    read(str,*,iostat=readStatus) IO_strAsInt
+  if (readStatus /= 0) call IO_error(111,'cannot represent "'//str//'" as integer')
-    if (readStatus /= 0) call IO_error(111,str)
-  else valid
-    IO_strAsInt = 0
-    call IO_error(111,str)
-  end if valid
 
 end function IO_strAsInt
 
@@ -402,27 +397,23 @@ real(pREAL) function IO_strAsReal(str)
 
   character(len=*), intent(in) :: str                                                               !< string for conversion to real value
 
-  integer                      :: readStatus
+  integer :: readStatus
-  character(len=*), parameter  :: VALIDCHARS = '0123456789eE.+- '
 
 
-  valid: if (verify(str,VALIDCHARS) == 0) then
+  read(str,*,iostat=readStatus) IO_strAsReal
-    read(str,*,iostat=readStatus) IO_strAsReal
+  if (readStatus /= 0) call IO_error(111,'cannot represent "'//str//'" as real')
-    if (readStatus /= 0) call IO_error(112,str)
-  else valid
-    IO_strAsReal = 0.0_pREAL
-    call IO_error(112,str)
-  end if valid
 
 end function IO_strAsReal
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief Return logical value from given string.
+!> @details: 'True' and 'true' are converted to .true.
+!> @details: 'False' and 'false' are converted to .false.
 !--------------------------------------------------------------------------------------------------
 logical function IO_strAsBool(str)
 
-  character(len=*), intent(in) :: str                                                               !< string for conversion to int value
+  character(len=*), intent(in) :: str                                                               !< string for conversion to boolean
 
 
   if     (trim(adjustl(str)) == 'True' .or.  trim(adjustl(str)) == 'true') then
@@ -430,8 +421,7 @@ logical function IO_strAsBool(str)
   elseif (trim(adjustl(str)) == 'False' .or. trim(adjustl(str)) == 'false') then
     IO_strAsBool = .false.
   else
-    IO_strAsBool = .false.
+    call IO_error(111,'cannot represent "'//str//'" as boolean')
-    call IO_error(113,str)
   end if
 
 end function IO_strAsBool
@@ -498,11 +488,7 @@ subroutine IO_error(error_ID,ext_msg,label1,ID1,label2,ID2)
     case (110)
       msg = 'invalid chunk selected'
     case (111)
-      msg = 'invalid character for int:'
+      msg = 'invalid string for conversion'
-    case (112)
-      msg = 'invalid character for real:'
-    case (113)
-      msg = 'invalid character for logical:'
     case (114)
       msg = 'cannot decode base64 string:'
 
@@ -742,6 +728,33 @@ pure function CRLF2LF(str)
 end function CRLF2LF
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief Fortran 2023 tokenize (first form).
+!--------------------------------------------------------------------------------------------------
+pure subroutine tokenize(string,set,tokens)
+
+  character(len=*), intent(in) :: string, set
+  character(len=:), dimension(:), allocatable, intent(out) :: tokens
+
+  integer, allocatable, dimension(:,:) :: pos
+  integer :: i, s, e
+
+
+  allocate(pos(2,0))
+  e = 0
+  do while (e < verify(string,set,back=.true.))
+    s = e + merge(verify(string(e+1:),set),1,scan(string(e+1:),set)/=0)
+    e = s + merge(scan(string(s:),set)-2,len(string(s:))-1,scan(string(s:),set)/=0)
+    pos = reshape([pos,[s,e]],[2,size(pos)/2+1])
+  end do
+  allocate(character(len=merge(maxval(pos(2,:)-pos(1,:))+1,0,size(pos)>0))::tokens(size(pos,2)))
+  do i = 1, size(pos,2)
+    tokens(i) = string(pos(1,i):pos(2,i))
+  end do
+
+end subroutine tokenize
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief Write statements to standard error.
 !--------------------------------------------------------------------------------------------------
@@ -808,6 +821,7 @@ subroutine IO_selfTest()
 
   integer, dimension(:), allocatable :: chunkPos
   character(len=:),      allocatable :: str,out
+  character(len=:), dimension(:), allocatable :: tokens
 
 
   if (dNeq(1.0_pREAL, IO_strAsReal('1.0')))          error stop 'IO_strAsReal'
@@ -887,6 +901,54 @@ subroutine IO_selfTest()
   if ('abc,'//IO_EOL//'xxdefg,'//IO_EOL//'xxhij' /= IO_wrapLines('abc,defg, hij',filler='xx',length=4)) &
                                                      error stop 'IO_wrapLines/7'
 
+  call tokenize('','$',tokens)
+  if (size(tokens) /= 0 .or. len(tokens) /=0) error stop 'tokenize empty'
+  call tokenize('abcd','dcba',tokens)
+  if (size(tokens) /= 0 .or. len(tokens) /=0) error stop 'tokenize only separators'
+
+  tokens=['a']
+  call test_tokenize('a','#',tokens)
+  call test_tokenize('#a','#',tokens)
+  call test_tokenize('a#','#',tokens)
+
+  tokens=['aa']
+  call test_tokenize('aa','#',tokens)
+  call test_tokenize('$aa','$',tokens)
+  call test_tokenize('aa$','$',tokens)
+
+  tokens=['a','b']
+  call test_tokenize('a$b','$',tokens)
+  call test_tokenize('@a@$b@','$@',tokens)
+
+  tokens=['aa','bb']
+  call test_tokenize('aa$bb','$',tokens)
+  call test_tokenize('aa$$bb','$',tokens)
+  call test_tokenize('aa$bb$','$',tokens)
+
+  tokens=['aa  ','bbb ','cccc']
+  call test_tokenize('aa$bbb$cccc','$',tokens)
+  call test_tokenize('$aa$bbb$cccc$','$',tokens)
+  call tokenize('#aa@@bbb!!!cccc#','#@!',tokens)
+
+
+  contains
+  subroutine test_tokenize(input,delimiter,solution)
+
+    character(len=*), intent(in) :: input, delimiter
+    character(len=*), dimension(:), intent(in) :: solution
+
+    character(len=:), dimension(:), allocatable :: tok
+    integer :: i
+
+
+    call tokenize(input,delimiter,tok)
+    do i = 1,size(tok)
+      !if (solution(i) /= tok(i)) error stop 'tokenize "'//solution(i)//'" vs. "'//tok(i)//'"'      ! requires 2018 standard
+      if (solution(i) /= tok(i)) error stop 'tokenize'
+    end do
+
+  end subroutine test_tokenize
+
 end subroutine IO_selfTest
 
 end module IO

src/Marc/materialpoint_Marc.f90

@@ -30,7 +30,8 @@ module materialpoint_Marc
   real(pREAL), dimension (:,:,:),   allocatable, private :: &
     materialpoint_cs                                                                                !< Cauchy stress
   real(pREAL), dimension (:,:,:,:), allocatable, private :: &
-    materialpoint_dcsdE                                                                             !< Cauchy stress tangent
+    materialpoint_dcsdE, &                                                                          !< Cauchy stress tangent
+    materialpoint_F                                                                                 !< deformation gradient
   real(pREAL), dimension (:,:,:,:), allocatable, private :: &
     materialpoint_dcsdE_knownGood                                                                   !< known good tangent
 
@@ -95,6 +96,7 @@ subroutine materialpoint_init()
 
   print'(/,1x,a)', '<<<+-  materialpoint init  -+>>>'; flush(IO_STDOUT)
 
+  allocate(materialpoint_F(              3,3,discretization_nIPs,discretization_Nelems), source= 0.0_pREAL)
   allocate(materialpoint_cs(               6,discretization_nIPs,discretization_Nelems), source= 0.0_pREAL)
   allocate(materialpoint_dcsdE(          6,6,discretization_nIPs,discretization_Nelems), source= 0.0_pREAL)
   allocate(materialpoint_dcsdE_knownGood(6,6,discretization_nIPs,discretization_Nelems), source= 0.0_pREAL)
@@ -140,10 +142,10 @@ subroutine materialpoint_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip,
   if (iand(mode, materialpoint_RESTOREJACOBIAN) /= 0) &
     materialpoint_dcsde = materialpoint_dcsde_knownGood
 
-  if (iand(mode, materialpoint_AGERESULTS) /= 0) call materialpoint_forward
+  if (iand(mode, materialpoint_AGERESULTS) /= 0) call materialpoint_forward()
 
-    homogenization_F0(1:3,1:3,ce) = ffn
     homogenization_F(1:3,1:3,ce) = ffn1
+    materialpoint_F(1:3,1:3,ip,elCP) = ffn1
 
   if (iand(mode, materialpoint_CALCRESULTS) /= 0) then
 
@@ -154,9 +156,11 @@ subroutine materialpoint_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip,
       materialpoint_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_eye(6)
 
     else validCalculation
-      call homogenization_mechanical_response(dt,(elCP-1)*discretization_nIPs + ip,(elCP-1)*discretization_nIPs + ip)
+      call homogenization_mechanical_response(dt,(elCP-1)*discretization_nIPs + ip, &
+                                                 (elCP-1)*discretization_nIPs + ip)
       if (.not. terminallyIll) &
-        call homogenization_mechanical_response2(dt,[ip,ip],[elCP,elCP])
+        call homogenization_mechanical_response2(dt,(elCP-1)*discretization_nIPs + ip, &
+                                                    (elCP-1)*discretization_nIPs + ip)
 
       terminalIllness: if (terminallyIll) then
 
@@ -168,17 +172,17 @@ subroutine materialpoint_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip,
       else terminalIllness
 
         ! translate from P to sigma
-        Kirchhoff = matmul(homogenization_P(1:3,1:3,ce), transpose(homogenization_F(1:3,1:3,ce)))
+        Kirchhoff = matmul(homogenization_P(1:3,1:3,ce), transpose(materialpoint_F(1:3,1:3,ip,elCP)))
-        J_inverse  = 1.0_pREAL / math_det33(homogenization_F(1:3,1:3,ce))
+        J_inverse  = 1.0_pREAL / math_det33(materialpoint_F(1:3,1:3,ip,elCP))
         materialpoint_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.)
 
         !  translate from dP/dF to dCS/dE
         H = 0.0_pREAL
         do i=1,3; do j=1,3; do k=1,3; do l=1,3; do m=1,3; do n=1,3
           H(i,j,k,l) = H(i,j,k,l) &
-                     +  homogenization_F(j,m,ce) * homogenization_F(l,n,ce) &
+                     +  materialpoint_F(j,m,ip,elCP) * materialpoint_F(l,n,ip,elCP) &
-                                                 * homogenization_dPdF(i,m,k,n,ce) &
+                                                     * homogenization_dPdF(i,m,k,n,ce) &
-                     -  math_delta(j,l) * homogenization_F(i,m,ce) * homogenization_P(k,m,ce) &
+                     -  math_delta(j,l) * materialpoint_F(i,m,ip,elCP) * homogenization_P(k,m,ce) &
                      +  0.5_pREAL * (  Kirchhoff(j,l)*math_delta(i,k) + Kirchhoff(i,k)*math_delta(j,l) &
                                      + Kirchhoff(j,k)*math_delta(i,l) + Kirchhoff(i,l)*math_delta(j,k))
         end do; end do; end do; end do; end do; end do

src/grid/DAMASK_grid.f90

@@ -23,6 +23,7 @@ program DAMASK_grid
   use materialpoint
   use material
   use spectral_utilities
+  use grid_mech_utilities
   use grid_mechanical_spectral_basic
   use grid_mechanical_spectral_polarization
   use grid_mechanical_FEM

src/grid/VTI.f90

@@ -201,26 +201,26 @@ subroutine cellsSizeOrigin(c,s,o,header)
   real(pREAL), dimension(3), intent(out) :: s,o
   character(len=*),          intent(in) :: header
 
-  character(len=:), allocatable :: temp
+  character(len=:), allocatable, dimension(:) :: temp
   real(pREAL), dimension(3) :: delta
   integer :: i
 
 
-  temp = getXMLValue(header,'Direction')
+  temp = [getXMLValue(header,'Direction')]
-  if (temp /= '1 0 0 0 1 0 0 0 1' .and. temp /= '') &                                               ! https://discourse.vtk.org/t/vti-specification/6526
+  if (temp(1) /= '1 0 0 0 1 0 0 0 1' .and. temp(1) /= '') &                                         ! https://discourse.vtk.org/t/vti-specification/6526
     call IO_error(error_ID = 844, ext_msg = 'coordinate order')
 
-  temp = getXMLValue(header,'WholeExtent')
+  call tokenize(getXMLValue(header,'WholeExtent'),' ',temp)
-  if (any([(IO_intValue(temp,IO_strPos(temp),i),i=1,5,2)] /= 0)) &
+  if (any([(IO_strAsInt(temp(i)),i=1,5,2)] /= 0)) &
     call IO_error(error_ID = 844, ext_msg = 'coordinate start')
-  c = [(IO_intValue(temp,IO_strPos(temp),i),i=2,6,2)]
+  c = [(IO_strAsInt(temp(i)),i=2,6,2)]
 
-  temp = getXMLValue(header,'Spacing')
+  call tokenize(getXMLValue(header,'Spacing'),' ',temp)
-  delta = [(IO_realValue(temp,IO_strPos(temp),i),i=1,3)]
+  delta = [(IO_strAsReal(temp(i)),i=1,3)]
   s = delta * real(c,pREAL)
 
-  temp = getXMLValue(header,'Origin')
+  call tokenize(getXMLValue(header,'Origin'),' ',temp)
-  o = [(IO_realValue(temp,IO_strPos(temp),i),i=1,3)]
+  o = [(IO_strAsReal(temp(i)),i=1,3)]
 
 end subroutine cellsSizeOrigin
 

src/grid/grid_damage_spectral.f90

@@ -44,7 +44,6 @@ module grid_damage_spectral
 
   type(tNumerics) :: num
 
-  type(tSolutionParams) :: params
 !--------------------------------------------------------------------------------------------------
 ! PETSc data
   SNES :: SNES_damage
@@ -57,7 +56,7 @@ module grid_damage_spectral
 ! reference diffusion tensor, mobility etc.
   integer                     :: totalIter = 0                                                      !< total iteration in current increment
   real(pREAL), dimension(3,3) :: K_ref
-  real(pREAL)                 :: mu_ref
+  real(pREAL)                 :: mu_ref, Delta_t_
 
   public :: &
     grid_damage_spectral_init, &
@@ -207,8 +206,7 @@ end subroutine grid_damage_spectral_init
 !--------------------------------------------------------------------------------------------------
 function grid_damage_spectral_solution(Delta_t) result(solution)
 
-  real(pREAL), intent(in) :: &
+  real(pREAL), intent(in) :: Delta_t                                                                !< increment in time for current solution
-    Delta_t                                                                                         !< increment in time for current solution
 
   type(tSolutionState) :: solution
   PetscInt  :: devNull
@@ -222,7 +220,7 @@ function grid_damage_spectral_solution(Delta_t) result(solution)
 
 !--------------------------------------------------------------------------------------------------
 ! set module wide availabe data
-  params%Delta_t = Delta_t
+  Delta_t_ = Delta_t
 
   call SNESSolve(SNES_damage,PETSC_NULL_VEC,phi_PETSc,err_PETSc)
   CHKERRQ(err_PETSc)
@@ -350,12 +348,12 @@ subroutine formResidual(residual_subdomain,x_scal,r,dummy,err_PETSc)
     ce = 0
     do k = 1, cells3;  do j = 1, cells(2);  do i = 1,cells(1)
       ce = ce + 1
-      r(i,j,k) = params%Delta_t*(r(i,j,k) + homogenization_f_phi(phi(i,j,k),ce)) &
+      r(i,j,k) = Delta_t_*(r(i,j,k) + homogenization_f_phi(phi(i,j,k),ce)) &
                + homogenization_mu_phi(ce)*(phi_lastInc(i,j,k) - phi(i,j,k)) &
                + mu_ref*phi(i,j,k)
     end do; end do; end do
 
-    r = max(min(utilities_GreenConvolution(r, K_ref, mu_ref, params%Delta_t),phi_lastInc),num%phi_min) &
+    r = max(min(utilities_GreenConvolution(r, K_ref, mu_ref, Delta_t_),phi_lastInc),num%phi_min) &
       - phi
   end associate
   err_PETSc = 0

src/grid/grid_mech_FEM.f90

@@ -23,6 +23,7 @@ module grid_mechanical_FEM
   use math
   use rotations
   use spectral_utilities
+  use grid_mech_utilities
   use config
   use homogenization
   use discretization
@@ -269,7 +270,6 @@ subroutine grid_mechanical_FEM_init(num_grid)
     F         = spread(spread(spread(math_I3,3,cells(1)),4,cells(2)),5,cells3)
   end if restartRead
 
-  homogenization_F0 = reshape(F_lastInc, [3,3,product(cells(1:2))*cells3])                          ! set starting condition for homogenization_mechanical_response
   call utilities_updateCoords(F)
   call utilities_constitutiveResponse(P_current,P_av,C_volAvg,devNull, &                            ! stress field, stress avg, global average of stiffness and (min+max)/2
                                       F, &                                                          ! target F
@@ -390,7 +390,6 @@ subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remai
 
     F_lastInc = F
 
-    homogenization_F0 = reshape(F, [3,3,product(cells(1:2))*cells3])
   end if
 
 !--------------------------------------------------------------------------------------------------

src/grid/grid_mech_spectral_basic.f90

@@ -23,6 +23,7 @@ module grid_mechanical_spectral_basic
   use math
   use rotations
   use spectral_utilities
+  use grid_mech_utilities
   use homogenization
   use discretization_grid
 
@@ -226,7 +227,6 @@ subroutine grid_mechanical_spectral_basic_init(num_grid)
     F = reshape(F_lastInc,[9,cells(1),cells(2),cells3])
   end if restartRead
 
-  homogenization_F0 = reshape(F_lastInc, [3,3,product(cells(1:2))*cells3])                          ! set starting condition for homogenization_mechanical_response
   call utilities_updateCoords(reshape(F,shape(F_lastInc)))
   call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, &                                ! stress field, stress avg, global average of stiffness and (min+max)/2
                                       reshape(F,shape(F_lastInc)), &                                ! target F
@@ -347,8 +347,6 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_
                                    F_lastInc,reshape(F,[3,3,cells(1),cells(2),cells3]),Delta_t_old, &
                                    rotation_BC%rotate(F_aimDot,active=.true.))
     F_lastInc = reshape(F,[3,3,cells(1),cells(2),cells3])
-
-    homogenization_F0 = reshape(F,[3,3,product(cells(1:2))*cells3])
   end if
 
 !--------------------------------------------------------------------------------------------------

src/grid/grid_mech_spectral_polarization.f90

@@ -23,6 +23,7 @@ module grid_mechanical_spectral_polarization
   use math
   use rotations
   use spectral_utilities
+  use grid_mech_utilities
   use config
   use homogenization
   use discretization_grid
@@ -255,7 +256,6 @@ subroutine grid_mechanical_spectral_polarization_init(num_grid)
     F_tau_lastInc = 2.0_pREAL*F_lastInc
   end if restartRead
 
-  homogenization_F0 = reshape(F_lastInc, [3,3,product(cells(1:2))*cells3])                          ! set starting condition for homogenization_mechanical_response
   call utilities_updateCoords(reshape(F,shape(F_lastInc)))
   call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, &                                ! stress field, stress avg, global average of stiffness and (min+max)/2
                                       reshape(F,shape(F_lastInc)), &                                ! target F
@@ -391,7 +391,6 @@ subroutine grid_mechanical_spectral_polarization_forward(cutBack,guess,Delta_t,D
     F_lastInc     = reshape(F,    [3,3,cells(1),cells(2),cells3])
     F_tau_lastInc = reshape(F_tau,[3,3,cells(1),cells(2),cells3])
 
-    homogenization_F0 = reshape(F,[3,3,product(cells(1:2))*cells3])
   end if
 
 !--------------------------------------------------------------------------------------------------

src/grid/grid_mech_utilities.f90

@@ -0,0 +1,253 @@
+!--------------------------------------------------------------------------------------------------
+!> @author Martin Diehl, KU Leuven
+!> @brief Utilities used by the mech grid solver variants
+!--------------------------------------------------------------------------------------------------
+module grid_mech_utilities
+#include <petsc/finclude/petscsys.h>
+  use PETScSys
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  use MPI_f08
+#endif
+
+  use prec
+  use parallelization
+  use math
+  use rotations
+  use IO
+  use discretization_grid
+  use discretization
+  use spectral_utilities
+  use homogenization
+
+
+#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
+  implicit none(type,external)
+#else
+  implicit none
+#endif
+  private
+
+!--------------------------------------------------------------------------------------------------
+! derived types
+  type, public :: tBoundaryCondition                                                                !< set of parameters defining a boundary condition
+    real(pREAL), dimension(3,3)   :: values = 0.0_pREAL
+    logical,     dimension(3,3)   :: mask   = .true.
+    character(len=:), allocatable :: myType
+  end type tBoundaryCondition
+
+  type, public :: tSolutionParams
+    real(pREAL), dimension(3,3) :: stress_BC
+    logical, dimension(3,3)     :: stress_mask
+    type(tRotation)             :: rotation_BC
+    real(pREAL) :: Delta_t
+  end type tSolutionParams
+
+  public :: &
+    utilities_maskedCompliance, &
+    utilities_constitutiveResponse, &
+    utilities_calculateRate, &
+    utilities_forwardTensorField
+
+contains
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Calculate masked compliance tensor used to adjust F to fullfill stress BC.
+!--------------------------------------------------------------------------------------------------
+function utilities_maskedCompliance(rot_BC,mask_stress,C)
+
+  real(pREAL),                dimension(3,3,3,3) :: utilities_maskedCompliance                      !< masked compliance
+  real(pREAL),    intent(in), dimension(3,3,3,3) :: C                                               !< current average stiffness
+  type(tRotation), intent(in)                    :: rot_BC                                          !< rotation of load frame
+  logical,        intent(in), dimension(3,3)     :: mask_stress                                     !< mask of stress BC
+
+  integer :: i, j
+  logical, dimension(9)   :: mask_stressVector
+  logical, dimension(9,9) :: mask
+  real(pREAL), dimension(9,9) :: temp99_real
+  integer :: size_reduced = 0
+  real(pREAL),              dimension(:,:), allocatable ::  &
+    s_reduced, &                                                                                    !< reduced compliance matrix (depending on number of stress BC)
+    c_reduced, &                                                                                    !< reduced stiffness (depending on number of stress BC)
+    sTimesC                                                                                         !< temp variable to check inversion
+  logical :: errmatinv
+  character(len=pSTRLEN):: formatString
+
+
+  mask_stressVector = .not. reshape(transpose(mask_stress), [9])
+  size_reduced = count(mask_stressVector)
+  if (size_reduced > 0) then
+    temp99_real = math_3333to99(rot_BC%rotate(C))
+
+    do i = 1,9; do j = 1,9
+      mask(i,j) = mask_stressVector(i) .and. mask_stressVector(j)
+    end do; end do
+    c_reduced = reshape(pack(temp99_Real,mask),[size_reduced,size_reduced])
+
+    allocate(s_reduced,mold = c_reduced)
+    call math_invert(s_reduced, errmatinv, c_reduced)                                               ! invert reduced stiffness
+    if (any(IEEE_is_NaN(s_reduced))) errmatinv = .true.
+
+!--------------------------------------------------------------------------------------------------
+! check if inversion was successful
+    sTimesC = matmul(c_reduced,s_reduced)
+    errmatinv = errmatinv .or. any(dNeq(sTimesC,math_eye(size_reduced),1.0e-12_pREAL))
+    if (errmatinv) then
+      write(formatString, '(i2)') size_reduced
+      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), 'C (load) ', transpose(c_reduced)
+      print trim(formatString), 'S (load) ', transpose(s_reduced)
+      if (errmatinv) error stop 'matrix inversion error'
+    end if
+    temp99_real = reshape(unpack(reshape(s_reduced,[size_reduced**2]),reshape(mask,[81]),0.0_pREAL),[9,9])
+  else
+    temp99_real = 0.0_pREAL
+  end if
+
+  utilities_maskedCompliance = math_99to3333(temp99_Real)
+
+end function utilities_maskedCompliance
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Calculate constitutive response.
+!--------------------------------------------------------------------------------------------------
+subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
+                                          F,Delta_t,rotation_BC)
+
+  real(pREAL),    intent(out), dimension(3,3,3,3)                   :: C_volAvg, C_minmaxAvg        !< average stiffness
+  real(pREAL),    intent(out), dimension(3,3)                       :: P_av                         !< average PK stress
+  real(pREAL),    intent(out), dimension(3,3,cells(1),cells(2),cells3) :: P                         !< PK stress
+  real(pREAL),    intent(in),  dimension(3,3,cells(1),cells(2),cells3) :: F                         !< deformation gradient target
+  real(pREAL),    intent(in)                                        :: Delta_t                      !< loading time
+  type(tRotation), intent(in),  optional                            :: rotation_BC                  !< rotation of load frame
+
+
+  integer :: i
+  integer(MPI_INTEGER_KIND) :: err_MPI
+  real(pREAL), dimension(3,3,3,3) :: dPdF_max,      dPdF_min
+  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'(/,1x,a)', '... evaluating constitutive response ......................................'
+  flush(IO_STDOUT)
+
+  homogenization_F  = reshape(F,[3,3,product(cells(1:2))*cells3])                                   ! set materialpoint target F to estimated field
+
+  call homogenization_mechanical_response(Delta_t,1,product(cells(1:2))*cells3)                     ! calculate P field
+  if (.not. terminallyIll) &
+    call homogenization_thermal_response(Delta_t,1,product(cells(1:2))*cells3)
+  if (.not. terminallyIll) &
+    call homogenization_mechanical_response2(Delta_t,1,product(cells(1:2))*cells3)
+
+  P = reshape(homogenization_P, [3,3,cells(1),cells(2),cells3])
+  P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt
+  call MPI_Allreduce(MPI_IN_PLACE,P_av,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
+  if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
+  if (present(rotation_BC)) then
+    if (any(dNeq(rotation_BC%asQuaternion(), real([1.0, 0.0, 0.0, 0.0],pREAL)))) &
+      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
+    P_av = rotation_BC%rotate(P_av)
+  end if
+  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)
+
+  dPdF_max = 0.0_pREAL
+  dPdF_norm_max = 0.0_pREAL
+  dPdF_min = huge(1.0_pREAL)
+  dPdF_norm_min = huge(1.0_pREAL)
+  do i = 1, product(cells(1:2))*cells3
+    if (dPdF_norm_max < sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)) then
+      dPdF_max = homogenization_dPdF(1:3,1:3,1:3,1:3,i)
+      dPdF_norm_max = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)
+    end if
+    if (dPdF_norm_min > sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)) then
+      dPdF_min = homogenization_dPdF(1:3,1:3,1:3,1:3,i)
+      dPdF_norm_min = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)
+    end if
+  end do
+
+  valueAndRank = [dPdF_norm_max,real(worldrank,pREAL)]
+  call MPI_Allreduce(MPI_IN_PLACE,valueAndRank,1_MPI_INTEGER_KIND,MPI_2DOUBLE_PRECISION,MPI_MAXLOC,MPI_COMM_WORLD,err_MPI)
+  if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
+  call MPI_Bcast(dPdF_max,81_MPI_INTEGER_KIND,MPI_DOUBLE,int(valueAndRank(2),MPI_INTEGER_KIND),MPI_COMM_WORLD,err_MPI)
+  if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
+
+  valueAndRank = [dPdF_norm_min,real(worldrank,pREAL)]
+  call MPI_Allreduce(MPI_IN_PLACE,valueAndRank,1_MPI_INTEGER_KIND,MPI_2DOUBLE_PRECISION,MPI_MINLOC,MPI_COMM_WORLD,err_MPI)
+  if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
+  call MPI_Bcast(dPdF_min,81_MPI_INTEGER_KIND,MPI_DOUBLE,int(valueAndRank(2),MPI_INTEGER_KIND),MPI_COMM_WORLD,err_MPI)
+  if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
+
+  C_minmaxAvg = 0.5_pREAL*(dPdF_max + dPdF_min)
+
+  C_volAvg = sum(homogenization_dPdF,dim=5)
+  call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
+  if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
+  C_volAvg = C_volAvg * wgt
+
+
+end subroutine utilities_constitutiveResponse
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Calculate forward rate, either as local guess or as homogeneous add on.
+!--------------------------------------------------------------------------------------------------
+pure function utilities_calculateRate(heterogeneous,field0,field,dt,avRate)
+
+  real(pREAL), intent(in), dimension(3,3) :: &
+    avRate                                                                                          !< homogeneous addon
+  real(pREAL), intent(in) :: &
+    dt                                                                                              !< Delta_t between field0 and field
+  logical, intent(in) :: &
+    heterogeneous                                                                                   !< calculate field of rates
+  real(pREAL), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: &
+    field0, &                                                                                       !< data of previous step
+    field                                                                                           !< data of current step
+  real(pREAL),             dimension(3,3,cells(1),cells(2),cells3) :: &
+    utilities_calculateRate
+
+
+  utilities_calculateRate = merge((field-field0) / dt, &
+                                  spread(spread(spread(avRate,3,cells(1)),4,cells(2)),5,cells3), &
+                                  heterogeneous)
+
+end function utilities_calculateRate
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief forwards a field with a pointwise given rate, if aim is given,
+!> ensures that the average matches the aim
+!--------------------------------------------------------------------------------------------------
+function utilities_forwardTensorField(Delta_t,field_lastInc,rate,aim)
+
+  real(pREAL), intent(in) :: &
+    Delta_t                                                                                         !< Delta_t of current step
+  real(pREAL), intent(in),           dimension(3,3,cells(1),cells(2),cells3) :: &
+    field_lastInc, &                                                                                !< initial field
+    rate                                                                                            !< rate by which to forward
+  real(pREAL), intent(in), optional, dimension(3,3) :: &
+    aim                                                                                             !< average field value aim
+
+  real(pREAL),                       dimension(3,3,cells(1),cells(2),cells3) :: &
+    utilities_forwardTensorField
+  real(pREAL),                       dimension(3,3) :: fieldDiff                                    !< <a + adot*t> - aim
+  integer(MPI_INTEGER_KIND) :: err_MPI
+
+
+  utilities_forwardTensorField = field_lastInc + rate*Delta_t
+  if (present(aim)) then                                                                            !< correct to match average
+    fieldDiff = sum(sum(sum(utilities_forwardTensorField,dim=5),dim=4),dim=3)*wgt
+    call MPI_Allreduce(MPI_IN_PLACE,fieldDiff,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
+    if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
+    fieldDiff = fieldDiff - aim
+    utilities_forwardTensorField = utilities_forwardTensorField &
+                                 - spread(spread(spread(fieldDiff,3,cells(1)),4,cells(2)),5,cells3)
+  end if
+
+end function utilities_forwardTensorField
+
+end module grid_mech_utilities

src/grid/grid_thermal_spectral.f90

@@ -43,7 +43,6 @@ module grid_thermal_spectral
 
   type(tNumerics) :: num
 
-  type(tSolutionParams) :: params
 !--------------------------------------------------------------------------------------------------
 ! PETSc data
   SNES :: SNES_thermal
@@ -56,7 +55,7 @@ module grid_thermal_spectral
 ! reference diffusion tensor, mobility etc.
   integer                     :: totalIter = 0                                                      !< total iteration in current increment
   real(pREAL), dimension(3,3) :: K_ref
-  real(pREAL)                 :: mu_ref
+  real(pREAL)                 :: mu_ref, Delta_t_
 
   public :: &
     grid_thermal_spectral_init, &
@@ -186,8 +185,7 @@ end subroutine grid_thermal_spectral_init
 !--------------------------------------------------------------------------------------------------
 function grid_thermal_spectral_solution(Delta_t) result(solution)
 
-  real(pREAL), intent(in) :: &
+  real(pREAL), intent(in) ::  Delta_t                                                               !< increment in time for current solution
-    Delta_t                                                                                         !< increment in time for current solution
 
   type(tSolutionState) :: solution
   PetscInt  :: devNull
@@ -201,7 +199,7 @@ function grid_thermal_spectral_solution(Delta_t) result(solution)
 
 !--------------------------------------------------------------------------------------------------
 ! set module wide availabe data
-  params%Delta_t = Delta_t
+  Delta_t_ = Delta_t
 
   call SNESSolve(SNES_thermal,PETSC_NULL_VEC,T_PETSc,err_PETSc)
   CHKERRQ(err_PETSc)
@@ -227,7 +225,7 @@ function grid_thermal_spectral_solution(Delta_t) result(solution)
   T_stagInc = T
 
   call homogenization_thermal_setField(reshape(T,[product(cells(1:2))*cells3]), &
-                                       reshape(T-T_lastInc,[product(cells(1:2))*cells3])/params%Delta_t)
+                                       reshape(T-T_lastInc,[product(cells(1:2))*cells3])/Delta_t_)
 
   call DMDAVecRestoreArrayF90(DM_thermal,T_PETSc,T,err_PETSc)
   CHKERRQ(err_PETSc)
@@ -264,7 +262,7 @@ subroutine grid_thermal_spectral_forward(cutBack)
     T = T_lastInc
     T_stagInc = T_lastInc
   else
-    dotT_lastInc = (T - T_lastInc)/params%Delta_t
+    dotT_lastInc = (T - T_lastInc)/Delta_t_
     T_lastInc = T
     call updateReference()
   end if
@@ -336,13 +334,13 @@ subroutine formResidual(residual_subdomain,x_scal,r,dummy,err_PETSc)
     ce = 0
     do k = 1, cells3;  do j = 1, cells(2);  do i = 1,cells(1)
       ce = ce + 1
-      r(i,j,k) = params%Delta_t*(r(i,j,k) + homogenization_f_T(ce)) &
+      r(i,j,k) = Delta_t_*(r(i,j,k) + homogenization_f_T(ce)) &
                + homogenization_mu_T(ce) * (T_lastInc(i,j,k) - T(i,j,k)) &
                + mu_ref*T(i,j,k)
     end do; end do; end do
 
     r = T &
-      - utilities_GreenConvolution(r, K_ref, mu_ref, params%Delta_t)
+      - utilities_GreenConvolution(r, K_ref, mu_ref, Delta_t_)
   end associate
   err_PETSc = 0
 

src/grid/spectral_utilities.f90

@@ -75,19 +75,6 @@ module spectral_utilities
        termIll           = .false.
   end type tSolutionState
 
-  type, public :: tBoundaryCondition                                                                !< set of parameters defining a boundary condition
-    real(pREAL), dimension(3,3)   :: values = 0.0_pREAL
-    logical,     dimension(3,3)   :: mask   = .true.
-    character(len=:), allocatable :: myType
-  end type tBoundaryCondition
-
-  type, public :: tSolutionParams
-    real(pREAL), dimension(3,3) :: stress_BC
-    logical, dimension(3,3)     :: stress_mask
-    type(tRotation)             :: rotation_BC
-    real(pREAL) :: Delta_t
-  end type tSolutionParams
-
   type :: tNumerics
     integer :: &
       divergence_correction                                                                         !< scale divergence/curl calculation: [0: no correction, 1: size scaled to 1, 2: size scaled to Npoints]
@@ -121,10 +108,6 @@ module spectral_utilities
     utilities_curlRMS, &
     utilities_scalarGradient, &
     utilities_vectorDivergence, &
-    utilities_maskedCompliance, &
-    utilities_constitutiveResponse, &
-    utilities_calculateRate, &
-    utilities_forwardTensorField, &
     utilities_updateCoords
 
 contains
@@ -653,65 +636,6 @@ real(pREAL) function utilities_curlRMS(tensorField)
 end function utilities_curlRMS
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief Calculate masked compliance tensor used to adjust F to fullfill stress BC.
-!--------------------------------------------------------------------------------------------------
-function utilities_maskedCompliance(rot_BC,mask_stress,C)
-
-  real(pREAL),                dimension(3,3,3,3) :: utilities_maskedCompliance                      !< masked compliance
-  real(pREAL),    intent(in), dimension(3,3,3,3) :: C                                               !< current average stiffness
-  type(tRotation), intent(in)                    :: rot_BC                                          !< rotation of load frame
-  logical,        intent(in), dimension(3,3)     :: mask_stress                                     !< mask of stress BC
-
-  integer :: i, j
-  logical, dimension(9)   :: mask_stressVector
-  logical, dimension(9,9) :: mask
-  real(pREAL), dimension(9,9) :: temp99_real
-  integer :: size_reduced = 0
-  real(pREAL),              dimension(:,:), allocatable ::  &
-    s_reduced, &                                                                                    !< reduced compliance matrix (depending on number of stress BC)
-    c_reduced, &                                                                                    !< reduced stiffness (depending on number of stress BC)
-    sTimesC                                                                                         !< temp variable to check inversion
-  logical :: errmatinv
-  character(len=pSTRLEN):: formatString
-
-
-  mask_stressVector = .not. reshape(transpose(mask_stress), [9])
-  size_reduced = count(mask_stressVector)
-  if (size_reduced > 0) then
-    temp99_real = math_3333to99(rot_BC%rotate(C))
-
-    do i = 1,9; do j = 1,9
-      mask(i,j) = mask_stressVector(i) .and. mask_stressVector(j)
-    end do; end do
-    c_reduced = reshape(pack(temp99_Real,mask),[size_reduced,size_reduced])
-
-    allocate(s_reduced,mold = c_reduced)
-    call math_invert(s_reduced, errmatinv, c_reduced)                                               ! invert reduced stiffness
-    if (any(IEEE_is_NaN(s_reduced))) errmatinv = .true.
-
-!--------------------------------------------------------------------------------------------------
-! check if inversion was successful
-    sTimesC = matmul(c_reduced,s_reduced)
-    errmatinv = errmatinv .or. any(dNeq(sTimesC,math_eye(size_reduced),1.0e-12_pREAL))
-    if (errmatinv) then
-      write(formatString, '(i2)') size_reduced
-      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), 'C (load) ', transpose(c_reduced)
-      print trim(formatString), 'S (load) ', transpose(s_reduced)
-      if (errmatinv) error stop 'matrix inversion error'
-    end if
-    temp99_real = reshape(unpack(reshape(s_reduced,[size_reduced**2]),reshape(mask,[81]),0.0_pREAL),[9,9])
-  else
-    temp99_real = 0.0_pREAL
-  end if
-
-  utilities_maskedCompliance = math_99to3333(temp99_Real)
-
-end function utilities_maskedCompliance
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief Calculate gradient of scalar field.
 !--------------------------------------------------------------------------------------------------
@@ -755,147 +679,6 @@ function utilities_vectorDivergence(field) result(div)
 end function utilities_vectorDivergence
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief calculate constitutive response from homogenization_F0 to F during Delta_t
-!--------------------------------------------------------------------------------------------------
-subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
-                                          F,Delta_t,rotation_BC)
-
-  real(pREAL),    intent(out), dimension(3,3,3,3)                   :: C_volAvg, C_minmaxAvg        !< average stiffness
-  real(pREAL),    intent(out), dimension(3,3)                       :: P_av                         !< average PK stress
-  real(pREAL),    intent(out), dimension(3,3,cells(1),cells(2),cells3) :: P                         !< PK stress
-  real(pREAL),    intent(in),  dimension(3,3,cells(1),cells(2),cells3) :: F                         !< deformation gradient target
-  real(pREAL),    intent(in)                                        :: Delta_t                      !< loading time
-  type(tRotation), intent(in),  optional                            :: rotation_BC                  !< rotation of load frame
-
-
-  integer :: i
-  integer(MPI_INTEGER_KIND) :: err_MPI
-  real(pREAL), dimension(3,3,3,3) :: dPdF_max,      dPdF_min
-  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'(/,1x,a)', '... evaluating constitutive response ......................................'
-  flush(IO_STDOUT)
-
-  homogenization_F  = reshape(F,[3,3,product(cells(1:2))*cells3])                                   ! set materialpoint target F to estimated field
-
-  call homogenization_mechanical_response(Delta_t,1,product(cells(1:2))*cells3)                     ! calculate P field
-  if (.not. terminallyIll) &
-    call homogenization_thermal_response(Delta_t,1,product(cells(1:2))*cells3)
-  if (.not. terminallyIll) &
-    call homogenization_mechanical_response2(Delta_t,[1,1],[1,product(cells(1:2))*cells3])
-
-  P = reshape(homogenization_P, [3,3,cells(1),cells(2),cells3])
-  P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt
-  call MPI_Allreduce(MPI_IN_PLACE,P_av,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
-  if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
-  if (present(rotation_BC)) then
-    if (any(dNeq(rotation_BC%asQuaternion(), real([1.0, 0.0, 0.0, 0.0],pREAL)))) &
-      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
-    P_av = rotation_BC%rotate(P_av)
-  end if
-  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)
-
-  dPdF_max = 0.0_pREAL
-  dPdF_norm_max = 0.0_pREAL
-  dPdF_min = huge(1.0_pREAL)
-  dPdF_norm_min = huge(1.0_pREAL)
-  do i = 1, product(cells(1:2))*cells3
-    if (dPdF_norm_max < sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)) then
-      dPdF_max = homogenization_dPdF(1:3,1:3,1:3,1:3,i)
-      dPdF_norm_max = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)
-    end if
-    if (dPdF_norm_min > sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)) then
-      dPdF_min = homogenization_dPdF(1:3,1:3,1:3,1:3,i)
-      dPdF_norm_min = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2)
-    end if
-  end do
-
-  valueAndRank = [dPdF_norm_max,real(worldrank,pREAL)]
-  call MPI_Allreduce(MPI_IN_PLACE,valueAndRank,1_MPI_INTEGER_KIND,MPI_2DOUBLE_PRECISION,MPI_MAXLOC,MPI_COMM_WORLD,err_MPI)
-  if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
-  call MPI_Bcast(dPdF_max,81_MPI_INTEGER_KIND,MPI_DOUBLE,int(valueAndRank(2),MPI_INTEGER_KIND),MPI_COMM_WORLD,err_MPI)
-  if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
-
-  valueAndRank = [dPdF_norm_min,real(worldrank,pREAL)]
-  call MPI_Allreduce(MPI_IN_PLACE,valueAndRank,1_MPI_INTEGER_KIND,MPI_2DOUBLE_PRECISION,MPI_MINLOC,MPI_COMM_WORLD,err_MPI)
-  if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
-  call MPI_Bcast(dPdF_min,81_MPI_INTEGER_KIND,MPI_DOUBLE,int(valueAndRank(2),MPI_INTEGER_KIND),MPI_COMM_WORLD,err_MPI)
-  if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
-
-  C_minmaxAvg = 0.5_pREAL*(dPdF_max + dPdF_min)
-
-  C_volAvg = sum(homogenization_dPdF,dim=5)
-  call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
-  if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
-  C_volAvg = C_volAvg * wgt
-
-
-end subroutine utilities_constitutiveResponse
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief Calculate forward rate, either as local guess or as homogeneous add on.
-!--------------------------------------------------------------------------------------------------
-pure function utilities_calculateRate(heterogeneous,field0,field,dt,avRate)
-
-  real(pREAL), intent(in), dimension(3,3) :: &
-    avRate                                                                                          !< homogeneous addon
-  real(pREAL), intent(in) :: &
-    dt                                                                                              !< Delta_t between field0 and field
-  logical, intent(in) :: &
-    heterogeneous                                                                                   !< calculate field of rates
-  real(pREAL), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: &
-    field0, &                                                                                       !< data of previous step
-    field                                                                                           !< data of current step
-  real(pREAL),             dimension(3,3,cells(1),cells(2),cells3) :: &
-    utilities_calculateRate
-
-
-  utilities_calculateRate = merge((field-field0) / dt, &
-                                  spread(spread(spread(avRate,3,cells(1)),4,cells(2)),5,cells3), &
-                                  heterogeneous)
-
-end function utilities_calculateRate
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief forwards a field with a pointwise given rate, if aim is given,
-!> ensures that the average matches the aim
-!--------------------------------------------------------------------------------------------------
-function utilities_forwardTensorField(Delta_t,field_lastInc,rate,aim)
-
-  real(pREAL), intent(in) :: &
-    Delta_t                                                                                         !< Delta_t of current step
-  real(pREAL), intent(in),           dimension(3,3,cells(1),cells(2),cells3) :: &
-    field_lastInc, &                                                                                !< initial field
-    rate                                                                                            !< rate by which to forward
-  real(pREAL), intent(in), optional, dimension(3,3) :: &
-    aim                                                                                             !< average field value aim
-
-  real(pREAL),                       dimension(3,3,cells(1),cells(2),cells3) :: &
-    utilities_forwardTensorField
-  real(pREAL),                       dimension(3,3) :: fieldDiff                                    !< <a + adot*t> - aim
-  integer(MPI_INTEGER_KIND) :: err_MPI
-
-
-  utilities_forwardTensorField = field_lastInc + rate*Delta_t
-  if (present(aim)) then                                                                            !< correct to match average
-    fieldDiff = sum(sum(sum(utilities_forwardTensorField,dim=5),dim=4),dim=3)*wgt
-    call MPI_Allreduce(MPI_IN_PLACE,fieldDiff,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
-    if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
-    fieldDiff = fieldDiff - aim
-    utilities_forwardTensorField = utilities_forwardTensorField &
-                                 - spread(spread(spread(fieldDiff,3,cells(1)),4,cells(2)),5,cells3)
-  end if
-
-end function utilities_forwardTensorField
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief Calculate Filter for Fourier convolution.
 !> @details this is the full operator to calculate derivatives, i.e. 2 \pi i k for the

src/homogenization.f90

@@ -52,7 +52,6 @@ module homogenization
 !--------------------------------------------------------------------------------------------------
 ! General variables for the homogenization at a  material point
   real(pREAL),   dimension(:,:,:),     allocatable, public :: &
-    homogenization_F0, &                                                                            !< def grad of IP at start of FE increment
     homogenization_F                                                                                !< def grad of IP to be reached at end of FE increment
   real(pREAL),   dimension(:,:,:),     allocatable, public :: & !, protected :: &                   Issue with ifort
     homogenization_P                                                                                !< first P--K stress of IP
@@ -274,6 +273,7 @@ subroutine homogenization_thermal_response(Delta_t,cell_start,cell_end)
   real(pREAL), intent(in) :: Delta_t                                                                !< time increment
   integer, intent(in) :: &
     cell_start, cell_end
+
   integer :: &
     co, ce, ho
 
@@ -297,37 +297,33 @@ end subroutine homogenization_thermal_response
 !--------------------------------------------------------------------------------------------------
 !> @brief
 !--------------------------------------------------------------------------------------------------
-subroutine homogenization_mechanical_response2(Delta_t,FEsolving_execIP,FEsolving_execElem)
+subroutine homogenization_mechanical_response2(Delta_t,cell_start,cell_end)
 
   real(pREAL), intent(in) :: Delta_t                                                                !< time increment
-  integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP
+  integer, intent(in) :: &
+    cell_start, cell_end
+
   integer :: &
-    ip, &                                                                                           !< integration point number
-    el, &                                                                                           !< element number
     co, ce, ho
 
 
-  !$OMP PARALLEL DO PRIVATE(ho,ce)
+  !$OMP PARALLEL DO PRIVATE(ho)
-  elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2)
+  do ce = cell_start, cell_end
-    IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2)
-      ce = (el-1)*discretization_nIPs + ip
       ho = material_ID_homogenization(ce)
       do co = 1, homogenization_Nconstituents(ho)
-        call crystallite_orientations(co,ip,el)
+        call crystallite_orientations(co,ce)
       end do
       call mechanical_homogenize(Delta_t,ce)
-    end do IpLooping3
+  end do
-  end do elementLooping3
   !$OMP END PARALLEL DO
 
-
 end subroutine homogenization_mechanical_response2
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief writes homogenization results to HDF5 output file
 !--------------------------------------------------------------------------------------------------
-subroutine homogenization_result
+subroutine homogenization_result()
 
   integer :: ho
   character(len=:), allocatable :: group_base,group
@@ -362,7 +358,7 @@ end subroutine homogenization_result
 !> @brief Forward data after successful increment.
 ! ToDo: Any guessing for the current states possible?
 !--------------------------------------------------------------------------------------------------
-subroutine homogenization_forward
+subroutine homogenization_forward()
 
   integer :: ho
 

src/homogenization_mechanical.f90

@@ -77,8 +77,7 @@ module subroutine mechanical_init()
   call parseMechanical()
 
   allocate(homogenization_dPdF(3,3,3,3,discretization_Ncells), source=0.0_pREAL)
-  homogenization_F0 = spread(math_I3,3,discretization_Ncells)
+  homogenization_F = spread(math_I3,3,discretization_Ncells)
-  homogenization_F = homogenization_F0
   allocate(homogenization_P(3,3,discretization_Ncells),source=0.0_pREAL)
 
   if (any(mechanical_type == MECHANICAL_PASS_ID))      call pass_init()

src/mesh/FEM_utilities.f90

@@ -131,7 +131,7 @@ subroutine utilities_constitutiveResponse(Delta_t,P_av,forwardData)
 
   call homogenization_mechanical_response(Delta_t,1,mesh_maxNips*mesh_NcpElems)                     ! calculate P field
   if (.not. terminallyIll) &
-    call homogenization_mechanical_response2(Delta_t,[1,mesh_maxNips],[1,mesh_NcpElems])
+    call homogenization_mechanical_response2(Delta_t,1,mesh_maxNips*mesh_NcpElems)
   cutBack = .false.
 
   P_av = sum(homogenization_P,dim=3) * wgt

src/mesh/mesh_mech_FEM.f90

@@ -686,7 +686,6 @@ subroutine FEM_mechanical_forward(guess,Delta_t,Delta_t_prev,mechBC)
 ! forward last inc
   if (guess .and. .not. cutBack) then
     ForwardData = .True.
-    homogenization_F0 = homogenization_F
     call SNESGetDM(mechanical_snes,dm_local,err_PETSc)                                              !< retrieve mesh info from mechanical_snes into dm_local
     CHKERRQ(err_PETSc)
     call DMGetSection(dm_local,section,err_PETSc)

src/phase.f90

@@ -326,11 +326,8 @@ module phase
       real(pREAL) :: f
     end function phase_f_T
 
-    module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,ip,el)
+    module subroutine plastic_nonlocal_updateCompatibility(orientation,ce)
-      integer, intent(in) :: &
+      integer, intent(in) :: ce
-        ph, &
-        ip, &
-        el
         type(tRotationContainer), dimension(:), intent(in) :: orientation
     end subroutine plastic_nonlocal_updateCompatibility
 
@@ -387,7 +384,7 @@ contains
 !--------------------------------------------------------------------------------------------------
 !> @brief Initialize constitutive models for individual physics
 !--------------------------------------------------------------------------------------------------
-subroutine phase_init
+subroutine phase_init()
 
   integer :: &
     ph, ce, co, ma
@@ -544,25 +541,16 @@ subroutine crystallite_init()
   integer :: &
     ce, &
     co, &                                                                                           !< counter in integration point component loop
-    ip, &                                                                                           !< counter in integration point loop
-    el, &                                                                                           !< counter in element loop
     en, ph
-  type(tDict), pointer :: &
-    num_phase, &
-    phases
 
-  phases => config_material%get_dict('phase')
 
-  !$OMP PARALLEL DO PRIVATE(ce,ph,en)
+  !$OMP PARALLEL DO PRIVATE(ph,en)
-  do el = 1, discretization_Nelems
+  do ce = 1, size(material_ID_homogenization)
-    do ip = 1, discretization_nIPs
+    do co = 1,homogenization_Nconstituents(material_ID_homogenization(ce))
-      ce = (el-1)*discretization_nIPs + ip
+      ph = material_ID_phase(co,ce)
-      do co = 1,homogenization_Nconstituents(material_ID_homogenization(ce))
+      en = material_entry_phase(co,ce)
-        en = material_entry_phase(co,ce)
+      call crystallite_orientations(co,ce)
-        ph = material_ID_phase(co,ce)
+      call plastic_dependentState(ph,en)                                                          ! update dependent state variables to be consistent with basic states
-        call crystallite_orientations(co,ip,el)
-        call plastic_dependentState(ph,en)                                                          ! update dependent state variables to be consistent with basic states
-     end do
     end do
   end do
   !$OMP END PARALLEL DO
@@ -572,32 +560,30 @@ end subroutine crystallite_init
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief calculates orientations
+!> @brief Update orientations and, if needed, compatibility.
 !--------------------------------------------------------------------------------------------------
-subroutine crystallite_orientations(co,ip,el)
+subroutine crystallite_orientations(co,ce)
 
   integer, intent(in) :: &
-    co, &                                                                                           !< counter in integration point component loop
+    co, &
-    ip, &                                                                                           !< counter in integration point loop
+    ce
-    el                                                                                              !< counter in element loop
 
   integer :: ph, en
 
 
-  ph = material_ID_phase(co,(el-1)*discretization_nIPs + ip)
+  ph = material_ID_phase(co,ce)
-  en = material_entry_phase(co,(el-1)*discretization_nIPs + ip)
+  en = material_entry_phase(co,ce)
 
   call phase_O(ph)%data(en)%fromMatrix(transpose(math_rotationalPart(mechanical_F_e(ph,en))))
 
-  if (plasticState(material_ID_phase(1,(el-1)*discretization_nIPs + ip))%nonlocal) &
+  if (plasticState(material_ID_phase(1,ce))%nonlocal) call plastic_nonlocal_updateCompatibility(phase_O,ce)
-    call plastic_nonlocal_updateCompatibility(phase_O,material_ID_phase(1,(el-1)*discretization_nIPs + ip),ip,el)
 
 
 end subroutine crystallite_orientations
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Map 2nd order tensor to reference config
+!> @brief Map 2nd order tensor to reference configuration.
 !--------------------------------------------------------------------------------------------------
 function crystallite_push33ToRef(co,ce, tensor33)
 
@@ -621,15 +607,17 @@ end function crystallite_push33ToRef
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief determines whether a point is converged
+!> @brief Determine whether a point is converged.
 !--------------------------------------------------------------------------------------------------
 logical pure function converged(residuum,state,atol)
 
-  real(pREAL), intent(in), dimension(:) ::&
+  real(pREAL), intent(in), dimension(:) :: &
     residuum, state, atol
+
   real(pREAL) :: &
     rTol
 
+
   rTol = num%rTol_crystalliteState
 
   converged = all(abs(residuum) <= max(atol, rtol*abs(state)))

src/phase_mechanical_plastic_nonlocal.f90

@@ -44,8 +44,7 @@ submodule(phase:plastic) nonlocal
   ! BEGIN DEPRECATED
   integer, dimension(:,:,:), allocatable :: &
     iRhoU, &                                                                                        !< state indices for unblocked density
-    iV, &                                                                                           !< state indices for dislocation velocities
+    iV                                                                                              !< state indices for dislocation velocities
-    iD                                                                                              !< state indices for stable dipole height
   !END DEPRECATED
 
   real(pREAL), dimension(:,:,:,:,:,:), allocatable :: &
@@ -124,7 +123,9 @@ submodule(phase:plastic) nonlocal
   type :: tNonlocalDependentState
     real(pREAL), allocatable, dimension(:,:) :: &
       tau_pass, &
-      tau_back
+      tau_back, &
+      rho_forest, &
+      max_dipole_height
     real(pREAL), allocatable, dimension(:,:,:,:,:) :: &
       compatibility
   end type tNonlocalDependentState
@@ -146,7 +147,6 @@ submodule(phase:plastic) nonlocal
         rhoDip, &
           rho_dip_edg, &
           rho_dip_scr, &
-        rho_forest, &
       gamma, &
       v, &
         v_edg_pos, &
@@ -177,7 +177,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
   integer :: &
     ph, &
     Nmembers, &
-    sizeState, sizeDotState, sizeDependentState, sizeDeltaState, &
+    sizeState, sizeDotState, sizeDeltaState, &
     s1, s2, &
     s, t, l
   real(pREAL), dimension(:,:), allocatable :: &
@@ -389,11 +389,9 @@ module function plastic_nonlocal_init() result(myPlasticity)
                             'rhoSglScrewPosImmobile','rhoSglScrewNegImmobile', &
                             'rhoDipEdge            ','rhoDipScrew           ', &
                             'gamma                 ' ]) * prm%sum_N_sl                              !< "basic" microstructural state variables that are independent from other state variables
-    sizeDependentState = size([ 'rhoForest   ']) * prm%sum_N_sl                                     !< microstructural state variables that depend on other state variables
+    sizeState          = sizeDotState &
-    sizeState          = sizeDotState + sizeDependentState &
                        + size([ 'velocityEdgePos     ','velocityEdgeNeg     ', &
-                                'velocityScrewPos    ','velocityScrewNeg    ', &
+                                'velocityScrewPos    ','velocityScrewNeg    ']) * prm%sum_N_sl      !< other dependent state variables that are not updated by microstructure
-                                'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl     !< other dependent state variables that are not updated by microstructure
     sizeDeltaState            = sizeDotState
 
     call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState,0)     ! ToDo: state structure does not follow convention
@@ -477,15 +475,17 @@ module function plastic_nonlocal_init() result(myPlasticity)
     if (any(plasticState(ph)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl) < 0.0_pREAL)) &
       extmsg = trim(extmsg)//' atol_gamma'
 
-    stt%rho_forest => plasticState(ph)%state                 (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:Nmembers)
+    stt%v          => plasticState(ph)%state                 (11*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:Nmembers)
-    stt%v          => plasticState(ph)%state                 (12*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nmembers)
+    st0%v          => plasticState(ph)%state0                (11*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:Nmembers)
-        stt%v_edg_pos  => plasticState(ph)%state             (12*prm%sum_N_sl + 1:13*prm%sum_N_sl,1:Nmembers)
+        stt%v_edg_pos  => plasticState(ph)%state             (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:Nmembers)
-        stt%v_edg_neg  => plasticState(ph)%state             (13*prm%sum_N_sl + 1:14*prm%sum_N_sl,1:Nmembers)
+        stt%v_edg_neg  => plasticState(ph)%state             (12*prm%sum_N_sl + 1:13*prm%sum_N_sl,1:Nmembers)
-        stt%v_scr_pos  => plasticState(ph)%state             (14*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:Nmembers)
+        stt%v_scr_pos  => plasticState(ph)%state             (13*prm%sum_N_sl + 1:14*prm%sum_N_sl,1:Nmembers)
-        stt%v_scr_neg  => plasticState(ph)%state             (15*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nmembers)
+        stt%v_scr_neg  => plasticState(ph)%state             (14*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:Nmembers)
 
     allocate(dst%tau_pass(prm%sum_N_sl,Nmembers),source=0.0_pREAL)
     allocate(dst%tau_back(prm%sum_N_sl,Nmembers),source=0.0_pREAL)
+    allocate(dst%rho_forest(prm%sum_N_sl,Nmembers),source=0.0_pREAL)
+    allocate(dst%max_dipole_height(2*prm%sum_N_sl,Nmembers),source=0.0_pREAL)                       ! edge and screw
     allocate(dst%compatibility(2,maxval(param%sum_N_sl),maxval(param%sum_N_sl),nIPneighbors,Nmembers),source=0.0_pREAL)
     end associate
 
@@ -503,7 +503,6 @@ module function plastic_nonlocal_init() result(myPlasticity)
 ! BEGIN DEPRECATED----------------------------------------------------------------------------------
   allocate(iRhoU(maxval(param%sum_N_sl),4,phases%length), source=0)
   allocate(iV(maxval(param%sum_N_sl),4,phases%length),    source=0)
-  allocate(iD(maxval(param%sum_N_sl),2,phases%length),    source=0)
 
   do ph = 1, phases%length
 
@@ -518,20 +517,14 @@ module function plastic_nonlocal_init() result(myPlasticity)
         iRhoU(s,t,ph) = l
       end do
     end do
-    l = l + (4+2+1+1)*param(ph)%sum_N_sl ! immobile(4), dipole(2), shear, forest
+    l = l + (4+2+1)*param(ph)%sum_N_sl ! immobile(4), dipole(2), shear
     do t = 1,4
       do s = 1,param(ph)%sum_N_sl
         l = l + 1
         iV(s,t,ph) = l
       end do
     end do
-    do t = 1,2
+    if (iV(param(ph)%sum_N_sl,4,ph) /= plasticState(ph)%sizeState) &
-      do s = 1,param(ph)%sum_N_sl
-        l = l + 1
-        iD(s,t,ph) = l
-      end do
-    end do
-    if (iD(param(ph)%sum_N_sl,2,ph) /= plasticState(ph)%sizeState) &
       error stop 'state indices not properly set (nonlocal)'
   end do
 
@@ -602,7 +595,7 @@ module subroutine nonlocal_dependentState(ph, en)
   nu = elastic_nu(ph,en,prm%isotropic_bound)
   rho = getRho(ph,en)
 
-  stt%rho_forest(:,en) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) &
+  dst%rho_forest(:,en) = matmul(prm%forestProjection_Edge, sum(abs(rho(:,edg)),2)) &
                        + matmul(prm%forestProjection_Screw,sum(abs(rho(:,scr)),2))
 
 
@@ -612,7 +605,7 @@ module subroutine nonlocal_dependentState(ph, en)
     myInteractionMatrix = prm%h_sl_sl &
                         * spread((  1.0_pREAL - prm%f_F &
                                    + prm%f_F &
-                                   * log(0.35_pREAL * prm%b_sl * sqrt(max(stt%rho_forest(:,en),prm%rho_significant))) &
+                                   * log(0.35_pREAL * prm%b_sl * sqrt(max(dst%rho_forest(:,en),prm%rho_significant))) &
                                    / log(0.35_pREAL * prm%b_sl * 1e6_pREAL))**2,2,prm%sum_N_sl)
   else
     myInteractionMatrix = prm%h_sl_sl
@@ -861,14 +854,14 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
     deltaDUpper                                                                                     ! change in maximum stable dipole distance for edges and screws
 
 
-  associate(prm => param(ph),dst => dependentState(ph),del => deltaState(ph))
+  associate(prm => param(ph),dst => dependentState(ph),del => deltaState(ph), stt=>state(ph))
 
   mu = elastic_mu(ph,en,prm%isotropic_bound)
   nu = elastic_nu(ph,en,prm%isotropic_bound)
 
   !*** 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)
+  v = reshape(stt%v(:,en),[prm%sum_N_sl,4])
-  forall (s = 1:prm%sum_N_sl, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,ph),en)
+  dUpperOld = reshape(dst%max_dipole_height(:,en),[prm%sum_N_sl,2])
 
   rho =  getRho(ph,en)
   rhoDip = rho(:,dip)
@@ -915,7 +908,7 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
                                        / (dUpperOld(s,c) - prm%minDipoleHeight(s,c))
 
   forall (t=1:4) deltaRhoDipole2SingleStress(:,t) = -0.5_pREAL * deltaRhoDipole2SingleStress(:,(t-1)/2+9)
-  forall (s = 1:prm%sum_N_sl, c = 1:2) plasticState(ph)%state(iD(s,c,ph),en) = dUpper(s,c)
+  dst%max_dipole_height(:,en) = pack(dUpper,.true.)
 
   plasticState(ph)%deltaState(:,en) = 0.0_pREAL
   del%rho(:,en) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*prm%sum_N_sl])
@@ -975,7 +968,8 @@ module subroutine nonlocal_dotState(Mp,timestep, &
     return
   end if
 
-  associate(prm => param(ph), dst => dependentState(ph), dot => dotState(ph), stt => state(ph))
+  associate(prm => param(ph), dst => dependentState(ph), dot => dotState(ph), &
+            stt => state(ph), st0 => state0(ph))
 
   mu = elastic_mu(ph,en,prm%isotropic_bound)
   nu = elastic_nu(ph,en,prm%isotropic_bound)
@@ -990,11 +984,10 @@ module subroutine nonlocal_dotState(Mp,timestep, &
   rho0 = getRho0(ph,en)
   my_rhoSgl0 = rho0(:,sgl)
 
-  forall (s = 1:prm%sum_N_sl, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en)
+  v = reshape(stt%v(:,en),[prm%sum_N_sl,4])
   dot_gamma = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
 
 
-
   ! limits for stable dipole height
   do s = 1,prm%sum_N_sl
     tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s)) + dst%tau_back(s,en)
@@ -1018,20 +1011,20 @@ module subroutine nonlocal_dotState(Mp,timestep, &
   isBCC: if (phase_lattice(ph) == 'cI') then
     forall (s = 1:prm%sum_N_sl, sum(abs(v(s,1:4))) > 0.0_pREAL)
       rhoDotMultiplication(s,1:2) = sum(abs(dot_gamma(s,3:4))) / prm%b_sl(s) &                      ! assuming double-cross-slip of screws to be decisive for multiplication
-                                  * sqrt(stt%rho_forest(s,en)) / prm%i_sl(s) ! &                    ! mean free path
+                                  * sqrt(dst%rho_forest(s,en)) / prm%i_sl(s) ! &                    ! mean free path
                                   ! * 2.0_pREAL * sum(abs(v(s,3:4))) / sum(abs(v(s,1:4)))           ! ratio of screw to overall velocity determines edge generation
       rhoDotMultiplication(s,3:4) = sum(abs(dot_gamma(s,3:4))) /prm%b_sl(s) &                       ! assuming double-cross-slip of screws to be decisive for multiplication
-                                  * sqrt(stt%rho_forest(s,en)) / prm%i_sl(s) ! &                    ! mean free path
+                                  * sqrt(dst%rho_forest(s,en)) / prm%i_sl(s) ! &                    ! mean free path
                                   ! * 2.0_pREAL * sum(abs(v(s,1:2))) / sum(abs(v(s,1:4)))           ! ratio of edge to overall velocity determines screw generation
     endforall
 
   else isBCC
     rhoDotMultiplication(:,1:4) = spread( &
           (sum(abs(dot_gamma(:,1:2)),2) * prm%f_ed_mult + sum(abs(dot_gamma(:,3:4)),2)) &
-        * sqrt(stt%rho_forest(:,en)) / prm%i_sl / prm%b_sl, 2, 4)                                   ! eq. 3.26
+        * sqrt(dst%rho_forest(:,en)) / prm%i_sl / prm%b_sl, 2, 4)                                   ! eq. 3.26
   end if isBCC
 
-  forall (s = 1:prm%sum_N_sl, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en)
+  v0 = reshape(st0%v(:,en),[prm%sum_N_sl,4])
 
 
   !****************************************************************************
@@ -1074,7 +1067,7 @@ module subroutine nonlocal_dotState(Mp,timestep, &
   if (phase_lattice(ph) == 'cF') &
     forall (s = 1:prm%sum_N_sl, prm%colinearSystem(s) > 0) &
       rhoDotAthermalAnnihilation(prm%colinearSystem(s),1:2) = - rhoDotAthermalAnnihilation(s,10) &
-        * 0.25_pREAL * sqrt(stt%rho_forest(s,en)) * (dUpper(s,2) + dLower(s,2)) * prm%f_ed
+        * 0.25_pREAL * sqrt(dst%rho_forest(s,en)) * (dUpper(s,2) + dLower(s,2)) * prm%f_ed
 
 
   ! thermally activated annihilation of edge dipoles by climb
@@ -1171,7 +1164,8 @@ function rhoDotFlux(timestep,ph,en)
 
   associate(prm => param(ph), &
             dst => dependentState(ph), &
-            stt => state(ph))
+            stt => state(ph), &
+            st0 => state0(ph))
   ns = prm%sum_N_sl
 
   dot_gamma = 0.0_pREAL
@@ -1181,10 +1175,10 @@ function rhoDotFlux(timestep,ph,en)
   rho0 = getRho0(ph,en)
   my_rhoSgl0 = rho0(:,sgl)
 
-  forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en)                         !ToDo: MD: I think we should use state0 here
+  v = reshape(stt%v(:,en),[prm%sum_N_sl,4])                                                         !ToDo: MD: I think we should use state0 here
   dot_gamma = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
 
-  forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en)
+  v0 = reshape(st0%v(:,en),[prm%sum_N_sl,4])
 
   !****************************************************************************
   !*** calculate dislocation fluxes (only for nonlocal plasticity)
@@ -1331,18 +1325,19 @@ end function rhoDotFlux
 ! plane normals and signed cosine of the angle between the slip directions. Only the largest values
 ! that sum up to a total of 1 are considered, all others are set to zero.
 !--------------------------------------------------------------------------------------------------
-module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,ip,el)
+module subroutine plastic_nonlocal_updateCompatibility(orientation,ce)
 
   type(tRotationContainer), dimension(:), intent(in) :: &
     orientation                                                                                     ! crystal orientation
   integer, intent(in) :: &
-    ph, &
+    ce
-    ip, &
-    el
 
   integer :: &
     n, &                                                                                            ! neighbor index
+    ph, &
     en, &
+    ip, &
+    el, &
     neighbor_e, &                                                                                   ! element index of my neighbor
     neighbor_i, &                                                                                   ! integration point index of my neighbor
     neighbor_me, &
@@ -1350,17 +1345,21 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,ip,el)
     ns, &                                                                                           ! number of active slip systems
     s1, &                                                                                           ! slip system index (en)
     s2                                                                                              ! slip system index (my neighbor)
-  real(pREAL), dimension(2,param(ph)%sum_N_sl,param(ph)%sum_N_sl,nIPneighbors) :: &
+  real(pREAL), dimension(2,param(material_ID_phase(1,ce))%sum_N_sl,param(material_ID_phase(1,ce))%sum_N_sl,nIPneighbors) :: &
     my_compatibility                                                                                ! my_compatibility for current element and ip
   real(pREAL) :: &
     my_compatibilitySum, &
     thresholdValue, &
     nThresholdValues
-  logical, dimension(param(ph)%sum_N_sl) :: &
+  logical, dimension(param(material_ID_phase(1,ce))%sum_N_sl) :: &
     belowThreshold
   type(tRotation) :: mis
 
 
+  ph = material_ID_phase(1,ce)
+  el = (ce-1)/discretization_nIPs + 1
+  ip = modulo(ce-1,discretization_nIPs) + 1
+
   associate(prm => param(ph))
     ns = prm%sum_N_sl
 
@@ -1486,7 +1485,7 @@ module subroutine plastic_nonlocal_result(ph,group)
           call result_writeDataset(stt%rho_dip_scr,group,trim(prm%output(ou)), &
                                    'screw dipole density','1/m²', prm%systems_sl)
         case('rho_f')
-          call result_writeDataset(stt%rho_forest,group,trim(prm%output(ou)), &
+          call result_writeDataset(dst%rho_forest,group,trim(prm%output(ou)), &
                                    'forest density','1/m²', prm%systems_sl)
         case('v_ed_pos')
           call result_writeDataset(stt%v_edg_pos,group,trim(prm%output(ou)), &