modifying storing of orientations;

needed for nonlocal, and also to remove the use if ip,el at homogenization level.
ip, el should be used only for looping eventually.

src/homogenization.f90

@@ -536,13 +536,15 @@ function damage_nonlocal_getDiffusion(ip,el)
     damage_nonlocal_getDiffusion
   integer :: &
     homog, &
-    grain
+    grain, &
+    ce
 
   homog  = material_homogenizationAt(el)
   damage_nonlocal_getDiffusion = 0.0_pReal
+  ce = (el-1)*discretization_nIPs + ip
   do grain = 1, homogenization_Nconstituents(homog)
     damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + &
-      crystallite_push33ToRef(grain,ip,el,lattice_D(1:3,1:3,material_phaseAt(grain,el)))
+      crystallite_push33ToRef(grain,ce,lattice_D(1:3,1:3,material_phaseAt2(grain,ce)))
   enddo
 
   damage_nonlocal_getDiffusion = &

src/homogenization_mechanical.f90

@@ -153,13 +153,13 @@ module subroutine mechanical_homogenize(dt,ip,el)
   chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
 
     case (HOMOGENIZATION_NONE_ID) chosenHomogenization
-        homogenization_P(1:3,1:3,ce)            = phase_mechanical_getP(1,ip,el)
+        homogenization_P(1:3,1:3,ce)            = phase_mechanical_getP(1,ce)
         homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(dt,1,ip,el)
 
     case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
       do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
         dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ip,el)
-        Ps(:,:,co) = phase_mechanical_getP(co,ip,el)
+        Ps(:,:,co) = phase_mechanical_getP(co,ce)
       enddo
       call mechanical_isostrain_averageStressAndItsTangent(&
         homogenization_P(1:3,1:3,ce), &
@@ -170,7 +170,7 @@ module subroutine mechanical_homogenize(dt,ip,el)
     case (HOMOGENIZATION_RGC_ID) chosenHomogenization
       do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
         dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ip,el)
-        Ps(:,:,co) = phase_mechanical_getP(co,ip,el)
+        Ps(:,:,co) = phase_mechanical_getP(co,ce)
       enddo
       call mechanical_RGC_averageStressAndItsTangent(&
         homogenization_P(1:3,1:3,ce), &
@@ -208,8 +208,8 @@ module function mechanical_updateState(subdt,subF,ce,ip,el) result(doneAndHappy)
   if (homogenization_type(material_homogenizationAt2(ce)) == HOMOGENIZATION_RGC_ID) then
       do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce))
         dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(subdt,co,ip,el)
-        Fs(:,:,co)        = phase_mechanical_getF(co,ip,el)
+        Fs(:,:,co)        = phase_mechanical_getF(co,ce)
-        Ps(:,:,co)        = phase_mechanical_getP(co,ip,el)
+        Ps(:,:,co)        = phase_mechanical_getP(co,ce)
       enddo
       doneAndHappy = mechanical_RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ce)
   else

src/homogenization_thermal.f90

@@ -107,16 +107,17 @@ module function thermal_conduction_getConductivity(ip,el) result(K)
   real(pReal), dimension(3,3) :: K
 
   integer :: &
-    co
+    co, &
-
+    ce
 
   K = 0.0_pReal
 
-  do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
+  ce = (el-1)*discretization_nIPs + ip
-    K = K + crystallite_push33ToRef(co,ip,el,lattice_K(:,:,material_phaseAt(co,el)))
+  do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce))
+    K = K + crystallite_push33ToRef(co,ce,lattice_K(:,:,material_phaseAt2(co,ce)))
   enddo
 
-  K = K / real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
+  K = K / real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal)
 
 end function thermal_conduction_getConductivity
 

src/material.f90

@@ -39,9 +39,6 @@ module material
   integer, dimension(:,:,:), allocatable, public, protected :: &                                    ! (constituent,IP,elem)
     material_phaseMemberAt                                                                          !< position of the element within its phase instance
 
-  type(Rotation), dimension(:,:,:), allocatable, public, protected :: &
-    material_orientation0                                                                           !< initial orientation of each grain,IP,element
-
   public :: &
     material_init
 
@@ -125,8 +122,6 @@ subroutine material_parseMaterial
   allocate(material_phaseAt2(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0)
   allocate(material_phaseMemberAt2(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0)
 
-  allocate(material_orientation0(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems))
-
   do el = 1, discretization_Nelems
     material     => materials%get(discretization_materialAt(el))
     constituents => material%get('constituents')
@@ -135,7 +130,7 @@ subroutine material_parseMaterial
     do ip = 1, discretization_nIPs
       ce = (el-1)*discretization_nIPs + ip
       counterHomogenization(material_homogenizationAt(el)) = counterHomogenization(material_homogenizationAt(el)) + 1
-      material_homogenizationMemberAt(ip,el)                = counterHomogenization(material_homogenizationAt(el))
+      material_homogenizationMemberAt(ip,el)               = counterHomogenization(material_homogenizationAt(el))
       material_homogenizationAt2(ce)       = material_homogenizationAt(el)
       material_homogenizationMemberAt2(ce) = material_homogenizationMemberAt(ip,el)
     enddo
@@ -153,7 +148,6 @@ subroutine material_parseMaterial
 
         material_phaseAt2(co,ce)       = material_phaseAt(co,el)
         material_phaseMemberAt2(co,ce) = material_phaseMemberAt(co,ip,el)
-        call material_orientation0(co,ip,el)%fromQuaternion(constituent%get_asFloats('O',requiredSize=4)) ! should be done in crystallite
       enddo
 
     enddo

src/phase.f90

@@ -19,6 +19,9 @@ module phase
   implicit none
   private
 
+  type(Rotation), dimension(:,:,:), allocatable, protected :: &
+    material_orientation0                                                                           !< initial orientation of each grain,IP,element
+
   type(rotation),            dimension(:,:,:),        allocatable :: &
     crystallite_orientation                                                                         !< current orientation
 
@@ -77,8 +80,8 @@ module phase
   interface
 
 ! == cleaned:begin =================================================================================
-    module subroutine mechanical_init(phases)
+    module subroutine mechanical_init(materials,phases)
-      class(tNode), pointer :: phases
+      class(tNode), pointer :: materials,phases
     end subroutine mechanical_init
 
     module subroutine damage_init
@@ -147,8 +150,8 @@ module phase
       real(pReal), dimension(3,3) :: L_p
     end function mechanical_L_p
 
-    module function phase_mechanical_getF(co,ip,el) result(F)
+    module function phase_mechanical_getF(co,ce) result(F)
-      integer, intent(in) :: co, ip, el
+      integer, intent(in) :: co, ce
       real(pReal), dimension(3,3) :: F
     end function phase_mechanical_getF
 
@@ -157,8 +160,8 @@ module phase
       real(pReal), dimension(3,3) :: F_e
     end function mechanical_F_e
 
-    module function phase_mechanical_getP(co,ip,el) result(P)
+    module function phase_mechanical_getP(co,ce) result(P)
-      integer, intent(in) :: co, ip, el
+      integer, intent(in) :: co, ce
       real(pReal), dimension(3,3) :: P
     end function phase_mechanical_getP
 
@@ -342,6 +345,7 @@ subroutine phase_init
     so                                                                                               !< counter in source loop
   class (tNode), pointer :: &
     debug_constitutive, &
+    materials, &
     phases
 
 
@@ -356,9 +360,10 @@ subroutine phase_init
   debugConstitutive%grain      =  config_debug%get_asInt('grain',defaultVal = 1)
 
 
+  materials => config_material%get('material')
   phases => config_material%get('phase')
 
-  call mechanical_init(phases)
+  call mechanical_init(materials,phases)
   call damage_init
   call thermal_init(phases)
 
@@ -624,19 +629,20 @@ end subroutine crystallite_orientations
 !--------------------------------------------------------------------------------------------------
 !> @brief Map 2nd order tensor to reference config
 !--------------------------------------------------------------------------------------------------
-function crystallite_push33ToRef(co,ip,el, tensor33)
+function crystallite_push33ToRef(co,ce, tensor33)
 
   real(pReal), dimension(3,3), intent(in) :: tensor33
   integer, intent(in):: &
-    el, &
+    co, &
-    ip, &
+    ce
-    co
   real(pReal), dimension(3,3) :: crystallite_push33ToRef
 
   real(pReal), dimension(3,3)             :: T
+  integer :: ph, me
 
-
+  ph = material_phaseAt2(co,ce)
-  T = matmul(material_orientation0(co,ip,el)%asMatrix(),transpose(math_inv33(phase_mechanical_getF(co,ip,el)))) ! ToDo: initial orientation correct?
+  me = material_phaseMemberAt2(co,ce)
+  T = matmul(material_orientation0(co,ph,me)%asMatrix(),transpose(math_inv33(phase_mechanical_getF(co,ce)))) ! ToDo: initial orientation correct?
 
   crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T))
 

src/phase_mechanical.f90

@@ -177,21 +177,26 @@ contains
 !> @brief Initialize mechanical field related constitutive models
 !> @details Initialize elasticity, plasticity and stiffness degradation models.
 !--------------------------------------------------------------------------------------------------
-module subroutine mechanical_init(phases)
+module subroutine mechanical_init(materials,phases)
 
   class(tNode), pointer :: &
+    materials, &
     phases
 
   integer :: &
     el, &
     ip, &
     co, &
+    ce, &
     ph, &
     me, &
     stiffDegradationCtr, &
     Nconstituents
   class(tNode), pointer :: &
     num_crystallite, &
+    material, &
+    constituents, &
+    constituent, &
     phase, &
     mech, &
     elastic, &
@@ -221,6 +226,8 @@ module subroutine mechanical_init(phases)
   allocate(phase_mechanical_P(phases%length))
   allocate(phase_mechanical_S0(phases%length))
 
+  allocate(material_orientation0(homogenization_maxNconstituents,phases%length,maxval(material_phaseMemberAt)))
+
   do ph = 1, phases%length
     Nconstituents = count(material_phaseAt == ph) * discretization_nIPs
 
@@ -271,14 +278,20 @@ module subroutine mechanical_init(phases)
     enddo
   endif
 
+  
   !$OMP PARALLEL DO PRIVATE(ph,me)
   do el = 1, size(material_phaseMemberAt,3); do ip = 1, size(material_phaseMemberAt,2)
     do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
+      material     => materials%get(discretization_materialAt(el))
+      constituents => material%get('constituents')
+      constituent => constituents%get(co)
 
       ph = material_phaseAt(co,el)
       me = material_phaseMemberAt(co,ip,el)
 
-      phase_mechanical_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)
+      call material_orientation0(co,ph,me)%fromQuaternion(constituent%get_asFloats('O',requiredSize=4))
+
+      phase_mechanical_Fp0(ph)%data(1:3,1:3,me) = material_orientation0(co,ph,me)%asMatrix()                         ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005)
       phase_mechanical_Fp0(ph)%data(1:3,1:3,me) = phase_mechanical_Fp0(ph)%data(1:3,1:3,me) &
                                                 / math_det33(phase_mechanical_Fp0(ph)%data(1:3,1:3,me))**(1.0_pReal/3.0_pReal)
       phase_mechanical_Fi0(ph)%data(1:3,1:3,me) = math_I3
@@ -1440,13 +1453,13 @@ end function mechanical_L_p
 !----------------------------------------------------------------------------------------------
 !< @brief Get deformation gradient (for use by homogenization)
 !----------------------------------------------------------------------------------------------
-module function phase_mechanical_getF(co,ip,el) result(F)
+module function phase_mechanical_getF(co,ce) result(F)
 
-  integer, intent(in) :: co, ip, el
+  integer, intent(in) :: co, ce
   real(pReal), dimension(3,3) :: F
 
 
-  F = phase_mechanical_F(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el))
+  F = phase_mechanical_F(material_phaseAt2(co,ce))%data(1:3,1:3,material_phaseMemberAt2(co,ce))
 
 end function phase_mechanical_getF
 
@@ -1469,13 +1482,13 @@ end function mechanical_F_e
 !----------------------------------------------------------------------------------------------
 !< @brief Get second Piola-Kichhoff stress (for use by homogenization)
 !----------------------------------------------------------------------------------------------
-module function phase_mechanical_getP(co,ip,el) result(P)
+module function phase_mechanical_getP(co,ce) result(P)
 
-  integer, intent(in) :: co, ip, el
+  integer, intent(in) :: co, ce
   real(pReal), dimension(3,3) :: P
 
 
-  P = phase_mechanical_P(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el))
+  P = phase_mechanical_P(material_phaseAt2(co,ce))%data(1:3,1:3,material_phaseMemberAt2(co,ce))
 
 end function phase_mechanical_getP
 

src/phase_mechanical_plastic_nonlocal.f90

@@ -1403,8 +1403,10 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
 
   integer :: &
     n, &                                                                                            ! neighbor index
+    me, &
     neighbor_e, &                                                                                   ! element index of my neighbor
     neighbor_i, &                                                                                   ! integration point index of my neighbor
+    neighbor_me, &                                                                                   
     neighbor_phase, &
     ns, &                                                                                           ! number of active slip systems
     s1, &                                                                                           ! slip system index (me)
@@ -1422,6 +1424,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
   associate(prm => param(ph))
   ns = prm%sum_N_sl
 
+  me = material_phaseMemberAt(1,i,e)
   !*** start out fully compatible
   my_compatibility = 0.0_pReal
   forall(s1 = 1:ns) my_compatibility(:,s1,s1,:) = 1.0_pReal
@@ -1429,7 +1432,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
   neighbors: do n = 1,nIPneighbors
     neighbor_e = IPneighborhood(1,n,i,e)
     neighbor_i = IPneighborhood(2,n,i,e)
-
+    neighbor_me = material_phaseMemberAt(1,neighbor_i,neighbor_e)
     neighbor_phase = material_phaseAt(1,neighbor_e)
 
     if (neighbor_e <= 0 .or. neighbor_i <= 0) then
@@ -1447,8 +1450,8 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
     elseif (prm%chi_GB >= 0.0_pReal) then
       !* GRAIN BOUNDARY !
       !* fixed transmissivity for adjacent ips with different texture (only if explicitly given in material.config)
-      if (any(dNeq(material_orientation0(1,i,e)%asQuaternion(), &
+      if (any(dNeq(material_orientation0(1,ph,me)%asQuaternion(), &
-                   material_orientation0(1,neighbor_i,neighbor_e)%asQuaternion())) .and. &
+                   material_orientation0(1,neighbor_phase,neighbor_me)%asQuaternion())) .and. &
           (.not. phase_localPlasticity(neighbor_phase))) &
         forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB)
     else