From 0b83fa0fb251d9ffb8c1dcab71b924ad4e093a98 Mon Sep 17 00:00:00 2001 From: Claudio Zambaldi Date: Wed, 28 Apr 2010 17:19:58 +0000 Subject: [PATCH] corrected (?) disorientation calc and introduced some new assisting functions --- code/crystallite.f90 | 70 ++++++++++---------- code/math.f90 | 150 +++++++++++++++++++++++++++++-------------- 2 files changed, 138 insertions(+), 82 deletions(-) diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 718642e90..1b1686ec7 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -60,7 +60,7 @@ real(pReal), dimension (:,:,:,:,:), allocatable :: & crystallite_partionedLp0,& ! plastic velocity grad at start of homog inc crystallite_subLp0,& ! plastic velocity grad at start of crystallite inc crystallite_P, & ! 1st Piola-Kirchhoff stress per grain - crystallite_misorientation ! misorientation between two neighboring ips (only calculated for single grain IPs) + crystallite_disorientation ! disorientation between two neighboring ips (only calculated for single grain IPs) real(pReal), dimension (:,:,:,:,:,:,:), allocatable :: & crystallite_dPdF, & ! individual dPdF per grain crystallite_fallbackdPdF ! dPdF fallback for non-converged grains (elastic prediction) @@ -340,7 +340,7 @@ subroutine crystallite_init(Temperature) write(6,'(a35,x,7(i5,x))') 'crystallite_orientation: ', shape(crystallite_orientation) write(6,'(a35,x,7(i5,x))') 'crystallite_orientation0: ', shape(crystallite_orientation0) write(6,'(a35,x,7(i5,x))') 'crystallite_rotation: ', shape(crystallite_rotation) - write(6,'(a35,x,7(i5,x))') 'crystallite_misorientation: ', shape(crystallite_misorientation) + write(6,'(a35,x,7(i5,x))') 'crystallite_disorientation: ', shape(crystallite_disorientation) write(6,'(a35,x,7(i5,x))') 'crystallite_dt: ', shape(crystallite_dt) write(6,'(a35,x,7(i5,x))') 'crystallite_subdt: ', shape(crystallite_subdt) write(6,'(a35,x,7(i5,x))') 'crystallite_subFrac: ', shape(crystallite_subFrac) @@ -621,7 +621,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) if (crystallite_todo(g,i,e)) then ! all undone crystallites call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & - crystallite_misorientation(:,:,g,i,e), crystallite_subdt(g,i,e), g, i, e) + crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), g, i, e) endif enddo; enddo; enddo !$OMPEND PARALLEL DO @@ -712,7 +712,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) if (crystallite_todo(g,i,e)) then ! all undone crystallites call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & - crystallite_misorientation(:,:,g,i,e), crystallite_subdt(g,i,e), g, i, e) + crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), g, i, e) delta_dotState1 = constitutive_dotState(g,i,e)%p - constitutive_previousDotState(g,i,e)%p delta_dotState2 = constitutive_previousDotState(g,i,e)%p - constitutive_previousDotState2(g,i,e)%p dot_prod12 = dot_product(delta_dotState1, delta_dotState2) @@ -869,7 +869,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) constitutive_dotState(g,i,e)%p = 0.0_pReal call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & - crystallite_misorientation(:,:,g,i,e), crystallite_subdt(g,i,e), & + crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), & g,i,e) stateConverged = crystallite_updateState(g,i,e) ! update state temperatureConverged = crystallite_updateTemperature(g,i,e) ! update temperature @@ -998,7 +998,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) if (crystallite_todo(g,i,e)) then call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), & - crystallite_misorientation(:,:,g,i,e), crystallite_subdt(g,i,e), & + crystallite_disorientation(:,:,g,i,e), crystallite_subdt(g,i,e), & g,i,e) ! collect dot state delta_dotState1 = constitutive_dotState(g,i,e)%p - constitutive_previousDotState(g,i,e)%p delta_dotState2 = constitutive_previousDotState(g,i,e)%p - constitutive_previousDotState2(g,i,e)%p @@ -1348,11 +1348,11 @@ endsubroutine ! inversion of Fp_current... invFp_current = math_inv3x3(Fp_current) if (all(invFp_current == 0.0_pReal)) then ! ... failed? - if (verboseDebugger) then + if (verboseDebugger .and. selectiveDebugger) then !$OMP CRITICAL (write2out) write(6,*) '::: integrateStress failed on invFp_current inversion',g,i,e write(6,*) - write(6,'(a11,3(i3,x),/,3(3(f12.7,x)/))') 'invFp_new at ',g,i,e,invFp_new + write(6,'(a11,i3,x,i2,x,i5,/,3(3(f12.7,x)/))') 'invFp_new at ',g,i,e,invFp_new !$OMPEND CRITICAL (write2out) endif return @@ -1406,10 +1406,10 @@ LpLoop: do if (tock < tick) debug_cumLpTicks = debug_cumLpTicks + maxticks if (verboseDebugger .and. selectiveDebugger) then !$OMP CRITICAL (write2out) - write(6,*) '::: integrateStress at ' ,g,i,e, ' ; iteration ', NiterationStress + write(6,'(a,i3,x,i2,x,i5,x,a,x,i3)') '::: integrateStress at ' ,g,i,e, ' ; iteration ', NiterationStress write(6,*) - write(6,'(a,/,3(3(f20.7,x)/))') 'Lp_constitutive', Lp_constitutive - write(6,'(a,/,3(3(f20.7,x)/))') 'Lpguess', Lpguess + write(6,'(a,/,3(3(e20.7,x)/))') 'Lp_constitutive', Lp_constitutive + write(6,'(a,/,3(3(e20.7,x)/))') 'Lpguess', Lpguess !$OMPEND CRITICAL (write2out) endif @@ -1430,7 +1430,7 @@ LpLoop: do if (any(residuum/=residuum) .and. leapfrog == 1.0) then if (verboseDebugger) then !$OMP CRITICAL (write2out) - write(6,*) '::: integrateStress encountered NaN at ',g,i,e,' ; iteration ', NiterationStress + write(6,'(a,i3,x,i2,x,i5,x,a,x,i3)') '::: integrateStress encountered NaN at ',g,i,e,' ; iteration ', NiterationStress !$OMPEND CRITICAL (write2out) endif return @@ -1462,14 +1462,14 @@ LpLoop: do invdRdLp = 0.0_pReal call math_invert(9,dRdLp,invdRdLp,dummy,error) ! invert dR/dLp --> dLp/dR if (error) then - if (verboseDebugger) then + if (verboseDebugger .and. selectiveDebugger) then !$OMP CRITICAL (write2out) - write(6,*) '::: integrateStress failed on dR/dLp inversion at ',g,i,e,' ; iteration ', NiterationStress + write(6,'(a,i3,x,i2,x,i5,x,a,x,i3)') '::: integrateStress failed on dR/dLp inversion at ',g,i,e,' ; iteration ', NiterationStress write(6,*) - write(6,'(a,/,9(9(f15.3,x)/))') 'dRdLp',dRdLp - write(6,'(a,/,9(9(f15.3,x)/))') 'dLpdT_constitutive',dLpdT_constitutive - write(6,'(a,/,3(3(f20.7,x)/))') 'Lp_constitutive',Lp_constitutive - write(6,'(a,/,3(3(f20.7,x)/))') 'Lpguess',Lpguess + write(6,'(a,/,9(9(e15.3,x)/))') 'dRdLp',dRdLp + write(6,'(a,/,9(9(e15.3,x)/))') 'dLpdT_constitutive',dLpdT_constitutive + write(6,'(a,/,3(3(e20.7,x)/))') 'Lp_constitutive',Lp_constitutive + write(6,'(a,/,3(3(e20.7,x)/))') 'Lpguess',Lpguess !$OMPEND CRITICAL (write2out) endif return @@ -1495,9 +1495,9 @@ LpLoop: do invFp_new = invFp_new/math_det3x3(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize by det call math_invert3x3(invFp_new,Fp_new,det,error) if (error) then - if (verboseDebugger) then + if (verboseDebugger .and. selectiveDebugger) then !$OMP CRITICAL (write2out) - write(6,*) '::: integrateStress failed on invFp_new inversion at ',g,i,e,' ; iteration ', NiterationStress + write(6,'(a,i3,x,i2,x,i5,x,a,x,i3)') '::: integrateStress failed on invFp_new inversion at ',g,i,e,' ; iteration ', NiterationStress write(6,*) write(6,'(a11,3(i3,x),/,3(3(f12.7,x)/))') 'invFp_new at ',g,i,e,invFp_new !$OMPEND CRITICAL (write2out) @@ -1523,7 +1523,7 @@ LpLoop: do crystallite_integrateStress = .true. if (verboseDebugger .and. selectiveDebugger) then !$OMP CRITICAL (write2out) - write(6,*) '::: integrateStress converged at ',g,i,e,' ; iteration ', NiterationStress + write(6,'(a,i3,x,i2,x,i5,x,a,x,i3)') '::: integrateStress converged at ',g,i,e,' ; iteration ', NiterationStress write(6,*) write(6,'(a,/,3(3(f12.7,x)/))') 'P / MPa',crystallite_P(:,:,g,i,e)/1e6 write(6,'(a,/,3(3(f12.7,x)/))') 'Lp',crystallite_Lp(:,:,g,i,e) @@ -1542,7 +1542,7 @@ LpLoop: do !******************************************************************** -! calculates orientations and misorientations (in case of single grain ips) +! calculates orientations and disorientations (in case of single grain ips) !******************************************************************** subroutine crystallite_orientations() @@ -1551,12 +1551,11 @@ use prec, only: pInt, & pReal use math, only: math_pDecomposition, & math_RtoQuaternion, & - math_misorientation, & + math_QuaternionDisorientation, & inDeg use FEsolving, only: FEsolving_execElem, & FEsolving_execIP use IO, only: IO_warning -!use lattice, only: lattice_symmetryTypes use material, only: material_phase, & homogenization_Ngrains, & phase_constitution, & @@ -1597,7 +1596,7 @@ logical error do g = 1,homogenization_Ngrains(mesh_element(3,e)) ! calculate orientation in terms of rotation matrix and euler angles - call math_pDecomposition(crystallite_Fe(:,:,g,i,e), U, R, error) ! polar decomposition of Fe + call math_pDecomposition(crystallite_Fe(:,:,g,i,e), U, R, error) ! polar decomposition of Fe if (error) then call IO_warning(650, e, i, g) crystallite_orientation(:,g,i,e) = (/1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal/) ! fake orientation @@ -1605,9 +1604,10 @@ logical error crystallite_orientation(:,g,i,e) = math_RtoQuaternion(R) endif - call math_misorientation( crystallite_rotation(:,g,i,e), & ! calculate grainrotation - crystallite_orientation(:,g,i,e), crystallite_orientation0(:,g,i,e), & - crystallite_symmetryID(g,i,e)) + crystallite_rotation(:,g,i,e) = & + math_QuaternionDisorientation( crystallite_orientation(:,g,i,e), & ! calculate grainrotation + crystallite_orientation0(:,g,i,e), & + crystallite_symmetryID(g,i,e) ) enddo enddo @@ -1631,17 +1631,17 @@ logical error neighboringPhase = material_phase(1,neighboring_i,neighboring_e) ! get my neighbor's crystal structure if (myPhase == neighboringPhase) then ! if my neighbor has same phase like me - call math_misorientation( crystallite_misorientation(:,n,1,i,e), & - crystallite_orientation(:,1,i,e), & - crystallite_orientation(:,1,neighboring_i,neighboring_e), & - crystallite_symmetryID(g,i,e)) ! calculate misorientation + crystallite_disorientation(:,n,1,i,e) = & + math_QuaternionDisorientation( crystallite_orientation(:,1,i,e), & + crystallite_orientation(:,1,neighboring_i,neighboring_e), & + crystallite_symmetryID(1,i,e)) ! calculate disorientation else ! for neighbor with different phase - crystallite_misorientation(:,n,1,i,e) = (/-1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal/) ! set misorientation to maximum + crystallite_disorientation(:,n,1,i,e) = (/1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal/) ! identity "rotation" endif else ! no existing neighbor - crystallite_misorientation(:,n,1,i,e) = (/1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal/) ! set misorientation to zero + crystallite_disorientation(:,n,1,i,e) = (/-1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal/) ! homomorphic identity endif enddo endif @@ -1729,7 +1729,7 @@ function crystallite_postResults(& crystallite_postResults(c+1) = constitutive_sizePostResults(g,i,e); c = c+1_pInt ! size of constitutive results crystallite_postResults(c+1:c+constitutive_sizePostResults(g,i,e)) = & constitutive_postResults(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Fe, crystallite_Fp, & - crystallite_Temperature(g,i,e), crystallite_misorientation(:,:,g,i,e), dt, & + crystallite_Temperature(g,i,e), crystallite_disorientation(:,:,g,i,e), dt, & crystallite_subdt(g,i,e), g, i, e) c = c + constitutive_sizePostResults(g,i,e) diff --git a/code/math.f90 b/code/math.f90 index fb6ddf52c..c540545c0 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -68,7 +68,8 @@ ! Symmetry operations as quaternions ! 24 for cubic, 12 for hexagonal = 36 -real(pReal), dimension(4,36), parameter :: symOperations = & +integer(pInt), dimension(2), parameter :: math_NsymOperations = (/24,12/) +real(pReal), dimension(4,36), parameter :: math_symOperations = & reshape((/& 1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal, & ! cubic symmetry operations 0.0_pReal, 0.0_pReal, 0.7071067811865476_pReal, 0.7071067811865476_pReal, & ! 2-fold symmetry @@ -146,48 +147,6 @@ real(pReal), dimension(4,36), parameter :: symOperations = & -!************************************************************************** -! calculates the misorientation for 2 orientations Q1 and Q2 (needs quaternions) -!************************************************************************** -subroutine math_misorientation(dQ, Q1, Q2, symmetryType) - - use prec, only: pReal, pInt - use IO, only: IO_warning - implicit none - - !*** input variables - real(pReal), dimension(4), intent(in) :: Q1, & ! 1st orientation - Q2 ! 2nd orientation - integer(pInt), intent(in) :: symmetryType ! Type of crystal symmetry; 1:cubic, 2:hexagonal - - !*** output variables - real(pReal), dimension(4), intent(out) :: dQ ! misorientation - - !*** local variables - real(pReal), dimension(4) :: this_dQ ! candidate for misorientation - integer(pInt) s - integer(pInt), dimension(2), parameter :: NsymOperations = (/24,12/) ! number of possible symmetry operations - real(pReal), dimension(:,:), allocatable :: mySymOperations ! symmetry Operations for my crystal symmetry - - dQ = 0.0_pReal - - if (symmetryType < 1_pInt .or. symmetryType > 2_pInt) then - dQ=NaN - !call IO_warning(700) - return - endif - - allocate(mySymOperations(4,NsymOperations(symmetryType))) - mySymOperations = symOperations(:,sum(NsymOperations(1:symmetryType-1))+1:sum(NsymOperations(1:symmetryType))) ! choose symmetry operations according to crystal symmetry - - dQ(1) = -1.0_pReal ! start with maximum misorientation angle - do s = 1,NsymOperations(symmetryType) ! loop ver symmetry operations - this_dQ = math_qMul( math_qConj(Q1), math_qMul(mySymOperations(:,s),Q2) ) ! misorientation candidate from Q1^-1*(sym*Q2) - if (this_dQ(1) > dQ(1)) dQ = this_dQ ! store if misorientation angle is smaller (cos is larger) than previous one - enddo - -endsubroutine - !************************************************************************** @@ -1446,6 +1405,25 @@ pure function math_transpose3x3(A) ENDFUNCTION +!******************************************************************** +! Rodrigues vector (x, y, z) from quaternion (w+ix+jy+kz) +!******************************************************************** + PURE FUNCTION math_QuaternionToRodrig(Q) + + use prec, only: pReal, pInt + implicit none + + real(pReal), dimension(4), intent(in) :: Q + real(pReal), dimension(3) :: math_QuaternionToRodrig + + if (Q(1) /= 0.0_pReal) then + math_QuaternionToRodrig = Q(2:4)/Q(1) + else + math_QuaternionToRodrig = NaN + endif + + ENDFUNCTION + !**************************************************************** ! rotation matrix from axis and angle (in radians) @@ -1489,23 +1467,101 @@ pure function math_transpose3x3(A) !************************************************************************** ! disorientation angle between two sets of Euler angles !************************************************************************** - pure function math_disorient(EulerA,EulerB) + pure function math_EulerMisorientation(EulerA,EulerB) use prec, only: pReal, pInt implicit none real(pReal), dimension(3), intent(in) :: EulerA,EulerB real(pReal), dimension(3,3) :: r - real(pReal) math_disorient, tr + real(pReal) math_EulerMisorientation, tr r = math_mul33x33(math_EulerToR(EulerB),transpose(math_EulerToR(EulerA))) tr = (r(1,1)+r(2,2)+r(3,3)-1.0_pReal)*0.4999999_pReal - math_disorient = abs(0.5_pReal*pi-asin(tr)) + math_EulerMisorientation = abs(0.5_pReal*pi-asin(tr)) return ENDFUNCTION +!************************************************************************** +! figures whether quat falls into stereographicc standard triangle +!************************************************************************** +pure function math_QuaternionInSST(Q, symmetryType) + + use prec, only: pReal, pInt + implicit none + + !*** input variables + real(pReal), dimension(4), intent(in) :: Q ! orientation + integer(pInt), intent(in) :: symmetryType ! Type of crystal symmetry; 1:cubic, 2:hexagonal + + !*** output variables + logical math_QuaternionInSST + + !*** local variables + real(pReal), dimension(3) :: Rodrig ! Rodrigues vector of Q + + Rodrig = math_QuaternionToRodrig(Q) + select case (symmetryType) + case (1) + math_QuaternionInSST = Rodrig(1) > Rodrig(2) .and. & + Rodrig(2) > Rodrig(3) .and. & + Rodrig(3) > 0.0_pReal + case (2) + math_QuaternionInSST = Rodrig(1) > dsqrt(3.0_pReal)*Rodrig(2) .and. & + Rodrig(2) > 0.0_pReal .and. & + Rodrig(3) > 0.0_pReal + case default + math_QuaternionInSST = .true. + end select + +endfunction + +!************************************************************************** +! calculates the disorientation for 2 orientations Q1 and Q2 (needs quaternions) +!************************************************************************** +pure function math_QuaternionDisorientation(Q1, Q2, symmetryType) + + use prec, only: pReal, pInt + use IO, only: IO_warning + implicit none + + !*** input variables + real(pReal), dimension(4), intent(in) :: Q1, & ! 1st orientation + Q2 ! 2nd orientation + integer(pInt), intent(in) :: symmetryType ! Type of crystal symmetry; 1:cubic, 2:hexagonal + + !*** output variables + real(pReal), dimension(4) :: math_QuaternionDisorientation ! disorientation + + !*** local variables + real(pReal), dimension(4) :: dQ,dQsymA,mis + integer(pInt) i,j,k,s + + dQ = math_qMul(math_qConj(Q1),Q2) + math_QuaternionDisorientation = dQ + + if (symmetryType > 0_pInt .and. symmetryType <= 2_pInt) then + s = sum(math_NsymOperations(1:symmetryType-1)) + do i = 1,2 + dQ = math_qConj(dQ) ! switch order of "from -- to" + do j = 1,math_NsymOperations(symmetryType) ! run through first crystals symmetries + dQsymA = math_qMul(math_symOperations(:,s+j),dQ) ! apply sym + do k = 1,math_NsymOperations(symmetryType) ! run through 2nd crystals symmetries + mis = math_qMul(dQsymA,math_symOperations(:,s+k)) ! apply sym + if (mis(1) < 0) & ! want positive angle + mis = -mis + if (mis(1)-math_QuaternionDisorientation(1) > -1e-8_pReal .and. & + math_QuaternionInSST(mis,symmetryType)) & + math_QuaternionDisorientation = mis ! found better one + enddo; enddo; enddo + endif + + +endfunction + + !******************************************************************** ! draw a random sample from Euler space @@ -1557,7 +1613,7 @@ endif disturb(1) = scatter * rnd(1) ! phi1 disturb(2) = sign(1.0_pReal,rnd(2))*acos(cosScatter+(1.0_pReal-cosScatter)*rnd(4)) ! Phi disturb(3) = scatter * rnd(2) ! phi2 - if (rnd(5) <= exp(-1.0_pReal*(math_disorient(origin,disturb)/scatter)**2)) exit + if (rnd(5) <= exp(-1.0_pReal*(math_EulerMisorientation(origin,disturb)/scatter)**2)) exit enddo math_sampleGaussOri = math_RtoEuler(math_mul33x33(math_EulerToR(disturb),math_EulerToR(center)))