Merge branch 'separating-physics' into 'development'

separate physics

See merge request damask/DAMASK!413

examples/config/phase/damage/isobrittle_generic.yaml

@@ -1,7 +1,5 @@
 type: isobrittle
 W_crit: 1400000.0
-m: 1.0
-isoBrittle_atol: 0.01
 
 output: [f_phi]
 

src/CPFEM.f90

@@ -189,9 +189,11 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
       CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_eye(6)
 
     else validCalculation
-      if (debugCPFEM%extensive) &
+      if (debugCPFEM%extensive)  print'(a,i8,1x,i2)', '<< CPFEM >> calculation for elFE ip ',elFE,ip
-        print'(a,i8,1x,i2)', '<< CPFEM >> calculation for elFE ip ',elFE,ip
+      call homogenization_mechanical_response(dt,[ip,ip],[elCP,elCP])
-      call materialpoint_stressAndItsTangent(dt,[ip,ip],[elCP,elCP])
+      if (.not. terminallyIll) &
+        call homogenization_mechanical_response2(dt,[ip,ip],[elCP,elCP])
+
 
       terminalIllness: if (terminallyIll) then
 

src/grid/grid_mech_FEM.f90

@@ -261,7 +261,7 @@ subroutine grid_mechanical_FEM_init
     F         = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3)
   endif restartRead
 
-  homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3])                            ! set starting condition for materialpoint_stressAndItsTangent
+  homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3])                            ! set starting condition for homogenization_mechanical_response
   call utilities_updateCoords(F)
   call utilities_constitutiveResponse(P_current,P_av,C_volAvg,devNull, &                            ! stress field, stress avg, global average of stiffness and (min+max)/2
                                       F, &                                                          ! target F
@@ -490,8 +490,8 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,i
     BCTol
 
   err_div = fnorm*sqrt(wgt)*geomSize(1)/scaledGeomSize(1)/detJ
-  divTol = max(maxval(abs(P_av))*num%eps_div_rtol   ,num%eps_div_atol)
+  divTol = max(maxval(abs(P_av))*num%eps_div_rtol, num%eps_div_atol)
-  BCTol  = max(maxval(abs(P_av))*num%eps_stress_rtol,num%eps_stress_atol)
+  BCTol  = max(maxval(abs(P_av))*num%eps_stress_rtol, num%eps_stress_atol)
 
   if ((totalIter >= num%itmin .and. all([err_div/divTol, err_BC/BCTol] < 1.0_pReal)) &
        .or. terminallyIll) then
@@ -502,8 +502,6 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,i
     reason = 0
   endif
 
-!--------------------------------------------------------------------------------------------------
-! report
   print'(1/,a)', ' ... reporting .............................................................'
   print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', &
           err_div/divTol,  ' (',err_div,' / m, tol = ',divTol,')'

src/grid/grid_mech_spectral_basic.f90

@@ -211,7 +211,7 @@ subroutine grid_mechanical_spectral_basic_init
     F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
   endif restartRead
 
-  homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3])                            ! set starting condition for materialpoint_stressAndItsTangent
+  homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3])                            ! set starting condition for homogenization_mechanical_response
   call utilities_updateCoords(reshape(F,shape(F_lastInc)))
   call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, &                                ! stress field, stress avg, global average of stiffness and (min+max)/2
                                       reshape(F,shape(F_lastInc)), &                                ! target F
@@ -429,8 +429,8 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
     divTol, &
     BCTol
 
-  divTol = max(maxval(abs(P_av))*num%eps_div_rtol   ,num%eps_div_atol)
+  divTol = max(maxval(abs(P_av))*num%eps_div_rtol, num%eps_div_atol)
-  BCTol  = max(maxval(abs(P_av))*num%eps_stress_rtol,num%eps_stress_atol)
+  BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol, num%eps_stress_atol)
 
   if ((totalIter >= num%itmin .and. all([err_div/divTol, err_BC/BCTol] < 1.0_pReal)) &
        .or. terminallyIll) then
@@ -441,8 +441,6 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
     reason = 0
   endif
 
-!--------------------------------------------------------------------------------------------------
-! report
   print'(1/,a)', ' ... reporting .............................................................'
   print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', &
           err_div/divTol,  ' (',err_div,' / m, tol = ',divTol,')'

src/grid/grid_mech_spectral_polarisation.f90

@@ -237,7 +237,7 @@ subroutine grid_mechanical_spectral_polarisation_init
     F_tau_lastInc = 2.0_pReal*F_lastInc
   endif restartRead
 
-  homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3])                            ! set starting condition for materialpoint_stressAndItsTangent
+  homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3])                            ! set starting condition for homogenization_mechanical_response
   call utilities_updateCoords(reshape(F,shape(F_lastInc)))
   call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, &                                ! stress field, stress avg, global average of stiffness and (min+max)/2
                                       reshape(F,shape(F_lastInc)), &                                ! target F
@@ -487,9 +487,9 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
     divTol, &
     BCTol
 
-  curlTol    = max(maxval(abs(F_aim-math_I3))*num%eps_curl_rtol  ,num%eps_curl_atol)
+  curlTol = max(maxval(abs(F_aim-math_I3))*num%eps_curl_rtol, num%eps_curl_atol)
-  divTol     = max(maxval(abs(P_av))         *num%eps_div_rtol   ,num%eps_div_atol)
+  divTol = max(maxval(abs(P_av))*num%eps_div_rtol, num%eps_div_atol)
-  BCTol      = max(maxval(abs(P_av))         *num%eps_stress_rtol,num%eps_stress_atol)
+  BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol, num%eps_stress_atol)
 
   if ((totalIter >= num%itmin .and. all([err_div/divTol, err_curl/curlTol, err_BC/BCTol] < 1.0_pReal)) &
        .or. terminallyIll) then
@@ -500,8 +500,6 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
     reason = 0
   endif
 
-!--------------------------------------------------------------------------------------------------
-! report
   print'(1/,a)', ' ... reporting .............................................................'
   print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', &
             err_div/divTol,  ' (',err_div, ' / m, tol = ',divTol,')'
@@ -542,13 +540,13 @@ subroutine formResidual(in, FandF_tau, &
 !---------------------------------------------------------------------------------------------------
 
   F              => FandF_tau(1:3,1:3,1,&
-                                 XG_RANGE,YG_RANGE,ZG_RANGE)
+                              XG_RANGE,YG_RANGE,ZG_RANGE)
   F_tau          => FandF_tau(1:3,1:3,2,&
-                                 XG_RANGE,YG_RANGE,ZG_RANGE)
+                              XG_RANGE,YG_RANGE,ZG_RANGE)
   residual_F     => residuum(1:3,1:3,1,&
-                                 X_RANGE, Y_RANGE, Z_RANGE)
+                             X_RANGE, Y_RANGE, Z_RANGE)
   residual_F_tau => residuum(1:3,1:3,2,&
-                                 X_RANGE, Y_RANGE, Z_RANGE)
+                             X_RANGE, Y_RANGE, Z_RANGE)
 
   F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt
   call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr)

src/grid/spectral_utilities.f90

@@ -815,10 +815,14 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
 
   homogenization_F  = reshape(F,[3,3,product(grid(1:2))*grid3])                                     ! set materialpoint target F to estimated field
 
-  call materialpoint_stressAndItsTangent(timeinc,[1,1],[1,product(grid(1:2))*grid3])                ! calculate P field
+  call homogenization_mechanical_response(timeinc,[1,1],[1,product(grid(1:2))*grid3])               ! calculate P field
+  if (.not. terminallyIll) &
+    call homogenization_thermal_response(timeinc,[1,1],[1,product(grid(1:2))*grid3])
+  if (.not. terminallyIll) &
+    call homogenization_mechanical_response2(timeinc,[1,1],[1,product(grid(1:2))*grid3])
 
   P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3])
-  P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt                                                   ! average of P
+  P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt
   call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr)
   if (debugRotation) print'(/,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &
     ' Piola--Kirchhoff stress (lab) / MPa =', transpose(P_av)*1.e-6_pReal

src/homogenization.f90

@@ -101,8 +101,8 @@ module homogenization
       integer,     intent(in) :: ce
     end subroutine damage_partition
 
-    module subroutine mechanical_homogenize(dt,ce)
+    module subroutine mechanical_homogenize(Delta_t,ce)
-     real(pReal), intent(in) :: dt
+     real(pReal), intent(in) :: Delta_t
      integer, intent(in) :: &
        ce                                                                                           !< cell
     end subroutine mechanical_homogenize
@@ -178,7 +178,9 @@ module homogenization
 
   public ::  &
     homogenization_init, &
-    materialpoint_stressAndItsTangent, &
+    homogenization_mechanical_response, &
+    homogenization_mechanical_response2, &
+    homogenization_thermal_response, &
     homogenization_mu_T, &
     homogenization_K_T, &
     homogenization_f_T, &
@@ -227,24 +229,24 @@ end subroutine homogenization_init
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief  parallelized calculation of stress and corresponding tangent at material points
+!> @brief
 !--------------------------------------------------------------------------------------------------
-subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execElem)
+subroutine homogenization_mechanical_response(Delta_t,FEsolving_execIP,FEsolving_execElem)
 
-  real(pReal), intent(in) :: dt                                                                     !< time increment
+  real(pReal), intent(in) :: Delta_t                                                                !< time increment
   integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP
   integer :: &
     NiterationMPstate, &
     ip, &                                                                                            !< integration point number
     el, &                                                                                            !< element number
-    co, ce, ho, en, ph
+    co, ce, ho, en
   logical :: &
     converged
   logical, dimension(2) :: &
     doneAndHappy
 
-  !$OMP PARALLEL
+
-  !$OMP DO PRIVATE(ce,en,ho,NiterationMPstate,converged,doneAndHappy)
+  !$OMP PARALLEL DO PRIVATE(ce,en,ho,NiterationMPstate,converged,doneAndHappy)
   do el = FEsolving_execElem(1),FEsolving_execElem(2)
 
     do ip = FEsolving_execIP(1),FEsolving_execIP(2)
@@ -266,65 +268,89 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
         NiterationMPstate = NiterationMPstate + 1
 
         call mechanical_partition(homogenization_F(1:3,1:3,ce),ce)
-        converged = .true.
+        converged = all([(phase_mechanical_constitutive(Delta_t,co,ip,el),co=1,homogenization_Nconstituents(ho))])
-        do co = 1, homogenization_Nconstituents(ho)
-          converged = converged .and. crystallite_stress(dt,co,ip,el)
-        enddo
-
         if (converged) then
-          doneAndHappy = mechanical_updateState(dt,homogenization_F(1:3,1:3,ce),ce)
+          doneAndHappy = mechanical_updateState(Delta_t,homogenization_F(1:3,1:3,ce),ce)
           converged = all(doneAndHappy)
         else
           doneAndHappy = [.true.,.false.]
         endif
-
       enddo convergenceLooping
+
+      converged = converged .and. all([(phase_damage_constitutive(Delta_t,co,ip,el),co=1,homogenization_Nconstituents(ho))])
+
       if (.not. converged) then
-        if (.not. terminallyIll) print*, ' Integration point ', ip,' at element ', el, ' terminally ill'
+        if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill'
         terminallyIll = .true.
       endif
     enddo
   enddo
-  !$OMP END DO
+  !$OMP END PARALLEL DO
 
-  if (.not. terminallyIll) then
+end subroutine homogenization_mechanical_response
-    !$OMP DO PRIVATE(ho,ph,ce)
+
-    do el = FEsolving_execElem(1),FEsolving_execElem(2)
+
-      if (terminallyIll) continue
+!--------------------------------------------------------------------------------------------------
-      do ip = FEsolving_execIP(1),FEsolving_execIP(2)
+!> @brief
-        ce = (el-1)*discretization_nIPs + ip
+!--------------------------------------------------------------------------------------------------
-        ho = material_homogenizationID(ce)
+subroutine homogenization_thermal_response(Delta_t,FEsolving_execIP,FEsolving_execElem)
-        call thermal_partition(ce)
+
-        do co = 1, homogenization_Nconstituents(ho)
+  real(pReal), intent(in) :: Delta_t                                                                !< time increment
-          ph = material_phaseID(co,ce)
+  integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP
-          if (.not. thermal_stress(dt,ph,material_phaseMemberAt(co,ip,el))) then
+  integer :: &
-            if (.not. terminallyIll) &                                                              ! so first signals terminally ill...
+    ip, &                                                                                           !< integration point number
-              print*, ' Integration point ', ip,' at element ', el, ' terminally ill'
+    el, &                                                                                           !< element number
-            terminallyIll = .true.                                                                  ! ...and kills all others
+    co, ce, ho
-          endif
+
-        enddo
+
+  !$OMP PARALLEL DO PRIVATE(ho,ce)
+  do el = FEsolving_execElem(1),FEsolving_execElem(2)
+    if (terminallyIll) continue
+    do ip = FEsolving_execIP(1),FEsolving_execIP(2)
+      ce = (el-1)*discretization_nIPs + ip
+      ho = material_homogenizationID(ce)
+      call thermal_partition(ce)
+      do co = 1, homogenization_Nconstituents(ho)
+        if (.not. phase_thermal_constitutive(Delta_t,material_phaseID(co,ce),material_phaseEntry(co,ce))) then
+          if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill'
+          terminallyIll = .true.
+        endif
       enddo
     enddo
-    !$OMP END DO
+  enddo
+  !$OMP END PARALLEL DO
 
-    !$OMP DO PRIVATE(ho,ce)
+end subroutine homogenization_thermal_response
-    elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2)
-      IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2)
-        ce = (el-1)*discretization_nIPs + ip
-        ho = material_homogenizationID(ce)
-        do co = 1, homogenization_Nconstituents(ho)
-          call crystallite_orientations(co,ip,el)
-        enddo
-        call mechanical_homogenize(dt,ce)
-      enddo IpLooping3
-    enddo elementLooping3
-    !$OMP END DO
-  else
-    print'(/,a,/)', ' << HOMOG >> Material Point terminally ill'
-  endif
-  !$OMP END PARALLEL
 
-end subroutine materialpoint_stressAndItsTangent
+
+!--------------------------------------------------------------------------------------------------
+!> @brief
+!--------------------------------------------------------------------------------------------------
+subroutine homogenization_mechanical_response2(Delta_t,FEsolving_execIP,FEsolving_execElem)
+
+  real(pReal), intent(in) :: Delta_t                                                                !< time increment
+  integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP
+  integer :: &
+    ip, &                                                                                           !< integration point number
+    el, &                                                                                           !< element number
+    co, ce, ho
+
+
+  !$OMP PARALLEL DO PRIVATE(ho,ce)
+  elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2)
+    IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2)
+      ce = (el-1)*discretization_nIPs + ip
+      ho = material_homogenizationID(ce)
+      do co = 1, homogenization_Nconstituents(ho)
+        call crystallite_orientations(co,ip,el)
+      enddo
+      call mechanical_homogenize(Delta_t,ce)
+    enddo IpLooping3
+  enddo elementLooping3
+  !$OMP END PARALLEL DO
+
+
+end subroutine homogenization_mechanical_response2
 
 
 !--------------------------------------------------------------------------------------------------

src/homogenization_mechanical.f90

@@ -123,21 +123,21 @@ end subroutine mechanical_partition
 !--------------------------------------------------------------------------------------------------
 !> @brief Average P and dPdF from the individual constituents.
 !--------------------------------------------------------------------------------------------------
-module subroutine mechanical_homogenize(dt,ce)
+module subroutine mechanical_homogenize(Delta_t,ce)
 
-  real(pReal), intent(in) :: dt
+  real(pReal), intent(in) :: Delta_t
   integer, intent(in) :: ce
 
   integer :: co
 
 
   homogenization_P(1:3,1:3,ce)            = phase_P(1,ce)
-  homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(dt,1,ce)
+  homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(Delta_t,1,ce)
   do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
     homogenization_P(1:3,1:3,ce)            = homogenization_P(1:3,1:3,ce) &
                                             + phase_P(co,ce)
     homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = homogenization_dPdF(1:3,1:3,1:3,1:3,ce) &
-                                            + phase_mechanical_dPdF(dt,co,ce)
+                                            + phase_mechanical_dPdF(Delta_t,co,ce)
   enddo
 
   homogenization_P(1:3,1:3,ce)            = homogenization_P(1:3,1:3,ce) &

src/mesh/FEM_utilities.f90

@@ -27,17 +27,17 @@ module FEM_utilities
   logical, public             :: cutBack = .false.                                                  !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
   integer, public, parameter  :: maxFields = 6
   integer, public             :: nActiveFields = 0
- 
+
 !--------------------------------------------------------------------------------------------------
 ! grid related information information
   real(pReal),   public       :: wgt                                                                !< weighting factor 1/Nelems
-  
+
 
 !--------------------------------------------------------------------------------------------------
 ! field labels information
   character(len=*),                         parameter,            public :: &
     FIELD_MECH_label     = 'mechanical'
- 
+
   enum, bind(c); enumerator :: &
     FIELD_UNDEFINED_ID, &
     FIELD_MECH_ID
@@ -48,7 +48,7 @@ module FEM_utilities
     COMPONENT_MECH_Y_ID, &
     COMPONENT_MECH_Z_ID
   end enum
- 
+
 !--------------------------------------------------------------------------------------------------
 ! variables controlling debugging
  logical :: &
@@ -57,23 +57,23 @@ module FEM_utilities
 !--------------------------------------------------------------------------------------------------
 ! derived types
   type, public :: tSolutionState                                                                    !< return type of solution from FEM solver variants
-    logical :: converged        = .true.   
+    logical :: converged        = .true.
-    logical :: stagConverged    = .true.   
+    logical :: stagConverged    = .true.
     integer :: iterationsNeeded = 0
   end type tSolutionState
- 
+
   type, public :: tComponentBC
     integer(kind(COMPONENT_UNDEFINED_ID))          :: ID
     real(pReal),                       allocatable, dimension(:) :: Value
-    logical,                           allocatable, dimension(:) :: Mask 
+    logical,                           allocatable, dimension(:) :: Mask
   end type tComponentBC
- 
+
   type, public :: tFieldBC
     integer(kind(FIELD_UNDEFINED_ID))  :: ID
     integer                            :: nComponents = 0
     type(tComponentBC),    allocatable :: componentBC(:)
   end type tFieldBC
- 
+
   type, public :: tLoadCase
     real(pReal)  :: time                   = 0.0_pReal                                              !< length of increment
     integer      :: incs                   = 0, &                                                   !< number of increments
@@ -83,7 +83,7 @@ module FEM_utilities
     integer,        allocatable, dimension(:) :: faceID
     type(tFieldBC), allocatable, dimension(:) :: fieldBC
   end type tLoadCase
-  
+
   public :: &
     FEM_utilities_init, &
     utilities_constitutiveResponse, &
@@ -94,14 +94,14 @@ module FEM_utilities
     COMPONENT_MECH_Y_ID, &
     COMPONENT_MECH_Z_ID
 
-contains 
+contains
 
 !ToDo: use functions in variable call
 !--------------------------------------------------------------------------------------------------
 !> @brief allocates all neccessary fields, sets debug flags
 !--------------------------------------------------------------------------------------------------
 subroutine FEM_utilities_init
-   
+
   character(len=pStringLen) :: petsc_optionsOrder
   class(tNode), pointer :: &
     num_mesh, &
@@ -113,7 +113,7 @@ subroutine FEM_utilities_init
   PetscErrorCode            :: ierr
 
   print'(/,a)',   ' <<<+-  FEM_utilities init  -+>>>'
- 
+
   num_mesh    => config_numerics%get('mesh',defaultVal=emptyDict)
   structOrder =  num_mesh%get_asInt('structOrder', defaultVal = 2)
 
@@ -141,7 +141,7 @@ subroutine FEM_utilities_init
   write(petsc_optionsOrder,'(a,i0)') '-mechFE_petscspace_degree ', structOrder
   call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsOrder),ierr)
   CHKERRQ(ierr)
-  
+
   wgt = 1.0/real(mesh_maxNips*mesh_NcpElemsGlobal,pReal)
 
 
@@ -152,20 +152,21 @@ end subroutine FEM_utilities_init
 !> @brief calculates constitutive response
 !--------------------------------------------------------------------------------------------------
 subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
- 
+
   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
-  
+
   PetscErrorCode :: ierr
 
   print'(/,a)', ' ... evaluating constitutive response ......................................'
 
-  call materialpoint_stressAndItsTangent(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems])                ! calculate P field
+  call homogenization_mechanical_response(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems])                ! calculate P field
+  if (.not. terminallyIll) &
+    call homogenization_mechanical_response2(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems])
+  cutBack = .false.
 
-  cutBack = .false.                                                                                 ! reset cutBack status
-  
   P_av = sum(homogenization_P,dim=3) * wgt
   call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr)
 
@@ -198,12 +199,12 @@ subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCVa
     do dof = offset+comp+1, offset+numDof, numComp
       localArray(dof) = localArray(dof) + BCValue + BCDotValue*timeinc
     enddo
-  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
 
 end module FEM_utilities

src/phase.f90

@@ -108,9 +108,13 @@ module phase
       logical, intent(in) :: includeL
     end subroutine mechanical_restore
 
+    module subroutine damage_restore(ce)
+      integer, intent(in) :: ce
+    end subroutine damage_restore
 
-    module function phase_mechanical_dPdF(dt,co,ce) result(dPdF)
+
-      real(pReal), intent(in) :: dt
+    module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
+      real(pReal), intent(in) :: Delta_t
       integer, intent(in) :: &
         co, &                                                                                       !< counter in constituent loop
         ce
@@ -164,8 +168,8 @@ module phase
       real(pReal) :: dot_T
     end function thermal_dot_T
 
-    module function damage_phi(ph,me) result(phi)
+    module function damage_phi(ph,en) result(phi)
-      integer, intent(in) :: ph,me
+      integer, intent(in) :: ph,en
       real(pReal) :: phi
     end function damage_phi
 
@@ -209,29 +213,27 @@ module phase
 
 ! == cleaned:end ===================================================================================
 
-    module function thermal_stress(Delta_t,ph,en) result(converged_)
+    module function phase_thermal_constitutive(Delta_t,ph,en) result(converged_)
 
       real(pReal), intent(in) :: Delta_t
       integer, intent(in) :: ph, en
       logical :: converged_
 
-    end function thermal_stress
+    end function phase_thermal_constitutive
 
-    module function integrateDamageState(dt,co,ce) result(broken)
+    module function phase_damage_constitutive(Delta_t,co,ip,el) result(converged_)
-      real(pReal), intent(in) :: dt
+      real(pReal), intent(in) :: Delta_t
-      integer, intent(in) :: &
-        ce, &
-        co
-      logical :: broken
-    end function integrateDamageState
-
-    module function crystallite_stress(dt,co,ip,el) result(converged_)
-      real(pReal), intent(in) :: dt
       integer, intent(in) :: co, ip, el
       logical :: converged_
-    end function crystallite_stress
+    end function phase_damage_constitutive
 
-    !ToDo: Try to merge the all stiffness functions
+    module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged_)
+      real(pReal), intent(in) :: Delta_t
+      integer, intent(in) :: co, ip, el
+      logical :: converged_
+    end function phase_mechanical_constitutive
+
+    !ToDo: Merge all the stiffness functions
     module function phase_homogenizedC(ph,en) result(C)
       integer, intent(in) :: ph, en
       real(pReal), dimension(6,6) :: C
@@ -271,8 +273,8 @@ module phase
     end subroutine plastic_dependentState
 
 
-    module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me)
+    module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,en)
-      integer, intent(in) :: ph, me
+      integer, intent(in) :: ph, en
       real(pReal),   intent(in),  dimension(3,3) :: &
         S
       real(pReal),   intent(out), dimension(3,3) :: &
@@ -315,14 +317,14 @@ module phase
     plastic_nonlocal_updateCompatibility, &
     converged, &
     crystallite_init, &
-    crystallite_stress, &
+    phase_mechanical_constitutive, &
-    thermal_stress, &
+    phase_thermal_constitutive, &
+    phase_damage_constitutive, &
     phase_mechanical_dPdF, &
     crystallite_orientations, &
     crystallite_push33ToRef, &
     phase_restartWrite, &
     phase_restartRead, &
-    integrateDamageState, &
     phase_thermal_setField, &
     phase_set_phi, &
     phase_P, &
@@ -435,17 +437,9 @@ subroutine phase_restore(ce,includeL)
   logical, intent(in) :: includeL
   integer, intent(in) :: ce
 
-  integer :: &
-    co
-
-
-  do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
-    if (damageState(material_phaseID(co,ce))%sizeState > 0) &
-    damageState(material_phaseID(co,ce))%state( :,material_phaseEntry(co,ce)) = &
-      damageState(material_phaseID(co,ce))%state0(:,material_phaseEntry(co,ce))
-  enddo
 
   call mechanical_restore(ce,includeL)
+  call damage_restore(ce)
 
 end subroutine phase_restore
 
@@ -545,10 +539,6 @@ subroutine crystallite_init()
 
   phases => config_material%get('phase')
 
-  do ph = 1, phases%length
-    if (damageState(ph)%sizeState > 0) allocate(damageState(ph)%subState0,source=damageState(ph)%state0)  ! ToDo: hack
-  enddo
-
   print'(a42,1x,i10)', '    # of elements:                       ', eMax
   print'(a42,1x,i10)', '    # of integration points/element:     ', iMax
   print'(a42,1x,i10)', 'max # of constituents/integration point: ', cMax

src/phase_damage.f90

@@ -39,8 +39,8 @@ submodule(phase) damage
     end function isobrittle_init
 
 
-    module subroutine isobrittle_deltaState(C, Fe, ph, me)
+    module subroutine isobrittle_deltaState(C, Fe, ph, en)
-      integer, intent(in) :: ph,me
+      integer, intent(in) :: ph,en
       real(pReal),  intent(in), dimension(3,3) :: &
         Fe
       real(pReal),  intent(in), dimension(6,6) :: &
@@ -48,8 +48,8 @@ submodule(phase) damage
     end subroutine isobrittle_deltaState
 
 
-    module subroutine anisobrittle_dotState(S, ph, me)
+    module subroutine anisobrittle_dotState(S, ph, en)
-      integer, intent(in) :: ph,me
+      integer, intent(in) :: ph,en
       real(pReal),  intent(in), dimension(3,3) :: &
         S
     end subroutine anisobrittle_dotState
@@ -125,6 +125,29 @@ module subroutine damage_init
 end subroutine damage_init
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief calculate stress (P)
+!--------------------------------------------------------------------------------------------------
+module function phase_damage_constitutive(Delta_t,co,ip,el) result(converged_)
+
+  real(pReal), intent(in) :: Delta_t
+  integer, intent(in) :: &
+    co, &
+    ip, &
+    el
+  logical :: converged_
+
+  integer :: &
+    ph, en
+
+  ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
+  en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
+
+  converged_ = .not. integrateDamageState(Delta_t,ph,en)
+
+end function phase_damage_constitutive
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief returns the degraded/modified elasticity matrix
 !--------------------------------------------------------------------------------------------------
@@ -144,6 +167,26 @@ module function phase_damage_C(C_homogenized,ph,en) result(C)
 end function phase_damage_C
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief Restore data after homog cutback.
+!--------------------------------------------------------------------------------------------------
+module subroutine damage_restore(ce)
+
+  integer, intent(in) :: ce
+
+  integer :: &
+    co
+
+
+  do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
+    if (damageState(material_phaseID(co,ce))%sizeState > 0) &
+    damageState(material_phaseID(co,ce))%state( :,material_phaseEntry(co,ce)) = &
+      damageState(material_phaseID(co,ce))%state0(:,material_phaseEntry(co,ce))
+  enddo
+
+end subroutine damage_restore
+
+
 !----------------------------------------------------------------------------------------------
 !< @brief returns local part of nonlocal damage driving force
 !----------------------------------------------------------------------------------------------
@@ -173,23 +216,20 @@ module function phase_f_phi(phi,co,ce) result(f)
 end function phase_f_phi
 
 
-
 !--------------------------------------------------------------------------------------------------
 !> @brief integrate stress, state with adaptive 1st order explicit Euler method
 !> using Fixed Point Iteration to adapt the stepsize
 !--------------------------------------------------------------------------------------------------
-module function integrateDamageState(dt,co,ce) result(broken)
+function integrateDamageState(Delta_t,ph,en) result(broken)
 
-  real(pReal), intent(in) :: dt
+  real(pReal), intent(in) :: Delta_t
   integer, intent(in) :: &
-    ce, &
+    ph, &
-    co
+    en
   logical :: broken
 
   integer :: &
     NiterationState, &                                                                              !< number of iterations in state loop
-    ph, &
-    me, &
     size_so
   real(pReal) :: &
     zeta
@@ -199,8 +239,6 @@ module function integrateDamageState(dt,co,ce) result(broken)
   logical :: &
     converged_
 
-  ph = material_phaseID(co,ce)
-  me = material_phaseEntry(co,ce)
 
   if (damageState(ph)%sizeState == 0) then
     broken = .false.
@@ -208,37 +246,37 @@ module function integrateDamageState(dt,co,ce) result(broken)
   endif
 
   converged_ = .true.
-  broken = phase_damage_collectDotState(ph,me)
+  broken = phase_damage_collectDotState(ph,en)
   if(broken) return
 
-    size_so = damageState(ph)%sizeDotState
+  size_so = damageState(ph)%sizeDotState
-    damageState(ph)%state(1:size_so,me) = damageState(ph)%subState0(1:size_so,me) &
+  damageState(ph)%state(1:size_so,en) = damageState(ph)%state0  (1:size_so,en) &
-                                        + damageState(ph)%dotState (1:size_so,me) * dt
+                                      + damageState(ph)%dotState(1:size_so,en) * Delta_t
-    source_dotState(1:size_so,2) = 0.0_pReal
+  source_dotState(1:size_so,2) = 0.0_pReal
 
   iteration: do NiterationState = 1, num%nState
 
     if(nIterationState > 1) source_dotState(1:size_so,2) = source_dotState(1:size_so,1)
-    source_dotState(1:size_so,1) = damageState(ph)%dotState(:,me)
+    source_dotState(1:size_so,1) = damageState(ph)%dotState(:,en)
 
-    broken = phase_damage_collectDotState(ph,me)
+    broken = phase_damage_collectDotState(ph,en)
     if(broken) exit iteration
 
 
-      zeta = damper(damageState(ph)%dotState(:,me),source_dotState(1:size_so,1),source_dotState(1:size_so,2))
+      zeta = damper(damageState(ph)%dotState(:,en),source_dotState(1:size_so,1),source_dotState(1:size_so,2))
-      damageState(ph)%dotState(:,me) = damageState(ph)%dotState(:,me) * zeta &
+      damageState(ph)%dotState(:,en) = damageState(ph)%dotState(:,en) * zeta &
                                      + source_dotState(1:size_so,1)* (1.0_pReal - zeta)
-      r(1:size_so) = damageState(ph)%state    (1:size_so,me)  &
+      r(1:size_so) = damageState(ph)%state   (1:size_so,en)  &
-                   - damageState(ph)%subState0(1:size_so,me)  &
+                   - damageState(ph)%State0  (1:size_so,en)  &
-                   - damageState(ph)%dotState (1:size_so,me) * dt
+                   - damageState(ph)%dotState(1:size_so,en) * Delta_t
-      damageState(ph)%state(1:size_so,me) = damageState(ph)%state(1:size_so,me) - r(1:size_so)
+      damageState(ph)%state(1:size_so,en) = damageState(ph)%state(1:size_so,en) - r(1:size_so)
       converged_ = converged_  .and. converged(r(1:size_so), &
-                                               damageState(ph)%state(1:size_so,me), &
+                                               damageState(ph)%state(1:size_so,en), &
                                                damageState(ph)%atol(1:size_so))
 
 
     if(converged_) then
-      broken = phase_damage_deltaState(mechanical_F_e(ph,me),ph,me)
+      broken = phase_damage_deltaState(mechanical_F_e(ph,en),ph,en)
       exit iteration
     endif
 
@@ -300,11 +338,11 @@ end subroutine damage_results
 !--------------------------------------------------------------------------------------------------
 !> @brief contains the constitutive equation for calculating the rate of change of microstructure
 !--------------------------------------------------------------------------------------------------
-function phase_damage_collectDotState(ph,me) result(broken)
+function phase_damage_collectDotState(ph,en) result(broken)
 
   integer, intent(in) :: &
     ph, &
-    me                                                                                         !< counter in source loop
+    en                                                                                         !< counter in source loop
   logical :: broken
 
 
@@ -315,11 +353,11 @@ function phase_damage_collectDotState(ph,me) result(broken)
     sourceType: select case (phase_damage(ph))
 
       case (DAMAGE_ANISOBRITTLE_ID) sourceType
-        call anisobrittle_dotState(mechanical_S(ph,me), ph,me) ! correct stress?
+        call anisobrittle_dotState(mechanical_S(ph,en), ph,en) ! correct stress?
 
     end select sourceType
 
-    broken = broken .or. any(IEEE_is_NaN(damageState(ph)%dotState(:,me)))
+    broken = broken .or. any(IEEE_is_NaN(damageState(ph)%dotState(:,en)))
 
   endif
 
@@ -358,11 +396,11 @@ end function phase_K_phi
 !> @brief for constitutive models having an instantaneous change of state
 !> will return false if delta state is not needed/supported by the constitutive model
 !--------------------------------------------------------------------------------------------------
-function phase_damage_deltaState(Fe, ph, me) result(broken)
+function phase_damage_deltaState(Fe, ph, en) result(broken)
 
   integer, intent(in) :: &
     ph, &
-    me
+    en
   real(pReal),   intent(in), dimension(3,3) :: &
     Fe                                                                                              !< elastic deformation gradient
   integer :: &
@@ -379,13 +417,13 @@ function phase_damage_deltaState(Fe, ph, me) result(broken)
    sourceType: select case (phase_damage(ph))
 
     case (DAMAGE_ISOBRITTLE_ID) sourceType
-      call isobrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me)
+      call isobrittle_deltaState(phase_homogenizedC(ph,en), Fe, ph,en)
-      broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,me)))
+      broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,en)))
       if(.not. broken) then
         myOffset = damageState(ph)%offsetDeltaState
         mySize   = damageState(ph)%sizeDeltaState
-        damageState(ph)%state(myOffset + 1: myOffset + mySize,me) = &
+        damageState(ph)%state(myOffset + 1: myOffset + mySize,en) = &
-        damageState(ph)%state(myOffset + 1: myOffset + mySize,me) + damageState(ph)%deltaState(1:mySize,me)
+        damageState(ph)%state(myOffset + 1: myOffset + mySize,en) + damageState(ph)%deltaState(1:mySize,en)
       endif
 
   end select sourceType
@@ -436,13 +474,13 @@ module subroutine phase_set_phi(phi,co,ce)
 end subroutine phase_set_phi
 
 
-module function damage_phi(ph,me) result(phi)
+module function damage_phi(ph,en) result(phi)
 
-  integer, intent(in) :: ph, me
+  integer, intent(in) :: ph, en
   real(pReal) :: phi
 
 
-  phi = current(ph)%phi(me)
+  phi = current(ph)%phi(en)
 
 end function damage_phi
 

src/phase_damage_anisobrittle.f90

@@ -40,7 +40,7 @@ module function anisobrittle_init() result(mySources)
     phase, &
     sources, &
     src
-  integer :: Nmembers,p
+  integer :: Nmembers,ph
   integer, dimension(:), allocatable :: N_cl
   character(len=pStringLen) :: extmsg = ''
 
@@ -56,12 +56,12 @@ module function anisobrittle_init() result(mySources)
   allocate(param(phases%length))
 
 
-  do p = 1, phases%length
+  do ph = 1, phases%length
-    if(mySources(p)) then
+    if(mySources(ph)) then
-    phase => phases%get(p)
+      phase => phases%get(ph)
-    sources => phase%get('damage')
+      sources => phase%get('damage')
 
-        associate(prm  => param(p))
+      associate(prm  => param(ph))
         src => sources%get(1)
 
         N_cl = src%get_as1dInt('N_cl',defaultVal=emptyIntArray)
@@ -92,17 +92,15 @@ module function anisobrittle_init() result(mySources)
         if (any(prm%g_crit <  0.0_pReal)) extmsg = trim(extmsg)//' g_crit'
         if (any(prm%s_crit <  0.0_pReal)) extmsg = trim(extmsg)//' s_crit'
 
-        Nmembers = count(material_phaseID==p)
+        Nmembers = count(material_phaseID==ph)
-        call phase_allocateState(damageState(p),Nmembers,1,1,0)
+        call phase_allocateState(damageState(ph),Nmembers,1,1,0)
-        damageState(p)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal)
+        damageState(ph)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal)
-        if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi'
+        if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi'
 
-        end associate
+      end associate
 
-!--------------------------------------------------------------------------------------------------
+      if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_anisoBrittle)')
-!  exit if any parameter is out of range
+    endif
-        if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_anisoBrittle)')
-      endif
 
   enddo
 
@@ -112,10 +110,10 @@ end function anisobrittle_init
 !--------------------------------------------------------------------------------------------------
 !> @brief
 !--------------------------------------------------------------------------------------------------
-module subroutine anisobrittle_dotState(S, ph,me)
+module subroutine anisobrittle_dotState(S, ph,en)
 
   integer, intent(in) :: &
-    ph,me
+    ph,en
   real(pReal),  intent(in), dimension(3,3) :: &
     S
 
@@ -126,15 +124,15 @@ module subroutine anisobrittle_dotState(S, ph,me)
 
 
   associate(prm => param(ph))
-    damageState(ph)%dotState(1,me) = 0.0_pReal
+    damageState(ph)%dotState(1,en) = 0.0_pReal
     do i = 1, prm%sum_N_cl
       traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i))
       traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i))
       traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i))
 
-      traction_crit = prm%g_crit(i)*damage_phi(ph,me)**2.0_pReal
+      traction_crit = prm%g_crit(i)*damage_phi(ph,en)**2.0_pReal
 
-      damageState(ph)%dotState(1,me) = damageState(ph)%dotState(1,me) &
+      damageState(ph)%dotState(1,en) = damageState(ph)%dotState(1,en) &
           + prm%dot_o / prm%s_crit(i) &
             * ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%q + &
                (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%q + &
@@ -171,10 +169,10 @@ end subroutine anisobrittle_results
 !--------------------------------------------------------------------------------------------------
 !> @brief  contains the constitutive equation for calculating the velocity gradient
 !--------------------------------------------------------------------------------------------------
-module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me)
+module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,en)
 
   integer, intent(in) :: &
-    ph,me
+    ph,en
   real(pReal),   intent(in),  dimension(3,3) :: &
     S
   real(pReal),   intent(out), dimension(3,3) :: &
@@ -193,7 +191,7 @@ module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me)
   dLd_dTstar = 0.0_pReal
   associate(prm => param(ph))
   do i = 1,prm%sum_N_cl
-    traction_crit = prm%g_crit(i)*damage_phi(ph,me)**2.0_pReal
+    traction_crit = prm%g_crit(i)*damage_phi(ph,en)**2.0_pReal
 
     traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i))
     if (abs(traction_d) > traction_crit + tol_math_check) then

src/phase_damage_isobrittle.f90

@@ -13,7 +13,15 @@ submodule(phase:damage) isobrittle
       output
   end type tParameters
 
-  type(tParameters), dimension(:), allocatable :: param                                             !< containers of constitutive parameters (len Ninstances)
+  type :: tIsobrittleState
+    real(pReal), pointer, dimension(:) :: &                                                         !< vectors along Nmembers
+      r_W                                                                                           !< ratio between actual and critical strain energy density
+  end type tIsobrittleState
+
+  type(tParameters),      allocatable, dimension(:) :: param                                        !< containers of constitutive parameters (len Ninstances)
+  type(tIsobrittleState), allocatable, dimension(:) :: &
+    deltaState, &
+    state
 
 contains
 
@@ -44,13 +52,15 @@ module function isobrittle_init() result(mySources)
 
   phases => config_material%get('phase')
   allocate(param(phases%length))
+  allocate(state(phases%length))
+  allocate(deltaState(phases%length))
 
   do ph = 1, phases%length
     if(mySources(ph)) then
-    phase => phases%get(ph)
+      phase => phases%get(ph)
-    sources => phase%get('damage')
+      sources => phase%get('damage')
 
-        associate(prm  => param(ph))
+      associate(prm => param(ph), dlt => deltaState(ph), stt => state(ph))
         src => sources%get(1)
 
         prm%W_crit = src%get_asFloat('W_crit')
@@ -69,12 +79,14 @@ module function isobrittle_init() result(mySources)
         damageState(ph)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal)
         if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi'
 
-        end associate
+        stt%r_W => damageState(ph)%state(1,:)
+        dlt%r_W => damageState(ph)%deltaState(1,:)
 
-!--------------------------------------------------------------------------------------------------
+      end associate
-!  exit if any parameter is out of range
+
-        if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isobrittle)')
+
-      endif
+      if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isobrittle)')
+    endif
 
   enddo
 
@@ -85,29 +97,27 @@ end function isobrittle_init
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates derived quantities from state
 !--------------------------------------------------------------------------------------------------
-module subroutine isobrittle_deltaState(C, Fe, ph,me)
+module subroutine isobrittle_deltaState(C, Fe, ph,en)
 
-  integer, intent(in) :: ph,me
+  integer, intent(in) :: ph,en
   real(pReal),  intent(in), dimension(3,3) :: &
     Fe
   real(pReal),  intent(in), dimension(6,6) :: &
     C
 
   real(pReal), dimension(6) :: &
-    strain
+    epsilon
   real(pReal) :: &
-    strainenergy
+    r_W
 
 
-  strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3)
+  epsilon = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3)
 
-  associate(prm => param(ph))
+  associate(prm => param(ph), stt => state(ph), dlt => deltaState(ph))
-    strainenergy = 2.0_pReal*sum(strain*matmul(C,strain))/prm%W_crit
+
-    ! ToDo: check strainenergy = 2.0_pReal*dot_product(strain,matmul(C,strain))/prm%W_crit
+    r_W = 2.0_pReal*dot_product(epsilon,matmul(C,epsilon))/prm%W_crit
+    dlt%r_W(en) = merge(r_W - stt%r_W(en), 0.0_pReal, r_W > stt%r_W(en))
 
-    damageState(ph)%deltaState(1,me) = merge(strainenergy - damageState(ph)%state(1,me), &
-                                             damageState(ph)%subState0(1,me) - damageState(ph)%state(1,me), &
-                                             strainenergy > damageState(ph)%subState0(1,me))
   end associate
 
 end subroutine isobrittle_deltaState
@@ -124,14 +134,15 @@ module subroutine isobrittle_results(phase,group)
   integer :: o
 
 
-  associate(prm => param(phase), &
+  associate(prm => param(phase), stt => damageState(phase)%state) ! point to state and output r_W (is scalar, not 1D vector)
-            stt => damageState(phase)%state)
+
     outputsLoop: do o = 1,size(prm%output)
       select case(trim(prm%output(o)))
         case ('f_phi')
-          call results_writeDataset(stt,group,trim(prm%output(o)),'driving force','J/m³')
+          call results_writeDataset(stt,group,trim(prm%output(o)),'driving force','J/m³') ! Wrong, this is dimensionless
       end select
     enddo outputsLoop
+
   end associate
 
 end subroutine isobrittle_results

src/phase_mechanical.f90

@@ -44,9 +44,9 @@ submodule(phase) mechanical
 
   interface
 
-    module subroutine eigendeformation_init(phases)
+    module subroutine eigen_init(phases)
       class(tNode), pointer :: phases
-    end subroutine eigendeformation_init
+    end subroutine eigen_init
 
     module subroutine elastic_init(phases)
       class(tNode), pointer :: phases
@@ -197,8 +197,7 @@ module subroutine mechanical_init(materials,phases)
     Nmembers
   class(tNode), pointer :: &
     num_crystallite, &
-    material, &
+
-    constituents, &
     phase, &
     mech
 
@@ -303,7 +302,7 @@ module subroutine mechanical_init(materials,phases)
   end select
 
 
-  call eigendeformation_init(phases)
+  call eigen_init(phases)
 
 
 end subroutine mechanical_init
@@ -977,9 +976,9 @@ end subroutine mechanical_forward
 !--------------------------------------------------------------------------------------------------
 !> @brief calculate stress (P)
 !--------------------------------------------------------------------------------------------------
-module function crystallite_stress(dt,co,ip,el) result(converged_)
+module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged_)
 
-  real(pReal), intent(in) :: dt
+  real(pReal), intent(in) :: Delta_t
   integer, intent(in) :: &
     co, &
     ip, &
@@ -1009,17 +1008,13 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
   subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,en)
   subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,en)
   subState0 = plasticState(ph)%State0(:,en)
-
-  if (damageState(ph)%sizeState > 0) &
-    damageState(ph)%subState0(:,en) = damageState(ph)%state0(:,en)
-
   subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,en)
   subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,en)
   subF0  = phase_mechanical_F0(ph)%data(1:3,1:3,en)
   subFrac = 0.0_pReal
-  subStep = 1.0_pReal/num%subStepSizeCryst
   todo = .true.
-  converged_ = .false.                                                      ! pretend failed step of 1/subStepSizeCryst
+  subStep = 1.0_pReal/num%subStepSizeCryst
+  converged_ = .false.                                                                              ! pretend failed step of 1/subStepSizeCryst
 
   todo = .true.
   cutbackLooping: do while (todo)
@@ -1038,9 +1033,6 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
         subFp0 = phase_mechanical_Fp(ph)%data(1:3,1:3,en)
         subFi0 = phase_mechanical_Fi(ph)%data(1:3,1:3,en)
         subState0 = plasticState(ph)%state(:,en)
-        if (damageState(ph)%sizeState > 0) &
-          damageState(ph)%subState0(:,en) = damageState(ph)%state(:,en)
-
       endif
 !--------------------------------------------------------------------------------------------------
 !  cut back (reduced time and restore)
@@ -1048,16 +1040,13 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
       subStep       = num%subStepSizeCryst * subStep
       phase_mechanical_Fp(ph)%data(1:3,1:3,en) = subFp0
       phase_mechanical_Fi(ph)%data(1:3,1:3,en) = subFi0
-      phase_mechanical_S(ph)%data(1:3,1:3,en) = phase_mechanical_S0(ph)%data(1:3,1:3,en)          ! why no subS0 ? is S0 of any use?
+      phase_mechanical_S(ph)%data(1:3,1:3,en) = phase_mechanical_S0(ph)%data(1:3,1:3,en)            ! why no subS0 ? is S0 of any use?
       if (subStep < 1.0_pReal) then                                                                 ! actual (not initial) cutback
         phase_mechanical_Lp(ph)%data(1:3,1:3,en) = subLp0
         phase_mechanical_Li(ph)%data(1:3,1:3,en) = subLi0
       endif
       plasticState(ph)%state(:,en) = subState0
-      if (damageState(ph)%sizeState > 0) &
+      todo = subStep > num%subStepMinCryst                                                          ! still on track or already done (beyond repair)
-        damageState(ph)%state(:,en) = damageState(ph)%subState0(:,en)
-
-      todo = subStep > num%subStepMinCryst                          ! still on track or already done (beyond repair)
     endif
 
 !--------------------------------------------------------------------------------------------------
@@ -1065,13 +1054,12 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
     if (todo) then
       subF = subF0 &
            + subStep * (phase_mechanical_F(ph)%data(1:3,1:3,en) - phase_mechanical_F0(ph)%data(1:3,1:3,en))
-      converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el)
+      converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * Delta_t,co,ip,el)
-      converged_ = converged_ .and. .not. integrateDamageState(subStep * dt,co,(el-1)*discretization_nIPs + ip)
     endif
 
   enddo cutbackLooping
 
-end function crystallite_stress
+end function phase_mechanical_constitutive
 
 
 !--------------------------------------------------------------------------------------------------
@@ -1104,12 +1092,13 @@ module subroutine mechanical_restore(ce,includeL)
 
 end subroutine mechanical_restore
 
+
 !--------------------------------------------------------------------------------------------------
 !> @brief Calculate tangent (dPdF).
 !--------------------------------------------------------------------------------------------------
-module function phase_mechanical_dPdF(dt,co,ce) result(dPdF)
+module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
 
-  real(pReal), intent(in) :: dt
+  real(pReal), intent(in) :: Delta_t
   integer, intent(in) :: &
     co, &                                                                                            !< counter in constituent loop
     ce
@@ -1159,11 +1148,11 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF)
     lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal
     do o=1,3; do p=1,3
       lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) &
-                            + matmul(invSubFi0,dLidFi(1:3,1:3,o,p)) * dt
+                            + matmul(invSubFi0,dLidFi(1:3,1:3,o,p)) * Delta_t
       lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) &
                             + invFi*invFi(p,o)
       rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) &
-                            - matmul(invSubFi0,dLidS(1:3,1:3,o,p)) * dt
+                            - matmul(invSubFi0,dLidS(1:3,1:3,o,p)) * Delta_t
     enddo; enddo
     call math_invert(temp_99,error,math_3333to99(lhs_3333))
     if (error) then
@@ -1192,7 +1181,7 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF)
     temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), invFi) &
                            + matmul(temp_33_3,dLidS(1:3,1:3,p,o))
   enddo; enddo
-  lhs_3333 = math_mul3333xx3333(dSdFe,temp_3333) * dt &
+  lhs_3333 = math_mul3333xx3333(dSdFe,temp_3333) * Delta_t &
            + math_mul3333xx3333(dSdFi,dFidS)
 
   call math_invert(temp_99,error,math_eye(9)+math_3333to99(lhs_3333))
@@ -1208,7 +1197,7 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF)
 ! calculate dFpinvdF
   temp_3333 = math_mul3333xx3333(dLpdS,dSdF)
   do o=1,3; do p=1,3
-    dFpinvdF(1:3,1:3,p,o) = - matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi)) * dt
+    dFpinvdF(1:3,1:3,p,o) = - matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi)) * Delta_t
   enddo; enddo
 
 !--------------------------------------------------------------------------------------------------

src/phase_mechanical_eigen.f90

@@ -32,7 +32,7 @@ submodule(phase:mechanical) eigen
 contains
 
 
-module subroutine eigendeformation_init(phases)
+module subroutine eigen_init(phases)
 
   class(tNode), pointer :: &
     phases
@@ -68,7 +68,7 @@ module subroutine eigendeformation_init(phases)
   where(damage_anisobrittle_init())  model_damage = KINEMATICS_cleavage_opening_ID
 
 
-end subroutine eigendeformation_init
+end subroutine eigen_init
 
 
 !--------------------------------------------------------------------------------------------------

src/phase_thermal.f90

@@ -216,7 +216,7 @@ module function phase_K_T(co,ce) result(K)
 end function phase_K_T
 
 
-module function thermal_stress(Delta_t,ph,en) result(converged_)           ! ?? why is this called "stress" when it seems closer to "updateState" ??
+module function phase_thermal_constitutive(Delta_t,ph,en) result(converged_)
 
   real(pReal), intent(in) :: Delta_t
   integer, intent(in) :: ph, en
@@ -225,7 +225,7 @@ module function thermal_stress(Delta_t,ph,en) result(converged_)           ! ??
 
   converged_ = .not. integrateThermalState(Delta_t,ph,en)
 
-end function thermal_stress
+end function phase_thermal_constitutive
 
 
 !--------------------------------------------------------------------------------------------------

src/prec.f90

@@ -38,16 +38,14 @@ module prec
       sizeDotState     = 0, &                                                                       !< size of dot state, i.e. state(1:sizeDot) follows time evolution by dotState rates
       offsetDeltaState = 0, &                                                                       !< index offset of delta state
       sizeDeltaState   = 0                                                                          !< size of delta state, i.e. state(offset+1:offset+sizeDelta) follows time evolution by deltaState increments
-    ! http://stackoverflow.com/questions/3948210
+    real(pReal), allocatable, dimension(:) :: &
-    real(pReal), pointer,     dimension(:), contiguous :: &
       atol
-    real(pReal), pointer,     dimension(:,:), contiguous :: &                                       ! a pointer is needed here because we might point to state/doState. However, they will never point to something, but are rather allocated and, hence, contiguous
+    ! http://stackoverflow.com/questions/3948210
+    real(pReal), pointer,     dimension(:,:), contiguous :: &                                       !< is basically an allocatable+target, but in a type needs to be pointer
       state0, &
       state, &                                                                                      !< state
       dotState, &                                                                                   !< rate of state change
       deltaState                                                                                    !< increment of state change
-    real(pReal), allocatable, dimension(:,:) :: &
-      subState0
   end type
 
   type, extends(tState) :: tPlasticState
@@ -61,12 +59,9 @@ module prec
   real(pReal), private, parameter :: PREAL_EPSILON = epsilon(0.0_pReal)                             !< minimum positive number such that 1.0 + EPSILON /= 1.0.
   real(pReal), private, parameter :: PREAL_MIN     = tiny(0.0_pReal)                                !< smallest normalized floating point number
 
-  integer,                   dimension(0), parameter :: &
+  integer,                   dimension(0), parameter :: emptyIntArray    = [integer::]
-    emptyIntArray    = [integer::]
+  real(pReal),               dimension(0), parameter :: emptyRealArray   = [real(pReal)::]
-  real(pReal),               dimension(0), parameter :: &
+  character(len=pStringLen), dimension(0), parameter :: emptyStringArray = [character(len=pStringLen)::]
-    emptyRealArray   = [real(pReal)::]
-  character(len=pStringLen), dimension(0), parameter :: &
-    emptyStringArray = [character(len=pStringLen)::]
 
   private :: &
     selfTest