223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge branch 'use-Voigt-notation' into 'development' 

use Voigt notation

See merge request damask/DAMASK!458
This commit is contained in:
Philip Eisenlohr 2021-11-23 17:49:47 +00:00
parent 1396ec6f56 da23c916ca
commit bf76d9f3a7
10 changed files with 193 additions and 200 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 277b4a693a57c1e4583bcf8ea2205f9b5d147198 Subproject commit 76bb51348de75207d483d369628670e5ae51dca9

6
src/homogenization_mechanical_RGC.f90
View File

@ -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 !> @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) :: & integer, intent(in) :: &
grainID,& co,&
ce ce
real(pReal), dimension(6,6) :: C 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') equivalentMu = lattice_equivalent_mu(C,'voigt')
end function equivalentMu end function equivalentMu

85
src/lattice.f90
View File

@ -405,7 +405,7 @@ module lattice
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Module initialization !> @brief module initialization
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine lattice_init 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) 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) function lattice_C66_twin(Ntwin,C66,lattice,CoverA)
@ -505,6 +505,7 @@ function lattice_C66_twin(Ntwin,C66,lattice,CoverA)
type(rotation) :: R type(rotation) :: R
integer :: i integer :: i
select case(lattice) select case(lattice)
case('cF') case('cF')
coordinateSystem = buildCoordinateSystem(Ntwin,FCC_NSLIPSYSTEM,FCC_SYSTEMTWIN,& coordinateSystem = buildCoordinateSystem(Ntwin,FCC_NSLIPSYSTEM,FCC_SYSTEMTWIN,&
@ -521,14 +522,14 @@ function lattice_C66_twin(Ntwin,C66,lattice,CoverA)
do i = 1, sum(Ntwin) do i = 1, sum(Ntwin)
call R%fromAxisAngle([coordinateSystem(1:3,2,i),PI],P=1) ! ToDo: Why always 180 deg? 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 enddo
end function lattice_C66_twin 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, & function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
cOverA_trans,a_bcc,a_fcc) 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)) call IO_error(137,ext_msg='lattice_C66_trans : '//trim(lattice_target))
endif endif
do i = 1, 6 do i = 1,6
if (abs(C_target_unrotated66(i,i))<tol_math_check) & if (abs(C_target_unrotated66(i,i))<tol_math_check) &
call IO_error(135,el=i,ext_msg='matrix diagonal "el"ement in transformation') call IO_error(135,el=i,ext_msg='matrix diagonal "el"ement in transformation')
enddo enddo
call buildTransformationSystem(Q,S,Ntrans,cOverA_trans,a_fcc,a_bcc) call buildTransformationSystem(Q,S,Ntrans,cOverA_trans,a_fcc,a_bcc)
do i = 1, sum(Ntrans) do i = 1,sum(Ntrans)
call R%fromMatrix(Q(1:3,1:3,i)) call R%fromMatrix(Q(1:3,1:3,i))
lattice_C66_trans(1:6,1:6,i) = R%rotTensor4sym(C_target_unrotated66) lattice_C66_trans(1:6,1:6,i) = math_3333toVoigt66(R%rotTensor4(math_Voigt66to3333(C_target_unrotated66)))
enddo enddo
end function lattice_C66_trans end function lattice_C66_trans
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @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) ! 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 ! 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 !> details only active slip systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(interactionMatrix) 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 !> details only active twin systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_TwinByTwin(Ntwin,interactionValues,lattice) result(interactionMatrix) 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 !> details only active trans systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_TransByTrans(Ntrans,interactionValues,lattice) result(interactionMatrix) 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 & 2,2,2,2,2,2,2,2,2,1,1,1 &
],shape(FCC_INTERACTIONTRANSTRANS)) !< Trans-trans interaction types for fcc ],shape(FCC_INTERACTIONTRANSTRANS)) !< Trans-trans interaction types for fcc
if(lattice == 'cF') then if (lattice == 'cF') then
interactionTypes = FCC_INTERACTIONTRANSTRANS interactionTypes = FCC_INTERACTIONTRANSTRANS
NtransMax = FCC_NTRANSSYSTEM NtransMax = FCC_NTRANSSYSTEM
else 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 !> details only active slip and twin systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) result(interactionMatrix) 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 !> details only active slip and trans systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,lattice) result(interactionMatrix) 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 !> details only active twin and slip systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,lattice) result(interactionMatrix) 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 !> details only active twin systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_SchmidMatrix_trans(Ntrans,lattice_target,cOverA,a_bcc,a_fcc) result(SchmidMatrix) 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) 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) 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) 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 !> details only active slip systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_labels_slip(Nslip,lattice) result(labels) 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) 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 !> @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) 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 !> details only active twin systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_labels_twin(Ntwin,lattice) result(labels) 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 !> @details: This projection is used to calculate forest hardening for edge dislocations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function slipProjection_transverse(Nslip,lattice,cOverA) result(projection) 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 !> @details: This projection is used to calculate forest hardening for screw dislocations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function slipProjection_direction(Nslip,lattice,cOverA) result(projection) 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) 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 !> @details Order: Direction, plane (normal), and common perpendicular
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function buildCoordinateSystem(active,potential,system,lattice,cOverA) 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. ! Needed to calculate Schmid matrix and rotated stiffness matrices.
! @details: set c/a = 0.0 for fcc -> bcc transformation ! @details: set c/a = 0.0 for fcc -> bcc transformation
! set a_Xcc = 0.0 for fcc -> hex 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 !> @details https://doi.org/10.1143/JPSJ.20.635
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_equivalent_nu(C,assumption) result(nu) function lattice_equivalent_nu(C,assumption) result(nu)
@ -2089,9 +2090,9 @@ function lattice_equivalent_nu(C,assumption) result(nu)
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))) & 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 / 9.0_pReal
elseif(IO_lc(assumption) == 'reuss') then elseif (IO_lc(assumption) == 'reuss') then
call math_invert(S,error,C) call math_invert(S,error,C)
if(error) error stop 'matrix inversion failed' if (error) error stop 'matrix inversion failed'
K = 1.0_pReal & 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))) / (S(1,1)+S(2,2)+S(3,3) +2.0_pReal*(S(1,2)+S(2,3)+S(1,3)))
else else
@ -2099,13 +2100,13 @@ function lattice_equivalent_nu(C,assumption) result(nu)
endif endif
mu = lattice_equivalent_mu(C,assumption) 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 end function lattice_equivalent_nu
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Equivalent shear modulus (μ) !> @brief equivalent shear modulus (μ)
!> @details https://doi.org/10.1143/JPSJ.20.635 !> @details https://doi.org/10.1143/JPSJ.20.635
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_equivalent_mu(C,assumption) result(mu) function lattice_equivalent_mu(C,assumption) result(mu)
@ -2121,9 +2122,9 @@ function lattice_equivalent_mu(C,assumption) result(mu)
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))) & 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 / 15.0_pReal
elseif(IO_lc(assumption) == 'reuss') then elseif (IO_lc(assumption) == 'reuss') then
call math_invert(S,error,C) call math_invert(S,error,C)
if(error) error stop 'matrix inversion failed' if (error) error stop 'matrix inversion failed'
mu = 15.0_pReal & 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))) / (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 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 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]) 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) 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 do i = 1, 10
call random_number(C) call random_number(C)
@ -2198,13 +2199,13 @@ subroutine selfTest
C(1,1) = C(1,1) + C(1,2) + 0.1_pReal 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(4,4) = 0.5_pReal * (C(1,1) - C(1,2))
C = lattice_symmetrize_C66(C,'cI') 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,'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,'reuss'),1.0e-12_pReal)) error stop 'equivalent_mu/reuss'
lambda = C(1,2) 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' 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' lattice_equivalent_nu(C,'reuss'),1.0e-12_pReal)) error stop 'equivalent_nu/reuss'
end subroutine selfTest end subroutine selfTest

119
src/math.f90
View File

@ -15,15 +15,15 @@ module math
implicit none implicit none
public public
#if __INTEL_COMPILER >= 1900 #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 :: & private :: &
IO, & IO, &
config config
#endif #endif
real(pReal), parameter :: PI = acos(-1.0_pReal) !< ratio of a circle's circumference to its diameter 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 :: INDEG = 180.0_pReal/PI !< conversion from radian to degree
real(pReal), parameter :: INRAD = PI/180.0_pReal !< conversion from degree into radian 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) complex(pReal), parameter :: TWOPIIMG = cmplx(0.0_pReal,2.0_pReal*PI) !< Re(0.0), Im(2xPi)
real(pReal), dimension(3,3), parameter :: & 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, intent(in),optional :: istart,iend, sortDim
integer :: ipivot,s,e,d integer :: ipivot,s,e,d
if(present(istart)) then if (present(istart)) then
s = istart s = istart
else else
s = lbound(a,2) s = lbound(a,2)
endif endif
if(present(iend)) then if (present(iend)) then
e = iend e = iend
else else
e = ubound(a,2) e = ubound(a,2)
endif endif
if(present(sortDim)) then if (present(sortDim)) then
d = sortDim d = sortDim
else else
d = 1 d = 1
@ -467,7 +467,7 @@ pure function math_inv33(A)
logical :: error logical :: error
call math_invert33(math_inv33,DetA,error,A) 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 end function math_inv33
@ -698,7 +698,7 @@ pure function math_sym33to6(m33,weighted)
integer :: i integer :: i
if(present(weighted)) then if (present(weighted)) then
w = merge(NRMMANDEL,1.0_pReal,weighted) w = merge(NRMMANDEL,1.0_pReal,weighted)
else else
w = NRMMANDEL w = NRMMANDEL
@ -725,7 +725,7 @@ pure function math_6toSym33(v6,weighted)
integer :: i integer :: i
if(present(weighted)) then if (present(weighted)) then
w = merge(INVNRMMANDEL,1.0_pReal,weighted) w = merge(INVNRMMANDEL,1.0_pReal,weighted)
else else
w = INVNRMMANDEL w = INVNRMMANDEL
@ -802,7 +802,7 @@ pure function math_sym3333to66(m3333,weighted)
integer :: i,j integer :: i,j
if(present(weighted)) then if (present(weighted)) then
w = merge(NRMMANDEL,1.0_pReal,weighted) w = merge(NRMMANDEL,1.0_pReal,weighted)
else else
w = NRMMANDEL 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 !> @details Weighted conversion (default) rearranges according to Nye and weights shear
! components according to Mandel. Advisable for matrix operations. ! components according to Mandel. Advisable for matrix operations.
! Unweighted conversion only rearranges order according to Nye ! Unweighted conversion only rearranges order according to Nye
@ -837,7 +837,7 @@ pure function math_66toSym3333(m66,weighted)
integer :: i,j integer :: i,j
if(present(weighted)) then if (present(weighted)) then
w = merge(INVNRMMANDEL,1.0_pReal,weighted) w = merge(INVNRMMANDEL,1.0_pReal,weighted)
else else
w = INVNRMMANDEL 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) pure function math_Voigt66to3333(m66)
real(pReal), dimension(3,3,3,3) :: math_Voigt66to3333 real(pReal), dimension(3,3,3,3) :: math_Voigt66to3333
real(pReal), dimension(6,6), intent(in) :: m66 !< 6x6 matrix real(pReal), dimension(6,6), intent(in) :: m66 !< 6x6 matrix
integer :: i,j integer :: i,j
@ -873,6 +874,31 @@ pure function math_Voigt66to3333(m66)
end function math_Voigt66to3333 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 !> @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), & v(2,2) + m(2, 3) * w(1), &
(m(1,1) - w(1)) * (m(2,2) - w(1)) - v(3,2)] (m(1,1) - w(1)) * (m(2,2) - w(1)) - v(3,2)]
norm = norm2(v(1:3, 1)) norm = norm2(v(1:3, 1))
fallback1: if(norm < threshold) then fallback1: if (norm < threshold) then
call math_eigh(w,v,error,m) call math_eigh(w,v,error,m)
else fallback1 else fallback1
v(1:3,1) = v(1:3, 1) / norm 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), & v(2,2) + m(2, 3) * w(2), &
(m(1,1) - w(2)) * (m(2,2) - w(2)) - v(3,2)] (m(1,1) - w(2)) * (m(2,2) - w(2)) - v(3,2)]
norm = norm2(v(1: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) call math_eigh(w,v,error,m)
else fallback2 else fallback2
v(1:3,2) = v(1:3, 2) / norm 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)))] 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) 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)) 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))) & 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) *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 & - 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) - 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) math_eigvalsh33 = math_eigvalsh(m)
else else
rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal 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(left)) math_clip = max(left,math_clip)
if (present(right)) math_clip = min(right,math_clip) if (present(right)) math_clip = min(right,math_clip)
if (present(left) .and. present(right)) then if (present(left) .and. present(right)) then
if(left>right) error stop 'left > right' if (left>right) error stop 'left > right'
endif endif
end function math_clip 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]' error stop 'math_expand [1,2] by [1,2,3] => [1,2,2,1,1,1]'
call math_sort(sort_in_,1,3,2) call math_sort(sort_in_,1,3,2)
if(any(sort_in_ /= sort_out_)) & if (any(sort_in_ /= sort_out_)) &
error stop 'math_sort' error stop 'math_sort'
if(any(math_range(5) /= range_out_)) & if (any(math_range(5) /= range_out_)) &
error stop 'math_range' 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)' 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)' error stop 'math_exp33(math_I3,128)'
call random_number(v9) 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' error stop 'math_33to9/math_9to33'
call random_number(t99) 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' error stop 'math_3333to99/math_99to3333'
call random_number(v6) 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' error stop 'math_sym33to6/math_6toSym33'
call random_number(t66) 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' 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) 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' error stop 'math_symmetric33'
call random_number(v3_1) call random_number(v3_1)
@ -1270,34 +1299,34 @@ subroutine selfTest
call random_number(v3_3) call random_number(v3_3)
call random_number(v3_4) 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)) & math_volTetrahedron(v3_1,v3_2,v3_3,v3_4),tol=1.0e-12_pReal)) &
error stop 'math_volTetrahedron' error stop 'math_volTetrahedron'
call random_number(t33) 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' 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' 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)' error stop 'math_inv33(math_I3)'
do while(abs(math_det33(t33))<1.0e-9_pReal) do while(abs(math_det33(t33))<1.0e-9_pReal)
call random_number(t33) call random_number(t33)
enddo 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' error stop 'math_inv33'
call math_invert33(t33_2,det,e,t33) 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' 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)' error stop 'math_invert33 (determinant)'
call math_invert(t33_2,e,t33) 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' error stop 'math_invert t33'
do while(math_det33(t33)<1.0e-2_pReal) ! O(det(F)) = 1 do while(math_det33(t33)<1.0e-2_pReal) ! O(det(F)) = 1
@ -1305,7 +1334,7 @@ subroutine selfTest
enddo enddo
t33_2 = math_rotationalPart(transpose(t33)) t33_2 = math_rotationalPart(transpose(t33))
t33 = math_rotationalPart(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' error stop 'math_rotationalPart'
call random_number(r) call random_number(r)
@ -1313,33 +1342,33 @@ subroutine selfTest
txx = math_eye(d) txx = math_eye(d)
allocate(txx_2(d,d)) allocate(txx_2(d,d))
call math_invert(txx_2,e,txx) 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' 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 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)' 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' error stop 'math_clip'
if(math_factorial(10) /= 3628800) & if (math_factorial(10) /= 3628800) &
error stop 'math_factorial' error stop 'math_factorial'
if(math_binomial(49,6) /= 13983816) & if (math_binomial(49,6) /= 13983816) &
error stop 'math_binomial' error stop 'math_binomial'
if(math_multinomial([1,2,3,4]) /= 12600) & if (math_multinomial([1,2,3,4]) /= 12600) &
error stop 'math_multinomial' error stop 'math_multinomial'
ijk = cshift([1,2,3],int(r*1.0e2_pReal)) 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)' error stop 'math_LeviCivita(even)'
ijk = cshift([3,2,1],int(r*2.0e2_pReal)) 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)' error stop 'math_LeviCivita(odd)'
ijk = cshift([2,2,1],int(r*2.0e2_pReal)) 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' error stop 'math_LeviCivita'
end subroutine selfTest end subroutine selfTest

14
src/phase.f90
View File

@ -230,15 +230,15 @@ module phase
end function phase_mechanical_constitutive end function phase_mechanical_constitutive
!ToDo: Merge all the stiffness functions !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 integer, intent(in) :: ph, en
real(pReal), dimension(6,6) :: C real(pReal), dimension(6,6) :: C
end function phase_homogenizedC end function phase_homogenizedC66
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(3,3,3,3), intent(in) :: C_homogenized real(pReal), dimension(6,6), intent(in) :: C66
integer, intent(in) :: ph,en integer, intent(in) :: ph,en
real(pReal), dimension(3,3,3,3) :: C real(pReal), dimension(6,6) :: C66_degraded
end function phase_damage_C end function phase_damage_C66
module function phase_f_phi(phi,co,ce) result(f) module function phase_f_phi(phi,co,ce) result(f)
integer, intent(in) :: ce,co integer, intent(in) :: ce,co
@ -299,7 +299,7 @@ module phase
public :: & public :: &
phase_init, & phase_init, &
phase_homogenizedC, & phase_homogenizedC66, &
phase_f_phi, & phase_f_phi, &
phase_f_T, & phase_f_T, &
phase_K_phi, & phase_K_phi, &

16
src/phase_damage.f90
View File

@ -139,6 +139,7 @@ module function phase_damage_constitutive(Delta_t,co,ip,el) result(converged_)
integer :: & integer :: &
ph, en ph, en
ph = material_phaseID(co,(el-1)*discretization_nIPs + ip) ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
en = material_phaseEntry(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 !> @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(3,3,3,3), intent(in) :: C_homogenized real(pReal), dimension(6,6), intent(in) :: C66
integer, intent(in) :: ph,en integer, intent(in) :: ph,en
real(pReal), dimension(3,3,3,3) :: C real(pReal), dimension(6,6) :: C66_degraded
damageType: select case (phase_damage(ph)) damageType: select case (phase_damage(ph))
case (DAMAGE_ISOBRITTLE_ID) damageType case (DAMAGE_ISOBRITTLE_ID) damageType
C = C_homogenized * damage_phi(ph,en)**2 C66_degraded = C66 * damage_phi(ph,en)**2
case default damageType case default damageType
C = C_homogenized C66_degraded = C66
end select damageType 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)) sourceType: select case (phase_damage(ph))
case (DAMAGE_ISOBRITTLE_ID) sourceType 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))) broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,en)))
if (.not. broken) then if (.not. broken) then
myOffset = damageState(ph)%offsetDeltaState myOffset = damageState(ph)%offsetDeltaState

6
src/phase_damage_isobrittle.f90
View File

@ -96,6 +96,7 @@ end function isobrittle_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state !> @brief calculates derived quantities from state
! ToDo: Use Voigt directly
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine isobrittle_deltaState(C, Fe, ph,en) module subroutine isobrittle_deltaState(C, Fe, ph,en)
@ -109,13 +110,16 @@ module subroutine isobrittle_deltaState(C, Fe, ph,en)
epsilon epsilon
real(pReal) :: & real(pReal) :: &
r_W 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) epsilon = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3)
associate(prm => param(ph), stt => state(ph), dlt => deltaState(ph)) 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)) dlt%r_W(en) = merge(r_W - stt%r_W(en), 0.0_pReal, r_W > stt%r_W(en))
end associate end associate

76
src/phase_mechanical_elastic.f90
View File

@ -15,7 +15,7 @@ submodule(phase:mechanical) elastic
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Initialize elasticity. !> @brief initialize elasticity
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine elastic_init(phases) module subroutine elastic_init(phases)
@ -62,52 +62,30 @@ end subroutine elastic_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Return 6x6 elasticity tensor. !> @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.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function elastic_C66(ph,en) result(C66) module function elastic_C66(ph,en) result(C66)
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
en en
real(pReal), dimension(6,6) :: & real(pReal), dimension(6,6) :: C66
C66
associate(prm => param(ph)) 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) if (any(phase_lattice(ph) == ['hP','tI'])) then
C66 = math_sym3333to66(math_Voigt66to3333(C66)) ! Literature data is in Voigt notation 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 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) module function elastic_mu(ph,en) result(mu)
@ -126,12 +104,13 @@ module function elastic_mu(ph,en) result(mu)
mu mu
mu = lattice_equivalent_mu(get_C66(ph,en),'voigt') mu = lattice_equivalent_mu(elastic_C66(ph,en),'voigt')
end function elastic_mu end function elastic_mu
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Return Poisson ratio. !> @brief return Poisson ratio
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function elastic_nu(ph,en) result(nu) module function elastic_nu(ph,en) result(nu)
@ -142,15 +121,16 @@ module function elastic_nu(ph,en) result(nu)
nu nu
nu = lattice_equivalent_nu(get_C66(ph,en),'voigt') nu = lattice_equivalent_nu(elastic_C66(ph,en),'voigt')
end function elastic_nu 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 !> 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, & module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
Fe, Fi, ph, en) Fe, Fi, ph, en)
@ -173,13 +153,12 @@ module subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
i, j i, j
C = math_66toSym3333(phase_homogenizedC(ph,en)) C = math_Voigt66to3333(phase_damage_C66(phase_homogenizedC66(ph,en),ph,en))
C = phase_damage_C(C,ph,en)
E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration 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 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_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 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 end do; end do
@ -188,10 +167,9 @@ end subroutine phase_hooke_SandItsTangents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns the homogenized elasticity matrix !> @brief Return the homogenized elasticity matrix.
!> ToDo: homogenizedC66 would be more consistent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module function phase_homogenizedC(ph,en) result(C) module function phase_homogenizedC66(ph,en) result(C)
real(pReal), dimension(6,6) :: C real(pReal), dimension(6,6) :: C
integer, intent(in) :: ph, en integer, intent(in) :: ph, en
@ -204,7 +182,7 @@ module function phase_homogenizedC(ph,en) result(C)
C = elastic_C66(ph,en) C = elastic_C66(ph,en)
end select plasticType end select plasticType
end function phase_homogenizedC end function phase_homogenizedC66
end submodule elastic end submodule elastic

34
src/phase_mechanical_plastic_dislotwin.f90
View File

@ -150,7 +150,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
myPlasticity = plastic_active('dislotwin') myPlasticity = plastic_active('dislotwin')
if(count(myPlasticity) == 0) return if (count(myPlasticity) == 0) return
print'(/,1x,a)', '<<<+- phase:mechanical:plastic:dislotwin init -+>>>' print'(/,1x,a)', '<<<+- phase:mechanical:plastic:dislotwin init -+>>>'
print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) 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 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)) 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 ! 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') prm%D = pl%get_asFloat('D')
if (prm%sum_N_tw + prm%sum_N_tr > 0) & 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,:) stt%gamma_sl=>plasticState(ph)%state(startIndex:endIndex,:)
dot%gamma_sl=>plasticState(ph)%dotState(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) 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 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw endIndex = endIndex + prm%sum_N_tw
@ -474,26 +474,26 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC)
ph, en ph, en
real(pReal), dimension(6,6) :: & real(pReal), dimension(6,6) :: &
homogenizedC, & homogenizedC, &
C66 C
real(pReal), dimension(:,:,:), allocatable :: & real(pReal), dimension(:,:,:), allocatable :: &
C66_tw, & C66_tw, &
C66_tr C66_tr
integer :: i integer :: i
real(pReal) :: f_unrotated 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 & f_unrotated = 1.0_pReal &
- sum(stt%f_tw(1:prm%sum_N_tw,en)) & - sum(stt%f_tw(1:prm%sum_N_tw,en)) &
- sum(stt%f_tr(1:prm%sum_N_tr,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 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 do i=1,prm%sum_N_tw
homogenizedC = homogenizedC & homogenizedC = homogenizedC &
+ stt%f_tw(i,en)*C66_tw(1:6,1:6,i) + stt%f_tw(i,en)*C66_tw(1:6,1:6,i)
@ -501,7 +501,7 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC)
end if twinActive end if twinActive
transActive: if (prm%sum_N_tr > 0) then 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 do i=1,prm%sum_N_tr
homogenizedC = homogenizedC & homogenizedC = homogenizedC &
+ stt%f_tr(i,en)*C66_tr(1:6,1:6,i) + 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 Lp = Lp * f_unrotated
dLp_dMp = dLp_dMp * 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) E_kB_T = prm%E_sb/(kB*T)
call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design? 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) 'threshold stress for twinning','Pa',prm%systems_tw)
case('f_tr') 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³') 'martensite volume fraction','m³/m³')
end select end select
@ -929,8 +929,8 @@ pure subroutine kinetics_sl(Mp,T,ph,en, &
end associate end associate
if(present(ddot_gamma_dtau_sl)) ddot_gamma_dtau_sl = ddot_gamma_dtau if (present(ddot_gamma_dtau_sl)) ddot_gamma_dtau_sl = ddot_gamma_dtau
if(present(tau_sl)) tau_sl = tau if (present(tau_sl)) tau_sl = tau
end subroutine kinetics_sl end subroutine kinetics_sl
@ -1000,7 +1000,7 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
end associate 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 end subroutine kinetics_tw
@ -1069,7 +1069,7 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
end associate 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 end subroutine kinetics_tr

23
src/rotations.f90
View File

@ -75,7 +75,6 @@ module rotations
procedure, public :: rotVector procedure, public :: rotVector
procedure, public :: rotTensor2 procedure, public :: rotTensor2
procedure, public :: rotTensor4 procedure, public :: rotTensor4
procedure, public :: rotTensor4sym
procedure, public :: misorientation procedure, public :: misorientation
procedure, public :: standardize procedure, public :: standardize
end type rotation end type rotation
@ -371,27 +370,6 @@ pure function rotTensor4(self,T,active) result(tRot)
end function rotTensor4 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 !> @brief misorientation
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
@ -400,6 +378,7 @@ pure elemental function misorientation(self,other)
type(rotation) :: misorientation type(rotation) :: misorientation
class(rotation), intent(in) :: self, other class(rotation), intent(in) :: self, other
misorientation%q = multiply_quaternion(other%q, conjugate_quaternion(self%q)) misorientation%q = multiply_quaternion(other%q, conjugate_quaternion(self%q))
end function misorientation end function misorientation