diff --git a/PRIVATE b/PRIVATE index 1766c855e..d2c229e6d 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 1766c855e54668d6225c9b22db536be7052605bc +Subproject commit d2c229e6dd81fb399b2c1b0658b750f7d309f175 diff --git a/python/damask/_configmaterial.py b/python/damask/_configmaterial.py index 725567b7f..7a06a2c69 100644 --- a/python/damask/_configmaterial.py +++ b/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: diff --git a/python/tests/test_ConfigMaterial.py b/python/tests/test_ConfigMaterial.py index d3ea986c5..f003e254e 100644 --- a/python/tests/test_ConfigMaterial.py +++ b/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): diff --git a/src/IO.f90 b/src/IO.f90 index 240c9e992..047d11add 100644 --- a/src/IO.f90 +++ b/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' diff --git a/src/material.f90 b/src/material.f90 index 2258c40a1..eafc411cd 100644 --- a/src/material.f90 +++ b/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') diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index fc9e40e02..88b86a8d9 100644 --- a/src/phase_mechanical.f90 +++ b/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_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_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 + phase_mechanical_F0(ph)%data = phase_mechanical_F(ph)%data + phase_mechanical_Fp0(ph)%data = phase_mechanical_Fp(ph)%data + phase_mechanical_Fi0(ph)%data = phase_mechanical_Fi(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 diff --git a/src/rotations.f90 b/src/rotations.f90 index 7aa81fcab..5367f2c07 100644 --- a/src/rotations.f90 +++ b/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 = matmul(matmul(self%asMatrix(),T),transpose(self%asMatrix())) - else - tRot = matmul(matmul(transpose(self%asMatrix()),T),self%asMatrix()) - endif + tRot = merge(matmul(matmul(self%asMatrix(),T),transpose(self%asMatrix())), & + matmul(matmul(transpose(self%asMatrix()),T),self%asMatrix()), & + passive) 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) - c = sqrt(PI) * c / sqrt(24.0_pReal) / XYZ(3) + s = c * PI/(24.0*XYZ(3)**2) + 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