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

corrected (?) disorientation calc and introduced some new assisting functions

This commit is contained in:
Claudio Zambaldi 2010-04-28 17:19:58 +00:00
parent 35cebfb132
commit 0b83fa0fb2
2 changed files with 138 additions and 82 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

70
code/crystallite.f90
View File

@ -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)

150
code/math.f90
View File

@ -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)))