From da23c916cadd4e0bdb18d2d086635c28464c546d Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Sun, 21 Nov 2021 15:49:04 -0500 Subject: [PATCH] polish --- src/lattice.f90 | 80 +++++++++--------- src/math.f90 | 94 +++++++++++----------- src/phase_mechanical_elastic.f90 | 10 +-- src/phase_mechanical_plastic_dislotwin.f90 | 28 +++---- 4 files changed, 105 insertions(+), 107 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 4b87fccbf..0e24d1b18 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -405,7 +405,7 @@ module lattice contains !-------------------------------------------------------------------------------------------------- -!> @brief Module initialization +!> @brief module initialization !-------------------------------------------------------------------------------------------------- subroutine lattice_init @@ -417,7 +417,7 @@ end subroutine lattice_init !-------------------------------------------------------------------------------------------------- -!> @brief Characteristic shear for twinning +!> @brief characteristic shear for twinning !-------------------------------------------------------------------------------------------------- function lattice_characteristicShear_Twin(Ntwin,lattice,CoverA) result(characteristicShear) @@ -491,7 +491,7 @@ end function lattice_characteristicShear_Twin !-------------------------------------------------------------------------------------------------- -!> @brief Rotated elasticity matrices for twinning in 66-vector notation +!> @brief rotated elasticity matrices for twinning in 6x6-matrix notation !-------------------------------------------------------------------------------------------------- function lattice_C66_twin(Ntwin,C66,lattice,CoverA) @@ -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 diff --git a/src/math.f90 b/src/math.f90 index 2515c64d7..5d2609237 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -15,15 +15,15 @@ module math implicit none public #if __INTEL_COMPILER >= 1900 - ! do not make use associated entities available to other modules + ! do not make use of associated entities available to other modules private :: & IO, & config #endif real(pReal), parameter :: PI = acos(-1.0_pReal) !< ratio of a circle's circumference to its diameter - real(pReal), parameter :: INDEG = 180.0_pReal/PI !< conversion from radian into degree - real(pReal), parameter :: INRAD = PI/180.0_pReal !< conversion from degree into radian + real(pReal), parameter :: INDEG = 180.0_pReal/PI !< conversion from radian to degree + real(pReal), parameter :: INRAD = PI/180.0_pReal !< conversion from degree to radian complex(pReal), parameter :: TWOPIIMG = cmplx(0.0_pReal,2.0_pReal*PI) !< Re(0.0), Im(2xPi) real(pReal), dimension(3,3), parameter :: & @@ -132,19 +132,19 @@ pure recursive subroutine math_sort(a, istart, iend, sortDim) integer, intent(in),optional :: istart,iend, sortDim integer :: ipivot,s,e,d - if(present(istart)) then + if (present(istart)) then s = istart else s = lbound(a,2) endif - if(present(iend)) then + if (present(iend)) then e = iend else e = ubound(a,2) endif - if(present(sortDim)) then + if (present(sortDim)) then d = sortDim else d = 1 @@ -467,7 +467,7 @@ pure function math_inv33(A) logical :: error call math_invert33(math_inv33,DetA,error,A) - if(error) math_inv33 = 0.0_pReal + if (error) math_inv33 = 0.0_pReal end function math_inv33 @@ -698,7 +698,7 @@ pure function math_sym33to6(m33,weighted) integer :: i - if(present(weighted)) then + if (present(weighted)) then w = merge(NRMMANDEL,1.0_pReal,weighted) else w = NRMMANDEL @@ -725,7 +725,7 @@ pure function math_6toSym33(v6,weighted) integer :: i - if(present(weighted)) then + if (present(weighted)) then w = merge(INVNRMMANDEL,1.0_pReal,weighted) else w = INVNRMMANDEL @@ -802,7 +802,7 @@ pure function math_sym3333to66(m3333,weighted) integer :: i,j - if(present(weighted)) then + if (present(weighted)) then w = merge(NRMMANDEL,1.0_pReal,weighted) else w = NRMMANDEL @@ -822,7 +822,7 @@ end function math_sym3333to66 !-------------------------------------------------------------------------------------------------- -!> @brief convert 66 matrix into symmetric 3x3x3x3 matrix +!> @brief convert 6x6 matrix into symmetric 3x3x3x3 matrix !> @details Weighted conversion (default) rearranges according to Nye and weights shear ! components according to Mandel. Advisable for matrix operations. ! Unweighted conversion only rearranges order according to Nye @@ -837,7 +837,7 @@ pure function math_66toSym3333(m66,weighted) integer :: i,j - if(present(weighted)) then + if (present(weighted)) then w = merge(INVNRMMANDEL,1.0_pReal,weighted) else w = INVNRMMANDEL @@ -854,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 diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index 5792ea3cd..1b39f51ee 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -15,7 +15,7 @@ submodule(phase:mechanical) elastic contains !-------------------------------------------------------------------------------------------------- -!> @brief Initialize elasticity. +!> @brief initialize elasticity !-------------------------------------------------------------------------------------------------- module subroutine elastic_init(phases) @@ -62,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 !-------------------------------------------------------------------------------------------------- diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index de027ae44..10522737c 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -150,7 +150,7 @@ module function plastic_dislotwin_init() result(myPlasticity) myPlasticity = plastic_active('dislotwin') - if(count(myPlasticity) == 0) return + if (count(myPlasticity) == 0) return print'(/,1x,a)', '<<<+- phase:mechanical:plastic:dislotwin init -+>>>' print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) @@ -173,7 +173,7 @@ module function plastic_dislotwin_init() result(myPlasticity) do ph = 1, phases%length - if(.not. myPlasticity(ph)) cycle + if (.not. myPlasticity(ph)) cycle associate(prm => param(ph), dot => dotState(ph), stt => state(ph), dst => dependentState(ph)) @@ -368,7 +368,7 @@ module function plastic_dislotwin_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! parameters required for several mechanisms and their interactions - if(prm%sum_N_sl + prm%sum_N_tw + prm%sum_N_tw > 0) & + if (prm%sum_N_sl + prm%sum_N_tw + prm%sum_N_tw > 0) & prm%D = pl%get_asFloat('D') if (prm%sum_N_tw + prm%sum_N_tr > 0) & @@ -425,7 +425,7 @@ module function plastic_dislotwin_init() result(myPlasticity) stt%gamma_sl=>plasticState(ph)%state(startIndex:endIndex,:) dot%gamma_sl=>plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) - if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' + if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tw @@ -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