Skip to content

Adding quaternion implementation #153

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
biancafreytag wants to merge 2 commits into
base: develop
Choose a base branch
from
Open

Adding quaternion implementation #153

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
0c66fc0
Adding quaternion implementation
Jun 13, 2018
38985d5
Change Coordinates_MaterialSystemCalculatedXdNu3D to double precision
Jun 21, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
169 changes: 102 additions & 67 deletions src/coordinate_routines.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4037,7 +4037,8 @@ SUBROUTINE Coordinates_MaterialSystemCalculatedXdNu3D(geometricInterpPointMetric
TYPE(VARYING_STRING), INTENT(OUT) :: error !<The error string
!Local Variables
REAL(DP) :: angles(3),det,dXdNu2(3,3),dXdNu3(3,3),dXdNuR(3,3),f(3),g(3),h(3), &
& Raf(3,3),Rbf(3,3),Rai(3,3),Rbi(3,3),Ras(3,3),Rbs(3,3)
& Raf(3,3),Rbf(3,3),Rai(3,3),Rbi(3,3),Ras(3,3),Rbs(3,3),a,b,c,d,s,quaternions(4), Rq(3,3)
!LOGICAL :: eulerAnglesFlag, quaternionsFlag

ENTERS("Coordinates_MaterialSystemCalculatedXdNu3D",err,error,*999)

Expand All @@ -4059,72 +4060,106 @@ SUBROUTINE Coordinates_MaterialSystemCalculatedXdNu3D(geometricInterpPointMetric
IF(diagnostics1) THEN
ENDIF

!FIBRE ANGLE(alpha) - angles(1)
!In order to rotate reference material CS by alpha(fibre angle) in anti-clockwise
!direction about its axis 3, following steps are performed.
!(a) first align reference material direction 3 with Z(spatial) axis by rotating the ref material CS.
!(b) then rotate the aligned material CS by alpha about Z axis in anti-clockwise direction
!(c) apply the inverse of step(a) to the CS in (b)
!It can be shown that steps (a),(b) and (c) are equivalent to post-multiplying
!rotation in (a) by rotation in (b). i.e. Ra*Rb

!The normalised reference material CS contains the transformation(rotation) between
!the spatial CS -> reference material CS. i.e. Raf
Raf=dXdNuR

!Initialise rotation matrix Rbf
CALL IdentityMatrix(Rbf,err,error,*999)
!Populate rotation matrix Rbf about axis 3 (Z)
Rbf(1,1)=COS(angles(1))
Rbf(1,2)=-1.0_DP*SIN(angles(1))
Rbf(2,1)=SIN(angles(1))
Rbf(2,2)=COS(angles(1))

CALL MatrixProduct(Raf,Rbf,dXdNu3,err,error,*999)

!IMBRICATION ANGLE (beta) - angles(2)
!In order to rotate alpha-rotated material CS by beta(imbrication angle) in anti-clockwise
!direction about its new axis 2, following steps are performed.
!(a) first align new material direction 2 with Y(spatial) axis by rotating the new material CS.
!(b) then rotate the aligned CS by beta about Y axis in anti-clockwise direction
!(c) apply the inverse of step(a) to the CS in (b)
!As mentioned above, (a),(b) and (c) are equivalent to post-multiplying
!rotation in (a) by rotation in (b). i.e. Rai*Rbi

!dXdNu3 contains the transformation(rotation) between
!the spatial CS -> alpha-rotated reference material CS. i.e. Rai
Rai=dXdNu3
!Initialise rotation matrix Rbi
CALL IdentityMatrix(Rbi,err,error,*999)
!Populate rotation matrix Rbi about axis 2 (Y). Note the sign change
Rbi(1,1)=COS(angles(2))
Rbi(1,3)=SIN(angles(2))
Rbi(3,1)=-1.0_DP*SIN(angles(2))
Rbi(3,3)=COS(angles(2))

CALL MatrixProduct(Rai,Rbi,dXdNu2,err,error,*999)

!SHEET ANGLE (gamma) - angles(3)
!In order to rotate alpha-beta-rotated material CS by gama(sheet angle) in anti-clockwise
!direction about its new axis 1, following steps are performed.
!(a) first align new material direction 1 with X(spatial) axis by rotating the new material CS.
!(b) then rotate the aligned CS by gama about X axis in anti-clockwise direction
!(c) apply the inverse of step(a) to the CS in (b)
!Again steps (a),(b) and (c) are equivalent to post-multiplying
!rotation in (a) by rotation in (b). i.e. Ras*Rbs

!dXdNu2 contains the transformation(rotation) between
!the spatial CS -> alpha-beta-rotated reference material CS. i.e. Ras
Ras=dXdNu2
!Initialise rotation matrix Rbs
CALL IdentityMatrix(Rbs,err,error,*999)
!Populate rotation matrix Rbs about axis 1 (X).
Rbs(2,2)=COS(angles(3))
Rbs(2,3)=-1.0_DP*SIN(angles(3))
Rbs(3,2)=SIN(angles(3))
Rbs(3,3)=COS(angles(3))

CALL MatrixProduct(Ras,Rbs,dXdNu,err,error,*999)

IF(SIZE(angle,1)==3) THEN


!FIBRE ANGLE(alpha) - angles(1)
!In order to rotate reference material CS by alpha(fibre angle) in anti-clockwise
!direction about its axis 3, following steps are performed.
!(a) first align reference material direction 3 with Z(spatial) axis by rotating the ref material CS.
!(b) then rotate the aligned material CS by alpha about Z axis in anti-clockwise direction
!(c) apply the inverse of step(a) to the CS in (b)
!It can be shown that steps (a),(b) and (c) are equivalent to post-multiplying
!rotation in (a) by rotation in (b). i.e. Ra*Rb

!The normalised reference material CS contains the transformation(rotation) between
!the spatial CS -> reference material CS. i.e. Raf
Raf=dXdNuR

!Initialise rotation matrix Rbf
CALL IdentityMatrix(Rbf,err,error,*999)
!Populate rotation matrix Rbf about axis 3 (Z)
Rbf(1,1)=COS(angles(1))
Rbf(1,2)=-1.0_DP*SIN(angles(1))
Rbf(2,1)=SIN(angles(1))
Rbf(2,2)=COS(angles(1))

CALL MatrixProduct(Raf,Rbf,dXdNu3,err,error,*999)

!IMBRICATION ANGLE (beta) - angles(2)
!In order to rotate alpha-rotated material CS by beta(imbrication angle) in anti-clockwise
!direction about its new axis 2, following steps are performed.
!(a) first align new material direction 2 with Y(spatial) axis by rotating the new material CS.
!(b) then rotate the aligned CS by beta about Y axis in anti-clockwise direction
!(c) apply the inverse of step(a) to the CS in (b)
!As mentioned above, (a),(b) and (c) are equivalent to post-multiplying
!rotation in (a) by rotation in (b). i.e. Rai*Rbi

!dXdNu3 contains the transformation(rotation) between
!the spatial CS -> alpha-rotated reference material CS. i.e. Rai
Rai=dXdNu3
!Initialise rotation matrix Rbi
CALL IdentityMatrix(Rbi,err,error,*999)
!Populate rotation matrix Rbi about axis 2 (Y). Note the sign change
Rbi(1,1)=COS(angles(2))
Rbi(1,3)=SIN(angles(2))
Rbi(3,1)=-1.0_DP*SIN(angles(2))
Rbi(3,3)=COS(angles(2))

CALL MatrixProduct(Rai,Rbi,dXdNu2,err,error,*999)

!SHEET ANGLE (gamma) - angles(3)
!In order to rotate alpha-beta-rotated material CS by gama(sheet angle) in anti-clockwise
!direction about its new axis 1, following steps are performed.
!(a) first align new material direction 1 with X(spatial) axis by rotating the new material CS.
!(b) then rotate the aligned CS by gama about X axis in anti-clockwise direction
!(c) apply the inverse of step(a) to the CS in (b)
!Again steps (a),(b) and (c) are equivalent to post-multiplying
!rotation in (a) by rotation in (b). i.e. Ras*Rbs

!dXdNu2 contains the transformation(rotation) between
!the spatial CS -> alpha-beta-rotated reference material CS. i.e. Ras
Ras=dXdNu2
!Initialise rotation matrix Rbs
CALL IdentityMatrix(Rbs,err,error,*999)
!Populate rotation matrix Rbs about axis 1 (X).
Rbs(2,2)=COS(angles(3))
Rbs(2,3)=-1.0_DP*SIN(angles(3))
Rbs(3,2)=SIN(angles(3))
Rbs(3,3)=COS(angles(3))

CALL MatrixProduct(Ras,Rbs,dXdNu,err,error,*999)

ELSE

quaternions=0.0_DP
quaternions(1:SIZE(angle,1))=angle(1:SIZE(angle,1))
a=quaternions(1)
b=quaternions(2)
c=quaternions(3)
d=quaternions(4)
s=1.0_DP/(a**2 + b**2 + c**2 + d**2)
CALL IdentityMatrix(Rq,err,error,*999)
!Populate rotation matrix Rbf about axis 3 (Z)

Rq(1,1)=(1.0_DP-2.0_DP*s*(c**2+d**2))
Rq(1,2)=(2.0_DP*s*(b*c-a*d))
Rq(1,3)=(2.0_DP*s*(b*d+a*c))
Rq(2,1)=(2.0_DP*s*(b*c+a*d))
Rq(2,2)=(1.0_DP-2.0_DP*s*(b**2+d**2))
Rq(2,3)=(2.0_DP*s*(c*d-a*b))
Rq(3,1)=(2.0_DP*s*(b*d-a*c))
Rq(3,2)=(2.0_DP*s*(c*d+a*b))
Rq(3,3)=(1.0_DP-2.0_DP*s*(b**2+c**2))

!WRITE(*,*) "USING QUATERNIONS"

CALL MatrixProduct(Rq,dXdNuR,dXdNu,err,error,*999)
ENDIF




IF(diagnostics1) THEN
CALL WriteString(DIAGNOSTIC_OUTPUT_TYPE,"",err,error,*999)
Expand Down