diff --git a/src/crystal.f90 b/src/crystal.f90 index ff3291010..455f7ed1a 100644 --- a/src/crystal.f90 +++ b/src/crystal.f90 @@ -404,6 +404,7 @@ module crystal crystal_slip_direction, & crystal_slip_transverse, & crystal_labels_slip, & + crystal_CorrespondenceMatrix_twin, & crystal_labels_twin contains @@ -2111,6 +2112,60 @@ function getlabels(active,potential,system) result(labels) end function getlabels +!-------------------------------------------------------------------------------------------------- +!> @brief correspondence matrix for twinning +!> details only active twin systems are considered +!-------------------------------------------------------------------------------------------------- +function crystal_CorrespondenceMatrix_twin(Ntwin,lattice,cOverA) result(CorrespondenceMatrix) + + integer, dimension(:), intent(in) :: Ntwin !< number of active twin systems per family + character(len=*), intent(in) :: lattice !< lattice structure + real(pReal), intent(in) :: cOverA !< c/a ratio + real(pReal), dimension(3,3,sum(Ntwin)) :: CorrespondenceMatrix + + real(pReal), dimension(3,3,sum(Ntwin)) :: coordinateSystem + real(pReal), dimension(sum(Ntwin)) :: characteristicShearTwin + real(pReal), dimension(3,3,sum(Ntwin)) :: SchmidMatrixTwin + real(pReal), dimension(:,:), allocatable :: twinSystems + integer, dimension(:), allocatable :: NtwinMax + integer :: i + + select case(lattice) + case('cF') + NtwinMax = CF_NTWINSYSTEM + twinSystems = CF_SYSTEMTWIN + case('cI') + NtwinMax = CI_NTWINSYSTEM + twinSystems = CI_SYSTEMTWIN + case('hP') + NtwinMax = HP_NTWINSYSTEM + twinSystems = HP_SYSTEMTWIN !< the twin system matrix is different from V2.0 + case default + allocate(NtwinMax(0)) + call IO_error(137,ext_msg='crystal_CorrespondenceMatrix_twin: '//trim(lattice)) + end select + + if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0)) & + call IO_error(145,ext_msg='Ntwin '//trim(lattice)) + if (any(Ntwin < 0)) & + call IO_error(144,ext_msg='Ntwin '//trim(lattice)) + + coordinateSystem = buildCoordinateSystem(Ntwin,NtwinMax,twinSystems,lattice,cOverA) + ! characteristicShearTwin = 0.0_pReal*lattice_characteristicShear_Twin(Ntwin,lattice,cOverA) ! for removing shear from CorrespondenceMatrix + characteristicShearTwin = crystal_characteristicShear_Twin(Ntwin,lattice,cOverA) + SchmidMatrixTwin = crystal_SchmidMatrix_twin(Ntwin,lattice,cOverA) + + !write(6,*)'coordinate system', coordinateSystem(1:3,2,1) + + !CorrespondenceMatrix(1:3,1:3,1) = math_axisAngleToR(coordinateSystem(1:3,2,6), 180.0_pReal*INRAD) ! delete this + + do i = 1, sum(Ntwin) + CorrespondenceMatrix(1:3,1:3,i) = matmul(math_axisAngleToR(coordinateSystem(1:3,2,i), & + 180.0_pReal*INRAD), MATH_I3 + characteristicShearTwin(i)* & + SchmidMatrixTwin(1:3,1:3,i)) + enddo + +end function crystal_CorrespondenceMatrix_twin !-------------------------------------------------------------------------------------------------- !> @brief Equivalent Poisson's ratio (ν) @@ -2307,4 +2362,7 @@ subroutine crystal_selfTest() end subroutine crystal_selfTest + + + end module crystal diff --git a/src/math.f90 b/src/math.f90 index ec97846c0..c119f323a 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1470,4 +1470,46 @@ subroutine math_selfTest() end subroutine math_selfTest +!-------------------------------------------------------------------------------------------------- +!> @brief rotation matrix from axis and angle (in radians) +!> @details rotation matrix is meant to represent a ACTIVE rotation +!> @details (see http://en.wikipedia.org/wiki/Euler_angles for definitions) +!> @details formula for active rotation taken from http://mathworld.wolfram.com/RodriguesRotationFormula.html +!> @details equivalent to eu2om (P=-1) from "D Rowenhorst et al. Consistent representations of and +!> @details conversions between 3D rotations, Model. Simul. Mater. Sci. Eng. 23-8 (2015)" +!-------------------------------------------------------------------------------------------------- +pure function math_axisAngleToR(axis,omega) + + implicit none + real(pReal), dimension(3,3) :: math_axisAngleToR + real(pReal), dimension(3), intent(in) :: axis + real(pReal), intent(in) :: omega + real(pReal), dimension(3) :: n + real(pReal) :: norm,s,c,c1 + + norm = norm2(axis) + wellDefined: if (norm > 1.0e-8_pReal) then + n = axis/norm ! normalize axis to be sure + + s = sin(omega) + c = cos(omega) + c1 = 1.0_pReal - c + + math_axisAngleToR(1,1) = c + c1*n(1)**2.0_pReal + math_axisAngleToR(1,2) = c1*n(1)*n(2) - s*n(3) + math_axisAngleToR(1,3) = c1*n(1)*n(3) + s*n(2) + + math_axisAngleToR(2,1) = c1*n(1)*n(2) + s*n(3) + math_axisAngleToR(2,2) = c + c1*n(2)**2.0_pReal + math_axisAngleToR(2,3) = c1*n(2)*n(3) - s*n(1) + + math_axisAngleToR(3,1) = c1*n(1)*n(3) - s*n(2) + math_axisAngleToR(3,2) = c1*n(2)*n(3) + s*n(1) + math_axisAngleToR(3,3) = c + c1*n(3)**2.0_pReal + else wellDefined + math_axisAngleToR = math_I3 + endif wellDefined + +end function math_axisAngleToR + end module math diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90 index c1767184b..034b41b84 100644 --- a/src/phase_mechanical_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanical_plastic_phenopowerlaw.f90 @@ -32,7 +32,8 @@ submodule(phase:plastic) phenopowerlaw P_sl, & P_tw, & P_nS_pos, & - P_nS_neg + P_nS_neg, & + CorrespondenceMatrix integer :: & sum_N_sl, & !< total number of active slip system sum_N_tw !< total number of active twin systems @@ -203,12 +204,13 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) prm%c_4 = math_expand(pl%get_as1dReal('c_4', requiredSize=size(N_tw), & defaultVal=misc_zeros(size(N_tw))), N_tw) + prm%CorrespondenceMatrix = crystal_CorrespondenceMatrix_twin(N_tw,phase_lattice(ph),phase_cOverA(ph)) prm%gamma_char = crystal_characteristicShear_twin(N_tw,phase_lattice(ph),phase_cOverA(ph)) prm%h_tw_tw = crystal_interaction_TwinByTwin(N_tw,pl%get_as1dReal('h_tw-tw'),phase_lattice(ph)) prm%P_tw = crystal_SchmidMatrix_twin(N_tw,phase_lattice(ph),phase_cOverA(ph)) prm%systems_tw = crystal_labels_twin(N_tw,phase_lattice(ph)) - + ! sanity checks if (any(prm%dot_gamma_0_tw <= 0.0_pREAL)) extmsg = trim(extmsg)//' dot_gamma_0_tw' if (any(prm%n_tw <= 0.0_pREAL)) extmsg = trim(extmsg)//' n_tw' @@ -225,6 +227,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) prm%h_0_tw_tw, & source=emptyRealArray) allocate(prm%h_tw_tw(0,0)) + allocate(prm%CorrespondenceMatrix(0,0,0)) end if twinActive !-------------------------------------------------------------------------------------------------- @@ -511,7 +514,7 @@ pure subroutine kinetics_tw(Mp,ph,en,& associate(prm => param(ph), stt => state(ph)) - tau_tw = [(math_tensordot(Mp,prm%P_tw(1:3,1:3,i)),i=1,prm%sum_N_tw)] + tau_tw = [(math_tensordot(Mp,prm%CorrespondenceMatrix(1:3,1:3,i)),i=1,prm%sum_N_tw)] where(tau_tw > 0.0_pREAL) dot_gamma_tw = (1.0_pREAL-sum(stt%gamma_tw(:,en)/prm%gamma_char)) & ! only twin in untwinned volume fraction diff --git a/src/prec.mod b/src/prec.mod new file mode 100644 index 000000000..5fe4cd9e8 Binary files /dev/null and b/src/prec.mod differ