Merge branch 'direct-R-from-F-calculation' into 'development'
Direct r from f calculation See merge request damask/DAMASK!194
This commit is contained in:
Franz Roters
commit
b7e03364c4
|
|
@ -104,7 +104,7 @@ def set_of_quaternions():
|
|||
[1.0,-1.0,-1.0,-1.0],
|
||||
])
|
||||
specials /= np.linalg.norm(specials,axis=1).reshape(-1,1)
|
||||
specials_scatter = specials + np.broadcast_to(np.random.rand(4)*scatter,specials.shape)
|
||||
specials_scatter = specials + np.broadcast_to((np.random.rand(4)*2.-1.)*scatter,specials.shape)
|
||||
specials_scatter /= np.linalg.norm(specials_scatter,axis=1).reshape(-1,1)
|
||||
specials_scatter[specials_scatter[:,0]<0]*=-1
|
||||
|
||||
|
|
|
|||
|
|
@ -130,21 +130,21 @@ def qu2ho(qu):
|
|||
|
||||
|
||||
#---------- Rotation matrix ----------
|
||||
def om2qu(a):
|
||||
trace = a[0,0] + a[1,1] + a[2,2]
|
||||
def om2qu(om):
|
||||
trace = om.trace()
|
||||
if trace > 0:
|
||||
s = 0.5 / np.sqrt(trace+ 1.0)
|
||||
qu = np.array([0.25 / s,( a[2,1] - a[1,2] ) * s,( a[0,2] - a[2,0] ) * s,( a[1,0] - a[0,1] ) * s])
|
||||
qu = np.array([0.25 / s,( om[2,1] - om[1,2] ) * s,( om[0,2] - om[2,0] ) * s,( om[1,0] - om[0,1] ) * s])
|
||||
else:
|
||||
if ( a[0,0] > a[1,1] and a[0,0] > a[2,2] ):
|
||||
s = 2.0 * np.sqrt( 1.0 + a[0,0] - a[1,1] - a[2,2])
|
||||
qu = np.array([ (a[2,1] - a[1,2]) / s,0.25 * s,(a[0,1] + a[1,0]) / s,(a[0,2] + a[2,0]) / s])
|
||||
elif (a[1,1] > a[2,2]):
|
||||
s = 2.0 * np.sqrt( 1.0 + a[1,1] - a[0,0] - a[2,2])
|
||||
qu = np.array([ (a[0,2] - a[2,0]) / s,(a[0,1] + a[1,0]) / s,0.25 * s,(a[1,2] + a[2,1]) / s])
|
||||
if ( om[0,0] > om[1,1] and om[0,0] > om[2,2] ):
|
||||
s = 2.0 * np.sqrt( 1.0 + om[0,0] - om[1,1] - om[2,2])
|
||||
qu = np.array([ (om[2,1] - om[1,2]) / s,0.25 * s,(om[0,1] + om[1,0]) / s,(om[0,2] + om[2,0]) / s])
|
||||
elif (om[1,1] > om[2,2]):
|
||||
s = 2.0 * np.sqrt( 1.0 + om[1,1] - om[0,0] - om[2,2])
|
||||
qu = np.array([ (om[0,2] - om[2,0]) / s,(om[0,1] + om[1,0]) / s,0.25 * s,(om[1,2] + om[2,1]) / s])
|
||||
else:
|
||||
s = 2.0 * np.sqrt( 1.0 + a[2,2] - a[0,0] - a[1,1] )
|
||||
qu = np.array([ (a[1,0] - a[0,1]) / s,(a[0,2] + a[2,0]) / s,(a[1,2] + a[2,1]) / s,0.25 * s])
|
||||
s = 2.0 * np.sqrt( 1.0 + om[2,2] - om[0,0] - om[1,1] )
|
||||
qu = np.array([ (om[1,0] - om[0,1]) / s,(om[0,2] + om[2,0]) / s,(om[1,2] + om[2,1]) / s,0.25 * s])
|
||||
if qu[0]<0: qu*=-1
|
||||
return qu*np.array([1.,_P,_P,_P])
|
||||
|
||||
|
|
@ -163,7 +163,6 @@ def om2eu(om):
|
|||
|
||||
def om2ax(om):
|
||||
"""Rotation matrix to axis angle pair."""
|
||||
#return qu2ax(om2qu(om)) # HOTFIX
|
||||
ax=np.empty(4)
|
||||
|
||||
# first get the rotation angle
|
||||
|
|
@ -446,11 +445,6 @@ def mul(me, other):
|
|||
other : numpy.ndarray or Rotation
|
||||
Vector, second or fourth order tensor, or rotation object that is rotated.
|
||||
|
||||
Todo
|
||||
----
|
||||
Document details active/passive)
|
||||
consider rotation of (3,3,3,3)-matrix
|
||||
|
||||
"""
|
||||
if me.quaternion.shape != (4,):
|
||||
raise NotImplementedError('Support for multiple rotations missing')
|
||||
|
|
|
|||
|
|
@ -2,11 +2,10 @@
|
|||
#include <stdio.h>
|
||||
#include <unistd.h>
|
||||
#include <dirent.h>
|
||||
#include <sys/types.h>
|
||||
#include <sys/stat.h>
|
||||
#include <stdio.h>
|
||||
#include <string.h>
|
||||
#include <signal.h>
|
||||
#include <sys/types.h>
|
||||
#include <sys/stat.h>
|
||||
|
||||
/* http://stackoverflow.com/questions/30279228/is-there-an-alternative-to-getcwd-in-fortran-2003-2008 */
|
||||
|
||||
|
|
|
|||
|
|
@ -588,7 +588,7 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
|
|||
shearBandingContribution: if(dNeq0(prm%sbVelocity)) then
|
||||
|
||||
BoltzmannRatio = prm%E_sb/(kB*T)
|
||||
call math_eigh33(Mp,eigValues,eigVectors) ! is Mp symmetric by design?
|
||||
call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design?
|
||||
|
||||
do i = 1,6
|
||||
P_sb = 0.5_pReal * math_outer(matmul(eigVectors,sb_sComposition(1:3,i)),&
|
||||
|
|
|
|||
126
src/math.f90
126
src/math.f90
|
|
@ -877,15 +877,14 @@ end function math_sampleGaussVar
|
|||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief eigenvalues and eigenvectors of symmetric matrix
|
||||
! ToDo: has wrong oder of arguments
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine math_eigh(m,w,v,error)
|
||||
subroutine math_eigh(w,v,error,m)
|
||||
|
||||
real(pReal), dimension(:,:), intent(in) :: m !< quadratic matrix to compute eigenvectors and values of
|
||||
real(pReal), dimension(size(m,1)), intent(out) :: w !< eigenvalues
|
||||
real(pReal), dimension(size(m,1),size(m,1)), intent(out) :: v !< eigenvectors
|
||||
|
||||
logical, intent(out) :: error
|
||||
logical, intent(out) :: error
|
||||
|
||||
integer :: ierr
|
||||
real(pReal), dimension(size(m,1)**2) :: work
|
||||
|
||||
|
|
@ -902,9 +901,8 @@ end subroutine math_eigh
|
|||
!> @author Joachim Kopp, Max-Planck-Institut für Kernphysik, Heidelberg (Copyright (C) 2006)
|
||||
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
|
||||
!> @details See http://arxiv.org/abs/physics/0610206 (DSYEVH3)
|
||||
! ToDo: has wrong oder of arguments
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine math_eigh33(m,w,v)
|
||||
subroutine math_eigh33(w,v,m)
|
||||
|
||||
real(pReal), dimension(3,3),intent(in) :: m !< 3x3 matrix to compute eigenvectors and values of
|
||||
real(pReal), dimension(3), intent(out) :: w !< eigenvalues
|
||||
|
|
@ -928,7 +926,7 @@ subroutine math_eigh33(m,w,v)
|
|||
(m(1,1) - w(1)) * (m(2,2) - w(1)) - v(3,2)]
|
||||
norm = norm2(v(1:3, 1))
|
||||
fallback1: if(norm < threshold) then
|
||||
call math_eigh(m,w,v,error)
|
||||
call math_eigh(w,v,error,m)
|
||||
else fallback1
|
||||
v(1:3,1) = v(1:3, 1) / norm
|
||||
v(1:3,2) = [ v(1,2) + m(1, 3) * w(2), &
|
||||
|
|
@ -936,7 +934,7 @@ subroutine math_eigh33(m,w,v)
|
|||
(m(1,1) - w(2)) * (m(2,2) - w(2)) - v(3,2)]
|
||||
norm = norm2(v(1:3, 2))
|
||||
fallback2: if(norm < threshold) then
|
||||
call math_eigh(m,w,v,error)
|
||||
call math_eigh(w,v,error,m)
|
||||
else fallback2
|
||||
v(1:3,2) = v(1:3, 2) / norm
|
||||
v(1:3,3) = math_cross(v(1:3,1),v(1:3,2))
|
||||
|
|
@ -946,87 +944,49 @@ subroutine math_eigh33(m,w,v)
|
|||
end subroutine math_eigh33
|
||||
|
||||
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief rotational part from polar decomposition of 3x3 tensor
|
||||
!> @brief Calculate rotational part of a deformation gradient
|
||||
!> @details https://www.jstor.org/stable/43637254
|
||||
!! https://www.jstor.org/stable/43637372
|
||||
!! https://doi.org/10.1023/A:1007407802076
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
function math_rotationalPart(m)
|
||||
pure function math_rotationalPart(F) result(R)
|
||||
|
||||
real(pReal), intent(in), dimension(3,3) :: m
|
||||
real(pReal), dimension(3,3) :: math_rotationalPart
|
||||
real(pReal), dimension(3,3) :: U , Uinv
|
||||
real(pReal), dimension(3,3), intent(in) :: &
|
||||
F ! deformation gradient
|
||||
real(pReal), dimension(3,3) :: &
|
||||
C, & ! right Cauchy-Green tensor
|
||||
R ! rotational part
|
||||
real(pReal), dimension(3) :: &
|
||||
lambda, & ! principal stretches
|
||||
I_C, & ! invariants of C
|
||||
I_U ! invariants of U
|
||||
real(pReal), dimension(2) :: &
|
||||
I_F ! first two invariants of F
|
||||
real(pReal) :: x,Phi
|
||||
|
||||
U = eigenvectorBasis(matmul(transpose(m),m))
|
||||
Uinv = math_inv33(U)
|
||||
C = matmul(transpose(F),F)
|
||||
I_C = math_invariantsSym33(C)
|
||||
I_F = [math_trace33(F), 0.5*(math_trace33(F)**2 - math_trace33(matmul(F,F)))]
|
||||
|
||||
inversionFailed: if (all(dEq0(Uinv))) then
|
||||
math_rotationalPart = math_I3
|
||||
call IO_warning(650)
|
||||
else inversionFailed
|
||||
math_rotationalPart = matmul(m,Uinv)
|
||||
endif inversionFailed
|
||||
x = math_clip(I_C(1)**2 -3.0_pReal*I_C(2),0.0_pReal)**(3.0_pReal/2.0_pReal)
|
||||
if(dNeq0(x)) then
|
||||
Phi = acos(math_clip((I_C(1)**3 -4.5_pReal*I_C(1)*I_C(2) +13.5_pReal*I_C(3))/x,-1.0_pReal,1.0_pReal))
|
||||
lambda = I_C(1) +(2.0_pReal * sqrt(math_clip(I_C(1)**2-3.0_pReal*I_C(2),0.0_pReal))) &
|
||||
*cos((Phi-2.0_pReal * PI*[1.0_pReal,2.0_pReal,3.0_pReal])/3.0_pReal)
|
||||
lambda = sqrt(math_clip(lambda,0.0_pReal)/3.0_pReal)
|
||||
else
|
||||
lambda = sqrt(I_C(1)/3.0_pReal)
|
||||
endif
|
||||
|
||||
contains
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief eigenvector basis of positive-definite 3x3 matrix
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
pure function eigenvectorBasis(m)
|
||||
I_U = [sum(lambda), lambda(1)*lambda(2)+lambda(2)*lambda(3)+lambda(3)*lambda(1), product(lambda)]
|
||||
|
||||
real(pReal), dimension(3,3) :: eigenvectorBasis
|
||||
real(pReal), dimension(3,3), intent(in) :: m !< positive-definite matrix of which the basis is computed
|
||||
|
||||
real(pReal), dimension(3) :: I, v
|
||||
real(pReal) :: P, Q, rho, phi
|
||||
real(pReal), parameter :: TOL=1.e-14_pReal
|
||||
real(pReal), dimension(3,3,3) :: N, EB
|
||||
|
||||
I = math_invariantsSym33(m)
|
||||
|
||||
P = I(2)-I(1)**2.0_pReal/3.0_pReal
|
||||
Q = -2.0_pReal/27.0_pReal*I(1)**3.0_pReal+product(I(1:2))/3.0_pReal-I(3)
|
||||
|
||||
threeSimilarEigVals: if(all(abs([P,Q]) < TOL)) then
|
||||
v = I(1)/3.0_pReal
|
||||
! this is not really correct, but at least the basis is correct
|
||||
EB = 0.0_pReal
|
||||
EB(1,1,1)=1.0_pReal
|
||||
EB(2,2,2)=1.0_pReal
|
||||
EB(3,3,3)=1.0_pReal
|
||||
else threeSimilarEigVals
|
||||
rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal
|
||||
phi=acos(math_clip(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal))
|
||||
v = 2.0_pReal*rho**(1.0_pReal/3.0_pReal)* [cos((phi )/3.0_pReal), &
|
||||
cos((phi+2.0_pReal*PI)/3.0_pReal), &
|
||||
cos((phi+4.0_pReal*PI)/3.0_pReal) &
|
||||
] + I(1)/3.0_pReal
|
||||
N(1:3,1:3,1) = m-v(1)*math_I3
|
||||
N(1:3,1:3,2) = m-v(2)*math_I3
|
||||
N(1:3,1:3,3) = m-v(3)*math_I3
|
||||
twoSimilarEigVals: if(abs(v(1)-v(2)) < TOL) then
|
||||
EB(1:3,1:3,3) = matmul(N(1:3,1:3,1),N(1:3,1:3,2))/((v(3)-v(1))*(v(3)-v(2)))
|
||||
EB(1:3,1:3,1) = math_I3-EB(1:3,1:3,3)
|
||||
EB(1:3,1:3,2) = 0.0_pReal
|
||||
elseif (abs(v(2)-v(3)) < TOL) then twoSimilarEigVals
|
||||
EB(1:3,1:3,1) = matmul(N(1:3,1:3,2),N(1:3,1:3,3))/((v(1)-v(2))*(v(1)-v(3)))
|
||||
EB(1:3,1:3,2) = math_I3-EB(1:3,1:3,1)
|
||||
EB(1:3,1:3,3) = 0.0_pReal
|
||||
elseif (abs(v(3)-v(1)) < TOL) then twoSimilarEigVals
|
||||
EB(1:3,1:3,2) = matmul(N(1:3,1:3,3),N(1:3,1:3,1))/((v(2)-v(3))*(v(2)-v(1)))
|
||||
EB(1:3,1:3,3) = math_I3-EB(1:3,1:3,2)
|
||||
EB(1:3,1:3,1) = 0.0_pReal
|
||||
else twoSimilarEigVals
|
||||
EB(1:3,1:3,1) = matmul(N(1:3,1:3,2),N(1:3,1:3,3))/((v(1)-v(2))*(v(1)-v(3)))
|
||||
EB(1:3,1:3,2) = matmul(N(1:3,1:3,3),N(1:3,1:3,1))/((v(2)-v(3))*(v(2)-v(1)))
|
||||
EB(1:3,1:3,3) = matmul(N(1:3,1:3,1),N(1:3,1:3,2))/((v(3)-v(1))*(v(3)-v(2)))
|
||||
endif twoSimilarEigVals
|
||||
endif threeSimilarEigVals
|
||||
|
||||
eigenvectorBasis = sqrt(v(1)) * EB(1:3,1:3,1) &
|
||||
+ sqrt(v(2)) * EB(1:3,1:3,2) &
|
||||
+ sqrt(v(3)) * EB(1:3,1:3,3)
|
||||
|
||||
end function eigenvectorBasis
|
||||
R = I_U(1)*I_F(2) * math_I3 &
|
||||
+(I_U(1)**2-I_U(2)) * F &
|
||||
- I_U(1)*I_F(1) * transpose(F) &
|
||||
+ I_U(1) * transpose(matmul(F,F)) &
|
||||
- matmul(F,C)
|
||||
R = R /(I_U(1)*I_U(2)-I_U(3))
|
||||
|
||||
end function math_rotationalPart
|
||||
|
||||
|
|
@ -1078,7 +1038,7 @@ function math_eigvalsh33(m)
|
|||
rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal
|
||||
phi=acos(math_clip(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal))
|
||||
math_eigvalsh33 = 2.0_pReal*rho**(1.0_pReal/3.0_pReal)* &
|
||||
[cos(phi/3.0_pReal), &
|
||||
[cos( phi /3.0_pReal), &
|
||||
cos((phi+2.0_pReal*PI)/3.0_pReal), &
|
||||
cos((phi+4.0_pReal*PI)/3.0_pReal) &
|
||||
] &
|
||||
|
|
|
|||
|
|
@ -34,7 +34,7 @@
|
|||
!> @details: rotation is internally stored as quaternion. It can be inialized from different
|
||||
!> representations and also returns itself in different representations.
|
||||
!
|
||||
! All methods and naming conventions based on Rowenhorst_etal2015
|
||||
! All methods and naming conventions based on Rowenhorst et al. 2015
|
||||
! Convention 1: coordinate frames are right-handed
|
||||
! Convention 2: a rotation angle ω is taken to be positive for a counterclockwise rotation
|
||||
! when viewing from the end point of the rotation axis towards the origin
|
||||
|
|
@ -566,7 +566,26 @@ pure function om2qu(om) result(qu)
|
|||
real(pReal), intent(in), dimension(3,3) :: om
|
||||
real(pReal), dimension(4) :: qu
|
||||
|
||||
qu = eu2qu(om2eu(om))
|
||||
real(pReal) :: trace,s
|
||||
trace = math_trace33(om)
|
||||
|
||||
if(trace > 0.0_pReal) then
|
||||
s = 0.5_pReal / sqrt(trace+1.0_pReal)
|
||||
qu = [0.25_pReal/s, (om(3,2)-om(2,3))*s,(om(1,3)-om(3,1))*s,(om(2,1)-om(1,2))*s]
|
||||
else
|
||||
if( om(1,1) > om(2,2) .and. om(1,1) > om(3,3) ) then
|
||||
s = 2.0_pReal * sqrt( 1.0_pReal + om(1,1) - om(2,2) - om(3,3))
|
||||
qu = [ (om(3,2) - om(2,3)) /s,0.25_pReal * s,(om(1,2) + om(2,1)) / s,(om(1,3) + om(3,1)) / s]
|
||||
elseif (om(2,2) > om(3,3)) then
|
||||
s = 2.0_pReal * sqrt( 1.0_pReal + om(2,2) - om(1,1) - om(3,3))
|
||||
qu = [ (om(1,3) - om(3,1)) /s,(om(1,2) + om(2,1)) / s,0.25_pReal * s,(om(2,3) + om(3,2)) / s]
|
||||
else
|
||||
s = 2.0_pReal * sqrt( 1.0_pReal + om(3,3) - om(1,1) - om(2,2) )
|
||||
qu = [ (om(2,1) - om(1,2)) /s,(om(1,3) + om(3,1)) / s,(om(2,3) + om(3,2)) / s,0.25_pReal * s]
|
||||
endif
|
||||
endif
|
||||
if(qu(1)<0._pReal) qu =-1.0_pReal * qu
|
||||
qu = qu*[1.0_pReal,P,P,P]
|
||||
|
||||
end function om2qu
|
||||
|
||||
|
|
@ -727,7 +746,7 @@ pure function eu2om(eu) result(om)
|
|||
om(3,2) = -c(1)*s(2)
|
||||
om(3,3) = c(2)
|
||||
|
||||
where(dEq0(om)) om = 0.0_pReal
|
||||
where(abs(om)<1.0e-12_pReal) om = 0.0_pReal
|
||||
|
||||
end function eu2om
|
||||
|
||||
|
|
@ -1386,49 +1405,37 @@ subroutine selfTest
|
|||
sin(2.0_pReal*PI*x(1))*A]
|
||||
if(qu(1)<0.0_pReal) qu = qu * (-1.0_pReal)
|
||||
endif
|
||||
#ifndef __PGI
|
||||
if(dNeq0(norm2(om2qu(qu2om(qu))-qu),1.0e-12_pReal)) msg = trim(msg)//'om2qu/qu2om,'
|
||||
if(dNeq0(norm2(eu2qu(qu2eu(qu))-qu),1.0e-12_pReal)) msg = trim(msg)//'eu2qu/qu2eu,'
|
||||
if(dNeq0(norm2(ax2qu(qu2ax(qu))-qu),1.0e-12_pReal)) msg = trim(msg)//'ax2qu/qu2ax,'
|
||||
if(dNeq0(norm2(ro2qu(qu2ro(qu))-qu),1.0e-12_pReal)) msg = trim(msg)//'ro2qu/qu2ro,'
|
||||
if(dNeq0(norm2(ho2qu(qu2ho(qu))-qu),1.0e-7_pReal)) msg = trim(msg)//'ho2qu/qu2ho,'
|
||||
if(dNeq0(norm2(cu2qu(qu2cu(qu))-qu),1.0e-7_pReal)) msg = trim(msg)//'cu2qu/qu2cu,'
|
||||
#endif
|
||||
if(.not. quaternion_equal(om2qu(qu2om(qu)),qu)) msg = trim(msg)//'om2qu/qu2om,'
|
||||
if(.not. quaternion_equal(eu2qu(qu2eu(qu)),qu)) msg = trim(msg)//'eu2qu/qu2eu,'
|
||||
if(.not. quaternion_equal(ax2qu(qu2ax(qu)),qu)) msg = trim(msg)//'ax2qu/qu2ax,'
|
||||
if(.not. quaternion_equal(ro2qu(qu2ro(qu)),qu)) msg = trim(msg)//'ro2qu/qu2ro,'
|
||||
if(.not. quaternion_equal(ho2qu(qu2ho(qu)),qu)) msg = trim(msg)//'ho2qu/qu2ho,'
|
||||
if(.not. quaternion_equal(cu2qu(qu2cu(qu)),qu)) msg = trim(msg)//'cu2qu/qu2cu,'
|
||||
|
||||
om = qu2om(qu)
|
||||
#ifndef __PGI
|
||||
if(dNeq0(norm2(om2qu(eu2om(om2eu(om)))-qu),1.0e-7_pReal)) msg = trim(msg)//'eu2om/om2eu,'
|
||||
if(dNeq0(norm2(om2qu(ax2om(om2ax(om)))-qu),1.0e-7_pReal)) msg = trim(msg)//'ax2om/om2ax,'
|
||||
if(dNeq0(norm2(om2qu(ro2om(om2ro(om)))-qu),1.0e-12_pReal)) msg = trim(msg)//'ro2om/om2ro,'
|
||||
if(dNeq0(norm2(om2qu(ho2om(om2ho(om)))-qu),1.0e-7_pReal)) msg = trim(msg)//'ho2om/om2ho,'
|
||||
if(dNeq0(norm2(om2qu(cu2om(om2cu(om)))-qu),1.0e-7_pReal)) msg = trim(msg)//'cu2om/om2cu,'
|
||||
#endif
|
||||
if(.not. quaternion_equal(om2qu(eu2om(om2eu(om))),qu)) msg = trim(msg)//'eu2om/om2eu,'
|
||||
if(.not. quaternion_equal(om2qu(ax2om(om2ax(om))),qu)) msg = trim(msg)//'ax2om/om2ax,'
|
||||
if(.not. quaternion_equal(om2qu(ro2om(om2ro(om))),qu)) msg = trim(msg)//'ro2om/om2ro,'
|
||||
if(.not. quaternion_equal(om2qu(ho2om(om2ho(om))),qu)) msg = trim(msg)//'ho2om/om2ho,'
|
||||
if(.not. quaternion_equal(om2qu(cu2om(om2cu(om))),qu)) msg = trim(msg)//'cu2om/om2cu,'
|
||||
|
||||
eu = qu2eu(qu)
|
||||
#ifndef __PGI
|
||||
if(dNeq0(norm2(eu2qu(ax2eu(eu2ax(eu)))-qu),1.0e-12_pReal)) msg = trim(msg)//'ax2eu/eu2ax,'
|
||||
if(dNeq0(norm2(eu2qu(ro2eu(eu2ro(eu)))-qu),1.0e-12_pReal)) msg = trim(msg)//'ro2eu/eu2ro,'
|
||||
if(dNeq0(norm2(eu2qu(ho2eu(eu2ho(eu)))-qu),1.0e-7_pReal)) msg = trim(msg)//'ho2eu/eu2ho,'
|
||||
if(dNeq0(norm2(eu2qu(cu2eu(eu2cu(eu)))-qu),1.0e-7_pReal)) msg = trim(msg)//'cu2eu/eu2cu,'
|
||||
#endif
|
||||
if(.not. quaternion_equal(eu2qu(ax2eu(eu2ax(eu))),qu)) msg = trim(msg)//'ax2eu/eu2ax,'
|
||||
if(.not. quaternion_equal(eu2qu(ro2eu(eu2ro(eu))),qu)) msg = trim(msg)//'ro2eu/eu2ro,'
|
||||
if(.not. quaternion_equal(eu2qu(ho2eu(eu2ho(eu))),qu)) msg = trim(msg)//'ho2eu/eu2ho,'
|
||||
if(.not. quaternion_equal(eu2qu(cu2eu(eu2cu(eu))),qu)) msg = trim(msg)//'cu2eu/eu2cu,'
|
||||
|
||||
ax = qu2ax(qu)
|
||||
#ifndef __PGI
|
||||
if(dNeq0(norm2(ax2qu(ro2ax(ax2ro(ax)))-qu),1.0e-12_pReal)) msg = trim(msg)//'ro2ax/ax2ro,'
|
||||
if(dNeq0(norm2(ax2qu(ho2ax(ax2ho(ax)))-qu),1.0e-7_pReal)) msg = trim(msg)//'ho2ax/ax2ho,'
|
||||
if(dNeq0(norm2(ax2qu(cu2ax(ax2cu(ax)))-qu),1.0e-7_pReal)) msg = trim(msg)//'cu2ax/ax2cu,'
|
||||
#endif
|
||||
if(.not. quaternion_equal(ax2qu(ro2ax(ax2ro(ax))),qu)) msg = trim(msg)//'ro2ax/ax2ro,'
|
||||
if(.not. quaternion_equal(ax2qu(ho2ax(ax2ho(ax))),qu)) msg = trim(msg)//'ho2ax/ax2ho,'
|
||||
if(.not. quaternion_equal(ax2qu(cu2ax(ax2cu(ax))),qu)) msg = trim(msg)//'cu2ax/ax2cu,'
|
||||
|
||||
ro = qu2ro(qu)
|
||||
#ifndef __PGI
|
||||
if(dNeq0(norm2(ro2qu(ho2ro(ro2ho(ro)))-qu),1.0e-7_pReal)) msg = trim(msg)//'ho2ro/ro2ho,'
|
||||
if(dNeq0(norm2(ro2qu(cu2ro(ro2cu(ro)))-qu),1.0e-7_pReal)) msg = trim(msg)//'cu2ro/ro2cu,'
|
||||
#endif
|
||||
if(.not. quaternion_equal(ro2qu(ho2ro(ro2ho(ro))),qu)) msg = trim(msg)//'ho2ro/ro2ho,'
|
||||
if(.not. quaternion_equal(ro2qu(cu2ro(ro2cu(ro))),qu)) msg = trim(msg)//'cu2ro/ro2cu,'
|
||||
|
||||
ho = qu2ho(qu)
|
||||
#ifndef __PGI
|
||||
if(dNeq0(norm2(ho2qu(cu2ho(ho2cu(ho)))-qu),1.0e-7_pReal)) msg = trim(msg)//'cu2ho/ho2cu,'
|
||||
#endif
|
||||
if(.not. quaternion_equal(ho2qu(cu2ho(ho2cu(ho))),qu)) msg = trim(msg)//'cu2ho/ho2cu,'
|
||||
|
||||
call R%fromMatrix(om)
|
||||
|
||||
|
|
@ -1447,6 +1454,18 @@ subroutine selfTest
|
|||
if(len_trim(msg) /= 0) call IO_error(0,ext_msg=msg)
|
||||
|
||||
enddo
|
||||
contains
|
||||
|
||||
function quaternion_equal(qu1,qu2) result(ok)
|
||||
|
||||
real(pReal), intent(in), dimension(4) :: qu1,qu2
|
||||
logical :: ok
|
||||
|
||||
ok = all(dEq(qu1,qu2,1.0e-7_pReal))
|
||||
if(dEq0(qu1(1),1.0e-12_pReal)) &
|
||||
ok = ok .or. all(dEq(-1.0_pReal*qu1,qu2,1.0e-7_pReal))
|
||||
|
||||
end function quaternion_equal
|
||||
|
||||
end subroutine selfTest
|
||||
|
||||
|
|
|
|||
|
|
@ -29,7 +29,7 @@ module system_routines
|
|||
|
||||
integer(C_INT) :: isDirectory_C
|
||||
character(kind=C_CHAR), dimension(pPathLen), intent(in) :: path ! C string is an array
|
||||
end function isDirectory_C
|
||||
end function isDirectory_C
|
||||
|
||||
subroutine getCurrentWorkDir_C(path, stat) bind(C)
|
||||
use, intrinsic :: ISO_C_Binding, only: &
|
||||
|
|
@ -40,7 +40,7 @@ module system_routines
|
|||
|
||||
character(kind=C_CHAR), dimension(pPathLen), intent(out) :: path ! C string is an array
|
||||
integer(C_INT), intent(out) :: stat
|
||||
end subroutine getCurrentWorkDir_C
|
||||
end subroutine getCurrentWorkDir_C
|
||||
|
||||
subroutine getHostName_C(str, stat) bind(C)
|
||||
use, intrinsic :: ISO_C_Binding, only: &
|
||||
|
|
@ -51,7 +51,7 @@ module system_routines
|
|||
|
||||
character(kind=C_CHAR), dimension(pStringLen), intent(out) :: str ! C string is an array
|
||||
integer(C_INT), intent(out) :: stat
|
||||
end subroutine getHostName_C
|
||||
end subroutine getHostName_C
|
||||
|
||||
function chdir_C(path) bind(C)
|
||||
use, intrinsic :: ISO_C_Binding, only: &
|
||||
|
|
|
|||
Loading…
Reference in New Issue