Merge branch 'development' into 42-new-coding-style-for-homogenization

.gitlab-ci.yml

@@ -3,8 +3,8 @@ stages:
   - prepareAll
   - preprocessing
   - postprocessing
-  - compileSpectralIntel
+  - compilePETScIntel
-  - compileSpectralGNU
+  - compilePETScGNU
   - prepareSpectral
   - spectral
   - compileMarc2017
@@ -186,8 +186,8 @@ Post_ParaviewRelated:
     - release
 
 ###################################################################################################
-Compile_Intel:
+Compile_Spectral_Intel:
-  stage: compileSpectralIntel
+  stage: compilePETScIntel
   script:
     - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel
     - SpectralAll_compile/test.py
@@ -195,9 +195,18 @@ Compile_Intel:
     - master
     - release
 
+Compile_FEM_Intel:
+  stage: compilePETScIntel
+  script:
+    - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel
+    - FEM_compile/test.py
+  except:
+    - master
+    - release
+
 ###################################################################################################
-Compile_GNU:
+Compile_Spectral_GNU:
-  stage: compileSpectralGNU
+  stage: compilePETScGNU
   script:
     - module load $GNUCompiler $MPICH_GNU $PETSc_MPICH_GNU
     - SpectralAll_compile/test.py
@@ -205,6 +214,15 @@ Compile_GNU:
     - master
     - release
 
+Compile_FEM_GNU:
+  stage: compilePETScGNU
+  script:
+    - module load $GNUCompiler $MPICH_GNU $PETSc_MPICH_GNU
+    - FEM_compile/test.py
+  except:
+    - master
+    - release
+
 ###################################################################################################
 Compile_Intel_Prepare:
   stage: prepareSpectral

PRIVATE

@@ -1 +1 @@
-Subproject commit dfd67ea44ba88ee1e0a33266a3986c64137908cf
+Subproject commit a764ade044735df35fac93a5204446291ee29abc

VERSION

@@ -1 +1 @@
-v2.0.2-401-ga3b472a7
+v2.0.2-474-g38fd517c

src/CMakeLists.txt

@@ -17,13 +17,7 @@ list(APPEND OBJECTFILES $<TARGET_OBJECTS:SYSTEM_ROUTINES>)
 add_library(PREC OBJECT "prec.f90")
 list(APPEND OBJECTFILES $<TARGET_OBJECTS:PREC>)
 
-if (PROJECT_NAME STREQUAL "DAMASK_spectral")
+add_library(DAMASK_INTERFACE OBJECT "DAMASK_interface.f90")
-      add_library(DAMASK_INTERFACE OBJECT "spectral_interface.f90")
-elseif (PROJECT_NAME STREQUAL "DAMASK_FEM")
-      add_library(DAMASK_INTERFACE OBJECT "FEM_interface.f90")
-else ()
-    message (FATAL_ERROR "Build target (PROJECT_NAME) is not defined")
-endif()
 add_dependencies(DAMASK_INTERFACE PREC SYSTEM_ROUTINES)
 list(APPEND OBJECTFILES $<TARGET_OBJECTS:DAMASK_INTERFACE>)
 
@@ -57,7 +51,7 @@ if (PROJECT_NAME STREQUAL "DAMASK_spectral")
       add_dependencies(MESH DAMASK_MATH)
       list(APPEND OBJECTFILES $<TARGET_OBJECTS:MESH>)
 elseif (PROJECT_NAME STREQUAL "DAMASK_FEM")
-      add_library(FEZoo OBJECT "FEZoo.f90")
+      add_library(FEZoo OBJECT "FEM_zoo.f90")
       add_dependencies(FEZoo DAMASK_MATH)
       list(APPEND OBJECTFILES $<TARGET_OBJECTS:FEZoo>)
       add_library(MESH  OBJECT "meshFEM.f90")
@@ -175,25 +169,24 @@ if (PROJECT_NAME STREQUAL "DAMASK_spectral")
       "spectral_mech_Basic.f90")
       add_dependencies(SPECTRAL_SOLVER SPECTRAL_UTILITIES)
       list(APPEND OBJECTFILES $<TARGET_OBJECTS:SPECTRAL_SOLVER>)
+
       if(NOT CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY")
         add_executable(DAMASK_spectral "DAMASK_spectral.f90" ${OBJECTFILES})
       else()
         add_library(DAMASK_spectral OBJECT "DAMASK_spectral.f90")
       endif()
+
       add_dependencies(DAMASK_spectral SPECTRAL_SOLVER)
 elseif (PROJECT_NAME STREQUAL "DAMASK_FEM")
       add_library(FEM_UTILITIES OBJECT "FEM_utilities.f90")
       add_dependencies(FEM_UTILITIES DAMASK_CPFE)
+      list(APPEND OBJECTFILES $<TARGET_OBJECTS:FEM_UTILITIES>)
 
       add_library(FEM_SOLVER OBJECT
-      "FEM_hydrogenflux.f90"
-      "FEM_porosity.f90"
-      "FEM_vacancyflux.f90"
-      "FEM_damage.f90"
-      "FEM_thermal.f90"
       "FEM_mech.f90")
       add_dependencies(FEM_SOLVER FEM_UTILITIES)
+      list(APPEND OBJECTFILES $<TARGET_OBJECTS:FEM_SOLVER>)
 
-      add_executable(DAMASK_FEM "DAMASK_FEM_driver.f90")
+      add_executable(DAMASK_FEM "DAMASK_FEM.f90" ${OBJECTFILES})
       add_dependencies(DAMASK_FEM FEM_SOLVER)
 endif()

src/CPFEM2.f90

@@ -50,8 +50,8 @@ subroutine CPFEM_initAll(el,ip)
    IO_init
  use DAMASK_interface
 #ifdef FEM
- use FEZoo, only: &
+ use FEM_Zoo, only: &
-   FEZoo_init
+   FEM_Zoo_init
 #endif
 
  implicit none
@@ -62,7 +62,7 @@ subroutine CPFEM_initAll(el,ip)
  call prec_init
  call IO_init
 #ifdef FEM
- call FEZoo_init
+ call FEM_Zoo_init
 #endif
  call numerics_init
  call debug_init
@@ -196,7 +196,7 @@ end subroutine CPFEM_init
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief perform initialization at first call, update variables and call the actual material model
+!> @brief forwards data after successful increment
 !--------------------------------------------------------------------------------------------------
 subroutine CPFEM_age()
  use prec, only: &
@@ -212,16 +212,6 @@ subroutine CPFEM_age()
    debug_levelSelective
  use FEsolving, only: &
    restartWrite
- use math, only: &
-   math_identity2nd, &
-   math_mul33x33, &
-   math_det33, &
-   math_transpose33, &
-   math_I3, &
-   math_Mandel3333to66, &
-   math_Mandel66to3333, &
-   math_Mandel33to6, &
-   math_Mandel6to33
  use material, only: &
    plasticState, &
    sourceState, &

src/DAMASK_FEM.f90

@@ -0,0 +1,654 @@
+!--------------------------------------------------------------------------------------------------
+!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
+!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
+!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
+!> @brief Driver controlling inner and outer load case looping of the FEM solver
+!> @details doing cutbacking, forwarding in case of restart, reporting statistics, writing
+!> results
+!--------------------------------------------------------------------------------------------------
+program DAMASK_FEM 
+ use, intrinsic :: &
+   iso_fortran_env                                                                                  ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment)
+ use prec, only: &
+   pInt, &
+   pReal, &
+   tol_math_check
+ use DAMASK_interface, only: &
+   DAMASK_interface_init, &
+   loadCaseFile, &
+   getSolverJobName
+ use IO, only: &
+   IO_read, &
+   IO_isBlank, &
+   IO_open_file, &
+   IO_stringPos, &
+   IO_stringValue, &
+   IO_floatValue, &
+   IO_intValue, &
+   IO_error, &
+   IO_lc, &
+   IO_intOut, &
+   IO_warning, &
+   IO_timeStamp, &
+   IO_EOF
+ use debug, only: &
+   debug_level, &
+   debug_spectral, &
+   debug_levelBasic
+ use math                                                                                           ! need to include the whole module for FFTW
+ use CPFEM2, only: &
+   CPFEM_initAll
+ use FEsolving, only: &
+   restartWrite, &
+   restartInc
+ use numerics, only: &
+   maxCutBack, &
+   stagItMax, &
+   worldrank
+ use mesh, only: &
+   mesh_Nboundaries, &
+   mesh_boundaries, &   
+   geomMesh
+ use FEM_Utilities, only: &
+   utilities_init, &
+   tSolutionState, &
+   tLoadCase, &
+   cutBack, &
+   maxFields, &
+   nActiveFields, &
+   FIELD_MECH_ID, &
+   FIELD_THERMAL_ID, &
+   FIELD_DAMAGE_ID, &
+   FIELD_SOLUTE_ID, &
+   FIELD_MGTWIN_ID, &
+   COMPONENT_MECH_X_ID, &
+   COMPONENT_MECH_Y_ID, &
+   COMPONENT_MECH_Z_ID, &
+   COMPONENT_THERMAL_T_ID, &
+   COMPONENT_DAMAGE_PHI_ID, &
+   COMPONENT_SOLUTE_CV_ID, &
+   COMPONENT_SOLUTE_CVPOT_ID, &
+   COMPONENT_SOLUTE_CH_ID, &
+   COMPONENT_SOLUTE_CHPOT_ID, &
+   COMPONENT_SOLUTE_CVaH_ID, &
+   COMPONENT_SOLUTE_CVaHPOT_ID, &
+   COMPONENT_MGTWIN_PHI_ID, &
+   FIELD_MECH_label, &
+   FIELD_THERMAL_label, &
+   FIELD_DAMAGE_label, &
+   FIELD_SOLUTE_label, &
+   FIELD_MGTWIN_label
+ use FEM_mech
+ 
+ implicit none
+#include <petsc/finclude/petsc.h>
+
+!--------------------------------------------------------------------------------------------------
+! variables related to information from load case and geom file
+ integer(pInt), parameter  :: FILEUNIT           = 234_pInt                                         !< file unit, DAMASK IO does not support newunit feature
+ integer(pInt), allocatable, dimension(:) :: chunkPos                                               ! this is longer than needed for geometry parsing
+ 
+ integer(pInt) :: &
+   N_def = 0_pInt                                                                                   !< # of rate of deformation specifiers found in load case file
+ character(len=65536) :: &
+   line
+
+!--------------------------------------------------------------------------------------------------
+! loop variables, convergence etc.
+
+ integer(pInt), parameter :: &
+   subStepFactor = 2_pInt                                                                           !< for each substep, divide the last time increment by 2.0
+ real(pReal) :: &
+   time = 0.0_pReal, &                                                                              !< elapsed time
+   time0 = 0.0_pReal, &                                                                             !< begin of interval
+   timeinc = 0.0_pReal, &                                                                           !< current time interval
+   timeIncOld = 0.0_pReal, &                                                                        !< previous time interval
+   remainingLoadCaseTime = 0.0_pReal                                                                !< remaining time of current load case
+ logical :: &
+   guess                                                                                            !< guess along former trajectory
+ integer(pInt) :: &
+   i, &
+   errorID, &
+   cutBackLevel = 0_pInt, &                                                                         !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$
+   stepFraction = 0_pInt                                                                            !< fraction of current time interval
+ integer(pInt) :: &
+   currentLoadcase = 0_pInt, &                                                                      !< current load case
+   currentFace = 0_pInt, &
+   inc, &                                                                                           !< current increment in current load case
+   totalIncsCounter = 0_pInt, &                                                                     !< total No. of increments
+   convergedCounter = 0_pInt, &                                                                     !< No. of converged increments
+   notConvergedCounter = 0_pInt, &                                                                  !< No. of non-converged increments
+   statUnit = 0_pInt, &                                                                             !< file unit for statistics output
+   lastRestartWritten = 0_pInt                                                                      !< total increment No. at which last restart information was written
+ integer(pInt) :: &
+   stagIter, &
+   component
+ logical :: &
+   stagIterate  
+ character(len=6)  :: loadcase_string
+ character(len=1024)  :: incInfo                                                                    !< string parsed to solution with information about current load case
+ type(tLoadCase), allocatable, dimension(:) :: loadCases                                            !< array of all load cases
+ type(tSolutionState), allocatable, dimension(:) :: solres
+ PetscInt :: faceSet, currentFaceSet
+ PetscInt :: field, dimPlex
+ PetscErrorCode :: ierr
+
+ external :: &
+   MPI_abort, &
+   DMGetDimension, &
+   DMGetLabelSize, &
+   DMGetLabelIdIS, &
+   ISDestroy, &
+   quit
+!--------------------------------------------------------------------------------------------------
+! init DAMASK (all modules)
+ call CPFEM_initAll(el = 1_pInt, ip = 1_pInt)
+   write(6,'(/,a)')   ' <<<+-  DAMASK_FEM init  -+>>>'
+   write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
+#include "compilation_info.f90"
+ 
+! reading basic information from load case file and allocate data structure containing load cases
+ call DMGetDimension(geomMesh,dimPlex,ierr)! CHKERRQ(ierr)                                            !< dimension of mesh (2D or 3D)
+ nActiveFields = 1
+ allocate(solres(nActiveFields))
+
+!--------------------------------------------------------------------------------------------------
+! reading basic information from load case file and allocate data structure containing load cases
+ call IO_open_file(FILEUNIT,trim(loadCaseFile))
+ rewind(FILEUNIT)
+ do
+   line = IO_read(FILEUNIT)
+   if (trim(line) == IO_EOF) exit
+   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
+   chunkPos = IO_stringPos(line)
+   do i = 1_pInt, chunkPos(1)                                                                      ! reading compulsory parameters for loadcase
+     select case (IO_lc(IO_stringValue(line,chunkPos,i)))
+       case('$loadcase')
+         N_def = N_def + 1_pInt
+     end select
+   enddo                                                                                            ! count all identifiers to allocate memory and do sanity check
+ enddo
+
+ allocate (loadCases(N_def))         
+
+ do i = 1, size(loadCases)
+   allocate(loadCases(i)%fieldBC(nActiveFields))
+   field = 1
+   loadCases(i)%fieldBC(field)%ID = FIELD_MECH_ID
+ enddo
+
+ do i = 1, size(loadCases)
+   do field = 1, nActiveFields
+     select case (loadCases(i)%fieldBC(field)%ID)
+       case(FIELD_MECH_ID)
+         loadCases(i)%fieldBC(field)%nComponents = dimPlex                                        !< X, Y (, Z) displacements
+         allocate(loadCases(i)%fieldBC(field)%componentBC(loadCases(i)%fieldBC(field)%nComponents))
+     end select
+     do component = 1, loadCases(i)%fieldBC(field)%nComponents
+       allocate(loadCases(i)%fieldBC(field)%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pReal)
+       allocate(loadCases(i)%fieldBC(field)%componentBC(component)%Mask (mesh_Nboundaries), source = .false.)
+     enddo
+   enddo       
+ enddo   
+
+!--------------------------------------------------------------------------------------------------
+! reading the load case and assign values to the allocated data structure
+ rewind(FILEUNIT)
+ do
+   line = IO_read(FILEUNIT)
+   if (trim(line) == IO_EOF) exit
+   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
+   chunkPos = IO_stringPos(line)
+   do i = 1_pInt, chunkPos(1)
+     select case (IO_lc(IO_stringValue(line,chunkPos,i)))
+!--------------------------------------------------------------------------------------------------
+! loadcase information
+       case('$loadcase')
+         currentLoadCase = IO_intValue(line,chunkPos,i+1_pInt)
+       case('face')
+         currentFace = IO_intValue(line,chunkPos,i+1_pInt)
+         currentFaceSet = -1_pInt
+         do faceSet = 1, mesh_Nboundaries
+           if (mesh_boundaries(faceSet) == currentFace) currentFaceSet = faceSet
+         enddo
+         if (currentFaceSet < 0_pInt) call IO_error(error_ID = errorID, ext_msg = 'invalid BC')  
+       case('t','time','delta')                                                                     ! increment time
+         loadCases(currentLoadCase)%time = IO_floatValue(line,chunkPos,i+1_pInt)
+       case('n','incs','increments','steps')                                                        ! number of increments
+         loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1_pInt)
+       case('logincs','logincrements','logsteps')                                                   ! number of increments (switch to log time scaling)
+         loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1_pInt)
+         loadCases(currentLoadCase)%logscale = 1_pInt
+       case('freq','frequency','outputfreq')                                                        ! frequency of result writings
+         loadCases(currentLoadCase)%outputfrequency = IO_intValue(line,chunkPos,i+1_pInt)                
+       case('r','restart','restartwrite')                                                           ! frequency of writing restart information
+         loadCases(currentLoadCase)%restartfrequency = &
+               max(0_pInt,IO_intValue(line,chunkPos,i+1_pInt))                
+       case('guessreset','dropguessing')
+         loadCases(currentLoadCase)%followFormerTrajectory = .false.                                ! do not continue to predict deformation along former trajectory
+
+!--------------------------------------------------------------------------------------------------
+! boundary condition information
+       case('x')                                                                                    ! X displacement field  
+         do field = 1, nActiveFields
+           if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then
+             do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
+               if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MECH_X_ID) then
+                 loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
+                     .true.
+                 loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
+                     IO_floatValue(line,chunkPos,i+1_pInt)
+               endif
+             enddo
+           endif
+         enddo  
+       case('y')                                                                                    ! Y displacement field                                                                            
+         do field = 1, nActiveFields
+           if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then
+             do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
+               if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MECH_Y_ID) then
+                 loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
+                     .true.
+                 loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
+                     IO_floatValue(line,chunkPos,i+1_pInt)
+               endif
+             enddo
+           endif
+         enddo  
+       case('z')                                                                                    ! Z displacement field                                                                       
+         do field = 1, nActiveFields
+           if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then
+             do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
+               if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MECH_Z_ID) then
+                 loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
+                     .true.
+                 loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
+                     IO_floatValue(line,chunkPos,i+1_pInt)
+               endif
+             enddo
+           endif
+         enddo  
+       case('temp','temperature')                                                                   ! thermal field                                                                      
+         do field = 1, nActiveFields
+           if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_THERMAL_ID) then
+             do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
+               if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_THERMAL_T_ID) then
+                 loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
+                     .true.
+                 loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
+                     IO_floatValue(line,chunkPos,i+1_pInt)
+               endif
+             enddo
+           endif
+         enddo  
+       case('mgtwin')                                                                               ! mgtwin field                                                                      
+         do field = 1, nActiveFields
+           if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MGTWIN_ID) then
+             do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
+               if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_MGTWIN_PHI_ID) then
+                 loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
+                     .true.
+                 loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
+                     IO_floatValue(line,chunkPos,i+1_pInt)
+               endif
+             enddo
+           endif
+         enddo 
+       case('damage')                                                                             
+         do field = 1, nActiveFields
+           if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_DAMAGE_ID) then
+             do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
+               if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_DAMAGE_PHI_ID) then
+                 loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
+                     .true.
+                 loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
+                     IO_floatValue(line,chunkPos,i+1_pInt)
+               endif
+             enddo
+           endif
+         enddo  
+       case('cv')                                                                             
+         do field = 1, nActiveFields
+           if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
+             do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
+               if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CV_ID) then
+                 loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
+                     .true.
+                 loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
+                     IO_floatValue(line,chunkPos,i+1_pInt)
+               endif
+             enddo
+           endif
+         enddo  
+       case('cvpot')                                                                             
+         do field = 1, nActiveFields
+           if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
+             do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
+               if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CVPOT_ID) then
+                 loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
+                     .true.
+                 loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
+                     IO_floatValue(line,chunkPos,i+1_pInt)
+               endif
+             enddo
+           endif
+         enddo  
+       case('ch')                                                                             
+         do field = 1, nActiveFields
+           if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
+             do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
+               if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CH_ID) then
+                 loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
+                     .true.
+                 loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
+                     IO_floatValue(line,chunkPos,i+1_pInt)
+               endif
+             enddo
+           endif
+         enddo  
+       case('chpot')                                                                             
+         do field = 1, nActiveFields
+           if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
+             do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
+               if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CHPOT_ID) then
+                 loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
+                     .true.
+                 loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
+                     IO_floatValue(line,chunkPos,i+1_pInt)
+               endif
+             enddo
+           endif
+         enddo  
+       case('cvah')                                                                             
+         do field = 1, nActiveFields
+           if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
+             do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
+               if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CVaH_ID) then
+                 loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
+                     .true.
+                 loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
+                     IO_floatValue(line,chunkPos,i+1_pInt)
+               endif
+             enddo
+           endif
+         enddo  
+       case('cvahpot')                                                                             
+         do field = 1, nActiveFields
+           if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_SOLUTE_ID) then
+             do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents
+               if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == COMPONENT_SOLUTE_CVaHPOT_ID) then
+                 loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = &
+                     .true.
+                 loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
+                     IO_floatValue(line,chunkPos,i+1_pInt)
+               endif
+             enddo
+           endif
+         enddo  
+
+     end select
+ enddo; enddo
+ close(FILEUNIT)
+ 
+!--------------------------------------------------------------------------------------------------
+! consistency checks and output of load case
+ loadCases(1)%followFormerTrajectory = .false.                                                      ! cannot guess along trajectory for first inc of first currentLoadCase
+ errorID = 0_pInt
+ if (worldrank == 0) then
+   checkLoadcases: do currentLoadCase = 1_pInt, size(loadCases)
+     write (loadcase_string, '(i6)' ) currentLoadCase
+     write(6,'(1x,a,i6)') 'load case: ', currentLoadCase
+     if (.not. loadCases(currentLoadCase)%followFormerTrajectory) &
+       write(6,'(2x,a)') 'drop guessing along trajectory'
+     do field = 1_pInt, nActiveFields
+       select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
+         case(FIELD_MECH_ID)
+           write(6,'(2x,a)') 'Field '//trim(FIELD_MECH_label)
+       
+         case(FIELD_THERMAL_ID)
+           write(6,'(2x,a)') 'Field '//trim(FIELD_THERMAL_label)
+              
+         case(FIELD_DAMAGE_ID)
+           write(6,'(2x,a)') 'Field '//trim(FIELD_DAMAGE_label)
+           
+         case(FIELD_MGTWIN_ID)
+           write(6,'(2x,a)') 'Field '//trim(FIELD_MGTWIN_label)  
+       
+         case(FIELD_SOLUTE_ID)
+           write(6,'(2x,a)') 'Field '//trim(FIELD_SOLUTE_label)
+       
+       end select
+       do faceSet = 1_pInt, mesh_Nboundaries
+          do component = 1_pInt, loadCases(currentLoadCase)%fieldBC(field)%nComponents
+            if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask(faceSet)) &
+              write(6,'(4x,a,i2,a,i2,a,f12.7)') 'Face  ', mesh_boundaries(faceSet), &
+                                                ' Component ', component, & 
+                                                ' Value ', loadCases(currentLoadCase)%fieldBC(field)% &
+                                                             componentBC(component)%Value(faceSet)
+          enddo
+        enddo       
+     enddo
+     write(6,'(2x,a,f12.6)') 'time:       ', loadCases(currentLoadCase)%time
+     if (loadCases(currentLoadCase)%incs < 1_pInt)             errorID = 835_pInt                     ! non-positive incs count
+     write(6,'(2x,a,i5)')    'increments: ', loadCases(currentLoadCase)%incs
+     if (loadCases(currentLoadCase)%outputfrequency < 1_pInt)  errorID = 836_pInt                     ! non-positive result frequency
+     write(6,'(2x,a,i5)')    'output  frequency:  ', &
+                loadCases(currentLoadCase)%outputfrequency
+     write(6,'(2x,a,i5,/)')    'restart frequency:  ', &
+                loadCases(currentLoadCase)%restartfrequency
+     if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string)               ! exit with error message
+   enddo checkLoadcases
+ endif
+
+!--------------------------------------------------------------------------------------------------
+! doing initialization depending on selected solver 
+ call Utilities_init()
+ do field = 1, nActiveFields
+   select case (loadCases(1)%fieldBC(field)%ID)
+     case(FIELD_MECH_ID)
+       call FEM_mech_init(loadCases(1)%fieldBC(field))
+   end select
+ enddo   
+
+!--------------------------------------------------------------------------------------------------
+! loopping over loadcases
+ loadCaseLooping: do currentLoadCase = 1_pInt, size(loadCases)
+   time0 = time                                                                                     ! currentLoadCase start time                
+   if (loadCases(currentLoadCase)%followFormerTrajectory) then
+     guess = .true.
+   else
+     guess = .false.                                                                                ! change of load case, homogeneous guess for the first inc
+   endif
+
+!--------------------------------------------------------------------------------------------------
+! loop oper incs defined in input file for current currentLoadCase
+   incLooping: do inc = 1_pInt, loadCases(currentLoadCase)%incs
+     totalIncsCounter = totalIncsCounter + 1_pInt                                                 
+
+!--------------------------------------------------------------------------------------------------
+! forwarding time
+     timeIncOld = timeinc
+     if (loadCases(currentLoadCase)%logscale == 0_pInt) then                                        ! linear scale
+       timeinc = loadCases(currentLoadCase)%time/loadCases(currentLoadCase)%incs                    ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used
+     else
+       if (currentLoadCase == 1_pInt) then                                                          ! 1st currentLoadCase of logarithmic scale            
+         if (inc == 1_pInt) then                                                                    ! 1st inc of 1st currentLoadCase of logarithmic scale
+           timeinc = loadCases(1)%time*(2.0_pReal**real(    1_pInt-loadCases(1)%incs ,pReal))       ! assume 1st inc is equal to 2nd 
+         else                                                                                       ! not-1st inc of 1st currentLoadCase of logarithmic scale
+           timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1_pInt-loadCases(1)%incs ,pReal))
+         endif
+       else                                                                                         ! not-1st currentLoadCase of logarithmic scale
+         timeinc = time0 * &
+              ( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real(          inc,pReal)/&
+                                                    real(loadCases(currentLoadCase)%incs ,pReal))&
+               -(1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( (inc-1_pInt),pReal)/&
+                                                     real(loadCases(currentLoadCase)%incs ,pReal)))
+       endif
+     endif
+     timeinc = timeinc / 2.0_pReal**real(cutBackLevel,pReal)                                        ! depending on cut back level, decrease time step
+
+     forwarding: if(totalIncsCounter >= restartInc) then
+       stepFraction = 0_pInt
+
+!--------------------------------------------------------------------------------------------------
+! loop over sub incs 
+       subIncLooping: do while (stepFraction/subStepFactor**cutBackLevel <1_pInt)
+         time = time + timeinc                                                                      ! forward time
+         stepFraction = stepFraction + 1_pInt 
+         remainingLoadCaseTime = time0 - time + loadCases(currentLoadCase)%time + timeInc
+           
+!--------------------------------------------------------------------------------------------------
+! report begin of new increment
+         if (worldrank == 0) then
+           write(6,'(/,a)') ' ###########################################################################'
+           write(6,'(1x,a,es12.5'//&
+                   ',a,'//IO_intOut(inc)//',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//&
+                   ',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//&
+                   ',a,'//IO_intOut(currentLoadCase)//',a,'//IO_intOut(size(loadCases))//')') &
+                   'Time', time, &
+                   's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,&
+                   '-', stepFraction, '/', subStepFactor**cutBackLevel,&
+                   ' of load case ', currentLoadCase,'/',size(loadCases)
+           flush(6)
+           write(incInfo,'(a,'//IO_intOut(totalIncsCounter)//',a,'//IO_intOut(sum(loadCases%incs))//&
+                 ',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//')') &
+                 'Increment ',totalIncsCounter,'/',sum(loadCases%incs),&
+                 '-',stepFraction, '/', subStepFactor**cutBackLevel
+         endif
+
+!--------------------------------------------------------------------------------------------------
+! forward fields
+         do field = 1, nActiveFields
+           select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
+             case(FIELD_MECH_ID)
+               call FEM_mech_forward (&
+                   guess,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field))
+
+          end select
+         enddo       
+           
+!--------------------------------------------------------------------------------------------------
+! solve fields
+         stagIter = 0_pInt
+         stagIterate = .true.
+         do while (stagIterate)
+           do field = 1, nActiveFields
+             select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
+               case(FIELD_MECH_ID)
+                 solres(field) = FEM_mech_solution (&
+                       incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field))
+
+             end select
+             if(.not. solres(field)%converged) exit                                                ! no solution found
+           enddo
+           stagIter = stagIter + 1_pInt
+           stagIterate = stagIter < stagItMax .and. &
+                         all(solres(:)%converged) .and. &
+                         .not. all(solres(:)%stagConverged)
+         enddo     
+           
+! check solution 
+         cutBack = .False.                                                                   
+         if(.not. all(solres(:)%converged .and. solres(:)%stagConverged)) then                      ! no solution found
+           if (cutBackLevel < maxCutBack) then                                                   ! do cut back
+             if (worldrank == 0) &
+               write(6,'(/,a)') ' cut back detected'
+             cutBack = .True.
+             stepFraction = (stepFraction - 1_pInt) * subStepFactor                                 ! adjust to new denominator
+             cutBackLevel = cutBackLevel + 1_pInt
+             time    = time - timeinc                                                               ! rewind time
+             timeinc = timeinc/2.0_pReal
+           else                                                                                     ! default behavior, exit if spectral solver does not converge                                
+             call IO_warning(850_pInt)
+             call quit(-1_pInt*(lastRestartWritten+1_pInt))                                         ! quit and provide information about last restart inc written (e.g. for regridding)                                                                                       ! continue from non-converged solution and start guessing after accepted (sub)inc
+           endif
+         else
+           guess = .true.                                                                           ! start guessing after first converged (sub)inc
+           timeIncOld = timeinc
+         endif
+         if (.not. cutBack) then
+           if (worldrank == 0)  write(statUnit,*) totalIncsCounter, time, cutBackLevel, &
+                             solres%converged, solres%iterationsNeeded                              ! write statistics about accepted solution
+         endif
+       enddo subIncLooping
+       cutBackLevel = max(0_pInt, cutBackLevel - 1_pInt)                                            ! try half number of subincs next inc
+       if(all(solres(:)%converged)) then                                                                    ! report converged inc
+         convergedCounter = convergedCounter + 1_pInt
+         if (worldrank == 0) then
+         write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') &
+                                     ' increment ', totalIncsCounter, ' converged'
+         endif
+       else
+         if (worldrank == 0) then
+         write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') &                                     ! report non-converged inc
+                                     ' increment ', totalIncsCounter, ' NOT converged'
+         endif
+         notConvergedCounter = notConvergedCounter + 1_pInt
+       endif; flush(6)
+       if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then                      ! at output frequency
+         if (worldrank == 0) then
+         write(6,'(1/,a)') ' ... writing results to file ......................................'
+         endif
+       endif
+       if( loadCases(currentLoadCase)%restartFrequency > 0_pInt .and. &                             ! at frequency of writing restart information set restart parameter for FEsolving 
+                      mod(inc,loadCases(currentLoadCase)%restartFrequency) == 0_pInt) then          ! ToDo first call to CPFEM_general will write? 
+         restartWrite = .true.
+         lastRestartWritten = inc
+       endif 
+     else forwarding
+       time = time + timeinc
+       guess = .true.
+     endif forwarding
+
+    enddo incLooping
+ enddo loadCaseLooping
+ 
+!--------------------------------------------------------------------------------------------------
+! report summary of whole calculation
+ if (worldrank == 0) then
+ write(6,'(/,a)') ' ###########################################################################'
+ write(6,'(1x,i6.6,a,i6.6,a,f5.1,a)') convergedCounter, ' out of ', &
+                                   notConvergedCounter + convergedCounter, ' (', &
+                                   real(convergedCounter, pReal)/&
+                                   real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, &
+                                   ' %) increments converged!'
+ endif
+ if (notConvergedCounter > 0_pInt) call quit(3_pInt)                                                ! error if some are not converged
+ call quit(0_pInt)                                                                                  ! no complains ;)
+
+end program DAMASK_FEM 
+
+
+!--------------------------------------------------------------------------------------------------
+!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
+!> @brief quit subroutine to mimic behavior of FEM solvers
+!> @details exits the Spectral solver and reports time and duration. Exit code 0 signals
+!> everything went fine. Exit code 1 signals an error, message according to IO_error. Exit code 
+!> 2 signals request for regridding, increment of last saved restart information is written to
+!> stderr. Exit code 3 signals no severe problems, but some increments did not converge
+!--------------------------------------------------------------------------------------------------
+subroutine quit(stop_id)
+ use prec, only: &
+   pInt
+   
+ implicit none
+ integer(pInt), intent(in) :: stop_id
+ integer, dimension(8) :: dateAndTime                                                               ! type default integer
+
+ call date_and_time(values = dateAndTime)
+ write(6,'(/,a)') 'DAMASK terminated on:'
+ write(6,'(a,2(i2.2,a),i4.4)') 'Date:               ',dateAndTime(3),'/',&
+                                                      dateAndTime(2),'/',&
+                                                      dateAndTime(1) 
+ write(6,'(a,2(i2.2,a),i2.2)') 'Time:               ',dateAndTime(5),':',&
+                                                      dateAndTime(6),':',&
+                                                      dateAndTime(7)  
+ if (stop_id == 0_pInt) stop 0                                                                      ! normal termination
+ if (stop_id <  0_pInt) then                                                                        ! trigger regridding
+   write(0,'(a,i6)') 'restart information available at ', stop_id*(-1_pInt)
+   stop 2
+ endif
+ if (stop_id == 3_pInt) stop 3                                                                      ! not all incs converged
+ stop 1                                                                                             ! error (message from IO_error)
+
+end subroutine quit

src/DAMASK_interface.f90

@@ -1,9 +1,11 @@
 !--------------------------------------------------------------------------------------------------
+!> @author   Jaeyong Jung, Max-Planck-Institut für Eisenforschung GmbH
+!> @author   Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
 !> @author   Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
 !> @author   Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief    Interfacing between the spectral solver and the material subroutines provided
+!> @brief    Interfacing between the PETSc-based solvers and the material subroutines provided
 !!           by DAMASK
-!> @details  Interfacing between the spectral solver and the material subroutines provided
+!> @details  Interfacing between the PETSc-based solvers and the material subroutines provided
 !>           by DAMASK. Interpretating the command line arguments to get load case, geometry file, 
 !>           and working directory.
 !--------------------------------------------------------------------------------------------------
@@ -13,12 +15,11 @@ module DAMASK_interface
 
  implicit none
  private
- logical,             public, protected :: appendToOutFile = .false.                                !< Append to existing spectralOut file (in case of restart, not in case of regridding)
+ integer(pInt),       public, protected :: &
- integer(pInt),       public, protected :: spectralRestartInc = 0_pInt                              !< Increment at which calculation starts
+   interface_restartInc = 0_pInt                                                                    !< Increment at which calculation starts
  character(len=1024), public, protected :: &
    geometryFile = '', &                                                                             !< parameter given for geometry file
    loadCaseFile = ''                                                                                !< parameter given for load case file
- character(len=1024), private           :: workingDirectory
 
  public :: &
    getSolverJobName, &
@@ -44,22 +45,36 @@ subroutine DAMASK_interface_init()
 #include <petsc/finclude/petscsys.h>
 #if PETSC_VERSION_MAJOR!=3 || PETSC_VERSION_MINOR!=9
 ===================================================================================================
+  3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x
+===================================================================================================
+=======   THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x   ===========================================
+==========   THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x   ========================================
+=============   THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x   =====================================
+================   THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x   ==================================
+===================   THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x   ===============================
+======================   THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x   ============================
 =========================   THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x   ========================= 
+============================   THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x   ======================
+===============================   THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x   ===================
+==================================   THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x   ================
+=====================================   THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x   =============
+========================================   THIS VERSION OF DAMASK REQUIRES PETSc 3.9.x   ==========
+===================================================================================================
+  3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x 3.9.x
 ===================================================================================================
 #endif
  use PETScSys
  use system_routines, only: &
-   getHostName
+   getHostName, &
+   getCWD
 
  implicit none
  character(len=1024) :: &
    commandLine, &                                                                                   !< command line call as string
-   loadcaseArg   = '', &                                                                            !< -l argument given to DAMASK_spectral.exe
+   loadcaseArg   = '', &                                                                            !< -l argument given to the executable
-   geometryArg   = '', &                                                                            !< -g argument given to DAMASK_spectral.exe
+   geometryArg   = '', &                                                                            !< -g argument given to the executable
-   workingDirArg = '', &                                                                            !< -w argument given to DAMASK_spectral.exe
+   workingDirArg = '', &                                                                            !< -w argument given to the executable
-   hostName, &                                                                                      !< name of machine on which DAMASK_spectral.exe is execute (might require export HOSTNAME)
+   userName                                                                                         !< name of user calling the executable
-   userName, &                                                                                      !< name of user calling DAMASK_spectral.exe
-   tag
  integer :: &
    i, &
 #ifdef _OPENMP
@@ -72,7 +87,6 @@ subroutine DAMASK_interface_init()
  integer, dimension(8) :: &
    dateAndTime                                                                                      ! type default integer
  PetscErrorCode :: ierr
- logical        :: error
  external :: &
    quit,&
    PETScErrorF, &                                                                                   ! is called in the CHKERRQ macro
@@ -110,7 +124,7 @@ subroutine DAMASK_interface_init()
  endif mainProcess
 
  call date_and_time(values = dateAndTime)
- write(6,'(/,a)') ' <<<+-  DAMASK_spectral  -+>>>'
+ write(6,'(/,a)') ' <<<+-  DAMASK_interface init  -+>>>'
  write(6,'(a,/)') ' Roters et al., Computational Materials Science, 2018'
  write(6,'(/,a)')              ' Version: '//DAMASKVERSION
  write(6,'(a,2(i2.2,a),i4.4)') ' Date:    ',dateAndTime(3),'/',&
@@ -120,7 +134,6 @@ subroutine DAMASK_interface_init()
                                             dateAndTime(6),':',&
                                             dateAndTime(7)  
  write(6,'(/,a,i4.1)') ' MPI processes: ',worldsize
- write(6,'(/,a)') ' <<<+-  DAMASK_interface init  -+>>>'
 #include "compilation_info.f90"
  
  call get_command(commandLine)
@@ -129,9 +142,8 @@ subroutine DAMASK_interface_init()
    select case(IIO_stringValue(commandLine,chunkPos,i))                                             ! extract key
      case ('-h','--help')
        write(6,'(a)')  ' #######################################################################'
-       write(6,'(a)')  ' DAMASK_spectral:'
+       write(6,'(a)')  ' DAMASK Command Line Interface:'
-       write(6,'(a)')  ' The spectral method boundary value problem solver for'
+       write(6,'(a)')  ' For PETSc-based solvers for the Düsseldorf Advanced Material Simulation Kit'
-       write(6,'(a)')  ' the Düsseldorf Advanced Material Simulation Kit'
        write(6,'(a,/)')' #######################################################################'
        write(6,'(a,/)')' Valid command line switches:'
        write(6,'(a)')  '    --geom         (-g, --geometry)'
@@ -141,23 +153,14 @@ subroutine DAMASK_interface_init()
        write(6,'(a)')  '    --help         (-h)'
        write(6,'(/,a)')' -----------------------------------------------------------------------'
        write(6,'(a)')  ' Mandatory arguments:'
-       write(6,'(/,a)')'   --geom PathToGeomFile/NameOfGeom.geom'
+       write(6,'(/,a)')'   --geom PathToGeomFile/NameOfGeom'
-       write(6,'(a)')  '        Specifies the location of the geometry definition file,' 
+       write(6,'(a)')  '        Specifies the location of the geometry definition file.'
-       write(6,'(a)')  '            if no extension is given, .geom will be appended.' 
+       write(6,'(/,a)')'   --load PathToLoadFile/NameOfLoadFile'
-       write(6,'(a)')  '        "PathToGeomFile" will be the working directory if not specified'
+       write(6,'(a)')  '        Specifies the location of the load case definition file.'
-       write(6,'(a)')  '            via --workingdir.'
-       write(6,'(a)')  '        Make sure the file "material.config" exists in the working'
-       write(6,'(a)')  '            directory.'   
-       write(6,'(a)')  '        For further configuration place "numerics.config"'
-       write(6,'(a)')'            and "numerics.config" in that directory.'
-       write(6,'(/,a)')'   --load PathToLoadFile/NameOfLoadFile.load'
-       write(6,'(a)')  '        Specifies the location of the load case definition file,' 
-       write(6,'(a)')  '            if no extension is given, .load will be appended.' 
        write(6,'(/,a)')' -----------------------------------------------------------------------'
        write(6,'(a)')  ' Optional arguments:'
        write(6,'(/,a)')'   --workingdirectory PathToWorkingDirectory'
-       write(6,'(a)')  '        Specifies the working directory and overwrites the default'
+       write(6,'(a)')  '        Specifies the working directory and overwrites the default ./'
-       write(6,'(a)')  '            "PathToGeomFile".'
        write(6,'(a)')  '        Make sure the file "material.config" exists in the working'
        write(6,'(a)')  '            directory.'   
        write(6,'(a)')  '        For further configuration place "numerics.config"'
@@ -166,7 +169,7 @@ subroutine DAMASK_interface_init()
        write(6,'(a)')  '        Reads in increment XX and continues with calculating'
        write(6,'(a)')  '            increment XX+1 based on this.'
        write(6,'(a)')  '        Appends to existing results file'
-       write(6,'(a)')  '            "NameOfGeom_NameOfLoadFile.spectralOut".'
+       write(6,'(a)')  '            "NameOfGeom_NameOfLoadFile".'
        write(6,'(a)')  '        Works only if the restart information for increment XX'
        write(6,'(a)')  '            is available in the working directory.'
        write(6,'(/,a)')' -----------------------------------------------------------------------'
@@ -182,8 +185,7 @@ subroutine DAMASK_interface_init()
        if (i < chunkPos(1)) workingDirArg = trim(IIO_stringValue(commandLine,chunkPos,i+1_pInt))
      case ('-r', '--rs', '--restart')
        if (i < chunkPos(1)) then
-         spectralRestartInc = IIO_IntValue(commandLine,chunkPos,i+1_pInt)
+         interface_restartInc = IIO_IntValue(commandLine,chunkPos,i+1_pInt)
-         appendToOutFile = .true.
        endif
    end select
  enddo
@@ -193,26 +195,25 @@ subroutine DAMASK_interface_init()
    call quit(1_pInt)
  endif
 
- workingDirectory = trim(setWorkingDirectory(trim(workingDirArg)))
+ if (len_trim(workingDirArg) > 0) call setWorkingDirectory(trim(workingDirArg))
  geometryFile = getGeometryFile(geometryArg)
  loadCaseFile = getLoadCaseFile(loadCaseArg)
 
  call get_environment_variable('USER',userName)
- error = getHostName(hostName)
+ ! ToDo: https://stackoverflow.com/questions/8953424/how-to-get-the-username-in-c-c-in-linux
- write(6,'(a,a)')      ' Host name:              ', trim(hostName)
+ write(6,'(a,a)')      ' Host name:              ', trim(getHostName())
  write(6,'(a,a)')      ' User name:              ', trim(userName)
  write(6,'(a,a)')      ' Command line call:      ', trim(commandLine)
  if (len(trim(workingDirArg)) > 0) &
    write(6,'(a,a)')    ' Working dir argument:   ', trim(workingDirArg)
  write(6,'(a,a)')      ' Geometry argument:      ', trim(geometryArg)
  write(6,'(a,a)')      ' Loadcase argument:      ', trim(loadcaseArg)
- write(6,'(a,a)')      ' Working directory:      ', trim(workingDirectory)
+ write(6,'(a,a)')      ' Working directory:      ', trim(getCWD())
  write(6,'(a,a)')      ' Geometry file:          ', trim(geometryFile)
  write(6,'(a,a)')      ' Loadcase file:          ', trim(loadCaseFile)
  write(6,'(a,a)')      ' Solver job name:        ', trim(getSolverJobName())
- if (SpectralRestartInc > 0_pInt) &
+ if (interface_restartInc > 0_pInt) &
-   write(6,'(a,i6.6)') ' Restart from increment: ', spectralRestartInc
+   write(6,'(a,i6.6)') ' Restart from increment: ', interface_restartInc
- write(6,'(a,l1,/)')   ' Append to result file:  ', appendToOutFile
 
 end subroutine DAMASK_interface_init
 
@@ -221,38 +222,32 @@ end subroutine DAMASK_interface_init
 !> @brief extract working directory from given argument or from location of geometry file,
 !!        possibly converting relative arguments to absolut path
 !--------------------------------------------------------------------------------------------------
-character(len=1024) function setWorkingDirectory(workingDirectoryArg)
+subroutine setWorkingDirectory(workingDirectoryArg)
  use system_routines, only: &
    getCWD, &
    setCWD
 
  implicit none
  character(len=*),  intent(in) :: workingDirectoryArg                                               !< working directory argument
- logical                       :: error
+ character(len=1024)           :: workingDirectory                                                  !< working directory argument
  external                      :: quit
+ logical                       :: error
 
- wdGiven: if (len(workingDirectoryArg)>0) then
  absolutePath: if (workingDirectoryArg(1:1) == '/') then
-     setWorkingDirectory = workingDirectoryArg
+   workingDirectory = workingDirectoryArg
  else absolutePath
-     error = getCWD(setWorkingDirectory)
+   workingDirectory = getCWD()
-     if (error) call quit(1_pInt)
+   workingDirectory = trim(workingDirectory)//'/'//workingDirectoryArg
-     setWorkingDirectory = trim(setWorkingDirectory)//'/'//workingDirectoryArg
  endif absolutePath
- else wdGiven
-   error = getCWD(setWorkingDirectory)                                                              ! relative path given as command line argument
-   if (error) call quit(1_pInt)
- endif wdGiven
 
- setWorkingDirectory = trim(rectifyPath(setWorkingDirectory))
+ workingDirectory = trim(rectifyPath(workingDirectory))
-
+ error = setCWD(trim(workingDirectory))
- error = setCWD(trim(setWorkingDirectory))
  if(error) then
-   write(6,'(a20,a,a16)') ' working directory "',trim(setWorkingDirectory),'" does not exist'
+   write(6,'(a20,a,a16)') ' working directory "',trim(workingDirectory),'" does not exist'
    call quit(1_pInt)
  endif
 
-end function setWorkingDirectory
+end subroutine setWorkingDirectory
 
 
 !--------------------------------------------------------------------------------------------------
@@ -284,22 +279,15 @@ end function getSolverJobName
 !> @brief basename of geometry file with extension from command line arguments
 !--------------------------------------------------------------------------------------------------
 character(len=1024) function getGeometryFile(geometryParameter)
+ use system_routines, only: &
+   getCWD
 
  implicit none
- character(len=1024), intent(in) :: &
+ character(len=1024), intent(in) :: geometryParameter
-   geometryParameter
- integer :: posExt, posSep
- external  :: quit
 
  getGeometryFile = trim(geometryParameter)
- posExt = scan(getGeometryFile,'.',back=.true.)
+ if (scan(getGeometryFile,'/') /= 1) getGeometryFile = trim(getCWD())//'/'//trim(getGeometryFile)
- posSep = scan(getGeometryFile,'/',back=.true.)
+ getGeometryFile = makeRelativePath(trim(getCWD()), getGeometryFile)
-
- if (posExt <= posSep) getGeometryFile = trim(getGeometryFile)//('.geom')
- if (scan(getGeometryFile,'/') /= 1) &
-   getGeometryFile = trim(workingDirectory)//'/'//trim(getGeometryFile)
-
- getGeometryFile = makeRelativePath(workingDirectory, getGeometryFile)
 
 
 end function getGeometryFile
@@ -309,22 +297,15 @@ end function getGeometryFile
 !> @brief relative path of loadcase from command line arguments
 !--------------------------------------------------------------------------------------------------
 character(len=1024) function getLoadCaseFile(loadCaseParameter)
+ use system_routines, only: &
+   getCWD
 
  implicit none
- character(len=1024), intent(in) :: &
+ character(len=1024), intent(in) :: loadCaseParameter
-   loadCaseParameter
- integer :: posExt, posSep
- external  :: quit
 
  getLoadCaseFile = trim(loadCaseParameter)
- posExt = scan(getLoadCaseFile,'.',back=.true.)
+ if (scan(getLoadCaseFile,'/') /= 1) getLoadCaseFile = trim(getCWD())//'/'//trim(getLoadCaseFile)
- posSep = scan(getLoadCaseFile,'/',back=.true.)
+ getLoadCaseFile = makeRelativePath(trim(getCWD()), getLoadCaseFile)
-
- if (posExt <= posSep) getLoadCaseFile = trim(getLoadCaseFile)//('.load')
- if (scan(getLoadCaseFile,'/') /= 1) &
-   getLoadCaseFile = trim(workingDirectory)//'/'//trim(getLoadCaseFile)
-
- getLoadCaseFile = makeRelativePath(workingDirectory, getLoadCaseFile)
 
 end function getLoadCaseFile
 
@@ -337,21 +318,20 @@ function rectifyPath(path)
 
  implicit none
  character(len=*) :: path
- character(len=len_trim(path)) :: rectifyPath
+ character(len=1024) :: rectifyPath
  integer :: i,j,k,l                                                                                 ! no pInt
 
 !--------------------------------------------------------------------------------------------------
 ! remove /./ from path
- l = len_trim(path)
+ rectifyPath = trim(path)
- rectifyPath = path
+ l = len_trim(rectifyPath)
  do i = l,3,-1
    if (rectifyPath(i-2:i) == '/./') rectifyPath(i-1:l) = rectifyPath(i+1:l)//'  '
  enddo
 
 !--------------------------------------------------------------------------------------------------
 ! remove // from path
- l = len_trim(path)
+ l = len_trim(rectifyPath)
- rectifyPath = path
  do i = l,2,-1
    if (rectifyPath(i-1:i) == '//') rectifyPath(i-1:l) = rectifyPath(i:l)//' '
  enddo

src/DAMASK_spectral.f90

@@ -20,14 +20,12 @@ program DAMASK_spectral
    pReal, &
    tol_math_check, &
    dNeq
- use system_routines, only: &
-   getCWD
  use DAMASK_interface, only: &
    DAMASK_interface_init, &
    loadCaseFile, &
    geometryFile, &
    getSolverJobName, &
-   appendToOutFile
+   interface_restartInc
  use IO, only: &
    IO_read, &
    IO_isBlank, &
@@ -383,8 +381,7 @@ program DAMASK_spectral
 !--------------------------------------------------------------------------------------------------
 ! write header of output file
  if (worldrank == 0) then
-   if (.not. appendToOutFile) then                                                                  ! after restart, append to existing results file
+   writeHeader: if (interface_restartInc < 1_pInt) then
-     if (getCWD(workingDir)) call IO_error(106_pInt,ext_msg=trim(workingDir))
      open(newunit=resUnit,file=trim(getSolverJobName())//&
                                  '.spectralOut',form='UNFORMATTED',status='REPLACE')
      write(resUnit) 'load:',       trim(loadCaseFile)                                               ! ... and write header
@@ -407,10 +404,10 @@ program DAMASK_spectral
      if (iand(debug_level(debug_spectral),debug_levelBasic) /= 0) &
        write(6,'(/,a)') ' header of result and statistics file written out'
      flush(6)
-   else                                                                                             ! open new files ...
+   else writeHeader
      open(newunit=statUnit,file=trim(getSolverJobName())//&
                                  '.sta',form='FORMATTED', position='APPEND', status='OLD')
-   endif
+   endif writeHeader
  endif
 
 !--------------------------------------------------------------------------------------------------
@@ -431,7 +428,7 @@ program DAMASK_spectral
  call MPI_file_seek (resUnit,fileOffset,MPI_SEEK_SET,ierr)
  if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_seek')
 
- if (.not. appendToOutFile) then                                                                    ! if not restarting, write 0th increment
+ writeUndeformed: if (interface_restartInc < 1_pInt) then
    write(6,'(1/,a)') ' ... writing initial configuration to file ........................'
    do i = 1, size(materialpoint_results,3)/(maxByteOut/(materialpoint_sizeResults*pReal))+1         ! slice the output of my process in chunks not exceeding the limit for one output
      outputIndex = int([(i-1_pInt)*((maxRealOut)/materialpoint_sizeResults)+1_pInt, &               ! QUESTION: why not starting i at 0 instead of murky 1?
@@ -443,7 +440,7 @@ program DAMASK_spectral
      if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_write')
    enddo
    fileOffset = fileOffset + sum(outputSize)                                                        ! forward to current file position
- endif
+ endif writeUndeformed
 !--------------------------------------------------------------------------------------------------
 ! looping over loadcases
  loadCaseLooping: do currentLoadCase = 1_pInt, size(loadCases)

src/FEM_mech.f90

@@ -0,0 +1,734 @@
+!--------------------------------------------------------------------------------------------------
+!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
+!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
+!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
+!> @brief FEM PETSc solver
+!--------------------------------------------------------------------------------------------------
+module FEM_mech
+#include <petsc/finclude/petsc.h>
+
+use PETScdmda
+use PETScsnes
+use PETScDM
+use PETScDMplex
+ use prec, only: & 
+   pInt, &
+   pReal
+ use math, only: &
+   math_I3
+ use FEM_utilities, only: &
+   tSolutionState, &
+   tFieldBC, &
+   tComponentBC
+ use numerics, only: &
+   worldrank, &
+   worldsize  
+ use mesh, only: &
+   mesh_Nboundaries, &
+   mesh_boundaries
+
+ implicit none
+ private
+   
+!--------------------------------------------------------------------------------------------------
+! derived types
+ type tSolutionParams 
+   type(tFieldBC)  :: fieldBC
+   real(pReal)     :: timeinc
+   real(pReal)     :: timeincOld
+ end type tSolutionParams
+ 
+ type(tSolutionParams),   private :: params
+
+!--------------------------------------------------------------------------------------------------
+! PETSc data
+ SNES,                           private :: mech_snes
+ Vec,                            private :: solution, solution_rate, solution_local
+ PetscInt,                       private :: dimPlex, cellDof, nQuadrature, nBasis
+ PetscReal, allocatable, target,dimension(:), private :: qPoints, qWeights
+ MatNullSpace,                   private :: matnull
+
+!--------------------------------------------------------------------------------------------------
+! stress, stiffness and compliance average etc.
+ character(len=1024), private :: incInfo   
+ real(pReal), private, dimension(3,3) :: &
+   P_av = 0.0_pReal
+ logical, private :: ForwardData
+ real(pReal), parameter, private :: eps = 1.0e-18_pReal
+
+ public :: &
+   FEM_mech_init, &
+   FEM_mech_solution ,&
+   FEM_mech_forward, &
+   FEM_mech_destroy
+
+ external :: &
+   MatZeroRowsColumnsLocalIS, &
+   PetscQuadratureCreate, &
+   PetscFECreateDefault, &
+   PetscFESetQuadrature, &
+   PetscFEGetDimension, &
+   PetscFEDestroy, &
+   PetscFEGetDualSpace, &
+   PetscQuadratureDestroy, &
+   PetscDSSetDiscretization, &
+   PetscDSGetTotalDimension, &
+   PetscDSGetDiscretization, &
+   PetscDualSpaceGetFunctional, &
+   DMGetLabelSize, &
+   DMSNESSetFunctionLocal, &
+   DMSNESSetJacobianLocal, &
+   SNESSetOptionsPrefix, &
+   SNESSetConvergenceTest, &
+   PetscObjectSetName
+   
+contains
+
+!--------------------------------------------------------------------------------------------------
+!> @brief allocates all neccessary fields and fills them with data, potentially from restart info
+!--------------------------------------------------------------------------------------------------
+subroutine FEM_mech_init(fieldBC)
+ use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment)
+ use IO, only: &
+   IO_timeStamp, &
+   IO_error
+ use DAMASK_interface, only: &
+   getSolverJobName
+ use mesh, only: &
+   geomMesh
+ use numerics, only: &
+   worldrank, &
+   itmax, &
+   integrationOrder
+ use FEM_Zoo, only: &
+   FEM_Zoo_nQuadrature, &
+   FEM_Zoo_QuadraturePoints, &
+   FEM_Zoo_QuadratureWeights  
+
+ implicit none
+ type(tFieldBC),             intent(in) :: fieldBC
+ DM                                     :: mech_mesh
+ PetscFE                                :: mechFE
+ PetscQuadrature                        :: mechQuad, functional
+ PetscDS                                :: mechDS
+ PetscDualSpace                         :: mechDualSpace
+ DMLabel                                :: BCLabel
+ PetscInt,          allocatable, target :: numComp(:), numDoF(:), bcField(:)
+ PetscInt,                      pointer :: pNumComp(:), pNumDof(:), pBcField(:), pBcPoint(:)
+ PetscInt                               :: numBC, bcSize
+ IS                                     :: bcPoint
+ IS,                allocatable, target :: bcComps(:), bcPoints(:)
+ IS,                            pointer :: pBcComps(:), pBcPoints(:)
+ PetscSection                           :: section
+ PetscInt                               :: field, faceSet, topologDim, nNodalPoints
+ PetscReal,                     pointer :: qPointsP(:), qWeightsP(:), &
+                                           nodalPointsP(:), nodalWeightsP(:)
+ PetscReal,         allocatable, target :: nodalPoints(:), nodalWeights(:)
+ PetscScalar,                   pointer :: px_scal(:)
+ PetscScalar,       allocatable, target ::  x_scal(:)
+ PetscReal                              :: detJ
+ PetscReal,         allocatable, target :: v0(:), cellJ(:), invcellJ(:), cellJMat(:,:)
+ PetscReal,                     pointer :: pV0(:), pCellJ(:), pInvcellJ(:)
+ PetscInt                               :: cellStart, cellEnd, cell, basis 
+ character(len=7)                       :: prefix = 'mechFE_'
+ PetscErrorCode                         :: ierr
+ 
+ write(6,'(/,a)') ' <<<+-  FEM_mech init  -+>>>'
+ write(6,'(a15,a)')   ' Current time: ',IO_timeStamp()
+#include "compilation_info.f90"
+
+!--------------------------------------------------------------------------------------------------
+! Setup FEM mech mesh
+ call DMClone(geomMesh,mech_mesh,ierr); CHKERRQ(ierr)
+ call DMGetDimension(mech_mesh,dimPlex,ierr); CHKERRQ(ierr)
+
+!--------------------------------------------------------------------------------------------------
+! Setup FEM mech discretization
+ allocate(qPoints(dimPlex*FEM_Zoo_nQuadrature(dimPlex,integrationOrder)))
+ allocate(qWeights(FEM_Zoo_nQuadrature(dimPlex,integrationOrder)))
+ qPoints = FEM_Zoo_QuadraturePoints(dimPlex,integrationOrder)%p
+ qWeights = FEM_Zoo_QuadratureWeights(dimPlex,integrationOrder)%p
+ nQuadrature = FEM_Zoo_nQuadrature(dimPlex,integrationOrder)
+ qPointsP => qPoints
+ qWeightsP => qWeights
+ call PetscQuadratureCreate(PETSC_COMM_SELF,mechQuad,ierr); CHKERRQ(ierr)
+ call PetscQuadratureSetData(mechQuad,dimPlex,nQuadrature,qPointsP,qWeightsP,ierr)
+ CHKERRQ(ierr)
+ call PetscFECreateDefault(mech_mesh,dimPlex,dimPlex,PETSC_TRUE,prefix, &
+                           integrationOrder,mechFE,ierr); CHKERRQ(ierr)
+ call PetscFESetQuadrature(mechFE,mechQuad,ierr); CHKERRQ(ierr)
+ call PetscFEGetDimension(mechFE,nBasis,ierr); CHKERRQ(ierr)
+ call DMGetDS(mech_mesh,mechDS,ierr); CHKERRQ(ierr)
+ call PetscDSAddDiscretization(mechDS,mechFE,ierr); CHKERRQ(ierr)
+ call PetscDSGetTotalDimension(mechDS,cellDof,ierr); CHKERRQ(ierr)
+ call PetscFEDestroy(mechFE,ierr); CHKERRQ(ierr)
+ call PetscQuadratureDestroy(mechQuad,ierr); CHKERRQ(ierr)
+
+!--------------------------------------------------------------------------------------------------
+! Setup FEM mech boundary conditions
+ call DMGetLabel(mech_mesh,'Face Sets',BCLabel,ierr); CHKERRQ(ierr)
+ call DMPlexLabelComplete(mech_mesh,BCLabel,ierr); CHKERRQ(ierr)
+ call DMGetDefaultSection(mech_mesh,section,ierr); CHKERRQ(ierr)
+ allocate(numComp(1), source=dimPlex); pNumComp => numComp
+ allocate(numDof(dimPlex+1), source = 0); pNumDof  => numDof
+ do topologDim = 0, dimPlex
+   call DMPlexGetDepthStratum(mech_mesh,topologDim,cellStart,cellEnd,ierr)
+   CHKERRQ(ierr)
+   call PetscSectionGetDof(section,cellStart,numDof(topologDim+1),ierr)
+   CHKERRQ(ierr)
+ enddo
+ numBC = 0
+ do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries
+   if (fieldBC%componentBC(field)%Mask(faceSet)) numBC = numBC + 1
+ enddo; enddo
+ allocate(bcField(numBC), source=0); pBcField => bcField 
+ allocate(bcComps(numBC)); pBcComps => bcComps
+ allocate(bcPoints(numBC)); pBcPoints => bcPoints
+ numBC = 0
+ do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries
+   if (fieldBC%componentBC(field)%Mask(faceSet)) then
+     numBC = numBC + 1
+     call ISCreateGeneral(PETSC_COMM_WORLD,1,[field-1],PETSC_COPY_VALUES,bcComps(numBC),ierr)
+     CHKERRQ(ierr)
+     call DMGetStratumSize(mech_mesh,'Face Sets',mesh_boundaries(faceSet),bcSize,ierr)
+     CHKERRQ(ierr)
+     if (bcSize > 0) then
+       call DMGetStratumIS(mech_mesh,'Face Sets',mesh_boundaries(faceSet),bcPoint,ierr)
+       CHKERRQ(ierr)
+       call ISGetIndicesF90(bcPoint,pBcPoint,ierr); CHKERRQ(ierr)
+       call ISCreateGeneral(PETSC_COMM_WORLD,bcSize,pBcPoint,PETSC_COPY_VALUES,bcPoints(numBC),ierr)
+       CHKERRQ(ierr)
+       call ISRestoreIndicesF90(bcPoint,pBcPoint,ierr); CHKERRQ(ierr)
+       call ISDestroy(bcPoint,ierr); CHKERRQ(ierr)
+     else
+       call ISCreateGeneral(PETSC_COMM_WORLD,0,[0],PETSC_COPY_VALUES,bcPoints(numBC),ierr)
+       CHKERRQ(ierr)
+     endif  
+   endif  
+ enddo; enddo
+ call DMPlexCreateSection(mech_mesh,dimPlex,1,pNumComp,pNumDof, &
+                          numBC,pBcField,pBcComps,pBcPoints,PETSC_NULL_IS, &
+                          section,ierr)
+ CHKERRQ(ierr)
+ call DMSetDefaultSection(mech_mesh,section,ierr); CHKERRQ(ierr)
+ do faceSet = 1, numBC
+   call ISDestroy(bcPoints(faceSet),ierr); CHKERRQ(ierr)
+ enddo
+ 
+!--------------------------------------------------------------------------------------------------
+! initialize solver specific parts of PETSc
+ call SNESCreate(PETSC_COMM_WORLD,mech_snes,ierr);CHKERRQ(ierr) 
+ call SNESSetOptionsPrefix(mech_snes,'mech_',ierr);CHKERRQ(ierr) 
+ call SNESSetDM(mech_snes,mech_mesh,ierr); CHKERRQ(ierr)                                             !< set the mesh for non-linear solver
+ call DMCreateGlobalVector(mech_mesh,solution        ,ierr); CHKERRQ(ierr)                           !< locally owned displacement Dofs
+ call DMCreateGlobalVector(mech_mesh,solution_rate   ,ierr); CHKERRQ(ierr)                           !< locally owned velocity Dofs to guess solution at next load step
+ call DMCreateLocalVector (mech_mesh,solution_local  ,ierr); CHKERRQ(ierr)                           !< locally owned velocity Dofs to guess solution at next load step
+ call DMSNESSetFunctionLocal(mech_mesh,FEM_mech_formResidual,PETSC_NULL_VEC,ierr)                 !< function to evaluate residual forces
+ CHKERRQ(ierr)
+ call DMSNESSetJacobianLocal(mech_mesh,FEM_mech_formJacobian,PETSC_NULL_VEC,ierr)                 !< function to evaluate stiffness matrix
+ CHKERRQ(ierr)
+ call SNESSetMaxLinearSolveFailures(mech_snes, huge(1), ierr); CHKERRQ(ierr)                         !< ignore linear solve failures 
+ call SNESSetConvergenceTest(mech_snes,FEM_mech_converged,PETSC_NULL_VEC,PETSC_NULL_FUNCTION,ierr)
+ CHKERRQ(ierr)
+ call SNESSetTolerances(mech_snes,1.0,0.0,0.0,itmax,itmax,ierr)
+ CHKERRQ(ierr)
+ call SNESSetFromOptions(mech_snes,ierr); CHKERRQ(ierr)
+  
+!--------------------------------------------------------------------------------------------------
+! init fields
+ call VecSet(solution        ,0.0,ierr); CHKERRQ(ierr)
+ call VecSet(solution_rate   ,0.0,ierr); CHKERRQ(ierr)
+ allocate(x_scal(cellDof))
+ allocate(nodalPoints (dimPlex))
+ allocate(nodalWeights(1))
+ nodalPointsP  => nodalPoints
+ nodalWeightsP => nodalWeights
+ allocate(v0(dimPlex))
+ allocate(cellJ(dimPlex*dimPlex))
+ allocate(invcellJ(dimPlex*dimPlex))
+ allocate(cellJMat(dimPlex,dimPlex))
+ pV0 => v0
+ pCellJ => cellJ
+ pInvcellJ => invcellJ
+ call DMGetDefaultSection(mech_mesh,section,ierr); CHKERRQ(ierr)
+ call DMGetDS(mech_mesh,mechDS,ierr); CHKERRQ(ierr)
+ call PetscDSGetDiscretization(mechDS,0,mechFE,ierr)
+ CHKERRQ(ierr)
+ call PetscFEGetDualSpace(mechFE,mechDualSpace,ierr); CHKERRQ(ierr)
+ call DMPlexGetHeightStratum(mech_mesh,0,cellStart,cellEnd,ierr)
+ CHKERRQ(ierr)
+ do cell = cellStart, cellEnd-1                                                                     !< loop over all elements 
+   x_scal = 0.0
+   call  DMPlexComputeCellGeometryAffineFEM(mech_mesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr) 
+   CHKERRQ(ierr)
+   cellJMat = reshape(pCellJ,shape=[dimPlex,dimPlex])
+   do basis = 0, nBasis-1
+     call PetscDualSpaceGetFunctional(mechDualSpace,basis,functional,ierr)
+     CHKERRQ(ierr)
+     call PetscQuadratureGetData(functional,dimPlex,nNodalPoints,nodalPointsP,nodalWeightsP,ierr)
+     CHKERRQ(ierr)
+     x_scal(basis*dimPlex+1:(basis+1)*dimPlex) = pV0 + matmul(transpose(cellJMat),nodalPointsP + 1.0)
+   enddo
+   px_scal => x_scal
+   call DMPlexVecSetClosure(mech_mesh,section,solution_local,cell,px_scal,INSERT_ALL_VALUES,ierr)
+   CHKERRQ(ierr) 
+ enddo 
+
+end subroutine FEM_mech_init
+  
+!--------------------------------------------------------------------------------------------------
+!> @brief solution for the FEM load step
+!--------------------------------------------------------------------------------------------------
+type(tSolutionState) function FEM_mech_solution( &
+             incInfoIn,timeinc,timeinc_old,fieldBC)
+ use numerics, only: &
+   itmax
+ use FEsolving, only: &
+   terminallyIll
+
+ implicit none
+!--------------------------------------------------------------------------------------------------
+! input data for solution
+ real(pReal), intent(in) :: &
+   timeinc, &                                                                                       !< increment in time for current solution
+   timeinc_old                                                                                      !< increment in time of last increment
+ type(tFieldBC),      intent(in) :: &
+   fieldBC
+ character(len=*), intent(in) :: &
+   incInfoIn
+ 
+!--------------------------------------------------------------------------------------------------
+! 
+ PetscErrorCode :: ierr   
+ SNESConvergedReason :: reason
+
+ incInfo = incInfoIn
+ FEM_mech_solution%converged =.false.
+!--------------------------------------------------------------------------------------------------
+! set module wide availabe data 
+ params%timeinc = timeinc
+ params%timeincOld = timeinc_old
+ params%fieldBC = fieldBC
+
+ call SNESSolve(mech_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr)                           ! solve mech_snes based on solution guess (result in solution)
+ call SNESGetConvergedReason(mech_snes,reason,ierr); CHKERRQ(ierr)                                  ! solution converged?
+ terminallyIll = .false.
+
+ if (reason < 1) then                                                                               ! 0: still iterating (will not occur), negative -> convergence error
+   FEM_mech_solution%converged = .false.
+   FEM_mech_solution%iterationsNeeded = itmax
+ else                                                                                               ! >= 1 proper convergence (or terminally ill)
+   FEM_mech_solution%converged = .true.
+   call SNESGetIterationNumber(mech_snes,FEM_mech_solution%iterationsNeeded,ierr)
+   CHKERRQ(ierr)
+ endif
+
+ write(6,'(/,a)') ' ==========================================================================='
+ flush(6) 
+
+end function FEM_mech_solution
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief forms the FEM residual vector
+!--------------------------------------------------------------------------------------------------
+subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
+ use numerics, only: &
+   BBarStabilisation
+ use FEM_utilities, only: &
+   utilities_projectBCValues, &
+   utilities_constitutiveResponse
+ use homogenization, only: &
+   materialpoint_F, &
+   materialpoint_P
+ use math, only: &
+   math_det33, &
+   math_inv33  
+ use FEsolving, only: &
+   terminallyIll
+
+ implicit none
+ DM                                 :: dm_local
+ PetscDS                            :: prob
+ Vec                                :: x_local, f_local, xx_local
+ PetscSection                       :: section
+ PetscScalar, dimension(:), pointer :: x_scal, pf_scal
+ PetscScalar,                target :: f_scal(cellDof) 
+ PetscReal                          :: detJ, IcellJMat(dimPlex,dimPlex)
+ PetscReal,                  target :: v0(dimPlex), cellJ(dimPlex*dimPlex), &
+                                                    invcellJ(dimPlex*dimPlex)
+ PetscReal,                 pointer :: pV0(:), pCellJ(:), pInvcellJ(:)
+ PetscReal,                 pointer :: basisField(:), basisFieldDer(:)
+ PetscInt                           :: cellStart, cellEnd, cell, field, face, &
+                                       qPt, basis, comp, cidx
+ PetscReal                          :: detFAvg
+ PetscReal                          :: BMat(dimPlex*dimPlex,cellDof)
+ PetscObject                        :: dummy
+ PetscInt                           :: bcSize
+ IS                                 :: bcPoints
+ PetscErrorCode                     :: ierr
+ 
+ pV0 => v0
+ pCellJ => cellJ
+ pInvcellJ => invcellJ
+ call DMGetDefaultSection(dm_local,section,ierr); CHKERRQ(ierr)
+ call DMGetDS(dm_local,prob,ierr); CHKERRQ(ierr)
+ call PetscDSGetTabulation(prob,0,basisField,basisFieldDer,ierr)
+ CHKERRQ(ierr)
+ call DMPlexGetHeightStratum(dm_local,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
+ call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
+ call VecWAXPY(x_local,1.0,xx_local,solution_local,ierr); CHKERRQ(ierr)
+ do field = 1, dimPlex; do face = 1, mesh_Nboundaries
+   if (params%fieldBC%componentBC(field)%Mask(face)) then
+     call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,ierr)
+     if (bcSize > 0) then
+       call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr)
+       CHKERRQ(ierr)
+       call utilities_projectBCValues(x_local,section,0,field-1,bcPoints, &
+                                      0.0,params%fieldBC%componentBC(field)%Value(face),params%timeinc)
+       call ISDestroy(bcPoints,ierr); CHKERRQ(ierr)
+     endif
+   endif
+ enddo; enddo
+
+!--------------------------------------------------------------------------------------------------
+! evaluate field derivatives
+ do cell = cellStart, cellEnd-1                                                                     !< loop over all elements 
+   call DMPlexVecGetClosure(dm_local,section,x_local,cell,x_scal,ierr)                              !< get Dofs belonging to element
+   CHKERRQ(ierr)
+   call  DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,ierr) 
+   CHKERRQ(ierr)
+   IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex])
+   do qPt = 0, nQuadrature-1
+     BMat = 0.0
+     do basis = 0, nBasis-1
+       do comp = 0, dimPlex-1
+         cidx = basis*dimPlex+comp
+         BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = &
+           matmul(IcellJMat,basisFieldDer((qPt*nBasis*dimPlex+cidx  )*dimPlex+1: &
+                                          (qPt*nBasis*dimPlex+cidx+1)*dimPlex   ))
+       enddo
+     enddo
+     materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1) = &
+       reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1])
+   enddo    
+   if (BBarStabilisation) then
+     detFAvg = math_det33(sum(materialpoint_F(1:3,1:3,1:nQuadrature,cell+1),dim=3)/real(nQuadrature))
+     do qPt = 1, nQuadrature
+       materialpoint_F(1:dimPlex,1:dimPlex,qPt,cell+1) = &
+         materialpoint_F(1:dimPlex,1:dimPlex,qPt,cell+1)* &
+         (detFAvg/math_det33(materialpoint_F(1:3,1:3,qPt,cell+1)))**(1.0/real(dimPlex))
+         
+     enddo  
+   endif
+   call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,ierr)
+   CHKERRQ(ierr) 
+ enddo 
+
+!--------------------------------------------------------------------------------------------------
+! evaluate constitutive response
+ call Utilities_constitutiveResponse(params%timeinc,P_av,ForwardData)
+ call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
+ ForwardData = .false.
+
+!--------------------------------------------------------------------------------------------------
+! integrating residual 
+ do cell = cellStart, cellEnd-1                                                                     !< loop over all elements 
+   call DMPlexVecGetClosure(dm_local,section,x_local,cell,x_scal,ierr)                              !< get Dofs belonging to element
+   CHKERRQ(ierr)
+   call  DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,ierr) 
+   CHKERRQ(ierr)
+   IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex])
+   f_scal = 0.0
+   do qPt = 0, nQuadrature-1
+     BMat = 0.0
+     do basis = 0, nBasis-1
+       do comp = 0, dimPlex-1
+         cidx = basis*dimPlex+comp
+         BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = &
+           matmul(IcellJMat,basisFieldDer((qPt*nBasis*dimPlex+cidx  )*dimPlex+1: &
+                                          (qPt*nBasis*dimPlex+cidx+1)*dimPlex   ))
+       enddo
+     enddo
+     f_scal = f_scal + &
+              matmul(transpose(BMat), &
+                     reshape(transpose(materialpoint_P(1:dimPlex,1:dimPlex,qPt+1,cell+1)), &
+                             shape=[dimPlex*dimPlex]))*qWeights(qPt+1)          
+   enddo    
+   f_scal = f_scal*abs(detJ)  
+   pf_scal => f_scal
+   call DMPlexVecSetClosure(dm_local,section,f_local,cell,pf_scal,ADD_VALUES,ierr)
+   CHKERRQ(ierr) 
+   call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,ierr)
+   CHKERRQ(ierr) 
+ enddo
+ call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
+ 
+end subroutine FEM_mech_formResidual
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief forms the FEM stiffness matrix
+!--------------------------------------------------------------------------------------------------
+subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
+ use numerics, only: &
+   BBarStabilisation
+ use homogenization, only: &
+   materialpoint_dPdF, &
+   materialpoint_F
+ use math, only: &
+   math_inv33, &
+   math_identity2nd, &
+   math_det33
+ use FEM_utilities, only: &
+   utilities_projectBCValues
+
+ implicit none
+
+ DM                                   :: dm_local
+ PetscDS                              :: prob
+ Vec                                  :: x_local, xx_local
+ Mat                                  :: Jac_pre, Jac
+ PetscSection                         :: section, gSection
+ PetscReal                            :: detJ, IcellJMat(dimPlex,dimPlex)
+ PetscReal,                    target :: v0(dimPlex), cellJ(dimPlex*dimPlex), &
+                                                      invcellJ(dimPlex*dimPlex)
+ PetscReal,                   pointer :: pV0(:), pCellJ(:), pInvcellJ(:)
+ PetscReal,   dimension(:),   pointer :: basisField, basisFieldDer
+ PetscInt                             :: cellStart, cellEnd, cell, field, face, &
+                                         qPt, basis, comp, cidx
+ PetscScalar,                  target :: K_e   (cellDof,cellDof), &
+                                         K_eA  (cellDof,cellDof), &  
+                                         K_eB  (cellDof,cellDof), &  
+                                         K_eVec(cellDof*cellDof)
+ PetscReal                            :: BMat   (dimPlex*dimPlex,cellDof), &
+                                         BMatAvg(dimPlex*dimPlex,cellDof), &
+                                         MatA   (dimPlex*dimPlex,cellDof), &
+                                         MatB   (1              ,cellDof)
+ PetscScalar, dimension(:),   pointer :: pK_e, x_scal
+ PetscReal,   dimension(3,3)          :: F = math_I3, FAvg, FInv
+ PetscObject                          :: dummy
+ PetscInt                             :: bcSize
+ IS                                   :: bcPoints
+ PetscErrorCode                       :: ierr
+ 
+ pV0 => v0
+ pCellJ => cellJ
+ pInvcellJ => invcellJ
+ call MatSetOption(Jac,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE,ierr); CHKERRQ(ierr)
+ call MatSetOption(Jac,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE,ierr); CHKERRQ(ierr)
+ call MatZeroEntries(Jac,ierr); CHKERRQ(ierr)
+ call DMGetDS(dm_local,prob,ierr); CHKERRQ(ierr)
+ call PetscDSGetTabulation(prob,0,basisField,basisFieldDer,ierr)
+ call DMGetDefaultSection(dm_local,section,ierr); CHKERRQ(ierr)
+ call DMGetDefaultGlobalSection(dm_local,gSection,ierr); CHKERRQ(ierr)
+ 
+ call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
+ call VecWAXPY(x_local,1.0,xx_local,solution_local,ierr); CHKERRQ(ierr)
+ do field = 1, dimPlex; do face = 1, mesh_Nboundaries
+   if (params%fieldBC%componentBC(field)%Mask(face)) then
+     call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,ierr)
+     if (bcSize > 0) then
+       call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr)
+       CHKERRQ(ierr)
+       call utilities_projectBCValues(x_local,section,0,field-1,bcPoints, &
+                                      0.0,params%fieldBC%componentBC(field)%Value(face),params%timeinc)
+       call ISDestroy(bcPoints,ierr); CHKERRQ(ierr)
+     endif
+   endif
+ enddo; enddo
+ call DMPlexGetHeightStratum(dm_local,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
+ do cell = cellStart, cellEnd-1                                                                     !< loop over all elements 
+   call DMPlexVecGetClosure(dm_local,section,x_local,cell,x_scal,ierr)                              !< get Dofs belonging to element
+   CHKERRQ(ierr)
+   call  DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,ierr) 
+   CHKERRQ(ierr)
+   IcellJMat = reshape(pInvcellJ, shape = [dimPlex,dimPlex])
+   K_eA = 0.0
+   K_eB = 0.0
+   MatB = 0.0
+   FAvg = 0.0
+   BMatAvg = 0.0
+   do qPt = 0, nQuadrature-1
+     BMat = 0.0
+     do basis = 0, nBasis-1
+       do comp = 0, dimPlex-1
+         cidx = basis*dimPlex+comp
+         BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = &
+           matmul(IcellJMat,basisFieldDer((qPt*nBasis*dimPlex+cidx  )*dimPlex+1: &
+                                          (qPt*nBasis*dimPlex+cidx+1)*dimPlex   ))
+       enddo
+     enddo
+     MatA = matmul(reshape(reshape(materialpoint_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,qPt+1,cell+1), &
+                                   shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), &
+                           shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1)        
+     if (BBarStabilisation) then
+       F(1:dimPlex,1:dimPlex) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex])
+       FInv = math_inv33(F)
+       K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0/real(dimPlex))        
+       K_eB = K_eB - &
+              matmul(transpose(matmul(reshape(materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1), &
+                                              shape=[dimPlex*dimPlex,1]), &
+                                      matmul(reshape(FInv(1:dimPlex,1:dimPlex), &
+                                                     shape=[1,dimPlex*dimPlex],order=[2,1]),BMat))),MatA)
+       MatB = MatB + &
+              matmul(reshape(materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1),shape=[1,dimPlex*dimPlex]),MatA)  
+       FAvg = FAvg + F      
+       BMatAvg = BMatAvg + BMat
+     else
+       K_eA = K_eA + matmul(transpose(BMat),MatA)
+     endif  
+   enddo  
+   if (BBarStabilisation) then
+     FInv = math_inv33(FAvg)
+     K_e = K_eA*math_det33(FAvg/real(nQuadrature))**(1.0/real(dimPlex)) + &
+           (matmul(matmul(transpose(BMatAvg), &
+                          reshape(FInv(1:dimPlex,1:dimPlex),shape=[dimPlex*dimPlex,1],order=[2,1])),MatB) + &
+            K_eB)/real(dimPlex)              
+     
+   else
+     K_e = K_eA
+   endif  
+   K_e = K_e + eps*math_identity2nd(cellDof)
+   K_eVec = reshape(K_e, [cellDof*cellDof])*abs(detJ) 
+   pK_e => K_eVec      
+   call DMPlexMatSetClosure(dm_local,section,gSection,Jac,cell,pK_e,ADD_VALUES,ierr)
+   CHKERRQ(ierr)
+   call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,ierr)
+   CHKERRQ(ierr) 
+ enddo 
+ call MatAssemblyBegin(Jac,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr)
+ call MatAssemblyEnd(Jac,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr)
+ call MatAssemblyBegin(Jac_pre,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr)
+ call MatAssemblyEnd(Jac_pre,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr)
+ call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
+ 
+!--------------------------------------------------------------------------------------------------
+! apply boundary conditions 
+ !call DMPlexCreateRigidBody(dm_local,matnull,ierr); CHKERRQ(ierr) MD: linker error
+ call MatSetNullSpace(Jac,matnull,ierr); CHKERRQ(ierr)
+ call MatSetNearNullSpace(Jac,matnull,ierr); CHKERRQ(ierr)
+ call MatNullSpaceDestroy(matnull,ierr); CHKERRQ(ierr)
+
+end subroutine FEM_mech_formJacobian
+
+!--------------------------------------------------------------------------------------------------
+!> @brief forwarding routine
+!--------------------------------------------------------------------------------------------------
+subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC)
+ use FEM_utilities, only: &
+   cutBack
+ use homogenization, only: &
+   materialpoint_F0, &
+   materialpoint_F
+ use FEM_utilities, only: &
+   utilities_projectBCValues
+
+ implicit none
+ type(tFieldBC), intent(in) :: &
+   fieldBC
+ real(pReal), intent(in) :: &
+   timeinc_old, &
+   timeinc 
+ logical, intent(in) :: &
+   guess
+ PetscInt                :: field, face
+ DM                      :: dm_local
+ Vec                     :: x_local
+ PetscSection            :: section
+ PetscInt                :: bcSize
+ IS                      :: bcPoints
+ PetscErrorCode          :: ierr   
+
+!--------------------------------------------------------------------------------------------------
+! forward last inc
+ if (guess .and. .not. cutBack) then 
+   ForwardData = .True.
+   materialpoint_F0 = materialpoint_F
+   call SNESGetDM(mech_snes,dm_local,ierr); CHKERRQ(ierr)                                                 !< retrieve mesh info from mech_snes into dm_local
+   call DMGetDefaultSection(dm_local,section,ierr); CHKERRQ(ierr)
+   call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
+   call VecSet(x_local,0.0,ierr); CHKERRQ(ierr)
+   call DMGlobalToLocalBegin(dm_local,solution,INSERT_VALUES,x_local,ierr)                       !< retrieve my partition of global solution vector
+   CHKERRQ(ierr)
+   call DMGlobalToLocalEnd(dm_local,solution,INSERT_VALUES,x_local,ierr)
+   CHKERRQ(ierr)
+   call VecAXPY(solution_local,1.0,x_local,ierr); CHKERRQ(ierr)
+   do field = 1, dimPlex; do face = 1, mesh_Nboundaries
+     if (fieldBC%componentBC(field)%Mask(face)) then
+       call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,ierr)
+       if (bcSize > 0) then
+         call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr)
+         CHKERRQ(ierr)
+         call utilities_projectBCValues(solution_local,section,0,field-1,bcPoints, &
+                                        0.0,fieldBC%componentBC(field)%Value(face),timeinc_old)
+         call ISDestroy(bcPoints,ierr); CHKERRQ(ierr)
+       endif
+     endif
+   enddo; enddo
+   call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
+
+!--------------------------------------------------------------------------------------------------
+! update rate and forward last inc
+   call VecCopy(solution,solution_rate,ierr); CHKERRQ(ierr)
+   call VecScale(solution_rate,1.0/timeinc_old,ierr); CHKERRQ(ierr)
+ endif
+ call VecCopy(solution_rate,solution,ierr); CHKERRQ(ierr)
+ call VecScale(solution,timeinc,ierr); CHKERRQ(ierr)
+ 
+end subroutine FEM_mech_forward
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief reporting
+!--------------------------------------------------------------------------------------------------
+subroutine FEM_mech_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr)
+ use numerics, only: &
+   err_struct_tolAbs, &
+   err_struct_tolRel
+ use IO, only: &
+   IO_intOut
+ use FEsolving, only: &
+   terminallyIll
+
+ implicit none
+ SNES :: snes_local
+ PetscInt :: PETScIter
+ PetscReal :: xnorm,snorm,fnorm,divTol
+ SNESConvergedReason :: reason
+ PetscObject :: dummy
+ PetscErrorCode :: ierr
+
+!--------------------------------------------------------------------------------------------------
+! report
+ divTol = max(maxval(abs(P_av(1:dimPlex,1:dimPlex)))*err_struct_tolRel,err_struct_tolAbs)
+ call SNESConvergedDefault(snes_local,PETScIter,xnorm,snorm,fnorm/divTol,reason,dummy,ierr)
+ CHKERRQ(ierr)
+ if (terminallyIll) reason = SNES_DIVERGED_FUNCTION_DOMAIN
+ if (worldrank == 0) then 
+   write(6,'(1/,1x,a,a,i0,a,i0,f0.3)') trim(incInfo), &
+                   ' @ Iteration ',PETScIter,' mechanical residual norm = ', &
+                                                   int(fnorm/divTol),fnorm/divTol-int(fnorm/divTol)
+   write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',&
+                                                       transpose(P_av)*1.e-6_pReal
+   flush(6) 
+ endif
+ 
+end subroutine FEM_mech_converged
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief destroy routine
+!--------------------------------------------------------------------------------------------------
+subroutine FEM_mech_destroy()
+
+ implicit none
+ PetscErrorCode :: ierr
+
+ call VecDestroy(solution,ierr); CHKERRQ(ierr)
+ call VecDestroy(solution_rate,ierr); CHKERRQ(ierr)
+ call SNESDestroy(mech_snes,ierr); CHKERRQ(ierr)
+
+end subroutine FEM_mech_destroy
+
+end module FEM_mech

src/FEM_utilities.f90

@@ -0,0 +1,751 @@
+!--------------------------------------------------------------------------------------------------
+!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
+!> @brief Utilities used by the FEM solver
+!--------------------------------------------------------------------------------------------------
+module FEM_utilities
+#include <petsc/finclude/petsc.h>
+ use prec, only: pReal, pInt
+
+use PETScdmda
+use PETScis
+
+ implicit none
+ private
+!--------------------------------------------------------------------------------------------------
+! 
+ logical,       public             :: cutBack = .false.                                             !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
+ integer(pInt), public, parameter  :: maxFields = 6_pInt
+ integer(pInt), public             :: nActiveFields = 0_pInt
+ 
+!--------------------------------------------------------------------------------------------------
+! grid related information information
+ real(pReal),   public                :: wgt                                                        !< weighting factor 1/Nelems
+ real(pReal),   public                :: wgtDof                                                     !< weighting factor 1/Nelems
+ real(pReal),   public                :: C_volAvg(3,3,3,3)
+ 
+!--------------------------------------------------------------------------------------------------
+! output data
+ PetscViewer,                    public :: resUnit
+ Vec,                            public :: coordinatesVec
+ Vec,               allocatable, public :: homogenizationResultsVec(:), &
+                                           crystalliteResultsVec(:,:), &
+                                           phaseResultsVec(:,:) 
+
+!--------------------------------------------------------------------------------------------------
+! field labels information
+ character(len=*),                         parameter,            public :: &
+   FIELD_MECH_label     = 'mechanical', &
+   FIELD_THERMAL_label  = 'thermal', &
+   FIELD_DAMAGE_label   = 'damage', &
+   FIELD_SOLUTE_label   = 'solute', &
+   FIELD_MGTWIN_label   = 'mgtwin'
+ 
+ enum, bind(c)
+   enumerator :: FIELD_UNDEFINED_ID, &
+                 FIELD_MECH_ID, &
+                 FIELD_THERMAL_ID, &
+                 FIELD_DAMAGE_ID, &
+                 FIELD_SOLUTE_ID, &
+                 FIELD_MGTWIN_ID
+ end enum
+ enum, bind(c)
+   enumerator :: COMPONENT_UNDEFINED_ID, &
+                 COMPONENT_MECH_X_ID, &
+                 COMPONENT_MECH_Y_ID, &
+                 COMPONENT_MECH_Z_ID, &
+                 COMPONENT_THERMAL_T_ID, &
+                 COMPONENT_DAMAGE_PHI_ID, &
+                 COMPONENT_SOLUTE_CV_ID, &
+                 COMPONENT_SOLUTE_CVPOT_ID, &
+                 COMPONENT_SOLUTE_CH_ID, &
+                 COMPONENT_SOLUTE_CHPOT_ID, &
+                 COMPONENT_SOLUTE_CVaH_ID, &
+                 COMPONENT_SOLUTE_CVaHPOT_ID, &
+                 COMPONENT_MGTWIN_PHI_ID
+ end enum
+ 
+!--------------------------------------------------------------------------------------------------
+! variables controlling debugging
+ logical, private :: &
+   debugGeneral, &                                                                                  !< general debugging of FEM solver
+   debugRotation, &                                                                                 !< also printing out results in lab frame
+   debugPETSc                                                                                       !< use some in debug defined options for more verbose PETSc solution
+
+!--------------------------------------------------------------------------------------------------
+! derived types
+ type, public :: tSolutionState                                                                     !< return type of solution from FEM solver variants
+   logical       :: converged         = .true.   
+   logical       :: stagConverged     = .true.   
+   logical       :: regrid            = .false.   
+   integer(pInt) :: iterationsNeeded  = 0_pInt
+ end type tSolutionState
+
+ type, public :: tComponentBC
+   integer(kind(COMPONENT_UNDEFINED_ID))          :: ID
+   real(pReal),                       allocatable :: Value(:)
+   logical,                           allocatable :: Mask(:)   
+ end type tComponentBC
+
+ type, public :: tFieldBC
+   integer(kind(FIELD_UNDEFINED_ID))              :: ID
+   integer(pInt)                                  :: nComponents = 0_pInt
+   type(tComponentBC),                allocatable :: componentBC(:)
+ end type tFieldBC
+
+ type, public :: tLoadCase
+   real(pReal)                       :: time                   = 0.0_pReal                               !< length of increment
+   integer(pInt)                     :: incs                   = 0_pInt, &                               !< number of increments
+                                        outputfrequency        = 1_pInt, &                               !< frequency of result writes
+                                        restartfrequency       = 0_pInt, &                               !< frequency of restart writes
+                                        logscale               = 0_pInt                                  !< linear/logarithmic time inc flag
+   logical                           :: followFormerTrajectory = .true.                                  !< follow trajectory of former loadcase
+   integer(pInt),        allocatable :: faceID(:)
+   type(tFieldBC),       allocatable :: fieldBC(:)
+ end type tLoadCase
+
+ type, public :: tFEMInterpolation
+   integer(pInt)                                :: n                             
+   real(pReal),   dimension(:,:)  , allocatable :: shapeFunc, shapeDerivReal, geomShapeDerivIso
+   real(pReal),   dimension(:,:,:), allocatable :: shapeDerivIso
+ end type tFEMInterpolation
+ 
+ type, public :: tQuadrature
+   integer(pInt)                            :: n
+   real(pReal), dimension(:)  , allocatable :: Weights
+   real(pReal), dimension(:,:), allocatable :: Points
+ end type tQuadrature   
+ 
+ public :: &
+   utilities_init, &
+   utilities_constitutiveResponse, &
+   utilities_indexBoundaryDofs, &
+   utilities_projectBCValues, &
+   utilities_indexActiveSet, &
+   utilities_destroy, &
+   FIELD_MECH_ID, &
+   FIELD_THERMAL_ID, &
+   FIELD_DAMAGE_ID, &
+   FIELD_SOLUTE_ID, &
+   FIELD_MGTWIN_ID, &
+   COMPONENT_MECH_X_ID, &
+   COMPONENT_MECH_Y_ID, &
+   COMPONENT_MECH_Z_ID, &
+   COMPONENT_THERMAL_T_ID, &
+   COMPONENT_DAMAGE_PHI_ID, &
+   COMPONENT_SOLUTE_CV_ID, &
+   COMPONENT_SOLUTE_CVPOT_ID, &
+   COMPONENT_SOLUTE_CH_ID, &
+   COMPONENT_SOLUTE_CHPOT_ID, &
+   COMPONENT_SOLUTE_CVaH_ID, &
+   COMPONENT_SOLUTE_CVaHPOT_ID, &
+   COMPONENT_MGTWIN_PHI_ID
+
+ external :: &
+   PetscOptionsInsertString, &
+   PetscObjectSetName, &
+   DMPlexGetHeightStratum, &
+   DMGetLabelIdIS, &
+   DMPlexGetChart, &
+   DMPlexLabelComplete, &
+   PetscViewerHDF5Open, &
+   PetscViewerHDF5PushGroup, &
+   PetscViewerHDF5PopGroup, &
+   PetscViewerDestroy
+
+contains 
+
+!--------------------------------------------------------------------------------------------------
+!> @brief allocates all neccessary fields, sets debug flags
+!--------------------------------------------------------------------------------------------------
+subroutine utilities_init()
+ use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment)
+ use DAMASK_interface, only: &
+  getSolverJobName
+ use IO, only: &
+   IO_error, &
+   IO_warning, &
+   IO_timeStamp, &
+   IO_open_file
+ use numerics, only: &                        
+   integrationOrder, &
+   worldsize, & 
+   worldrank, &
+   petsc_defaultOptions, &
+   petsc_options
+ use debug, only: &
+   debug_level, &
+   debug_SPECTRAL, &
+   debug_LEVELBASIC, &
+   debug_SPECTRALPETSC, &
+   debug_SPECTRALROTATION
+ use debug, only: &
+   PETSCDEBUG
+ use math                                                                                           ! must use the whole module for use of FFTW
+ use mesh, only: &
+   mesh_NcpElemsGlobal, &
+   mesh_maxNips, &
+   geomMesh
+ use material, only: &
+   material_homog
+
+ implicit none
+
+ character(len=1024)                :: petsc_optionsPhysics, grainStr
+ integer(pInt)                      :: dimPlex
+ integer(pInt)                      :: headerID = 205_pInt
+ PetscInt,    dimension(:), pointer :: points
+ PetscInt,    allocatable           :: nEntities(:), nOutputCells(:), nOutputNodes(:), mappingCells(:)
+ PetscInt                           :: cellStart, cellEnd, cell, ip, dim, ctr, qPt
+ PetscInt,              allocatable :: connectivity(:,:)
+ Vec                                :: connectivityVec
+ PetscErrorCode                     :: ierr
+
+ write(6,'(/,a)')   ' <<<+-  DAMASK_FEM_utilities init  -+>>>'
+ write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
+#include "compilation_info.f90"
+ 
+!--------------------------------------------------------------------------------------------------
+! set debugging parameters
+ debugGeneral    = iand(debug_level(debug_SPECTRAL),debug_LEVELBASIC)         /= 0
+ debugRotation   = iand(debug_level(debug_SPECTRAL),debug_SPECTRALROTATION)   /= 0
+ debugPETSc      = iand(debug_level(debug_SPECTRAL),debug_SPECTRALPETSC)      /= 0
+ if(debugPETSc) write(6,'(3(/,a),/)') &
+                ' Initializing PETSc with debug options: ', &
+                trim(PETScDebug), &
+                ' add more using the PETSc_Options keyword in numerics.config '
+ flush(6)
+ call PetscOptionsClear(PETSC_NULL_OPTIONS,ierr)
+ CHKERRQ(ierr)
+ if(debugPETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),ierr)
+ CHKERRQ(ierr)
+ call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_defaultOptions),ierr)
+ call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr)
+ CHKERRQ(ierr)
+ !write(petsc_optionsPhysics,'(a,i0)') '-mechFE_petscspace_order '   , structOrder
+ call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsPhysics),ierr)
+ CHKERRQ(ierr)
+ 
+ wgt = 1.0/real(mesh_maxNips*mesh_NcpElemsGlobal,pReal)
+
+ call PetscViewerHDF5Open(PETSC_COMM_WORLD, trim(getSolverJobName())//'.h5', &
+                          FILE_MODE_WRITE, resUnit, ierr); CHKERRQ(ierr)
+ call PetscViewerHDF5PushGroup(resUnit, '/', ierr); CHKERRQ(ierr)  
+ call DMGetDimension(geomMesh,dimPlex,ierr); CHKERRQ(ierr) 
+ allocate(nEntities(dimPlex+1), source=0)
+ allocate(nOutputNodes(worldsize), source = 0)
+ allocate(nOutputCells(worldsize), source = 0)
+ do dim = 0, dimPlex
+   call DMGetStratumSize(geomMesh,'depth',dim,nEntities(dim+1),ierr)
+   CHKERRQ(ierr)
+ enddo
+ select case (integrationOrder)
+   case(1_pInt)
+     nOutputNodes(worldrank+1) = nEntities(1)
+   case(2_pInt)  
+     nOutputNodes(worldrank+1) = sum(nEntities)
+   case default
+     nOutputNodes(worldrank+1) = mesh_maxNips*nEntities(dimPlex+1)
+ end select    
+ nOutputCells(worldrank+1) = count(material_homog > 0_pInt)
+ call MPI_Allreduce(MPI_IN_PLACE,nOutputNodes,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)
+ call MPI_Allreduce(MPI_IN_PLACE,nOutputCells,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)
+ if (worldrank == 0_pInt) then
+   open(unit=headerID, file=trim(getSolverJobName())//'.header', &
+        form='FORMATTED', status='REPLACE')
+   write(headerID, '(a,i0)') 'dimension : ',   dimPlex                                             
+   write(headerID, '(a,i0)') 'number of nodes : ', sum(nOutputNodes)
+   write(headerID, '(a,i0)') 'number of cells : ', sum(nOutputCells)
+ endif 
+
+ allocate(connectivity(2**dimPlex,nOutputCells(worldrank+1)))
+ call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr)
+ CHKERRQ(ierr)
+ ctr = 0
+ select case (integrationOrder)
+   case(1_pInt)
+     do cell = cellStart, cellEnd-1                                                                     !< loop over all elements 
+       call DMPlexGetTransitiveClosure(geomMesh,cell,PETSC_TRUE,points,ierr)
+       CHKERRQ(ierr)
+       if (dimPlex == 2) then
+         connectivity(:,ctr+1) =  [points( 9), points(11), points(13), points(13)]  - nEntities(dimPlex+1)
+         ctr = ctr + 1
+       else
+         connectivity(:,ctr+1) =  [points(23), points(25), points(27), points(27), &
+                                   points(29), points(29), points(29), points(29)] - nEntities(dimPlex+1)
+         ctr = ctr + 1
+       endif 
+     enddo 
+   
+   case(2_pInt)  
+     do cell = cellStart, cellEnd-1                                                                     !< loop over all elements 
+       call DMPlexGetTransitiveClosure(geomMesh,cell,PETSC_TRUE,points,ierr)
+       CHKERRQ(ierr)
+       if (dimPlex == 2) then
+         connectivity(:,ctr+1) = [points(9 ), points(3), points(1), points(7)]
+         connectivity(:,ctr+2) = [points(11), points(5), points(1), points(3)]
+         connectivity(:,ctr+3) = [points(13), points(7), points(1), points(5)]
+         ctr = ctr + 3
+       else
+         connectivity(:,ctr+1) = [points(23), points(11), points(3), points(15), points(17), points(5), points(1), points(7)]
+         connectivity(:,ctr+2) = [points(25), points(13), points(3), points(11), points(19), points(9), points(1), points(5)]
+         connectivity(:,ctr+3) = [points(27), points(15), points(3), points(13), points(21), points(7), points(1), points(9)]
+         connectivity(:,ctr+4) = [points(29), points(17), points(7), points(21), points(19), points(5), points(1), points(9)]
+         ctr = ctr + 4_pInt
+       endif  
+     enddo 
+   
+   case default
+     do cell = cellStart, cellEnd-1; do ip = 0, mesh_maxNips-1
+       connectivity(:,ctr+1) = cell*mesh_maxNips + ip
+       ctr = ctr + 1
+     enddo; enddo        
+   
+ end select 
+ connectivity = connectivity + sum(nOutputNodes(1:worldrank))
+
+ call VecCreateMPI(PETSC_COMM_WORLD,dimPlex*nOutputNodes(worldrank+1),dimPlex*sum(nOutputNodes), &
+                   coordinatesVec,ierr);CHKERRQ(ierr)
+ call PetscObjectSetName(coordinatesVec, 'NodalCoordinates',ierr)
+ call VecSetFromOptions(coordinatesVec, ierr); CHKERRQ(ierr)
+ 
+ !allocate(mappingCells(worldsize), source = 0)
+ !do homog = 1, material_Nhomogenization
+ !  mappingCells = 0_pInt; mappingCells(worldrank+1) = homogOutput(homog)%sizeIpCells
+ !  call MPI_Allreduce(MPI_IN_PLACE,mappingCells,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)
+ !  call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1),sum(mappingCells), &
+ !                    homogenizationResultsVec(homog),ierr);CHKERRQ(ierr)
+ !  if (sum(mappingCells) > 0) then
+ !    call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1)*2**dimPlex,sum(mappingCells)*2**dimPlex, &
+ !                      connectivityVec,ierr);CHKERRQ(ierr)
+ !    call PetscObjectSetName(connectivityVec,'mapping_'//trim(homogenization_name(homog)),ierr)
+ !    CHKERRQ(ierr)
+ !    call VecGetArrayF90(connectivityVec,results,ierr); CHKERRQ(ierr) 
+ !    results = 0.0_pReal; ctr = 1_pInt
+ !    do cell = cellStart, cellEnd-1; do qPt = 1, mesh_maxNips
+ !      if (material_homog(qPt,cell+1) == homog) then
+ !        results(ctr:ctr+2**dimPlex-1) = real(reshape(connectivity(1:2**dimPlex,mesh_maxNips*cell+qPt), &
+ !                                                     shape=[2**dimPlex]))
+ !        ctr = ctr + 2**dimPlex
+ !      endif  
+ !    enddo; enddo
+ !    call VecRestoreArrayF90(connectivityVec, results, ierr); CHKERRQ(ierr)
+ !    call VecAssemblyBegin(connectivityVec, ierr); CHKERRQ(ierr)
+ !    call VecAssemblyEnd  (connectivityVec, ierr); CHKERRQ(ierr)
+ !    call VecView(connectivityVec, resUnit, ierr); CHKERRQ(ierr)
+ !    call VecDestroy(connectivityVec, ierr); CHKERRQ(ierr)
+ !  endif
+ !enddo  
+ !do cryst = 1, material_Ncrystallite; do grain = 1, homogenization_maxNgrains
+ !  mappingCells              = 0_pInt
+ !  mappingCells(worldrank+1) = crystalliteOutput(cryst,grain)%sizeIpCells
+ !  call MPI_Allreduce(MPI_IN_PLACE,mappingCells,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)
+ !  call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1),sum(mappingCells), &
+ !                    crystalliteResultsVec(cryst,grain),ierr);CHKERRQ(ierr)
+ !  if (sum(mappingCells) > 0) then
+ !    call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1)*2**dimPlex,sum(mappingCells)*2**dimPlex, &
+ !                      connectivityVec,ierr);CHKERRQ(ierr)
+ !    write(grainStr,'(a,i0)') 'Grain',grain
+ !    call PetscObjectSetName(connectivityVec,'mapping_'// &
+ !                                            trim(crystallite_name(cryst))//'_'// &
+ !                                            trim(grainStr),ierr)
+ !    CHKERRQ(ierr)
+ !    call VecGetArrayF90(connectivityVec, results, ierr); CHKERRQ(ierr) 
+ !    results = 0.0_pReal; ctr = 1_pInt
+ !    do cell = cellStart, cellEnd-1; do qPt = 1, mesh_maxNips
+ !      if (homogenization_Ngrains    (mesh_element(3,cell+1)) >= grain .and. &
+ !          microstructure_crystallite(mesh_element(4,cell+1)) == cryst) then
+ !        results(ctr:ctr+2**dimPlex-1) = real(reshape(connectivity(1:2**dimPlex,mesh_maxNips*cell+qPt), &
+ !                                                     shape=[2**dimPlex]))
+ !        ctr = ctr + 2**dimPlex
+ !      endif  
+ !    enddo; enddo
+ !    call VecRestoreArrayF90(connectivityVec, results, ierr); CHKERRQ(ierr)
+ !    call VecAssemblyBegin(connectivityVec, ierr); CHKERRQ(ierr)
+ !    call VecAssemblyEnd  (connectivityVec, ierr); CHKERRQ(ierr)
+ !    call VecView(connectivityVec, resUnit, ierr); CHKERRQ(ierr)
+ !    call VecDestroy(connectivityVec, ierr); CHKERRQ(ierr)
+ !  endif
+ !enddo; enddo
+ !do phase = 1, material_Nphase; do grain = 1, homogenization_maxNgrains
+ !  mappingCells              = 0_pInt
+ !  mappingCells(worldrank+1) = phaseOutput(phase,grain)%sizeIpCells
+ !  call MPI_Allreduce(MPI_IN_PLACE,mappingCells,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)
+ !  call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1),sum(mappingCells), &
+ !                    phaseResultsVec(phase,grain),ierr);CHKERRQ(ierr)
+ !  if (sum(mappingCells) > 0) then
+ !    call VecCreateMPI(PETSC_COMM_WORLD,mappingCells(worldrank+1)*2**dimPlex,sum(mappingCells)*2**dimPlex, &
+ !                      connectivityVec,ierr);CHKERRQ(ierr)
+ !    write(grainStr,'(a,i0)') 'Grain',grain
+ !    call PetscObjectSetName(connectivityVec,&
+ !                            'mapping_'//trim(phase_name(phase))//'_'// &
+ !                            trim(grainStr),ierr)
+ !    CHKERRQ(ierr)
+ !    call VecGetArrayF90(connectivityVec, results, ierr)
+ !    CHKERRQ(ierr) 
+ !    results = 0.0_pReal; ctr = 1_pInt
+ !    do cell = cellStart, cellEnd-1; do qPt = 1, mesh_maxNips
+ !      if (material_phase(grain,qPt,cell+1) == phase) then
+ !        results(ctr:ctr+2**dimPlex-1) = real(reshape(connectivity(1:2**dimPlex,mesh_maxNips*cell+qPt), &
+ !                                                     shape=[2**dimPlex]))
+ !        ctr = ctr + 2**dimPlex
+ !      endif  
+ !    enddo; enddo
+ !    call VecRestoreArrayF90(connectivityVec, results, ierr)
+ !    CHKERRQ(ierr)
+ !    call VecAssemblyBegin(connectivityVec, ierr);CHKERRQ(ierr)
+ !    call VecAssemblyEnd  (connectivityVec, ierr);CHKERRQ(ierr)
+ !    call VecView(connectivityVec, resUnit, ierr);CHKERRQ(ierr)
+ !    call VecDestroy(connectivityVec, ierr); CHKERRQ(ierr)
+ !  endif
+ !enddo; enddo  
+ !if (worldrank == 0_pInt) then
+ !  do homog = 1, material_Nhomogenization
+ !    call VecGetSize(homogenizationResultsVec(homog),mappingCells(1),ierr)
+ !    CHKERRQ(ierr)
+ !    if (mappingCells(1) > 0) &
+ !      write(headerID, '(a,i0)') 'number of homog_'// &
+ !                                trim(homogenization_name(homog))//'_'// &
+ !                                'cells : ', mappingCells(1)
+ !  enddo
+ !  do cryst = 1, material_Ncrystallite; do grain = 1, homogenization_maxNgrains
+ !    call VecGetSize(crystalliteResultsVec(cryst,grain),mappingCells(1),ierr)
+ !    CHKERRQ(ierr)
+ !    write(grainStr,'(a,i0)') 'Grain',grain
+ !    if (mappingCells(1) > 0) &
+ !      write(headerID, '(a,i0)') 'number of cryst_'// &
+ !                                trim(crystallite_name(cryst))//'_'// &
+ !                                trim(grainStr)//'_'// &
+ !                                'cells : ', mappingCells(1)
+ !  enddo; enddo
+ !  do phase = 1, material_Nphase; do grain = 1, homogenization_maxNgrains
+ !    call VecGetSize(phaseResultsVec(phase,grain),mappingCells(1),ierr)
+ !    CHKERRQ(ierr)
+ !    write(grainStr,'(a,i0)') 'Grain',grain
+ !    if (mappingCells(1) > 0) &
+ !      write(headerID, '(a,i0)') 'number of phase_'// &
+ !                                trim(phase_name(phase))//'_'//trim(grainStr)//'_'// &
+ !                                'cells : ', mappingCells(1)
+ !  enddo; enddo  
+ !  close(headerID) 
+ !endif 
+
+end subroutine utilities_init
+
+!--------------------------------------------------------------------------------------------------
+!> @brief calculates constitutive response
+!--------------------------------------------------------------------------------------------------
+subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
+ use debug, only: &
+   debug_reset, &
+   debug_info
+ use numerics, only: &
+   worldrank
+ use math, only: &
+   math_transpose33, &
+   math_rotate_forward33, &
+   math_det33
+ use FEsolving, only: &
+   restartWrite
+ use homogenization, only: &
+   materialpoint_F0, &
+   materialpoint_F, &
+   materialpoint_P, &
+   materialpoint_dPdF, &
+   materialpoint_stressAndItsTangent
+ use mesh, only: &
+   mesh_NcpElems  
+ 
+ implicit none
+ real(pReal), intent(in)                                         :: timeinc                         !< loading time
+ logical,     intent(in)                                         :: forwardData                     !< age results
+ 
+ real(pReal),intent(out), dimension(3,3)                         :: P_av                            !< average PK stress
+ 
+ logical :: &
+   age
+
+ integer(pInt) :: &
+   j
+ real(pReal)   :: defgradDetMin, defgradDetMax, defgradDet
+ PetscErrorCode :: ierr
+
+ write(6,'(/,a)') ' ... evaluating constitutive response ......................................'
+
+ age = .False.
+ if (forwardData) then                                                                              ! aging results
+   age = .True.
+ endif
+ if (cutBack) then                                                                                  ! restore saved variables
+   age = .False.
+ endif
+ call debug_reset()
+
+!--------------------------------------------------------------------------------------------------
+! calculate bounds of det(F) and report
+ if(debugGeneral) then
+   defgradDetMax = -huge(1.0_pReal)
+   defgradDetMin = +huge(1.0_pReal)
+   do j = 1_pInt, mesh_NcpElems
+     defgradDet = math_det33(materialpoint_F(1:3,1:3,1,j))
+     defgradDetMax = max(defgradDetMax,defgradDet)
+     defgradDetMin = min(defgradDetMin,defgradDet) 
+   end do
+   write(6,'(a,1x,es11.4)') ' max determinant of deformation =', defgradDetMax
+   write(6,'(a,1x,es11.4)') ' min determinant of deformation =', defgradDetMin
+   flush(6)
+ endif
+  
+ call materialpoint_stressAndItsTangent(.true.,timeinc)                                             ! calculate P field
+
+ call debug_info()
+ 
+ restartWrite = .false.                                                                             ! reset restartWrite status
+ cutBack = .false.                                                                                  ! reset cutBack status
+ 
+ P_av = sum(sum(materialpoint_P,dim=4),dim=3) * wgt                                                    ! average of P 
+ C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) * wgt
+ call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
+ call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD, ierr)
+
+end subroutine utilities_constitutiveResponse
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Create index sets of boundary dofs (in local and global numbering)
+!--------------------------------------------------------------------------------------------------
+subroutine utilities_indexBoundaryDofs(dm_local,nFaceSets,numFields,local2global,section,localIS,globalIS)
+
+ implicit none
+ 
+ DM                                    :: dm_local
+ ISLocalToGlobalMapping                :: local2global
+ PetscSection                          :: section
+ PetscInt                              :: nFaceSets, numFields, nDof
+ IS,    dimension(nFaceSets,numFields) :: localIS, globalIS
+ PetscInt                              :: field, faceSet, point, dof, offset
+ PetscInt                              :: localSize, storageSize, ISSize
+ PetscInt, dimension(:)  , allocatable :: localIndices
+ IS                                    :: faceSetIS, BC_IS, dummyIS
+ PetscInt, dimension(:)  ,     pointer :: pFaceSets, pBCvertex, pBCvertexlc
+ DMLabel                               :: BCLabel
+ PetscErrorCode                        :: ierr
+
+ call DMGetLabel(dm_local,'Face Sets',BCLabel,ierr); CHKERRQ(ierr)
+ call DMPlexLabelComplete(dm_local,BCLabel,ierr); CHKERRQ(ierr)
+ call PetscSectionGetStorageSize(section,storageSize,ierr); CHKERRQ(ierr)
+ call DMGetLabelIdIS(dm_local,'Face Sets',faceSetIS,ierr); CHKERRQ(ierr)
+ call ISGetIndicesF90(faceSetIS,pFaceSets,ierr); CHKERRQ(ierr)
+ allocate(localIndices (storageSize))
+ do faceSet = 1, nFaceSets
+   call DMGetStratumSize(dm_local,'Face Sets',pFaceSets(faceSet),ISSize,ierr)
+   CHKERRQ(ierr)
+   call DMGetStratumIS(dm_local,'Face Sets',pFaceSets(faceSet),BC_IS,ierr)
+   CHKERRQ(ierr)
+   if (ISSize > 0) call ISGetIndicesF90(BC_IS,pBCvertex,ierr)
+   do field = 1, numFields
+     localSize  = 0
+     do point = 1, ISSize
+       call PetscSectionGetFieldDof(section,pBCvertex(point),field-1,nDof,ierr)
+       CHKERRQ(ierr)
+       call PetscSectionGetFieldOffset(section,pBCvertex(point),field-1,offset,ierr)
+       CHKERRQ(ierr)
+       do dof = 1, nDof
+         localSize = localSize + 1
+         localIndices(localSize) = offset + dof - 1
+       enddo  
+     enddo
+     call ISCreateGeneral(PETSC_COMM_SELF,localSize,localIndices,PETSC_COPY_VALUES, &
+                          localIS(faceSet,field),ierr)
+     CHKERRQ(ierr)
+     call ISLocalToGlobalMappingApplyIS(local2global,localIS(faceSet,field), &
+                                        globalIS(faceSet,field),ierr)
+     CHKERRQ(ierr)
+   enddo
+   if (ISSize > 0) call ISRestoreIndicesF90(BC_IS,pBCvertex,ierr)
+   call ISDestroy(BC_IS,ierr); CHKERRQ(ierr)
+ enddo  
+ call ISRestoreIndicesF90(faceSetIS,pFaceSets,ierr); CHKERRQ(ierr)
+ call ISDestroy(faceSetIS,ierr); CHKERRQ(ierr)
+ 
+ do faceSet = 1, nFaceSets; do field = 1, numFields
+   call ISGetSize(globalIS(faceSet,field),ISSize,ierr); CHKERRQ(ierr)
+   if (ISSize > 0) then
+     call ISGetIndicesF90(localIS(faceSet,field),pBCvertexlc,ierr); CHKERRQ(ierr)
+     call ISGetIndicesF90(globalIS(faceSet,field),pBCvertex,ierr); CHKERRQ(ierr)
+   endif  
+   localSize = 0
+   do point = 1, ISSize
+     if (pBCvertex(point) >= 0) then
+       localSize = localSize + 1
+       localIndices(localSize) = pBCvertexlc(point)
+     endif  
+   enddo
+   if (ISSize > 0) then
+     call ISRestoreIndicesF90(localIS(faceSet,field),pBCvertexlc,ierr); CHKERRQ(ierr)
+     call ISRestoreIndicesF90(globalIS(faceSet,field),pBCvertex,ierr); CHKERRQ(ierr)
+   endif  
+   call ISDestroy(globalIS(faceSet,field),ierr); CHKERRQ(ierr) 
+   call ISCreateGeneral(PETSC_COMM_SELF,localSize,localIndices,PETSC_COPY_VALUES, &
+                        globalIS(faceSet,field),ierr) 
+   CHKERRQ(ierr)
+   if (ISSize > 0) then
+     call ISDuplicate(localIS(faceSet,field),dummyIS,ierr); CHKERRQ(ierr)
+     call ISDestroy(localIS(faceSet,field),ierr); CHKERRQ(ierr)
+     call ISDifference(dummyIS,globalIS(faceSet,field),localIS(faceSet,field),ierr)
+     CHKERRQ(ierr)
+     call ISDestroy(dummyIS,ierr); CHKERRQ(ierr)
+   endif
+ enddo; enddo                       
+ deallocate(localIndices)
+ 
+end subroutine utilities_indexBoundaryDofs
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Project BC values to local vector
+!--------------------------------------------------------------------------------------------------
+subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCValue,BCDotValue,timeinc)
+
+ implicit none
+ 
+ Vec                  :: localVec
+ PetscInt             :: field, comp, nBcPoints, point, dof, numDof, numComp, offset
+ PetscSection         :: section
+ IS                   :: bcPointsIS
+ PetscInt,    pointer :: bcPoints(:)
+ PetscScalar, pointer :: localArray(:)
+ PetscScalar          :: BCValue,BCDotValue,timeinc
+ PetscErrorCode       :: ierr
+
+ call PetscSectionGetFieldComponents(section,field,numComp,ierr); CHKERRQ(ierr)
+ call ISGetSize(bcPointsIS,nBcPoints,ierr); CHKERRQ(ierr)
+ if (nBcPoints > 0) call ISGetIndicesF90(bcPointsIS,bcPoints,ierr)
+ call VecGetArrayF90(localVec,localArray,ierr); CHKERRQ(ierr)
+ do point = 1, nBcPoints
+   call PetscSectionGetFieldDof(section,bcPoints(point),field,numDof,ierr)
+   CHKERRQ(ierr)
+   call PetscSectionGetFieldOffset(section,bcPoints(point),field,offset,ierr)
+   CHKERRQ(ierr)
+   do dof = offset+comp+1, offset+numDof, numComp
+     localArray(dof) = localArray(dof) + BCValue + BCDotValue*timeinc
+   enddo
+ enddo    
+ call VecRestoreArrayF90(localVec,localArray,ierr); CHKERRQ(ierr)
+ call VecAssemblyBegin(localVec, ierr); CHKERRQ(ierr)
+ call VecAssemblyEnd  (localVec, ierr); CHKERRQ(ierr)
+ if (nBcPoints > 0) call ISRestoreIndicesF90(bcPointsIS,bcPoints,ierr)
+ 
+end subroutine utilities_projectBCValues
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Create index sets of boundary dofs (in local and global numbering)
+!--------------------------------------------------------------------------------------------------
+subroutine utilities_indexActiveSet(field,section,x_local,f_local,localIS,globalIS)
+ use mesh, only: &
+   geomMesh
+
+ implicit none
+ 
+ ISLocalToGlobalMapping                   :: local2global
+ PetscSection                             :: section
+ Vec                                      :: x_local, f_local
+ PetscInt                                 :: field
+ IS                                       :: localIS, globalIS, dummyIS
+ PetscScalar, dimension(:)  ,     pointer :: x_scal, f_scal
+ PetscInt                                 :: ISSize
+ PetscInt                                 :: chart, chartStart, chartEnd, nDof, dof, offset
+ PetscInt                                 :: localSize
+ PetscInt,    dimension(:)  , allocatable :: localIndices
+ PetscInt,    dimension(:)  ,     pointer :: pBCvertex, pBCvertexlc
+ PetscErrorCode                           :: ierr
+
+ call DMGetLocalToGlobalMapping(geomMesh,local2global,ierr)
+ CHKERRQ(ierr)
+ call DMPlexGetChart(geomMesh,chartStart,chartEnd,ierr)
+ CHKERRQ(ierr)
+ call VecGetArrayF90(x_local,x_scal,ierr); CHKERRQ(ierr)
+ call VecGetArrayF90(f_local,f_scal,ierr); CHKERRQ(ierr)
+ localSize = 0
+ do chart = chartStart, chartEnd-1
+   call PetscSectionGetFieldDof(section,chart,field-1,nDof,ierr); CHKERRQ(ierr)
+   call PetscSectionGetFieldOffset(section,chart,field-1,offset,ierr); CHKERRQ(ierr)
+   do dof = offset+1, offset+nDof
+     if (((x_scal(dof) <       1.0e-8) .and. (f_scal(dof) > 0.0)) .or. &
+         ((x_scal(dof) > 1.0 - 1.0e-8) .and. (f_scal(dof) < 0.0))) localSize = localSize + 1
+   enddo
+ enddo
+ allocate(localIndices(localSize))        
+ localSize = 0
+ do chart = chartStart, chartEnd-1
+   call PetscSectionGetFieldDof(section,chart,field-1,nDof,ierr); CHKERRQ(ierr)
+   call PetscSectionGetFieldOffset(section,chart,field-1,offset,ierr); CHKERRQ(ierr)
+   do dof = offset+1, offset+nDof
+     if (((x_scal(dof) <       1.0e-8) .and. (f_scal(dof) > 0.0)) .or. &
+         ((x_scal(dof) > 1.0 - 1.0e-8) .and. (f_scal(dof) < 0.0))) then
+       localSize = localSize + 1
+       localIndices(localSize) = dof-1
+     endif  
+   enddo
+ enddo
+ call VecRestoreArrayF90(x_local,x_scal,ierr); CHKERRQ(ierr)
+ call VecRestoreArrayF90(f_local,f_scal,ierr); CHKERRQ(ierr)
+ call ISCreateGeneral(PETSC_COMM_SELF,localSize,localIndices,PETSC_COPY_VALUES,localIS,ierr)
+ CHKERRQ(ierr)
+ call ISLocalToGlobalMappingApplyIS(local2global,localIS,globalIS,ierr)
+ CHKERRQ(ierr)
+ call ISGetSize(globalIS,ISSize,ierr); CHKERRQ(ierr)
+ if (ISSize > 0) then
+   call ISGetIndicesF90(localIS,pBCvertexlc,ierr); CHKERRQ(ierr)
+   call ISGetIndicesF90(globalIS,pBCvertex,ierr); CHKERRQ(ierr)
+ endif  
+ localSize = 0
+ do chart = 1, ISSize
+   if (pBCvertex(chart) >= 0) then
+     localSize = localSize + 1
+     localIndices(localSize) = pBCvertexlc(chart)
+   endif  
+ enddo
+ if (ISSize > 0) then
+   call ISRestoreIndicesF90(localIS,pBCvertexlc,ierr); CHKERRQ(ierr)
+   call ISRestoreIndicesF90(globalIS,pBCvertex,ierr); CHKERRQ(ierr)
+ endif  
+ call ISDestroy(globalIS,ierr); CHKERRQ(ierr) 
+ call ISCreateGeneral(PETSC_COMM_SELF,localSize,localIndices,PETSC_COPY_VALUES,globalIS,ierr) 
+ CHKERRQ(ierr)
+ if (ISSize > 0) then
+   call ISDuplicate(localIS,dummyIS,ierr); CHKERRQ(ierr)
+   call ISDestroy(localIS,ierr); CHKERRQ(ierr)
+   call ISDifference(dummyIS,globalIS,localIS,ierr)
+   CHKERRQ(ierr)
+   call ISDestroy(dummyIS,ierr); CHKERRQ(ierr)
+ endif                      
+ deallocate(localIndices)
+ 
+end subroutine utilities_indexActiveSet
+
+!--------------------------------------------------------------------------------------------------
+!> @brief cleans up
+!--------------------------------------------------------------------------------------------------
+subroutine utilities_destroy()
+ !use material, only: &
+ !  homogenization_Ngrains
+
+ !implicit none
+ !PetscInt       :: homog, cryst, grain, phase 
+ !PetscErrorCode :: ierr
+
+ !call PetscViewerHDF5PopGroup(resUnit, ierr); CHKERRQ(ierr)
+ !call VecDestroy(coordinatesVec,ierr); CHKERRQ(ierr)
+ !do homog = 1, material_Nhomogenization
+ !  call VecDestroy(homogenizationResultsVec(homog),ierr);CHKERRQ(ierr)
+ !  do cryst = 1, material_Ncrystallite; do grain = 1, homogenization_Ngrains(homog)
+ !    call VecDestroy(crystalliteResultsVec(cryst,grain),ierr);CHKERRQ(ierr)
+ !  enddo; enddo
+ !  do phase = 1, material_Nphase; do grain = 1, homogenization_Ngrains(homog)
+ !    call VecDestroy(phaseResultsVec(phase,grain),ierr);CHKERRQ(ierr)
+ !  enddo; enddo  
+ !enddo      
+ !call PetscViewerDestroy(resUnit, ierr); CHKERRQ(ierr)
+
+end subroutine utilities_destroy
+
+
+end module FEM_utilities

src/FEM_zoo.f90

@@ -0,0 +1,350 @@
+!--------------------------------------------------------------------------------------------------
+!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
+!> @brief Interpolation data used by the FEM solver
+!--------------------------------------------------------------------------------------------------
+module FEM_Zoo     
+ use prec, only: pReal, pInt, group_float
+
+ implicit none
+ private
+ integer(pInt), parameter, public:: &
+   maxOrder = 5                                                                            !< current max interpolation set at cubic (intended to be arbitrary)
+ real(pReal),   dimension(2,3),       private, protected :: &
+   triangle    = reshape([-1.0_pReal, -1.0_pReal, &
+                           1.0_pReal, -1.0_pReal, &
+                          -1.0_pReal,  1.0_pReal], shape=[2,3])
+ real(pReal),   dimension(3,4),       private, protected :: &
+   tetrahedron = reshape([-1.0_pReal, -1.0_pReal, -1.0_pReal, &
+                           1.0_pReal, -1.0_pReal, -1.0_pReal, &
+                          -1.0_pReal,  1.0_pReal, -1.0_pReal, &
+                          -1.0_pReal, -1.0_pReal,  1.0_pReal], shape=[3,4])
+ integer(pInt), dimension(3,maxOrder), public, protected :: &
+   FEM_Zoo_nQuadrature                                                                       !< number of quadrature points for a given spatial dimension(1-3) and interpolation order(1-maxOrder)
+ type(group_float),   dimension(3,maxOrder), public, protected :: &
+   FEM_Zoo_QuadratureWeights, &                                                              !< quadrature weights for each quadrature rule 
+   FEM_Zoo_QuadraturePoints                                                                  !< quadrature point coordinates (in simplical system) for each quadrature rule
+
+ public :: &
+   FEM_Zoo_init
+
+contains
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief initializes FEM interpolation data
+!--------------------------------------------------------------------------------------------------
+subroutine FEM_Zoo_init
+#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
+ use, intrinsic :: iso_fortran_env, only: &
+   compiler_version, &
+   compiler_options
+#endif
+ use IO, only: &
+   IO_timeStamp
+   
+ implicit none
+  
+ write(6,'(/,a)')   ' <<<+-  FEM_Zoo init  -+>>>'
+ write(6,'(a)')     ' $Id: FEM_Zoo.f90 4354 2015-08-04 15:04:53Z MPIE\p.shanthraj $'
+ write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
+#include "compilation_info.f90"
+
+!--------------------------------------------------------------------------------------------------
+! 2D linear
+ FEM_Zoo_nQuadrature(2,1) = 1
+ allocate(FEM_Zoo_QuadratureWeights(2,1)%p(1))
+ allocate(FEM_Zoo_QuadraturePoints (2,1)%p(2))
+ FEM_Zoo_QuadratureWeights(2,1)%p(1) = 1.0_pReal
+ call FEM_Zoo_permutationStar3([1.0_pReal/3.0_pReal], &
+                              FEM_Zoo_QuadraturePoints(2,1)%p(1:2))
+ 
+!--------------------------------------------------------------------------------------------------
+! 2D quadratic
+ FEM_Zoo_nQuadrature(2,2) = 3
+ allocate(FEM_Zoo_QuadratureWeights(2,2)%p(3))
+ allocate(FEM_Zoo_QuadraturePoints (2,2)%p(6))
+ FEM_Zoo_QuadratureWeights(2,2)%p(1:3) = 1.0_pReal/3.0_pReal
+ call FEM_Zoo_permutationStar21([1.0_pReal/6.0_pReal], &
+                               FEM_Zoo_QuadraturePoints(2,2)%p(1:6))
+ 
+!--------------------------------------------------------------------------------------------------
+! 2D cubic
+ FEM_Zoo_nQuadrature(2,3) = 6
+ allocate(FEM_Zoo_QuadratureWeights(2,3)%p(6 ))
+ allocate(FEM_Zoo_QuadraturePoints (2,3)%p(12))
+ FEM_Zoo_QuadratureWeights(2,3)%p(1:3) = 0.22338158967801146570_pReal
+ call FEM_Zoo_permutationStar21([0.44594849091596488632_pReal], &
+                               FEM_Zoo_QuadraturePoints(2,3)%p(1:6))
+ FEM_Zoo_QuadratureWeights(2,3)%p(4:6) = 0.10995174365532186764_pReal
+ call FEM_Zoo_permutationStar21([0.091576213509770743460_pReal], &
+                               FEM_Zoo_QuadraturePoints(2,3)%p(7:12))
+
+!--------------------------------------------------------------------------------------------------
+! 2D quartic
+ FEM_Zoo_nQuadrature(2,4) = 12
+ allocate(FEM_Zoo_QuadratureWeights(2,4)%p(12))
+ allocate(FEM_Zoo_QuadraturePoints (2,4)%p(24))
+ FEM_Zoo_QuadratureWeights(2,4)%p(1:3) = 0.11678627572638_pReal
+ call FEM_Zoo_permutationStar21([0.24928674517091_pReal], &
+                               FEM_Zoo_QuadraturePoints(2,4)%p(1:6))
+ FEM_Zoo_QuadratureWeights(2,4)%p(4:6) = 0.05084490637021_pReal
+ call FEM_Zoo_permutationStar21([0.06308901449150_pReal], &
+                               FEM_Zoo_QuadraturePoints(2,4)%p(7:12))
+ FEM_Zoo_QuadratureWeights(2,4)%p(7:12) = 0.08285107561837_pReal
+ call FEM_Zoo_permutationStar111([0.31035245103378_pReal, 0.63650249912140_pReal], &
+                                FEM_Zoo_QuadraturePoints(2,4)%p(13:24))
+
+!--------------------------------------------------------------------------------------------------
+! 2D order 5
+ FEM_Zoo_nQuadrature(2,5) = 16
+ allocate(FEM_Zoo_QuadratureWeights(2,5)%p(16))
+ allocate(FEM_Zoo_QuadraturePoints (2,5)%p(32))
+ FEM_Zoo_QuadratureWeights(2,5)%p(1  ) = 0.14431560767779_pReal
+ call FEM_Zoo_permutationStar3([0.33333333333333_pReal], &
+                              FEM_Zoo_QuadraturePoints(2,5)%p(1:2))
+ FEM_Zoo_QuadratureWeights(2,5)%p(2:4) = 0.09509163426728_pReal
+ call FEM_Zoo_permutationStar21([0.45929258829272_pReal], &
+                               FEM_Zoo_QuadraturePoints(2,5)%p(3:8))
+ FEM_Zoo_QuadratureWeights(2,5)%p(5:7) = 0.10321737053472_pReal
+ call FEM_Zoo_permutationStar21([0.17056930775176_pReal], &
+                               FEM_Zoo_QuadraturePoints(2,5)%p(9:14))
+ FEM_Zoo_QuadratureWeights(2,5)%p(8:10) = 0.03245849762320_pReal
+ call FEM_Zoo_permutationStar21([0.05054722831703_pReal], &
+                               FEM_Zoo_QuadraturePoints(2,5)%p(15:20))
+ FEM_Zoo_QuadratureWeights(2,5)%p(11:16) = 0.02723031417443_pReal
+ call FEM_Zoo_permutationStar111([0.26311282963464_pReal, 0.72849239295540_pReal], &
+                                FEM_Zoo_QuadraturePoints(2,5)%p(21:32))
+
+!--------------------------------------------------------------------------------------------------
+! 3D linear
+ FEM_Zoo_nQuadrature(3,1) = 1
+ allocate(FEM_Zoo_QuadratureWeights(3,1)%p(1))
+ allocate(FEM_Zoo_QuadraturePoints (3,1)%p(3))
+ FEM_Zoo_QuadratureWeights(3,1)%p(1) = 1.0_pReal
+ call FEM_Zoo_permutationStar4([0.25_pReal], &
+                              FEM_Zoo_QuadraturePoints(3,1)%p(1:3))
+ 
+!--------------------------------------------------------------------------------------------------
+! 3D quadratic
+ FEM_Zoo_nQuadrature(3,2) = 4
+ allocate(FEM_Zoo_QuadratureWeights(3,2)%p(4 ))
+ allocate(FEM_Zoo_QuadraturePoints (3,2)%p(12))
+ FEM_Zoo_QuadratureWeights(3,2)%p(1:4) = 0.25_pReal
+ call FEM_Zoo_permutationStar31([0.13819660112501051518_pReal], &
+                               FEM_Zoo_QuadraturePoints(3,2)%p(1:12))
+ 
+!--------------------------------------------------------------------------------------------------
+! 3D cubic
+ FEM_Zoo_nQuadrature(3,3) = 14
+ allocate(FEM_Zoo_QuadratureWeights(3,3)%p(14))
+ allocate(FEM_Zoo_QuadraturePoints (3,3)%p(42))
+ FEM_Zoo_QuadratureWeights(3,3)%p(1:4) = 0.073493043116361949544_pReal
+ call FEM_Zoo_permutationStar31([0.092735250310891226402_pReal], &
+                               FEM_Zoo_QuadraturePoints(3,3)%p(1:12))
+ FEM_Zoo_QuadratureWeights(3,3)%p(5:8) = 0.11268792571801585080_pReal
+ call FEM_Zoo_permutationStar31([0.31088591926330060980_pReal], &
+                               FEM_Zoo_QuadraturePoints(3,3)%p(13:24))
+ FEM_Zoo_QuadratureWeights(3,3)%p(9:14) = 0.042546020777081466438_pReal
+ call FEM_Zoo_permutationStar22([0.045503704125649649492_pReal], &
+                               FEM_Zoo_QuadraturePoints(3,3)%p(25:42))
+
+!--------------------------------------------------------------------------------------------------
+! 3D quartic
+ FEM_Zoo_nQuadrature(3,4) = 35
+ allocate(FEM_Zoo_QuadratureWeights(3,4)%p(35))
+ allocate(FEM_Zoo_QuadraturePoints (3,4)%p(105))
+ FEM_Zoo_QuadratureWeights(3,4)%p(1:4) = 0.0021900463965388_pReal
+ call FEM_Zoo_permutationStar31([0.0267367755543735_pReal], &
+                               FEM_Zoo_QuadraturePoints(3,4)%p(1:12))
+ FEM_Zoo_QuadratureWeights(3,4)%p(5:16) = 0.0143395670177665_pReal
+ call FEM_Zoo_permutationStar211([0.0391022406356488_pReal, 0.7477598884818090_pReal], &
+                                FEM_Zoo_QuadraturePoints(3,4)%p(13:48))
+ FEM_Zoo_QuadratureWeights(3,4)%p(17:22) = 0.0250305395686746_pReal
+ call FEM_Zoo_permutationStar22([0.4547545999844830_pReal], &
+                               FEM_Zoo_QuadraturePoints(3,4)%p(49:66))
+ FEM_Zoo_QuadratureWeights(3,4)%p(23:34) = 0.0479839333057554_pReal
+ call FEM_Zoo_permutationStar211([0.2232010379623150_pReal, 0.0504792790607720_pReal], &
+                                FEM_Zoo_QuadraturePoints(3,4)%p(67:102))
+ FEM_Zoo_QuadratureWeights(3,4)%p(35) = 0.0931745731195340_pReal
+ call FEM_Zoo_permutationStar4([0.25_pReal], &
+                              FEM_Zoo_QuadraturePoints(3,4)%p(103:105))
+!--------------------------------------------------------------------------------------------------
+! 3D quintic
+ FEM_Zoo_nQuadrature(3,5) = 56
+ allocate(FEM_Zoo_QuadratureWeights(3,5)%p(56))
+ allocate(FEM_Zoo_QuadraturePoints (3,5)%p(168))
+ FEM_Zoo_QuadratureWeights(3,5)%p(1:4) = 0.0010373112336140_pReal
+ call FEM_Zoo_permutationStar31([0.0149520651530592_pReal], &
+                               FEM_Zoo_QuadraturePoints(3,5)%p(1:12))
+ FEM_Zoo_QuadratureWeights(3,5)%p(5:16) = 0.0096016645399480_pReal
+ call FEM_Zoo_permutationStar211([0.0340960211962615_pReal, 0.1518319491659370_pReal], &
+                                FEM_Zoo_QuadraturePoints(3,5)%p(13:48))
+ FEM_Zoo_QuadratureWeights(3,5)%p(17:28) = 0.0164493976798232_pReal
+ call FEM_Zoo_permutationStar211([0.0462051504150017_pReal, 0.3549340560639790_pReal], &
+                                FEM_Zoo_QuadraturePoints(3,5)%p(49:84))
+ FEM_Zoo_QuadratureWeights(3,5)%p(29:40) = 0.0153747766513310_pReal
+ call FEM_Zoo_permutationStar211([0.2281904610687610_pReal, 0.0055147549744775_pReal], &
+                                FEM_Zoo_QuadraturePoints(3,5)%p(85:120))
+ FEM_Zoo_QuadratureWeights(3,5)%p(41:52) = 0.0293520118375230_pReal
+ call FEM_Zoo_permutationStar211([0.3523052600879940_pReal, 0.0992057202494530_pReal], &
+                                FEM_Zoo_QuadraturePoints(3,5)%p(121:156))
+ FEM_Zoo_QuadratureWeights(3,5)%p(53:56) = 0.0366291366405108_pReal
+ call FEM_Zoo_permutationStar31([0.1344783347929940_pReal], &
+                               FEM_Zoo_QuadraturePoints(3,5)%p(157:168))
+
+end subroutine FEM_Zoo_init
+
+!--------------------------------------------------------------------------------------------------
+!> @brief star 3 permutation of input 
+!--------------------------------------------------------------------------------------------------
+subroutine FEM_Zoo_permutationStar3(point,qPt)
+
+ implicit none
+ real(pReal)                  :: point(1), qPt(2,1), temp(3,1)
+ 
+ temp(:,1) = [point(1), point(1), point(1)]
+ qPt = matmul(triangle, temp)
+  
+end subroutine FEM_Zoo_permutationStar3
+
+!--------------------------------------------------------------------------------------------------
+!> @brief star 21 permutation of input 
+!--------------------------------------------------------------------------------------------------
+subroutine FEM_Zoo_permutationStar21(point,qPt)
+
+ implicit none
+ real(pReal)                  :: point(1), qPt(2,3), temp(3,3)
+ 
+ temp(:,1) = [point(1), point(1), 1.0_pReal - 2.0_pReal*point(1)]
+ temp(:,2) = [point(1), 1.0_pReal - 2.0_pReal*point(1), point(1)]
+ temp(:,3) = [1.0_pReal - 2.0_pReal*point(1), point(1), point(1)]
+ qPt = matmul(triangle, temp)
+  
+end subroutine FEM_Zoo_permutationStar21
+
+!--------------------------------------------------------------------------------------------------
+!> @brief star 111 permutation of input 
+!--------------------------------------------------------------------------------------------------
+subroutine FEM_Zoo_permutationStar111(point,qPt)
+
+ implicit none
+ real(pReal)                  :: point(2), qPt(2,6), temp(3,6)
+ 
+ temp(:,1) = [point(1), point(2), 1.0_pReal - point(1) - point(2)]
+ temp(:,2) = [point(1), 1.0_pReal - point(1) - point(2), point(2)]
+ temp(:,4) = [point(2), 1.0_pReal - point(1) - point(2), point(1)]
+ temp(:,5) = [1.0_pReal - point(1) - point(2), point(2), point(1)]
+ temp(:,6) = [1.0_pReal - point(1) - point(2), point(1), point(2)]
+ qPt = matmul(triangle, temp)
+  
+end subroutine FEM_Zoo_permutationStar111
+
+!--------------------------------------------------------------------------------------------------
+!> @brief star 4 permutation of input 
+!--------------------------------------------------------------------------------------------------
+subroutine FEM_Zoo_permutationStar4(point,qPt)
+
+ implicit none
+ real(pReal)                  :: point(1), qPt(3,1), temp(4,1)
+ 
+ temp(:,1) = [point(1), point(1), point(1), point(1)]
+ qPt = matmul(tetrahedron, temp)
+  
+end subroutine FEM_Zoo_permutationStar4
+
+!--------------------------------------------------------------------------------------------------
+!> @brief star 31 permutation of input 
+!--------------------------------------------------------------------------------------------------
+subroutine FEM_Zoo_permutationStar31(point,qPt)
+
+ implicit none
+ real(pReal)                  :: point(1), qPt(3,4), temp(4,4)
+ 
+ temp(:,1) = [point(1), point(1), point(1), 1.0_pReal - 3.0_pReal*point(1)]
+ temp(:,2) = [point(1), point(1), 1.0_pReal - 3.0_pReal*point(1), point(1)]
+ temp(:,3) = [point(1), 1.0_pReal - 3.0_pReal*point(1), point(1), point(1)]
+ temp(:,4) = [1.0_pReal - 3.0_pReal*point(1), point(1), point(1), point(1)]
+ qPt = matmul(tetrahedron, temp)
+  
+end subroutine FEM_Zoo_permutationStar31
+
+!--------------------------------------------------------------------------------------------------
+!> @brief star 22 permutation of input 
+!--------------------------------------------------------------------------------------------------
+subroutine FEM_Zoo_permutationStar22(point,qPt)
+
+ implicit none
+ real(pReal)                  :: point(1), qPt(3,6), temp(4,6)
+ 
+ temp(:,1) = [point(1), point(1), 0.5_pReal - point(1), 0.5_pReal - point(1)]
+ temp(:,2) = [point(1), 0.5_pReal - point(1), point(1), 0.5_pReal - point(1)]
+ temp(:,3) = [0.5_pReal - point(1), point(1), point(1), 0.5_pReal - point(1)]
+ temp(:,4) = [0.5_pReal - point(1), point(1), 0.5_pReal - point(1), point(1)]
+ temp(:,5) = [0.5_pReal - point(1), 0.5_pReal - point(1), point(1), point(1)]
+ temp(:,6) = [point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), point(1)]
+ qPt = matmul(tetrahedron, temp)
+  
+end subroutine FEM_Zoo_permutationStar22
+
+!--------------------------------------------------------------------------------------------------
+!> @brief star 211 permutation of input 
+!--------------------------------------------------------------------------------------------------
+subroutine FEM_Zoo_permutationStar211(point,qPt)
+
+ implicit none
+ real(pReal)                  :: point(2), qPt(3,12), temp(4,12)
+ 
+ temp(:,1 ) = [point(1), point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2)]
+ temp(:,2 ) = [point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2)]
+ temp(:,3 ) = [point(1), point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2)]
+ temp(:,4 ) = [point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1)]
+ temp(:,5 ) = [point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2)]
+ temp(:,6 ) = [point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1)]
+ temp(:,7 ) = [point(2), point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2)]
+ temp(:,8 ) = [point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1)]
+ temp(:,9 ) = [point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1)]
+ temp(:,10) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1), point(2)]
+ temp(:,11) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2), point(1)]
+ temp(:,12) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1), point(1)]
+ qPt = matmul(tetrahedron, temp)
+  
+end subroutine FEM_Zoo_permutationStar211
+
+!--------------------------------------------------------------------------------------------------
+!> @brief star 1111 permutation of input 
+!--------------------------------------------------------------------------------------------------
+subroutine FEM_Zoo_permutationStar1111(point,qPt)
+
+ implicit none
+ real(pReal)                  :: point(3), qPt(3,24), temp(4,24)
+ 
+ temp(:,1 ) = [point(1), point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3)]
+ temp(:,2 ) = [point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3)]
+ temp(:,3 ) = [point(1), point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3)]
+ temp(:,4 ) = [point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2)]
+ temp(:,5 ) = [point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(3)]
+ temp(:,6 ) = [point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(2)]
+ temp(:,7 ) = [point(2), point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3)]
+ temp(:,8 ) = [point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3)]
+ temp(:,9 ) = [point(2), point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3)]
+ temp(:,10) = [point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1)]
+ temp(:,11) = [point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(3)]
+ temp(:,12) = [point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(1)]
+ temp(:,13) = [point(3), point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3)]
+ temp(:,14) = [point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2)]
+ temp(:,15) = [point(3), point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3)]
+ temp(:,16) = [point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1)]
+ temp(:,17) = [point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(2)]
+ temp(:,18) = [point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(1)]
+ temp(:,19) = [1.0_pReal - point(1) - point(2)- point(3), point(1), point(2), point(3)]
+ temp(:,20) = [1.0_pReal - point(1) - point(2)- point(3), point(1), point(3), point(2)]
+ temp(:,21) = [1.0_pReal - point(1) - point(2)- point(3), point(2), point(1), point(3)]
+ temp(:,22) = [1.0_pReal - point(1) - point(2)- point(3), point(2), point(3), point(1)]
+ temp(:,23) = [1.0_pReal - point(1) - point(2)- point(3), point(3), point(1), point(2)]
+ temp(:,24) = [1.0_pReal - point(1) - point(2)- point(3), point(3), point(2), point(1)]
+ qPt = matmul(tetrahedron, temp)
+  
+end subroutine FEM_Zoo_permutationStar1111
+
+
+end module FEM_Zoo

src/FEsolving.f90

@@ -81,20 +81,13 @@ subroutine FE_init
  modelName = getSolverJobName()
 
 #if defined(Spectral) || defined(FEM)
-
+ restartInc = interface_RestartInc
-#ifdef Spectral
- restartInc = spectralRestartInc
-#endif
-#ifdef FEM
- restartInc = FEMRestartInc
-#endif
 
  if(restartInc < 0_pInt) then
    call IO_warning(warning_ID=34_pInt)
    restartInc = 0_pInt
  endif
  restartRead = restartInc > 0_pInt                                                                  ! only read in if "true" restart requested
-
 #else
  call IO_open_inputFile(FILEUNIT,modelName)
  rewind(FILEUNIT)

src/IO.f90

@@ -22,6 +22,7 @@ module IO
  public :: &
    IO_init, &
    IO_read, &
+   IO_recursiveRead, &
    IO_checkAndRewind, &
    IO_open_file_stat, &
    IO_open_jobFile_stat, &
@@ -35,10 +36,6 @@ module IO
    IO_hybridIA, &
    IO_isBlank, &
    IO_getTag, &
-   IO_countSections, &
-   IO_countTagInPart, &
-   IO_spotTagInPart, &
-   IO_globalTagInPart, &
    IO_stringPos, &
    IO_stringValue, &
    IO_fixedStringValue ,&
@@ -100,6 +97,7 @@ end subroutine IO_init
 !--------------------------------------------------------------------------------------------------
 !> @brief recursively reads a line from a text file.
 !!        Recursion is triggered by "{path/to/inputfile}" in a line
+!> @details unstable and buggy
 !--------------------------------------------------------------------------------------------------
 recursive function IO_read(fileUnit,reset) result(line)
 
@@ -151,7 +149,7 @@ recursive function IO_read(fileUnit,reset) result(line)
    pathOn(stack) = path(1:scan(path,SEP,.true.))//input                                             ! glue include to current file's dir
  endif
 
- open(newunit=unitOn(stack),iostat=myStat,file=pathOn(stack),action='read')                         ! open included file
+ open(newunit=unitOn(stack),iostat=myStat,file=pathOn(stack),action='read',status='old',position='rewind')                         ! open included file
  if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=pathOn(stack))
 
  line = IO_read(fileUnit)
@@ -170,6 +168,80 @@ recursive function IO_read(fileUnit,reset) result(line)
 
 end function IO_read
 
+!--------------------------------------------------------------------------------------------------
+!> @brief recursively reads a text file.
+!!        Recursion is triggered by "{path/to/inputfile}" in a line
+!--------------------------------------------------------------------------------------------------
+recursive function IO_recursiveRead(fileName,cnt) result(fileContent)
+
+  implicit none
+  character(len=*),   intent(in)                :: fileName
+  integer(pInt),      intent(in), optional      :: cnt                                              !< recursion counter
+  character(len=256), dimension(:), allocatable :: fileContent                                      !< file content, separated per lines
+  character(len=256), dimension(:), allocatable :: includedContent
+  character(len=256)                            :: line
+  character(len=256), parameter                 :: dummy = 'https://damask.mpie.de'                 !< to fill up remaining array
+  character(len=:),                 allocatable :: rawData
+  integer(pInt) ::  &
+    fileLength, &
+    fileUnit, &
+    startPos, endPos, &
+    myTotalLines, &                                                                                 !< # lines read from file without include statements
+    includedLines, &                                                                                !< # lines included from other file(s)
+    missingLines, &                                                                                 !< # lines missing from current file
+    l,i, &
+    myStat
+
+  if (merge(cnt,0_pInt,present(cnt))>10_pInt) call IO_error(106_pInt,ext_msg=trim(fileName))
+
+!--------------------------------------------------------------------------------------------------
+! read data as stream
+  inquire(file = fileName, size=fileLength)
+  open(newunit=fileUnit, file=fileName, access='stream',&
+       status='old', position='rewind', action='read',iostat=myStat)
+  if(myStat /= 0_pInt) call IO_error(100_pInt,ext_msg=trim(fileName))
+  allocate(character(len=fileLength)::rawData)
+  read(fileUnit) rawData
+  close(fileUnit)
+
+!--------------------------------------------------------------------------------------------------
+! count lines to allocate string array
+  myTotalLines = 0_pInt
+  do l=1_pInt, len(rawData)
+    if (rawData(l:l) == new_line('')) myTotalLines = myTotalLines+1
+  enddo
+  allocate(fileContent(myTotalLines))
+
+!--------------------------------------------------------------------------------------------------
+! split raw data at end of line and handle includes
+  startPos = 1_pInt
+  endPos = 0_pInt
+
+  includedLines=0_pInt
+  l=0_pInt
+  do while (startPos <= len(rawData))
+    l = l + 1_pInt
+    endPos = endPos + scan(rawData(startPos:),new_line(''))
+    if(endPos - startPos >256) call IO_error(107_pInt,ext_msg=trim(fileName))
+    line = rawData(startPos:endPos-1_pInt)
+    startPos = endPos + 1_pInt
+
+    recursion: if(scan(trim(line),'{') < scan(trim(line),'}')) then
+      myTotalLines    = myTotalLines - 1_pInt
+      includedContent = IO_recursiveRead(trim(line(scan(line,'{')+1_pInt:scan(line,'}')-1_pInt)), &
+                        merge(cnt,1_pInt,present(cnt)))                                             ! to track recursion depth
+      includedLines   = includedLines + size(includedContent)
+      missingLines    = myTotalLines + includedLines - size(fileContent(1:l-1)) -size(includedContent)
+      fileContent     = [ fileContent(1:l-1_pInt), includedContent, [(dummy,i=1,missingLines)] ]    ! add content and grow array
+      l = l - 1_pInt + size(includedContent)
+    else recursion
+      fileContent(l) = line
+    endif recursion
+
+  enddo
+
+end function IO_recursiveRead
+
 
 !--------------------------------------------------------------------------------------------------
 !> @brief checks if unit is opened for reading, if true rewinds. Otherwise stops with
@@ -203,7 +275,7 @@ subroutine IO_open_file(fileUnit,path)
 
  integer(pInt)                  :: myStat
 
- open(fileUnit,status='old',iostat=myStat,file=path)
+ open(fileUnit,status='old',iostat=myStat,file=path,action='read',position='rewind')
  if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=path)
 
 end subroutine IO_open_file
@@ -222,7 +294,8 @@ logical function IO_open_file_stat(fileUnit,path)
 
  integer(pInt)                  :: myStat
 
- open(fileUnit,status='old',iostat=myStat,file=path)
+ open(fileUnit,status='old',iostat=myStat,file=path,action='read',position='rewind')
+ if (myStat /= 0_pInt) close(fileUnit)
  IO_open_file_stat = (myStat == 0_pInt)
 
 end function IO_open_file_stat
@@ -246,7 +319,7 @@ subroutine IO_open_jobFile(fileUnit,ext)
  character(len=1024)            :: path
 
  path = trim(getSolverJobName())//'.'//ext
- open(fileUnit,status='old',iostat=myStat,file=path)
+ open(fileUnit,status='old',iostat=myStat,file=path,action='read',position='rewind')
  if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=path)
 
 end subroutine IO_open_jobFile
@@ -270,7 +343,8 @@ logical function IO_open_jobFile_stat(fileUnit,ext)
  character(len=1024)            :: path
 
  path = trim(getSolverJobName())//'.'//ext
- open(fileUnit,status='old',iostat=myStat,file=path)
+ open(fileUnit,status='old',iostat=myStat,file=path,action='read',position='rewind')
+ if (myStat /= 0_pInt) close(fileUnit)
  IO_open_jobFile_stat = (myStat == 0_pInt)
 
 end function IO_open_JobFile_stat
@@ -296,11 +370,11 @@ subroutine IO_open_inputFile(fileUnit,modelName)
 
  fileType = 1_pInt                                                                                  ! assume .pes
  path = trim(modelName)//inputFileExtension(fileType)                                               ! attempt .pes, if it exists: it should be used
- open(fileUnit+1,status='old',iostat=myStat,file=path)
+ open(fileUnit+1,status='old',iostat=myStat,file=path,action='read',position='rewind')
  if(myStat /= 0_pInt) then                                                                          ! if .pes does not work / exist; use conventional extension, i.e.".inp"
     fileType = 2_pInt
     path = trim(modelName)//inputFileExtension(fileType)
-    open(fileUnit+1,status='old',iostat=myStat,file=path)
+    open(fileUnit+1,status='old',iostat=myStat,file=path,action='read',position='rewind')
  endif
  if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=path)
 
@@ -335,7 +409,7 @@ subroutine IO_open_logFile(fileUnit)
  character(len=1024)            :: path
 
  path = trim(getSolverJobName())//LogFileExtension
- open(fileUnit,status='old',iostat=myStat,file=path)
+ open(fileUnit,status='old',iostat=myStat,file=path,action='read',position='rewind')
  if (myStat /= 0_pInt) call IO_error(100_pInt,el=myStat,ext_msg=path)
 
 end subroutine IO_open_logFile
@@ -755,16 +829,22 @@ pure function IO_getTag(string,openChar,closeChar)
  character(len=*), intent(in)  :: string                                                            !< string to check for tag
  character(len=len_trim(string)) :: IO_getTag
 
- character(len=*), intent(in)  :: openChar, &                                                       !< indicates beginning of tag
+ character, intent(in)  :: openChar, &                                                              !< indicates beginning of tag
                            closeChar                                                                !< indicates end of tag
 
  character(len=*), parameter   :: SEP=achar(32)//achar(9)//achar(10)//achar(13)                     ! whitespaces
-
  integer :: left,right                                                                              ! no pInt
 
  IO_getTag = ''
+
+
+ if (openChar /= closeChar) then
    left  = scan(string,openChar)
    right = scan(string,closeChar)
+ else
+   left  = scan(string,openChar)
+   right = left + merge(scan(string(left+1:),openChar),0_pInt,len(string) > left)
+ endif
 
  if (left == verify(string,SEP) .and. right > left) &                                               ! openChar is first and closeChar occurs
    IO_getTag = string(left+1:right-1)
@@ -772,173 +852,6 @@ pure function IO_getTag(string,openChar,closeChar)
 end function IO_getTag
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief count number of [sections] in <part> for given file handle
-!--------------------------------------------------------------------------------------------------
-integer(pInt) function IO_countSections(fileUnit,part)
-
- implicit none
- integer(pInt),      intent(in) :: fileUnit                                                         !< file handle
- character(len=*),   intent(in) :: part                                                             !< part name in which sections are counted
-
- character(len=65536)           :: line
-
- line = ''
- IO_countSections = 0_pInt
- rewind(fileUnit)
-
- do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= part)                       ! search for part
-   line = IO_read(fileUnit)
- enddo
-
- do while (trim(line) /= IO_EOF)
-   line = IO_read(fileUnit)
-   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
-   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
-     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
-     exit
-   endif
-   if (IO_getTag(line,'[',']') /= '') &                                                             ! found [section] identifier
-     IO_countSections = IO_countSections + 1_pInt
- enddo
-
-end function IO_countSections
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief returns array of tag counts within <part> for at most N [sections]
-!--------------------------------------------------------------------------------------------------
-function IO_countTagInPart(fileUnit,part,tag,Nsections)
-
- implicit none
- integer(pInt),   intent(in)                :: Nsections                                            !< maximum number of sections in which tag is searched for
- integer(pInt),   dimension(Nsections)      :: IO_countTagInPart
- integer(pInt),   intent(in)                :: fileUnit                                             !< file handle
- character(len=*),intent(in)                :: part, &                                              !< part in which tag is searched for
-                                               tag                                                  !< tag to search for
-
-
- integer(pInt),   dimension(Nsections)      :: counter
- integer(pInt), allocatable, dimension(:)   :: chunkPos
- integer(pInt)                              :: section
- character(len=65536)                       :: line
-
- line = ''
- counter = 0_pInt
- section = 0_pInt
-
- rewind(fileUnit)
- do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= part)                       ! search for part
-   line = IO_read(fileUnit)
- enddo
-
- do while (trim(line) /= IO_EOF)
-   line = IO_read(fileUnit)
-   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
-   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
-     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
-     exit
-   endif
-   if (IO_getTag(line,'[',']') /= '') section = section + 1_pInt                                    ! found [section] identifier
-   if (section > 0) then
-     chunkPos = IO_stringPos(line)
-     if (tag == trim(IO_lc(IO_stringValue(line,chunkPos,1_pInt)))) &                                ! match
-       counter(section) = counter(section) + 1_pInt
-   endif
- enddo
-
- IO_countTagInPart = counter
-
-end function IO_countTagInPart
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief returns array of tag presence within <part> for at most N [sections]
-!--------------------------------------------------------------------------------------------------
-function IO_spotTagInPart(fileUnit,part,tag,Nsections)
-
- implicit none
- integer(pInt),   intent(in)                :: Nsections                                            !< maximum number of sections in which tag is searched for
- logical,         dimension(Nsections)      :: IO_spotTagInPart
- integer(pInt),   intent(in)                :: fileUnit                                             !< file handle
- character(len=*),intent(in)                :: part, &                                              !< part in which tag is searched for
-                                               tag                                                  !< tag to search for
-
-
- integer(pInt), allocatable, dimension(:) :: chunkPos
- integer(pInt)                            :: section
- character(len=65536)                     :: line
-
- IO_spotTagInPart = .false.                                                                         ! assume to nowhere spot tag
- section = 0_pInt
- line = ''
-
- rewind(fileUnit)
- do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= part)                       ! search for part
-   line = IO_read(fileUnit)
- enddo
-
- do while (trim(line) /= IO_EOF)
-   line = IO_read(fileUnit)
-   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
-   foundNextPart: if (IO_getTag(line,'<','>') /= '') then
-     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
-     exit
-   endif foundNextPart
-   if (IO_getTag(line,'[',']') /= '') section = section + 1_pInt                                    ! found [section] identifier
-   if (section > 0_pInt) then
-     chunkPos = IO_stringPos(line)
-     if (tag == trim(IO_lc(IO_stringValue(line,chunkPos,1_pInt)))) &                                ! match
-       IO_spotTagInPart(section) = .true.
-   endif
- enddo
-
- end function IO_spotTagInPart
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief return logical whether tag is present within <part> before any [sections]
-!--------------------------------------------------------------------------------------------------
-logical function IO_globalTagInPart(fileUnit,part,tag)
-
- implicit none
- integer(pInt),   intent(in)                :: fileUnit                                             !< file handle
- character(len=*),intent(in)                :: part, &                                              !< part in which tag is searched for
-                                               tag                                                  !< tag to search for
-
- integer(pInt), allocatable, dimension(:) :: chunkPos
- character(len=65536)                     :: line
-
- IO_globalTagInPart = .false.                                                                       ! assume to nowhere spot tag
- line =''
-
- rewind(fileUnit)
- do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= part)                       ! search for part
-   line = IO_read(fileUnit)
- enddo
-
- do while (trim(line) /= IO_EOF)
-   line = IO_read(fileUnit)
-   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
-   foundNextPart: if (IO_getTag(line,'<','>') /= '') then
-     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
-     exit
-   endif foundNextPart
-   foundFirstSection: if (IO_getTag(line,'[',']') /= '') then
-     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
-     exit
-   endif foundFirstSection
-   chunkPos = IO_stringPos(line)
-   match: if (tag == trim(IO_lc(IO_stringValue(line,chunkPos,1_pInt)))) then
-     IO_globalTagInPart = .true.
-     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
-     exit
-   endif match
- enddo
-
-end function IO_globalTagInPart
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief locates all space-separated chunks in given string and returns array containing number
 !! them and the left/right position to be used by IO_xxxVal
@@ -1513,6 +1426,8 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
    msg = 'unknown output:'
  case (106_pInt)
    msg = 'working directory does not exist:'
+ case (107_pInt)
+   msg = 'line length exceeds limit of 256'
 
 !--------------------------------------------------------------------------------------------------
 ! lattice error messages

src/config.f90

@@ -20,12 +20,17 @@ module config
  type, public :: tPartitionedStringList
    type(tPartitionedString)               :: string
    type(tPartitionedStringList),  pointer :: next => null()
-
    contains
      procedure :: add            => add
      procedure :: show           => show
      procedure :: free           => free
 
+! currently, a finalize is needed for all shapes of tPartitionedStringList.
+! with Fortran 2015, we can define one recursive elemental function
+! https://software.intel.com/en-us/forums/intel-visual-fortran-compiler-for-windows/topic/543326
+     final     :: finalize, &
+                  finalizeArray 
+
      procedure :: keyExists      => keyExists
      procedure :: countKeys      => countKeys
 
@@ -37,11 +42,10 @@ module config
      procedure :: getInts        => getInts
      procedure :: getStrings     => getStrings
 
+
  end type tPartitionedStringList
 
- type(tPartitionedStringList), public :: emptyList
+ type(tPartitionedStringList), public, protected, allocatable, dimension(:) :: &
-
- type(tPartitionedStringList), public, protected, allocatable, dimension(:) :: &                    ! QUESTION: rename to config_XXX?
    config_phase, &
    config_microstructure, &
    config_homogenization, &
@@ -76,7 +80,6 @@ module config
    MATERIAL_configFile         = 'material.config', &                                               !< generic name for material configuration file
    MATERIAL_localFileExt       = 'materialConfig'                                                   !< extension of solver job name depending material configuration file
 
-
  public :: &
    config_init, &
    config_deallocate
@@ -92,12 +95,14 @@ subroutine config_init()
    compiler_version, &
    compiler_options
 #endif
+ use prec, only: &
+   pStringLen
+ use DAMASK_interface, only: &
+   getSolverJobName
  use IO, only: &
    IO_error, &
-   IO_open_file, &
-   IO_read, &
    IO_lc, &
-   IO_open_jobFile_stat, &
+   IO_recursiveRead, &
    IO_getTag, &
    IO_timeStamp, &
    IO_EOF
@@ -107,12 +112,13 @@ subroutine config_init()
    debug_levelBasic
 
  implicit none
- integer(pInt), parameter :: FILEUNIT = 200_pInt
+ integer(pInt) :: myDebug,i
- integer(pInt)            :: myDebug
 
- character(len=65536) :: &                                                                          
+ character(len=pStringLen) :: &
   line, &
   part
+ character(len=pStringLen), dimension(:), allocatable :: fileContent
+ logical :: fileExists
 
  write(6,'(/,a)') ' <<<+-  config init  -+>>>'
  write(6,'(a15,a)')   ' Current time: ',IO_timeStamp()
@@ -120,39 +126,40 @@ subroutine config_init()
 
  myDebug = debug_level(debug_material)
 
- if (.not. IO_open_jobFile_stat(FILEUNIT,material_localFileExt)) &                                  ! no local material configuration present...
+ inquire(file=trim(getSolverJobName())//'.'//material_localFileExt,exist=fileExists)
-   call IO_open_file(FILEUNIT,material_configFile)                                                  ! ...open material.config file
+ if(fileExists) then
+   fileContent = IO_recursiveRead(trim(getSolverJobName())//'.'//material_localFileExt)
+ else
+   inquire(file='material.config',exist=fileExists)
+   if(.not. fileExists) call IO_error(100_pInt,ext_msg='material.config')
+   fileContent = IO_recursiveRead('material.config')
+ endif
 
- rewind(fileUnit)
+ do i = 1_pInt, size(fileContent)
- line        = ''                                                                                   ! to have it initialized
+   line = trim(fileContent(i))
- do while (trim(line) /= IO_EOF)
    part = IO_lc(IO_getTag(line,'<','>'))
-
    select case (trim(part))
     
      case (trim(material_partPhase))
-       call parseFile(line,phase_name,config_phase,FILEUNIT)
+       call parseFile(line,phase_name,config_phase,fileContent(i+1:))
        if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Phase parsed'; flush(6)
     
      case (trim(material_partMicrostructure))
-       call parseFile(line,microstructure_name,config_microstructure,FILEUNIT)
+       call parseFile(line,microstructure_name,config_microstructure,fileContent(i+1:))
        if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Microstructure parsed'; flush(6)
     
      case (trim(material_partCrystallite))
-       call parseFile(line,crystallite_name,config_crystallite,FILEUNIT)
+       call parseFile(line,crystallite_name,config_crystallite,fileContent(i+1:))
        if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Crystallite parsed'; flush(6)
     
      case (trim(material_partHomogenization))
-       call parseFile(line,homogenization_name,config_homogenization,FILEUNIT)
+       call parseFile(line,homogenization_name,config_homogenization,fileContent(i+1:))
        if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Homogenization parsed'; flush(6)
     
      case (trim(material_partTexture))
-       call parseFile(line,texture_name,config_texture,FILEUNIT)
+       call parseFile(line,texture_name,config_texture,fileContent(i+1:))
        if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Texture parsed'; flush(6)
 
-     case default
-       line = IO_read(fileUnit)
-
    end select
 
  enddo
@@ -173,107 +180,81 @@ end subroutine config_init
 !--------------------------------------------------------------------------------------------------
 !> @brief parses the material.config file
 !--------------------------------------------------------------------------------------------------
-subroutine parseFile(line,&
+subroutine parseFile(line,sectionNames,part,&
-                     sectionNames,part,fileUnit)
+                     fileContent)
+ use prec, only: &
+   pStringLen
  use IO, only: &
-   IO_read, &
    IO_error, &
-   IO_lc, &
+   IO_getTag
-   IO_getTag, &
-   IO_isBlank, &
-   IO_stringValue, &
-   IO_stringPos, &
-   IO_EOF
 
  implicit none
- integer(pInt),    intent(in) :: fileUnit
+ character(len=pStringLen),                               intent(out) :: line
- character(len=*),  dimension(:), allocatable, intent(inout)  :: sectionNames
+ character(len=64),            allocatable, dimension(:), intent(out) :: sectionNames
- type(tPartitionedStringList), allocatable, dimension(:), intent(inout) :: part
+ type(tPartitionedStringList), allocatable, dimension(:), intent(out) :: part
- character(len=65536),intent(out) :: line
+ character(len=pStringLen),                 dimension(:), intent(in)  :: fileContent
 
- integer(pInt), allocatable, dimension(:) :: chunkPos
+ integer(pInt),                allocatable, dimension(:)              :: partPosition               ! position of [] tags + last line in section
- integer(pInt)        :: s
+ integer(pInt)        :: i, j
- character(len=65536) :: devNull
- character(len=64)    :: tag
  logical              :: echo
 
  echo = .false. 
- allocate(part(0))
+ allocate(partPosition(0))
  
- s = 0_pInt
+ do i = 1_pInt, size(fileContent)
- do while (trim(line) /= IO_EOF)                                                                    ! read through sections of material part
+   line = trim(fileContent(i))
-   line = IO_read(fileUnit)
+   if (IO_getTag(line,'<','>') /= '') exit
-   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
-   foundNextPart: if (IO_getTag(line,'<','>') /= '') then
-     devNull = IO_read(fileUnit, .true.)                                                            ! reset IO_read to close any recursively included files
-     exit
-   endif foundNextPart
    nextSection: if (IO_getTag(line,'[',']') /= '') then
-     s = s + 1_pInt
+     partPosition = [partPosition, i]
-     part = [part, emptyList]
-     tag = IO_getTag(line,'[',']')
-     GfortranBug86033: if (.not. allocated(sectionNames)) then
-       allocate(sectionNames(1),source=tag)
-     else GfortranBug86033
-       sectionNames  = [sectionNames,tag]
-     endif GfortranBug86033
      cycle
    endif nextSection
-   chunkPos = IO_stringPos(line)
+   if (size(partPosition) < 1_pInt) &
-   tag = IO_lc(IO_stringValue(trim(line),chunkPos,1_pInt))                                          ! extract key
+     echo = (trim(IO_getTag(line,'/','/')) == 'echo') .or. echo
-   inSection: if (s > 0_pInt) then
-     call part(s)%add(IO_lc(trim(line)))
-   else inSection
-     echo = (trim(tag) == '/echo/')
-   endif inSection
  enddo
 
+ allocate(sectionNames(size(partPosition)))
+ allocate(part(size(partPosition)))
+
+ partPosition = [partPosition, i]                                                                   ! needed when actually storing content
+
+ do i = 1_pInt, size(partPosition) -1_pInt
+   sectionNames(i) = trim(adjustl(fileContent(partPosition(i))))
+   do j = partPosition(i) + 1_pInt,  partPosition(i+1) -1_pInt
+     call part(i)%add(trim(adjustl(fileContent(j))))
+   enddo
    if (echo) then
-   do s = 1, size(sectionNames)
+     write(6,*) 'section',i, '"'//trim(sectionNames(i))//'"'
-     call part(s)%show()
+     call part(i)%show()
-   end do
+   endif
- end if 
+ enddo
 
 end subroutine parseFile
 
+!--------------------------------------------------------------------------------------------------
+!> @brief deallocates the linked lists that store the content of the configuration files
+!--------------------------------------------------------------------------------------------------
 subroutine config_deallocate(what)
  use IO, only: &
    IO_error
 
  implicit none
  character(len=*), intent(in) :: what
- integer(pInt) :: i
 
- select case(what)
+ select case(trim(what))
 
    case('material.config/phase')
-     do i=1, size(config_phase)
-       call config_phase(i)%free
-     enddo
      deallocate(config_phase)
 
    case('material.config/microstructure')
-     do i=1, size(config_microstructure)
-       call config_microstructure(i)%free
-     enddo
      deallocate(config_microstructure)
 
    case('material.config/crystallite')
-     do i=1, size(config_crystallite)
-       call config_crystallite(i)%free
-     enddo
      deallocate(config_crystallite)
 
    case('material.config/homogenization')
-     do i=1, size(config_homogenization)
-       call config_homogenization(i)%free
-     enddo
      deallocate(config_homogenization)
 
    case('material.config/texture')
-     do i=1, size(config_texture)
-       call config_texture(i)%free
-     enddo
      deallocate(config_texture)
 
    case default
@@ -284,11 +265,17 @@ subroutine config_deallocate(what)
 end subroutine config_deallocate
 
 
+!##################################################################################################
+! The folowing functions are part of the tPartitionedStringList object
+!##################################################################################################
+
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief add element
 !> @details Adds a string together with the start/end position of chunks in this string. The new 
 !! element is added at the end of the list. Empty strings are not added. All strings are converted
-!! to lower case
+!! to lower case. The data is not stored in the new element but in the current.
 !--------------------------------------------------------------------------------------------------
 subroutine add(this,string)
   use IO, only: &
@@ -299,19 +286,18 @@ subroutine add(this,string)
   implicit none
   class(tPartitionedStringList),  target, intent(in) :: this
   character(len=*),                       intent(in) :: string
-  type(tPartitionedStringList),   pointer            :: new, item
+  type(tPartitionedStringList),   pointer            :: new, temp
 
   if (IO_isBlank(string)) return
 
   allocate(new)
-  new%string%val = IO_lc       (trim(string))
+  temp => this
-  new%string%pos = IO_stringPos(trim(string))
+  do while (associated(temp%next))
-
+    temp => temp%next
-  item => this
-  do while (associated(item%next))
-    item => item%next
   enddo
-  item%next => new
+  temp%string%val = IO_lc       (trim(string))
+  temp%string%pos = IO_stringPos(trim(string))
+  temp%next => new
 
 end subroutine add
 
@@ -323,12 +309,12 @@ end subroutine add
 subroutine show(this)
 
  implicit none
- class(tPartitionedStringList) :: this
+ class(tPartitionedStringList), target, intent(in) :: this
  type(tPartitionedStringList),  pointer            :: item
 
- item => this%next
+ item => this
- do while (associated(item))
+ do while (associated(item%next))
-   write(6,'(a)') trim(item%string%val)
+   write(6,'(a)') ' '//trim(item%string%val)
    item => item%next
  end do
 
@@ -336,28 +322,55 @@ end subroutine show
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief cleans entire list
+!> @brief empties list and frees associated memory
-!> @details list head is remains alive
+!> @details explicit interface to reset list. Triggers final statement (and following chain reaction)
 !--------------------------------------------------------------------------------------------------
 subroutine free(this)
 
   implicit none
-  class(tPartitionedStringList),  target, intent(in) :: this
+  class(tPartitionedStringList),  intent(inout) :: this
-  type(tPartitionedStringList),   pointer            :: new, item
 
-  if (.not. associated(this%next)) return
+  if(associated(this%next)) deallocate(this%next)
-
-  item => this%next
-  do while (associated(item%next))
-    new => item
-    deallocate(item)
-    item => new%next
-  enddo
-  deallocate(item)
 
 end subroutine free
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief empties list and frees associated memory
+!> @details called when variable goes out of scope. Triggers chain reaction for list
+!--------------------------------------------------------------------------------------------------
+recursive subroutine finalize(this)
+
+  implicit none
+  type(tPartitionedStringList),  intent(inout) :: this
+
+  if(associated(this%next)) deallocate(this%next)
+
+end subroutine finalize
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief cleans entire array of linke lists
+!> @details called when variable goes out of scope.
+!--------------------------------------------------------------------------------------------------
+subroutine finalizeArray(this)
+
+  implicit none
+  integer :: i
+  type(tPartitionedStringList),  intent(inout), dimension(:) :: this
+  type(tPartitionedStringList),  pointer :: temp ! bug in Gfortran?
+
+  do i=1, size(this)
+    !if (associated(this(i)%next)) then
+      temp => this(i)%next
+      !deallocate(this(i)) !internal compiler error: in gfc_build_final_call, at fortran/trans.c:975
+      deallocate(temp)
+    !endif
+  enddo
+
+end subroutine finalizeArray
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief reports wether a given key (string value at first position) exists in the list
 !--------------------------------------------------------------------------------------------------
@@ -366,14 +379,14 @@ logical function keyExists(this,key)
    IO_stringValue
 
  implicit none
- class(tPartitionedStringList), intent(in) :: this
+ class(tPartitionedStringList), target, intent(in) :: this
  character(len=*),                      intent(in) :: key
  type(tPartitionedStringList),  pointer            :: item
 
  keyExists = .false.
 
- item => this%next
+ item => this
- do while (associated(item) .and. .not. keyExists)
+ do while (associated(item%next) .and. .not. keyExists)
    keyExists = trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)
    item => item%next
  end do
@@ -391,14 +404,14 @@ integer(pInt) function countKeys(this,key)
 
  implicit none
 
- class(tPartitionedStringList), intent(in) :: this
+ class(tPartitionedStringList), target, intent(in) :: this
  character(len=*),                      intent(in) :: key
  type(tPartitionedStringList),  pointer            :: item
 
  countKeys = 0_pInt
 
- item => this%next
+ item => this
- do while (associated(item))
+ do while (associated(item%next))
    if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) &
      countKeys = countKeys + 1_pInt
    item => item%next
@@ -419,7 +432,7 @@ real(pReal) function getFloat(this,key,defaultVal)
    IO_FloatValue
 
  implicit none
- class(tPartitionedStringList), intent(in)           :: this
+ class(tPartitionedStringList), target, intent(in)           :: this
  character(len=*),                      intent(in)           :: key
  real(pReal),                           intent(in), optional :: defaultVal
  type(tPartitionedStringList), pointer                       :: item
@@ -428,8 +441,8 @@ real(pReal) function getFloat(this,key,defaultVal)
  found = present(defaultVal)
  if (found) getFloat = defaultVal
  
- item => this%next
+ item => this
- do while (associated(item))
+ do while (associated(item%next))
    if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
      found = .true.
      if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key)
@@ -455,7 +468,7 @@ integer(pInt) function getInt(this,key,defaultVal)
    IO_IntValue
 
  implicit none
- class(tPartitionedStringList), intent(in)           :: this
+ class(tPartitionedStringList), target, intent(in)           :: this
  character(len=*),                      intent(in)           :: key
  integer(pInt),                         intent(in), optional :: defaultVal
  type(tPartitionedStringList), pointer                       :: item
@@ -464,8 +477,8 @@ integer(pInt) function getInt(this,key,defaultVal)
  found = present(defaultVal)
  if (found) getInt = defaultVal
  
- item => this%next
+ item => this
- do while (associated(item))
+ do while (associated(item%next))
    if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
      found = .true.
      if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key)
@@ -491,7 +504,7 @@ character(len=65536) function getString(this,key,defaultVal,raw)
    IO_stringValue
 
  implicit none
- class(tPartitionedStringList), intent(in)           :: this
+ class(tPartitionedStringList), target, intent(in)           :: this
  character(len=*),                      intent(in)           :: key
  character(len=65536),                  intent(in), optional :: defaultVal
  logical,                               intent(in), optional :: raw
@@ -506,8 +519,8 @@ character(len=65536) function getString(this,key,defaultVal,raw)
    if (len_trim(getString) /= len_trim(defaultVal)) call IO_error(0_pInt,ext_msg='getString')
  endif
 
- item => this%next
+ item => this
- do while (associated(item))
+ do while (associated(item%next))
    if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
      found = .true.
      if (item%string%pos(1) < 2_pInt) call IO_error(143_pInt,ext_msg=key)
@@ -539,7 +552,7 @@ function getFloats(this,key,defaultVal,requiredShape)
 
  implicit none
  real(pReal),     dimension(:), allocatable          :: getFloats
- class(tPartitionedStringList), intent(in)           :: this
+ class(tPartitionedStringList), target, intent(in)   :: this
  character(len=*),              intent(in)           :: key
  real(pReal),   dimension(:),   intent(in), optional :: defaultVal
  integer(pInt), dimension(:),   intent(in), optional :: requiredShape
@@ -553,8 +566,8 @@ function getFloats(this,key,defaultVal,requiredShape)
 
  allocate(getFloats(0))
 
- item => this%next
+ item => this
- do while (associated(item))
+ do while (associated(item%next))
    if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
      found = .true.
      if (.not. cumulative) getFloats = [real(pReal)::]
@@ -586,7 +599,7 @@ function getInts(this,key,defaultVal,requiredShape)
 
  implicit none
  integer(pInt), dimension(:), allocatable            :: getInts
- class(tPartitionedStringList), intent(in)           :: this
+ class(tPartitionedStringList), target, intent(in)   :: this
  character(len=*),              intent(in)           :: key
  integer(pInt), dimension(:),   intent(in), optional :: defaultVal, &
                                                         requiredShape
@@ -600,8 +613,8 @@ function getInts(this,key,defaultVal,requiredShape)
 
  allocate(getInts(0))
 
- item => this%next
+ item => this
- do while (associated(item))
+ do while (associated(item%next))
    if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
      found = .true.
      if (.not. cumulative) getInts = [integer(pInt)::]
@@ -633,7 +646,7 @@ function getStrings(this,key,defaultVal,requiredShape,raw)
 
  implicit none
  character(len=65536),dimension(:), allocatable           :: getStrings
- class(tPartitionedStringList),      intent(in)           :: this
+ class(tPartitionedStringList), target, intent(in)        :: this
  character(len=*),                   intent(in)           :: key
  character(len=65536),dimension(:),  intent(in), optional :: defaultVal
  integer(pInt),       dimension(:),  intent(in), optional :: requiredShape
@@ -649,8 +662,8 @@ function getStrings(this,key,defaultVal,requiredShape,raw)
  whole = merge(raw,.false.,present(raw))
  found = .false.
 
- item => this%next
+ item => this
- do while (associated(item))
+ do while (associated(item%next))
    if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) then
      found = .true.
      if (allocated(getStrings) .and. .not. cumulative) deallocate(getStrings)

src/constitutive.f90

@@ -58,14 +58,15 @@ subroutine constitutive_init()
    IO_write_jobIntFile, &
    IO_timeStamp
  use config, only: &
-   config_deallocate
+   config_phase
  use mesh, only: &
    FE_geomtype
  use config, only: &
    material_Nphase, &
    material_localFileExt, &
    phase_name, &
-   material_configFile
+   material_configFile, &
+   config_deallocate
  use material, only: &
    material_phase, &
    phase_plasticity, &
@@ -138,7 +139,7 @@ subroutine constitutive_init()
  use kinematics_hydrogen_strain
 
  implicit none
- integer(pInt), parameter :: FILEUNIT = 200_pInt
+ integer(pInt), parameter :: FILEUNIT = 204_pInt
  integer(pInt) :: &
    o, &                                                                                             !< counter in output loop
    ph, &                                                                                            !< counter in phase loop

src/crystallite.f90

@@ -172,9 +172,9 @@ subroutine crystallite_init
    IO_error
  use material
  use config, only: &
+  config_deallocate, &
   config_crystallite, &
-  crystallite_name, &
+  crystallite_name
-  config_deallocate
  use constitutive, only: &
    constitutive_initialFi, &
    constitutive_microstructure                                                                     ! derived (shortcut) quantities of given state

src/debug.f90

@@ -102,7 +102,7 @@ subroutine debug_init
    IO_EOF
 
  implicit none
- integer(pInt), parameter                 :: FILEUNIT    = 300_pInt
+ integer(pInt), parameter                 :: FILEUNIT    = 330_pInt
 
  integer(pInt)                            :: i, what
  integer(pInt), allocatable, dimension(:) :: chunkPos

src/homogenization.f90

@@ -6,9 +6,6 @@
 !--------------------------------------------------------------------------------------------------
 module homogenization
  use prec, only: &
-#ifdef FEM
-   tOutputData, &
-#endif
    pInt, &
    pReal
 
@@ -22,16 +19,8 @@ module homogenization
    materialpoint_P                                                                                  !< first P--K stress of IP
  real(pReal),   dimension(:,:,:,:,:,:), allocatable, public ::  &
    materialpoint_dPdF                                                                               !< tangent of first P--K stress at IP
-#ifdef FEM
- type(tOutputData), dimension(:),       allocatable, public :: &
-   homogOutput
- type(tOutputData), dimension(:,:),     allocatable, public :: &
-   crystalliteOutput, &
-   phaseOutput
-#else
  real(pReal),   dimension(:,:,:),       allocatable, public :: &
    materialpoint_results                                                                            !< results array of material point
-#endif
  integer(pInt),                                      public, protected  :: &
    materialpoint_sizeResults, &
    homogenization_maxSizePostResults, &
@@ -90,20 +79,15 @@ subroutine homogenization_init
    mesh_element, &
    FE_Nips, &
    FE_geomtype
-#ifdef FEM
- use crystallite, only: &
-   crystallite_sizePostResults
-#else
  use constitutive, only: &
    constitutive_plasticity_maxSizePostResults, &
    constitutive_source_maxSizePostResults
  use crystallite, only: &
    crystallite_maxSizePostResults
-#endif
  use config, only: &
-  config_deallocate, &
   material_configFile, &
   material_localFileExt, &
+  config_deallocate, &
   config_homogenization, &
   homogenization_name
  use material
@@ -412,33 +396,6 @@ subroutine homogenization_init
    hydrogenflux_maxSizePostResults   = max(hydrogenflux_maxSizePostResults  ,hydrogenfluxState(p)%sizePostResults)
  enddo
 
-#ifdef FEM
- allocate(homogOutput      (material_Nhomogenization                         ))
- allocate(crystalliteOutput(material_Ncrystallite,  homogenization_maxNgrains))
- allocate(phaseOutput      (material_Nphase,        homogenization_maxNgrains))
- do p = 1, material_Nhomogenization
-   homogOutput(p)%sizeResults = homogState       (p)%sizePostResults + &
-                                thermalState     (p)%sizePostResults + &
-                                damageState      (p)%sizePostResults + &
-                                vacancyfluxState (p)%sizePostResults + &
-                                porosityState    (p)%sizePostResults + &
-                                hydrogenfluxState(p)%sizePostResults
-   homogOutput(p)%sizeIpCells = count(material_homog==p)
-   allocate(homogOutput(p)%output(homogOutput(p)%sizeResults,homogOutput(p)%sizeIpCells))
- enddo
- do p = 1, material_Ncrystallite; do e = 1, homogenization_maxNgrains
-   crystalliteOutput(p,e)%sizeResults = crystallite_sizePostResults(p)
-   crystalliteOutput(p,e)%sizeIpCells = count(microstructure_crystallite(mesh_element(4,:)) == p .and. &
-                                              homogenization_Ngrains    (mesh_element(3,:)) >= e)*mesh_maxNips
-   allocate(crystalliteOutput(p,e)%output(crystalliteOutput(p,e)%sizeResults,crystalliteOutput(p,e)%sizeIpCells))
- enddo; enddo
- do p = 1, material_Nphase; do e = 1, homogenization_maxNgrains
-   phaseOutput(p,e)%sizeResults = plasticState    (p)%sizePostResults + &
-                                  sum(sourceState (p)%p(:)%sizePostResults)
-   phaseOutput(p,e)%sizeIpCells = count(material_phase(e,:,:) == p)
-   allocate(phaseOutput(p,e)%output(phaseOutput(p,e)%sizeResults,phaseOutput(p,e)%sizeIpCells))
- enddo; enddo
-#else
  materialpoint_sizeResults = 1 &                                                                    ! grain count
                            + 1 + homogenization_maxSizePostResults &                                ! homogSize & homogResult
                                + thermal_maxSizePostResults        &
@@ -450,7 +407,6 @@ subroutine homogenization_init
                                                         + 1 + constitutive_plasticity_maxSizePostResults &     ! constitutive size & constitutive results
                                                             + constitutive_source_maxSizePostResults)
  allocate(materialpoint_results(materialpoint_sizeResults,mesh_maxNips,mesh_NcpElems))
-#endif
 
  write(6,'(/,a)')   ' <<<+-  homogenization init  -+>>>'
  write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
@@ -474,9 +430,6 @@ subroutine homogenization_init
    write(6,'(a32,1x,7(i8,1x))')   'materialpoint_requested:        ', shape(materialpoint_requested)
    write(6,'(a32,1x,7(i8,1x))')   'materialpoint_converged:        ', shape(materialpoint_converged)
    write(6,'(a32,1x,7(i8,1x),/)') 'materialpoint_doneAndHappy:     ', shape(materialpoint_doneAndHappy)
-#ifndef FEM
-   write(6,'(a32,1x,7(i8,1x),/)') 'materialpoint_results:          ', shape(materialpoint_results)
-#endif
    write(6,'(a32,1x,7(i8,1x))')   'maxSizePostResults: ', homogenization_maxSizePostResults
  endif
  flush(6)
@@ -905,33 +858,18 @@ subroutine materialpoint_postResults
    mesh_element
  use material, only: &
    mappingHomogenization, &
-#ifdef FEM
-   phaseAt, phasememberAt, &
-   homogenization_maxNgrains, &
-   material_Ncrystallite, &
-   material_Nphase, &
-#else
    homogState, &
    thermalState, &
    damageState, &
    vacancyfluxState, &
    porosityState, &
    hydrogenfluxState, &
-#endif
    plasticState, &
    sourceState, &
    material_phase, &
    homogenization_Ngrains, &
    microstructure_crystallite
-#ifdef FEM
- use constitutive, only: &
-   constitutive_plasticity_maxSizePostResults, &
-   constitutive_source_maxSizePostResults
-#endif
  use crystallite, only: &
-#ifdef FEM
-   crystallite_maxSizePostResults, &
-#endif
    crystallite_sizePostResults, &
    crystallite_postResults
 
@@ -944,55 +882,6 @@ subroutine materialpoint_postResults
    g, &                                                                                             !< grain number
    i, &                                                                                             !< integration point number
    e                                                                                                !< element number
-#ifdef FEM
- integer(pInt) :: &
-   myHomog, &
-   myPhase, &
-   crystalliteCtr(material_Ncrystallite,  homogenization_maxNgrains), &
-   phaseCtr      (material_Nphase,        homogenization_maxNgrains)
- real(pReal), dimension(1+crystallite_maxSizePostResults + &
-                        1+constitutive_plasticity_maxSizePostResults + &
-                          constitutive_source_maxSizePostResults) :: &
-   crystalliteResults
-
-
-
- crystalliteCtr = 0_pInt; phaseCtr = 0_pInt
- elementLooping: do e = FEsolving_execElem(1),FEsolving_execElem(2)
-   myNgrains = homogenization_Ngrains(mesh_element(3,e))
-   myCrystallite = microstructure_crystallite(mesh_element(4,e))
-   IpLooping: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
-     myHomog = mappingHomogenization(2,i,e)
-     thePos =  mappingHomogenization(1,i,e)
-     homogOutput(myHomog)%output(1: &
-                                 homogOutput(myHomog)%sizeResults, &
-                                 thePos) = homogenization_postResults(i,e)
-
-     grainLooping :do g = 1,myNgrains
-       myPhase = phaseAt(g,i,e)
-       crystalliteResults(1:1+crystallite_sizePostResults(myCrystallite) + &
-                            1+plasticState(myPhase)%sizePostResults + &
-                              sum(sourceState(myPhase)%p(:)%sizePostResults)) = crystallite_postResults(g,i,e)
-       if (microstructure_crystallite(mesh_element(4,e)) == myCrystallite .and. &
-           homogenization_Ngrains    (mesh_element(3,e)) >= g) then
-         crystalliteCtr(myCrystallite,g) = crystalliteCtr(myCrystallite,g) + 1_pInt
-         crystalliteOutput(myCrystallite,g)% &
-           output(1:crystalliteOutput(myCrystallite,g)%sizeResults,crystalliteCtr(myCrystallite,g)) = &
-             crystalliteResults(2:1+crystalliteOutput(myCrystallite,g)%sizeResults)
-       endif
-       if (material_phase(g,i,e) == myPhase) then
-         phaseCtr(myPhase,g) = phaseCtr(myPhase,g) + 1_pInt
-         phaseOutput(myPhase,g)% &
-           output(1:phaseOutput(myPhase,g)%sizeResults,phaseCtr(myPhase,g)) = &
-             crystalliteResults(3 + crystalliteOutput(myCrystallite,g)%sizeResults: &
-                                1 + crystalliteOutput(myCrystallite,g)%sizeResults + &
-                                1 + plasticState    (myphase)%sizePostResults + &
-                                    sum(sourceState(myphase)%p(:)%sizePostResults))
-       endif
-     enddo grainLooping
-   enddo IpLooping
- enddo elementLooping
-#else
 
  !$OMP PARALLEL DO PRIVATE(myNgrains,myCrystallite,thePos,theSize)
    elementLooping: do e = FEsolving_execElem(1),FEsolving_execElem(2)
@@ -1028,7 +917,6 @@ subroutine materialpoint_postResults
      enddo IpLooping
    enddo elementLooping
  !$OMP END PARALLEL DO
-#endif
 
 end subroutine materialpoint_postResults
 

src/lattice.f90

@@ -1278,7 +1278,7 @@ subroutine lattice_init
  integer(pInt) :: Nphases
  character(len=65536) :: &
    tag  = ''
- integer(pInt) :: section = 0_pInt,i,p
+ integer(pInt) :: i,p
  real(pReal),  dimension(:), allocatable :: &
    temp, &
    CoverA, &                                                                                        !< c/a ratio for low symmetry type lattice
@@ -1388,9 +1388,9 @@ subroutine lattice_init
    tag = config_phase(p)%getString('trans_lattice_structure',defaultVal=tag)
    select case(trim(tag))
      case('bcc')
-        trans_lattice_structure(section) = LATTICE_bcc_ID
+        trans_lattice_structure(p) = LATTICE_bcc_ID
      case('hex','hexagonal')
-        trans_lattice_structure(section) = LATTICE_hex_ID
+        trans_lattice_structure(p) = LATTICE_hex_ID
    end select
 
    lattice_C66(1,1,p) = config_phase(p)%getFloat('c11',defaultVal=0.0_pReal)

src/material.f90

@@ -16,7 +16,7 @@ module material
    tSourceState, &
    tHomogMapping, &
    tPhaseMapping, &
-   group_scalar, &
+   group_float, &
    group_int
 
  implicit none
@@ -268,7 +268,7 @@ module material
    porosityMapping, &                                                                               !< mapping for porosity state/fields
    hydrogenfluxMapping                                                                              !< mapping for hydrogen conc state/fields
 
- type(group_scalar),  allocatable, dimension(:), public :: &
+ type(group_float),  allocatable, dimension(:), public :: &
    temperature, &                                                                                   !< temperature field
    damage, &                                                                                        !< damage field
    vacancyConc, &                                                                                   !< vacancy conc field
@@ -360,8 +360,7 @@ subroutine material_init()
    homogenization_name, &
    microstructure_name, &
    phase_name, &
-   texture_name, &
+   texture_name
-   config_deallocate
  use mesh, only: &
    mesh_maxNips, &
    mesh_NcpElems, &
@@ -370,7 +369,7 @@ subroutine material_init()
    FE_geomtype
 
  implicit none
- integer(pInt), parameter :: FILEUNIT = 200_pInt
+ integer(pInt), parameter :: FILEUNIT = 210_pInt
  integer(pInt)            :: m,c,h, myDebug, myPhase, myHomog
  integer(pInt) :: &
   g, &                                                                                              !< grain number
@@ -469,7 +468,6 @@ subroutine material_init()
  endif debugOut
 
  call material_populateGrains
- call config_deallocate('material.config/microstructure')
 
  allocate(phaseAt                   (  homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems),source=0_pInt)
  allocate(phasememberAt             (  homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems),source=0_pInt)
@@ -921,8 +919,8 @@ subroutine material_parseTexture
    IO_floatValue, &
    IO_stringValue
  use config, only: &
-   config_texture, &
+   config_deallocate, &
-   config_deallocate
+   config_texture
  use math, only: &
    inRad, &
    math_sampleRandomOri, &
@@ -1093,6 +1091,7 @@ subroutine material_populateGrains
  use config, only: &
    config_homogenization, &
    config_microstructure, &
+   config_deallocate, &
    homogenization_name, &
    microstructure_name 
  use IO, only: &
@@ -1429,6 +1428,7 @@ subroutine material_populateGrains
  deallocate(texture_transformation)
  deallocate(Nelems)
  deallocate(elemsOfHomogMicro)
+ call config_deallocate('material.config/microstructure')
 
 end subroutine material_populateGrains
 

src/math.f90

@@ -12,7 +12,7 @@ module math
 
  implicit none
  private
- real(pReal),    parameter, public :: PI = 3.141592653589793_pReal                                  !< ratio of a circle's circumference to its diameter
+ real(pReal),    parameter, public :: PI = acos(-1.0_pReal)                                         !< ratio of a circle's circumference to its diameter
  real(pReal),    parameter, public :: INDEG = 180.0_pReal/PI                                        !< conversion from radian into degree
  real(pReal),    parameter, public :: INRAD = PI/180.0_pReal                                        !< conversion from degree into radian
  complex(pReal), parameter, public :: TWOPIIMG = (0.0_pReal,2.0_pReal)*(PI,0.0_pReal)               !< Re(0.0), Im(2xPi)

src/mesh.f90

@@ -95,9 +95,11 @@ module mesh
  integer(pInt), dimension(:,:), allocatable, private :: &
    mesh_cellnodeParent                                                                              !< cellnode's parent element ID, cellnode's intra-element ID
 
+#if defined(Marc4DAMASK) || defined(Abaqus)
  integer(pInt), dimension(:,:), allocatable, target, private :: &
    mesh_mapFEtoCPelem, &                                                                            !< [sorted FEid, corresponding CPid]
    mesh_mapFEtoCPnode                                                                               !< [sorted FEid, corresponding CPid]
+#endif
 
  integer(pInt),dimension(:,:,:), allocatable, private :: &
    mesh_cell                                                                                        !< cell connectivity for each element,ip/cell
@@ -402,7 +404,9 @@ module mesh
 
  public :: &
    mesh_init, &
+#if defined(Marc4DAMASK) || defined(Abaqus)
    mesh_FEasCP, &
+#endif
    mesh_build_cellnodes, &
    mesh_build_ipVolumes, &
    mesh_build_ipCoordinates, &
@@ -420,7 +424,6 @@ module mesh
 #ifdef Spectral
    mesh_spectral_getHomogenization, &
    mesh_spectral_count, &
-   mesh_spectral_mapNodesAndElems, &
    mesh_spectral_count_cpSizes, &
    mesh_spectral_build_nodes, &
    mesh_spectral_build_elements, &
@@ -552,8 +555,6 @@ subroutine mesh_init(ip,el)
  if (myDebug) write(6,'(a)') ' Grid partitioned'; flush(6)
  call mesh_spectral_count()
  if (myDebug) write(6,'(a)') ' Counted nodes/elements'; flush(6)
- call mesh_spectral_mapNodesAndElems
- if (myDebug) write(6,'(a)') ' Mapped nodes and elements'; flush(6)
  call mesh_spectral_count_cpSizes
  if (myDebug) write(6,'(a)') ' Built CP statistics'; flush(6)
  call mesh_spectral_build_nodes()
@@ -659,12 +660,16 @@ subroutine mesh_init(ip,el)
 
  allocate(calcMode(mesh_maxNips,mesh_NcpElems))
  calcMode = .false.                                                                                 ! pretend to have collected what first call is asking (F = I)
+#if defined(Marc4DAMASK) || defined(Abaqus)
  calcMode(ip,mesh_FEasCP('elem',el)) = .true.                                                       ! first ip,el needs to be already pingponged to "calc"
-
+#else
+ calcMode(ip,el) = .true.                                                                           ! first ip,el needs to be already pingponged to "calc"
+#endif
 
 end subroutine mesh_init
 
 
+#if defined(Marc4DAMASK) || defined(Abaqus)
 !--------------------------------------------------------------------------------------------------
 !> @brief Gives the FE to CP ID mapping by binary search through lookup array
 !! valid questions (what) are 'elem', 'node'
@@ -713,7 +718,7 @@ integer(pInt) function mesh_FEasCP(what,myID)
  enddo binarySearch
 
 end function mesh_FEasCP
-
+#endif
 
 !--------------------------------------------------------------------------------------------------
 !> @brief Split CP elements into cells.
@@ -1188,24 +1193,6 @@ subroutine mesh_spectral_count()
 end subroutine mesh_spectral_count
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief fake map node from FE ID to internal (consecutive) representation for node and element
-!! Allocates global array 'mesh_mapFEtoCPnode' and 'mesh_mapFEtoCPelem'
-!--------------------------------------------------------------------------------------------------
-subroutine mesh_spectral_mapNodesAndElems
- use math, only: &
-   math_range
-
- implicit none
- allocate (mesh_mapFEtoCPnode(2_pInt,mesh_Nnodes), source = 0_pInt)
- allocate (mesh_mapFEtoCPelem(2_pInt,mesh_NcpElems), source = 0_pInt)
-
- mesh_mapFEtoCPnode = spread(math_range(mesh_Nnodes),1,2)
- mesh_mapFEtoCPelem = spread(math_range(mesh_NcpElems),1,2)
-
-end subroutine mesh_spectral_mapNodesAndElems
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief Gets maximum count of nodes, IPs, IP neighbors, and subNodes among cpElements.
 !! Sets global values 'mesh_maxNnodes', 'mesh_maxNips', 'mesh_maxNipNeighbors',

src/meshFEM.f90

@@ -0,0 +1,356 @@
+!--------------------------------------------------------------------------------------------------
+!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
+!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
+!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
+!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
+!> @brief Driver controlling inner and outer load case looping of the FEM solver
+!> @details doing cutbacking, forwarding in case of restart, reporting statistics, writing
+!> results
+!--------------------------------------------------------------------------------------------------
+module mesh     
+#include <petsc/finclude/petscis.h>
+#include <petsc/finclude/petscdmda.h>
+ use prec, only: pReal, pInt
+
+use PETScdmda
+use PETScis
+
+ implicit none
+ private
+
+ integer(pInt), public, protected :: &
+   mesh_Nboundaries, &
+   mesh_NcpElems, &                                                                                 !< total number of CP elements in mesh
+   mesh_NcpElemsGlobal, &
+   mesh_Nnodes, &                                                                                   !< total number of nodes in mesh
+   mesh_maxNnodes, &                                                                                !< max number of nodes in any CP element
+   mesh_maxNips, &                                                                                  !< max number of IPs in any CP element
+   mesh_maxNipNeighbors, &
+   mesh_Nelems                                                                                      !< total number of elements in mesh
+
+ real(pReal), public, protected :: charLength
+ 
+ integer(pInt), dimension(:,:), allocatable, public, protected :: &
+   mesh_element                                                                                     !< FEid, type(internal representation), material, texture, node indices as CP IDs
+
+ real(pReal), dimension(:,:), allocatable, public :: &
+   mesh_node                                                                                        !< node x,y,z coordinates (after deformation! ONLY FOR MARC!!!)
+ 
+ real(pReal), dimension(:,:), allocatable, public, protected :: &
+   mesh_ipVolume, &                                                                                 !< volume associated with IP (initially!)
+   mesh_node0                                                                                       !< node x,y,z coordinates (initially!)
+ 
+ real(pReal), dimension(:,:,:), allocatable, public :: &
+   mesh_ipCoordinates                                                                               !< IP x,y,z coordinates (after deformation!)
+
+ real(pReal), dimension(:,:,:), allocatable, public, protected :: &
+   mesh_ipArea                                                                                      !< area of interface to neighboring IP (initially!)
+
+ real(pReal),dimension(:,:,:,:), allocatable, public, protected :: & 
+   mesh_ipAreaNormal                                                                                !< area normal of interface to neighboring IP (initially!)
+ 
+ integer(pInt), dimension(:,:,:,:), allocatable, public, protected :: &
+   mesh_ipNeighborhood                                                                              !< 6 or less neighboring IPs as [element_num, IP_index, neighbor_index that points to me]
+
+ logical, dimension(3), public, protected :: mesh_periodicSurface                                   !< flag indicating periodic outer surfaces (used for fluxes)
+
+ DM, public :: geomMesh
+ 
+ integer(pInt), dimension(:), allocatable, public, protected :: &
+   mesh_boundaries
+
+
+ integer(pInt), parameter, public :: &
+   FE_Nelemtypes =  1_pInt, &
+   FE_Ngeomtypes =  1_pInt, &
+   FE_Ncelltypes =  1_pInt, &
+   FE_maxNnodes  =  1_pInt, &
+   FE_maxNips    = 14_pInt
+                      
+ integer(pInt), dimension(FE_Nelemtypes), parameter, public :: FE_geomtype = &                      !< geometry type of particular element type
+ int([1],pInt)
+
+ integer(pInt), dimension(FE_Ngeomtypes), parameter, public  :: FE_celltype = &                     !< cell type that is used by each geometry type
+ int([1],pInt)
+
+ integer(pInt), dimension(FE_Nelemtypes), parameter, public :: FE_Nnodes = &                        !< number of nodes that constitute a specific type of element
+ int([0],pInt) 
+
+ integer(pInt), dimension(FE_Ngeomtypes),            public :: FE_Nips = &                          !< number of IPs in a specific type of element
+ int([0],pInt)
+
+ integer(pInt), dimension(FE_Ncelltypes), parameter, public :: FE_NipNeighbors = &                  !< number of ip neighbors / cell faces in a specific cell type
+ int([6],pInt)
+
+
+ public :: &
+   mesh_init, &
+   mesh_FEM_build_ipVolumes, &
+   mesh_FEM_build_ipCoordinates, &
+   mesh_cellCenterCoordinates
+
+ external :: &
+   DMPlexCreateFromFile, &
+   DMPlexDistribute, &
+   DMPlexCopyCoordinates, &
+   DMGetStratumSize, &
+   DMPlexGetHeightStratum, &
+   DMPlexGetLabelValue, &
+   DMPlexSetLabelValue
+   
+contains
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief initializes the mesh by calling all necessary private routines the mesh module
+!! Order and routines strongly depend on type of solver
+!--------------------------------------------------------------------------------------------------
+subroutine mesh_init(ip,el)
+ use DAMASK_interface
+ use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
+ use IO, only: &
+   IO_timeStamp, &
+   IO_error, &
+   IO_open_file, &
+   IO_stringPos, &
+   IO_intValue, &
+   IO_EOF, &
+   IO_read, &
+   IO_isBlank
+ use debug, only: &
+   debug_e, &
+   debug_i
+ use numerics, only: &
+   usePingPong, &
+   integrationOrder, &
+   worldrank, &
+   worldsize
+ use FEsolving, only: &
+   FEsolving_execElem, &
+   FEsolving_execIP, &
+   calcMode
+ use FEM_Zoo, only: &
+   FEM_Zoo_nQuadrature, &
+   FEM_Zoo_QuadraturePoints
+ 
+ implicit none
+ integer(pInt), parameter :: FILEUNIT = 222_pInt
+ integer(pInt), intent(in) :: el, ip
+ integer(pInt) :: j
+ integer(pInt), allocatable, dimension(:) :: chunkPos
+ integer :: dimPlex
+ character(len=512) :: &
+   line
+ logical :: flag
+ PetscSF :: sf
+ DM :: globalMesh
+ PetscInt :: face, nFaceSets
+ PetscInt, pointer :: pFaceSets(:)
+ IS :: faceSetIS 
+ PetscErrorCode :: ierr
+
+ 
+ write(6,'(/,a)')   ' <<<+-  mesh init  -+>>>'
+ write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
+#include "compilation_info.f90"
+
+ call DMPlexCreateFromFile(PETSC_COMM_WORLD,geometryFile,PETSC_TRUE,globalMesh,ierr)
+ CHKERRQ(ierr)
+ call DMGetDimension(globalMesh,dimPlex,ierr)
+ CHKERRQ(ierr)
+ call DMGetStratumSize(globalMesh,'depth',dimPlex,mesh_NcpElemsGlobal,ierr)
+ CHKERRQ(ierr)
+ call DMGetLabelSize(globalMesh,'Face Sets',mesh_Nboundaries,ierr)
+ CHKERRQ(ierr)
+ call MPI_Bcast(mesh_Nboundaries,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
+ call MPI_Bcast(mesh_NcpElemsGlobal,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
+ call MPI_Bcast(dimPlex,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
+
+ allocate(mesh_boundaries(mesh_Nboundaries), source = 0_pInt)
+ call DMGetLabelSize(globalMesh,'Face Sets',nFaceSets,ierr)
+ CHKERRQ(ierr)
+ call DMGetLabelIdIS(globalMesh,'Face Sets',faceSetIS,ierr)
+ CHKERRQ(ierr)
+ if (nFaceSets > 0) call ISGetIndicesF90(faceSetIS,pFaceSets,ierr)
+ do face = 1, nFaceSets
+   mesh_boundaries(face) = pFaceSets(face)
+ enddo  
+ if (nFaceSets > 0) call ISRestoreIndicesF90(faceSetIS,pFaceSets,ierr)
+ call MPI_Bcast(mesh_boundaries,mesh_Nboundaries,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
+ 
+ if (worldrank == 0) then
+   j = 0
+   flag = .false.
+   call IO_open_file(FILEUNIT,trim(geometryFile))
+   do
+     read(FILEUNIT,'(a512)') line
+     if (trim(line) == IO_EOF) exit                                                                    ! skip empty lines
+     if (trim(line) == '$Elements') then
+       read(FILEUNIT,'(a512)') line
+       read(FILEUNIT,'(a512)') line
+       flag = .true.  
+     endif
+     if (trim(line) == '$EndElements') exit
+     if (flag) then
+       chunkPos = IO_stringPos(line)
+       if (chunkPos(1) == 3+IO_intValue(line,chunkPos,3)+dimPlex+1) then
+         call DMSetLabelValue(globalMesh,'material',j,IO_intValue(line,chunkPos,4),ierr)
+         CHKERRQ(ierr)
+         j = j + 1  
+       endif                                                                                           ! count all identifiers to allocate memory and do sanity check
+     endif
+   enddo
+   close (FILEUNIT)
+ endif
+ 
+ if (worldsize > 1) then
+   call DMPlexDistribute(globalMesh,0,sf,geomMesh,ierr)
+   CHKERRQ(ierr)
+ else 
+   call DMClone(globalMesh,geomMesh,ierr)
+   CHKERRQ(ierr)
+ endif  
+ call DMDestroy(globalMesh,ierr); CHKERRQ(ierr)
+ 
+ call DMGetStratumSize(geomMesh,'depth',dimPlex,mesh_Nelems,ierr)
+ CHKERRQ(ierr)
+ call DMGetStratumSize(geomMesh,'depth',0,mesh_Nnodes,ierr)
+ CHKERRQ(ierr)
+ mesh_NcpElems = mesh_Nelems
+
+ FE_Nips(FE_geomtype(1_pInt)) = FEM_Zoo_nQuadrature(dimPlex,integrationOrder)
+ mesh_maxNnodes = FE_Nnodes(1_pInt)
+ mesh_maxNips = FE_Nips(1_pInt)
+ call mesh_FEM_build_ipCoordinates(dimPlex,FEM_Zoo_QuadraturePoints(dimPlex,integrationOrder)%p)
+ call mesh_FEM_build_ipVolumes(dimPlex)
+ 
+ allocate (mesh_element (4_pInt+mesh_maxNnodes,mesh_NcpElems)); mesh_element = 0_pInt
+ do j = 1, mesh_NcpElems
+   mesh_element( 1,j) = j
+   mesh_element( 2,j) = 1_pInt                                                                      ! elem type
+   mesh_element( 3,j) = 1_pInt                                                                      ! homogenization
+   call DMGetLabelValue(geomMesh,'material',j-1,mesh_element(4,j),ierr)
+   CHKERRQ(ierr)
+ end do 
+
+ if (usePingPong .and. (mesh_Nelems /= mesh_NcpElems)) &
+   call IO_error(600_pInt)                                                                          ! ping-pong must be disabled when having non-DAMASK elements
+ if (debug_e < 1 .or. debug_e > mesh_NcpElems) &
+   call IO_error(602_pInt,ext_msg='element')                                                        ! selected element does not exist
+ if (debug_i < 1 .or. debug_i > FE_Nips(FE_geomtype(mesh_element(2_pInt,debug_e)))) &
+   call IO_error(602_pInt,ext_msg='IP')                                                             ! selected element does not have requested IP
+ 
+ FEsolving_execElem = [ 1_pInt,mesh_NcpElems ]                                                      ! parallel loop bounds set to comprise all DAMASK elements
+ if (allocated(FEsolving_execIP)) deallocate(FEsolving_execIP)
+ allocate(FEsolving_execIP(2_pInt,mesh_NcpElems)); FEsolving_execIP = 1_pInt                        ! parallel loop bounds set to comprise from first IP...
+ forall (j = 1_pInt:mesh_NcpElems) FEsolving_execIP(2,j) = FE_Nips(FE_geomtype(mesh_element(2,j)))  ! ...up to own IP count for each element
+ 
+ if (allocated(calcMode)) deallocate(calcMode)
+ allocate(calcMode(mesh_maxNips,mesh_NcpElems))
+ calcMode = .false.                                                                                 ! pretend to have collected what first call is asking (F = I)
+ calcMode(ip,el) = .true.                                                                           ! first ip,el needs to be already pingponged to "calc"
+
+end subroutine mesh_init
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Calculates cell center coordinates.
+!--------------------------------------------------------------------------------------------------
+pure function mesh_cellCenterCoordinates(ip,el)
+ 
+ implicit none
+ integer(pInt), intent(in) :: el, &                                                                  !< element number
+                              ip                                                                     !< integration point number
+ real(pReal), dimension(3) :: mesh_cellCenterCoordinates                                             !< x,y,z coordinates of the cell center of the requested IP cell
+
+ end function mesh_cellCenterCoordinates
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Calculates IP volume. Allocates global array 'mesh_ipVolume'
+!> @details The IP volume is calculated differently depending on the cell type.
+!> 2D cells assume an element depth of one in order to calculate the volume.
+!> For the hexahedral cell we subdivide the cell into subvolumes of pyramidal
+!> shape with a cell face as basis and the central ip at the tip. This subvolume is
+!> calculated as an average of four tetrahedals with three corners on the cell face 
+!> and one corner at the central ip.
+!--------------------------------------------------------------------------------------------------
+subroutine mesh_FEM_build_ipVolumes(dimPlex)
+ use math, only: &
+   math_I3, &
+   math_det33
+ 
+ implicit none
+ PetscInt           :: dimPlex
+ PetscReal          :: vol
+ PetscReal,  target :: cent(dimPlex), norm(dimPlex)
+ PetscReal, pointer :: pCent(:), pNorm(:)
+ PetscInt           :: cellStart, cellEnd, cell
+ PetscErrorCode     :: ierr
+
+ if (.not. allocated(mesh_ipVolume)) then
+   allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems))
+   mesh_ipVolume = 0.0_pReal 
+ endif
+
+ call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
+ pCent => cent
+ pNorm => norm
+ do cell = cellStart, cellEnd-1
+   call  DMPlexComputeCellGeometryFVM(geomMesh,cell,vol,pCent,pNorm,ierr) 
+   CHKERRQ(ierr)
+   mesh_ipVolume(:,cell+1) = vol/real(mesh_maxNips,pReal)
+ enddo  
+
+end subroutine mesh_FEM_build_ipVolumes
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Calculates IP Coordinates. Allocates global array 'mesh_ipCoordinates'
+! Called by all solvers in mesh_init in order to initialize the ip coordinates.
+! Later on the current ip coordinates are directly prvided by the spectral solver and by Abaqus,
+! so no need to use this subroutine anymore; Marc however only provides nodal displacements,
+! so in this case the ip coordinates are always calculated on the basis of this subroutine.
+! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! FOR THE MOMENT THIS SUBROUTINE ACTUALLY CALCULATES THE CELL CENTER AND NOT THE IP COORDINATES,
+! AS THE IP IS NOT (ALWAYS) LOCATED IN THE CENTER OF THE IP VOLUME. 
+! HAS TO BE CHANGED IN A LATER VERSION. 
+! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!--------------------------------------------------------------------------------------------------
+subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints)
+
+ implicit none
+ PetscInt,      intent(in) :: dimPlex
+ PetscReal,     intent(in) :: qPoints(mesh_maxNips*dimPlex)
+ PetscReal,         target :: v0(dimPlex), cellJ(dimPlex*dimPlex), invcellJ(dimPlex*dimPlex)
+ PetscReal,        pointer :: pV0(:), pCellJ(:), pInvcellJ(:)
+ PetscReal                 :: detJ
+ PetscInt                  :: cellStart, cellEnd, cell, qPt, dirI, dirJ, qOffset
+ PetscErrorCode            :: ierr
+
+
+ allocate(mesh_ipCoordinates(3,mesh_maxNips,mesh_NcpElems),source=0.0_pReal)
+
+ pV0 => v0
+ pCellJ => cellJ
+ pInvcellJ => invcellJ
+ call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
+ do cell = cellStart, cellEnd-1                                                                     !< loop over all elements 
+   call DMPlexComputeCellGeometryAffineFEM(geomMesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
+   CHKERRQ(ierr)
+   qOffset = 0
+   do qPt = 1, mesh_maxNips
+     do dirI = 1, dimPlex
+       mesh_ipCoordinates(dirI,qPt,cell+1) = pV0(dirI)
+       do dirJ = 1, dimPlex
+         mesh_ipCoordinates(dirI,qPt,cell+1) = mesh_ipCoordinates(dirI,qPt,cell+1) + &
+                                               pCellJ((dirI-1)*dimPlex+dirJ)*(qPoints(qOffset+dirJ) + 1.0)
+       enddo
+     enddo
+     qOffset = qOffset + dimPlex
+   enddo
+ enddo  
+
+end subroutine mesh_FEM_build_ipCoordinates
+
+end module mesh

src/plastic_isotropic.f90

@@ -109,11 +109,9 @@ use IO
  type(tParameters), pointer :: prm
  
  integer(pInt) :: &
-   o, &
    phase, & 
    instance, &
    maxNinstance, &
-   mySize, &
    sizeDotState, &
    sizeState, &
    sizeDeltaState
@@ -136,7 +134,6 @@ use IO
           plastic_isotropic_output = ''
  allocate(plastic_isotropic_Noutput(maxNinstance),                              source=0_pInt)
 
-! inernal variable 
  allocate(param(maxNinstance))                                                                      ! one container of parameters per instance
  allocate(state(maxNinstance))                                                                      ! internal state aliases
  allocate(dotState(maxNinstance))

src/plastic_phenopowerlaw.f90

@@ -241,29 +241,29 @@ subroutine plastic_phenopowerlaw_init
      select case(outputs(i))
         case ('resistance_slip')
           outputID = resistance_slip_ID
-          outputSize = sum(prm%Nslip)
+          outputSize = prm%totalNslip
         case ('accumulatedshear_slip')
           outputID = accumulatedshear_slip_ID
-          outputSize = sum(prm%Nslip)
+          outputSize = prm%totalNslip
         case ('shearrate_slip')
           outputID = shearrate_slip_ID
-          outputSize = sum(prm%Nslip)
+          outputSize = prm%totalNslip
         case ('resolvedstress_slip')
           outputID = resolvedstress_slip_ID
-          outputSize = sum(prm%Nslip)
+          outputSize = prm%totalNslip
 
         case ('resistance_twin')
           outputID = resistance_twin_ID
-          outputSize = sum(prm%Ntwin)
+          outputSize = prm%totalNtwin
         case ('accumulatedshear_twin')
           outputID = accumulatedshear_twin_ID
-          outputSize = sum(prm%Ntwin)
+          outputSize = prm%totalNtwin
         case ('shearrate_twin')
           outputID = shearrate_twin_ID
-          outputSize = sum(prm%Ntwin)
+          outputSize = prm%totalNtwin
         case ('resolvedstress_twin')
           outputID = resolvedstress_twin_ID
-          outputSize = sum(prm%Ntwin)
+          outputSize = prm%totalNtwin
 
         case ('totalvolfrac_twin')
           outputID = totalvolfrac_twin_ID

src/prec.f90

@@ -7,6 +7,7 @@
 !> @brief    setting precision for real and int type
 !--------------------------------------------------------------------------------------------------
 module prec
+! ToDo:    use, intrinsic :: iso_fortran_env, only : I8 => int64, WP => real64 
  implicit none
  private 
 #if (FLOAT==8)
@@ -23,14 +24,15 @@ module prec
  NO SUITABLE PRECISION FOR INTEGER SELECTED, STOPPING COMPILATION
 #endif
 
+ integer,     parameter, public :: pStringLen = 256                                                 !< default string lenth
  integer,     parameter, public :: pLongInt  = 8                                                    !< integer representation 64 bit (was selected_int_kind(12), number with at least up to +- 1e12)
  real(pReal), parameter, public :: tol_math_check = 1.0e-8_pReal                                    !< tolerance for internal math self-checks (rotation)
 
  integer(pInt), allocatable, dimension(:) :: realloc_lhs_test
  
- type, public :: group_scalar                                                                       !< variable length datatype used for storage of state
+ type, public :: group_float                                                                       !< variable length datatype used for storage of state
    real(pReal), dimension(:), pointer :: p
- end type group_scalar
+ end type group_float
 
  type, public :: group_int
    integer(pInt), dimension(:), pointer :: p

src/system_routines.f90

@@ -78,28 +78,30 @@ end function isDirectory
 !--------------------------------------------------------------------------------------------------
 !> @brief gets the current working directory
 !--------------------------------------------------------------------------------------------------
-logical function getCWD(str)
+character(len=1024) function getCWD()
   use, intrinsic :: ISO_C_Binding, only: &
     C_INT, &
     C_CHAR, &
     C_NULL_CHAR
 
   implicit none
-  character(len=*), intent(out) :: str
+  character(kind=C_CHAR), dimension(1024) :: charArray                                              ! C string is an array
-  character(kind=C_CHAR), dimension(1024) :: strFixedLength                                         ! C string is an array
   integer(C_INT) :: stat
   integer :: i
 
-  str = repeat('',len(str))
+  call getCurrentWorkDir_C(charArray,stat)
-  call getCurrentWorkDir_C(strFixedLength,stat)
+  if (stat /= 0_C_INT) then
-  do i=1,1024                                                                                       ! copy array components until Null string is found
+    getCWD = 'Error occured when getting currend working directory'
-    if (strFixedLength(i) /= C_NULL_CHAR) then
+  else
-      str(i:i)=strFixedLength(i)
+    getCWD = repeat('',len(getCWD))
+    arrayToString: do i=1,len(getCWD)
+      if (charArray(i) /= C_NULL_CHAR) then
+        getCWD(i:i)=charArray(i)
       else
         exit
       endif
-  enddo
+    enddo arrayToString
-  getCWD=merge(.True.,.False.,stat /= 0_C_INT)
+  endif
 
 end function getCWD
 
@@ -107,28 +109,30 @@ end function getCWD
 !--------------------------------------------------------------------------------------------------
 !> @brief gets the current host name
 !--------------------------------------------------------------------------------------------------
-logical function getHostName(str)
+character(len=1024) function getHostName()
   use, intrinsic :: ISO_C_Binding, only: &
     C_INT, &
     C_CHAR, &
     C_NULL_CHAR
 
   implicit none
-  character(len=*), intent(out) :: str
+  character(kind=C_CHAR), dimension(1024) :: charArray                                              ! C string is an array
-  character(kind=C_CHAR), dimension(1024) :: strFixedLength                                         ! C string is an array
   integer(C_INT) :: stat
   integer :: i
 
-  str = repeat('',len(str))
+  call getHostName_C(charArray,stat)
-  call getHostName_C(strFixedLength,stat)
+  if (stat /= 0_C_INT) then
-  do i=1,1024                                                                                       ! copy array components until Null string is found
+    getHostName = 'Error occured when getting host name'
-    if (strFixedLength(i) /= C_NULL_CHAR) then
+  else
-      str(i:i)=strFixedLength(i)
+    getHostName = repeat('',len(getHostName))
+    arrayToString: do i=1,len(getHostName)
+      if (charArray(i) /= C_NULL_CHAR) then
+        getHostName(i:i)=charArray(i)
       else
         exit
       endif
-  enddo
+    enddo arrayToString
-  getHostName=merge(.True.,.False.,stat /= 0_C_INT)
+  endif
 
 end function getHostName
 

src/vacancyflux_cahnhilliard.f90

@@ -7,7 +7,7 @@ module vacancyflux_cahnhilliard
  use prec, only: &
    pReal, &
    pInt, &
-   group_scalar
+   group_float
 
  implicit none
  private
@@ -26,7 +26,7 @@ module vacancyflux_cahnhilliard
  real(pReal),                         dimension(:),           allocatable,        private :: &
    vacancyflux_cahnhilliard_flucAmplitude
 
- type(group_scalar),                  dimension(:),           allocatable,        private :: &
+ type(group_float),                  dimension(:),           allocatable,        private :: &
    vacancyflux_cahnhilliard_thermalFluc
 
  real(pReal),                                                 parameter,           private :: &