Merge branch 'initial-V_e' into 'development'

Draf: Initial v e

See merge request damask/DAMASK!589

PRIVATE

@@ -1 +1 @@
-Subproject commit 1766c855e54668d6225c9b22db536be7052605bc
+Subproject commit d2c229e6dd81fb399b2c1b0658b750f7d309f175

python/damask/_configmaterial.py

@@ -482,7 +482,7 @@ class ConfigMaterial(Config):
         """
         N,n,shaped = 1,1,{}
 
-        map_dim = {'O':-1,'F_i':-2}
+        map_dim = {'O':-1,'V_e':-2}
         for k,v in kwargs.items():
             shaped[k] = np.array(v)
             s = shaped[k].shape[:map_dim.get(k,None)]
@@ -494,12 +494,12 @@ class ConfigMaterial(Config):
         if 'v' not in kwargs:
             shaped['v'] = np.broadcast_to(1/n,(N,n))
 
-        map_shape = {'O':(N,n,4),'F_i':(N,n,3,3)}
+        map_shape = {'O':(N,n,4),'V_e':(N,n,3,3)}
         for k,v in shaped.items():
             target = map_shape.get(k,(N,n))
             obj = np.broadcast_to(v.reshape(util.shapeshifter(v.shape, target, mode = 'right')), target)
             for i in range(N):
-                if k in ['phase','O','v','F_i']:
+                if k in ['phase','O','v','V_e']:
                     for j in range(n):
                         mat[i]['constituents'][j][k] = obj[i,j].item() if isinstance(obj[i,j],np.generic) else obj[i,j]
                 else:

python/tests/test_ConfigMaterial.py

@@ -113,15 +113,15 @@ class TestConfigMaterial:
     @pytest.mark.parametrize('N,n,kw',[
                                         (1,1,{'phase':'Gold',
                                               'O':[1,0,0,0],
-                                              'F_i':np.eye(3),
+                                              'V_e':np.eye(3),
                                               'homogenization':'SX'}),
                                         (3,1,{'phase':'Gold',
                                               'O':Rotation.from_random(3),
-                                              'F_i':np.broadcast_to(np.eye(3),(3,3,3)),
+                                              'V_e':np.broadcast_to(np.eye(3),(3,3,3)),
                                               'homogenization':'SX'}),
                                         (2,3,{'phase':np.broadcast_to(['a','b','c'],(2,3)),
                                               'O':Rotation.from_random((2,3)),
-                                              'F_i':np.broadcast_to(np.eye(3),(2,3,3,3)),
+                                              'V_e':np.broadcast_to(np.eye(3),(2,3,3,3)),
                                               'homogenization':['SX','PX']}),
                                         ])
     def test_material_add(self,kw,N,n):

src/IO.f90

@@ -428,6 +428,8 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
       msg = 'too many systems requested'
     case (146)
       msg = 'number of values does not match'
+    case (147)
+      msg = 'V_e needs to be symmetric'
     case (148)
       msg = 'Nconstituents mismatch between homogenization and material'
 

src/material.f90

@@ -27,7 +27,7 @@ module material
 
 
   type(tRotationContainer), dimension(:), allocatable, public, protected :: material_O_0
-  type(tTensorContainer),   dimension(:), allocatable, public, protected :: material_F_i_0
+  type(tTensorContainer),   dimension(:), allocatable, public, protected :: material_V_e_0
 
   integer, dimension(:), allocatable, public, protected :: &
     homogenization_Nconstituents                                                                    !< number of grains in each homogenization
@@ -133,7 +133,7 @@ subroutine parse()
   allocate(material_v(homogenization_maxNconstituents,discretization_Ncells),source=0.0_pReal)
 
   allocate(material_O_0(materials%length))
-  allocate(material_F_i_0(materials%length))
+  allocate(material_V_e_0(materials%length))
 
   allocate(ho_of(materials%length))
   allocate(ph_of(materials%length,homogenization_maxNconstituents),source=-1)
@@ -154,7 +154,7 @@ subroutine parse()
         if (constituents%length /= homogenization%get_asInt('N_constituents')) call IO_error(148)
 
         allocate(material_O_0(ma)%data(constituents%length))
-        allocate(material_F_i_0(ma)%data(1:3,1:3,constituents%length))
+        allocate(material_V_e_0(ma)%data(1:3,1:3,constituents%length))
 
         do co = 1, constituents%length
           constituent => constituents%get(co)
@@ -162,7 +162,9 @@ subroutine parse()
           ph_of(ma,co) = phases%getIndex(constituent%get_asString('phase'))
 
           call material_O_0(ma)%data(co)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4))
-          material_F_i_0(ma)%data(1:3,1:3,co) = constituent%get_as2dFloat('F_i',defaultVal=math_I3,requiredShape=[3,3])
+          material_V_e_0(ma)%data(1:3,1:3,co) = constituent%get_as2dFloat('V_e',defaultVal=math_I3,requiredShape=[3,3])
+          if (any(dNeq(material_V_e_0(ma)%data(1:3,1:3,co),transpose(material_V_e_0(ma)%data(1:3,1:3,co))))) &
+            call IO_error(147)
 
         end do
         if (dNeq(sum(v_of(ma,:)),1.0_pReal,1.e-9_pReal)) call IO_error(153,ext_msg='constituent')

src/phase_mechanical.f90

@@ -210,7 +210,6 @@ module subroutine mechanical_init(phases)
     Nmembers
   class(tNode), pointer :: &
     num_crystallite, &
-
     phase, &
     mech
 
@@ -239,11 +238,8 @@ module subroutine mechanical_init(phases)
 
     allocate(phase_mechanical_Fe(ph)%data(3,3,Nmembers))
     allocate(phase_mechanical_Fi(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_Fi0(ph)%data(3,3,Nmembers))
     allocate(phase_mechanical_Fp(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_Fp0(ph)%data(3,3,Nmembers))
     allocate(phase_mechanical_F(ph)%data(3,3,Nmembers))
-    allocate(phase_mechanical_F0(ph)%data(3,3,Nmembers))
     allocate(phase_mechanical_Li(ph)%data(3,3,Nmembers),source=0.0_pReal)
     allocate(phase_mechanical_Li0(ph)%data(3,3,Nmembers),source=0.0_pReal)
     allocate(phase_mechanical_Lp(ph)%data(3,3,Nmembers),source=0.0_pReal)
@@ -265,23 +261,19 @@ module subroutine mechanical_init(phases)
     ma = discretization_materialAt((ce-1)/discretization_nIPs+1)
     do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
       ph = material_phaseID(co,ce)
-      phase_mechanical_Fi0(ph)%data(1:3,1:3,material_phaseEntry(co,ce)) = material_F_i_0(ma)%data(1:3,1:3,co)
+      en = material_phaseEntry(co,ce)
+      phase_mechanical_F(ph)%data(1:3,1:3,en)  = math_I3
+      phase_mechanical_Fp(ph)%data(1:3,1:3,en) = material_O_0(ma)%data(co)%asMatrix()              ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005)
+      phase_mechanical_Fe(ph)%data(1:3,1:3,en) = matmul(material_V_e_0(ma)%data(1:3,1:3,co), &
+                                                        transpose(phase_mechanical_Fp(ph)%data(1:3,1:3,en)))
+      phase_mechanical_Fi(ph)%data(1:3,1:3,en) = material_O_0(ma)%data(co)%rotate(math_inv33(material_V_e_0(ma)%data(1:3,1:3,co)))
     enddo
   enddo
 
   do ph = 1, phases%length
-    do en = 1, count(material_phaseID == ph)
+    phase_mechanical_F0(ph)%data  = phase_mechanical_F(ph)%data
-
+    phase_mechanical_Fp0(ph)%data = phase_mechanical_Fp(ph)%data
-      phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = phase_O_0(ph)%data(en)%asMatrix()                 ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005)
+    phase_mechanical_Fi0(ph)%data = phase_mechanical_Fi(ph)%data
-      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_F0(ph)%data(1:3,1:3,en) = math_I3
-      phase_mechanical_Fe(ph)%data(1:3,1:3,en) = math_inv33(matmul(phase_mechanical_Fi0(ph)%data(1:3,1:3,en), &
-                                                                   phase_mechanical_Fp0(ph)%data(1:3,1:3,en)))  ! assuming that euler angles are given in internal strain free configuration
-    enddo
-    phase_mechanical_Fp(ph)%data = phase_mechanical_Fp0(ph)%data
-    phase_mechanical_Fi(ph)%data = phase_mechanical_Fi0(ph)%data
-    phase_mechanical_F(ph)%data  = phase_mechanical_F0(ph)%data
   enddo
 
 
@@ -318,7 +310,6 @@ module subroutine mechanical_init(phases)
 
   end select
 
-
   call eigen_init(phases)
 
 
@@ -1244,6 +1235,9 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
 end function phase_mechanical_dPdF
 
 
+!--------------------------------------------------------------------------------------------------
+!< @brief Write restart information to file.
+!--------------------------------------------------------------------------------------------------
 module subroutine mechanical_restartWrite(groupHandle,ph)
 
   integer(HID_T), intent(in) :: groupHandle
@@ -1251,36 +1245,41 @@ module subroutine mechanical_restartWrite(groupHandle,ph)
 
 
   call HDF5_write(plasticState(ph)%state,groupHandle,'omega_plastic')
-  call HDF5_write(phase_mechanical_Fi(ph)%data,groupHandle,'F_i')
-  call HDF5_write(phase_mechanical_Li(ph)%data,groupHandle,'L_i')
-  call HDF5_write(phase_mechanical_Lp(ph)%data,groupHandle,'L_p')
-  call HDF5_write(phase_mechanical_Fp(ph)%data,groupHandle,'F_p')
   call HDF5_write(phase_mechanical_S(ph)%data,groupHandle,'S')
   call HDF5_write(phase_mechanical_F(ph)%data,groupHandle,'F')
+  call HDF5_write(phase_mechanical_Fp(ph)%data,groupHandle,'F_p')
+  call HDF5_write(phase_mechanical_Fi(ph)%data,groupHandle,'F_i')
+  call HDF5_write(phase_mechanical_Lp(ph)%data,groupHandle,'L_p')
+  call HDF5_write(phase_mechanical_Li(ph)%data,groupHandle,'L_i')
 
 end subroutine mechanical_restartWrite
 
 
+!--------------------------------------------------------------------------------------------------
+!< @brief Read restart information from file.
+!--------------------------------------------------------------------------------------------------
 module subroutine mechanical_restartRead(groupHandle,ph)
 
   integer(HID_T), intent(in) :: groupHandle
   integer, intent(in) :: ph
 
+  integer :: en
+
 
   call HDF5_read(plasticState(ph)%state0,groupHandle,'omega_plastic')
-  call HDF5_read(phase_mechanical_Fi0(ph)%data,groupHandle,'F_i')
-  call HDF5_read(phase_mechanical_Li0(ph)%data,groupHandle,'L_i')
-  call HDF5_read(phase_mechanical_Lp0(ph)%data,groupHandle,'L_p')
-  call HDF5_read(phase_mechanical_Fp0(ph)%data,groupHandle,'F_p')
   call HDF5_read(phase_mechanical_S0(ph)%data,groupHandle,'S')
   call HDF5_read(phase_mechanical_F0(ph)%data,groupHandle,'F')
+  call HDF5_read(phase_mechanical_Fp0(ph)%data,groupHandle,'F_p')
+  call HDF5_read(phase_mechanical_Fi0(ph)%data,groupHandle,'F_i')
+  call HDF5_read(phase_mechanical_Lp0(ph)%data,groupHandle,'L_p')
+  call HDF5_read(phase_mechanical_Li0(ph)%data,groupHandle,'L_i')
 
 end subroutine mechanical_restartRead
 
 
-!----------------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------
-!< @brief Get first Piola-Kichhoff stress (for use by non-mech physics)
+!< @brief Get first Piola-Kichhoff stress (for use by non-mech physics).
-!----------------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------
 module function mechanical_S(ph,en) result(S)
 
   integer, intent(in) :: ph,en
@@ -1292,9 +1291,9 @@ module function mechanical_S(ph,en) result(S)
 end function mechanical_S
 
 
-!----------------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------
-!< @brief Get plastic velocity gradient (for use by non-mech physics)
+!< @brief Get plastic velocity gradient (for use by non-mech physics).
-!----------------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------
 module function mechanical_L_p(ph,en) result(L_p)
 
   integer, intent(in) :: ph,en
@@ -1306,9 +1305,9 @@ module function mechanical_L_p(ph,en) result(L_p)
 end function mechanical_L_p
 
 
-!----------------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------
-!< @brief Get elastic deformation gradient (for use by non-mech physics)
+!< @brief Get elastic deformation gradient (for use by non-mech physics).
-!----------------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------
 module function mechanical_F_e(ph,en) result(F_e)
 
   integer, intent(in) :: ph,en
@@ -1320,9 +1319,9 @@ module function mechanical_F_e(ph,en) result(F_e)
 end function mechanical_F_e
 
 
-!----------------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------
-!< @brief Get second Piola-Kichhoff stress (for use by homogenization)
+!< @brief Get second Piola-Kichhoff stress (for use by homogenization).
-!----------------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------
 module function phase_P(co,ce) result(P)
 
   integer, intent(in) :: co, ce
@@ -1334,9 +1333,9 @@ module function phase_P(co,ce) result(P)
 end function phase_P
 
 
-!----------------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------
-!< @brief Get deformation gradient (for use by homogenization)
+!< @brief Get deformation gradient (for use by homogenization).
-!----------------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------
 module function phase_F(co,ce) result(F)
 
   integer, intent(in) :: co, ce
@@ -1348,9 +1347,9 @@ module function phase_F(co,ce) result(F)
 end function phase_F
 
 
-!----------------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------
-!< @brief Set deformation gradient (for use by homogenization)
+!< @brief Set deformation gradient (for use by homogenization).
-!----------------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------
 module subroutine phase_set_F(F,co,ce)
 
   real(pReal), dimension(3,3), intent(in) :: F

src/rotations.f90

@@ -323,17 +323,16 @@ pure function rotTensor2(self,T,active) result(tRot)
 
   logical           :: passive
 
+
   if (present(active)) then
     passive = .not. active
   else
     passive = .true.
   endif
 
-  if (passive) then
+  tRot = merge(matmul(matmul(self%asMatrix(),T),transpose(self%asMatrix())), &
-    tRot = matmul(matmul(self%asMatrix(),T),transpose(self%asMatrix()))
+               matmul(matmul(transpose(self%asMatrix()),T),self%asMatrix()), &
-  else
+               passive)
-    tRot = matmul(matmul(transpose(self%asMatrix()),T),self%asMatrix())
-  endif
 
 end function rotTensor2
 
@@ -450,6 +449,7 @@ pure function qu2om(qu) result(om)
   om(1,3) = 2.0_pReal*(qu(2)*qu(4)+qu(1)*qu(3))
 
   if (sign(1.0_pReal,P) < 0.0_pReal) om = transpose(om)
+  om = om/math_det33(om)**(1.0_pReal/3.0_pReal)
 
 end function qu2om
 
@@ -575,7 +575,7 @@ end function qu2cu
 
 !---------------------------------------------------------------------------------------------------
 !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief convert rotation matrix to cubochoric
+!> @brief convert rotation matrix to unit quaternion
 !> @details the original formulation (direct conversion) had (numerical?) issues
 !---------------------------------------------------------------------------------------------------
 pure function om2qu(om) result(qu)
@@ -602,14 +602,15 @@ pure function om2qu(om) result(qu)
       endif
   endif
   if(sign(1.0_pReal,qu(1))<0.0_pReal) qu =-1.0_pReal * qu
-  qu = qu*[1.0_pReal,P,P,P]
+  qu(2:4) = merge(qu(2:4),qu(2:4)*P,dEq0(qu(2:4)))
+  qu = qu/norm2(qu)
 
 end function om2qu
 
 
 !---------------------------------------------------------------------------------------------------
 !> @author Marc De Graef, Carnegie Mellon University
-!> @brief orientation matrix to Euler angles
+!> @brief convert orientation matrix to Euler angles
 !> @details Two step check for special cases to avoid invalid operations (not needed for python)
 !---------------------------------------------------------------------------------------------------
 pure function om2eu(om) result(eu)
@@ -1334,8 +1335,8 @@ pure function cu2ho(cu) result(ho)
       ! transform to sphere grid (inverse Lambert)
       ! [note that there is no need to worry about dividing by zero, since XYZ(3) can not become zero]
       c = sum(T**2)
-      s = PI * c/(24.0*XYZ(3)**2)
+      s = c * PI/(24.0*XYZ(3)**2)
-      c = sqrt(PI) * c / sqrt(24.0_pReal) / XYZ(3)
+      c = c * sqrt(PI/24.0_pReal) / XYZ(3)
       q = sqrt( 1.0 - s )
       LamXYZ = [ T(order(2)) * q, T(order(1)) * q, PREF * XYZ(3) - c ]
     end if special