natural data structure

src/discretization.f90

@@ -9,20 +9,20 @@ module discretization
 
   implicit none
   private
-  
+
   integer,     public, protected :: &
     discretization_nIPs, &
     discretization_Nelems
-    
-  integer,     public, protected, dimension(:),   allocatable :: &
-    discretization_materialAt                                                                       !ToDo: discretization_materialID
 
-  real(pReal), public, protected, dimension(:,:), allocatable :: & 
+  integer,     public, protected, dimension(:),   allocatable :: &
+    discretization_materialAt                                                                       !ToDo: discretization_ID_material
+
+  real(pReal), public, protected, dimension(:,:), allocatable :: &
     discretization_IPcoords0, &
     discretization_IPcoords, &
     discretization_NodeCoords0, &
     discretization_NodeCoords
-    
+
   integer :: &
     discretization_sharedNodesBegin
 
@@ -33,7 +33,7 @@ module discretization
     discretization_setNodeCoords
 
 contains
-  
+
 !--------------------------------------------------------------------------------------------------
 !> @brief stores the relevant information in globally accesible variables
 !--------------------------------------------------------------------------------------------------
@@ -54,20 +54,20 @@ subroutine discretization_init(materialAt,&
   discretization_Nelems = size(materialAt,1)
   discretization_nIPs   = size(IPcoords0,2)/discretization_Nelems
 
-  discretization_materialAt = materialAt  
+  discretization_materialAt = materialAt
 
   discretization_IPcoords0   = IPcoords0
   discretization_IPcoords    = IPcoords0
 
   discretization_NodeCoords0 = NodeCoords0
   discretization_NodeCoords  = NodeCoords0
-  
+
   if(present(sharedNodesBegin)) then
     discretization_sharedNodesBegin = sharedNodesBegin
   else
     discretization_sharedNodesBegin = size(discretization_NodeCoords0,2)
   endif
-  
+
 end subroutine discretization_init
 
 
@@ -77,13 +77,13 @@ end subroutine discretization_init
 subroutine discretization_results
 
   real(pReal), dimension(:,:), allocatable :: u
-  
+
   call results_closeGroup(results_addGroup('current/geometry'))
-  
+
   u = discretization_NodeCoords (1:3,:discretization_sharedNodesBegin) &
     - discretization_NodeCoords0(1:3,:discretization_sharedNodesBegin)
   call results_writeDataset('current/geometry',u,'u_n','displacements of the nodes','m')
-  
+
   u = discretization_IPcoords &
     - discretization_IPcoords0
   call results_writeDataset('current/geometry',u,'u_p','displacements of the materialpoints (cell centers)','m')
@@ -97,7 +97,7 @@ end subroutine discretization_results
 subroutine discretization_setIPcoords(IPcoords)
 
   real(pReal), dimension(:,:), intent(in) :: IPcoords
-  
+
   discretization_IPcoords = IPcoords
 
 end subroutine discretization_setIPcoords
@@ -109,7 +109,7 @@ end subroutine discretization_setIPcoords
 subroutine discretization_setNodeCoords(NodeCoords)
 
   real(pReal), dimension(:,:), intent(in) :: NodeCoords
-  
+
   discretization_NodeCoords = NodeCoords
 
 end subroutine discretization_setNodeCoords

src/material.f90

@@ -16,6 +16,12 @@ module material
   implicit none
   private
 
+  type :: tRotationContainer
+    type(Rotation), dimension(:),  allocatable :: data
+  end type
+
+  type(tRotationContainer), dimension(:), allocatable :: material_orientation0_2
+
   integer, dimension(:), allocatable, public, protected :: &
     homogenization_Nconstituents                                                                    !< number of grains in each homogenization
   integer, public, protected :: &
@@ -39,6 +45,8 @@ module material
     material_phaseMemberAt                                                                          !< position of the element within its phase instance
 
   public :: &
+    tRotationContainer, &
+    material_orientation0_2, &
     material_init
 
 contains
@@ -88,7 +96,7 @@ subroutine parse()
   real(pReal) :: &
     frac
   integer :: &
-    el, ip, co, &
+    el, ip, co, ma, &
     h, ce
 
   materials       => config_material%get('material')
@@ -153,6 +161,17 @@ subroutine parse()
 
   enddo
 
+  allocate(material_orientation0_2(materials%length))
+
+  do ma = 1, materials%length
+    material     => materials%get(ma)
+    constituents => material%get('constituents')
+    allocate(material_orientation0_2(ma)%data(constituents%length))
+    do co = 1, constituents%length
+      call material_orientation0_2(ma)%data(co)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4))
+    enddo
+ enddo
+
 end subroutine parse
 
 

src/phase.f90

@@ -19,8 +19,9 @@ module phase
   implicit none
   private
 
-  type(Rotation), dimension(:,:,:), allocatable :: &
+  type(tRotationContainer), dimension(:), allocatable :: &
-    material_orientation0                                                                           !< initial orientation of each grain,IP,element
+    phase_orientation0, &
+    phase_orientation
 
   type(rotation),            dimension(:,:,:),        allocatable :: &
     crystallite_orientation                                                                         !< current orientation
@@ -346,7 +347,7 @@ contains
 subroutine phase_init
 
   integer :: &
-    ph
+    ph, ce, co, ma
   class (tNode), pointer :: &
     debug_constitutive, &
     materials, &
@@ -367,6 +368,25 @@ subroutine phase_init
   materials => config_material%get('material')
   phases    => config_material%get('phase')
 
+
+  allocate(phase_orientation0(phases%length))
+  do ph = 1,phases%length
+    allocate(phase_orientation0(ph)%data(count(material_phaseID==ph)))
+  enddo
+
+  do ce = 1, size(material_phaseID,2)
+    ma = discretization_materialAt((ce-1)/discretization_nIPs+1)
+    do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
+      ph = material_phaseID(co,ce)
+      phase_orientation0(ph)%data(material_phaseEntry(co,ce)) = material_orientation0_2(ma)%data(co)
+    enddo
+  enddo
+
+  allocate(phase_orientation(phases%length))
+  do ph = 1,phases%length
+    phase_orientation(ph)%data = phase_orientation0(ph)%data
+  enddo
+
   call mechanical_init(materials,phases)
   call damage_init
   call thermal_init(phases)
@@ -602,7 +622,7 @@ function crystallite_push33ToRef(co,ce, tensor33)
 
   ph = material_phaseID(co,ce)
   en = material_phaseEntry(co,ce)
-  T = matmul(material_orientation0(co,ph,en)%asMatrix(),transpose(math_inv33(phase_F(co,ce)))) ! ToDo: initial orientation correct?
+  T = matmul(phase_orientation0(ph)%data(en)%asMatrix(),transpose(math_inv33(phase_F(co,ce))))      ! ToDo: initial orientation correct?
 
   crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T))
 

src/phase_mechanical.f90

@@ -68,7 +68,7 @@ submodule(phase) mechanical
         dS_dFe, &                                                                                   !< derivative of 2nd P-K stress with respect to elastic deformation gradient
         dS_dFi                                                                                      !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
     end subroutine phase_hooke_SandItsTangents
-    
+
     module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en)
       real(pReal), dimension(3,3),     intent(out) :: &
         Li                                                                                          !< inleastic velocity gradient
@@ -204,7 +204,6 @@ module subroutine mechanical_init(materials,phases)
     num_crystallite, &
     material, &
     constituents, &
-    constituent, &
     phase, &
     mech
 
@@ -228,8 +227,6 @@ module subroutine mechanical_init(materials,phases)
   allocate(phase_mechanical_P(phases%length))
   allocate(phase_mechanical_S0(phases%length))
 
-  allocate(material_orientation0(homogenization_maxNconstituents,phases%length,maxVal(material_phaseEntry)))
-
   do ph = 1, phases%length
     Nmembers = count(material_phaseID == ph)
 
@@ -260,15 +257,11 @@ module subroutine mechanical_init(materials,phases)
   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_phaseID(co,(el-1)*discretization_nIPs + ip)
       en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
 
-      call material_orientation0(co,ph,en)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4))
+      phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = phase_orientation0(ph)%data(en)%asMatrix()        ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005)
-
-      phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = material_orientation0(co,ph,en)%asMatrix()                         ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005)
       phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) &
                                                 / math_det33(phase_mechanical_Fp0(ph)%data(1:3,1:3,en))**(1.0_pReal/3.0_pReal)
       phase_mechanical_Fi0(ph)%data(1:3,1:3,en) = math_I3

src/phase_mechanical_plastic_nonlocal.f90

@@ -1450,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,ph,en)%asQuaternion(), &
+      if (any(dNeq(phase_orientation0(ph)%data(en)%asQuaternion(), &
-                   material_orientation0(1,neighbor_phase,neighbor_me)%asQuaternion())) .and. &
+                   phase_orientation0(neighbor_phase)%data(neighbor_me)%asQuaternion())) .and. &
           (.not. phase_localPlasticity(neighbor_phase))) &
         forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB)
     else