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

fixed some errors in QuaternionToEuler, RToQuaternion to get the correct texture representation in the output.

This commit is contained in:
Denny Tjahjanto 2010-05-26 15:52:54 +00:00
parent e5f0af638e
commit d114a600c3
4 changed files with 70 additions and 49 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
code/IO.f90
View File

@ -262,6 +262,7 @@ endsubroutine
IO_hybridIA(1,i) = deltas(1)*(mod(bin/(steps(3)*steps(2)),steps(1))+center) ! phi1
IO_hybridIA(2,i) = deltas(2)*(mod(bin/ steps(3) ,steps(2))+center) ! Phi
IO_hybridIA(3,i) = deltas(3)*(mod(bin ,steps(3))+center) ! phi2
! write(6,'(i6,x,a3,i5,x,a5,x,3(f8.5,x))') i,'bin',bin,'euler',IO_hybridIA(:,i)
binSet(j) = binSet(i)
enddo
close(999)

9
code/crystallite.f90
View File

@ -86,7 +86,8 @@ subroutine crystallite_init(Temperature)
use debug, only: debug_info, &
debug_reset
use math, only: math_I3, &
math_EulerToR
math_EulerToR, &
math_RToEuler
use FEsolving, only: FEsolving_execElem, &
FEsolving_execIP
use mesh, only: mesh_element, &
@ -1043,7 +1044,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
k = kl(1,mod(cycleCounter/iJacoStiffness,9)+1,g,i,e)
l = kl(2,mod(cycleCounter/iJacoStiffness,9)+1,g,i,e)
crystallite_dPdF(:,:,k,l,g,i,e) = (crystallite_P(:,:,g,i,e) - storedP(:,:,g,i,e))/pert_Fg ! ... use tangent dP_ij/dFg_kl
elseif (.not. storedConvergenceFlag(g,i,e)) then ! if crystallite didnÕt converge before...
elseif (.not. storedConvergenceFlag(g,i,e)) then ! if crystallite didn<EFBFBD>t converge before...
crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e) ! ... use (elastic) fallback
endif
enddo; enddo; enddo
@ -1549,6 +1550,8 @@ use prec, only: pInt, &
pReal
use math, only: math_pDecomposition, &
math_RtoQuaternion, &
math_RtoEuler, &
math_QuaternionToEuler, &
math_QuaternionDisorientation, &
inDeg, &
math_qConj
@ -1557,6 +1560,7 @@ use FEsolving, only: FEsolving_execElem, &
use IO, only: IO_warning
use material, only: material_phase, &
homogenization_Ngrains, &
material_EulerAngles, &
phase_constitution
use mesh, only: mesh_element, &
mesh_ipNeighborhood, &
@ -1588,6 +1592,7 @@ logical error
!$OMP PARALLEL DO
n = 0_pInt
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1,homogenization_Ngrains(mesh_element(3,e))

5
code/material.f90
View File

@ -667,6 +667,9 @@ subroutine material_populateGrains()
constituentGrain = constituentGrain + myNorientations
endif
! do j = 1,myNorientations
! write(6,'(i6,x,a5,x,3(f8.5,x))') j,'euler',orientationOfGrain(:,grain+j)
! enddo
! ----------------------------------------------------------------------------
symExtension = texture_symmetry(textureID) - 1_pInt
if (symExtension > 0_pInt) then ! sample symmetry
@ -687,6 +690,7 @@ subroutine material_populateGrains()
enddo ! constituent
! ----------------------------------------------------------------------------
if (.not. microstructure_elemhomo(micro)) then ! unless element homogeneous, reshuffle grains
do i=1,myNgrains-1 ! walk thru grains
call random_number(rnd)
@ -720,6 +724,7 @@ subroutine material_populateGrains()
material_phase(g,i,e) = phaseOfGrain(grain+(i-1)*dGrains+g)
material_EulerAngles(:,g,i,e) = orientationOfGrain(:,grain+(i-1)*dGrains+g)
end forall
grain = grain + FE_Nips(mesh_element(2,e)) * dGrains ! wind forward by Nips*NgrainsPerIP
endif
! write (6,*) e

102
code/math.f90
View File

@ -431,7 +431,7 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = &
!**************************************************************************
! matrix multiplication 33x33 = 3x3
!**************************************************************************´
!**************************************************************************<EFBFBD>
pure function math_mul33x33(A,B)
use prec, only: pReal, pInt
@ -1256,9 +1256,9 @@ pure function math_transpose3x3(A)
if(val > 1.0_pReal) val = 1.0_pReal
if(val < -1.0_pReal) val = -1.0_pReal
math_RtoEuler(2) = acos(val)
math_RtoEuler(2) = dacos(val)
if(math_RtoEuler(2) < 1.0e-30_pReal) then
if(math_RtoEuler(2) < 1.0e-8_pReal) then
! calculate phi2
math_RtoEuler(3) = 0.0_pReal
! calculate phi1
@ -1266,7 +1266,7 @@ pure function math_transpose3x3(A)
if(val > 1.0_pReal) val = 1.0_pReal
if(val < -1.0_pReal) val = -1.0_pReal
math_RtoEuler(1) = acos(val)
math_RtoEuler(1) = dacos(val)
if(R(2,1) > 0.0_pReal) math_RtoEuler(1) = 2.0_pReal*pi-math_RtoEuler(1)
else
! calculate phi2
@ -1274,14 +1274,14 @@ pure function math_transpose3x3(A)
if(val > 1.0_pReal) val = 1.0_pReal
if(val < -1.0_pReal) val = -1.0_pReal
math_RtoEuler(3) = acos(val)
math_RtoEuler(3) = dacos(val)
if(R(1,3) < 0.0) math_RtoEuler(3) = 2.0_pReal*pi-math_RtoEuler(3)
! calculate phi1
val=-R(3,2)/sin(math_RtoEuler(2))
if(val > 1.0_pReal) val = 1.0_pReal
if(val < -1.0_pReal) val = -1.0_pReal
math_RtoEuler(1) = acos(val)
math_RtoEuler(1) = dacos(val)
if(R(3,1) < 0.0) math_RtoEuler(1) = 2.0_pReal*pi-math_RtoEuler(1)
end if
return
@ -1297,43 +1297,48 @@ pure function math_transpose3x3(A)
use prec, only: pReal, pInt
implicit none
real(pReal), dimension (3,3), parameter :: evenPerm = reshape((/1,2,3,2,3,1,3,1,2/),(/3,3/))
real(pReal), dimension (3,3), intent(in) :: R
real(pReal), dimension(4) :: math_RtoQuaternion
real(pReal) largestDiagElem,root
integer(pInt) i,largest
real(pReal), dimension(4) :: absQ,math_RtoQuaternion
integer(pInt), dimension(1) :: largest
! math adopted from http://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation#Quaternions_and_other_representations_of_rotations
! assuming that above corresponds to active rotation, we negate rotation angle
! math adopted from http://code.google.com/p/mtex/source/browse/trunk/geometry/geometry_tools/mat2quat.m
largestDiagElem = 0.0_pReal
largest = 0_pInt
math_RtoQuaternion = 0.0_pReal
do i = 1,3
if (dabs(R(i,i)) >= largestDiagElem) then
largestDiagElem = dabs(R(i,i))
largest = i
endif
enddo
absQ(1) = 0.5_pReal * dsqrt(1_pReal+R(1,1)+R(2,2)+R(3,3))
absQ(2) = 0.5_pReal * dsqrt(1_pReal+R(1,1)-R(2,2)-R(3,3))
absQ(3) = 0.5_pReal * dsqrt(1_pReal-R(1,1)+R(2,2)-R(3,3))
absQ(4) = 0.5_pReal * dsqrt(1_pReal-R(1,1)-R(2,2)+R(3,3))
root = dsqrt( max(0.0_pReal, &
1.0_pReal - (R(1,1)+R(2,2)+R(3,3)) + 2.0_pReal*R(largest,largest)) )
if (root < 1.0e-8_pReal) then
math_RtoQuaternion = 0.0_pReal
math_RtoQuaternion(1) = 1.0_pReal
else
math_RtoQuaternion(1) = 0.5_pReal / root * &
( R(evenPerm(2,largest),evenPerm(3,largest)) - &
R(evenPerm(3,largest),evenPerm(2,largest)) )
math_RtoQuaternion(1+evenPerm(1,largest)) = 0.5_pReal * root
math_RtoQuaternion(1+evenPerm(2,largest)) = 0.5_pReal / root * &
( R(evenPerm(1,largest),evenPerm(2,largest)) + &
R(evenPerm(2,largest),evenPerm(1,largest)) )
math_RtoQuaternion(1+evenPerm(3,largest)) = 0.5_pReal / root * &
( R(evenPerm(3,largest),evenPerm(1,largest)) + &
R(evenPerm(1,largest),evenPerm(3,largest)) )
endif
largest = maxloc(absQ)
select case(largest(1))
case (1)
math_RtoQuaternion(2) = R(2,3)-R(3,2)
math_RtoQuaternion(3) = R(3,1)-R(1,3)
math_RtoQuaternion(4) = R(1,2)-R(2,1)
case (2)
math_RtoQuaternion(1) = R(2,3)-R(3,2)
math_RtoQuaternion(3) = R(1,2)+R(2,1)
math_RtoQuaternion(4) = R(3,1)+R(1,3)
case (3)
math_RtoQuaternion(1) = R(3,1)-R(1,3)
math_RtoQuaternion(2) = R(1,2)+R(2,1)
math_RtoQuaternion(4) = R(2,3)+R(3,2)
case (4)
math_RtoQuaternion (1) = R(1,2)-R(2,1)
math_RtoQuaternion (2) = R(3,1)+R(1,3)
math_RtoQuaternion (3) = R(3,2)+R(2,3)
end select
math_RtoQuaternion = math_RtoQuaternion*0.25_pReal/maxval(absQ)
math_RtoQuaternion(largest(1)) = maxval(absQ)
return
@ -1515,19 +1520,24 @@ pure function math_transpose3x3(A)
real(pReal), dimension(4), intent(in) :: Q
real(pReal), dimension(3) :: math_QuaternionToEuler
math_QuaternionToEuler(2) = dacos(1.0_pReal-2.0_pReal*(Q(2)*Q(2)+Q(3)*Q(3)))
math_QuaternionToEuler(1) = atan2(Q(1)*Q(3)+Q(2)*Q(4), Q(1)*Q(2)-Q(3)*Q(4))
if (math_QuaternionToEuler(1) < 0.0_pReal) &
math_QuaternionToEuler(1) = math_QuaternionToEuler(1) + 2.0_pReal * pi
if (dabs(math_QuaternionToEuler(2)) < 1.0e-3_pReal) then
math_QuaternionToEuler(1) = 2.0_pReal*dacos(Q(1))
math_QuaternionToEuler(3) = 0.0_pReal
else
math_QuaternionToEuler(1) = datan2(Q(1)*Q(3)+Q(2)*Q(4), Q(1)*Q(2)-Q(3)*Q(4))
if (math_QuaternionToEuler(1) < 0.0_pReal) &
math_QuaternionToEuler(1) = math_QuaternionToEuler(1) + 2.0_pReal * pi
math_QuaternionToEuler(3) = datan2(-Q(1)*Q(3)+Q(2)*Q(4), Q(1)*Q(2)+Q(3)*Q(4))
if (math_QuaternionToEuler(3) < 0.0_pReal) &
math_QuaternionToEuler(3) = math_QuaternionToEuler(3) + 2.0_pReal * pi
endif
math_QuaternionToEuler(2) = acos(1.0_pReal-2.0_pReal*(Q(2)*Q(2)+Q(3)*Q(3)))
if (math_QuaternionToEuler(2) < 0.0_pReal) &
math_QuaternionToEuler(2) = math_QuaternionToEuler(2) + pi
math_QuaternionToEuler(3) = atan2(-Q(1)*Q(3)+Q(2)*Q(4), Q(1)*Q(2)+Q(3)*Q(4))
if (math_QuaternionToEuler(3) < 0.0_pReal) &
math_QuaternionToEuler(3) = math_QuaternionToEuler(3) + 2.0_pReal * pi
return
endfunction