Merge commit 'v2.0.2-443-gd765cf28'

.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 81fd7109fea8456b8eecaaef0eec041edcce7792
+Subproject commit a764ade044735df35fac93a5204446291ee29abc

VERSION

@@ -1 +1 @@
-v2.0.2-409-gac011684
+v2.0.2-442-gb11666ef

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/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,16 +79,11 @@ 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, &
@@ -411,33 +395,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        &
@@ -449,7 +406,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()
@@ -473,9 +429,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)
@@ -904,33 +857,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
 
@@ -943,55 +881,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)
@@ -1027,7 +916,6 @@ subroutine materialpoint_postResults
      enddo IpLooping
    enddo elementLooping
  !$OMP END PARALLEL DO
-#endif
 
 end subroutine materialpoint_postResults
 

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

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

@@ -28,9 +28,9 @@ module prec
 
  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 :: &