avoiding public variables

src/constitutive.f90

@@ -45,12 +45,8 @@ module constitutive
 
   type(rotation),            dimension(:,:,:),        allocatable :: &
     crystallite_orientation                                                                         !< current orientation
-  real(pReal),               dimension(:,:,:,:,:),    allocatable :: &
-    crystallite_F0                                                                                  !< def grad at start of FE inc
   real(pReal),               dimension(:,:,:,:,:),    allocatable, public :: &
-    crystallite_P, &                                                                                !< 1st Piola-Kirchhoff stress per grain
+    crystallite_P                                                                                   !< 1st Piola-Kirchhoff stress per grain
-    crystallite_partitionedF0, &                                                                    !< def grad at start of homog inc
-    crystallite_F                                                                                   !< def grad to be reached at end of homog inc
 
   type :: tTensorContainer
     real(pReal), dimension(:,:,:), allocatable :: data
@@ -61,18 +57,21 @@ module constitutive
     constitutive_mech_Fe, &
     constitutive_mech_Fi, &
     constitutive_mech_Fp, &
+    constitutive_mech_F, &
     constitutive_mech_Li, &
     constitutive_mech_Lp, &
     constitutive_mech_S, &
     ! converged value at end of last solver increment
     constitutive_mech_Fi0, &
     constitutive_mech_Fp0, &
+    constitutive_mech_F0, &
     constitutive_mech_Li0, &
     constitutive_mech_Lp0, &
     constitutive_mech_S0, &
     ! converged value at end of last homogenization increment (RGC only)
     constitutive_mech_partitionedFi0, &
     constitutive_mech_partitionedFp0, &
+    constitutive_mech_partitionedF0, &
     constitutive_mech_partitionedLi0, &
     constitutive_mech_partitionedLp0, &
     constitutive_mech_partitionedS0
@@ -339,13 +338,11 @@ module constitutive
     end subroutine constitutive_plastic_LpAndItsTangents
 
 
-    module subroutine constitutive_plastic_dependentState(F, co, ip, el)
+    module subroutine constitutive_plastic_dependentState(co,ip,el)
       integer, intent(in) :: &
         co, &                                                                                      !< component-ID of integration point
         ip, &                                                                                       !< integration point
         el                                                                                          !< element
-      real(pReal),   intent(in), dimension(3,3) :: &
-        F                                                                                           !< elastic deformation gradient
     end subroutine constitutive_plastic_dependentState
 
 
@@ -394,6 +391,7 @@ module constitutive
     integrateSourceState, &
     constitutive_mech_setF, &
     constitutive_mech_getLp, &
+    constitutive_mech_getF, &
     constitutive_mech_getS, &
     crystallite_restartRead, &
     constitutive_initializeRestorationPoints, &
@@ -789,15 +787,14 @@ end subroutine constitutive_restore
 !--------------------------------------------------------------------------------------------------
 subroutine constitutive_forward
 
-  integer :: i, j
+  integer :: ph, so
 
-  crystallite_F0  = crystallite_F
   
   call constitutive_mech_forward()
 
-  do i = 1, size(sourceState)
+  do ph = 1, size(sourceState)
-    do j = 1,phase_Nsources(i)
+    do so = 1,phase_Nsources(ph)
-      sourceState(i)%p(j)%state0 = sourceState(i)%p(j)%state
+      sourceState(ph)%p(so)%state0 = sourceState(ph)%p(so)%state
   enddo; enddo
 
 end subroutine constitutive_forward
@@ -862,12 +859,6 @@ subroutine crystallite_init
   eMax = discretization_Nelems
 
   allocate(crystallite_P(3,3,cMax,iMax,eMax),source=0.0_pReal)
-
-  allocate(crystallite_F0, &
-           crystallite_partitionedF0,&
-           crystallite_F, &
-           source = crystallite_P)
-
   allocate(crystallite_orientation(cMax,iMax,eMax))
 
 
@@ -911,6 +902,9 @@ subroutine crystallite_init
   allocate(constitutive_mech_Fp(phases%length))
   allocate(constitutive_mech_Fp0(phases%length))
   allocate(constitutive_mech_partitionedFp0(phases%length))
+  allocate(constitutive_mech_F(phases%length))
+  allocate(constitutive_mech_F0(phases%length))
+  allocate(constitutive_mech_partitionedF0(phases%length))
   allocate(constitutive_mech_Li(phases%length))
   allocate(constitutive_mech_Li0(phases%length))
   allocate(constitutive_mech_partitionedLi0(phases%length))
@@ -939,6 +933,9 @@ subroutine crystallite_init
     allocate(constitutive_mech_S(ph)%data(3,3,Nconstituents),source=0.0_pReal)
     allocate(constitutive_mech_S0(ph)%data(3,3,Nconstituents),source=0.0_pReal)
     allocate(constitutive_mech_partitionedS0(ph)%data(3,3,Nconstituents),source=0.0_pReal)
+    allocate(constitutive_mech_F(ph)%data(3,3,Nconstituents))
+    allocate(constitutive_mech_F0(ph)%data(3,3,Nconstituents))
+    allocate(constitutive_mech_partitionedF0(ph)%data(3,3,Nconstituents))
     do so = 1, phase_Nsources(ph)
       allocate(sourceState(ph)%p(so)%subState0,source=sourceState(ph)%p(so)%state0)                 ! ToDo: hack
     enddo
@@ -955,28 +952,27 @@ subroutine crystallite_init
 
       ph = material_phaseAt(co,el)
       me = material_phaseMemberAt(co,ip,el)
+      
       constitutive_mech_Fp0(ph)%data(1:3,1:3,me) = material_orientation0(co,ip,el)%asMatrix()                      ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005)
       constitutive_mech_Fp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) &
                                                  / math_det33(constitutive_mech_Fp0(ph)%data(1:3,1:3,me))**(1.0_pReal/3.0_pReal)
       constitutive_mech_Fi0(ph)%data(1:3,1:3,me) = math_I3
-
+      constitutive_mech_F0(ph)%data(1:3,1:3,me)  = math_I3
-      crystallite_F0(1:3,1:3,co,ip,el)  = math_I3
       
       constitutive_mech_Fe(ph)%data(1:3,1:3,me) = math_inv33(matmul(constitutive_mech_Fi0(ph)%data(1:3,1:3,me), &
                                                                     constitutive_mech_Fp0(ph)%data(1:3,1:3,me)))           ! assuming that euler angles are given in internal strain free configuration
       constitutive_mech_Fp(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me)
       constitutive_mech_Fi(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me)
+      constitutive_mech_F(ph)%data(1:3,1:3,me)  = constitutive_mech_F0(ph)%data(1:3,1:3,me)
 
       constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me)
       constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me)
+      constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me)  = constitutive_mech_F0(ph)%data(1:3,1:3,me)
 
     enddo
   enddo; enddo
   !$OMP END PARALLEL DO
 
-  crystallite_partitionedF0  = crystallite_F0
-  crystallite_F   = crystallite_F0
-
 
   !$OMP PARALLEL DO PRIVATE(ph,me)
   do el = 1, size(material_phaseMemberAt,3)
@@ -985,7 +981,7 @@ subroutine crystallite_init
         ph = material_phaseAt(co,el)
         me = material_phaseMemberAt(co,ip,el)
         call crystallite_orientations(co,ip,el)
-        call constitutive_plastic_dependentState(crystallite_partitionedF0(1:3,1:3,co,ip,el),co,ip,el)  ! update dependent state variables to be consistent with basic states
+        call constitutive_plastic_dependentState(co,ip,el)                                          ! update dependent state variables to be consistent with basic states
      enddo
     enddo
   enddo
@@ -1010,13 +1006,11 @@ subroutine constitutive_initializeRestorationPoints(ip,el)
   do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
     ph = material_phaseAt(co,el)
     me = material_phaseMemberAt(co,ip,el)
-    crystallite_partitionedF0(1:3,1:3,co,ip,el)  = crystallite_F0(1:3,1:3,co,ip,el)
 
     call mech_initializeRestorationPoints(ph,me)
 
     do so = 1, phase_Nsources(material_phaseAt(co,el))
-      sourceState(material_phaseAt(co,el))%p(so)%partitionedState0(:,material_phasememberAt(co,ip,el)) = &
+      sourceState(ph)%p(so)%partitionedState0(:,me) = sourceState(ph)%p(so)%state0(:,me)
-      sourceState(material_phaseAt(co,el))%p(so)%state0(           :,material_phasememberAt(co,ip,el))
     enddo
   enddo
 
@@ -1040,7 +1034,6 @@ subroutine constitutive_windForward(ip,el)
   do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
     ph = material_phaseAt(co,el)
     me = material_phaseMemberAt(co,ip,el)
-    crystallite_partitionedF0 (1:3,1:3,co,ip,el) = crystallite_F (1:3,1:3,co,ip,el)
 
     call constitutive_mech_windForward(ph,me)
     do so = 1, phase_Nsources(material_phaseAt(co,el))
@@ -1132,8 +1125,8 @@ function crystallite_stressTangent(dt,co,ip,el) result(dPdF)
 !--------------------------------------------------------------------------------------------------
 ! calculate dSdF
   temp_33_1 = transpose(matmul(invFp,invFi))
-  temp_33_2 = matmul(crystallite_F(1:3,1:3,co,ip,el),invSubFp0)
+  temp_33_2 = matmul(constitutive_mech_F(ph)%data(1:3,1:3,me),invSubFp0)
-  temp_33_3 = matmul(matmul(crystallite_F(1:3,1:3,co,ip,el),invFp), invSubFi0)
+  temp_33_3 = matmul(matmul(constitutive_mech_F(ph)%data(1:3,1:3,me),invFp), invSubFi0)
 
   do o=1,3; do p=1,3
     rhs_3333(p,o,1:3,1:3)  = matmul(dSdFe(p,o,1:3,1:3),temp_33_1)
@@ -1162,7 +1155,7 @@ function crystallite_stressTangent(dt,co,ip,el) result(dPdF)
 !--------------------------------------------------------------------------------------------------
 ! assemble dPdF
   temp_33_1 = matmul(constitutive_mech_S(ph)%data(1:3,1:3,me),transpose(invFp))
-  temp_33_2 = matmul(crystallite_F(1:3,1:3,co,ip,el),invFp)
+  temp_33_2 = matmul(constitutive_mech_F(ph)%data(1:3,1:3,me),invFp)
   temp_33_3 = matmul(temp_33_2,constitutive_mech_S(ph)%data(1:3,1:3,me))
 
   dPdF = 0.0_pReal
@@ -1171,7 +1164,7 @@ function crystallite_stressTangent(dt,co,ip,el) result(dPdF)
   enddo
   do o=1,3; do p=1,3
     dPdF(1:3,1:3,p,o) = dPdF(1:3,1:3,p,o) &
-                      + matmul(matmul(crystallite_F(1:3,1:3,co,ip,el),dFpinvdF(1:3,1:3,p,o)),temp_33_1) &
+                      + matmul(matmul(constitutive_mech_F(ph)%data(1:3,1:3,me),dFpinvdF(1:3,1:3,p,o)),temp_33_1) &
                       + matmul(matmul(temp_33_2,dSdF(1:3,1:3,p,o)),transpose(invFp)) &
                       + matmul(temp_33_3,transpose(dFpinvdF(1:3,1:3,p,o)))
   enddo; enddo
@@ -1207,17 +1200,17 @@ end subroutine crystallite_orientations
 function crystallite_push33ToRef(co,ip,el, tensor33)
 
   real(pReal), dimension(3,3), intent(in) :: tensor33
-  real(pReal), dimension(3,3)             :: T
   integer, intent(in):: &
     el, &
     ip, &
     co
-
   real(pReal), dimension(3,3) :: crystallite_push33ToRef
   
+  real(pReal), dimension(3,3)             :: T
+
 
   T = matmul(material_orientation0(co,ip,el)%asMatrix(), &                                         ! ToDo: initial orientation correct?
-             transpose(math_inv33(crystallite_F(1:3,1:3,co,ip,el))))
+             transpose(math_inv33(constitutive_mech_F(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)))))
   crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T))
 
 end function crystallite_push33ToRef
@@ -1360,8 +1353,6 @@ subroutine crystallite_restartWrite
   write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5'
   fileHandle = HDF5_openFile(fileName,'a')
 
-  call HDF5_write(fileHandle,crystallite_F,'F')
-
   groupHandle = HDF5_addGroup(fileHandle,'phase')
   do ph = 1,size(material_name_phase)
     write(datasetName,'(i0,a)') ph,'_omega'
@@ -1376,6 +1367,8 @@ subroutine crystallite_restartWrite
     call HDF5_write(groupHandle,constitutive_mech_Fp(ph)%data,datasetName)
     write(datasetName,'(i0,a)') ph,'_S'
     call HDF5_write(groupHandle,constitutive_mech_S(ph)%data,datasetName)
+    write(datasetName,'(i0,a)') ph,'_F'
+    call HDF5_write(groupHandle,constitutive_mech_F(ph)%data,datasetName)
   enddo
   call HDF5_closeGroup(groupHandle)
 
@@ -1406,8 +1399,6 @@ subroutine crystallite_restartRead
   write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5'
   fileHandle = HDF5_openFile(fileName)
 
-  call HDF5_read(fileHandle,crystallite_F0, 'F')
-
   groupHandle = HDF5_openGroup(fileHandle,'phase')
   do ph = 1,size(material_name_phase)
     write(datasetName,'(i0,a)') ph,'_omega'
@@ -1422,6 +1413,8 @@ subroutine crystallite_restartRead
     call HDF5_read(groupHandle,constitutive_mech_Fp0(ph)%data,datasetName)
     write(datasetName,'(i0,a)') ph,'_S'
     call HDF5_read(groupHandle,constitutive_mech_S0(ph)%data,datasetName)
+    write(datasetName,'(i0,a)') ph,'_F'
+    call HDF5_read(groupHandle,constitutive_mech_F0(ph)%data,datasetName)
   enddo
   call HDF5_closeGroup(groupHandle)
 
@@ -1461,6 +1454,18 @@ function constitutive_mech_getLp(co,ip,el) result(Lp)
 end function constitutive_mech_getLp
 
 
+! getter for non-mech (e.g. thermal)
+function constitutive_mech_getF(co,ip,el) result(F)
+
+  integer, intent(in) :: co, ip, el
+  real(pReal), dimension(3,3) :: F
+
+
+  F = constitutive_mech_F(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el))
+
+end function constitutive_mech_getF
+
+
 ! setter for homogenization
 subroutine constitutive_mech_setF(F,co,ip,el)
 
@@ -1468,8 +1473,7 @@ subroutine constitutive_mech_setF(F,co,ip,el)
   integer, intent(in) :: co, ip, el
 
 
-  crystallite_F(1:3,1:3,co,ip,el) = F
+  constitutive_mech_F(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) = F
-  !constitutive_mech_F(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) = F
 
 end subroutine constitutive_mech_setF
 

src/constitutive_mech.f90

@@ -184,12 +184,9 @@ submodule(constitutive) constitutive_mech
         of
     end subroutine plastic_disloTungsten_dotState
 
-    module subroutine plastic_nonlocal_dotState(Mp, F, Temperature,timestep, &
+    module subroutine plastic_nonlocal_dotState(Mp,Temperature,timestep,instance,of,ip,el)
-                                                        instance,of,ip,el)
       real(pReal), dimension(3,3), intent(in) :: &
         Mp                                                                                          !< MandelStress
-      real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: &
-        F                                                                                           !< deformation gradient
       real(pReal), intent(in) :: &
         Temperature, &                                                                              !< temperature
         timestep                                                                                    !< substepped crystallite time increment
@@ -215,9 +212,7 @@ submodule(constitutive) constitutive_mech
         of
     end subroutine plastic_dislotungsten_dependentState
 
-    module subroutine plastic_nonlocal_dependentState(F, instance, of, ip, el)
+    module subroutine plastic_nonlocal_dependentState(instance, of, ip, el)
-      real(pReal), dimension(3,3), intent(in) :: &
-        F                                                                                           !< deformation gradient
       integer, intent(in) :: &
         instance, &
         of, &
@@ -480,32 +475,35 @@ end subroutine constitutive_hooke_SandItsTangents
 !--------------------------------------------------------------------------------------------------
 !> @brief calls microstructure function of the different plasticity constitutive models
 !--------------------------------------------------------------------------------------------------
-module subroutine constitutive_plastic_dependentState(F, co, ip, el)
+module subroutine constitutive_plastic_dependentState(co, ip, el)
 
   integer, intent(in) :: &
     co, &                                                                                           !< component-ID of integration point
     ip, &                                                                                           !< integration point
     el                                                                                              !< element
-  real(pReal),   intent(in), dimension(3,3) :: &
-    F                                                                                               !< deformation gradient
 
   integer :: &
     ho, &                                                                                           !< homogenization
     tme, &                                                                                          !< thermal member position
-    instance, of
+    instance, me
+
 
   ho  = material_homogenizationAt(el)
   tme = material_homogenizationMemberAt(ip,el)
-  of  = material_phasememberAt(co,ip,el)
+  me  = material_phasememberAt(co,ip,el)
   instance = phase_plasticityInstance(material_phaseAt(co,el))
 
   plasticityType: select case (phase_plasticity(material_phaseAt(co,el)))
+    
     case (PLASTICITY_DISLOTWIN_ID) plasticityType
-      call plastic_dislotwin_dependentState(temperature(ho)%p(tme),instance,of)
+      call plastic_dislotwin_dependentState(temperature(ho)%p(tme),instance,me)
+    
     case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType
-      call plastic_dislotungsten_dependentState(instance,of)
+      call plastic_dislotungsten_dependentState(instance,me)
+    
     case (PLASTICITY_NONLOCAL_ID) plasticityType
-      call plastic_nonlocal_dependentState (F,instance,of,ip,el)
+      call plastic_nonlocal_dependentState(instance,me,ip,el)
+  
   end select plasticityType
 
 end subroutine constitutive_plastic_dependentState
@@ -539,13 +537,13 @@ module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
     ho, &                                                                                           !< homogenization
     tme                                                                                             !< thermal member position
   integer :: &
-    i, j, instance, of
+    i, j, instance, me
 
   ho = material_homogenizationAt(el)
   tme = material_homogenizationMemberAt(ip,el)
 
   Mp = matmul(matmul(transpose(Fi),Fi),S)
-  of = material_phasememberAt(co,ip,el)
+  me = material_phasememberAt(co,ip,el)
   instance = phase_plasticityInstance(material_phaseAt(co,el))
 
   plasticityType: select case (phase_plasticity(material_phaseAt(co,el)))
@@ -555,22 +553,22 @@ module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
       dLp_dMp = 0.0_pReal
 
     case (PLASTICITY_ISOTROPIC_ID) plasticityType
-      call plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
+      call plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me)
 
     case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
-      call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
+      call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me)
 
     case (PLASTICITY_KINEHARDENING_ID) plasticityType
-      call plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
+      call plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,me)
 
     case (PLASTICITY_NONLOCAL_ID) plasticityType
-      call plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, temperature(ho)%p(tme),instance,of,ip,el)
+      call plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp,Mp, temperature(ho)%p(tme),instance,me,ip,el)
 
     case (PLASTICITY_DISLOTWIN_ID) plasticityType
-      call plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of)
+      call plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,me)
 
     case (PLASTICITY_DISLOTUNGSTEN_ID) plasticityType
-      call plastic_dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of)
+      call plastic_dislotungsten_LpAndItsTangent(Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,me)
 
   end select plasticityType
 
@@ -586,7 +584,7 @@ end subroutine constitutive_plastic_LpAndItsTangents
 !--------------------------------------------------------------------------------------------------
 !> @brief contains the constitutive equation for calculating the rate of change of microstructure
 !--------------------------------------------------------------------------------------------------
-function mech_collectDotState(subdt, co, ip, el,ph,of) result(broken)
+function mech_collectDotState(subdt,co,ip,el,ph,of) result(broken)
 
   integer, intent(in) :: &
     co, &                                                                                          !< component-ID of integration point
@@ -601,9 +599,9 @@ function mech_collectDotState(subdt, co, ip, el,ph,of) result(broken)
   integer :: &
     ho, &                                                                                           !< homogenization
     tme, &                                                                                          !< thermal member position
-    i, &                                                                                            !< counter in source loop
     instance
   logical :: broken
+  
   ho = material_homogenizationAt(el)
   tme = material_homogenizationMemberAt(ip,el)
   instance = phase_plasticityInstance(ph)
@@ -629,8 +627,7 @@ function mech_collectDotState(subdt, co, ip, el,ph,of) result(broken)
       call plastic_disloTungsten_dotState(Mp,temperature(ho)%p(tme),instance,of)
 
     case (PLASTICITY_NONLOCAL_ID) plasticityType
-      call plastic_nonlocal_dotState(Mp,crystallite_partitionedF0,temperature(ho)%p(tme),subdt, &
+      call plastic_nonlocal_dotState(Mp,temperature(ho)%p(tme),subdt,instance,of,ip,el)
-                                          instance,of,ip,el)
   end select plasticityType
   broken = any(IEEE_is_NaN(plasticState(ph)%dotState(:,of)))
 
@@ -798,13 +795,14 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
                                       jacoCounterLi                                                 ! counters to check for Jacobian update
   logical :: error,broken
 
-  broken = .true.
 
-  call constitutive_plastic_dependentState(crystallite_F(1:3,1:3,co,ip,el),co,ip,el)
+  broken = .true.
 
   ph = material_phaseAt(co,el)
   me = material_phaseMemberAt(co,ip,el)
   
+  call constitutive_plastic_dependentState(co,ip,el)
+
   Lpguess = constitutive_mech_Lp(ph)%data(1:3,1:3,me)                                              ! take as first guess
   Liguess = constitutive_mech_Li(ph)%data(1:3,1:3,me)                                              ! take as first guess
 
@@ -1289,8 +1287,7 @@ subroutine crystallite_results(group,ph)
 
       select case (output_constituent(ph)%label(ou))
         case('F')
-          selected_tensors = select_tensors(crystallite_F,ph)
+          call results_writeDataset(group//'/mechanics/',constitutive_mech_F(ph)%data,output_constituent(ph)%label(ou),&
-          call results_writeDataset(group//'/mechanics/',selected_tensors,output_constituent(ph)%label(ou),&
                                    'deformation gradient','1')
         case('F_e')
           call results_writeDataset(group//'/mechanics/',constitutive_mech_Fe(ph)%data,output_constituent(ph)%label(ou),&
@@ -1405,6 +1402,7 @@ module subroutine mech_initializeRestorationPoints(ph,me)
 
   constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me)
   constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me)
+  constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me)  = constitutive_mech_F0(ph)%data(1:3,1:3,me)
   constitutive_mech_partitionedLi0(ph)%data(1:3,1:3,me) = constitutive_mech_Li0(ph)%data(1:3,1:3,me)
   constitutive_mech_partitionedLp0(ph)%data(1:3,1:3,me) = constitutive_mech_Lp0(ph)%data(1:3,1:3,me)
   constitutive_mech_partitionedS0(ph)%data(1:3,1:3,me)  = constitutive_mech_S0(ph)%data(1:3,1:3,me)
@@ -1424,6 +1422,7 @@ module subroutine constitutive_mech_windForward(ph,me)
 
   constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp(ph)%data(1:3,1:3,me)
   constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi(ph)%data(1:3,1:3,me)
+  constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me)  = constitutive_mech_F(ph)%data(1:3,1:3,me)
   constitutive_mech_partitionedLi0(ph)%data(1:3,1:3,me) = constitutive_mech_Li(ph)%data(1:3,1:3,me)
   constitutive_mech_partitionedLp0(ph)%data(1:3,1:3,me) = constitutive_mech_Lp(ph)%data(1:3,1:3,me)
   constitutive_mech_partitionedS0(ph)%data(1:3,1:3,me)  = constitutive_mech_S(ph)%data(1:3,1:3,me)
@@ -1445,6 +1444,7 @@ module subroutine constitutive_mech_forward()
   do ph = 1, size(plasticState)
     constitutive_mech_Fi0(ph) = constitutive_mech_Fi(ph)
     constitutive_mech_Fp0(ph) = constitutive_mech_Fp(ph)
+    constitutive_mech_F0(ph)  = constitutive_mech_F(ph)
     constitutive_mech_Li0(ph) = constitutive_mech_Li(ph)
     constitutive_mech_Lp0(ph) = constitutive_mech_Lp(ph)
     constitutive_mech_S0(ph)  = constitutive_mech_S(ph)
@@ -1519,7 +1519,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
   enddo
   subFp0 = constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me)
   subFi0 = constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me)
-  subF0  = crystallite_partitionedF0(1:3,1:3,co,ip,el)
+  subF0  = constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me)
   subFrac = 0.0_pReal
   subStep = 1.0_pReal/num%subStepSizeCryst
   todo = .true.
@@ -1569,7 +1569,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
 !  prepare for integration
     if (todo) then
       subF = subF0 &
-           + subStep * (crystallite_F(1:3,1:3,co,ip,el) - crystallite_partitionedF0(1:3,1:3,co,ip,el))
+           + subStep * (constitutive_mech_F(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me))
       constitutive_mech_Fe(ph)%data(1:3,1:3,me) = matmul(subF,math_inv33(matmul(constitutive_mech_Fi(ph)%data(1:3,1:3,me), &
                                                                                 constitutive_mech_Fp(ph)%data(1:3,1:3,me))))
       converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el)

src/constitutive_plastic_nonlocal.f90

@@ -552,10 +552,8 @@ end function plastic_nonlocal_init
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates quantities characterizing the microstructure
 !--------------------------------------------------------------------------------------------------
-module subroutine plastic_nonlocal_dependentState(F, instance, of, ip, el)
+module subroutine plastic_nonlocal_dependentState(instance, of, ip, el)
 
-  real(pReal), dimension(3,3), intent(in) :: &
-    F
   integer, intent(in) :: &
     instance, &
     of, &
@@ -647,7 +645,7 @@ module subroutine plastic_nonlocal_dependentState(F, instance, of, ip, el)
     ph = material_phaseAt(1,el)
     me = material_phaseMemberAt(1,ip,el)
     invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me))
-    invFe = matmul(constitutive_mech_Fp(ph)%data(1:3,1:3,me),math_inv33(F))
+    invFe = math_inv33(constitutive_mech_Fe(ph)%data(1:3,1:3,me))
 
     rho_edg_delta = rho0(:,mob_edg_pos) - rho0(:,mob_edg_neg)
     rho_scr_delta = rho0(:,mob_scr_pos) - rho0(:,mob_scr_neg)
@@ -976,13 +974,11 @@ end subroutine plastic_nonlocal_deltaState
 !---------------------------------------------------------------------------------------------------
 !> @brief calculates the rate of change of microstructure
 !---------------------------------------------------------------------------------------------------
-module subroutine plastic_nonlocal_dotState(Mp, F, Temperature,timestep, &
+module subroutine plastic_nonlocal_dotState(Mp, Temperature,timestep, &
                                             instance,of,ip,el)
 
   real(pReal), dimension(3,3), intent(in) :: &
     Mp                                                                                              !< MandelStress
-  real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: &
-    F                                                                                               !< Deformation gradient
   real(pReal), intent(in) :: &
     Temperature, &                                                                                  !< temperature
     timestep                                                                                        !< substepped crystallite time increment
@@ -1149,7 +1145,7 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Temperature,timestep, &
                                          - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) &
                                                                   - rhoDotSingle2DipoleGlide(s,9))  ! make sure that we do not annihilate more dipoles than we have
 
-  rhoDot = rhoDotFlux(F,timestep,  instance,of,ip,el) &
+  rhoDot = rhoDotFlux(timestep, instance,of,ip,el) &
          + rhoDotMultiplication &
          + rhoDotSingle2DipoleGlide &
          + rhoDotAthermalAnnihilation &
@@ -1178,10 +1174,8 @@ end subroutine plastic_nonlocal_dotState
 !---------------------------------------------------------------------------------------------------
 !> @brief calculates the rate of change of microstructure
 !---------------------------------------------------------------------------------------------------
-function rhoDotFlux(F,timestep,  instance,of,ip,el)
+function rhoDotFlux(timestep,instance,of,ip,el)
 
-  real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: &
-    F                                                                                               !< Deformation gradient
   real(pReal), intent(in) :: &
     timestep                                                                                        !< substepped crystallite time increment
   integer, intent(in) :: &
@@ -1293,7 +1287,7 @@ function rhoDotFlux(F,timestep,  instance,of,ip,el)
     m(1:3,:,3) = -prm%slip_transverse
     m(1:3,:,4) =  prm%slip_transverse
 
-    my_F = F(1:3,1:3,1,ip,el)
+    my_F = constitutive_mech_F(ph)%data(1:3,1:3,of)
     my_Fe = matmul(my_F, math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,of)))
 
     neighbors: do n = 1,nIPneighbors
@@ -1311,7 +1305,7 @@ function rhoDotFlux(F,timestep,  instance,of,ip,el)
 
       if (neighbor_n > 0) then                                                                      ! if neighbor exists, average deformation gradient
         neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el))
-        neighbor_F = F(1:3,1:3,1,neighbor_ip,neighbor_el)
+        neighbor_F = constitutive_mech_F(np)%data(1:3,1:3,no)
         neighbor_Fe = matmul(neighbor_F, math_inv33(constitutive_mech_Fp(np)%data(1:3,1:3,no)))
         Favg = 0.5_pReal * (my_F + neighbor_F)
       else                                                                                          ! if no neighbor, take my value as average

src/homogenization_mech.f90

@@ -202,15 +202,17 @@ module function mech_updateState(subdt,subF,ip,el) result(doneAndHappy)
 
   integer :: co
   real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
+  real(pReal) :: Fs(3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
 
 
   if (homogenization_type(material_homogenizationAt(el)) == HOMOGENIZATION_RGC_ID) then
       do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
         dPdFs(:,:,:,:,co) = crystallite_stressTangent(subdt,co,ip,el)
+        Fs(:,:,co)        = constitutive_mech_getF(co,ip,el)
       enddo
       doneAndHappy = &
           mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
-                        crystallite_F(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
+                               Fs, &
                                subF,&
                                subdt, &
                                dPdFs, &

src/thermal_conduction.f90

@@ -112,14 +112,16 @@ function thermal_conduction_getConductivity(ip,el)
     el                                                                                              !< element number
   real(pReal), dimension(3,3) :: &
     thermal_conduction_getConductivity
+  
   integer :: &
-    grain
+    co
 
 
   thermal_conduction_getConductivity = 0.0_pReal
-  do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el))
+  
+  do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
     thermal_conduction_getConductivity = thermal_conduction_getConductivity + &
-     crystallite_push33ToRef(grain,ip,el,lattice_K(:,:,material_phaseAt(grain,el)))
+     crystallite_push33ToRef(co,ip,el,lattice_K(:,:,material_phaseAt(co,el)))
   enddo
 
   thermal_conduction_getConductivity = thermal_conduction_getConductivity &
@@ -138,14 +140,16 @@ function thermal_conduction_getSpecificHeat(ip,el)
     el                                                                                              !< element number
   real(pReal) :: &
     thermal_conduction_getSpecificHeat
+
   integer :: &
-    grain
+    co
+
 
   thermal_conduction_getSpecificHeat = 0.0_pReal
 
-  do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el))
+  do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
     thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat &
-                                       + lattice_c_p(material_phaseAt(grain,el))
+                                       + lattice_c_p(material_phaseAt(co,el))
   enddo
 
   thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat &
@@ -164,15 +168,16 @@ function thermal_conduction_getMassDensity(ip,el)
     el                                                                                              !< element number
   real(pReal) :: &
     thermal_conduction_getMassDensity
+  
   integer :: &
-    grain
+    co
+
 
   thermal_conduction_getMassDensity = 0.0_pReal
 
-
+  do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
-  do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el))
     thermal_conduction_getMassDensity = thermal_conduction_getMassDensity &
-                                      + lattice_rho(material_phaseAt(grain,el))
+                                      + lattice_rho(material_phaseAt(co,el))
   enddo
 
   thermal_conduction_getMassDensity = thermal_conduction_getMassDensity &