From d59051f57670770c2653bdf4f6b3c6e3418ea2eb Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Wed, 7 Apr 2021 22:41:49 +0200
Subject: [PATCH 01/24] systematic names

---
 src/grid/grid_thermal_spectral.f90 |  8 +++---
 src/homogenization.f90             | 24 ++++++++---------
 src/homogenization_thermal.f90     | 36 ++++++++++---------------
 src/phase.f90                      |  9 +++----
 src/phase_thermal.f90              | 42 +++++++++++++-----------------
 src/phase_thermal_dissipation.f90  | 10 +++----
 src/phase_thermal_externalheat.f90 | 12 ++++-----
 7 files changed, 61 insertions(+), 80 deletions(-)

diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90
index 14e2affb9..7d4878b12 100644
--- a/src/grid/grid_thermal_spectral.f90
+++ b/src/grid/grid_thermal_spectral.f90
@@ -258,7 +258,6 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
   PetscObject :: dummy
   PetscErrorCode :: ierr
   integer :: i, j, k, ce
-  real(pReal)   :: Tdot
 
   T_current = x_scal
 !--------------------------------------------------------------------------------------------------
@@ -271,7 +270,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
   ce = 0
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
     ce = ce + 1
-    vectorField_real(1:3,i,j,k) = matmul(thermal_conduction_getConductivity(ce) - K_ref, &
+    vectorField_real(1:3,i,j,k) = matmul(homogenization_K(ce) - K_ref, &
                                          vectorField_real(1:3,i,j,k))
   enddo; enddo; enddo
   call utilities_FFTvectorForward
@@ -280,8 +279,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
   ce = 0
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
     ce = ce + 1
-    call thermal_conduction_getSource(Tdot,1,ce)
-    scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + Tdot) &
+    scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + homogenization_f_T(ce)) &
                             + homogenization_thermal_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) &
                             + mu_ref*T_current(i,j,k)
   enddo; enddo; enddo
@@ -311,7 +309,7 @@ subroutine updateReference
   mu_ref = 0.0_pReal
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
     ce = ce + 1
-    K_ref  = K_ref  + thermal_conduction_getConductivity(ce)
+    K_ref  = K_ref  + homogenization_K(ce)
     mu_ref = mu_ref + homogenization_thermal_mu_T(ce)
   enddo; enddo; enddo
   K_ref = K_ref*wgt
diff --git a/src/homogenization.f90 b/src/homogenization.f90
index 4ce75feca..478efad53 100644
--- a/src/homogenization.f90
+++ b/src/homogenization.f90
@@ -136,10 +136,10 @@ module homogenization
     end function mechanical_updateState
 
 
-    module function thermal_conduction_getConductivity(ce) result(K)
+    module function homogenization_K(ce) result(K)
       integer, intent(in) :: ce
       real(pReal), dimension(3,3) :: K
-    end function thermal_conduction_getConductivity
+    end function homogenization_K
 
     module function homogenization_thermal_mu_T(ce) result(mu_T)
       integer, intent(in) :: ce
@@ -151,17 +151,15 @@ module homogenization
       real(pReal),   intent(in) :: T, dot_T
     end subroutine homogenization_thermal_setField
 
-    module function homogenization_thermal_T(ce) result(T)
+    module function homogenization_T(ce) result(T)
       integer, intent(in) :: ce
       real(pReal) :: T
-    end function homogenization_thermal_T
+    end function homogenization_T
 
-    module subroutine thermal_conduction_getSource(Tdot, ip, el)
-      integer, intent(in) :: &
-        ip, &
-        el
-      real(pReal), intent(out) :: Tdot
-    end subroutine thermal_conduction_getSource
+    module function homogenization_f_T(ce) result(f_T)
+      integer, intent(in) :: ce
+      real(pReal) :: f_T
+    end function homogenization_f_T
 
     module function damage_nonlocal_getMobility(ce) result(M)
       integer, intent(in) :: ce
@@ -188,13 +186,13 @@ module homogenization
     homogenization_init, &
     materialpoint_stressAndItsTangent, &
     homogenization_thermal_mu_T, &
-    thermal_conduction_getConductivity, &
-    thermal_conduction_getSource, &
+    homogenization_K, &
+    homogenization_f_T, &
     damage_nonlocal_getMobility, &
     damage_nonlocal_getSourceAndItsTangent, &
     homogenization_set_phi, &
     homogenization_thermal_setfield, &
-    homogenization_thermal_T, &
+    homogenization_T, &
     homogenization_forward, &
     homogenization_results, &
     homogenization_restartRead, &
diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90
index 8dcddda7a..a372a7494 100644
--- a/src/homogenization_thermal.f90
+++ b/src/homogenization_thermal.f90
@@ -109,7 +109,7 @@ end subroutine thermal_homogenize
 !--------------------------------------------------------------------------------------------------
 !> @brief return homogenized thermal conductivity in reference configuration
 !--------------------------------------------------------------------------------------------------
-module function thermal_conduction_getConductivity(ce) result(K)
+module function homogenization_K(ce) result(K)
 
   integer, intent(in) :: ce
   real(pReal), dimension(3,3) :: K
@@ -125,11 +125,11 @@ module function thermal_conduction_getConductivity(ce) result(K)
 
   K = K / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
 
-end function thermal_conduction_getConductivity
+end function homogenization_K
 
 
 module function homogenization_thermal_mu_T(ce) result(mu_T)
-  
+
   integer, intent(in) :: ce
   real(pReal) :: mu_T
 
@@ -220,43 +220,35 @@ module subroutine thermal_results(ho,group)
 end subroutine thermal_results
 
 
-module function homogenization_thermal_T(ce) result(T)
+module function homogenization_T(ce) result(T)
 
   integer, intent(in) :: ce
   real(pReal) :: T
 
   T = current(material_homogenizationID(ce))%T(material_homogenizationEntry(ce))
 
-end function homogenization_thermal_T
+end function homogenization_T
 
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief return heat generation rate
 !--------------------------------------------------------------------------------------------------
-module subroutine thermal_conduction_getSource(Tdot, ip, el)
+module function homogenization_f_T(ce) result(f_T)
 
-  integer, intent(in) :: &
-    ip, &
-    el
-  real(pReal), intent(out) :: &
-    Tdot
+  integer, intent(in) :: ce
+  real(pReal) :: f_T
 
-  integer :: co, ho,ph,me
-  real(pReal) :: dot_T_temp
+  integer :: co
 
-  ho = material_homogenizationAt(el)
-  Tdot = 0.0_pReal
-  do co = 1, homogenization_Nconstituents(ho)
-     ph = material_phaseAt(co,el)
-     me = material_phasememberAt(co,ip,el)
-     call phase_thermal_getRate(dot_T_temp, ph,me)
-     Tdot = Tdot + dot_T_temp
+  f_T = phase_f_T(material_phaseID(1,ce),material_phaseEntry(1,ce))
+  do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
+    f_T = f_T + phase_f_T(material_phaseID(co,ce),material_phaseEntry(co,ce))
   enddo
 
-  Tdot = Tdot/real(homogenization_Nconstituents(ho),pReal)
+  f_T = f_T/real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
 
-end subroutine thermal_conduction_getSource
+end function homogenization_f_T
 
 
 end submodule thermal
diff --git a/src/phase.f90 b/src/phase.f90
index 1adba03d9..dce63a5ac 100644
--- a/src/phase.f90
+++ b/src/phase.f90
@@ -235,11 +235,10 @@ module phase
         phi_dot
     end function phase_damage_phi_dot
 
-    module subroutine phase_thermal_getRate(TDot, ph,me)
+    module function phase_f_T(ph,me) result(f_T)
       integer, intent(in) :: ph, me
-      real(pReal), intent(out) :: &
-        TDot
-    end subroutine phase_thermal_getRate
+      real(pReal) :: f_T
+    end function phase_f_T
 
     module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
       integer, intent(in) :: &
@@ -301,7 +300,7 @@ module phase
     phase_init, &
     phase_homogenizedC, &
     phase_damage_phi_dot, &
-    phase_thermal_getRate, &
+    phase_f_T, &
     phase_results, &
     phase_allocateState, &
     phase_forward, &
diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90
index 1e157f338..e7f9abdc5 100644
--- a/src/phase_thermal.f90
+++ b/src/phase_thermal.f90
@@ -21,7 +21,7 @@ submodule(phase) thermal
   integer(kind(THERMAL_UNDEFINED_ID)),  dimension(:,:), allocatable :: &
     thermal_source
 
-  type(tDataContainer), dimension(:), allocatable :: current          ! ?? not very telling name. Better: "field" ??
+  type(tDataContainer), dimension(:), allocatable :: current          ! ?? not very telling name. Better: "field" ?? MD: current(ho)%T(me) reads quite good
 
   integer :: thermal_source_maxSizeDotState
 
@@ -45,21 +45,19 @@ submodule(phase) thermal
         me
     end subroutine externalheat_dotState
 
-    module subroutine dissipation_getRate(TDot, ph,me)
+    module function dissipation_f_T(ph,me) result(f_T)
       integer, intent(in) :: &
         ph, &
         me
-      real(pReal),  intent(out) :: &
-        TDot
-    end subroutine dissipation_getRate
+      real(pReal) :: f_T
+    end function dissipation_f_T
 
-    module subroutine externalheat_getRate(TDot, ph,me)
+    module function externalheat_f_T(ph,me)  result(f_T)
       integer, intent(in) :: &
         ph, &
         me
-      real(pReal),  intent(out) :: &
-        TDot
-    end subroutine externalheat_getRate
+      real(pReal) :: f_T
+    end function externalheat_f_T
 
  end interface
 
@@ -123,35 +121,31 @@ end subroutine thermal_init
 !----------------------------------------------------------------------------------------------
 !< @brief calculates thermal dissipation rate
 !----------------------------------------------------------------------------------------------
-module subroutine phase_thermal_getRate(TDot, ph,me)
+module function phase_f_T(ph,me) result(f_T)
 
   integer, intent(in) :: ph, me
-  real(pReal), intent(out) :: &
-    TDot
-
-  real(pReal) :: &
-    my_Tdot
-  integer :: &
-    so
+  real(pReal) :: f_T
 
 
-  TDot = 0.0_pReal
+  integer :: so
+
+
+  f_T = 0.0_pReal
 
   do so = 1, thermal_Nsources(ph)
    select case(thermal_source(so,ph))
+
      case (THERMAL_DISSIPATION_ID)
-       call dissipation_getRate(my_Tdot, ph,me)
+       f_T = f_T + dissipation_f_T(ph,me)
 
      case (THERMAL_EXTERNALHEAT_ID)
-       call externalheat_getRate(my_Tdot, ph,me)
+       f_T = f_T + externalheat_f_T(ph,me)
 
-     case default
-       my_Tdot = 0.0_pReal
    end select
-   Tdot = Tdot + my_Tdot
+
   enddo
 
-end subroutine phase_thermal_getRate
+end function phase_f_T
 
 
 !--------------------------------------------------------------------------------------------------
diff --git a/src/phase_thermal_dissipation.f90 b/src/phase_thermal_dissipation.f90
index b16ed3cf7..3a4ee651a 100644
--- a/src/phase_thermal_dissipation.f90
+++ b/src/phase_thermal_dissipation.f90
@@ -69,17 +69,17 @@ end function dissipation_init
 !--------------------------------------------------------------------------------------------------
 !> @brief Ninstancess dissipation rate
 !--------------------------------------------------------------------------------------------------
-module subroutine dissipation_getRate(TDot, ph,me)
+module function dissipation_f_T(ph,me) result(f_T)
 
   integer, intent(in) :: ph, me
-  real(pReal),  intent(out) :: &
-    TDot
+  real(pReal) :: &
+    f_T
 
 
   associate(prm => param(ph))
-    TDot = prm%kappa*sum(abs(mechanical_S(ph,me)*mechanical_L_p(ph,me)))
+    f_T = prm%kappa*sum(abs(mechanical_S(ph,me)*mechanical_L_p(ph,me)))
   end associate
 
-end subroutine dissipation_getRate
+end function dissipation_f_T
 
 end submodule dissipation
diff --git a/src/phase_thermal_externalheat.f90 b/src/phase_thermal_externalheat.f90
index d5bbc7c38..6d4403ab8 100644
--- a/src/phase_thermal_externalheat.f90
+++ b/src/phase_thermal_externalheat.f90
@@ -100,13 +100,13 @@ end subroutine externalheat_dotState
 !--------------------------------------------------------------------------------------------------
 !> @brief returns local heat generation rate
 !--------------------------------------------------------------------------------------------------
-module subroutine externalheat_getRate(TDot, ph, me)
+module function externalheat_f_T(ph,me) result(f_T)
 
   integer, intent(in) :: &
     ph, &
     me
-  real(pReal),  intent(out) :: &
-    TDot
+  real(pReal) :: &
+    f_T
 
   integer :: &
     so, interval
@@ -122,12 +122,12 @@ module subroutine externalheat_getRate(TDot, ph, me)
       if (     (frac_time <  0.0_pReal .and. interval == 1) &
           .or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) &
           .or. (frac_time >= 0.0_pReal .and. frac_time < 1.0_pReal) ) &
-        TDot = prm%f_T(interval  ) * (1.0_pReal - frac_time) + &
-               prm%f_T(interval+1) * frac_time                                                      ! interpolate heat rate between segment boundaries...
+        f_T = prm%f_T(interval  ) * (1.0_pReal - frac_time) + &
+              prm%f_T(interval+1) * frac_time                                                      ! interpolate heat rate between segment boundaries...
                                                                                                     ! ...or extrapolate if outside me bounds
     enddo
   end associate
 
-end subroutine externalheat_getRate
+end function externalheat_f_T
 
 end submodule externalheat

From c4765d3742076187d4742d387b18b3f7d64b037b Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Thu, 8 Apr 2021 13:31:21 +0200
Subject: [PATCH 02/24] following paper

---
 src/grid/grid_damage_spectral.f90 | 13 ++++++------
 src/homogenization.f90            | 34 +++++++++++++++----------------
 src/homogenization_damage.f90     | 15 +++++++-------
 src/phase.f90                     | 29 +++++++++++---------------
 src/phase_damage.f90              | 18 ++++------------
 5 files changed, 45 insertions(+), 64 deletions(-)

diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90
index 002a15708..7e0d4112a 100644
--- a/src/grid/grid_damage_spectral.f90
+++ b/src/grid/grid_damage_spectral.f90
@@ -259,7 +259,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
   PetscObject :: dummy
   PetscErrorCode :: ierr
   integer :: i, j, k, ce
-  real(pReal)   :: phiDot, mobility
+  real(pReal) :: mobility
 
   phi_current = x_scal
 !--------------------------------------------------------------------------------------------------
@@ -272,7 +272,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
   ce = 0
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
     ce = ce + 1
-    vectorField_real(1:3,i,j,k) = matmul(damage_nonlocal_getDiffusion(ce) - K_ref, &
+    vectorField_real(1:3,i,j,k) = matmul(homogenization_K_phi(ce) - K_ref, &
                                          vectorField_real(1:3,i,j,k))
   enddo; enddo; enddo
   call utilities_FFTvectorForward
@@ -281,9 +281,8 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
   ce = 0
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
     ce = ce + 1
-    call damage_nonlocal_getSourceAndItsTangent(phiDot, phi_current(i,j,k),ce)
-    mobility = damage_nonlocal_getMobility(ce)
-    scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + phiDot) &
+    mobility = homogenization_mu_phi(ce)
+    scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + homogenization_f_phi(phi_current(i,j,k),ce)) &
                             + mobility*(phi_lastInc(i,j,k) - phi_current(i,j,k)) &
                             + mu_ref*phi_current(i,j,k)
   enddo; enddo; enddo
@@ -317,8 +316,8 @@ subroutine updateReference
   mu_ref = 0.0_pReal
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
     ce = ce + 1
-    K_ref  = K_ref  + damage_nonlocal_getDiffusion(ce)
-    mu_ref = mu_ref + damage_nonlocal_getMobility(ce)
+    K_ref  = K_ref  + homogenization_K_phi(ce)
+    mu_ref = mu_ref + homogenization_mu_phi(ce)
   enddo; enddo; enddo
   K_ref = K_ref*wgt
   call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
diff --git a/src/homogenization.f90 b/src/homogenization.f90
index 478efad53..9fe78dd77 100644
--- a/src/homogenization.f90
+++ b/src/homogenization.f90
@@ -161,18 +161,16 @@ module homogenization
       real(pReal) :: f_T
     end function homogenization_f_T
 
-    module function damage_nonlocal_getMobility(ce) result(M)
+    module function homogenization_mu_phi(ce) result(M)
       integer, intent(in) :: ce
       real(pReal) :: M
-    end function damage_nonlocal_getMobility
+    end function homogenization_mu_phi
 
-    module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, phi, ce)
+    module function homogenization_f_phi(phi,ce) result(f_phi)
       integer, intent(in) :: ce
-      real(pReal), intent(in) :: &
-        phi
-      real(pReal), intent(out) :: &
-        phiDot
-    end subroutine damage_nonlocal_getSourceAndItsTangent
+      real(pReal), intent(in) :: phi
+      real(pReal) :: f_phi
+    end function homogenization_f_phi
 
     module subroutine homogenization_set_phi(phi,ce)
       integer, intent(in) :: ce
@@ -188,8 +186,8 @@ module homogenization
     homogenization_thermal_mu_T, &
     homogenization_K, &
     homogenization_f_T, &
-    damage_nonlocal_getMobility, &
-    damage_nonlocal_getSourceAndItsTangent, &
+    homogenization_mu_phi, &
+    homogenization_f_phi, &
     homogenization_set_phi, &
     homogenization_thermal_setfield, &
     homogenization_T, &
@@ -201,7 +199,7 @@ module homogenization
     DAMAGE_NONLOCAL_ID
 
   public :: &
-    damage_nonlocal_getDiffusion
+    homogenization_K_phi
 
 contains
 
@@ -448,27 +446,27 @@ end subroutine homogenization_restartRead
 !--------------------------------------------------------------------------------------------------
 !> @brief returns homogenized non local damage diffusion tensor in reference configuration
 !--------------------------------------------------------------------------------------------------
-function damage_nonlocal_getDiffusion(ce)
+function homogenization_K_phi(ce)
 
   integer, intent(in) :: ce
   real(pReal), dimension(3,3) :: &
-    damage_nonlocal_getDiffusion
+    homogenization_K_phi
   integer :: &
     ho, &
     co
 
   ho  = material_homogenizationID(ce)
-  damage_nonlocal_getDiffusion = 0.0_pReal
+  homogenization_K_phi = 0.0_pReal
 
   do co = 1, homogenization_Nconstituents(ho)
-    damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + &
+    homogenization_K_phi = homogenization_K_phi + &
       crystallite_push33ToRef(co,ce,lattice_D(1:3,1:3,material_phaseID(co,ce)))
   enddo
 
-  damage_nonlocal_getDiffusion = &
-    num_damage%charLength**2*damage_nonlocal_getDiffusion/real(homogenization_Nconstituents(ho),pReal)
+  homogenization_K_phi = &
+    num_damage%charLength**2*homogenization_K_phi/real(homogenization_Nconstituents(ho),pReal)
 
-end function damage_nonlocal_getDiffusion
+end function homogenization_K_phi
 
 
 !--------------------------------------------------------------------------------------------------
diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90
index a067dde43..19a5b0f83 100644
--- a/src/homogenization_damage.f90
+++ b/src/homogenization_damage.f90
@@ -100,7 +100,7 @@ module subroutine damage_partition(ce)
 
   if(damageState_h(material_homogenizationID(ce))%sizeState < 1) return
   phi = damagestate_h(material_homogenizationID(ce))%state(1,material_homogenizationEntry(ce))
-  call phase_damage_set_phi(phi,1,ce)
+  call phase_set_phi(phi,1,ce)
 
 end subroutine damage_partition
 
@@ -109,30 +109,29 @@ end subroutine damage_partition
 !--------------------------------------------------------------------------------------------------
 !> @brief Returns homogenized nonlocal damage mobility
 !--------------------------------------------------------------------------------------------------
-module function damage_nonlocal_getMobility(ce) result(M)
+module function homogenization_mu_phi(ce) result(M)
 
   integer, intent(in) :: ce
   real(pReal) :: M
 
   M = lattice_M(material_phaseID(1,ce))
 
-end function damage_nonlocal_getMobility
+end function homogenization_mu_phi
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief  calculates homogenized damage driving forces
 !--------------------------------------------------------------------------------------------------
-module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, phi, ce)
+module function homogenization_f_phi(phi,ce) result(f_phi)
 
   integer, intent(in) :: ce
   real(pReal), intent(in) :: &
     phi
-  real(pReal), intent(out) :: &
-    phiDot
+  real(pReal) :: f_phi
 
-  phiDot = phase_damage_phi_dot(phi, 1, ce)
+  f_phi = phase_f_phi(phi, 1, ce)
 
-end subroutine damage_nonlocal_getSourceAndItsTangent
+end function homogenization_f_phi
 
 
 !--------------------------------------------------------------------------------------------------
diff --git a/src/phase.f90 b/src/phase.f90
index dce63a5ac..57c83c4b4 100644
--- a/src/phase.f90
+++ b/src/phase.f90
@@ -145,26 +145,22 @@ module phase
       real(pReal), dimension(3,3) :: L_p
     end function mechanical_L_p
 
-    module function phase_F(co,ce) result(F)
-      integer, intent(in) :: co, ce
-      real(pReal), dimension(3,3) :: F
-    end function phase_F
-
     module function mechanical_F_e(ph,me) result(F_e)
       integer, intent(in) :: ph,me
       real(pReal), dimension(3,3) :: F_e
     end function mechanical_F_e
 
+
+    module function phase_F(co,ce) result(F)
+      integer, intent(in) :: co, ce
+      real(pReal), dimension(3,3) :: F
+    end function phase_F
+
     module function phase_P(co,ce) result(P)
       integer, intent(in) :: co, ce
       real(pReal), dimension(3,3) :: P
     end function phase_P
 
-    module function phase_damage_get_phi(co,ip,el) result(phi)
-      integer, intent(in) :: co, ip, el
-      real(pReal) :: phi
-    end function phase_damage_get_phi
-
     module function thermal_T(ph,me) result(T)
       integer, intent(in) :: ph,me
       real(pReal) :: T
@@ -191,10 +187,10 @@ module phase
       integer, intent(in) :: ce, co
     end subroutine phase_thermal_setField
 
-    module subroutine phase_damage_set_phi(phi,co,ce)
+    module subroutine phase_set_phi(phi,co,ce)
       real(pReal), intent(in) :: phi
       integer, intent(in) :: co, ce
-    end subroutine phase_damage_set_phi
+    end subroutine phase_set_phi
 
 ! == cleaned:end ===================================================================================
 
@@ -227,13 +223,13 @@ module phase
     end function phase_homogenizedC
 
 
-    module function phase_damage_phi_dot(phi,co,ce) result(phi_dot)
+    module function phase_f_phi(phi,co,ce) result(phi_dot)
       integer, intent(in) :: ce,co
       real(pReal), intent(in) :: &
         phi                                                                                         !< damage parameter
       real(pReal) :: &
         phi_dot
-    end function phase_damage_phi_dot
+    end function phase_f_phi
 
     module function phase_f_T(ph,me) result(f_T)
       integer, intent(in) :: ph, me
@@ -299,7 +295,7 @@ module phase
   public :: &
     phase_init, &
     phase_homogenizedC, &
-    phase_damage_phi_dot, &
+    phase_f_phi, &
     phase_f_T, &
     phase_results, &
     phase_allocateState, &
@@ -317,8 +313,7 @@ module phase
     phase_restartRead, &
     integrateDamageState, &
     phase_thermal_setField, &
-    phase_damage_set_phi, &
-    phase_damage_get_phi, &
+    phase_set_phi, &
     phase_P, &
     phase_set_F, &
     phase_F
diff --git a/src/phase_damage.f90 b/src/phase_damage.f90
index 62ce70368..a6d97032f 100644
--- a/src/phase_damage.f90
+++ b/src/phase_damage.f90
@@ -142,7 +142,7 @@ end subroutine damage_init
 !----------------------------------------------------------------------------------------------
 !< @brief returns local part of nonlocal damage driving force
 !----------------------------------------------------------------------------------------------
-module function phase_damage_phi_dot(phi,co,ce) result(phi_dot)
+module function phase_f_phi(phi,co,ce) result(phi_dot)
 
   integer, intent(in) :: ce,co
   real(pReal), intent(in) :: &
@@ -165,7 +165,7 @@ module function phase_damage_phi_dot(phi,co,ce) result(phi_dot)
       phi_dot = 0.0_pReal
   end select
 
-end function phase_damage_phi_dot
+end function phase_f_phi
 
 
 
@@ -411,7 +411,7 @@ end function source_active
 !----------------------------------------------------------------------------------------------
 !< @brief Set damage parameter
 !----------------------------------------------------------------------------------------------
-module subroutine phase_damage_set_phi(phi,co,ce)
+module subroutine phase_set_phi(phi,co,ce)
 
   real(pReal), intent(in) :: phi
   integer, intent(in) :: ce, co
@@ -419,17 +419,7 @@ module subroutine phase_damage_set_phi(phi,co,ce)
 
   current(material_phaseID(co,ce))%phi(material_phaseEntry(co,ce)) = phi
 
-end subroutine phase_damage_set_phi
-
-
-module function phase_damage_get_phi(co,ip,el) result(phi)
-
-  integer, intent(in) :: co, ip, el
-  real(pReal) :: phi
-
-  phi = current(material_phaseAt(co,el))%phi(material_phaseMemberAt(co,ip,el))
-
-end function phase_damage_get_phi
+end subroutine phase_set_phi
 
 
 module function damage_phi(ph,me) result(phi)

From 1b89032086e1b41090b2c5a45fcaf81714810637 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Thu, 8 Apr 2021 23:40:20 +0200
Subject: [PATCH 03/24] names as in DAMASK paper

---
 src/grid/grid_thermal_spectral.f90 |  8 ++++----
 src/homogenization.f90             | 30 ++++++++++++++----------------
 src/homogenization_damage.f90      | 12 ++++++------
 src/homogenization_thermal.f90     | 22 +++++++++++-----------
 src/phase.f90                      |  8 ++++----
 src/phase_damage.f90               | 10 +++++-----
 src/phase_thermal.f90              | 10 +++++-----
 7 files changed, 49 insertions(+), 51 deletions(-)

diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90
index 7d4878b12..ea6caec0e 100644
--- a/src/grid/grid_thermal_spectral.f90
+++ b/src/grid/grid_thermal_spectral.f90
@@ -270,7 +270,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
   ce = 0
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
     ce = ce + 1
-    vectorField_real(1:3,i,j,k) = matmul(homogenization_K(ce) - K_ref, &
+    vectorField_real(1:3,i,j,k) = matmul(homogenization_K_T(ce) - K_ref, &
                                          vectorField_real(1:3,i,j,k))
   enddo; enddo; enddo
   call utilities_FFTvectorForward
@@ -280,7 +280,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
     ce = ce + 1
     scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + homogenization_f_T(ce)) &
-                            + homogenization_thermal_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) &
+                            + homogenization_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) &
                             + mu_ref*T_current(i,j,k)
   enddo; enddo; enddo
 
@@ -309,8 +309,8 @@ subroutine updateReference
   mu_ref = 0.0_pReal
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
     ce = ce + 1
-    K_ref  = K_ref  + homogenization_K(ce)
-    mu_ref = mu_ref + homogenization_thermal_mu_T(ce)
+    K_ref  = K_ref  + homogenization_K_T(ce)
+    mu_ref = mu_ref + homogenization_mu_T(ce)
   enddo; enddo; enddo
   K_ref = K_ref*wgt
   call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
diff --git a/src/homogenization.f90 b/src/homogenization.f90
index 9fe78dd77..34e71a0f2 100644
--- a/src/homogenization.f90
+++ b/src/homogenization.f90
@@ -136,15 +136,15 @@ module homogenization
     end function mechanical_updateState
 
 
-    module function homogenization_K(ce) result(K)
+    module function homogenization_K_T(ce) result(K)
       integer, intent(in) :: ce
       real(pReal), dimension(3,3) :: K
-    end function homogenization_K
+    end function homogenization_K_T
 
-    module function homogenization_thermal_mu_T(ce) result(mu_T)
+    module function homogenization_mu_T(ce) result(mu)
       integer, intent(in) :: ce
-      real(pReal) :: mu_T
-    end function homogenization_thermal_mu_T
+      real(pReal) :: mu
+    end function homogenization_mu_T
 
     module subroutine homogenization_thermal_setField(T,dot_T, ce)
       integer, intent(in) :: ce
@@ -156,20 +156,20 @@ module homogenization
       real(pReal) :: T
     end function homogenization_T
 
-    module function homogenization_f_T(ce) result(f_T)
+    module function homogenization_f_T(ce) result(f)
       integer, intent(in) :: ce
-      real(pReal) :: f_T
+      real(pReal) :: f
     end function homogenization_f_T
 
-    module function homogenization_mu_phi(ce) result(M)
+    module function homogenization_mu_phi(ce) result(mu)
       integer, intent(in) :: ce
-      real(pReal) :: M
+      real(pReal) :: mu
     end function homogenization_mu_phi
 
-    module function homogenization_f_phi(phi,ce) result(f_phi)
+    module function homogenization_f_phi(phi,ce) result(f)
       integer, intent(in) :: ce
       real(pReal), intent(in) :: phi
-      real(pReal) :: f_phi
+      real(pReal) :: f
     end function homogenization_f_phi
 
     module subroutine homogenization_set_phi(phi,ce)
@@ -183,9 +183,10 @@ module homogenization
   public ::  &
     homogenization_init, &
     materialpoint_stressAndItsTangent, &
-    homogenization_thermal_mu_T, &
-    homogenization_K, &
+    homogenization_mu_T, &
+    homogenization_K_T, &
     homogenization_f_T, &
+    homogenization_K_phi, &
     homogenization_mu_phi, &
     homogenization_f_phi, &
     homogenization_set_phi, &
@@ -198,9 +199,6 @@ module homogenization
     THERMAL_CONDUCTION_ID, &
     DAMAGE_NONLOCAL_ID
 
-  public :: &
-    homogenization_K_phi
-
 contains
 
 
diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90
index 19a5b0f83..f0f13bdbb 100644
--- a/src/homogenization_damage.f90
+++ b/src/homogenization_damage.f90
@@ -109,12 +109,12 @@ end subroutine damage_partition
 !--------------------------------------------------------------------------------------------------
 !> @brief Returns homogenized nonlocal damage mobility
 !--------------------------------------------------------------------------------------------------
-module function homogenization_mu_phi(ce) result(M)
+module function homogenization_mu_phi(ce) result(mu)
 
   integer, intent(in) :: ce
-  real(pReal) :: M
+  real(pReal) :: mu
 
-  M = lattice_M(material_phaseID(1,ce))
+  mu = lattice_M(material_phaseID(1,ce))
 
 end function homogenization_mu_phi
 
@@ -122,14 +122,14 @@ end function homogenization_mu_phi
 !--------------------------------------------------------------------------------------------------
 !> @brief  calculates homogenized damage driving forces
 !--------------------------------------------------------------------------------------------------
-module function homogenization_f_phi(phi,ce) result(f_phi)
+module function homogenization_f_phi(phi,ce) result(f)
 
   integer, intent(in) :: ce
   real(pReal), intent(in) :: &
     phi
-  real(pReal) :: f_phi
+  real(pReal) :: f
 
-  f_phi = phase_f_phi(phi, 1, ce)
+  f = phase_f_phi(phi, 1, ce)
 
 end function homogenization_f_phi
 
diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90
index a372a7494..8b4eafc3a 100644
--- a/src/homogenization_thermal.f90
+++ b/src/homogenization_thermal.f90
@@ -109,7 +109,7 @@ end subroutine thermal_homogenize
 !--------------------------------------------------------------------------------------------------
 !> @brief return homogenized thermal conductivity in reference configuration
 !--------------------------------------------------------------------------------------------------
-module function homogenization_K(ce) result(K)
+module function homogenization_K_T(ce) result(K)
 
   integer, intent(in) :: ce
   real(pReal), dimension(3,3) :: K
@@ -125,17 +125,17 @@ module function homogenization_K(ce) result(K)
 
   K = K / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
 
-end function homogenization_K
+end function homogenization_K_T
 
 
-module function homogenization_thermal_mu_T(ce) result(mu_T)
+module function homogenization_mu_T(ce) result(mu)
 
   integer, intent(in) :: ce
-  real(pReal) :: mu_T
+  real(pReal) :: mu
 
-  mu_T = c_P(ce) * rho(ce)
+  mu = c_P(ce) * rho(ce)
 
-end function homogenization_thermal_mu_T
+end function homogenization_mu_T
 
 
 !--------------------------------------------------------------------------------------------------
@@ -234,19 +234,19 @@ end function homogenization_T
 !--------------------------------------------------------------------------------------------------
 !> @brief return heat generation rate
 !--------------------------------------------------------------------------------------------------
-module function homogenization_f_T(ce) result(f_T)
+module function homogenization_f_T(ce) result(f)
 
   integer, intent(in) :: ce
-  real(pReal) :: f_T
+  real(pReal) :: f
 
   integer :: co
 
-  f_T = phase_f_T(material_phaseID(1,ce),material_phaseEntry(1,ce))
+  f = phase_f_T(material_phaseID(1,ce),material_phaseEntry(1,ce))
   do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
-    f_T = f_T + phase_f_T(material_phaseID(co,ce),material_phaseEntry(co,ce))
+    f = f + phase_f_T(material_phaseID(co,ce),material_phaseEntry(co,ce))
   enddo
 
-  f_T = f_T/real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
+  f = f/real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
 
 end function homogenization_f_T
 
diff --git a/src/phase.f90 b/src/phase.f90
index 57c83c4b4..346772ce4 100644
--- a/src/phase.f90
+++ b/src/phase.f90
@@ -223,17 +223,17 @@ module phase
     end function phase_homogenizedC
 
 
-    module function phase_f_phi(phi,co,ce) result(phi_dot)
+    module function phase_f_phi(phi,co,ce) result(f)
       integer, intent(in) :: ce,co
       real(pReal), intent(in) :: &
         phi                                                                                         !< damage parameter
       real(pReal) :: &
-        phi_dot
+        f
     end function phase_f_phi
 
-    module function phase_f_T(ph,me) result(f_T)
+    module function phase_f_T(ph,me) result(f)
       integer, intent(in) :: ph, me
-      real(pReal) :: f_T
+      real(pReal) :: f
     end function phase_f_T
 
     module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
diff --git a/src/phase_damage.f90 b/src/phase_damage.f90
index a6d97032f..189a5eabf 100644
--- a/src/phase_damage.f90
+++ b/src/phase_damage.f90
@@ -142,13 +142,13 @@ end subroutine damage_init
 !----------------------------------------------------------------------------------------------
 !< @brief returns local part of nonlocal damage driving force
 !----------------------------------------------------------------------------------------------
-module function phase_f_phi(phi,co,ce) result(phi_dot)
+module function phase_f_phi(phi,co,ce) result(f)
 
   integer, intent(in) :: ce,co
   real(pReal), intent(in) :: &
     phi                                                                                             !< damage parameter
   real(pReal) :: &
-    phi_dot
+    f
 
   integer :: &
     ph, &
@@ -159,10 +159,10 @@ module function phase_f_phi(phi,co,ce) result(phi_dot)
 
   select case(phase_source(ph))
     case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ISODUCTILE_ID,DAMAGE_ANISOBRITTLE_ID,DAMAGE_ANISODUCTILE_ID)
-      phi_dot = 1.0_pReal &
-              - phi*damageState(ph)%state(1,en)
+      f = 1.0_pReal &
+        - phi*damageState(ph)%state(1,en)
     case default
-      phi_dot = 0.0_pReal
+      f = 0.0_pReal
   end select
 
 end function phase_f_phi
diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90
index e7f9abdc5..807e1f655 100644
--- a/src/phase_thermal.f90
+++ b/src/phase_thermal.f90
@@ -121,25 +121,25 @@ end subroutine thermal_init
 !----------------------------------------------------------------------------------------------
 !< @brief calculates thermal dissipation rate
 !----------------------------------------------------------------------------------------------
-module function phase_f_T(ph,me) result(f_T)
+module function phase_f_T(ph,me) result(f)
 
   integer, intent(in) :: ph, me
-  real(pReal) :: f_T
+  real(pReal) :: f
 
 
   integer :: so
 
 
-  f_T = 0.0_pReal
+  f = 0.0_pReal
 
   do so = 1, thermal_Nsources(ph)
    select case(thermal_source(so,ph))
 
      case (THERMAL_DISSIPATION_ID)
-       f_T = f_T + dissipation_f_T(ph,me)
+       f = f + dissipation_f_T(ph,me)
 
      case (THERMAL_EXTERNALHEAT_ID)
-       f_T = f_T + externalheat_f_T(ph,me)
+       f = f + externalheat_f_T(ph,me)
 
    end select
 

From 5f608ed572d33122f21c7dd3f0c2e27e6429f904 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Fri, 9 Apr 2021 08:25:30 +0200
Subject: [PATCH 04/24] only 2 and 3 dimension can be 1

---
 src/grid/discretization_grid.f90 | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90
index bb6fa9d8d..88216f0aa 100644
--- a/src/grid/discretization_grid.f90
+++ b/src/grid/discretization_grid.f90
@@ -74,6 +74,8 @@ subroutine discretization_grid_init(restart)
     allocate(materialAt_global(0))                                                                  ! needed for IntelMPI
   endif
 
+  if (grid(1) < 2) call IO_error(844, ext_msg='cells(1) must be larger than  1')
+
   call MPI_Bcast(grid,3,MPI_INTEGER,0,PETSC_COMM_WORLD, ierr)
   if (ierr /= 0) error stop 'MPI error'
   call MPI_Bcast(geomSize,3,MPI_DOUBLE,0,PETSC_COMM_WORLD, ierr)

From 4d9949547c43a917688c47cb621e4b968a9b7507 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sat, 10 Apr 2021 20:52:46 +0200
Subject: [PATCH 05/24] more systematic name

---
 src/CPFEM.f90                    |  2 +-
 src/DAMASK_Marc.f90              | 12 ++++++------
 src/Marc/discretization_Marc.f90 | 20 ++++++++++----------
 3 files changed, 17 insertions(+), 17 deletions(-)

diff --git a/src/CPFEM.f90 b/src/CPFEM.f90
index 4961b99d4..60dc78d2e 100644
--- a/src/CPFEM.f90
+++ b/src/CPFEM.f90
@@ -156,7 +156,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
                                                      ODD_JACOBIAN  = 1e50_pReal                     !< return value for jacobian if terminallyIll
 
 
-  elCP = mesh_FEM2DAMASK_elem(elFE)
+  elCP = discretization_Marc_FEM2DAMASK_elem(elFE)
 
   ma = (elCP-1) * discretization_nIPs + ip
 
diff --git a/src/DAMASK_Marc.f90 b/src/DAMASK_Marc.f90
index e417be2fa..a1ce8fa0d 100644
--- a/src/DAMASK_Marc.f90
+++ b/src/DAMASK_Marc.f90
@@ -266,7 +266,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
       computationMode = CPFEM_RESTOREJACOBIAN
   elseif (lovl == 6) then                                                                           ! stress requested by marc
     computationMode = CPFEM_CALCRESULTS
-    cp_en = mesh_FEM2DAMASK_elem(m(1))
+    cp_en = discretization_Marc_FEM2DAMASK_elem(m(1))
     if (cptim > theTime .or. inc /= theInc) then                                                    ! reached "convergence"
       terminallyIll = .false.
       cycleCounter = -1                                                                             ! first calc step increments this to cycle = 0
@@ -370,11 +370,11 @@ subroutine uedinc(inc,incsub)
 
 
   if (inc > inc_written) then
-    allocate(d_n(3,count(mesh_FEM2DAMASK_node /= -1)))
-    do n = lbound(mesh_FEM2DAMASK_node,1), ubound(mesh_FEM2DAMASK_node,1)
-      if (mesh_FEM2DAMASK_node(n) /= -1) then
-        call nodvar(1,n,d_n(1:3,mesh_FEM2DAMASK_node(n)),nqncomp,nqdatatype)
-        if(nqncomp == 2) d_n(3,mesh_FEM2DAMASK_node(n)) = 0.0_pReal
+    allocate(d_n(3,count(discretization_Marc_FEM2DAMASK_node /= -1)))
+    do n = lbound(discretization_Marc_FEM2DAMASK_node,1), ubound(discretization_Marc_FEM2DAMASK_node,1)
+      if (discretization_Marc_FEM2DAMASK_node(n) /= -1) then
+        call nodvar(1,n,d_n(1:3,discretization_Marc_FEM2DAMASK_node(n)),nqncomp,nqdatatype)
+        if(nqncomp == 2) d_n(3,discretization_Marc_FEM2DAMASK_node(n)) = 0.0_pReal
       endif
     enddo
 
diff --git a/src/Marc/discretization_Marc.f90 b/src/Marc/discretization_Marc.f90
index d92623215..677cf1b52 100644
--- a/src/Marc/discretization_Marc.f90
+++ b/src/Marc/discretization_Marc.f90
@@ -24,8 +24,8 @@ module discretization_marc
     mesh_unitlength                                                                                 !< physical length of one unit in mesh MD: needs systematic_name
 
   integer,  dimension(:), allocatable, public, protected :: &
-    mesh_FEM2DAMASK_elem, &                                                                         !< DAMASK element ID for Marc element ID MD: Needs systematic name
-    mesh_FEM2DAMASK_node                                                                            !< DAMASK node ID for Marc node ID MD: needs systematic_name
+    discretization_Marc_FEM2DAMASK_elem, &                                                          !< DAMASK element ID for Marc element ID
+    discretization_Marc_FEM2DAMASK_node                                                             !< DAMASK node ID for Marc node ID
 
 
   type tCellNodeDefinition
@@ -127,7 +127,7 @@ subroutine discretization_marc_UpdateNodeAndIpCoords(d_n)
   real(pReal), dimension(:,:), allocatable :: node_cell
 
  
-  node_cell = buildCellNodes(discretization_NodeCoords0(1:3,1:maxval(mesh_FEM2DAMASK_node)) + d_n)
+  node_cell = buildCellNodes(discretization_NodeCoords0(1:3,1:maxval(discretization_Marc_FEM2DAMASK_node)) + d_n)
 
   call discretization_setNodeCoords(node_cell)
   call discretization_setIPcoords(buildIPcoordinates(node_cell))
@@ -219,10 +219,10 @@ subroutine inputRead(elem,node0_elem,connectivity_elem,materialAt)
   call inputRead_elemType(elem, &
                           nElems,inputFile)
 
-  call inputRead_mapElems(mesh_FEM2DAMASK_elem,&
+  call inputRead_mapElems(discretization_Marc_FEM2DAMASK_elem,&
                           nElems,elem%nNodes,inputFile)
 
-  call inputRead_mapNodes(mesh_FEM2DAMASK_node,&
+  call inputRead_mapNodes(discretization_Marc_FEM2DAMASK_node,&
                           nNodes,inputFile)
 
   call inputRead_elemNodes(node0_elem, &
@@ -532,7 +532,7 @@ subroutine inputRead_elemNodes(nodes, &
     if(IO_lc(IO_stringValue(fileContent(l),chunkPos,1)) == 'coordinates') then
       chunkPos = [4,1,10,11,30,31,50,51,70]
       do i=1,nNode
-        m = mesh_FEM2DAMASK_node(IO_intValue(fileContent(l+1+i),chunkPos,1))
+        m = discretization_Marc_FEM2DAMASK_node(IO_intValue(fileContent(l+1+i),chunkPos,1))
         do j = 1,3
           nodes(j,m) = mesh_unitlength * IO_floatValue(fileContent(l+1+i),chunkPos,j+1)
         enddo
@@ -653,11 +653,11 @@ function inputRead_connectivityElem(nElem,nNodes,fileContent)
       j = 0
       do i = 1,nElem
         chunkPos = IO_stringPos(fileContent(l+1+i+j))
-        e = mesh_FEM2DAMASK_elem(IO_intValue(fileContent(l+1+i+j),chunkPos,1))
+        e = discretization_Marc_FEM2DAMASK_elem(IO_intValue(fileContent(l+1+i+j),chunkPos,1))
         if (e /= 0) then                                                                            ! disregard non CP elems
           do k = 1,chunkPos(1)-2
             inputRead_connectivityElem(k,e) = &
-              mesh_FEM2DAMASK_node(IO_IntValue(fileContent(l+1+i+j),chunkPos,k+2))
+              discretization_Marc_FEM2DAMASK_node(IO_IntValue(fileContent(l+1+i+j),chunkPos,k+2))
           enddo
           nNodesAlreadyRead = chunkPos(1) - 2
           do while(nNodesAlreadyRead < nNodes)                                                      ! read on if not all nodes in one line
@@ -665,7 +665,7 @@ function inputRead_connectivityElem(nElem,nNodes,fileContent)
             chunkPos = IO_stringPos(fileContent(l+1+i+j))
             do k = 1,chunkPos(1)
               inputRead_connectivityElem(nNodesAlreadyRead+k,e) = &
-                mesh_FEM2DAMASK_node(IO_IntValue(fileContent(l+1+i+j),chunkPos,k))
+                discretization_Marc_FEM2DAMASK_node(IO_IntValue(fileContent(l+1+i+j),chunkPos,k))
             enddo
             nNodesAlreadyRead = nNodesAlreadyRead + chunkPos(1)
           enddo
@@ -718,7 +718,7 @@ subroutine inputRead_material(materialAt,&
           if (initialcondTableStyle == 2) m = m + 2
           contInts = continuousIntValues(fileContent(l+k+m+1:),nElem,nameElemSet,mapElemSet,size(nameElemSet)) ! get affected elements
           do i = 1,contInts(1)
-            e = mesh_FEM2DAMASK_elem(contInts(1+i))
+            e = discretization_Marc_FEM2DAMASK_elem(contInts(1+i))
             materialAt(e) = myVal
           enddo
           if (initialcondTableStyle == 0) m = m + 1

From 690777ac8808a7573b0806b4471c175ccbc02a93 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sat, 10 Apr 2021 21:10:17 +0200
Subject: [PATCH 06/24] base access on cell numbers

DAMASK does not care about elem, IP, etc..
---
 src/Marc/discretization_Marc.f90 | 36 ++++++++++++++++++++++++--------
 1 file changed, 27 insertions(+), 9 deletions(-)

diff --git a/src/Marc/discretization_Marc.f90 b/src/Marc/discretization_Marc.f90
index 677cf1b52..c9f099758 100644
--- a/src/Marc/discretization_Marc.f90
+++ b/src/Marc/discretization_Marc.f90
@@ -39,8 +39,9 @@ module discretization_marc
     connectivity_cell                                                                               !< cell connectivity for each element,ip/cell
 
   public :: &
-    discretization_marc_init, &
-    discretization_marc_UpdateNodeAndIpCoords
+    discretization_Marc_init, &
+    discretization_Marc_updateNodeAndIpCoords, &
+    discretization_Marc_FEM2DAMASK_cell
 
 contains
 
@@ -48,7 +49,7 @@ contains
 !> @brief initializes the mesh by calling all necessary private routines the mesh module
 !! Order and routines strongly depend on type of solver
 !--------------------------------------------------------------------------------------------------
-subroutine discretization_marc_init
+subroutine discretization_Marc_init
 
   real(pReal), dimension(:,:),     allocatable :: &
    node0_elem, &                                                                                    !< node x,y,z coordinates (initially!)
@@ -96,7 +97,7 @@ subroutine discretization_marc_init
   call buildCells(connectivity_cell,cellNodeDefinition,&
                   elem,connectivity_elem)
   node0_cell  = buildCellNodes(node0_elem)
-  
+
   IP_reshaped = buildIPcoordinates(node0_cell)
 
   call discretization_init(materialAt, IP_reshaped, node0_cell)
@@ -114,25 +115,42 @@ subroutine discretization_marc_init
   call geometry_plastic_nonlocal_setIPneighborhood(IPneighborhood(elem))
   call geometry_plastic_nonlocal_results
 
-end subroutine discretization_marc_init
+end subroutine discretization_Marc_init
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief Calculate and set current nodal and IP positions (including cell nodes)
 !--------------------------------------------------------------------------------------------------
-subroutine discretization_marc_UpdateNodeAndIpCoords(d_n)
- 
+subroutine discretization_Marc_updateNodeAndIpCoords(d_n)
+
   real(pReal), dimension(:,:), intent(in)  :: d_n
 
   real(pReal), dimension(:,:), allocatable :: node_cell
 
- 
+
   node_cell = buildCellNodes(discretization_NodeCoords0(1:3,1:maxval(discretization_Marc_FEM2DAMASK_node)) + d_n)
 
   call discretization_setNodeCoords(node_cell)
   call discretization_setIPcoords(buildIPcoordinates(node_cell))
 
-end subroutine discretization_marc_UpdateNodeAndIpCoords
+end subroutine discretization_Marc_updateNodeAndIpCoords
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Calculate and set current nodal and IP positions (including cell nodes)
+!--------------------------------------------------------------------------------------------------
+function discretization_marc_FEM2DAMASK_cell(IP_FEM,elem_FEM) result(cell)
+
+  integer, intent(in) :: IP_FEM, elem_FEM
+  integer :: cell
+
+  real(pReal), dimension(:,:), allocatable :: node_cell
+
+
+  cell = (discretization_Marc_FEM2DAMASK_elem(elem_FEM)-1)*discretization_nIPs + IP_FEM
+
+
+end function discretization_marc_FEM2DAMASK_cell
 
 
 !--------------------------------------------------------------------------------------------------

From 97d426718ad7b6aa3a515645daedae5e287c1808 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sat, 10 Apr 2021 21:16:57 +0200
Subject: [PATCH 07/24] following renames and access pattern

---
 src/DAMASK_Marc.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/DAMASK_Marc.f90 b/src/DAMASK_Marc.f90
index a1ce8fa0d..43ec9b084 100644
--- a/src/DAMASK_Marc.f90
+++ b/src/DAMASK_Marc.f90
@@ -344,8 +344,8 @@ subroutine flux(f,ts,n,time)
   real(pReal), dimension(2),           intent(out) :: &
     f
 
+  f(1) = homogenization_f_T(discretization_Marc_FEM2DAMASK_cell(n(3),n(1)))
   f(2) = 0.0_pReal
-  call thermal_conduction_getSource(f(1), n(3),mesh_FEM2DAMASK_elem(n(1)))
 
  end subroutine flux
 

From 1133090b6cb1bfb2d8fa62d7f18c5e92c8f52ecb Mon Sep 17 00:00:00 2001
From: Martin Diehl <m.diehl@mpie.de>
Date: Sat, 10 Apr 2021 22:47:00 +0200
Subject: [PATCH 08/24] logfile does not contain much valuable information

Marc automatically creates .out
---
 PRIVATE                       |  2 +-
 python/damask/solver/_marc.py | 29 +++++++++--------------------
 2 files changed, 10 insertions(+), 21 deletions(-)

diff --git a/PRIVATE b/PRIVATE
index 129812414..1490f9741 160000
--- a/PRIVATE
+++ b/PRIVATE
@@ -1 +1 @@
-Subproject commit 1298124143e7e2901d0b9c2e79ab6388cb78a1e3
+Subproject commit 1490f97417664d6ae11d7ceafb6b98c9cc2dade1
diff --git a/python/damask/solver/_marc.py b/python/damask/solver/_marc.py
index 14e21144a..26823911d 100644
--- a/python/damask/solver/_marc.py
+++ b/python/damask/solver/_marc.py
@@ -41,13 +41,9 @@ class Marc:
         return path_tools
 
 
-    def submit_job(self,
-                   model,
-                   job          = 'job1',
-                   logfile      = False,
+    def submit_job(self, model, job,
                    compile      = False,
-                   optimization = '',
-                  ):
+                   optimization = ''):
 
         usersub = Path(os.environ['DAMASK_ROOT'])/'src/DAMASK_Marc'
         usersub = usersub.parent/(usersub.name + ('.f90' if compile else '.marc'))
@@ -64,22 +60,15 @@ class Marc:
                ' -prog ' + str(usersub.with_suffix(''))
         print(cmd)
 
-        if logfile is not None:
-            try:
-                f = open(logfile,'w+',newline='\n')
-            except TypeError:
-                f = logfile
-        else:
-            f = io.StringIO()
-
-        proc = subprocess.Popen(shlex.split(cmd),stdout=f,stderr=subprocess.STDOUT)
-        proc.wait()
-        f.seek(0)
+        ret = subprocess.run(shlex.split(cmd),capture_output=True)
 
         try:
-            v = int(re.search('Exit number ([0-9]+)',''.join(f.readlines())).group(1))
+            if 3004 != int(re.search('Exit number ([0-9]+)',ret.stderr.decode()).group(1)):
+                print(ret.stderr.decode())
+                print(ret.stdout.decode())
+                raise RuntimeError(f'Marc simulation failed ({v})')
         except (AttributeError,ValueError):
+            print(ret.stderr.decode())
+            print(ret.stdout.decode())
             raise RuntimeError('Marc simulation failed (unknown return value)')
 
-        if v != 3004:
-            raise RuntimeError(f'Marc simulation failed ({v})')

From cfbb2d416f1f6638734a864858660ec7aa8567f9 Mon Sep 17 00:00:00 2001
From: Martin Diehl <m.diehl@mpie.de>
Date: Sat, 10 Apr 2021 23:52:59 +0200
Subject: [PATCH 09/24] better example

---
 examples/.gitignore             |   3 +
 examples/Marc/material.config   | 429 --------------------------------
 examples/Marc/material.yaml     |  31 +++
 examples/Marc/rotation.mud      | Bin 20700 -> 0 bytes
 examples/Marc/sheet_r-value.dat | 406 ++++++++++++++++++++++++++++++
 python/damask/solver/_marc.py   |   4 +-
 6 files changed, 442 insertions(+), 431 deletions(-)
 delete mode 100644 examples/Marc/material.config
 create mode 100644 examples/Marc/material.yaml
 delete mode 100644 examples/Marc/rotation.mud
 create mode 100644 examples/Marc/sheet_r-value.dat

diff --git a/examples/.gitignore b/examples/.gitignore
index 93d78295b..b832fb271 100644
--- a/examples/.gitignore
+++ b/examples/.gitignore
@@ -3,3 +3,6 @@
 *.xdmf
 *.sta
 *.vt*
+*.out
+*.sts
+*.t16
diff --git a/examples/Marc/material.config b/examples/Marc/material.config
deleted file mode 100644
index 46ea44367..000000000
--- a/examples/Marc/material.config
+++ /dev/null
@@ -1,429 +0,0 @@
-#-------------------#
-<homogenization>
-#-------------------#
-
-{../ConfigFiles/Homogenization_None_Dummy.config}
-
-#-------------------#
-<microstructure>
-#-------------------#
-
-[Grain001]
-(constituent)  phase 1   texture   1   fraction 1.0
-[Grain002]
-(constituent)  phase 1   texture   2   fraction 1.0
-[Grain003]
-(constituent)  phase 1   texture   3   fraction 1.0
-[Grain004]
-(constituent)  phase 1   texture   4   fraction 1.0
-[Grain005]
-(constituent)  phase 1   texture   5   fraction 1.0
-[Grain006]
-(constituent)  phase 1   texture   6   fraction 1.0
-[Grain007]
-(constituent)  phase 1   texture   7   fraction 1.0
-[Grain008]
-(constituent)  phase 1   texture   8   fraction 1.0
-[Grain009]
-(constituent)  phase 1   texture   9   fraction 1.0
-[Grain010]
-(constituent)  phase 1   texture  10   fraction 1.0
-[Grain011]
-(constituent)  phase 1   texture  11   fraction 1.0
-[Grain012]
-(constituent)  phase 1   texture  12   fraction 1.0
-[Grain013]
-(constituent)  phase 1   texture  13   fraction 1.0
-[Grain014]
-(constituent)  phase 1   texture  14   fraction 1.0
-[Grain015]
-(constituent)  phase 1   texture  15   fraction 1.0
-[Grain016]
-(constituent)  phase 1   texture  16   fraction 1.0
-[Grain017]
-(constituent)  phase 1   texture  17   fraction 1.0
-[Grain018]
-(constituent)  phase 1   texture  18   fraction 1.0
-[Grain019]
-(constituent)  phase 1   texture  19   fraction 1.0
-[Grain020]
-(constituent)  phase 1   texture  20   fraction 1.0
-[Grain021]
-(constituent)  phase 1   texture  21   fraction 1.0
-[Grain022]
-(constituent)  phase 1   texture  22   fraction 1.0
-[Grain023]
-(constituent)  phase 1   texture  23   fraction 1.0
-[Grain024]
-(constituent)  phase 1   texture  24   fraction 1.0
-[Grain025]
-(constituent)  phase 1   texture  25   fraction 1.0
-[Grain026]
-(constituent)  phase 1   texture  26   fraction 1.0
-[Grain027]
-(constituent)  phase 1   texture  27   fraction 1.0
-[Grain028]
-(constituent)  phase 1   texture  28   fraction 1.0
-[Grain029]
-(constituent)  phase 1   texture  29   fraction 1.0
-[Grain030]
-(constituent)  phase 1   texture  30   fraction 1.0
-[Grain031]
-(constituent)  phase 1   texture  31   fraction 1.0
-[Grain032]
-(constituent)  phase 1   texture  32   fraction 1.0
-[Grain033]
-(constituent)  phase 1   texture  33   fraction 1.0
-[Grain034]
-(constituent)  phase 1   texture  34   fraction 1.0
-[Grain035]
-(constituent)  phase 1   texture  35   fraction 1.0
-[Grain036]
-(constituent)  phase 1   texture  36   fraction 1.0
-[Grain037]
-(constituent)  phase 1   texture  37   fraction 1.0
-[Grain038]
-(constituent)  phase 1   texture  38   fraction 1.0
-[Grain039]
-(constituent)  phase 1   texture  39   fraction 1.0
-[Grain040]
-(constituent)  phase 1   texture  40   fraction 1.0
-[Grain041]
-(constituent)  phase 1   texture  41   fraction 1.0
-[Grain042]
-(constituent)  phase 1   texture  42   fraction 1.0
-[Grain043]
-(constituent)  phase 1   texture  43   fraction 1.0
-[Grain044]
-(constituent)  phase 1   texture  44   fraction 1.0
-[Grain045]
-(constituent)  phase 1   texture  45   fraction 1.0
-[Grain046]
-(constituent)  phase 1   texture  46   fraction 1.0
-[Grain047]
-(constituent)  phase 1   texture  47   fraction 1.0
-[Grain048]
-(constituent)  phase 1   texture  48   fraction 1.0
-[Grain049]
-(constituent)  phase 1   texture  49   fraction 1.0
-[Grain050]
-(constituent)  phase 1   texture  50   fraction 1.0
-[Grain051]
-(constituent)  phase 1   texture  51   fraction 1.0
-[Grain052]
-(constituent)  phase 1   texture  52   fraction 1.0
-[Grain053]
-(constituent)  phase 1   texture  53   fraction 1.0
-[Grain054]
-(constituent)  phase 1   texture  54   fraction 1.0
-[Grain055]
-(constituent)  phase 1   texture  55   fraction 1.0
-[Grain056]
-(constituent)  phase 1   texture  56   fraction 1.0
-[Grain057]
-(constituent)  phase 1   texture  57   fraction 1.0
-[Grain058]
-(constituent)  phase 1   texture  58   fraction 1.0
-[Grain059]
-(constituent)  phase 1   texture  59   fraction 1.0
-[Grain060]
-(constituent)  phase 1   texture  60   fraction 1.0
-[Grain061]
-(constituent)  phase 1   texture  61   fraction 1.0
-[Grain062]
-(constituent)  phase 1   texture  62   fraction 1.0
-[Grain063]
-(constituent)  phase 1   texture  63   fraction 1.0
-[Grain064]
-(constituent)  phase 1   texture  64   fraction 1.0
-[Grain065]
-(constituent)  phase 1   texture  65   fraction 1.0
-[Grain066]
-(constituent)  phase 1   texture  66   fraction 1.0
-[Grain067]
-(constituent)  phase 1   texture  67   fraction 1.0
-[Grain068]
-(constituent)  phase 1   texture  68   fraction 1.0
-[Grain069]
-(constituent)  phase 1   texture  69   fraction 1.0
-[Grain070]
-(constituent)  phase 1   texture  70   fraction 1.0
-[Grain071]
-(constituent)  phase 1   texture  71   fraction 1.0
-[Grain072]
-(constituent)  phase 1   texture  72   fraction 1.0
-[Grain073]
-(constituent)  phase 1   texture  73   fraction 1.0
-[Grain074]
-(constituent)  phase 1   texture  74   fraction 1.0
-[Grain075]
-(constituent)  phase 1   texture  75   fraction 1.0
-[Grain076]
-(constituent)  phase 1   texture  76   fraction 1.0
-[Grain077]
-(constituent)  phase 1   texture  77   fraction 1.0
-[Grain078]
-(constituent)  phase 1   texture  78   fraction 1.0
-[Grain079]
-(constituent)  phase 1   texture  79   fraction 1.0
-[Grain080]
-(constituent)  phase 1   texture  80   fraction 1.0
-[Grain081]
-(constituent)  phase 1   texture  81   fraction 1.0
-[Grain082]
-(constituent)  phase 1   texture  82   fraction 1.0
-[Grain083]
-(constituent)  phase 1   texture  83   fraction 1.0
-[Grain084]
-(constituent)  phase 1   texture  84   fraction 1.0
-[Grain085]
-(constituent)  phase 1   texture  85   fraction 1.0
-[Grain086]
-(constituent)  phase 1   texture  86   fraction 1.0
-[Grain087]
-(constituent)  phase 1   texture  87   fraction 1.0
-[Grain088]
-(constituent)  phase 1   texture  88   fraction 1.0
-[Grain089]
-(constituent)  phase 1   texture  89   fraction 1.0
-[Grain090]
-(constituent)  phase 1   texture  90   fraction 1.0
-[Grain091]
-(constituent)  phase 1   texture  91   fraction 1.0
-[Grain092]
-(constituent)  phase 1   texture  92   fraction 1.0
-[Grain093]
-(constituent)  phase 1   texture  93   fraction 1.0
-[Grain094]
-(constituent)  phase 1   texture  94   fraction 1.0
-[Grain095]
-(constituent)  phase 1   texture  95   fraction 1.0
-[Grain096]
-(constituent)  phase 1   texture  96   fraction 1.0
-[Grain097]
-(constituent)  phase 1   texture  97   fraction 1.0
-[Grain098]
-(constituent)  phase 1   texture  98   fraction 1.0
-[Grain099]
-(constituent)  phase 1   texture  99   fraction 1.0
-[Grain100]
-(constituent)  phase 1   texture 100   fraction 1.0
-[cubeGrain]
-(constituent)  phase 1   texture 101   fraction 1.0
-
-#-------------------#
-<texture>
-#-------------------#
-
-[Grain001]
-(gauss)  phi1 359.121452   Phi  82.319471   Phi2 347.729535
-[Grain002]
-(gauss)  phi1 269.253967   Phi 105.379919   Phi2 173.029284
-[Grain003]
-(gauss)  phi1  26.551535   Phi 171.606752   Phi2 124.949264
-[Grain004]
-(gauss)  phi1 123.207774   Phi 124.339577   Phi2  47.937748
-[Grain005]
-(gauss)  phi1 324.188825   Phi 103.089216   Phi2 160.373624
-[Grain006]
-(gauss)  phi1 238.295585   Phi 165.416882   Phi2 234.307741
-[Grain007]
-(gauss)  phi1 232.707177   Phi 110.733726   Phi2 308.049265
-[Grain008]
-(gauss)  phi1 144.463291   Phi 125.891441   Phi2 348.674207
-[Grain009]
-(gauss)  phi1 215.423832   Phi  69.759502   Phi2 164.477632
-[Grain010]
-(gauss)  phi1 118.805444   Phi 143.057031   Phi2 271.963190
-[Grain011]
-(gauss)  phi1 218.049576   Phi  64.017550   Phi2 323.040457
-[Grain012]
-(gauss)  phi1 236.962483   Phi 134.312093   Phi2 220.433366
-[Grain013]
-(gauss)  phi1 352.317686   Phi   3.356527   Phi2  92.447275
-[Grain014]
-(gauss)  phi1 198.311545   Phi  71.452240   Phi2 199.441849
-[Grain015]
-(gauss)  phi1 351.993635   Phi  36.500987   Phi2 236.852886
-[Grain016]
-(gauss)  phi1 262.389063   Phi 101.249950   Phi2 334.305959
-[Grain017]
-(gauss)  phi1  53.220668   Phi  69.570254   Phi2 277.061151
-[Grain018]
-(gauss)  phi1 122.156119   Phi 140.207051   Phi2 221.172906
-[Grain019]
-(gauss)  phi1 295.422170   Phi  26.595511   Phi2 263.206315
-[Grain020]
-(gauss)  phi1 179.137406   Phi 104.500977   Phi2 151.742108
-[Grain021]
-(gauss)  phi1 199.045094   Phi   5.228899   Phi2 356.542109
-[Grain022]
-(gauss)  phi1 268.671476   Phi  24.835403   Phi2  33.578889
-[Grain023]
-(gauss)  phi1 264.248527   Phi  59.766630   Phi2 340.865462
-[Grain024]
-(gauss)  phi1 254.223491   Phi  51.125301   Phi2 201.094027
-[Grain025]
-(gauss)  phi1  22.214008   Phi  92.248774   Phi2 215.168318
-[Grain026]
-(gauss)  phi1  49.511491   Phi  79.933539   Phi2 187.188575
-[Grain027]
-(gauss)  phi1 318.916204   Phi 113.102650   Phi2 241.076629
-[Grain028]
-(gauss)  phi1 239.378433   Phi  89.578655   Phi2  94.167043
-[Grain029]
-(gauss)  phi1  27.561421   Phi 142.892093   Phi2 197.735666
-[Grain030]
-(gauss)  phi1 135.210581   Phi 165.859834   Phi2 285.449561
-[Grain031]
-(gauss)  phi1 223.515916   Phi  56.824378   Phi2 343.289074
-[Grain032]
-(gauss)  phi1  41.127974   Phi 111.289145   Phi2 214.855145
-[Grain033]
-(gauss)  phi1  17.335045   Phi 140.496745   Phi2  77.747371
-[Grain034]
-(gauss)  phi1  36.206421   Phi 148.574232   Phi2  88.870226
-[Grain035]
-(gauss)  phi1 159.618336   Phi 125.680504   Phi2 204.119403
-[Grain036]
-(gauss)  phi1   8.752464   Phi  99.173166   Phi2 143.227089
-[Grain037]
-(gauss)  phi1 351.570753   Phi  67.343218   Phi2   1.779612
-[Grain038]
-(gauss)  phi1  46.771572   Phi 155.018674   Phi2 302.319987
-[Grain039]
-(gauss)  phi1 244.255976   Phi  80.566566   Phi2 264.069331
-[Grain040]
-(gauss)  phi1  41.775388   Phi  47.109507   Phi2 300.598550
-[Grain041]
-(gauss)  phi1 268.753103   Phi  46.654050   Phi2 190.382041
-[Grain042]
-(gauss)  phi1 239.574480   Phi  62.517793   Phi2 147.817535
-[Grain043]
-(gauss)  phi1 128.059775   Phi  61.916743   Phi2 169.674359
-[Grain044]
-(gauss)  phi1 166.545156   Phi  58.709099   Phi2 252.885391
-[Grain045]
-(gauss)  phi1  92.867691   Phi  28.906456   Phi2 164.197290
-[Grain046]
-(gauss)  phi1 291.056147   Phi  35.145174   Phi2 250.155599
-[Grain047]
-(gauss)  phi1  79.015862   Phi  44.772479   Phi2 267.982808
-[Grain048]
-(gauss)  phi1 108.400702   Phi  69.883075   Phi2 222.737053
-[Grain049]
-(gauss)  phi1 348.326500   Phi  11.339714   Phi2 121.682346
-[Grain050]
-(gauss)  phi1 331.476209   Phi 108.775043   Phi2 335.139671
-[Grain051]
-(gauss)  phi1 196.750278   Phi  93.955106   Phi2  63.689075
-[Grain052]
-(gauss)  phi1 136.077875   Phi 130.508342   Phi2 128.468976
-[Grain053]
-(gauss)  phi1 239.643513   Phi  76.284643   Phi2 168.821008
-[Grain054]
-(gauss)  phi1 113.850670   Phi 117.531757   Phi2  71.971648
-[Grain055]
-(gauss)  phi1 149.554071   Phi  16.543098   Phi2 195.556172
-[Grain056]
-(gauss)  phi1  46.626579   Phi  52.447846   Phi2 304.495569
-[Grain057]
-(gauss)  phi1 255.251821   Phi  86.678048   Phi2 238.982712
-[Grain058]
-(gauss)  phi1 324.266133   Phi  28.075458   Phi2  41.191295
-[Grain059]
-(gauss)  phi1 312.000332   Phi  74.648725   Phi2  87.403581
-[Grain060]
-(gauss)  phi1  57.742481   Phi 163.241519   Phi2  68.491438
-[Grain061]
-(gauss)  phi1 112.442447   Phi  51.735320   Phi2 206.538656
-[Grain062]
-(gauss)  phi1 297.453842   Phi 115.283041   Phi2  57.785319
-[Grain063]
-(gauss)  phi1 119.132681   Phi 117.923565   Phi2 196.121206
-[Grain064]
-(gauss)  phi1 199.267314   Phi 163.091476   Phi2  53.549301
-[Grain065]
-(gauss)  phi1  37.765215   Phi  76.795488   Phi2 146.264753
-[Grain066]
-(gauss)  phi1 324.550183   Phi  27.665150   Phi2  56.383148
-[Grain067]
-(gauss)  phi1 337.305377   Phi 136.807151   Phi2 133.661586
-[Grain068]
-(gauss)  phi1 115.744041   Phi  64.536978   Phi2 262.694800
-[Grain069]
-(gauss)  phi1 136.293403   Phi  48.862462   Phi2 343.319175
-[Grain070]
-(gauss)  phi1 111.030931   Phi  80.823213   Phi2  84.041594
-[Grain071]
-(gauss)  phi1 303.985249   Phi 118.929631   Phi2 302.307709
-[Grain072]
-(gauss)  phi1 193.556259   Phi  75.928015   Phi2 176.696899
-[Grain073]
-(gauss)  phi1 102.543259   Phi 121.929923   Phi2 234.496773
-[Grain074]
-(gauss)  phi1 218.581323   Phi 101.753894   Phi2 305.566089
-[Grain075]
-(gauss)  phi1 229.542114   Phi 118.839215   Phi2 129.179156
-[Grain076]
-(gauss)  phi1 202.258840   Phi 139.205956   Phi2 352.248979
-[Grain077]
-(gauss)  phi1 137.954289   Phi  63.806918   Phi2 128.975049
-[Grain078]
-(gauss)  phi1 327.557366   Phi  84.987420   Phi2 345.483143
-[Grain079]
-(gauss)  phi1 334.610243   Phi  74.535474   Phi2 106.419231
-[Grain080]
-(gauss)  phi1  62.906243   Phi  46.752029   Phi2 222.692276
-[Grain081]
-(gauss)  phi1 254.121439   Phi 121.005485   Phi2 287.265977
-[Grain082]
-(gauss)  phi1 140.765045   Phi 141.268031   Phi2 271.327656
-[Grain083]
-(gauss)  phi1  10.726984   Phi  66.339177   Phi2 189.073212
-[Grain084]
-(gauss)  phi1 270.921536   Phi  72.821127   Phi2 313.590515
-[Grain085]
-(gauss)  phi1 299.059668   Phi  23.884874   Phi2  80.016277
-[Grain086]
-(gauss)  phi1 208.617406   Phi  11.031834   Phi2 302.388247
-[Grain087]
-(gauss)  phi1  62.929967   Phi  65.223261   Phi2 108.558265
-[Grain088]
-(gauss)  phi1   9.014959   Phi  33.542169   Phi2 247.970366
-[Grain089]
-(gauss)  phi1 272.432808   Phi  30.065174   Phi2  19.803570
-[Grain090]
-(gauss)  phi1 179.621980   Phi 151.763475   Phi2  61.871794
-[Grain091]
-(gauss)  phi1 247.810321   Phi 112.752980   Phi2 264.668469
-[Grain092]
-(gauss)  phi1 270.780630   Phi 102.037858   Phi2  31.602610
-[Grain093]
-(gauss)  phi1  17.626672   Phi  56.032415   Phi2 245.079600
-[Grain094]
-(gauss)  phi1 112.165186   Phi  87.390459   Phi2 182.086729
-[Grain095]
-(gauss)  phi1 157.869381   Phi  79.905131   Phi2 107.037081
-[Grain096]
-(gauss)  phi1 106.163846   Phi 148.477084   Phi2 350.980466
-[Grain097]
-(gauss)  phi1 262.138550   Phi  58.923588   Phi2 111.303439
-[Grain098]
-(gauss)  phi1  88.739397   Phi 119.092789   Phi2 222.502594
-[Grain099]
-(gauss)  phi1 337.603765   Phi  10.145102   Phi2  80.934916
-[Grain100]
-(gauss)  phi1 341.022242   Phi  45.927285   Phi2 252.045476
-
-[cube]
-(gauss)  phi1 0   Phi  0   phi2 0
-
-
-#-------------------#
-<phase>
-#-------------------#
-
-{../ConfigFiles/Phase_Phenopowerlaw_Aluminum.config}
-
-{../ConfigFiles/Phase_Isotropic_AluminumIsotropic.config}
diff --git a/examples/Marc/material.yaml b/examples/Marc/material.yaml
new file mode 100644
index 000000000..5c893ac52
--- /dev/null
+++ b/examples/Marc/material.yaml
@@ -0,0 +1,31 @@
+---
+homogenization:
+  SX:
+    N_constituents: 1
+    mechanical: {type: pass}
+
+phase:
+  Aluminum:
+    lattice: cF
+    mechanical:
+      output: [F, P, F_e, F_p, L_p, O]
+      elastic: {type: Hooke, C_11: 106.75e9, C_12: 60.41e9, C_44: 28.34e9}
+      plastic:
+        type: phenopowerlaw
+        N_sl: [12]
+        a_sl: 2.25
+        atol_xi: 1.0
+        dot_gamma_0_sl: 0.001
+        h_0_sl_sl: 75e6
+        h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4]
+        n_sl: 20
+        output: [xi_sl]
+        xi_0_sl: [31e6]
+        xi_inf_sl: [63e6]
+
+material:
+  - homogenization: SX
+    constituents:
+      - phase: Aluminum
+        v: 1.0
+        O: [0.9330127018922194, 0.25, 0.06698729810778066, 0.25]
diff --git a/examples/Marc/rotation.mud b/examples/Marc/rotation.mud
deleted file mode 100644
index 21eff797487e1eb32b16d7463f7543ecb34797f0..0000000000000000000000000000000000000000
GIT binary patch
literal 0
HcmV?d00001

literal 20700
zcmeI434Baf`^RrlVqYq)y@Sxi5}8OMiOQJ}vBa(vZ4t>JWRnDmSVnE32<7!!+t_Mf
zYprr)DOJ)cs<o7eAh9J{wX~M^_sltGW-`enq2B(NkMsF_=Q;a+&+^=x^wdX1heU*H
zyfqyIqk=RY_2DsrF&bAFt(V4M9~d({O0Ve>7&KI))i`>%4G4*GQd$hq4{D)-fmc}I
zppYOZ<zz8MX`v4fZeb2|)i@%nkeCpCw3E_Eu?p9^xN984BZ4WA%xlKqynm2)F7phg
zNG@OLk<DzSDJr9wWH#eg+i7y=CDSNM?l486y1kw&f0KvJG{sF(%=s5ow-UM1$h5ZG
zV`nps?e^ptTvwN77~D{O7zWs4{3zeyfF!k_+C%jtH;pfi7xh2&8HJNye&dc{lgUwj
zGp&clF)%VRG!CVEBM#+O7L*4SKz+~z$Pz5&BlwpCWT-|fz-pj^{onvN2)+S_K??W<
zoCo&+H|Ab`QTM_HWY(BwU*7rH=>`Nw>qVjF?lMmmQ;^A05jNY8hACQsmb~(y;j`?_
z0y8OO1YTBDL`=WfezYp|i!~*b)hMXRxu0xP7KOIRvMdw)V@-}Tl`5N3P?Pgod`y`X
z((?RG&N=ldJ<W*17+#jeD33Cwl1(Y7$@HZEUGoA>ro{n<H=UZBS02VFRzEymj@8T4
z1*=zBV2nO0Brp^USQ(L>%(kX*(}HFX<n#8Wa8p;60+ogPO~|*vATSh!f-n#cB0wY<
z2BJVThyk=rj{qZq`BBSeyJ}6k*>%^Q0v4s5Gl06KajP0P(>$uN8Mn75JvhAC=kc1m
zCp!-4BJy$5p%@r9oeYn!ys^#YQ4`tS!B-NtKTP+L;apbGahWR8SJCmsayOfW0J+#x
zo8rL)FcHXpvXnnyrEd)n+J>ml>HtST=MmaCS^z4X3eo{z&=Rx-?Ld3b7Yqc{zk@*t
z7y_*On)-DdFo5ykV?cGA0usP<kOY>3RbU-h4`^KXf_*>*$zVS?01krhz%g(f`~jYU
zdao&p6KDoJL01q8B0waF0x@7Dhy$ZQJQxkef^onA#)FT+1TYay0+YcMkN~EFY2Xtu
z9n1hTK_Zw1W`j9k8CVWhfR*42unMdOYrtBt4lKgtT?{h84R8x&gS+52@H_YeJOY1$
zzrbVg8V3D!P#ly1cAyj}2g-vApdxq^R0GvP4Nw!j4eA01P!H4x4M7v&0<^#txB*Y#
z1-wCXpaXuOHSh;*Ks(S8bON107tjOr27SN>;6u<C3;;nO80f(e5Dp?hBp3#w!AKAb
z;=m{{4om<O!6YylB!KB)2AByF!CWvOEC36^BCrH31<Sy4umY?B8^A`e32X*i!49w!
z>;k*Nhj>UL!7va7qCpH82}Xf2zyKzI$zUp&4idq9uoN^ZjWGoNWfa97i<k#!20Vcm
z@CMDnyPyT|0SxGXFK7u`0YA_h_=7f}EocYYg8<M0bOfD1XV3+71>Hb*@E+&^J_11?
z3`_%`fazccm<bZWTrdyJ2MfSLummK58CXMRf<!P2%m$x=IbbgM7nleB4d#OdU?Erp
z7K0_=Gw?Y`0!zU%upF!aE5R3F6<7_{fVE&9SPwRUjbIbl47Px+;7hO#YzI5QPOuB?
z248_a;A^lK>;o!D2K&JQz>P+yo4!4^<|cbT?v`tN6)XP=$+x=F?y!YPDKgM;SzJgv
zA+wpx*IEuFQ)!qf9#1dJwU)+ZnW7XBU+xvw=?m-ryXgxqNA?5X%Z}2P17ka#+#?H$
zC&Oj_(*4P~(QDVPy}@`sa|>BsOo8!a{;wqUfbH>n)jzNL=hgi9r=Jh&eL}je<v(4o
zg36QiD=6N7FZ{;zB(L5+PKSl`Pdae(X+Ya;UbfuU>)W2aKcBrc`g{v<S5xJk%{rcR
z*K9hg^{Tiu#^rBf-;w*Uvrfj7GT-<HkLQc8_&jHPbHC?X52e%d<(t~&`Qj&cd%pM<
z`&0h)Hn&eP8hMVh)8YDMG^AW5Q-lknqE%RpLi$TT>9!Ldu8`zM`Dg$)=5pzdPUc-e
zx(!D_F9#ZchM*B>44eTecU^!MP`*up8*m36pc(K4UV!fT%>g&&^8UBEU$vi)2#q@#
zs0)n7K}b2DaK$tl8BY2yGSZD}tZ|-^7M131&U0}p8ypv@4-KSQ!(T^SGurr%%!A4z
zJFb7Gd?|(ur+mpycin%=l(a0!QlqHVHucmv4v825Wg;nCbcjfLRGLt@gBz8I)gbGD
zE}$DAh1O6I2ExH)Fa=Bn(?BAa2Nr^j;445%s{P;?I1YXUKY^bCy&*XZZh<U7imLnI
z0eA@Lt;?U_38;fQ)&&lL6j@%NBj^N3k@X?y3;KgV5CrsK5Eu+diA72)Qdp4^>QgWm
zB!Q)X6i_q?AdGz=2b0{Ay!>QXVWdnirDdAD{4B!@W0!eIDbrZHrN5N2-4)@n3+*LC
zVU&&@2nwfpL5ngyNHl+Gj?hFWod7M)6+lH$3A_PlF}4TQK}}E#ybUy<HlPQz9;grA
z0rUXdTAby=ObaqC^esSM3fNZg_XGZ*4F~`oKu6FCbOtn+-v_-wZ}1`L2ZBH_(1SrB
z6hwm<FdU2ku^<kN0;EJ91I7Xa7!N)M6M%d!=$V)W=$V)g765uS7J<cJ3HTf&fu&#>
zSPoWzm0&ek1J;7|U<07%Y7^KDwtz3eHn1J+0J{M_cVC0Oz?OpeFv5<2Z^2RU9rzyn
z08RjUPEUeU;57INoB{H=JqP=FZ~^=ZE`m!S6<h{s;0j0wSHU%K9b|wka2Nar?twqR
zePI3E{{_Ft;0bsNo&lPeuYn?<C@2PA2gN}NU<XQqQlK;_1ImJOpggDmxZ#_+VWP7y
zcURU~`qL%M`UM3KchDJ?Q3Lr455ts$BAsF-eYjhfsFO+=Ddh_dV}1zF`Qk$yHikS*
zHii#*$n#Mgg!I$NdKtqN>v~&<`-ybfrWY+K$H`<fe%O|8lVv`#oZcdxn<+h&r%y={
z{l&x0{da;<IA#Esy4pIU6EASMH23GQgTs4{85Zs%{d3BPG9Q)u%k*aX-j$V$4kc;$
zI25Q9R7@{|a+YK3P?yT8Z>d!wVB`#{$L8>hTYTmC$V*L*e01Y<+^txROH1KzZnqqN
z9^V|!-PZMy^GW(^O!<lUynir%yr;XQ&5Ubpqd#SOThdG9EBe>aRk$k)Df>?<AEm#S
zPQHlCd1*@*>8#6>^NIVL=Z$qdTW&=+ZN{kb&Fr`PI9F_YcSP*+7791KJd1JQ-ZG!Q
zJlx_^lQUe#D_F{T@S^EoRyx_=Qu2I*bzEA{8$01IWx?&_`MmsiyuN3}^Znxa7ZVMo
zy}XvN)3RO!6a^Gs0+a%DMyUYa0B?e-paysw)B*JXojn=?XP^acfX*D=fF_v^v;zL1
z9iTG>oh7(wsWVE$2rf<Q^ozI8nBPm7uP-|vFPC8bAf73HsmIAsFBjH0AYw^CqhnXt
zME~V@<&i(l1Um60fLq`;pi@dV(4nDpdl(Po4MN@=mWfUvA*?y1wI8CHj0~gM+Cw^E
zr`O}sG?5m_aTfNMZIw^a%g*=zN`8Dr>D;=zwwG4v+@noTt#Mj6i42Sx+$1KViJxyr
z-|hkZ!tk*xK5}(&9_Zy_Zd#6OwW~{0jbpT)K5>lEhehh6j9<g{lOH~sQqdq=BLnH{
zJ>wT`(bC7bq>;lii43&(Hi~};XHNFIh$E9pH+fqmi?b@NxK%RVJ>1hG*(j@I^i-3P
zrO|pU4Nu0uzm$H)C4iOzN<d|hq(pKNhc#bgT#)*O1Vw}go7+t53$5}pGbyce;)$(l
zM2p$quTf|aM$lX%x-ZinL^YAXYBXUXK~WL?v@Y7*Q=|+c$mmjvR5`wz*1Eewl@l5f
z7_2dt>ts~6R7cdeKwaPfXrrqSoB(ZhjR8&Krl1Gt3FuNh5YQq|%V8J@2N56=3<D#<
zGO!$M0CXWe2x#17Ep15|M@n1q3i78fil`1`C@$qE=L}78n!hv|X)4f&(g2hLwA9l;
zR0MB;%D^5_U(unUCZMHI1FZY)9r!r{>c0kn`py%00Xitq={o>)03AVR&;@h_q&({m
zdV@aT1MngE2=oR0Kz|Sj27n+C4D?_nSPj+z+25OB-weJ2d%%8h02~D0fFs~ra1<N|
zDd2lRm-t`6GhlgHw!A3QCAku?2US3IPy^HiwZPjz18M_O#?=S!07u{i8i1z24Y&gj
zK*~CA&>XxAS^yuwfDW_<{vZH!08+8n8TKxK6nx!4ckmvd?XoBM5PSrz6@S6-8wf(c
z5HJ+b*KN@t1`G!y0Db8e3#`9(8w0<wKz{Ky0rrVt5||98fT>^_pfB9Gkq5qZrVFP$
zpxcr>or8|cRB<8SZFMnsu)0`UyL*fDZbq3M@=({*XbxKC)7w-pUY?&wCogo_;;WB|
zjEA?g3YTx$q?GkGx0|OsA=1luBHS>L=NBms%BBOeOkX%@QvqS2aF7a;9EuhTg^nB#
zUSHl0>Nn{xWzGw`Sus&K1%w<I>vlPqF3z%iK3=~wF0Y9GvrZ@TlTwDe2?dk<T1P5i
zB6xo*?Q~LcA={O|%oOFD=ZkfGnIG@x2BvwI^Mb4q<;!|ztGvDK3Q$<&Y~OP$3X<P9
zDcoeedATpH2$1c{nV+oZd)Do=F@AO;o|H1Zl%{^Nk^a9Q&fAy!g3h0V7$*CH_ot|j
z&OFbgzwxlg9r%?CWLcBX^I2%e{{LV!XhiY<;fBze@_7x7AD!~Z(5Z?H`O&G1+;m=U
z1rOTAjmFSZC#9S@Jb>>|=I@&-@OV;X%p(>0w<uy#9h3lH-(Gx{<+zG@(BDi9(no~p
zW1@_IP}mxAUk7zST~H4+0dxeWx?4gS{|nHukR~o28R_`v3t9qsJR1o6AP@$^K?H~d
zbX1N4(I5s42O|I-pJPEB7zN_NXfOth2Ok4&azE~JTjSwY+iM>BBK|m=^|tHqX(xT$
z$<xbaNYQb(JO&gbWxCw;<M|ZiuhY?(XpNC+2(R&E#$Or%I<)6yD2$c@l4Np>bu$|4
z(2($<xOdArDl;N)UK}N(<T49#TN<j43@aw=Bxy3pbe0m{Ft!@3BvZ+*JQ99oE+!8P
zi)k*yFRzj*ttC8Z{w342=EkyJ(oaeXBgo`cvJ97!$IQQEHy+6%N0Wzzh1Npr-loS>
z_LF5#%XkI#w!0~R%E_F{<ZZN2$y5p%*=Ea<EqhwVE2v}-o*jNI+hy_;77DZOZOf8n
zFH6aj<=@++s^*4fzO6P}hRMDmncw_2=Vhw~(%;TF*XcgY4ebf$)Y6wrs)_7XTQ#vB
z0U6)+DCDdO-O9P4RUp5b$kDJ>6FCx=lIJW-v?Y0;XkA1_xGKU{RV+Kqa!N`&rL!de
zD)6O8YhLw81Cn3SvVUc9a%NB(>#-<Q(dK23jvV<FEe9a4qAfdJ#<M*Ilw`|2!Bwku
z(KvGDpp$XupoY-oqy3!hR4WJJhRCHNssz6`KxJSL-UOsls}1S`Zfp(z%k#HTB+EjL
zb*DK~(d5G|wtCtsS}v!Y05YHyEY>PBTm2!^yevs=FtsB&N=Nzf`kAz)<YhE8Yx3#`
zs{e}^>wcgIhs<;Kh0d7hCH~~-6CJ-Y^152i`KGS)FEYY{`^5wf2&I3M5gkH@Qf_1x
zvTN~@!=EBJj?gyIdb%`d>a2Csx_h~LI4L8v&;!eOmgzr$gCH$tt}d46ensB&vq`~c
zW#nU_=y-rUPF4_ZdHk#*+%<&M2w6u+2O;YV=_q6aAsY$lEToH&S|MG9Y$~LiknTcy
z2-!?XPa(a8^cJ$YknalFLP#GWnUFdz4MHlu!Y!oIQn-awS_!w1il1-`sk9buAr*h&
z7E);=+(IgCg<D8yg%mOdspV@g+(Ie=!Y!oILAZrfItsUt@*8$TC+@Fw7U4oFUAS8(
zq(R7JA(gHoTuAvDO0sY(-9)&M1|gG$RJx0BAq_$%3#q&(!i6*lnJlEzLxc-y5HeXv
zrKboN(lpocc42*J{nD3aGcUjMhk9LH%Cdv}eUBVl%VyW?GC1-47S{f0(AG1Zcd%7O
z7tQH^;45Zt7uM>utIFPaIC^&IiT!MW{!w7*;)hs~x$4Mj-iO(tGKpuxLceA0o%gRe
zyZAfy$@Fc`uE&nErG4CQ1U&ejE$eXi>oK!VF#juUekrx$N48Sm^}A*}PO{N=cWAQ?
zpJFqYwEV96nbYiA(S-2TSAJr#`e{`SS!bA=fB%dbkA7y^{ZCil``TG%I9asXqY~%X
z;SU}~G%j<VEnC?7RO@mV*uAQA>pGYFmDyEo5qh%pMb`M{sjIb0Tw>i%Y<ckP*(LUL
z&b}Jj_*8Z`qNjJ&KT=s0rP!&BVVBvqMYG0bq+e$5I9&hz>pp3$<7X+)>im$#_I|!N
zv0}?BY~%QQ>&Jb0g|!KY92Z<SolR}sv*XU$=`6c}|LBOP=}dPuefGqVtL*LLlico{
zyvka=|6{`y-q+Z>%|Ze{OS;DXI5PF}>0;N}loPAZvq9I{z6lkV|9s>+`)K=^O>^sK
zunL~*mqbj+VBLG&YZ7=dgIykcrqh6CnXFTfsU<^aX0p>wCe0gpC6mQ(D)ru2j~lGi
z+*RMkO}W9wUER4n=*$gvu2}kB_qsP(!o)>YZ-m`s-|d@zWz>$F>?=*x=(NW-*?rC7
z_SJoFu`*}ZR?$qj#p-#T+i~;2Ew&=%Xu~+=Hv78V_O%C_-)4<mpUgQNbDJF=IcfNq
zb+=hqX6e%Dmu|CtQ^xw-uXu+gUFhJr+4l|`e|vlRauIh}i7AzamRWd*&AL!Mb=848
z?EG+Nzh5)%uzmWei`SINVmA*>tyRi7i@g(?*5q-AEcWlr9|n#I$zruOw4L_Jge>;g
zAg50q7G<%Pi$7ZBzd4Kj*{hgC_JJ%`QZswb+o!Tvzuz;qo=VMP<)_3)mcEt6qUZIx
z_07F3w!DvTY>~gRSej<XzmuP4u_yFDh@NE`)dN8h(TX&*h7(1EhjHIPgxvc^W9_S*
z18=-^yrW+29K0W1?HtB_(B<Aalr9Sw81C0qm2by9tZ~1UstoCqQZJ^bs<{1Y!pwU?
zs*+%TCZJJ-s`xz^;T}0cH7pOUd3(xORo8Ij>K=C{sD^8wUp?nBRW*G4M@o&FiK-^G
z(f8d>e5xkfKUp5vf1cX^FRvqed={$xe;fYn*pScEr1vL}Zy2>yJ#qKMzVSaTSG%;F
zUfQeQ7pkVx_`S8tuTlG7Y3cGLeVv+gyzkLPsT)<tT?58mQnsje7j`;UUa?iJ^0-#h
zW4*Vj3DXne(`xTflQx$d=J9x^n(Dph&1D(8)e04d6n}PYkE;3VNOZjid)3jmx4$)_
zOtLzEONkzp>+e^4*=NPS-}ZnSdhu?iWAs5aWAf3i!@fAAhK9K9t&;YQx~<yjnZ@fL
zQQeMJnOHmOTebf*ubrI_A5~-hyq1orb4+d9?P2x!$;Z_b!@33g{+6P~K51KHgZ>BA
zbJ5~2x<sB(-I|woH)wxU=atQN?0Nr3we7f{91re0seW2--JG4%PO1IX%<*;moK`n(
zJAL?mqo34`e|>#dQ}m45t9`r9voD-cw~jh^!(;c)Y9+g&iOPas)Lqk#AH5KNR;~0_
zt<2ql=Tw6}NtxXKysE4!mg($qL4ADVjPEjsU)6n&ZbhZoUsN5ZoZYab<Rx`=QcxCl
zN;RRzcthmqR5kcq<GpYGnX3A)55Alcds%hzi75H&!^`TeEswl2<I>ayKb+MUeVnGI
zU;8<@_P8r*;s<e0J&UHR@q3S7WeMr(?V-)~G%R;jUA)$>>~C|gs#~j%=$%mGn)=NQ
z`?BdPuBrRV+D-nm@pbij!)0&m-+5g<8Z_oXqm~(Jqj{YMHciP;!zS)cINB>y4SKi5
z==AhVwOsJk?s3CzsQ0>es5<894b>+rFeUwyo9aBrfRL3{Z>jNhd@g;t>XusNT;J-i
zd)!t%0~QyXdFZy<;PH$xt-9Y)!&AHJtEJsh7aE4lFEK1j4eGh%E_M!;<~ekp`oJjT
zjHe&Dc@z~s#?6E85~eu9BKH{DXG#fj#xdO;gm=F83sW4jSYDFc-CgLP5aK_Eg*q7{
KQNdI{9_c@eI>7D#

diff --git a/examples/Marc/sheet_r-value.dat b/examples/Marc/sheet_r-value.dat
new file mode 100644
index 000000000..529f4d31a
--- /dev/null
+++ b/examples/Marc/sheet_r-value.dat
@@ -0,0 +1,406 @@
+title               r-value
+$....MARC input file produced by Marc Mentat 2019 (64bit)
+$...................................
+$....input file using extended precision
+extended
+$...................................
+sizing                                 0        80       165         0
+alloc                       25
+elements                     7
+version                     14         1         0         1
+table                        0         0         2         1         1         0         0         1
+processor                    1         1         1         0
+$no list
+large stra                   2         1         0         0         0         0         0
+all points
+no echo                      1         2         3         4
+state vars                   3
+end
+$...................
+solver
+         8         1         0         0         0         0         0         0         0         0         0         0         0         0         0         0
+optimize                    11
+connectivity
+         0         0         1         0         1         1         0         0         0
+         1         7         2         5        20        17         1         4        19        16
+         2         7         3         6        21        18         2         5        20        17
+         3         7         5         8        23        20         4         7        22        19
+         4         7         6         9        24        21         5         8        23        20
+         5         7         8        11        26        23         7        10        25        22
+         6         7         9        12        27        24         8        11        26        23
+         7         7        11        14        29        26        10        13        28        25
+         8         7        12        15        30        27        11        14        29        26
+         9         7        17        20        35        32        16        19        34        31
+        10         7        18        21        36        33        17        20        35        32
+        11         7        20        23        38        35        19        22        37        34
+        12         7        21        24        39        36        20        23        38        35
+        13         7        23        26        41        38        22        25        40        37
+        14         7        24        27        42        39        23        26        41        38
+        15         7        26        29        44        41        25        28        43        40
+        16         7        27        30        45        42        26        29        44        41
+        17         7        32        35        50        47        31        34        49        46
+        18         7        33        36        51        48        32        35        50        47
+        19         7        35        38        53        50        34        37        52        49
+        20         7        36        39        54        51        35        38        53        50
+        21         7        38        41        56        53        37        40        55        52
+        22         7        39        42        57        54        38        41        56        53
+        23         7        41        44        59        56        40        43        58        55
+        24         7        42        45        60        57        41        44        59        56
+        25         7        47        50        65        62        46        49        64        61
+        26         7        48        51        66        63        47        50        65        62
+        27         7        50        53        68        65        49        52        67        64
+        28         7        51        54        69        66        50        53        68        65
+        29         7        53        56        71        68        52        55        70        67
+        30         7        54        57        72        69        53        56        71        68
+        31         7        56        59        74        71        55        58        73        70
+        32         7        57        60        75        72        56        59        74        71
+        33         7        62        65        80        77        61        64        79        76
+        34         7        63        66        81        78        62        65        80        77
+        35         7        65        68        83        80        64        67        82        79
+        36         7        66        69        84        81        65        68        83        80
+        37         7        68        71        86        83        67        70        85        82
+        38         7        69        72        87        84        68        71        86        83
+        39         7        71        74        89        86        70        73        88        85
+        40         7        72        75        90        87        71        74        89        86
+        41         7        77        80        95        92        76        79        94        91
+        42         7        78        81        96        93        77        80        95        92
+        43         7        80        83        98        95        79        82        97        94
+        44         7        81        84        99        96        80        83        98        95
+        45         7        83        86       101        98        82        85       100        97
+        46         7        84        87       102        99        83        86       101        98
+        47         7        86        89       104       101        85        88       103       100
+        48         7        87        90       105       102        86        89       104       101
+        49         7        92        95       110       107        91        94       109       106
+        50         7        93        96       111       108        92        95       110       107
+        51         7        95        98       113       110        94        97       112       109
+        52         7        96        99       114       111        95        98       113       110
+        53         7        98       101       116       113        97       100       115       112
+        54         7        99       102       117       114        98       101       116       113
+        55         7       101       104       119       116       100       103       118       115
+        56         7       102       105       120       117       101       104       119       116
+        57         7       107       110       125       122       106       109       124       121
+        58         7       108       111       126       123       107       110       125       122
+        59         7       110       113       128       125       109       112       127       124
+        60         7       111       114       129       126       110       113       128       125
+        61         7       113       116       131       128       112       115       130       127
+        62         7       114       117       132       129       113       116       131       128
+        63         7       116       119       134       131       115       118       133       130
+        64         7       117       120       135       132       116       119       134       131
+        65         7       122       125       140       137       121       124       139       136
+        66         7       123       126       141       138       122       125       140       137
+        67         7       125       128       143       140       124       127       142       139
+        68         7       126       129       144       141       125       128       143       140
+        69         7       128       131       146       143       127       130       145       142
+        70         7       129       132       147       144       128       131       146       143
+        71         7       131       134       149       146       130       133       148       145
+        72         7       132       135       150       147       131       134       149       146
+        73         7       137       140       155       152       136       139       154       151
+        74         7       138       141       156       153       137       140       155       152
+        75         7       140       143       158       155       139       142       157       154
+        76         7       141       144       159       156       140       143       158       155
+        77         7       143       146       161       158       142       145       160       157
+        78         7       144       147       162       159       143       146       161       158
+        79         7       146       149       164       161       145       148       163       160
+        80         7       147       150       165       162       146       149       164       161
+coordinates
+         3       165         0         1
+         1-4.000000000000000+1-1.000000000000000+1-5.000000000000000-1
+         2-4.000000000000000+1-1.000000000000000+1 0.000000000000000+0
+         3-4.000000000000000+1-1.000000000000000+1 5.000000000000000-1
+         4-4.000000000000000+1-5.000000000000000+0-5.000000000000000-1
+         5-4.000000000000000+1-5.000000000000000+0 0.000000000000000+0
+         6-4.000000000000000+1-5.000000000000000+0 5.000000000000000-1
+         7-4.000000000000000+1 0.000000000000000+0-5.000000000000000-1
+         8-4.000000000000000+1 0.000000000000000+0 0.000000000000000+0
+         9-4.000000000000000+1 0.000000000000000+0 5.000000000000000-1
+        10-4.000000000000000+1 5.000000000000000+0-5.000000000000000-1
+        11-4.000000000000000+1 5.000000000000000+0 0.000000000000000+0
+        12-4.000000000000000+1 5.000000000000000+0 5.000000000000000-1
+        13-4.000000000000000+1 9.999999999999996+0-5.000000000000000-1
+        14-4.000000000000000+1 9.999999999999996+0 0.000000000000000+0
+        15-4.000000000000000+1 9.999999999999996+0 5.000000000000000-1
+        16-3.200000000000000+1-1.000000000000000+1-5.000000000000000-1
+        17-3.200000000000000+1-1.000000000000000+1 0.000000000000000+0
+        18-3.200000000000000+1-1.000000000000000+1 5.000000000000000-1
+        19-3.200000000000000+1-5.000000000000000+0-5.000000000000000-1
+        20-3.200000000000000+1-5.000000000000000+0 0.000000000000000+0
+        21-3.200000000000000+1-5.000000000000000+0 5.000000000000000-1
+        22-3.200000000000000+1 0.000000000000000+0-5.000000000000000-1
+        23-3.200000000000000+1 0.000000000000000+0 0.000000000000000+0
+        24-3.200000000000000+1 0.000000000000000+0 5.000000000000000-1
+        25-3.200000000000000+1 5.000000000000000+0-5.000000000000000-1
+        26-3.200000000000000+1 5.000000000000000+0 0.000000000000000+0
+        27-3.200000000000000+1 5.000000000000000+0 5.000000000000000-1
+        28-3.200000000000000+1 1.000000000000000+1-5.000000000000000-1
+        29-3.200000000000000+1 1.000000000000000+1 0.000000000000000+0
+        30-3.200000000000000+1 1.000000000000000+1 5.000000000000000-1
+        31-2.400000000000001+1-1.000000000000000+1-5.000000000000000-1
+        32-2.400000000000001+1-1.000000000000000+1 0.000000000000000+0
+        33-2.400000000000001+1-1.000000000000000+1 5.000000000000000-1
+        34-2.400000000000000+1-5.000000000000000+0-5.000000000000000-1
+        35-2.400000000000000+1-5.000000000000000+0 0.000000000000000+0
+        36-2.400000000000000+1-5.000000000000000+0 5.000000000000000-1
+        37-2.400000000000001+1 0.000000000000000+0-5.000000000000000-1
+        38-2.400000000000001+1 0.000000000000000+0 0.000000000000000+0
+        39-2.400000000000001+1 0.000000000000000+0 5.000000000000000-1
+        40-2.400000000000001+1 5.000000000000002+0-5.000000000000000-1
+        41-2.400000000000001+1 5.000000000000002+0 0.000000000000000+0
+        42-2.400000000000001+1 5.000000000000002+0 5.000000000000000-1
+        43-2.400000000000001+1 9.999999999999996+0-5.000000000000000-1
+        44-2.400000000000001+1 9.999999999999996+0 0.000000000000000+0
+        45-2.400000000000001+1 9.999999999999996+0 5.000000000000000-1
+        46-1.600000000000000+1-1.000000000000000+1-5.000000000000000-1
+        47-1.600000000000000+1-1.000000000000000+1 0.000000000000000+0
+        48-1.600000000000000+1-1.000000000000000+1 5.000000000000000-1
+        49-1.600000000000000+1-5.000000000000000+0-5.000000000000000-1
+        50-1.600000000000000+1-5.000000000000000+0 0.000000000000000+0
+        51-1.600000000000000+1-5.000000000000000+0 5.000000000000000-1
+        52-1.600000000000000+1 0.000000000000000+0-5.000000000000000-1
+        53-1.600000000000000+1 0.000000000000000+0 0.000000000000000+0
+        54-1.600000000000000+1 0.000000000000000+0 5.000000000000000-1
+        55-1.600000000000000+1 4.999999999999998+0-5.000000000000000-1
+        56-1.600000000000000+1 4.999999999999998+0 0.000000000000000+0
+        57-1.600000000000000+1 4.999999999999998+0 5.000000000000000-1
+        58-1.600000000000000+1 9.999999999999996+0-5.000000000000000-1
+        59-1.600000000000000+1 9.999999999999996+0 0.000000000000000+0
+        60-1.600000000000000+1 9.999999999999996+0 5.000000000000000-1
+        61-8.000000000000000+0-1.000000000000000+1-5.000000000000000-1
+        62-8.000000000000000+0-1.000000000000000+1 0.000000000000000+0
+        63-8.000000000000000+0-1.000000000000000+1 5.000000000000000-1
+        64-7.999999999999994+0-5.000000000000000+0-5.000000000000000-1
+        65-7.999999999999994+0-5.000000000000000+0 0.000000000000000+0
+        66-7.999999999999994+0-5.000000000000000+0 5.000000000000000-1
+        67-8.000000000000000+0 0.000000000000000+0-5.000000000000000-1
+        68-8.000000000000000+0 0.000000000000000+0 0.000000000000000+0
+        69-8.000000000000000+0 0.000000000000000+0 5.000000000000000-1
+        70-8.000000000000000+0 5.000000000000000+0-5.000000000000000-1
+        71-8.000000000000000+0 5.000000000000000+0 0.000000000000000+0
+        72-8.000000000000000+0 5.000000000000000+0 5.000000000000000-1
+        73-8.000000000000000+0 9.999999999999996+0-5.000000000000000-1
+        74-8.000000000000000+0 9.999999999999996+0 0.000000000000000+0
+        75-8.000000000000000+0 9.999999999999996+0 5.000000000000000-1
+        76 0.000000000000000+0-1.000000000000000+1-5.000000000000000-1
+        77 0.000000000000000+0-1.000000000000000+1 0.000000000000000+0
+        78 0.000000000000000+0-1.000000000000000+1 5.000000000000000-1
+        79 0.000000000000000+0-5.000000000000000+0-5.000000000000000-1
+        80 0.000000000000000+0-5.000000000000000+0 0.000000000000000+0
+        81 0.000000000000000+0-5.000000000000000+0 5.000000000000000-1
+        82 0.000000000000000+0 0.000000000000000+0-5.000000000000000-1
+        83 0.000000000000000+0 0.000000000000000+0 0.000000000000000+0
+        84 0.000000000000000+0 0.000000000000000+0 5.000000000000000-1
+        85 0.000000000000000+0 5.000000000000000+0-5.000000000000000-1
+        86 0.000000000000000+0 5.000000000000000+0 0.000000000000000+0
+        87 0.000000000000000+0 5.000000000000000+0 5.000000000000000-1
+        88 0.000000000000000+0 9.999999999999996+0-5.000000000000000-1
+        89 0.000000000000000+0 9.999999999999996+0 0.000000000000000+0
+        90 0.000000000000000+0 9.999999999999996+0 5.000000000000000-1
+        91 8.000000000000000+0-1.000000000000000+1-5.000000000000000-1
+        92 8.000000000000000+0-1.000000000000000+1 0.000000000000000+0
+        93 8.000000000000000+0-1.000000000000000+1 5.000000000000000-1
+        94 8.000000000000000+0-5.000000000000000+0-5.000000000000000-1
+        95 8.000000000000000+0-5.000000000000000+0 0.000000000000000+0
+        96 8.000000000000000+0-5.000000000000000+0 5.000000000000000-1
+        97 8.000000000000000+0 0.000000000000000+0-5.000000000000000-1
+        98 8.000000000000000+0 0.000000000000000+0 0.000000000000000+0
+        99 8.000000000000000+0 0.000000000000000+0 5.000000000000000-1
+       100 8.000000000000000+0 5.000000000000000+0-5.000000000000000-1
+       101 8.000000000000000+0 5.000000000000000+0 0.000000000000000+0
+       102 8.000000000000000+0 5.000000000000000+0 5.000000000000000-1
+       103 7.999999999999994+0 9.999999999999996+0-5.000000000000000-1
+       104 7.999999999999994+0 9.999999999999996+0 0.000000000000000+0
+       105 7.999999999999994+0 9.999999999999996+0 5.000000000000000-1
+       106 1.600000000000000+1-1.000000000000000+1-5.000000000000000-1
+       107 1.600000000000000+1-1.000000000000000+1 0.000000000000000+0
+       108 1.600000000000000+1-1.000000000000000+1 5.000000000000000-1
+       109 1.600000000000000+1-5.000000000000000+0-5.000000000000000-1
+       110 1.600000000000000+1-5.000000000000000+0 0.000000000000000+0
+       111 1.600000000000000+1-5.000000000000000+0 5.000000000000000-1
+       112 1.600000000000000+1 0.000000000000000+0-5.000000000000000-1
+       113 1.600000000000000+1 0.000000000000000+0 0.000000000000000+0
+       114 1.600000000000000+1 0.000000000000000+0 5.000000000000000-1
+       115 1.600000000000000+1 5.000000000000000+0-5.000000000000000-1
+       116 1.600000000000000+1 5.000000000000000+0 0.000000000000000+0
+       117 1.600000000000000+1 5.000000000000000+0 5.000000000000000-1
+       118 1.600000000000000+1 9.999999999999996+0-5.000000000000000-1
+       119 1.600000000000000+1 9.999999999999996+0 0.000000000000000+0
+       120 1.600000000000000+1 9.999999999999996+0 5.000000000000000-1
+       121 2.400000000000000+1-1.000000000000000+1-5.000000000000000-1
+       122 2.400000000000000+1-1.000000000000000+1 0.000000000000000+0
+       123 2.400000000000000+1-1.000000000000000+1 5.000000000000000-1
+       124 2.400000000000001+1-5.000000000000000+0-5.000000000000000-1
+       125 2.400000000000001+1-5.000000000000000+0 0.000000000000000+0
+       126 2.400000000000001+1-5.000000000000000+0 5.000000000000000-1
+       127 2.400000000000000+1 0.000000000000000+0-5.000000000000000-1
+       128 2.400000000000000+1 0.000000000000000+0 0.000000000000000+0
+       129 2.400000000000000+1 0.000000000000000+0 5.000000000000000-1
+       130 2.400000000000001+1 5.000000000000000+0-5.000000000000000-1
+       131 2.400000000000001+1 5.000000000000000+0 0.000000000000000+0
+       132 2.400000000000001+1 5.000000000000000+0 5.000000000000000-1
+       133 2.400000000000000+1 9.999999999999996+0-5.000000000000000-1
+       134 2.400000000000000+1 9.999999999999996+0 0.000000000000000+0
+       135 2.400000000000000+1 9.999999999999996+0 5.000000000000000-1
+       136 3.200000000000000+1-1.000000000000000+1-5.000000000000000-1
+       137 3.200000000000000+1-1.000000000000000+1 0.000000000000000+0
+       138 3.200000000000000+1-1.000000000000000+1 5.000000000000000-1
+       139 3.199999999999999+1-5.000000000000000+0-5.000000000000000-1
+       140 3.199999999999999+1-5.000000000000000+0 0.000000000000000+0
+       141 3.199999999999999+1-5.000000000000000+0 5.000000000000000-1
+       142 3.200000000000000+1 0.000000000000000+0-5.000000000000000-1
+       143 3.200000000000000+1 0.000000000000000+0 0.000000000000000+0
+       144 3.200000000000000+1 0.000000000000000+0 5.000000000000000-1
+       145 3.200000000000000+1 5.000000000000000+0-5.000000000000000-1
+       146 3.200000000000000+1 5.000000000000000+0 0.000000000000000+0
+       147 3.200000000000000+1 5.000000000000000+0 5.000000000000000-1
+       148 3.199999999999999+1 9.999999999999996+0-5.000000000000000-1
+       149 3.199999999999999+1 9.999999999999996+0 0.000000000000000+0
+       150 3.199999999999999+1 9.999999999999996+0 5.000000000000000-1
+       151 3.999999999999999+1-1.000000000000000+1-5.000000000000000-1
+       152 3.999999999999999+1-1.000000000000000+1 0.000000000000000+0
+       153 3.999999999999999+1-1.000000000000000+1 5.000000000000000-1
+       154 3.999999999999999+1-5.000000000000000+0-5.000000000000000-1
+       155 3.999999999999999+1-5.000000000000000+0 0.000000000000000+0
+       156 3.999999999999999+1-5.000000000000000+0 5.000000000000000-1
+       157 3.999999999999999+1 0.000000000000000+0-5.000000000000000-1
+       158 3.999999999999999+1 0.000000000000000+0 0.000000000000000+0
+       159 3.999999999999999+1 0.000000000000000+0 5.000000000000000-1
+       160 3.999999999999999+1 5.000000000000000+0-5.000000000000000-1
+       161 3.999999999999999+1 5.000000000000000+0 0.000000000000000+0
+       162 3.999999999999999+1 5.000000000000000+0 5.000000000000000-1
+       163 3.999999999999999+1 9.999999999999996+0-5.000000000000000-1
+       164 3.999999999999999+1 9.999999999999996+0 0.000000000000000+0
+       165 3.999999999999999+1 9.999999999999996+0 5.000000000000000-1
+define              element             set                 Material_Nummer_elements
+           1          to          80
+define              node                set                 unten_y_nodes
+           2           5           8          11          14
+define              node                set                 oben_y_nodes
+         152         155         158         161         164
+define              node                set                 unten_fest_nodes
+           1          to          15
+define              node                set                 oben_ziehen_nodes
+         151          to         165
+define              node                set                 unten_z_nodes
+           7          to           9
+define              node                set                 oben_z_nodes
+         157          to         159
+define              element             set                 texture_elements
+           1          to          80
+hypoelastic
+
+         1         0         1         0         1TKS                              0
+ 1.000000000000000+0 0.000000000000000+0 0.000000000000000+0 0.000000000000000+0 0.000000000000000+0 0.000000000000000+0 0.000000000000000+0
+         0         0         0         0         0         0         0
+
+mat color
+
+         1         1       230         0         0
+table               weg_x
+         1         1         0         0         2
+         1         2         2         0         0         2         0         0         2         0         0         2
+ 0.000000000000000+0 0.000000000000000+0
+ 2.000000000000000+2 1.600000000000000+1
+geometry
+         0         0         2
+         1         9         1       230         0         0
+r-value-sample
+ 0.000000000000000+0 0.000000000000000+0 0.000000000000000+0 0.000000000000000+0 0.000000000000000+0 0.000000000000000+0 0.000000000000000+0
+
+usdata                       1
+fixed disp
+
+         1         0         0         0         1         0unten_z
+ 0.000000000000000+0 0.000000000000000+0
+         0         0
+         1         2
+         2
+unten_z_nodes
+         1         0         0         0         1         0unten_y
+ 0.000000000000000+0 0.000000000000000+0
+         0         0
+         1         3
+         2
+unten_y_nodes
+         1         0         0         0         1         0oben_z
+ 1.000000000000000+0 0.000000000000000+0
+         1         0
+         1         2
+         2
+oben_z_nodes
+         1         0         0         0         1         0oben_y
+ 1.000000000000000+0 0.000000000000000+0
+         1         0
+         1         3
+         2
+oben_y_nodes
+         1         0         0         0         1         0unten_fest
+ 0.000000000000000+0
+         0
+         1
+         2
+unten_fest_nodes
+         1         0         0         0         1         0oben_ziehen
+ 1.000000000000000+0
+         1
+         1
+         2
+oben_ziehen_nodes
+initial state
+
+         2         6         1         0         0         0Material_Nummer
+ 1.000000000000000+0
+         0
+         1
+Material_Nummer_elements
+initial state
+
+         3         6         1         0         0         0texture
+ 1.000000000000000+0
+         0
+         1
+texture_elements
+loadcase            r-value
+         5
+Material_Nummer
+texture
+unten_z
+unten_y
+unten_fest
+no print
+post
+         6        16        17         0         0        19        20         0         1         0         0         0         0         0         0         0
+parameters
+ 1.000000000000000+0 1.000000000000000+9 1.000000000000000+2 1.000000000000000+6 2.500000000000000-1 5.000000000000000-1 1.500000000000000+0-5.000000000000000-1
+ 8.625000000000000+0 2.000000000000000+1 1.000000000000000-4 1.000000000000000-6 1.000000000000000+0 1.000000000000000-4
+ 8.314000000000000+0 2.731500000000000+2 5.000000000000000-1 0.000000000000000+0 5.670510000000000-8 1.438769000000000-2 2.997900000000000+8 1.00000000000000+30
+ 0.000000000000000+0 0.000000000000000+0 1.000000000000000+2 0.000000000000000+0 1.000000000000000+0-2.000000000000000+0 1.000000000000000+6 3.000000000000000+0
+ 0.000000000000000+0 0.000000000000000+0 1.256637061000000-6 8.85418781700000-12 1.200000000000000+2 1.000000000000000-3 1.600000000000000+2 0.000000000000000+0
+ 3.000000000000000+0 4.000000000000000-1
+end option
+$...................
+$....start of loadcase Tensile
+title               Tensile
+loadcase            Tensile
+         6
+unten_z
+unten_y
+oben_z
+oben_y
+unten_fest
+oben_ziehen
+control
+     99999        10         0         0         0         1         0         0         1         0         0         0         0         0         0
+ 1.000000000000000-1 0.000000000000000+0 0.000000000000000+0 0.000000000000000+0 0.000000000000000+0 0.000000000000000+0-1.000000000000000+0 0.000000000000000+0
+parameters
+ 1.000000000000000+0 1.000000000000000+9 1.000000000000000+2 1.000000000000000+6 2.500000000000000-1 5.000000000000000-1 1.500000000000000+0-5.000000000000000-1
+ 8.625000000000000+0 2.000000000000000+1 1.000000000000000-4 1.000000000000000-6 1.000000000000000+0 1.000000000000000-4
+ 8.314000000000000+0 2.731500000000000+2 5.000000000000000-1 0.000000000000000+0 5.670510000000000-8 1.438769000000000-2 2.997900000000000+8 1.00000000000000+30
+ 0.000000000000000+0 0.000000000000000+0 1.000000000000000+2 0.000000000000000+0 1.000000000000000+0-1.000000000000000+0 1.000000000000000+6 3.000000000000000+0
+ 0.000000000000000+0 0.000000000000000+0 1.256637061000000-6 8.85418781700000-12 1.200000000000000+2 1.000000000000000-3 1.600000000000000+2 0.000000000000000+0
+ 3.000000000000000+0 4.000000000000000-1
+auto load
+       100         0        10         0         0
+time step
+ 2.000000000000000+0
+continue
+$....end of loadcase Tensile
+$...................
diff --git a/python/damask/solver/_marc.py b/python/damask/solver/_marc.py
index 26823911d..9fb07fc1c 100644
--- a/python/damask/solver/_marc.py
+++ b/python/damask/solver/_marc.py
@@ -1,7 +1,6 @@
 import subprocess
 import shlex
 import re
-import io
 import os
 from pathlib import Path
 
@@ -63,7 +62,8 @@ class Marc:
         ret = subprocess.run(shlex.split(cmd),capture_output=True)
 
         try:
-            if 3004 != int(re.search('Exit number ([0-9]+)',ret.stderr.decode()).group(1)):
+            v = int(re.search('Exit number ([0-9]+)',ret.stderr.decode()).group(1))
+            if 3004 != v:
                 print(ret.stderr.decode())
                 print(ret.stdout.decode())
                 raise RuntimeError(f'Marc simulation failed ({v})')

From d488f1708a476acb6aa8d7f8fc068eea4192254c Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sun, 11 Apr 2021 07:35:43 +0200
Subject: [PATCH 10/24] consistent naming

---
 src/homogenization.f90         |  2 +-
 src/homogenization_damage.f90  |  2 +-
 src/homogenization_thermal.f90 |  7 +++----
 src/lattice.f90                | 30 +++++++++++++++---------------
 4 files changed, 20 insertions(+), 21 deletions(-)

diff --git a/src/homogenization.f90 b/src/homogenization.f90
index 34e71a0f2..2408462c9 100644
--- a/src/homogenization.f90
+++ b/src/homogenization.f90
@@ -458,7 +458,7 @@ function homogenization_K_phi(ce)
 
   do co = 1, homogenization_Nconstituents(ho)
     homogenization_K_phi = homogenization_K_phi + &
-      crystallite_push33ToRef(co,ce,lattice_D(1:3,1:3,material_phaseID(co,ce)))
+      crystallite_push33ToRef(co,ce,lattice_K_phi(1:3,1:3,material_phaseID(co,ce)))
   enddo
 
   homogenization_K_phi = &
diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90
index f0f13bdbb..c7c82331b 100644
--- a/src/homogenization_damage.f90
+++ b/src/homogenization_damage.f90
@@ -114,7 +114,7 @@ module function homogenization_mu_phi(ce) result(mu)
   integer, intent(in) :: ce
   real(pReal) :: mu
 
-  mu = lattice_M(material_phaseID(1,ce))
+  mu = lattice_mu_phi(material_phaseID(1,ce))
 
 end function homogenization_mu_phi
 
diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90
index 8b4eafc3a..eb9f796ce 100644
--- a/src/homogenization_thermal.f90
+++ b/src/homogenization_thermal.f90
@@ -117,10 +117,9 @@ module function homogenization_K_T(ce) result(K)
   integer :: &
     co
 
-  K = 0.0_pReal
-
-  do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
-    K = K + crystallite_push33ToRef(co,ce,lattice_K(:,:,material_phaseID(co,ce)))
+  K = crystallite_push33ToRef(co,1,lattice_K_T(:,:,material_phaseID(1,ce)))
+  do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
+    K = K + crystallite_push33ToRef(co,ce,lattice_K_T(:,:,material_phaseID(co,ce)))
   enddo
 
   K = K / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
diff --git a/src/lattice.f90 b/src/lattice.f90
index 5ec365fa7..69ee65857 100644
--- a/src/lattice.f90
+++ b/src/lattice.f90
@@ -394,13 +394,13 @@ module lattice
 ! SHOULD NOT BE PART OF LATTICE BEGIN
   real(pReal),                        dimension(:),     allocatable, public, protected :: &
     lattice_mu, lattice_nu, &
-    lattice_M, &
+    lattice_mu_phi, &
     lattice_rho, &
     lattice_c_p
    real(pReal),                       dimension(:,:,:), allocatable, public, protected :: &
     lattice_C66, &
-    lattice_K, &
-    lattice_D
+    lattice_K_T, &
+    lattice_K_phi
  integer(kind(lattice_UNDEFINED_ID)), dimension(:),     allocatable, public, protected :: &
     lattice_structure
 ! SHOULD NOT BE PART OF LATTICE END
@@ -470,10 +470,10 @@ subroutine lattice_init
   allocate(lattice_structure(Nphases),source = lattice_UNDEFINED_ID)
   allocate(lattice_C66(6,6,Nphases),  source=0.0_pReal)
 
-  allocate(lattice_K  (3,3,Nphases), source=0.0_pReal)
-  allocate(lattice_D  (3,3,Nphases), source=0.0_pReal)
+  allocate(lattice_K_T  (3,3,Nphases), source=0.0_pReal)
+  allocate(lattice_K_phi  (3,3,Nphases), source=0.0_pReal)
 
-  allocate(lattice_M,&
+  allocate(lattice_mu_phi,&
            lattice_rho,lattice_c_p, &
            lattice_mu, lattice_nu,&
            source=[(0.0_pReal,i=1,Nphases)])
@@ -527,10 +527,10 @@ subroutine lattice_init
 
     if (phase%contains('thermal')) then
       thermal  => phase%get('thermal')
-      lattice_K(1,1,ph) = thermal%get_asFloat('K_11',defaultVal=0.0_pReal)
-      lattice_K(2,2,ph) = thermal%get_asFloat('K_22',defaultVal=0.0_pReal)
-      lattice_K(3,3,ph) = thermal%get_asFloat('K_33',defaultVal=0.0_pReal)
-      lattice_K(1:3,1:3,ph) = lattice_applyLatticeSymmetry33(lattice_K(1:3,1:3,ph), &
+      lattice_K_T(1,1,ph) = thermal%get_asFloat('K_11',defaultVal=0.0_pReal)
+      lattice_K_T(2,2,ph) = thermal%get_asFloat('K_22',defaultVal=0.0_pReal)
+      lattice_K_T(3,3,ph) = thermal%get_asFloat('K_33',defaultVal=0.0_pReal)
+      lattice_K_T(1:3,1:3,ph) = lattice_applyLatticeSymmetry33(lattice_K_T(1:3,1:3,ph), &
                                                              phase%get_asString('lattice'))
       lattice_c_p(ph) = thermal%get_asFloat('c_p', defaultVal=0.0_pReal)
     endif
@@ -539,13 +539,13 @@ subroutine lattice_init
     if (phase%contains('damage')) then
       damage => phase%get('damage')
       damage => damage%get(1)
-      lattice_D(1,1,ph) = damage%get_asFloat('D_11',defaultVal=0.0_pReal)
-      lattice_D(2,2,ph) = damage%get_asFloat('D_22',defaultVal=0.0_pReal)
-      lattice_D(3,3,ph) = damage%get_asFloat('D_33',defaultVal=0.0_pReal)
-      lattice_D(1:3,1:3,ph) = lattice_applyLatticeSymmetry33(lattice_D(1:3,1:3,ph), &
+      lattice_K_phi(1,1,ph) = damage%get_asFloat('D_11',defaultVal=0.0_pReal)
+      lattice_K_phi(2,2,ph) = damage%get_asFloat('D_22',defaultVal=0.0_pReal)
+      lattice_K_phi(3,3,ph) = damage%get_asFloat('D_33',defaultVal=0.0_pReal)
+      lattice_K_phi(1:3,1:3,ph) = lattice_applyLatticeSymmetry33(lattice_K_phi(1:3,1:3,ph), &
                                                           phase%get_asString('lattice'))
 
-      lattice_M(ph) = damage%get_asFloat('M',defaultVal=0.0_pReal)
+      lattice_mu_phi(ph) = damage%get_asFloat('M',defaultVal=0.0_pReal)
     endif
     ! SHOULD NOT BE PART OF LATTICE END
 

From a386b82f748c92430fd4f4292bc5802f29cc3639 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sun, 11 Apr 2021 07:48:54 +0200
Subject: [PATCH 11/24] distributing responsibility

---
 src/homogenization.f90         | 31 +++++--------------------------
 src/homogenization_damage.f90  | 15 +++++++++++++++
 src/homogenization_thermal.f90 |  1 -
 3 files changed, 20 insertions(+), 27 deletions(-)

diff --git a/src/homogenization.f90 b/src/homogenization.f90
index 2408462c9..195ede741 100644
--- a/src/homogenization.f90
+++ b/src/homogenization.f90
@@ -166,6 +166,11 @@ module homogenization
       real(pReal) :: mu
     end function homogenization_mu_phi
 
+    module function homogenization_K_phi(ce) result(K)
+      integer, intent(in) :: ce
+      real(pReal), dimension(3,3) :: K
+    end function homogenization_K_phi
+
     module function homogenization_f_phi(phi,ce) result(f)
       integer, intent(in) :: ce
       real(pReal), intent(in) :: phi
@@ -441,32 +446,6 @@ subroutine homogenization_restartRead(fileHandle)
 end subroutine homogenization_restartRead
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief returns homogenized non local damage diffusion tensor in reference configuration
-!--------------------------------------------------------------------------------------------------
-function homogenization_K_phi(ce)
-
-  integer, intent(in) :: ce
-  real(pReal), dimension(3,3) :: &
-    homogenization_K_phi
-  integer :: &
-    ho, &
-    co
-
-  ho  = material_homogenizationID(ce)
-  homogenization_K_phi = 0.0_pReal
-
-  do co = 1, homogenization_Nconstituents(ho)
-    homogenization_K_phi = homogenization_K_phi + &
-      crystallite_push33ToRef(co,ce,lattice_K_phi(1:3,1:3,material_phaseID(co,ce)))
-  enddo
-
-  homogenization_K_phi = &
-    num_damage%charLength**2*homogenization_K_phi/real(homogenization_Nconstituents(ho),pReal)
-
-end function homogenization_K_phi
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief parses the homogenization part from the material configuration
 !--------------------------------------------------------------------------------------------------
diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90
index c7c82331b..79906b8f0 100644
--- a/src/homogenization_damage.f90
+++ b/src/homogenization_damage.f90
@@ -154,6 +154,21 @@ module subroutine homogenization_set_phi(phi,ce)
 end subroutine homogenization_set_phi
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief returns homogenized non local damage diffusion tensor in reference configuration
+!--------------------------------------------------------------------------------------------------
+module function homogenization_K_phi(ce) result(K)
+
+  integer, intent(in) :: ce
+  real(pReal), dimension(3,3) :: K
+
+
+  K = crystallite_push33ToRef(1,ce,lattice_K_phi(1:3,1:3,material_phaseID(1,ce))) \
+    * num_damage%charLength**2
+
+end function homogenization_K_phi
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief writes results to HDF5 output file
 !--------------------------------------------------------------------------------------------------
diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90
index eb9f796ce..6e47d8854 100644
--- a/src/homogenization_thermal.f90
+++ b/src/homogenization_thermal.f90
@@ -229,7 +229,6 @@ module function homogenization_T(ce) result(T)
 end function homogenization_T
 
 
-
 !--------------------------------------------------------------------------------------------------
 !> @brief return heat generation rate
 !--------------------------------------------------------------------------------------------------

From 4b89e2f40cb6955d134bc389cc4926c10ab1f592 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sun, 11 Apr 2021 08:32:13 +0200
Subject: [PATCH 12/24] sorted and documented

---
 src/homogenization.f90         |  25 +++----
 src/homogenization_damage.f90  |  81 +++++++++++-----------
 src/homogenization_thermal.f90 | 122 +++++++++++++++------------------
 3 files changed, 108 insertions(+), 120 deletions(-)

diff --git a/src/homogenization.f90 b/src/homogenization.f90
index 195ede741..a1aaff1db 100644
--- a/src/homogenization.f90
+++ b/src/homogenization.f90
@@ -135,32 +135,26 @@ module homogenization
       logical, dimension(2) :: doneAndHappy
     end function mechanical_updateState
 
+    module function homogenization_mu_T(ce) result(mu)
+      integer, intent(in) :: ce
+      real(pReal) :: mu
+    end function homogenization_mu_T
 
     module function homogenization_K_T(ce) result(K)
       integer, intent(in) :: ce
       real(pReal), dimension(3,3) :: K
     end function homogenization_K_T
 
-    module function homogenization_mu_T(ce) result(mu)
+    module function homogenization_f_T(ce) result(f)
       integer, intent(in) :: ce
-      real(pReal) :: mu
-    end function homogenization_mu_T
+      real(pReal) :: f
+    end function homogenization_f_T
 
     module subroutine homogenization_thermal_setField(T,dot_T, ce)
       integer, intent(in) :: ce
       real(pReal),   intent(in) :: T, dot_T
     end subroutine homogenization_thermal_setField
 
-    module function homogenization_T(ce) result(T)
-      integer, intent(in) :: ce
-      real(pReal) :: T
-    end function homogenization_T
-
-    module function homogenization_f_T(ce) result(f)
-      integer, intent(in) :: ce
-      real(pReal) :: f
-    end function homogenization_f_T
-
     module function homogenization_mu_phi(ce) result(mu)
       integer, intent(in) :: ce
       real(pReal) :: mu
@@ -191,12 +185,11 @@ module homogenization
     homogenization_mu_T, &
     homogenization_K_T, &
     homogenization_f_T, &
-    homogenization_K_phi, &
+    homogenization_thermal_setfield, &
     homogenization_mu_phi, &
+    homogenization_K_phi, &
     homogenization_f_phi, &
     homogenization_set_phi, &
-    homogenization_thermal_setfield, &
-    homogenization_T, &
     homogenization_forward, &
     homogenization_results, &
     homogenization_restartRead, &
diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90
index 79906b8f0..8a4e596a8 100644
--- a/src/homogenization_damage.f90
+++ b/src/homogenization_damage.f90
@@ -26,8 +26,8 @@ submodule(homogenization) damage
   type(tparameters),             dimension(:), allocatable :: &
     param
 
-contains
 
+contains
 
 !--------------------------------------------------------------------------------------------------
 !> @brief Allocate variables and set parameters.
@@ -105,57 +105,22 @@ module subroutine damage_partition(ce)
 end subroutine damage_partition
 
 
-
 !--------------------------------------------------------------------------------------------------
-!> @brief Returns homogenized nonlocal damage mobility
+!> @brief Homogenized damage viscosity.
 !--------------------------------------------------------------------------------------------------
 module function homogenization_mu_phi(ce) result(mu)
 
   integer, intent(in) :: ce
   real(pReal) :: mu
 
+
   mu = lattice_mu_phi(material_phaseID(1,ce))
 
 end function homogenization_mu_phi
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief  calculates homogenized damage driving forces
-!--------------------------------------------------------------------------------------------------
-module function homogenization_f_phi(phi,ce) result(f)
-
-  integer, intent(in) :: ce
-  real(pReal), intent(in) :: &
-    phi
-  real(pReal) :: f
-
-  f = phase_f_phi(phi, 1, ce)
-
-end function homogenization_f_phi
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief updated nonlocal damage field with solution from damage phase field PDE
-!--------------------------------------------------------------------------------------------------
-module subroutine homogenization_set_phi(phi,ce)
-
-  integer, intent(in) :: ce
-  real(pReal),   intent(in) :: &
-    phi
-  integer :: &
-    ho, &
-    en
-
-  ho = material_homogenizationID(ce)
-  en = material_homogenizationEntry(ce)
-  damagestate_h(ho)%state(1,en) = phi
-  current(ho)%phi(en) = phi
-
-end subroutine homogenization_set_phi
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief returns homogenized non local damage diffusion tensor in reference configuration
+!> @brief Homogenized damage conductivity/diffusivity in reference configuration.
 !--------------------------------------------------------------------------------------------------
 module function homogenization_K_phi(ce) result(K)
 
@@ -169,6 +134,44 @@ module function homogenization_K_phi(ce) result(K)
 end function homogenization_K_phi
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief Homogenized damage driving force.
+!--------------------------------------------------------------------------------------------------
+module function homogenization_f_phi(phi,ce) result(f)
+
+  integer, intent(in) :: ce
+  real(pReal), intent(in) :: &
+    phi
+  real(pReal) :: f
+
+
+  f = phase_f_phi(phi, 1, ce)
+
+end function homogenization_f_phi
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Set damage field.
+!--------------------------------------------------------------------------------------------------
+module subroutine homogenization_set_phi(phi,ce)
+
+  integer, intent(in) :: ce
+  real(pReal),   intent(in) :: &
+    phi
+
+  integer :: &
+    ho, &
+    en
+
+
+  ho = material_homogenizationID(ce)
+  en = material_homogenizationEntry(ce)
+  damagestate_h(ho)%state(1,en) = phi
+  current(ho)%phi(en) = phi
+
+end subroutine homogenization_set_phi
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief writes results to HDF5 output file
 !--------------------------------------------------------------------------------------------------
diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90
index 6e47d8854..9c58fc6a8 100644
--- a/src/homogenization_thermal.f90
+++ b/src/homogenization_thermal.f90
@@ -107,15 +107,29 @@ end subroutine thermal_homogenize
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief return homogenized thermal conductivity in reference configuration
+!> @brief Homogenized thermal viscosity.
+!--------------------------------------------------------------------------------------------------
+module function homogenization_mu_T(ce) result(mu)
+
+  integer, intent(in) :: ce
+  real(pReal) :: mu
+
+
+  mu = c_P(ce) * rho(ce)
+
+end function homogenization_mu_T
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Homogenized thermal conductivity in reference configuration.
 !--------------------------------------------------------------------------------------------------
 module function homogenization_K_T(ce) result(K)
 
   integer, intent(in) :: ce
   real(pReal), dimension(3,3) :: K
 
-  integer :: &
-    co
+  integer :: co
+
 
   K = crystallite_push33ToRef(co,1,lattice_K_T(:,:,material_phaseID(1,ce)))
   do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
@@ -127,61 +141,29 @@ module function homogenization_K_T(ce) result(K)
 end function homogenization_K_T
 
 
-module function homogenization_mu_T(ce) result(mu)
+!--------------------------------------------------------------------------------------------------
+!> @brief Homogenized heat generation rate.
+!--------------------------------------------------------------------------------------------------
+module function homogenization_f_T(ce) result(f)
 
   integer, intent(in) :: ce
-  real(pReal) :: mu
-
-  mu = c_P(ce) * rho(ce)
-
-end function homogenization_mu_T
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief returns homogenized specific heat capacity
-!--------------------------------------------------------------------------------------------------
-function c_P(ce)
-
-  integer, intent(in) :: ce
-  real(pReal) :: c_P
+  real(pReal) :: f
 
   integer :: co
 
 
-  c_P = lattice_c_p(material_phaseID(1,ce))
+  f = phase_f_T(material_phaseID(1,ce),material_phaseEntry(1,ce))
   do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
-    c_P = c_P + lattice_c_p(material_phaseID(co,ce))
+    f = f + phase_f_T(material_phaseID(co,ce),material_phaseEntry(co,ce))
   enddo
 
-  c_P = c_P / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
+  f = f/real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
 
-end function c_P
+end function homogenization_f_T
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief returns homogenized mass density
-!--------------------------------------------------------------------------------------------------
-function rho(ce)
-
-  integer, intent(in) :: ce
-  real(pReal) :: rho
-
-  integer :: co
-
-
-  rho = lattice_rho(material_phaseID(1,ce))
-  do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
-    rho = rho + lattice_rho(material_phaseID(co,ce))
-  enddo
-
-  rho = rho / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
-
-end function rho
-
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief Set thermal field and its rate (T and dot_T)
+!> @brief Set thermal field and its rate (T and dot_T).
 !--------------------------------------------------------------------------------------------------
 module subroutine homogenization_thermal_setField(T,dot_T, ce)
 
@@ -196,7 +178,6 @@ module subroutine homogenization_thermal_setField(T,dot_T, ce)
 end subroutine homogenization_thermal_setField
 
 
-
 !--------------------------------------------------------------------------------------------------
 !> @brief writes results to HDF5 output file
 !--------------------------------------------------------------------------------------------------
@@ -219,34 +200,45 @@ module subroutine thermal_results(ho,group)
 end subroutine thermal_results
 
 
-module function homogenization_T(ce) result(T)
+!--------------------------------------------------------------------------------------------------
+!> @brief Homogenize specific heat capacity.
+!--------------------------------------------------------------------------------------------------
+function c_P(ce)
 
   integer, intent(in) :: ce
-  real(pReal) :: T
-
-  T = current(material_homogenizationID(ce))%T(material_homogenizationEntry(ce))
-
-end function homogenization_T
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief return heat generation rate
-!--------------------------------------------------------------------------------------------------
-module function homogenization_f_T(ce) result(f)
-
-  integer, intent(in) :: ce
-  real(pReal) :: f
+  real(pReal) :: c_P
 
   integer :: co
 
-  f = phase_f_T(material_phaseID(1,ce),material_phaseEntry(1,ce))
+
+  c_P = lattice_c_p(material_phaseID(1,ce))
   do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
-    f = f + phase_f_T(material_phaseID(co,ce),material_phaseEntry(co,ce))
+    c_P = c_P + lattice_c_p(material_phaseID(co,ce))
   enddo
 
-  f = f/real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
+  c_P = c_P / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
 
-end function homogenization_f_T
+end function c_P
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief Homogenize mass density.
+!--------------------------------------------------------------------------------------------------
+function rho(ce)
+
+  integer, intent(in) :: ce
+  real(pReal) :: rho
+
+  integer :: co
+
+
+  rho = lattice_rho(material_phaseID(1,ce))
+  do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
+    rho = rho + lattice_rho(material_phaseID(co,ce))
+  enddo
+
+  rho = rho / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
+
+end function rho
+
 end submodule thermal

From 34bb4c65a9b5b95732969bbb287c1a900ab4c546 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sun, 11 Apr 2021 08:46:31 +0200
Subject: [PATCH 13/24] distributing responsibility

---
 src/phase_damage.f90  | 14 +++++++++++++-
 src/phase_thermal.f90 | 11 +++++++++++
 2 files changed, 24 insertions(+), 1 deletion(-)

diff --git a/src/phase_damage.f90 b/src/phase_damage.f90
index 189a5eabf..c0c0768e7 100644
--- a/src/phase_damage.f90
+++ b/src/phase_damage.f90
@@ -2,6 +2,12 @@
 !> @brief internal microstructure state for all damage sources and kinematics constitutive models
 !----------------------------------------------------------------------------------------------------
 submodule(phase) damage
+
+  type :: tDamageParameters
+    real(pReal) ::                 mu = 0.0_pReal                                                   !< viscosity
+    real(pReal), dimension(3,3) :: K  = 0.0_pReal                                                   !< conductivity/diffusivity
+  end type tDamageParameters
+
   enum, bind(c); enumerator :: &
     DAMAGE_UNDEFINED_ID, &
     DAMAGE_ISOBRITTLE_ID, &
@@ -16,13 +22,15 @@ submodule(phase) damage
   end type tDataContainer
 
   integer(kind(DAMAGE_UNDEFINED_ID)),     dimension(:), allocatable :: &
-    phase_source                                                                                !< active sources mechanisms of each phase
+    phase_source                                                                                    !< active sources mechanisms of each phase
 
   integer, dimension(:), allocatable :: &
     phase_Nsources
 
   type(tDataContainer), dimension(:), allocatable :: current
 
+  type(tDamageParameters), dimension(:), allocatable :: param
+
   interface
 
     module function anisobrittle_init() result(mySources)
@@ -111,6 +119,7 @@ module subroutine damage_init
 
   allocate(damageState (phases%length))
   allocate(phase_Nsources(phases%length),source = 0)
+  allocate(param(phases%length))
 
   do ph = 1,phases%length
 
@@ -121,6 +130,9 @@ module subroutine damage_init
 
     phase => phases%get(ph)
     sources => phase%get('damage',defaultVal=emptyList)
+    param(ph)%K(1,1) = sources%get_asFloat('K_11',defaultVal=0.0_pReal)
+    param(ph)%K(2,2) = sources%get_asFloat('K_22',defaultVal=0.0_pReal)
+    param(ph)%K(3,3) = sources%get_asFloat('K_33',defaultVal=0.0_pReal)
     if (sources%length > 1) error stop
     phase_Nsources(ph) = sources%length
 
diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90
index 807e1f655..c2b4b34b2 100644
--- a/src/phase_thermal.f90
+++ b/src/phase_thermal.f90
@@ -3,6 +3,11 @@
 !----------------------------------------------------------------------------------------------------
 submodule(phase) thermal
 
+  type :: tThermalParameters
+    real(pReal) ::                 C_p = 0.0_pReal                                                  !< heat capacity
+    real(pReal), dimension(3,3) :: K   = 0.0_pReal                                                  !< thermal conductivity
+  end type tThermalParameters
+
   integer, dimension(:), allocatable :: &
     thermal_Nsources
 
@@ -23,6 +28,8 @@ submodule(phase) thermal
 
   type(tDataContainer), dimension(:), allocatable :: current          ! ?? not very telling name. Better: "field" ?? MD: current(ho)%T(me) reads quite good
 
+  type(tThermalParameters), dimension(:), allocatable :: param
+
   integer :: thermal_source_maxSizeDotState
 
 
@@ -85,6 +92,7 @@ module subroutine thermal_init(phases)
 
   allocate(thermalState(phases%length))
   allocate(thermal_Nsources(phases%length),source = 0)
+  allocate(param(phases%length))
 
   do ph = 1, phases%length
     Nmembers = count(material_phaseID == ph)
@@ -93,6 +101,9 @@ module subroutine thermal_init(phases)
     phase => phases%get(ph)
     thermal => phase%get('thermal',defaultVal=emptyDict)
     sources => thermal%get('source',defaultVal=emptyList)
+    param(ph)%K(1,1) = thermal%get_asFloat('K_11',defaultVal=0.0_pReal)
+    param(ph)%K(2,2) = thermal%get_asFloat('K_22',defaultVal=0.0_pReal)
+    param(ph)%K(3,3) = thermal%get_asFloat('K_33',defaultVal=0.0_pReal)
     thermal_Nsources(ph) = sources%length
     allocate(thermalstate(ph)%p(thermal_Nsources(ph)))
   enddo

From b55e721ec486e34d34de493e6625cd638ae410c9 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sun, 11 Apr 2021 09:25:45 +0200
Subject: [PATCH 14/24] standard name

---
 src/CPFEM.f90 | 19 +++++++++----------
 1 file changed, 9 insertions(+), 10 deletions(-)

diff --git a/src/CPFEM.f90 b/src/CPFEM.f90
index 60dc78d2e..d8f4a0c3c 100644
--- a/src/CPFEM.f90
+++ b/src/CPFEM.f90
@@ -150,15 +150,14 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
                                                       H
 
   integer(pInt)                                       elCP, &                                       ! crystal plasticity element number
-                                                      i, j, k, l, m, n, ph, homog, mySource,ma
+                                                      i, j, k, l, m, n, ph, homog, mySource,ce
 
   real(pReal), parameter ::                          ODD_STRESS    = 1e15_pReal, &                  !< return value for stress if terminallyIll
                                                      ODD_JACOBIAN  = 1e50_pReal                     !< return value for jacobian if terminallyIll
 
 
   elCP = discretization_Marc_FEM2DAMASK_elem(elFE)
-
-  ma = (elCP-1) * discretization_nIPs + ip
+  ce   = discretization_Marc_FEM2DAMASK_cell(ip,elFE)
 
   if (debugCPFEM%basic .and. elCP == debugCPFEM%element .and. ip == debugCPFEM%ip) then
     print'(/,a)', '#############################################'
@@ -183,8 +182,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
     !    temperature(material_homogenizationAt(elCP))%p(material_homogenizationMemberAt(ip,elCP)) = &
     !      temperature_inp
     !end select chosenThermal1
-    homogenization_F0(1:3,1:3,ma) = ffn
-    homogenization_F(1:3,1:3,ma) = ffn1
+    homogenization_F0(1:3,1:3,ce) = ffn
+    homogenization_F(1:3,1:3,ce) = ffn1
 
   if (iand(mode, CPFEM_CALCRESULTS) /= 0_pInt) then
 
@@ -209,17 +208,17 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
       else terminalIllness
 
         ! translate from P to sigma
-        Kirchhoff = matmul(homogenization_P(1:3,1:3,ma), transpose(homogenization_F(1:3,1:3,ma)))
-        J_inverse  = 1.0_pReal / math_det33(homogenization_F(1:3,1:3,ma))
+        Kirchhoff = matmul(homogenization_P(1:3,1:3,ce), transpose(homogenization_F(1:3,1:3,ce)))
+        J_inverse  = 1.0_pReal / math_det33(homogenization_F(1:3,1:3,ce))
         CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.)
 
         !  translate from dP/dF to dCS/dE
         H = 0.0_pReal
         do i=1,3; do j=1,3; do k=1,3; do l=1,3; do m=1,3; do n=1,3
           H(i,j,k,l) = H(i,j,k,l) &
-                     +  homogenization_F(j,m,ma) * homogenization_F(l,n,ma) &
-                                                 * homogenization_dPdF(i,m,k,n,ma) &
-                     -  math_delta(j,l) * homogenization_F(i,m,ma) * homogenization_P(k,m,ma) &
+                     +  homogenization_F(j,m,ce) * homogenization_F(l,n,ce) &
+                                                 * homogenization_dPdF(i,m,k,n,ce) &
+                     -  math_delta(j,l) * homogenization_F(i,m,ce) * homogenization_P(k,m,ce) &
                      +  0.5_pReal * (  Kirchhoff(j,l)*math_delta(i,k) + Kirchhoff(i,k)*math_delta(j,l) &
                                      + Kirchhoff(j,k)*math_delta(i,l) + Kirchhoff(i,l)*math_delta(j,k))
         enddo; enddo; enddo; enddo; enddo; enddo

From 547f2ffa69fcb2b2548968d37189f74de1f113c6 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sun, 11 Apr 2021 09:46:11 +0200
Subject: [PATCH 15/24] cleaning

---
 src/DAMASK_Marc.f90                |  3 +--
 src/grid/discretization_grid.f90   |  2 +-
 src/grid/grid_damage_spectral.f90  | 14 ++++++--------
 src/grid/grid_thermal_spectral.f90 | 10 +++++-----
 4 files changed, 13 insertions(+), 16 deletions(-)

diff --git a/src/DAMASK_Marc.f90 b/src/DAMASK_Marc.f90
index 43ec9b084..910ca86c0 100644
--- a/src/DAMASK_Marc.f90
+++ b/src/DAMASK_Marc.f90
@@ -218,7 +218,7 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
   logical :: cutBack
   real(pReal), dimension(6) ::   stress
   real(pReal), dimension(6,6) :: ddsdde
-  integer :: computationMode, i, cp_en, node, CPnodeID
+  integer :: computationMode, i, node, CPnodeID
   integer(pI32) :: defaultNumThreadsInt                                                             !< default value set by Marc
 
   integer(pInt), save :: &
@@ -266,7 +266,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
       computationMode = CPFEM_RESTOREJACOBIAN
   elseif (lovl == 6) then                                                                           ! stress requested by marc
     computationMode = CPFEM_CALCRESULTS
-    cp_en = discretization_Marc_FEM2DAMASK_elem(m(1))
     if (cptim > theTime .or. inc /= theInc) then                                                    ! reached "convergence"
       terminallyIll = .false.
       cycleCounter = -1                                                                             ! first calc step increments this to cycle = 0
diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90
index 88216f0aa..14391fd23 100644
--- a/src/grid/discretization_grid.f90
+++ b/src/grid/discretization_grid.f90
@@ -74,7 +74,7 @@ subroutine discretization_grid_init(restart)
     allocate(materialAt_global(0))                                                                  ! needed for IntelMPI
   endif
 
-  if (grid(1) < 2) call IO_error(844, ext_msg='cells(1) must be larger than  1')
+  if (grid(1) < 2) call IO_error(844, ext_msg='cells(1) must be larger than 1')
 
   call MPI_Bcast(grid,3,MPI_INTEGER,0,PETSC_COMM_WORLD, ierr)
   if (ierr /= 0) error stop 'MPI error'
diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90
index 7e0d4112a..7fe9a4a25 100644
--- a/src/grid/grid_damage_spectral.f90
+++ b/src/grid/grid_damage_spectral.f90
@@ -259,7 +259,6 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
   PetscObject :: dummy
   PetscErrorCode :: ierr
   integer :: i, j, k, ce
-  real(pReal) :: mobility
 
   phi_current = x_scal
 !--------------------------------------------------------------------------------------------------
@@ -281,9 +280,8 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
   ce = 0
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
     ce = ce + 1
-    mobility = homogenization_mu_phi(ce)
     scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + homogenization_f_phi(phi_current(i,j,k),ce)) &
-                            + mobility*(phi_lastInc(i,j,k) - phi_current(i,j,k)) &
+                            + homogenization_mu_phi(ce)*(phi_lastInc(i,j,k) - phi_current(i,j,k)) &
                             + mu_ref*phi_current(i,j,k)
   enddo; enddo; enddo
 
@@ -309,16 +307,16 @@ end subroutine formResidual
 !--------------------------------------------------------------------------------------------------
 subroutine updateReference
 
-  integer :: i,j,k,ce,ierr
+  integer :: ce,ierr
+
 
-  ce = 0
   K_ref = 0.0_pReal
   mu_ref = 0.0_pReal
-  do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
-    ce = ce + 1
+  do ce = 1, product(grid(1:2))*grid3
     K_ref  = K_ref  + homogenization_K_phi(ce)
     mu_ref = mu_ref + homogenization_mu_phi(ce)
-  enddo; enddo; enddo
+  enddo
+
   K_ref = K_ref*wgt
   call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
   mu_ref = mu_ref*wgt
diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90
index ea6caec0e..60d8d99db 100644
--- a/src/grid/grid_thermal_spectral.f90
+++ b/src/grid/grid_thermal_spectral.f90
@@ -302,16 +302,16 @@ end subroutine formResidual
 !--------------------------------------------------------------------------------------------------
 subroutine updateReference
 
-  integer :: i,j,k,ce,ierr
+  integer :: ce,ierr
+
 
-  ce = 0
   K_ref = 0.0_pReal
   mu_ref = 0.0_pReal
-  do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
-    ce = ce + 1
+  do ce = 1, product(grid(1:2))*grid3
     K_ref  = K_ref  + homogenization_K_T(ce)
     mu_ref = mu_ref + homogenization_mu_T(ce)
-  enddo; enddo; enddo
+  enddo
+
   K_ref = K_ref*wgt
   call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
   mu_ref = mu_ref*wgt

From b67e8375483a351d548ad9c9f4410cad6fc14981 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sun, 11 Apr 2021 09:58:40 +0200
Subject: [PATCH 16/24] needs broadcast, otherwise rank>0 fail

---
 src/grid/discretization_grid.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90
index 14391fd23..ad614d8c5 100644
--- a/src/grid/discretization_grid.f90
+++ b/src/grid/discretization_grid.f90
@@ -74,10 +74,10 @@ subroutine discretization_grid_init(restart)
     allocate(materialAt_global(0))                                                                  ! needed for IntelMPI
   endif
 
-  if (grid(1) < 2) call IO_error(844, ext_msg='cells(1) must be larger than 1')
 
   call MPI_Bcast(grid,3,MPI_INTEGER,0,PETSC_COMM_WORLD, ierr)
   if (ierr /= 0) error stop 'MPI error'
+  if (grid(1) < 2) call IO_error(844, ext_msg='cells(1) must be larger than 1')
   call MPI_Bcast(geomSize,3,MPI_DOUBLE,0,PETSC_COMM_WORLD, ierr)
   if (ierr /= 0) error stop 'MPI error'
   call MPI_Bcast(origin,3,MPI_DOUBLE,0,PETSC_COMM_WORLD, ierr)

From 8b7f77718627c88789802838a0690fc5a0d333e8 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sun, 11 Apr 2021 11:17:52 +0200
Subject: [PATCH 17/24] probably not needed

---
 src/homogenization.f90         |  5 -----
 src/homogenization_thermal.f90 | 12 ------------
 2 files changed, 17 deletions(-)

diff --git a/src/homogenization.f90 b/src/homogenization.f90
index a1aaff1db..3fadf1e96 100644
--- a/src/homogenization.f90
+++ b/src/homogenization.f90
@@ -100,10 +100,6 @@ module homogenization
       integer,     intent(in) :: ce
     end subroutine damage_partition
 
-    module subroutine thermal_homogenize(ip,el)
-      integer, intent(in) :: ip,el
-    end subroutine thermal_homogenize
-
     module subroutine mechanical_homogenize(dt,ce)
      real(pReal), intent(in) :: dt
      integer, intent(in) :: &
@@ -306,7 +302,6 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
             terminallyIll = .true.                                                                  ! ...and kills all others
           endif
         enddo
-        call thermal_homogenize(ip,el)
       enddo
     enddo
     !$OMP END DO
diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90
index 9c58fc6a8..40c466ab5 100644
--- a/src/homogenization_thermal.f90
+++ b/src/homogenization_thermal.f90
@@ -94,18 +94,6 @@ module subroutine thermal_partition(ce)
 end subroutine thermal_partition
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief Homogenize temperature rates
-!--------------------------------------------------------------------------------------------------
-module subroutine thermal_homogenize(ip,el)
-
-  integer, intent(in) :: ip,el
-
-  !call phase_thermal_getRate(homogenization_dot_T((el-1)*discretization_nIPs+ip), ip,el)
-
-end subroutine thermal_homogenize
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief Homogenized thermal viscosity.
 !--------------------------------------------------------------------------------------------------

From 071c1a5ad453f21d2ba56716147f25a3fa9ea81f Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sun, 11 Apr 2021 12:18:26 +0200
Subject: [PATCH 18/24] avoid repetition

---
 src/homogenization_mechanical.f90           | 61 +++------------
 src/homogenization_mechanical_RGC.f90       | 25 +------
 src/homogenization_mechanical_isostrain.f90 | 83 +++------------------
 src/homogenization_mechanical_pass.f90      | 21 +++---
 4 files changed, 36 insertions(+), 154 deletions(-)

diff --git a/src/homogenization_mechanical.f90 b/src/homogenization_mechanical.f90
index 7d1c64445..7ea1f74e9 100644
--- a/src/homogenization_mechanical.f90
+++ b/src/homogenization_mechanical.f90
@@ -32,25 +32,6 @@ submodule(homogenization) mechanical
     end subroutine RGC_partitionDeformation
 
 
-    module subroutine isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho)
-      real(pReal),   dimension (3,3),       intent(out) :: avgP                                     !< average stress at material point
-      real(pReal),   dimension (3,3,3,3),   intent(out) :: dAvgPdAvgF                               !< average stiffness at material point
-
-      real(pReal),   dimension (:,:,:),     intent(in)  :: P                                        !< partitioned stresses
-      real(pReal),   dimension (:,:,:,:,:), intent(in)  :: dPdF                                     !< partitioned stiffnesses
-      integer,                              intent(in)  :: ho
-    end subroutine isostrain_averageStressAndItsTangent
-
-    module subroutine RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho)
-      real(pReal),   dimension (3,3),       intent(out) :: avgP                                     !< average stress at material point
-      real(pReal),   dimension (3,3,3,3),   intent(out) :: dAvgPdAvgF                               !< average stiffness at material point
-
-      real(pReal),   dimension (:,:,:),     intent(in)  :: P                                        !< partitioned stresses
-      real(pReal),   dimension (:,:,:,:,:), intent(in)  :: dPdF                                     !< partitioned stiffnesses
-      integer,                              intent(in)  :: ho
-    end subroutine RGC_averageStressAndItsTangent
-
-
     module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
       logical, dimension(2) :: doneAndHappy
       real(pReal), dimension(:,:,:),     intent(in)    :: &
@@ -148,39 +129,21 @@ module subroutine mechanical_homogenize(dt,ce)
   integer, intent(in) :: ce
 
   integer :: co
-  real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationID(ce)))
-  real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationID(ce)))
 
 
-  chosenHomogenization: select case(homogenization_type(material_homogenizationID(ce)))
+  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)
+  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)
+  enddo
 
-    case (HOMOGENIZATION_NONE_ID) chosenHomogenization
-        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)
-
-    case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
-      do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
-        dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ce)
-        Ps(:,:,co) = phase_P(co,ce)
-      enddo
-      call isostrain_averageStressAndItsTangent(&
-        homogenization_P(1:3,1:3,ce), &
-        homogenization_dPdF(1:3,1:3,1:3,1:3,ce),&
-        Ps,dPdFs, &
-        material_homogenizationID(ce))
-
-    case (HOMOGENIZATION_RGC_ID) chosenHomogenization
-      do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
-        dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ce)
-        Ps(:,:,co) = phase_P(co,ce)
-      enddo
-      call RGC_averageStressAndItsTangent(&
-        homogenization_P(1:3,1:3,ce), &
-        homogenization_dPdF(1:3,1:3,1:3,1:3,ce),&
-        Ps,dPdFs, &
-        material_homogenizationID(ce))
-
-  end select chosenHomogenization
+  homogenization_P(1:3,1:3,ce)            = homogenization_P(1:3,1:3,ce) &
+                                          / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
+  homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = homogenization_dPdF(1:3,1:3,1:3,1:3,ce) &
+                                          / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
 
 end subroutine mechanical_homogenize
 
diff --git a/src/homogenization_mechanical_RGC.f90 b/src/homogenization_mechanical_RGC.f90
index 772d1c4f5..745c266d4 100644
--- a/src/homogenization_mechanical_RGC.f90
+++ b/src/homogenization_mechanical_RGC.f90
@@ -89,7 +89,8 @@ module subroutine RGC_init(num_homogMech)
 
   print'(/,a)', ' <<<+-  homogenization:mechanical:RGC init  -+>>>'
 
-  print'(a,i2)', ' # instances: ',count(homogenization_type == HOMOGENIZATION_RGC_ID); flush(IO_STDOUT)
+  print'(a,i2)', ' # instances: ',count(homogenization_type == HOMOGENIZATION_RGC_ID)
+  flush(IO_STDOUT)
 
   print*, 'D.D. Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009'
   print*, 'https://doi.org/10.1007/s12289-009-0619-1'//IO_EOL
@@ -207,7 +208,7 @@ module subroutine RGC_partitionDeformation(F,avgF,ce)
   integer ::  iGrain,iFace,i,j,ho,en
 
   associate(prm => param(material_homogenizationID(ce)))
- 
+
   ho = material_homogenizationID(ce)
   en = material_homogenizationEntry(ce)
 !--------------------------------------------------------------------------------------------------
@@ -700,24 +701,6 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy)
 end function RGC_updateState
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief derive average stress and stiffness from constituent quantities
-!--------------------------------------------------------------------------------------------------
-module subroutine RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho)
-
-  real(pReal), dimension (3,3),        intent(out) :: avgP                                          !< average stress at material point
-  real(pReal), dimension (3,3,3,3),    intent(out) :: dAvgPdAvgF                                    !< average stiffness at material point
-
-  real(pReal), dimension (:,:,:),      intent(in)  :: P                                             !< partitioned stresses
-  real(pReal), dimension (:,:,:,:,:),  intent(in)  :: dPdF                                          !< partitioned stiffnesses
-  integer,                             intent(in)  :: ho
-
-  avgP       = sum(P,3)   /real(product(param(ho)%N_constituents),pReal)
-  dAvgPdAvgF = sum(dPdF,5)/real(product(param(ho)%N_constituents),pReal)
-
-end subroutine RGC_averageStressAndItsTangent
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief writes results to HDF5 output file
 !--------------------------------------------------------------------------------------------------
@@ -802,7 +785,7 @@ pure function interfaceNormal(intFace,ho,en)
   interfaceNormal(nPos) = real(intFace(1)/abs(intFace(1)),pReal)                                    ! get the normal vector w.r.t. cluster axis
 
   interfaceNormal = matmul(dst%orientation(1:3,1:3,en),interfaceNormal)                             ! map the normal vector into sample coordinate system (basis)
-  
+
   end associate
 
 end function interfaceNormal
diff --git a/src/homogenization_mechanical_isostrain.f90 b/src/homogenization_mechanical_isostrain.f90
index 9a3704575..7b114d04f 100644
--- a/src/homogenization_mechanical_isostrain.f90
+++ b/src/homogenization_mechanical_isostrain.f90
@@ -6,21 +6,6 @@
 !--------------------------------------------------------------------------------------------------
 submodule(homogenization:mechanical) isostrain
 
-  enum, bind(c); enumerator :: &
-    parallel_ID, &
-    average_ID
-  end enum
-
-  type :: tParameters                                                                               !< container type for internal constitutive parameters
-    integer :: &
-      N_constituents
-    integer(kind(average_ID)) :: &
-      mapping
-  end type
-
-  type(tParameters), dimension(:), allocatable :: param                                             !< containers of constitutive parameters (len Ninstances)
-
-
 contains
 
 !--------------------------------------------------------------------------------------------------
@@ -29,42 +14,21 @@ contains
 module subroutine isostrain_init
 
   integer :: &
-    h, &
+    ho, &
     Nmaterialpoints
-  class(tNode), pointer :: &
-    material_homogenization, &
-    homog, &
-    homogMech
 
   print'(/,a)', ' <<<+-  homogenization:mechanical:isostrain init  -+>>>'
 
-  print'(a,i2)', ' # instances: ',count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID); flush(IO_STDOUT)
+  print'(a,i2)', ' # instances: ',count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)
+  flush(IO_STDOUT)
 
-  material_homogenization => config_material%get('homogenization')
-  allocate(param(material_homogenization%length))                                                   ! one container of parameters per homog
+  do ho = 1, size(homogenization_type)
+    if (homogenization_type(ho) /= HOMOGENIZATION_ISOSTRAIN_ID) cycle
 
-  do h = 1, size(homogenization_type)
-    if (homogenization_type(h) /= HOMOGENIZATION_ISOSTRAIN_ID) cycle
-    homog => material_homogenization%get(h)
-    homogMech => homog%get('mechanical')
-    associate(prm => param(h))
-
-    prm%N_constituents = homogenization_Nconstituents(h)
-    select case(homogMech%get_asString('mapping',defaultVal = 'sum'))
-      case ('sum')
-        prm%mapping = parallel_ID
-      case ('avg')
-        prm%mapping = average_ID
-      case default
-        call IO_error(211,ext_msg='sum'//' (isostrain)')
-    end select
-
-    Nmaterialpoints = count(material_homogenizationAt == h)
-    homogState(h)%sizeState       = 0
-    allocate(homogState(h)%state0   (0,Nmaterialpoints))
-    allocate(homogState(h)%state    (0,Nmaterialpoints))
-
-    end associate
+    Nmaterialpoints = count(material_homogenizationAt == ho)
+    homogState(ho)%sizeState = 0
+    allocate(homogState(ho)%state0(0,Nmaterialpoints))
+    allocate(homogState(ho)%state (0,Nmaterialpoints))
 
   enddo
 
@@ -80,36 +44,9 @@ module subroutine isostrain_partitionDeformation(F,avgF)
 
   real(pReal),   dimension (3,3),   intent(in)  :: avgF                                             !< average deformation gradient at material point
 
+
   F = spread(avgF,3,size(F,3))
 
 end subroutine isostrain_partitionDeformation
 
-
-!--------------------------------------------------------------------------------------------------
-!> @brief derive average stress and stiffness from constituent quantities
-!--------------------------------------------------------------------------------------------------
-module subroutine isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho)
-
-  real(pReal),   dimension (3,3),       intent(out) :: avgP                                         !< average stress at material point
-  real(pReal),   dimension (3,3,3,3),   intent(out) :: dAvgPdAvgF                                   !< average stiffness at material point
-
-  real(pReal),   dimension (:,:,:),     intent(in)  :: P                                            !< partitioned stresses
-  real(pReal),   dimension (:,:,:,:,:), intent(in)  :: dPdF                                         !< partitioned stiffnesses
-  integer,                              intent(in)  :: ho
-
-  associate(prm => param(ho))
-
-  select case (prm%mapping)
-    case (parallel_ID)
-      avgP       = sum(P,3)
-      dAvgPdAvgF = sum(dPdF,5)
-    case (average_ID)
-      avgP       = sum(P,3)   /real(prm%N_constituents,pReal)
-      dAvgPdAvgF = sum(dPdF,5)/real(prm%N_constituents,pReal)
-  end select
-
-  end associate
-
-end subroutine isostrain_averageStressAndItsTangent
-
 end submodule isostrain
diff --git a/src/homogenization_mechanical_pass.f90 b/src/homogenization_mechanical_pass.f90
index 23fe74f44..e2e44658a 100644
--- a/src/homogenization_mechanical_pass.f90
+++ b/src/homogenization_mechanical_pass.f90
@@ -14,25 +14,24 @@ contains
 module subroutine pass_init
 
   integer :: &
-    Ninstances, &
-    h, &
+    ho, &
     Nmaterialpoints
 
   print'(/,a)', ' <<<+-  homogenization:mechanical:pass init  -+>>>'
 
-  Ninstances = count(homogenization_type == HOMOGENIZATION_NONE_ID)
-  print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
+  print'(a,i2)', ' # instances: ',count(homogenization_type == HOMOGENIZATION_NONE_ID)
+  flush(IO_STDOUT)
 
-  do h = 1, size(homogenization_type)
-    if(homogenization_type(h) /= HOMOGENIZATION_NONE_ID) cycle
+  do ho = 1, size(homogenization_type)
+    if(homogenization_type(ho) /= HOMOGENIZATION_NONE_ID) cycle
 
-    if(homogenization_Nconstituents(h) /= 1) &
+    if(homogenization_Nconstituents(ho) /= 1) &
       call IO_error(211,ext_msg='N_constituents (pass)')
 
-    Nmaterialpoints = count(material_homogenizationAt == h)
-    homogState(h)%sizeState = 0
-    allocate(homogState(h)%state0   (0,Nmaterialpoints))
-    allocate(homogState(h)%state    (0,Nmaterialpoints))
+    Nmaterialpoints = count(material_homogenizationAt == ho)
+    homogState(ho)%sizeState = 0
+    allocate(homogState(ho)%state0(0,Nmaterialpoints))
+    allocate(homogState(ho)%state (0,Nmaterialpoints))
 
   enddo
 

From 081d51f224d64e65a8628d7d375ac29529746f6c Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sun, 11 Apr 2021 12:37:46 +0200
Subject: [PATCH 19/24] multiple damage mechanisms currently not supported

---
 src/phase_damage.f90 | 59 +++++++++++++++++++-------------------------
 1 file changed, 25 insertions(+), 34 deletions(-)

diff --git a/src/phase_damage.f90 b/src/phase_damage.f90
index c0c0768e7..51a3ddac5 100644
--- a/src/phase_damage.f90
+++ b/src/phase_damage.f90
@@ -24,9 +24,6 @@ submodule(phase) damage
   integer(kind(DAMAGE_UNDEFINED_ID)),     dimension(:), allocatable :: &
     phase_source                                                                                    !< active sources mechanisms of each phase
 
-  integer, dimension(:), allocatable :: &
-    phase_Nsources
-
   type(tDataContainer), dimension(:), allocatable :: current
 
   type(tDamageParameters), dimension(:), allocatable :: param
@@ -109,7 +106,7 @@ module subroutine damage_init
    phases, &
    phase, &
    sources
-
+  logical:: damage_active
 
   print'(/,a)', ' <<<+-  phase:damage init  -+>>>'
 
@@ -118,9 +115,9 @@ module subroutine damage_init
   allocate(current(phases%length))
 
   allocate(damageState (phases%length))
-  allocate(phase_Nsources(phases%length),source = 0)
   allocate(param(phases%length))
 
+  damage_active = .false.
   do ph = 1,phases%length
 
     Nmembers = count(material_phaseID == ph)
@@ -134,14 +131,13 @@ module subroutine damage_init
     param(ph)%K(2,2) = sources%get_asFloat('K_22',defaultVal=0.0_pReal)
     param(ph)%K(3,3) = sources%get_asFloat('K_33',defaultVal=0.0_pReal)
     if (sources%length > 1) error stop
-    phase_Nsources(ph) = sources%length
+    if (sources%length == 1) damage_active = .true.
 
   enddo
 
   allocate(phase_source(phases%length), source = DAMAGE_UNDEFINED_ID)
 
-! initialize source mechanisms
-  if(maxval(phase_Nsources) /= 0) then
+  if (damage_active) then
     where(isobrittle_init()  ) phase_source = DAMAGE_ISOBRITTLE_ID
     where(isoductile_init()  ) phase_source = DAMAGE_ISODUCTILE_ID
     where(anisobrittle_init()) phase_source = DAMAGE_ANISOBRITTLE_ID
@@ -288,30 +284,25 @@ module subroutine damage_results(group,ph)
   character(len=*), intent(in) :: group
   integer,          intent(in) :: ph
 
-  integer :: so
-
-  sourceLoop: do so = 1, phase_Nsources(ph)
 
   if (phase_source(ph) /= DAMAGE_UNDEFINED_ID) &
     call results_closeGroup(results_addGroup(group//'damage'))
 
-    sourceType: select case (phase_source(ph))
+  sourceType: select case (phase_source(ph))
 
-      case (DAMAGE_ISOBRITTLE_ID) sourceType
-        call isobrittle_results(ph,group//'damage/')
+    case (DAMAGE_ISOBRITTLE_ID) sourceType
+      call isobrittle_results(ph,group//'damage/')
 
-      case (DAMAGE_ISODUCTILE_ID) sourceType
-        call isoductile_results(ph,group//'damage/')
+    case (DAMAGE_ISODUCTILE_ID) sourceType
+      call isoductile_results(ph,group//'damage/')
 
-      case (DAMAGE_ANISOBRITTLE_ID) sourceType
-        call anisobrittle_results(ph,group//'damage/')
+    case (DAMAGE_ANISOBRITTLE_ID) sourceType
+      call anisobrittle_results(ph,group//'damage/')
 
-      case (DAMAGE_ANISODUCTILE_ID) sourceType
-        call anisoductile_results(ph,group//'damage/')
+    case (DAMAGE_ANISODUCTILE_ID) sourceType
+      call anisoductile_results(ph,group//'damage/')
 
-    end select sourceType
-
-  enddo SourceLoop
+  end select sourceType
 
 end subroutine damage_results
 
@@ -374,19 +365,19 @@ function phase_damage_deltaState(Fe, ph, me) result(broken)
 
   if (damageState(ph)%sizeState == 0) return
 
-     sourceType: select case (phase_source(ph))
+   sourceType: select case (phase_source(ph))
 
-      case (DAMAGE_ISOBRITTLE_ID) sourceType
-        call isobrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me)
-        broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,me)))
-        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,me) + damageState(ph)%deltaState(1:mySize,me)
-        endif
+    case (DAMAGE_ISOBRITTLE_ID) sourceType
+      call isobrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me)
+      broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,me)))
+      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,me) + damageState(ph)%deltaState(1:mySize,me)
+      endif
 
-    end select sourceType
+  end select sourceType
 
 
 end function phase_damage_deltaState

From 9d90ed525e442f4a843db05cf22be33f30ae6ea9 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sun, 11 Apr 2021 12:42:57 +0200
Subject: [PATCH 20/24] YAML structure allows more than one mechanism

---
 src/phase_damage.f90 | 15 +++++++++------
 1 file changed, 9 insertions(+), 6 deletions(-)

diff --git a/src/phase_damage.f90 b/src/phase_damage.f90
index 51a3ddac5..5f6d65ab8 100644
--- a/src/phase_damage.f90
+++ b/src/phase_damage.f90
@@ -105,7 +105,8 @@ module subroutine damage_init
   class(tNode), pointer :: &
    phases, &
    phase, &
-   sources
+   sources, &
+   source
   logical:: damage_active
 
   print'(/,a)', ' <<<+-  phase:damage init  -+>>>'
@@ -113,7 +114,6 @@ module subroutine damage_init
   phases => config_material%get('phase')
 
   allocate(current(phases%length))
-
   allocate(damageState (phases%length))
   allocate(param(phases%length))
 
@@ -127,11 +127,14 @@ module subroutine damage_init
 
     phase => phases%get(ph)
     sources => phase%get('damage',defaultVal=emptyList)
-    param(ph)%K(1,1) = sources%get_asFloat('K_11',defaultVal=0.0_pReal)
-    param(ph)%K(2,2) = sources%get_asFloat('K_22',defaultVal=0.0_pReal)
-    param(ph)%K(3,3) = sources%get_asFloat('K_33',defaultVal=0.0_pReal)
     if (sources%length > 1) error stop
-    if (sources%length == 1) damage_active = .true.
+    if (sources%length == 1) then
+      damage_active = .true.
+      source => sources%get(1)
+      param(ph)%K(1,1) = source%get_asFloat('K_11',defaultVal=0.0_pReal)
+      param(ph)%K(2,2) = source%get_asFloat('K_22',defaultVal=0.0_pReal)
+      param(ph)%K(3,3) = source%get_asFloat('K_33',defaultVal=0.0_pReal)
+    endif
 
   enddo
 

From 8ec2d3a9ced292501d5ad16adea143bd15ede788 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sun, 11 Apr 2021 12:49:28 +0200
Subject: [PATCH 21/24] bugfix: mixed up order

---
 src/homogenization_thermal.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90
index 40c466ab5..27988dfb2 100644
--- a/src/homogenization_thermal.f90
+++ b/src/homogenization_thermal.f90
@@ -119,7 +119,7 @@ module function homogenization_K_T(ce) result(K)
   integer :: co
 
 
-  K = crystallite_push33ToRef(co,1,lattice_K_T(:,:,material_phaseID(1,ce)))
+  K = crystallite_push33ToRef(1,ce,lattice_K_T(:,:,material_phaseID(1,ce)))
   do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
     K = K + crystallite_push33ToRef(co,ce,lattice_K_T(:,:,material_phaseID(co,ce)))
   enddo

From b2292986f472777cc7552e49e156c778a5ad78b2 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sun, 11 Apr 2021 13:21:27 +0200
Subject: [PATCH 22/24] distributing responsibilities

---
 src/homogenization_damage.f90  |  4 +--
 src/homogenization_thermal.f90 | 55 +++++++---------------------------
 src/phase.f90                  | 28 ++++++++++++++++-
 src/phase_damage.f90           | 27 +++++++++++++++++
 src/phase_thermal.f90          | 29 ++++++++++++++++++
 5 files changed, 95 insertions(+), 48 deletions(-)

diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90
index 8a4e596a8..b8aecccae 100644
--- a/src/homogenization_damage.f90
+++ b/src/homogenization_damage.f90
@@ -114,7 +114,7 @@ module function homogenization_mu_phi(ce) result(mu)
   real(pReal) :: mu
 
 
-  mu = lattice_mu_phi(material_phaseID(1,ce))
+  mu = phase_mu_phi(1,ce)
 
 end function homogenization_mu_phi
 
@@ -128,7 +128,7 @@ module function homogenization_K_phi(ce) result(K)
   real(pReal), dimension(3,3) :: K
 
 
-  K = crystallite_push33ToRef(1,ce,lattice_K_phi(1:3,1:3,material_phaseID(1,ce))) \
+  K = phase_K_phi(1,ce) &
     * num_damage%charLength**2
 
 end function homogenization_K_phi
diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90
index 27988dfb2..0effe09ed 100644
--- a/src/homogenization_thermal.f90
+++ b/src/homogenization_thermal.f90
@@ -102,8 +102,15 @@ module function homogenization_mu_T(ce) result(mu)
   integer, intent(in) :: ce
   real(pReal) :: mu
 
+  integer :: co
 
-  mu = c_P(ce) * rho(ce)
+
+  mu = phase_mu_T(1,ce)
+  do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
+    mu = mu + phase_mu_T(co,ce)
+  enddo
+
+  mu = mu / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
 
 end function homogenization_mu_T
 
@@ -119,9 +126,9 @@ module function homogenization_K_T(ce) result(K)
   integer :: co
 
 
-  K = crystallite_push33ToRef(1,ce,lattice_K_T(:,:,material_phaseID(1,ce)))
+  K = phase_K_T(1,ce)
   do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
-    K = K + crystallite_push33ToRef(co,ce,lattice_K_T(:,:,material_phaseID(co,ce)))
+    K = K + phase_K_T(co,ce)
   enddo
 
   K = K / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
@@ -187,46 +194,4 @@ module subroutine thermal_results(ho,group)
 
 end subroutine thermal_results
 
-
-!--------------------------------------------------------------------------------------------------
-!> @brief Homogenize specific heat capacity.
-!--------------------------------------------------------------------------------------------------
-function c_P(ce)
-
-  integer, intent(in) :: ce
-  real(pReal) :: c_P
-
-  integer :: co
-
-
-  c_P = lattice_c_p(material_phaseID(1,ce))
-  do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
-    c_P = c_P + lattice_c_p(material_phaseID(co,ce))
-  enddo
-
-  c_P = c_P / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
-
-end function c_P
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief Homogenize mass density.
-!--------------------------------------------------------------------------------------------------
-function rho(ce)
-
-  integer, intent(in) :: ce
-  real(pReal) :: rho
-
-  integer :: co
-
-
-  rho = lattice_rho(material_phaseID(1,ce))
-  do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
-    rho = rho + lattice_rho(material_phaseID(co,ce))
-  enddo
-
-  rho = rho / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
-
-end function rho
-
 end submodule thermal
diff --git a/src/phase.f90 b/src/phase.f90
index 346772ce4..54116f3c9 100644
--- a/src/phase.f90
+++ b/src/phase.f90
@@ -184,7 +184,7 @@ module phase
 
     module subroutine phase_thermal_setField(T,dot_T, co,ce)
       real(pReal), intent(in) :: T, dot_T
-      integer, intent(in) :: ce, co
+      integer, intent(in) :: co, ce
     end subroutine phase_thermal_setField
 
     module subroutine phase_set_phi(phi,co,ce)
@@ -192,6 +192,28 @@ module phase
       integer, intent(in) :: co, ce
     end subroutine phase_set_phi
 
+
+    module function phase_mu_phi(co,ce) result(mu)
+      integer, intent(in) :: co, ce
+      real(pReal) :: mu
+    end function phase_mu_phi
+
+    module function phase_K_phi(co,ce) result(K)
+      integer, intent(in) :: co, ce
+      real(pReal), dimension(3,3) :: K
+    end function phase_K_phi
+
+
+    module function phase_mu_T(co,ce) result(mu)
+      integer, intent(in) :: co, ce
+      real(pReal) :: mu
+    end function phase_mu_T
+
+    module function phase_K_T(co,ce) result(K)
+      integer, intent(in) :: co, ce
+      real(pReal), dimension(3,3) :: K
+    end function phase_K_T
+
 ! == cleaned:end ===================================================================================
 
     module function thermal_stress(Delta_t,ph,me) result(converged_)
@@ -297,6 +319,10 @@ module phase
     phase_homogenizedC, &
     phase_f_phi, &
     phase_f_T, &
+    phase_K_phi, &
+    phase_K_T, &
+    phase_mu_phi, &
+    phase_mu_T, &
     phase_results, &
     phase_allocateState, &
     phase_forward, &
diff --git a/src/phase_damage.f90 b/src/phase_damage.f90
index 5f6d65ab8..9a20030d2 100644
--- a/src/phase_damage.f90
+++ b/src/phase_damage.f90
@@ -345,6 +345,33 @@ function phase_damage_collectDotState(ph,me) result(broken)
 end function phase_damage_collectDotState
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief Damage viscosity.
+!--------------------------------------------------------------------------------------------------
+module function phase_mu_phi(co,ce) result(mu)
+
+  integer, intent(in) :: co, ce
+  real(pReal) :: mu
+
+
+  mu = lattice_mu_phi(material_phaseID(co,ce))
+
+end function phase_mu_phi
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Damage conductivity/diffusivity in reference configuration.
+!--------------------------------------------------------------------------------------------------
+module function phase_K_phi(co,ce) result(K)
+
+  integer, intent(in) :: co, ce
+  real(pReal), dimension(3,3) :: K
+
+
+  K = crystallite_push33ToRef(co,ce,lattice_K_phi(1:3,1:3,material_phaseID(co,ce)))
+
+end function phase_K_phi
+
 
 !--------------------------------------------------------------------------------------------------
 !> @brief for constitutive models having an instantaneous change of state
diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90
index c2b4b34b2..a0e6adb60 100644
--- a/src/phase_thermal.f90
+++ b/src/phase_thermal.f90
@@ -184,6 +184,35 @@ function phase_thermal_collectDotState(ph,me) result(broken)
 end function phase_thermal_collectDotState
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief Damage viscosity.
+!--------------------------------------------------------------------------------------------------
+module function phase_mu_T(co,ce) result(mu)
+
+  integer, intent(in) :: co, ce
+  real(pReal) :: mu
+
+
+  mu = lattice_c_p(material_phaseID(co,ce)) &
+     * lattice_rho(material_phaseID(co,ce))
+
+end function phase_mu_T
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Damage conductivity/diffusivity in reference configuration.
+!--------------------------------------------------------------------------------------------------
+module function phase_K_T(co,ce) result(K)
+
+  integer, intent(in) :: co, ce
+  real(pReal), dimension(3,3) :: K
+
+
+  K = crystallite_push33ToRef(co,ce,lattice_K_T(1:3,1:3,material_phaseID(co,ce)))
+
+end function phase_K_T
+
+
 module function thermal_stress(Delta_t,ph,me) result(converged_)           ! ?? why is this called "stress" when it seems closer to "updateState" ??
 
   real(pReal), intent(in) :: Delta_t

From f655a6fe5cadd6f012726339786b698e6981deb1 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sun, 11 Apr 2021 13:55:30 +0200
Subject: [PATCH 23/24] store data locally

---
 src/lattice.f90       | 43 ++++++-------------------------------------
 src/phase_damage.f90  | 12 +++++++-----
 src/phase_thermal.f90 | 13 +++++++++----
 3 files changed, 22 insertions(+), 46 deletions(-)

diff --git a/src/lattice.f90 b/src/lattice.f90
index 69ee65857..06986b42b 100644
--- a/src/lattice.f90
+++ b/src/lattice.f90
@@ -398,9 +398,7 @@ module lattice
     lattice_rho, &
     lattice_c_p
    real(pReal),                       dimension(:,:,:), allocatable, public, protected :: &
-    lattice_C66, &
-    lattice_K_T, &
-    lattice_K_phi
+    lattice_C66
  integer(kind(lattice_UNDEFINED_ID)), dimension(:),     allocatable, public, protected :: &
     lattice_structure
 ! SHOULD NOT BE PART OF LATTICE END
@@ -458,23 +456,19 @@ subroutine lattice_init
     phases, &
     phase, &
     mech, &
-    elasticity, &
-    thermal, &
-    damage
+    elasticity
 
   print'(/,a)', ' <<<+-  lattice init  -+>>>'; flush(IO_STDOUT)
 
+! SHOULD NOT BE PART OF LATTICE BEGIN
+
   phases => config_material%get('phase')
   Nphases = phases%length
 
   allocate(lattice_structure(Nphases),source = lattice_UNDEFINED_ID)
   allocate(lattice_C66(6,6,Nphases),  source=0.0_pReal)
 
-  allocate(lattice_K_T  (3,3,Nphases), source=0.0_pReal)
-  allocate(lattice_K_phi  (3,3,Nphases), source=0.0_pReal)
-
-  allocate(lattice_mu_phi,&
-           lattice_rho,lattice_c_p, &
+  allocate(lattice_rho, &
            lattice_mu, lattice_nu,&
            source=[(0.0_pReal,i=1,Nphases)])
 
@@ -522,32 +516,7 @@ subroutine lattice_init
     enddo
 
     lattice_rho(ph) = phase%get_asFloat('rho', defaultVal=0.0_pReal)
-
-    ! SHOULD NOT BE PART OF LATTICE BEGIN
-
-    if (phase%contains('thermal')) then
-      thermal  => phase%get('thermal')
-      lattice_K_T(1,1,ph) = thermal%get_asFloat('K_11',defaultVal=0.0_pReal)
-      lattice_K_T(2,2,ph) = thermal%get_asFloat('K_22',defaultVal=0.0_pReal)
-      lattice_K_T(3,3,ph) = thermal%get_asFloat('K_33',defaultVal=0.0_pReal)
-      lattice_K_T(1:3,1:3,ph) = lattice_applyLatticeSymmetry33(lattice_K_T(1:3,1:3,ph), &
-                                                             phase%get_asString('lattice'))
-      lattice_c_p(ph) = thermal%get_asFloat('c_p', defaultVal=0.0_pReal)
-    endif
-
-
-    if (phase%contains('damage')) then
-      damage => phase%get('damage')
-      damage => damage%get(1)
-      lattice_K_phi(1,1,ph) = damage%get_asFloat('D_11',defaultVal=0.0_pReal)
-      lattice_K_phi(2,2,ph) = damage%get_asFloat('D_22',defaultVal=0.0_pReal)
-      lattice_K_phi(3,3,ph) = damage%get_asFloat('D_33',defaultVal=0.0_pReal)
-      lattice_K_phi(1:3,1:3,ph) = lattice_applyLatticeSymmetry33(lattice_K_phi(1:3,1:3,ph), &
-                                                          phase%get_asString('lattice'))
-
-      lattice_mu_phi(ph) = damage%get_asFloat('M',defaultVal=0.0_pReal)
-    endif
-    ! SHOULD NOT BE PART OF LATTICE END
+! SHOULD NOT BE PART OF LATTICE END
 
     call selfTest
 
diff --git a/src/phase_damage.f90 b/src/phase_damage.f90
index 9a20030d2..8df2cd5e4 100644
--- a/src/phase_damage.f90
+++ b/src/phase_damage.f90
@@ -131,9 +131,11 @@ module subroutine damage_init
     if (sources%length == 1) then
       damage_active = .true.
       source => sources%get(1)
-      param(ph)%K(1,1) = source%get_asFloat('K_11',defaultVal=0.0_pReal)
-      param(ph)%K(2,2) = source%get_asFloat('K_22',defaultVal=0.0_pReal)
-      param(ph)%K(3,3) = source%get_asFloat('K_33',defaultVal=0.0_pReal)
+      param(ph)%mu     = source%get_asFloat('M',defaultVal=0.0_pReal)
+      param(ph)%K(1,1) = source%get_asFloat('D_11',defaultVal=0.0_pReal)
+      param(ph)%K(2,2) = source%get_asFloat('D_22',defaultVal=0.0_pReal)
+      param(ph)%K(3,3) = source%get_asFloat('D_33',defaultVal=0.0_pReal)
+      param(ph)%K = lattice_applyLatticeSymmetry33(param(ph)%K,phase%get_asString('lattice'))
     endif
 
   enddo
@@ -354,7 +356,7 @@ module function phase_mu_phi(co,ce) result(mu)
   real(pReal) :: mu
 
 
-  mu = lattice_mu_phi(material_phaseID(co,ce))
+  mu = param(material_phaseID(co,ce))%mu
 
 end function phase_mu_phi
 
@@ -368,7 +370,7 @@ module function phase_K_phi(co,ce) result(K)
   real(pReal), dimension(3,3) :: K
 
 
-  K = crystallite_push33ToRef(co,ce,lattice_K_phi(1:3,1:3,material_phaseID(co,ce)))
+  K = crystallite_push33ToRef(co,ce,param(material_phaseID(co,ce))%K)
 
 end function phase_K_phi
 
diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90
index a0e6adb60..f6d046266 100644
--- a/src/phase_thermal.f90
+++ b/src/phase_thermal.f90
@@ -100,12 +100,17 @@ module subroutine thermal_init(phases)
     allocate(current(ph)%dot_T(Nmembers),source=0.0_pReal)
     phase => phases%get(ph)
     thermal => phase%get('thermal',defaultVal=emptyDict)
-    sources => thermal%get('source',defaultVal=emptyList)
+    param(ph)%C_p = thermal%get_asFloat('c_p',defaultVal=0.0_pReal)
+    if (param(ph)%C_p <= 0) param(ph)%C_p = thermal%get_asFloat('C_p',defaultVal=0.0_pReal)
     param(ph)%K(1,1) = thermal%get_asFloat('K_11',defaultVal=0.0_pReal)
     param(ph)%K(2,2) = thermal%get_asFloat('K_22',defaultVal=0.0_pReal)
     param(ph)%K(3,3) = thermal%get_asFloat('K_33',defaultVal=0.0_pReal)
+    param(ph)%K = lattice_applyLatticeSymmetry33(param(ph)%K,phase%get_asString('lattice'))
+
+    sources => thermal%get('source',defaultVal=emptyList)
     thermal_Nsources(ph) = sources%length
     allocate(thermalstate(ph)%p(thermal_Nsources(ph)))
+
   enddo
 
   allocate(thermal_source(maxval(thermal_Nsources),phases%length), source = THERMAL_UNDEFINED_ID)
@@ -193,8 +198,8 @@ module function phase_mu_T(co,ce) result(mu)
   real(pReal) :: mu
 
 
-  mu = lattice_c_p(material_phaseID(co,ce)) &
-     * lattice_rho(material_phaseID(co,ce))
+  mu = lattice_rho(material_phaseID(co,ce)) &
+     *  param(material_phaseID(co,ce))%C_p
 
 end function phase_mu_T
 
@@ -208,7 +213,7 @@ module function phase_K_T(co,ce) result(K)
   real(pReal), dimension(3,3) :: K
 
 
-  K = crystallite_push33ToRef(co,ce,lattice_K_T(1:3,1:3,material_phaseID(co,ce)))
+  K = crystallite_push33ToRef(co,ce,param(material_phaseID(co,ce))%K)
 
 end function phase_K_T
 

From 887524bcc1cd53bc6a842f0da5799d9082656ede Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sun, 11 Apr 2021 15:32:17 +0200
Subject: [PATCH 24/24] polishing

---
 src/grid/grid_damage_spectral.f90   | 27 ++++++++++----------
 src/grid/grid_thermal_spectral.f90  | 34 +++++++++++--------------
 src/homogenization_damage.f90       | 39 +++++++++++------------------
 src/homogenization_damage_pass.f90  |  5 ++--
 src/homogenization_thermal.f90      |  4 +--
 src/homogenization_thermal_pass.f90 |  5 ++--
 6 files changed, 52 insertions(+), 62 deletions(-)

diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90
index 7fe9a4a25..967ddb73a 100644
--- a/src/grid/grid_damage_spectral.f90
+++ b/src/grid/grid_damage_spectral.f90
@@ -15,8 +15,8 @@ module grid_damage_spectral
   use IO
   use spectral_utilities
   use discretization_grid
-  use YAML_types
   use homogenization
+  use YAML_types
   use config
 
   implicit none
@@ -61,7 +61,7 @@ contains
 !> @brief allocates all neccessary fields and fills them with data
 ! ToDo: Restart not implemented
 !--------------------------------------------------------------------------------------------------
-subroutine grid_damage_spectral_init
+subroutine grid_damage_spectral_init()
 
   PetscInt, dimension(0:worldsize-1) :: localK
   DM :: damage_grid
@@ -88,10 +88,10 @@ subroutine grid_damage_spectral_init
   num_generic => config_numerics%get('generic',defaultVal=emptyDict)
   num%residualStiffness = num_generic%get_asFloat('residualStiffness', defaultVal=1.0e-6_pReal)
 
-  if (num%residualStiffness < 0.0_pReal)   call IO_error(301,ext_msg='residualStiffness')
-  if (num%itmax <= 1)                      call IO_error(301,ext_msg='itmax')
-  if (num%eps_damage_atol <= 0.0_pReal)    call IO_error(301,ext_msg='eps_damage_atol')
-  if (num%eps_damage_rtol <= 0.0_pReal)    call IO_error(301,ext_msg='eps_damage_rtol')
+  if (num%residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness')
+  if (num%itmax <= 1)                    call IO_error(301,ext_msg='itmax')
+  if (num%eps_damage_atol <= 0.0_pReal)  call IO_error(301,ext_msg='eps_damage_atol')
+  if (num%eps_damage_rtol <= 0.0_pReal)  call IO_error(301,ext_msg='eps_damage_rtol')
 
 !--------------------------------------------------------------------------------------------------
 ! set default and user defined options for PETSc
@@ -146,6 +146,7 @@ subroutine grid_damage_spectral_init
   allocate(phi_current(grid(1),grid(2),grid3), source=1.0_pReal)
   allocate(phi_lastInc(grid(1),grid(2),grid3), source=1.0_pReal)
   allocate(phi_stagInc(grid(1),grid(2),grid3), source=1.0_pReal)
+
   call VecSet(solution_vec,1.0_pReal,ierr); CHKERRQ(ierr)
 
   call updateReference
@@ -216,21 +217,21 @@ end function grid_damage_spectral_solution
 subroutine grid_damage_spectral_forward(cutBack)
 
   logical, intent(in) :: cutBack
-  integer                               :: i, j, k, ce
+  integer :: i, j, k, ce
   DM :: dm_local
-  PetscScalar,  dimension(:,:,:), pointer     :: x_scal
-  PetscErrorCode                              :: ierr
+  PetscScalar,  dimension(:,:,:), pointer :: x_scal
+  PetscErrorCode :: ierr
 
   if (cutBack) then
     phi_current = phi_lastInc
     phi_stagInc = phi_lastInc
 !--------------------------------------------------------------------------------------------------
 ! reverting damage field state
-    ce = 0
     call SNESGetDM(damage_snes,dm_local,ierr); CHKERRQ(ierr)
     call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr)                       !< get the data out of PETSc to work with
     x_scal(xstart:xend,ystart:yend,zstart:zend) = phi_current
     call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr)
+    ce = 0
     do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
       ce = ce + 1
       call homogenization_set_phi(phi_current(i,j,k),ce)
@@ -271,8 +272,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
   ce = 0
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
     ce = ce + 1
-    vectorField_real(1:3,i,j,k) = matmul(homogenization_K_phi(ce) - K_ref, &
-                                         vectorField_real(1:3,i,j,k))
+    vectorField_real(1:3,i,j,k) = matmul(homogenization_K_phi(ce) - K_ref, vectorField_real(1:3,i,j,k))
   enddo; enddo; enddo
   call utilities_FFTvectorForward
   call utilities_fourierVectorDivergence                                                            !< calculate damage divergence in fourier field
@@ -290,6 +290,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
   call utilities_FFTscalarForward
   call utilities_fourierGreenConvolution(K_ref, mu_ref, params%timeinc)
   call utilities_FFTscalarBackward
+
   where(scalarField_real(1:grid(1),1:grid(2),1:grid3) > phi_lastInc) &
         scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_lastInc
   where(scalarField_real(1:grid(1),1:grid(2),1:grid3) < num%residualStiffness) &
@@ -305,7 +306,7 @@ end subroutine formResidual
 !--------------------------------------------------------------------------------------------------
 !> @brief update reference viscosity and conductivity
 !--------------------------------------------------------------------------------------------------
-subroutine updateReference
+subroutine updateReference()
 
   integer :: ce,ierr
 
diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90
index 60d8d99db..1ccb3b901 100644
--- a/src/grid/grid_thermal_spectral.f90
+++ b/src/grid/grid_thermal_spectral.f90
@@ -18,7 +18,6 @@ module grid_thermal_spectral
   use homogenization
   use YAML_types
   use config
-  use material
 
   implicit none
   private
@@ -46,7 +45,7 @@ module grid_thermal_spectral
 
 !--------------------------------------------------------------------------------------------------
 ! reference diffusion tensor, mobility etc.
-  integer                     :: totalIter = 0                                             !< total iteration in current increment
+  integer                     :: totalIter = 0                                                      !< total iteration in current increment
   real(pReal), dimension(3,3) :: K_ref
   real(pReal)                 :: mu_ref
 
@@ -68,7 +67,7 @@ subroutine grid_thermal_spectral_init(T_0)
   PetscInt, dimension(0:worldsize-1) :: localK
   integer :: i, j, k, ce
   DM :: thermal_grid
-  PetscScalar,  dimension(:,:,:), pointer :: x_scal
+  PetscScalar, dimension(:,:,:), pointer :: T_PETSc
   PetscErrorCode :: ierr
   class(tNode), pointer :: &
     num_grid
@@ -85,9 +84,9 @@ subroutine grid_thermal_spectral_init(T_0)
   num%eps_thermal_atol = num_grid%get_asFloat ('eps_thermal_atol',defaultVal=1.0e-2_pReal)
   num%eps_thermal_rtol = num_grid%get_asFloat ('eps_thermal_rtol',defaultVal=1.0e-6_pReal)
 
-  if (num%itmax <= 1)                         call IO_error(301,ext_msg='itmax')
-  if (num%eps_thermal_atol <= 0.0_pReal)      call IO_error(301,ext_msg='eps_thermal_atol')
-  if (num%eps_thermal_rtol <= 0.0_pReal)      call IO_error(301,ext_msg='eps_thermal_rtol')
+  if (num%itmax <= 1)                    call IO_error(301,ext_msg='itmax')
+  if (num%eps_thermal_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_thermal_atol')
+  if (num%eps_thermal_rtol <= 0.0_pReal) call IO_error(301,ext_msg='eps_thermal_rtol')
 
 !--------------------------------------------------------------------------------------------------
 ! set default and user defined options for PETSc
@@ -130,6 +129,7 @@ subroutine grid_thermal_spectral_init(T_0)
   allocate(T_current(grid(1),grid(2),grid3), source=0.0_pReal)
   allocate(T_lastInc(grid(1),grid(2),grid3), source=0.0_pReal)
   allocate(T_stagInc(grid(1),grid(2),grid3), source=0.0_pReal)
+
   ce = 0
   do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
     ce = ce + 1
@@ -138,9 +138,10 @@ subroutine grid_thermal_spectral_init(T_0)
     T_stagInc(i,j,k) = T_current(i,j,k)
     call homogenization_thermal_setField(T_0,0.0_pReal,ce)
   enddo; enddo; enddo
-  call DMDAVecGetArrayF90(thermal_grid,solution_vec,x_scal,ierr); CHKERRQ(ierr)                     !< get the data out of PETSc to work with
-  x_scal(xstart:xend,ystart:yend,zstart:zend) = T_current
-  call DMDAVecRestoreArrayF90(thermal_grid,solution_vec,x_scal,ierr); CHKERRQ(ierr)
+
+  call DMDAVecGetArrayF90(thermal_grid,solution_vec,T_PETSc,ierr); CHKERRQ(ierr)
+  T_PETSc(xstart:xend,ystart:yend,zstart:zend) = T_current
+  call DMDAVecRestoreArrayF90(thermal_grid,solution_vec,T_PETSc,ierr); CHKERRQ(ierr)
 
   call updateReference
 
@@ -190,9 +191,7 @@ function grid_thermal_spectral_solution(timeinc) result(solution)
   ce = 0
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
     ce = ce + 1
-    call homogenization_thermal_setField(T_current(i,j,k), &
-                                                     (T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, &
-                                                     ce)
+    call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc,ce)
   enddo; enddo; enddo
 
   call VecMin(solution_vec,devNull,T_min,ierr); CHKERRQ(ierr)
@@ -223,16 +222,14 @@ subroutine grid_thermal_spectral_forward(cutBack)
 
 !--------------------------------------------------------------------------------------------------
 ! reverting thermal field state
-    ce = 0
     call SNESGetDM(thermal_snes,dm_local,ierr); CHKERRQ(ierr)
     call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr)                       !< get the data out of PETSc to work with
     x_scal(xstart:xend,ystart:yend,zstart:zend) = T_current
     call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr)
+    ce = 0
     do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
       ce = ce + 1
-      call homogenization_thermal_setField(T_current(i,j,k), &
-                                          (T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc, &
-                                                       ce)
+      call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc,ce)
     enddo; enddo; enddo
   else
     T_lastInc = T_current
@@ -270,8 +267,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
   ce = 0
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
     ce = ce + 1
-    vectorField_real(1:3,i,j,k) = matmul(homogenization_K_T(ce) - K_ref, &
-                                         vectorField_real(1:3,i,j,k))
+    vectorField_real(1:3,i,j,k) = matmul(homogenization_K_T(ce) - K_ref, vectorField_real(1:3,i,j,k))
   enddo; enddo; enddo
   call utilities_FFTvectorForward
   call utilities_fourierVectorDivergence                                                            !< calculate temperature divergence in fourier field
@@ -300,7 +296,7 @@ end subroutine formResidual
 !--------------------------------------------------------------------------------------------------
 !> @brief update reference viscosity and conductivity
 !--------------------------------------------------------------------------------------------------
-subroutine updateReference
+subroutine updateReference()
 
   integer :: ce,ierr
 
diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90
index b8aecccae..0546834fd 100644
--- a/src/homogenization_damage.f90
+++ b/src/homogenization_damage.f90
@@ -38,14 +38,12 @@ module subroutine damage_init()
     configHomogenizations, &
     configHomogenization, &
     configHomogenizationDamage, &
-    num_generic, &
-    material_homogenization
-  integer :: ho
-  integer :: Ninstances,Nmaterialpoints,h
+    num_generic
+  integer :: ho,Nmaterialpoints
 
 
   print'(/,a)', ' <<<+-  homogenization:damage init  -+>>>'
-  print'(/,a)', ' <<<+-  homogenization:damage:pass init  -+>>>'
+
 
   configHomogenizations => config_material%get('homogenization')
   allocate(param(configHomogenizations%length))
@@ -62,6 +60,10 @@ module subroutine damage_init()
 #else
         prm%output = configHomogenizationDamage%get_as1dString('output',defaultVal=emptyStringArray)
 #endif
+        Nmaterialpoints = count(material_homogenizationAt == ho)
+        damageState_h(ho)%sizeState = 1
+        allocate(damageState_h(ho)%state0(1,Nmaterialpoints), source=1.0_pReal)
+        allocate(damageState_h(ho)%state (1,Nmaterialpoints), source=1.0_pReal)
       else
         prm%output = emptyStringArray
       endif
@@ -73,18 +75,7 @@ module subroutine damage_init()
   num_generic => config_numerics%get('generic',defaultVal= emptyDict)
   num_damage%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal)
 
-  Ninstances = count(damage_type == DAMAGE_nonlocal_ID)
-
-  material_homogenization => config_material%get('homogenization')
-  do h = 1, material_homogenization%length
-    if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle
-
-    Nmaterialpoints = count(material_homogenizationAt == h)
-    damageState_h(h)%sizeState = 1
-    allocate(damageState_h(h)%state0   (1,Nmaterialpoints), source=1.0_pReal)
-    allocate(damageState_h(h)%state    (1,Nmaterialpoints), source=1.0_pReal)
-
-  enddo
+  call pass_init()
 
 end subroutine damage_init
 
@@ -183,13 +174,13 @@ module subroutine damage_results(ho,group)
   integer :: o
 
   associate(prm => param(ho))
-  outputsLoop: do o = 1,size(prm%output)
-    select case(prm%output(o))
-      case ('phi')
-        call results_writeDataset(group,damagestate_h(ho)%state(1,:),prm%output(o),&
-                                  'damage indicator','-')
-    end select
-  enddo outputsLoop
+      outputsLoop: do o = 1,size(prm%output)
+        select case(prm%output(o))
+          case ('phi')
+            call results_writeDataset(group,damagestate_h(ho)%state(1,:),prm%output(o),&
+                                      'damage indicator','-')
+        end select
+      enddo outputsLoop
   end associate
 
 end subroutine damage_results
diff --git a/src/homogenization_damage_pass.f90 b/src/homogenization_damage_pass.f90
index cbb7ec79f..bdb44b822 100644
--- a/src/homogenization_damage_pass.f90
+++ b/src/homogenization_damage_pass.f90
@@ -6,8 +6,9 @@ submodule(homogenization:damage) damage_pass
 
 contains
 
-module subroutine pass_init
-
+module subroutine pass_init()
+  
+  print'(/,a)', ' <<<+-  homogenization:damage:pass init  -+>>>'
 
 end subroutine pass_init
 
diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90
index 0effe09ed..f681036de 100644
--- a/src/homogenization_thermal.f90
+++ b/src/homogenization_thermal.f90
@@ -45,8 +45,6 @@ module subroutine thermal_init()
 
 
   print'(/,a)', ' <<<+-  homogenization:thermal init  -+>>>'
-  print'(/,a)', ' <<<+-  homogenization:thermal:pass init  -+>>>'
-
 
 
   configHomogenizations => config_material%get('homogenization')
@@ -71,6 +69,8 @@ module subroutine thermal_init()
     end associate
   enddo
 
+  call pass_init()
+
 end subroutine thermal_init
 
 
diff --git a/src/homogenization_thermal_pass.f90 b/src/homogenization_thermal_pass.f90
index 940a15024..87020d8d2 100644
--- a/src/homogenization_thermal_pass.f90
+++ b/src/homogenization_thermal_pass.f90
@@ -6,8 +6,9 @@ submodule(homogenization:thermal) thermal_pass
 
 contains
 
-module subroutine pass_init
-
+module subroutine pass_init()
+  
+  print'(/,a)', ' <<<+-  homogenization:thermal:pass init  -+>>>'
 
 end subroutine pass_init