diff --git a/PRIVATE b/PRIVATE index 277b4a693..76bb51348 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 277b4a693a57c1e4583bcf8ea2205f9b5d147198 +Subproject commit 76bb51348de75207d483d369628670e5ae51dca9 diff --git a/src/homogenization_mechanical_RGC.f90 b/src/homogenization_mechanical_RGC.f90 index 71d1542c8..76676726a 100644 --- a/src/homogenization_mechanical_RGC.f90 +++ b/src/homogenization_mechanical_RGC.f90 @@ -643,16 +643,16 @@ module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) !------------------------------------------------------------------------------------------------- !> @brief compute the equivalent shear and bulk moduli from the elasticity tensor !------------------------------------------------------------------------------------------------- - real(pReal) function equivalentMu(grainID,ce) + real(pReal) function equivalentMu(co,ce) integer, intent(in) :: & - grainID,& + co,& ce real(pReal), dimension(6,6) :: C - C = phase_homogenizedC(material_phaseID(grainID,ce),material_phaseEntry(grainID,ce)) + C = phase_homogenizedC66(material_phaseID(co,ce),material_phaseEntry(co,ce)) ! damage not included! equivalentMu = lattice_equivalent_mu(C,'voigt') end function equivalentMu diff --git a/src/lattice.f90 b/src/lattice.f90 index 11d1dca0f..0e24d1b18 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -405,7 +405,7 @@ module lattice contains !-------------------------------------------------------------------------------------------------- -!> @brief Module initialization +!> @brief module initialization !-------------------------------------------------------------------------------------------------- subroutine lattice_init @@ -417,7 +417,7 @@ end subroutine lattice_init !-------------------------------------------------------------------------------------------------- -!> @brief Characteristic shear for twinning +!> @brief characteristic shear for twinning !-------------------------------------------------------------------------------------------------- function lattice_characteristicShear_Twin(Ntwin,lattice,CoverA) result(characteristicShear) @@ -491,7 +491,7 @@ end function lattice_characteristicShear_Twin !-------------------------------------------------------------------------------------------------- -!> @brief Rotated elasticity matrices for twinning in 66-vector notation +!> @brief rotated elasticity matrices for twinning in 6x6-matrix notation !-------------------------------------------------------------------------------------------------- function lattice_C66_twin(Ntwin,C66,lattice,CoverA) @@ -505,6 +505,7 @@ function lattice_C66_twin(Ntwin,C66,lattice,CoverA) type(rotation) :: R integer :: i + select case(lattice) case('cF') coordinateSystem = buildCoordinateSystem(Ntwin,FCC_NSLIPSYSTEM,FCC_SYSTEMTWIN,& @@ -521,14 +522,14 @@ function lattice_C66_twin(Ntwin,C66,lattice,CoverA) do i = 1, sum(Ntwin) call R%fromAxisAngle([coordinateSystem(1:3,2,i),PI],P=1) ! ToDo: Why always 180 deg? - lattice_C66_twin(1:6,1:6,i) = R%rotTensor4sym(C66) + lattice_C66_twin(1:6,1:6,i) = math_3333toVoigt66(R%rotTensor4(math_Voigt66to3333(C66))) enddo end function lattice_C66_twin !-------------------------------------------------------------------------------------------------- -!> @brief Rotated elasticity matrices for transformation in 66-vector notation +!> @brief rotated elasticity matrices for transformation in 6x6-matrix notation !-------------------------------------------------------------------------------------------------- function lattice_C66_trans(Ntrans,C_parent66,lattice_target, & cOverA_trans,a_bcc,a_fcc) @@ -571,23 +572,23 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, & call IO_error(137,ext_msg='lattice_C66_trans : '//trim(lattice_target)) endif - do i = 1, 6 + do i = 1,6 if (abs(C_target_unrotated66(i,i)) @brief Non-schmid projections for bcc with up to 6 coefficients +!> @brief non-Schmid projections for bcc with up to 6 coefficients ! Koester et al. 2012, Acta Materialia 60 (2012) 3894–3901, eq. (17) ! Gröger et al. 2008, Acta Materialia 56 (2008) 5412–5425, table 1 !-------------------------------------------------------------------------------------------------- @@ -634,7 +635,7 @@ end function lattice_nonSchmidMatrix !-------------------------------------------------------------------------------------------------- -!> @brief Slip-slip interaction matrix +!> @brief slip-slip interaction matrix !> details only active slip systems are considered !-------------------------------------------------------------------------------------------------- function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(interactionMatrix) @@ -882,7 +883,7 @@ end function lattice_interaction_SlipBySlip !-------------------------------------------------------------------------------------------------- -!> @brief Twin-twin interaction matrix +!> @brief twin-twin interaction matrix !> details only active twin systems are considered !-------------------------------------------------------------------------------------------------- function lattice_interaction_TwinByTwin(Ntwin,interactionValues,lattice) result(interactionMatrix) @@ -980,7 +981,7 @@ end function lattice_interaction_TwinByTwin !-------------------------------------------------------------------------------------------------- -!> @brief Trans-trans interaction matrix +!> @brief trans-trans interaction matrix !> details only active trans systems are considered !-------------------------------------------------------------------------------------------------- function lattice_interaction_TransByTrans(Ntrans,interactionValues,lattice) result(interactionMatrix) @@ -1009,7 +1010,7 @@ function lattice_interaction_TransByTrans(Ntrans,interactionValues,lattice) resu 2,2,2,2,2,2,2,2,2,1,1,1 & ],shape(FCC_INTERACTIONTRANSTRANS)) !< Trans-trans interaction types for fcc - if(lattice == 'cF') then + if (lattice == 'cF') then interactionTypes = FCC_INTERACTIONTRANSTRANS NtransMax = FCC_NTRANSSYSTEM else @@ -1022,7 +1023,7 @@ end function lattice_interaction_TransByTrans !-------------------------------------------------------------------------------------------------- -!> @brief Slip-twin interaction matrix +!> @brief slip-twin interaction matrix !> details only active slip and twin systems are considered !-------------------------------------------------------------------------------------------------- function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) result(interactionMatrix) @@ -1185,7 +1186,7 @@ end function lattice_interaction_SlipByTwin !-------------------------------------------------------------------------------------------------- -!> @brief Slip-trans interaction matrix +!> @brief slip-trans interaction matrix !> details only active slip and trans systems are considered !-------------------------------------------------------------------------------------------------- function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,lattice) result(interactionMatrix) @@ -1238,7 +1239,7 @@ function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,lattice) !-------------------------------------------------------------------------------------------------- -!> @brief Twin-slip interaction matrix +!> @brief twin-slip interaction matrix !> details only active twin and slip systems are considered !-------------------------------------------------------------------------------------------------- function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,lattice) result(interactionMatrix) @@ -1410,7 +1411,7 @@ end function lattice_SchmidMatrix_twin !-------------------------------------------------------------------------------------------------- -!> @brief Schmid matrix for twinning +!> @brief Schmid matrix for transformation !> details only active twin systems are considered !-------------------------------------------------------------------------------------------------- function lattice_SchmidMatrix_trans(Ntrans,lattice_target,cOverA,a_bcc,a_fcc) result(SchmidMatrix) @@ -1482,7 +1483,7 @@ end function lattice_SchmidMatrix_cleavage !-------------------------------------------------------------------------------------------------- -!> @brief Slip direction of slip systems (|| b) +!> @brief slip direction of slip systems (|| b) !-------------------------------------------------------------------------------------------------- function lattice_slip_direction(Nslip,lattice,cOverA) result(d) @@ -1500,7 +1501,7 @@ end function lattice_slip_direction !-------------------------------------------------------------------------------------------------- -!> @brief Normal direction of slip systems (|| n) +!> @brief normal direction of slip systems (|| n) !-------------------------------------------------------------------------------------------------- function lattice_slip_normal(Nslip,lattice,cOverA) result(n) @@ -1518,7 +1519,7 @@ end function lattice_slip_normal !-------------------------------------------------------------------------------------------------- -!> @brief Transverse direction of slip systems (|| t = b x n) +!> @brief transverse direction of slip systems (|| t = b x n) !-------------------------------------------------------------------------------------------------- function lattice_slip_transverse(Nslip,lattice,cOverA) result(t) @@ -1536,7 +1537,7 @@ end function lattice_slip_transverse !-------------------------------------------------------------------------------------------------- -!> @brief Labels for slip systems +!> @brief labels of slip systems !> details only active slip systems are considered !-------------------------------------------------------------------------------------------------- function lattice_labels_slip(Nslip,lattice) result(labels) @@ -1577,7 +1578,7 @@ end function lattice_labels_slip !-------------------------------------------------------------------------------------------------- -!> @brief Return 3x3 tensor with symmetry according to given Bravais lattice +!> @brief return 3x3 tensor with symmetry according to given Bravais lattice !-------------------------------------------------------------------------------------------------- pure function lattice_symmetrize_33(T,lattice) result(T_sym) @@ -1604,7 +1605,7 @@ end function lattice_symmetrize_33 !-------------------------------------------------------------------------------------------------- -!> @brief Return stiffness matrix in 6x6 notation with symmetry according to given Bravais lattice +!> @brief return stiffness matrix in 6x6 notation with symmetry according to given Bravais lattice !> @details J. A. Rayne and B. S. Chandrasekhar Phys. Rev. 120, 1658 Erratum Phys. Rev. 122, 1962 !-------------------------------------------------------------------------------------------------- pure function lattice_symmetrize_C66(C66,lattice) result(C66_sym) @@ -1650,7 +1651,7 @@ end function lattice_symmetrize_C66 !-------------------------------------------------------------------------------------------------- -!> @brief Labels for twin systems +!> @brief labels of twin systems !> details only active twin systems are considered !-------------------------------------------------------------------------------------------------- function lattice_labels_twin(Ntwin,lattice) result(labels) @@ -1688,7 +1689,7 @@ end function lattice_labels_twin !-------------------------------------------------------------------------------------------------- -!> @brief Projection of the transverse direction onto the slip plane +!> @brief projection of the transverse direction onto the slip plane !> @details: This projection is used to calculate forest hardening for edge dislocations !-------------------------------------------------------------------------------------------------- function slipProjection_transverse(Nslip,lattice,cOverA) result(projection) @@ -1712,7 +1713,7 @@ end function slipProjection_transverse !-------------------------------------------------------------------------------------------------- -!> @brief Projection of the slip direction onto the slip plane +!> @brief projection of the slip direction onto the slip plane !> @details: This projection is used to calculate forest hardening for screw dislocations !-------------------------------------------------------------------------------------------------- function slipProjection_direction(Nslip,lattice,cOverA) result(projection) @@ -1778,7 +1779,7 @@ end function coordinateSystem_slip !-------------------------------------------------------------------------------------------------- -!> @brief Populate reduced interaction matrix +!> @brief populate reduced interaction matrix !-------------------------------------------------------------------------------------------------- function buildInteraction(reacting_used,acting_used,reacting_max,acting_max,values,matrix) @@ -1821,7 +1822,7 @@ end function buildInteraction !-------------------------------------------------------------------------------------------------- -!> @brief Build a local coordinate system on slip, twin, trans, cleavage systems +!> @brief build a local coordinate system on slip, twin, trans, cleavage systems !> @details Order: Direction, plane (normal), and common perpendicular !-------------------------------------------------------------------------------------------------- function buildCoordinateSystem(active,potential,system,lattice,cOverA) @@ -1888,7 +1889,7 @@ end function buildCoordinateSystem !-------------------------------------------------------------------------------------------------- -!> @brief Helper function to define transformation systems +!> @brief helper function to define transformation systems ! Needed to calculate Schmid matrix and rotated stiffness matrices. ! @details: set c/a = 0.0 for fcc -> bcc transformation ! set a_Xcc = 0.0 for fcc -> hex transformation @@ -2072,7 +2073,7 @@ end function getlabels !-------------------------------------------------------------------------------------------------- -!> @brief Equivalent Poisson's ratio (ν) +!> @brief equivalent Poisson's ratio (ν) !> @details https://doi.org/10.1143/JPSJ.20.635 !-------------------------------------------------------------------------------------------------- function lattice_equivalent_nu(C,assumption) result(nu) @@ -2086,12 +2087,12 @@ function lattice_equivalent_nu(C,assumption) result(nu) real(pReal), dimension(6,6) :: S - if (IO_lc(assumption) == 'voigt') then + if (IO_lc(assumption) == 'voigt') then K = (C(1,1)+C(2,2)+C(3,3) +2.0_pReal*(C(1,2)+C(2,3)+C(1,3))) & / 9.0_pReal - elseif(IO_lc(assumption) == 'reuss') then + elseif (IO_lc(assumption) == 'reuss') then call math_invert(S,error,C) - if(error) error stop 'matrix inversion failed' + if (error) error stop 'matrix inversion failed' K = 1.0_pReal & / (S(1,1)+S(2,2)+S(3,3) +2.0_pReal*(S(1,2)+S(2,3)+S(1,3))) else @@ -2099,13 +2100,13 @@ function lattice_equivalent_nu(C,assumption) result(nu) endif mu = lattice_equivalent_mu(C,assumption) - nu = (1.5_pReal*K -mu)/(3.0_pReal*K+mu) + nu = (1.5_pReal*K-mu)/(3.0_pReal*K+mu) end function lattice_equivalent_nu !-------------------------------------------------------------------------------------------------- -!> @brief Equivalent shear modulus (μ) +!> @brief equivalent shear modulus (μ) !> @details https://doi.org/10.1143/JPSJ.20.635 !-------------------------------------------------------------------------------------------------- function lattice_equivalent_mu(C,assumption) result(mu) @@ -2118,12 +2119,12 @@ function lattice_equivalent_mu(C,assumption) result(mu) real(pReal), dimension(6,6) :: S - if (IO_lc(assumption) == 'voigt') then + if (IO_lc(assumption) == 'voigt') then mu = (1.0_pReal*(C(1,1)+C(2,2)+C(3,3)) -1.0_pReal*(C(1,2)+C(2,3)+C(1,3)) +3.0_pReal*(C(4,4)+C(5,5)+C(6,6))) & / 15.0_pReal - elseif(IO_lc(assumption) == 'reuss') then + elseif (IO_lc(assumption) == 'reuss') then call math_invert(S,error,C) - if(error) error stop 'matrix inversion failed' + if (error) error stop 'matrix inversion failed' mu = 15.0_pReal & / (4.0_pReal*(S(1,1)+S(2,2)+S(3,3)) -4.0_pReal*(S(1,2)+S(2,3)+S(1,3)) +3.0_pReal*(S(4,4)+S(5,5)+S(6,6))) else @@ -2134,7 +2135,7 @@ end function lattice_equivalent_mu !-------------------------------------------------------------------------------------------------- -!> @brief Check correctness of some lattice functions. +!> @brief check correctness of some lattice functions !-------------------------------------------------------------------------------------------------- subroutine selfTest @@ -2152,7 +2153,7 @@ subroutine selfTest system = reshape([1.0_pReal+r(1),0.0_pReal,0.0_pReal, 0.0_pReal,1.0_pReal+r(2),0.0_pReal],[6,1]) CoSy = buildCoordinateSystem([1],[1],system,'cF',0.0_pReal) - if(any(dNeq(CoSy(1:3,1:3,1),math_I3))) error stop 'buildCoordinateSystem' + if (any(dNeq(CoSy(1:3,1:3,1),math_I3))) error stop 'buildCoordinateSystem' do i = 1, 10 call random_number(C) @@ -2198,13 +2199,13 @@ subroutine selfTest C(1,1) = C(1,1) + C(1,2) + 0.1_pReal C(4,4) = 0.5_pReal * (C(1,1) - C(1,2)) C = lattice_symmetrize_C66(C,'cI') - if(dNeq(C(4,4),lattice_equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/voigt' - if(dNeq(C(4,4),lattice_equivalent_mu(C,'reuss'),1.0e-12_pReal)) error stop 'equivalent_mu/reuss' + if (dNeq(C(4,4),lattice_equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/voigt' + if (dNeq(C(4,4),lattice_equivalent_mu(C,'reuss'),1.0e-12_pReal)) error stop 'equivalent_mu/reuss' lambda = C(1,2) - if(dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'voigt')), & + if (dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'voigt')), & lattice_equivalent_nu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_nu/voigt' - if(dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'reuss')), & + if (dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'reuss')), & lattice_equivalent_nu(C,'reuss'),1.0e-12_pReal)) error stop 'equivalent_nu/reuss' end subroutine selfTest diff --git a/src/math.f90 b/src/math.f90 index 31bd3b952..5d2609237 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -15,15 +15,15 @@ module math implicit none public #if __INTEL_COMPILER >= 1900 - ! do not make use associated entities available to other modules + ! do not make use of associated entities available to other modules private :: & IO, & config #endif real(pReal), parameter :: PI = acos(-1.0_pReal) !< ratio of a circle's circumference to its diameter - real(pReal), parameter :: INDEG = 180.0_pReal/PI !< conversion from radian into degree - real(pReal), parameter :: INRAD = PI/180.0_pReal !< conversion from degree into radian + real(pReal), parameter :: INDEG = 180.0_pReal/PI !< conversion from radian to degree + real(pReal), parameter :: INRAD = PI/180.0_pReal !< conversion from degree to radian complex(pReal), parameter :: TWOPIIMG = cmplx(0.0_pReal,2.0_pReal*PI) !< Re(0.0), Im(2xPi) real(pReal), dimension(3,3), parameter :: & @@ -132,19 +132,19 @@ pure recursive subroutine math_sort(a, istart, iend, sortDim) integer, intent(in),optional :: istart,iend, sortDim integer :: ipivot,s,e,d - if(present(istart)) then + if (present(istart)) then s = istart else s = lbound(a,2) endif - if(present(iend)) then + if (present(iend)) then e = iend else e = ubound(a,2) endif - if(present(sortDim)) then + if (present(sortDim)) then d = sortDim else d = 1 @@ -467,7 +467,7 @@ pure function math_inv33(A) logical :: error call math_invert33(math_inv33,DetA,error,A) - if(error) math_inv33 = 0.0_pReal + if (error) math_inv33 = 0.0_pReal end function math_inv33 @@ -698,7 +698,7 @@ pure function math_sym33to6(m33,weighted) integer :: i - if(present(weighted)) then + if (present(weighted)) then w = merge(NRMMANDEL,1.0_pReal,weighted) else w = NRMMANDEL @@ -725,7 +725,7 @@ pure function math_6toSym33(v6,weighted) integer :: i - if(present(weighted)) then + if (present(weighted)) then w = merge(INVNRMMANDEL,1.0_pReal,weighted) else w = INVNRMMANDEL @@ -802,7 +802,7 @@ pure function math_sym3333to66(m3333,weighted) integer :: i,j - if(present(weighted)) then + if (present(weighted)) then w = merge(NRMMANDEL,1.0_pReal,weighted) else w = NRMMANDEL @@ -822,7 +822,7 @@ end function math_sym3333to66 !-------------------------------------------------------------------------------------------------- -!> @brief convert 66 matrix into symmetric 3x3x3x3 matrix +!> @brief convert 6x6 matrix into symmetric 3x3x3x3 matrix !> @details Weighted conversion (default) rearranges according to Nye and weights shear ! components according to Mandel. Advisable for matrix operations. ! Unweighted conversion only rearranges order according to Nye @@ -837,7 +837,7 @@ pure function math_66toSym3333(m66,weighted) integer :: i,j - if(present(weighted)) then + if (present(weighted)) then w = merge(INVNRMMANDEL,1.0_pReal,weighted) else w = INVNRMMANDEL @@ -854,12 +854,13 @@ end function math_66toSym3333 !-------------------------------------------------------------------------------------------------- -!> @brief convert 66 Voigt matrix into symmetric 3x3x3x3 matrix +!> @brief convert 6x6 Voigt matrix into symmetric 3x3x3x3 matrix !-------------------------------------------------------------------------------------------------- pure function math_Voigt66to3333(m66) real(pReal), dimension(3,3,3,3) :: math_Voigt66to3333 real(pReal), dimension(6,6), intent(in) :: m66 !< 6x6 matrix + integer :: i,j @@ -873,6 +874,31 @@ pure function math_Voigt66to3333(m66) end function math_Voigt66to3333 +!-------------------------------------------------------------------------------------------------- +!> @brief convert symmetric 3x3x3x3 matrix into 6x6 Voigt matrix +!-------------------------------------------------------------------------------------------------- +pure function math_3333toVoigt66(m3333) + + real(pReal), dimension(6,6) :: math_3333toVoigt66 + real(pReal), dimension(3,3,3,3), intent(in) :: m3333 !< symmetric 3x3x3x3 matrix (no internal check) + + integer :: i,j + + +#ifndef __INTEL_COMPILER + do concurrent(i=1:6, j=1:6) + math_3333toVoigt66(i,j) = m3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) + end do +#else + do i=1,6; do j=1,6 + math_3333toVoigt66(i,j) = m3333(MAPVOIGT(1,i),MAPVOIGT(2,i),MAPVOIGT(1,j),MAPVOIGT(2,j)) + end do; end do +#endif + +end function math_3333toVoigt66 + + + !-------------------------------------------------------------------------------------------------- !> @brief draw a random sample from Gauss variable !-------------------------------------------------------------------------------------------------- @@ -958,7 +984,7 @@ subroutine math_eigh33(w,v,m) v(2,2) + m(2, 3) * w(1), & (m(1,1) - w(1)) * (m(2,2) - w(1)) - v(3,2)] norm = norm2(v(1:3, 1)) - fallback1: if(norm < threshold) then + fallback1: if (norm < threshold) then call math_eigh(w,v,error,m) else fallback1 v(1:3,1) = v(1:3, 1) / norm @@ -966,7 +992,7 @@ subroutine math_eigh33(w,v,m) v(2,2) + m(2, 3) * w(2), & (m(1,1) - w(2)) * (m(2,2) - w(2)) - v(3,2)] norm = norm2(v(1:3, 2)) - fallback2: if(norm < threshold) then + fallback2: if (norm < threshold) then call math_eigh(w,v,error,m) else fallback2 v(1:3,2) = v(1:3, 2) / norm @@ -1003,7 +1029,7 @@ pure function math_rotationalPart(F) result(R) I_F = [math_trace33(F), 0.5*(math_trace33(F)**2 - math_trace33(matmul(F,F)))] x = math_clip(I_C(1)**2 -3.0_pReal*I_C(2),0.0_pReal)**(3.0_pReal/2.0_pReal) - if(dNeq0(x)) then + if (dNeq0(x)) then Phi = acos(math_clip((I_C(1)**3 -4.5_pReal*I_C(1)*I_C(2) +13.5_pReal*I_C(3))/x,-1.0_pReal,1.0_pReal)) lambda = I_C(1) +(2.0_pReal * sqrt(math_clip(I_C(1)**2-3.0_pReal*I_C(2),0.0_pReal))) & *cos((Phi-2.0_pReal * PI*[1.0_pReal,2.0_pReal,3.0_pReal])/3.0_pReal) @@ -1065,7 +1091,7 @@ function math_eigvalsh33(m) - 2.0_pReal/27.0_pReal*I(1)**3.0_pReal & - I(3) ! different from http://arxiv.org/abs/physics/0610206 (this formulation was in DAMASK) - if(all(abs([P,Q]) < TOL)) then + if (all(abs([P,Q]) < TOL)) then math_eigvalsh33 = math_eigvalsh(m) else rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal @@ -1188,7 +1214,7 @@ real(pReal) pure elemental function math_clip(a, left, right) if (present(left)) math_clip = max(left,math_clip) if (present(right)) math_clip = min(right,math_clip) if (present(left) .and. present(right)) then - if(left>right) error stop 'left > right' + if (left>right) error stop 'left > right' endif end function math_clip @@ -1234,35 +1260,38 @@ subroutine selfTest error stop 'math_expand [1,2] by [1,2,3] => [1,2,2,1,1,1]' call math_sort(sort_in_,1,3,2) - if(any(sort_in_ /= sort_out_)) & + if (any(sort_in_ /= sort_out_)) & error stop 'math_sort' - if(any(math_range(5) /= range_out_)) & + if (any(math_range(5) /= range_out_)) & error stop 'math_range' - if(any(dNeq(math_exp33(math_I3,0),math_I3))) & + if (any(dNeq(math_exp33(math_I3,0),math_I3))) & error stop 'math_exp33(math_I3,1)' - if(any(dNeq(math_exp33(math_I3,128),exp(1.0_pReal)*math_I3))) & + if (any(dNeq(math_exp33(math_I3,128),exp(1.0_pReal)*math_I3))) & error stop 'math_exp33(math_I3,128)' call random_number(v9) - if(any(dNeq(math_33to9(math_9to33(v9)),v9))) & + if (any(dNeq(math_33to9(math_9to33(v9)),v9))) & error stop 'math_33to9/math_9to33' call random_number(t99) - if(any(dNeq(math_3333to99(math_99to3333(t99)),t99))) & + if (any(dNeq(math_3333to99(math_99to3333(t99)),t99))) & error stop 'math_3333to99/math_99to3333' call random_number(v6) - if(any(dNeq(math_sym33to6(math_6toSym33(v6)),v6))) & + if (any(dNeq(math_sym33to6(math_6toSym33(v6)),v6))) & error stop 'math_sym33to6/math_6toSym33' call random_number(t66) - if(any(dNeq(math_sym3333to66(math_66toSym3333(t66)),t66,1.0e-15_pReal))) & + if (any(dNeq(math_sym3333to66(math_66toSym3333(t66)),t66,1.0e-15_pReal))) & error stop 'math_sym3333to66/math_66toSym3333' + if (any(dNeq(math_3333toVoigt66(math_Voigt66to3333(t66)),t66,1.0e-15_pReal))) & + error stop 'math_3333toVoigt66/math_Voigt66to3333' + call random_number(v6) - if(any(dNeq0(math_6toSym33(v6) - math_symmetric33(math_6toSym33(v6))))) & + if (any(dNeq0(math_6toSym33(v6) - math_symmetric33(math_6toSym33(v6))))) & error stop 'math_symmetric33' call random_number(v3_1) @@ -1270,34 +1299,34 @@ subroutine selfTest call random_number(v3_3) call random_number(v3_4) - if(dNeq(abs(dot_product(math_cross(v3_1-v3_4,v3_2-v3_4),v3_3-v3_4))/6.0, & + if (dNeq(abs(dot_product(math_cross(v3_1-v3_4,v3_2-v3_4),v3_3-v3_4))/6.0, & math_volTetrahedron(v3_1,v3_2,v3_3,v3_4),tol=1.0e-12_pReal)) & error stop 'math_volTetrahedron' call random_number(t33) - if(dNeq(math_det33(math_symmetric33(t33)),math_detSym33(math_symmetric33(t33)),tol=1.0e-12_pReal)) & + if (dNeq(math_det33(math_symmetric33(t33)),math_detSym33(math_symmetric33(t33)),tol=1.0e-12_pReal)) & error stop 'math_det33/math_detSym33' - if(any(dNeq(t33+transpose(t33),math_mul3333xx33(math_identity4th(),t33+transpose(t33))))) & + if (any(dNeq(t33+transpose(t33),math_mul3333xx33(math_identity4th(),t33+transpose(t33))))) & error stop 'math_mul3333xx33/math_identity4th' - if(any(dNeq0(math_eye(3),math_inv33(math_I3)))) & + if (any(dNeq0(math_eye(3),math_inv33(math_I3)))) & error stop 'math_inv33(math_I3)' do while(abs(math_det33(t33))<1.0e-9_pReal) call random_number(t33) enddo - if(any(dNeq0(matmul(t33,math_inv33(t33)) - math_eye(3),tol=1.0e-9_pReal))) & + if (any(dNeq0(matmul(t33,math_inv33(t33)) - math_eye(3),tol=1.0e-9_pReal))) & error stop 'math_inv33' call math_invert33(t33_2,det,e,t33) - if(any(dNeq0(matmul(t33,t33_2) - math_eye(3),tol=1.0e-9_pReal)) .or. e) & + if (any(dNeq0(matmul(t33,t33_2) - math_eye(3),tol=1.0e-9_pReal)) .or. e) & error stop 'math_invert33: T:T^-1 != I' - if(dNeq(det,math_det33(t33),tol=1.0e-12_pReal)) & + if (dNeq(det,math_det33(t33),tol=1.0e-12_pReal)) & error stop 'math_invert33 (determinant)' call math_invert(t33_2,e,t33) - if(any(dNeq0(matmul(t33,t33_2) - math_eye(3),tol=1.0e-9_pReal)) .or. e) & + if (any(dNeq0(matmul(t33,t33_2) - math_eye(3),tol=1.0e-9_pReal)) .or. e) & error stop 'math_invert t33' do while(math_det33(t33)<1.0e-2_pReal) ! O(det(F)) = 1 @@ -1305,7 +1334,7 @@ subroutine selfTest enddo t33_2 = math_rotationalPart(transpose(t33)) t33 = math_rotationalPart(t33) - if(any(dNeq0(matmul(t33_2,t33) - math_I3,tol=1.0e-10_pReal))) & + if (any(dNeq0(matmul(t33_2,t33) - math_I3,tol=1.0e-10_pReal))) & error stop 'math_rotationalPart' call random_number(r) @@ -1313,33 +1342,33 @@ subroutine selfTest txx = math_eye(d) allocate(txx_2(d,d)) call math_invert(txx_2,e,txx) - if(any(dNeq0(txx_2,txx) .or. e)) & + if (any(dNeq0(txx_2,txx) .or. e)) & error stop 'math_invert(txx)/math_eye' call math_invert(t99_2,e,t99) ! not sure how likely it is that we get a singular matrix - if(any(dNeq0(matmul(t99_2,t99)-math_eye(9),tol=1.0e-9_pReal)) .or. e) & + if (any(dNeq0(matmul(t99_2,t99)-math_eye(9),tol=1.0e-9_pReal)) .or. e) & error stop 'math_invert(t99)' - if(any(dNeq(math_clip([4.0_pReal,9.0_pReal],5.0_pReal,6.5_pReal),[5.0_pReal,6.5_pReal]))) & + if (any(dNeq(math_clip([4.0_pReal,9.0_pReal],5.0_pReal,6.5_pReal),[5.0_pReal,6.5_pReal]))) & error stop 'math_clip' - if(math_factorial(10) /= 3628800) & + if (math_factorial(10) /= 3628800) & error stop 'math_factorial' - if(math_binomial(49,6) /= 13983816) & + if (math_binomial(49,6) /= 13983816) & error stop 'math_binomial' - if(math_multinomial([1,2,3,4]) /= 12600) & + if (math_multinomial([1,2,3,4]) /= 12600) & error stop 'math_multinomial' ijk = cshift([1,2,3],int(r*1.0e2_pReal)) - if(dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),+1.0_pReal)) & + if (dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),+1.0_pReal)) & error stop 'math_LeviCivita(even)' ijk = cshift([3,2,1],int(r*2.0e2_pReal)) - if(dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),-1.0_pReal)) & + if (dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),-1.0_pReal)) & error stop 'math_LeviCivita(odd)' ijk = cshift([2,2,1],int(r*2.0e2_pReal)) - if(dNeq0(math_LeviCivita(ijk(1),ijk(2),ijk(3)))) & + if (dNeq0(math_LeviCivita(ijk(1),ijk(2),ijk(3)))) & error stop 'math_LeviCivita' end subroutine selfTest diff --git a/src/phase.f90 b/src/phase.f90 index a36fc26b9..22c35416b 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -230,15 +230,15 @@ module phase end function phase_mechanical_constitutive !ToDo: Merge all the stiffness functions - module function phase_homogenizedC(ph,en) result(C) + module function phase_homogenizedC66(ph,en) result(C) integer, intent(in) :: ph, en real(pReal), dimension(6,6) :: C - end function phase_homogenizedC - module function phase_damage_C(C_homogenized,ph,en) result(C) - real(pReal), dimension(3,3,3,3), intent(in) :: C_homogenized - integer, intent(in) :: ph,en - real(pReal), dimension(3,3,3,3) :: C - end function phase_damage_C + end function phase_homogenizedC66 + module function phase_damage_C66(C66,ph,en) result(C66_degraded) + real(pReal), dimension(6,6), intent(in) :: C66 + integer, intent(in) :: ph,en + real(pReal), dimension(6,6) :: C66_degraded + end function phase_damage_C66 module function phase_f_phi(phi,co,ce) result(f) integer, intent(in) :: ce,co @@ -299,7 +299,7 @@ module phase public :: & phase_init, & - phase_homogenizedC, & + phase_homogenizedC66, & phase_f_phi, & phase_f_T, & phase_K_phi, & diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index 49b9f9528..3b3ace8ab 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -139,6 +139,7 @@ module function phase_damage_constitutive(Delta_t,co,ip,el) result(converged_) integer :: & ph, en + ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip) @@ -150,20 +151,21 @@ end function phase_damage_constitutive !-------------------------------------------------------------------------------------------------- !> @brief returns the degraded/modified elasticity matrix !-------------------------------------------------------------------------------------------------- -module function phase_damage_C(C_homogenized,ph,en) result(C) +module function phase_damage_C66(C66,ph,en) result(C66_degraded) + + real(pReal), dimension(6,6), intent(in) :: C66 + integer, intent(in) :: ph,en + real(pReal), dimension(6,6) :: C66_degraded - real(pReal), dimension(3,3,3,3), intent(in) :: C_homogenized - integer, intent(in) :: ph,en - real(pReal), dimension(3,3,3,3) :: C damageType: select case (phase_damage(ph)) case (DAMAGE_ISOBRITTLE_ID) damageType - C = C_homogenized * damage_phi(ph,en)**2 + C66_degraded = C66 * damage_phi(ph,en)**2 case default damageType - C = C_homogenized + C66_degraded = C66 end select damageType -end function phase_damage_C +end function phase_damage_C66 !-------------------------------------------------------------------------------------------------- @@ -417,7 +419,7 @@ function phase_damage_deltaState(Fe, ph, en) result(broken) sourceType: select case (phase_damage(ph)) case (DAMAGE_ISOBRITTLE_ID) sourceType - call isobrittle_deltaState(phase_homogenizedC(ph,en), Fe, ph,en) + call isobrittle_deltaState(phase_homogenizedC66(ph,en), Fe, ph,en) broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,en))) if (.not. broken) then myOffset = damageState(ph)%offsetDeltaState diff --git a/src/phase_damage_isobrittle.f90 b/src/phase_damage_isobrittle.f90 index 25ce1d3e9..9d5b4f508 100644 --- a/src/phase_damage_isobrittle.f90 +++ b/src/phase_damage_isobrittle.f90 @@ -96,6 +96,7 @@ end function isobrittle_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state +! ToDo: Use Voigt directly !-------------------------------------------------------------------------------------------------- module subroutine isobrittle_deltaState(C, Fe, ph,en) @@ -109,13 +110,16 @@ module subroutine isobrittle_deltaState(C, Fe, ph,en) epsilon real(pReal) :: & r_W + real(pReal), dimension(6,6) :: & + C_sym + C_sym = math_sym3333to66(math_Voigt66to3333(C)) epsilon = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3) associate(prm => param(ph), stt => state(ph), dlt => deltaState(ph)) - r_W = 2.0_pReal*dot_product(epsilon,matmul(C,epsilon))/prm%W_crit + r_W = 2.0_pReal*dot_product(epsilon,matmul(C_sym,epsilon))/prm%W_crit dlt%r_W(en) = merge(r_W - stt%r_W(en), 0.0_pReal, r_W > stt%r_W(en)) end associate diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index 9c9a159a0..1b39f51ee 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -15,7 +15,7 @@ submodule(phase:mechanical) elastic contains !-------------------------------------------------------------------------------------------------- -!> @brief Initialize elasticity. +!> @brief initialize elasticity !-------------------------------------------------------------------------------------------------- module subroutine elastic_init(phases) @@ -62,52 +62,30 @@ end subroutine elastic_init !-------------------------------------------------------------------------------------------------- -!> @brief Return 6x6 elasticity tensor. -!-------------------------------------------------------------------------------------------------- -function get_C66(ph,en) - - integer, intent(in) :: & - ph, & - en - real(pReal), dimension(6,6) :: get_C66 - - - associate(prm => param(ph)) - get_C66 = 0.0_pReal - get_C66(1,1) = prm%C_11 - get_C66(1,2) = prm%C_12 - get_C66(4,4) = prm%C_44 - - if (any(phase_lattice(ph) == ['hP','tI'])) then - get_C66(1,3) = prm%C_13 - get_C66(3,3) = prm%C_33 - end if - - if (phase_lattice(ph) == 'tI') get_C66(6,6) = prm%C_66 - - get_C66 = lattice_symmetrize_C66(get_C66,phase_lattice(ph)) - - end associate - -end function get_C66 - - -!-------------------------------------------------------------------------------------------------- -!> @brief Return 6x6 elasticity tensor. +!> @brief return 6x6 elasticity tensor !-------------------------------------------------------------------------------------------------- module function elastic_C66(ph,en) result(C66) integer, intent(in) :: & ph, & en - real(pReal), dimension(6,6) :: & - C66 + real(pReal), dimension(6,6) :: C66 associate(prm => param(ph)) + C66 = 0.0_pReal + C66(1,1) = prm%C_11 + C66(1,2) = prm%C_12 + C66(4,4) = prm%C_44 - C66 = get_C66(ph,en) - C66 = math_sym3333to66(math_Voigt66to3333(C66)) ! Literature data is in Voigt notation + if (any(phase_lattice(ph) == ['hP','tI'])) then + C66(1,3) = prm%C_13 + C66(3,3) = prm%C_33 + end if + + if (phase_lattice(ph) == 'tI') C66(6,6) = prm%C_66 + + C66 = lattice_symmetrize_C66(C66,phase_lattice(ph)) end associate @@ -115,7 +93,7 @@ end function elastic_C66 !-------------------------------------------------------------------------------------------------- -!> @brief Return shear modulus. +!> @brief return shear modulus !-------------------------------------------------------------------------------------------------- module function elastic_mu(ph,en) result(mu) @@ -126,12 +104,13 @@ module function elastic_mu(ph,en) result(mu) mu - mu = lattice_equivalent_mu(get_C66(ph,en),'voigt') + mu = lattice_equivalent_mu(elastic_C66(ph,en),'voigt') end function elastic_mu + !-------------------------------------------------------------------------------------------------- -!> @brief Return Poisson ratio. +!> @brief return Poisson ratio !-------------------------------------------------------------------------------------------------- module function elastic_nu(ph,en) result(nu) @@ -142,15 +121,16 @@ module function elastic_nu(ph,en) result(nu) nu - nu = lattice_equivalent_nu(get_C66(ph,en),'voigt') + nu = lattice_equivalent_nu(elastic_C66(ph,en),'voigt') end function elastic_nu !-------------------------------------------------------------------------------------------------- -!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to +!> @brief return the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to !> the elastic and intermediate deformation gradients using Hooke's law +! ToDo: Use Voigt matrix directly !-------------------------------------------------------------------------------------------------- module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & Fe, Fi, ph, en) @@ -173,13 +153,12 @@ module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & i, j - C = math_66toSym3333(phase_homogenizedC(ph,en)) - C = phase_damage_C(C,ph,en) + C = math_Voigt66to3333(phase_damage_C66(phase_homogenizedC66(ph,en),ph,en)) E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration - do i =1, 3;do j=1,3 + do i =1,3; do j=1,3 dS_dFe(i,j,1:3,1:3) = matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn end do; end do @@ -188,10 +167,9 @@ end subroutine phase_hooke_SandItsTangents !-------------------------------------------------------------------------------------------------- -!> @brief returns the homogenized elasticity matrix -!> ToDo: homogenizedC66 would be more consistent +!> @brief Return the homogenized elasticity matrix. !-------------------------------------------------------------------------------------------------- -module function phase_homogenizedC(ph,en) result(C) +module function phase_homogenizedC66(ph,en) result(C) real(pReal), dimension(6,6) :: C integer, intent(in) :: ph, en @@ -204,7 +182,7 @@ module function phase_homogenizedC(ph,en) result(C) C = elastic_C66(ph,en) end select plasticType -end function phase_homogenizedC +end function phase_homogenizedC66 end submodule elastic diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index 62119c819..10522737c 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -150,7 +150,7 @@ module function plastic_dislotwin_init() result(myPlasticity) myPlasticity = plastic_active('dislotwin') - if(count(myPlasticity) == 0) return + if (count(myPlasticity) == 0) return print'(/,1x,a)', '<<<+- phase:mechanical:plastic:dislotwin init -+>>>' print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) @@ -173,7 +173,7 @@ module function plastic_dislotwin_init() result(myPlasticity) do ph = 1, phases%length - if(.not. myPlasticity(ph)) cycle + if (.not. myPlasticity(ph)) cycle associate(prm => param(ph), dot => dotState(ph), stt => state(ph), dst => dependentState(ph)) @@ -368,7 +368,7 @@ module function plastic_dislotwin_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! parameters required for several mechanisms and their interactions - if(prm%sum_N_sl + prm%sum_N_tw + prm%sum_N_tw > 0) & + if (prm%sum_N_sl + prm%sum_N_tw + prm%sum_N_tw > 0) & prm%D = pl%get_asFloat('D') if (prm%sum_N_tw + prm%sum_N_tr > 0) & @@ -425,7 +425,7 @@ module function plastic_dislotwin_init() result(myPlasticity) stt%gamma_sl=>plasticState(ph)%state(startIndex:endIndex,:) dot%gamma_sl=>plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) - if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' + if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tw @@ -474,34 +474,34 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC) ph, en real(pReal), dimension(6,6) :: & homogenizedC, & - C66 + C real(pReal), dimension(:,:,:), allocatable :: & C66_tw, & C66_tr - integer :: i real(pReal) :: f_unrotated - associate(prm => param(ph), stt => state(ph)) - C66 = elastic_C66(ph,en) + C = elastic_C66(ph,en) + + associate(prm => param(ph), stt => state(ph)) f_unrotated = 1.0_pReal & - sum(stt%f_tw(1:prm%sum_N_tw,en)) & - sum(stt%f_tr(1:prm%sum_N_tr,en)) - homogenizedC = f_unrotated * C66 + homogenizedC = f_unrotated * C twinActive: if (prm%sum_N_tw > 0) then - C66_tw = lattice_C66_twin(prm%N_tw,C66,phase_lattice(ph),phase_cOverA(ph)) + C66_tw = lattice_C66_twin(prm%N_tw,C,phase_lattice(ph),phase_cOverA(ph)) do i=1,prm%sum_N_tw homogenizedC = homogenizedC & + stt%f_tw(i,en)*C66_tw(1:6,1:6,i) end do - end if twinActive - + end if twinActive + transActive: if (prm%sum_N_tr > 0) then - C66_tr = lattice_C66_trans(prm%N_tr,C66,prm%lattice_tr,0.0_pReal,prm%a_cI,prm%a_cF) + C66_tr = lattice_C66_trans(prm%N_tr,C,prm%lattice_tr,0.0_pReal,prm%a_cI,prm%a_cF) do i=1,prm%sum_N_tr homogenizedC = homogenizedC & + stt%f_tr(i,en)*C66_tr(1:6,1:6,i) @@ -595,7 +595,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en) Lp = Lp * f_unrotated dLp_dMp = dLp_dMp * f_unrotated - shearBandingContribution: if(dNeq0(prm%v_sb)) then + shearBandingContribution: if (dNeq0(prm%v_sb)) then E_kB_T = prm%E_sb/(kB*T) call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design? @@ -846,7 +846,7 @@ module subroutine plastic_dislotwin_results(ph,group) 'threshold stress for twinning','Pa',prm%systems_tw) case('f_tr') - if(prm%sum_N_tr>0) call results_writeDataset(stt%f_tr,group,trim(prm%output(ou)), & + if (prm%sum_N_tr>0) call results_writeDataset(stt%f_tr,group,trim(prm%output(ou)), & 'martensite volume fraction','m³/m³') end select @@ -929,8 +929,8 @@ pure subroutine kinetics_sl(Mp,T,ph,en, & end associate - if(present(ddot_gamma_dtau_sl)) ddot_gamma_dtau_sl = ddot_gamma_dtau - if(present(tau_sl)) tau_sl = tau + if (present(ddot_gamma_dtau_sl)) ddot_gamma_dtau_sl = ddot_gamma_dtau + if (present(tau_sl)) tau_sl = tau end subroutine kinetics_sl @@ -1000,7 +1000,7 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,& end associate - if(present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw = ddot_gamma_dtau + if (present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw = ddot_gamma_dtau end subroutine kinetics_tw @@ -1069,7 +1069,7 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,& end associate - if(present(ddot_gamma_dtau_tr)) ddot_gamma_dtau_tr = ddot_gamma_dtau + if (present(ddot_gamma_dtau_tr)) ddot_gamma_dtau_tr = ddot_gamma_dtau end subroutine kinetics_tr diff --git a/src/rotations.f90 b/src/rotations.f90 index 68f9273a3..e73fb2782 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -75,7 +75,6 @@ module rotations procedure, public :: rotVector procedure, public :: rotTensor2 procedure, public :: rotTensor4 - procedure, public :: rotTensor4sym procedure, public :: misorientation procedure, public :: standardize end type rotation @@ -371,27 +370,6 @@ pure function rotTensor4(self,T,active) result(tRot) end function rotTensor4 -!--------------------------------------------------------------------------------------------------- -!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @brief rotate a symmetric rank-4 tensor stored as (6,6) passively (default) or actively -!! ToDo: Need to check active/passive !!! -!--------------------------------------------------------------------------------------------------- -pure function rotTensor4sym(self,T,active) result(tRot) - - real(pReal), dimension(6,6) :: tRot - class(rotation), intent(in) :: self - real(pReal), intent(in), dimension(6,6) :: T - logical, intent(in), optional :: active - - if (present(active)) then - tRot = math_sym3333to66(rotTensor4(self,math_66toSym3333(T),active)) - else - tRot = math_sym3333to66(rotTensor4(self,math_66toSym3333(T))) - endif - -end function rotTensor4sym - - !--------------------------------------------------------------------------------------------------- !> @brief misorientation !--------------------------------------------------------------------------------------------------- @@ -400,6 +378,7 @@ pure elemental function misorientation(self,other) type(rotation) :: misorientation class(rotation), intent(in) :: self, other + misorientation%q = multiply_quaternion(other%q, conjugate_quaternion(self%q)) end function misorientation