polish
This commit is contained in:
Philip Eisenlohr
parent
021d614daf
commit
da23c916ca
|
|
@ -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)
|
||||
|
||||
|
|
@ -529,7 +529,7 @@ 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)
|
||||
|
|
@ -588,7 +588,7 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
|
|||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @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
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
@ -635,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)
|
||||
|
|
@ -883,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)
|
||||
|
|
@ -981,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)
|
||||
|
|
@ -1010,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
|
||||
|
|
@ -1023,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)
|
||||
|
|
@ -1186,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)
|
||||
|
|
@ -1239,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)
|
||||
|
|
@ -1411,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)
|
||||
|
|
@ -1483,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)
|
||||
|
||||
|
|
@ -1501,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)
|
||||
|
||||
|
|
@ -1519,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)
|
||||
|
||||
|
|
@ -1537,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)
|
||||
|
|
@ -1578,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)
|
||||
|
||||
|
|
@ -1605,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)
|
||||
|
|
@ -1651,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)
|
||||
|
|
@ -1689,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)
|
||||
|
|
@ -1713,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)
|
||||
|
|
@ -1779,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)
|
||||
|
||||
|
|
@ -1822,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)
|
||||
|
|
@ -1889,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
|
||||
|
|
@ -2073,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)
|
||||
|
|
@ -2087,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
|
||||
|
|
@ -2100,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)
|
||||
|
|
@ -2119,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
|
||||
|
|
@ -2135,7 +2135,7 @@ end function lattice_equivalent_mu
|
|||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Check correctness of some lattice functions.
|
||||
!> @brief check correctness of some lattice functions
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine selfTest
|
||||
|
||||
|
|
@ -2153,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)
|
||||
|
|
@ -2199,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
|
||||
|
|
|
|||
94
src/math.f90
94
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,7 +854,7 @@ end function math_66toSym3333
|
|||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Convert 6x6 Voigt matrix into symmetric 3x3x3x3 matrix.
|
||||
!> @brief convert 6x6 Voigt matrix into symmetric 3x3x3x3 matrix
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
pure function math_Voigt66to3333(m66)
|
||||
|
||||
|
|
@ -875,7 +875,7 @@ end function math_Voigt66to3333
|
|||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Convert symmetric 3x3x3x3 matrix into 6x6 Voigt matrix.
|
||||
!> @brief convert symmetric 3x3x3x3 matrix into 6x6 Voigt matrix
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
pure function math_3333toVoigt66(m3333)
|
||||
|
||||
|
|
@ -984,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
|
||||
|
|
@ -992,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
|
||||
|
|
@ -1029,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)
|
||||
|
|
@ -1091,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
|
||||
|
|
@ -1214,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
|
||||
|
|
@ -1260,38 +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))) &
|
||||
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)
|
||||
|
|
@ -1299,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
|
||||
|
|
@ -1334,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)
|
||||
|
|
@ -1342,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
|
||||
|
|
|
|||
|
|
@ -15,7 +15,7 @@ submodule(phase:mechanical) elastic
|
|||
contains
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Initialize elasticity.
|
||||
!> @brief initialize elasticity
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
module subroutine elastic_init(phases)
|
||||
|
||||
|
|
@ -62,7 +62,7 @@ end subroutine elastic_init
|
|||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Return 6x6 elasticity tensor.
|
||||
!> @brief return 6x6 elasticity tensor
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
module function elastic_C66(ph,en) result(C66)
|
||||
|
||||
|
|
@ -93,7 +93,7 @@ end function elastic_C66
|
|||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Return shear modulus.
|
||||
!> @brief return shear modulus
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
module function elastic_mu(ph,en) result(mu)
|
||||
|
||||
|
|
@ -110,7 +110,7 @@ end function elastic_mu
|
|||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Return Poisson ratio.
|
||||
!> @brief return Poisson ratio
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
module function elastic_nu(ph,en) result(nu)
|
||||
|
||||
|
|
@ -128,7 +128,7 @@ 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
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
@ -473,9 +473,7 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC)
|
|||
integer, intent(in) :: &
|
||||
ph, en
|
||||
real(pReal), dimension(6,6) :: &
|
||||
homogenizedC
|
||||
|
||||
real(pReal), dimension(6,6) :: &
|
||||
homogenizedC, &
|
||||
C
|
||||
real(pReal), dimension(:,:,:), allocatable :: &
|
||||
C66_tw, &
|
||||
|
|
@ -500,8 +498,8 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC)
|
|||
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,C,prm%lattice_tr,0.0_pReal,prm%a_cI,prm%a_cF)
|
||||
do i=1,prm%sum_N_tr
|
||||
|
|
@ -597,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?
|
||||
|
|
@ -848,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
|
||||
|
|
@ -931,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
|
||||
|
||||
|
|
@ -1002,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
|
||||
|
||||
|
|
@ -1071,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
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue