2012-07-19 00:09:59 +05:30
!--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @author Christoph Kords, Max-Planck-Institut für Eisenforschung GmbH
2016-01-12 16:30:23 +05:30
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
2012-07-19 00:09:59 +05:30
!> @brief Mathematical library, including random number generation and tensor represenations
!--------------------------------------------------------------------------------------------------
2012-08-25 17:16:36 +05:30
module math
2013-01-31 21:58:08 +05:30
use prec , only : &
pReal , &
pInt
2012-08-25 17:16:36 +05:30
2007-03-20 19:25:22 +05:30
implicit none
2013-06-06 14:55:27 +05:30
private
2016-05-25 11:22:56 +05:30
real ( pReal ) , parameter , public :: PI = 3.141592653589793_pReal !< ratio of a circle's circumference to its diameter
2013-03-07 01:04:30 +05:30
real ( pReal ) , parameter , public :: INDEG = 18 0.0_pReal / PI !< conversion from radian into degree
real ( pReal ) , parameter , public :: INRAD = PI / 18 0.0_pReal !< conversion from degree into radian
2016-05-25 11:22:56 +05:30
complex ( pReal ) , parameter , public :: TWOPIIMG = ( 0.0_pReal , 2.0_pReal ) * ( PI , 0.0_pReal ) !< Re(0.0), Im(2xPi)
2011-12-01 17:31:13 +05:30
2012-03-30 01:24:31 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , parameter , public :: &
2013-09-19 01:24:39 +05:30
MATH_I3 = reshape ( [ &
2012-03-30 01:24:31 +05:30
1.0_pReal , 0.0_pReal , 0.0_pReal , &
0.0_pReal , 1.0_pReal , 0.0_pReal , &
0.0_pReal , 0.0_pReal , 1.0_pReal &
2012-07-19 00:09:59 +05:30
] , [ 3 , 3 ] ) !< 3x3 Identity
2008-02-15 18:12:27 +05:30
2012-03-30 01:24:31 +05:30
integer ( pInt ) , dimension ( 2 , 6 ) , parameter , private :: &
mapMandel = reshape ( [ &
1_pInt , 1_pInt , &
2_pInt , 2_pInt , &
3_pInt , 3_pInt , &
1_pInt , 2_pInt , &
2_pInt , 3_pInt , &
1_pInt , 3_pInt &
2012-08-25 17:16:36 +05:30
] , [ 2 , 6 ] ) !< arrangement in Mandel notation
2012-03-30 01:24:31 +05:30
real ( pReal ) , dimension ( 6 ) , parameter , private :: &
nrmMandel = [ &
1.0_pReal , 1.0_pReal , 1.0_pReal , &
2013-03-05 20:05:26 +05:30
1.414213562373095_pReal , 1.414213562373095_pReal , 1.414213562373095_pReal ] !< weighting for Mandel notation (forward)
2012-08-25 17:16:36 +05:30
2012-03-30 01:24:31 +05:30
real ( pReal ) , dimension ( 6 ) , parameter , public :: &
invnrmMandel = [ &
1.0_pReal , 1.0_pReal , 1.0_pReal , &
2013-03-05 20:05:26 +05:30
0.7071067811865476_pReal , 0.7071067811865476_pReal , 0.7071067811865476_pReal ] !< weighting for Mandel notation (backward)
2008-02-15 18:12:27 +05:30
2012-03-30 01:24:31 +05:30
integer ( pInt ) , dimension ( 2 , 6 ) , parameter , private :: &
mapVoigt = reshape ( [ &
1_pInt , 1_pInt , &
2_pInt , 2_pInt , &
3_pInt , 3_pInt , &
2_pInt , 3_pInt , &
1_pInt , 3_pInt , &
1_pInt , 2_pInt &
2012-08-25 17:16:36 +05:30
] , [ 2 , 6 ] ) !< arrangement in Voigt notation
2012-03-30 01:24:31 +05:30
2012-08-25 17:16:36 +05:30
real ( pReal ) , dimension ( 6 ) , parameter , private :: &
nrmVoigt = 1.0_pReal , & !< weighting for Voigt notation (forward)
invnrmVoigt = 1.0_pReal !< weighting for Voigt notation (backward)
2008-02-15 18:12:27 +05:30
2012-03-30 01:24:31 +05:30
integer ( pInt ) , dimension ( 2 , 9 ) , parameter , private :: &
mapPlain = reshape ( [ &
1_pInt , 1_pInt , &
1_pInt , 2_pInt , &
1_pInt , 3_pInt , &
2_pInt , 1_pInt , &
2_pInt , 2_pInt , &
2_pInt , 3_pInt , &
3_pInt , 1_pInt , &
3_pInt , 2_pInt , &
3_pInt , 3_pInt &
2012-08-25 17:16:36 +05:30
] , [ 2 , 9 ] ) !< arrangement in Plain notation
2008-02-15 18:12:27 +05:30
2013-01-31 21:58:08 +05:30
public :: &
math_init , &
2013-02-01 12:56:21 +05:30
math_qsort , &
2013-01-31 21:58:08 +05:30
math_range , &
math_identity2nd , &
math_identity4th , &
math_civita , &
math_delta , &
2016-01-10 19:04:26 +05:30
math_crossproduct , &
math_tensorproduct33 , &
2013-01-31 21:58:08 +05:30
math_mul3x3 , &
math_mul6x6 , &
math_mul33xx33 , &
math_mul3333xx33 , &
math_mul3333xx3333 , &
math_mul33x33 , &
math_mul66x66 , &
math_mul99x99 , &
math_mul33x3 , &
math_mul33x3_complex , &
math_mul66x6 , &
math_exp33 , &
math_transpose33 , &
math_inv33 , &
math_invert33 , &
math_invSym3333 , &
math_invert , &
math_symmetric33 , &
math_symmetric66 , &
math_skew33 , &
2016-01-09 00:27:37 +05:30
math_spherical33 , &
2013-01-31 21:58:08 +05:30
math_deviatoric33 , &
math_equivStrain33 , &
2015-01-13 15:13:05 +05:30
math_equivStress33 , &
2013-01-31 21:58:08 +05:30
math_trace33 , &
math_det33 , &
math_Plain33to9 , &
math_Plain9to33 , &
math_Mandel33to6 , &
math_Mandel6to33 , &
math_Plain3333to99 , &
math_Plain99to3333 , &
math_Mandel66toPlain66 , &
math_Plain66toMandel66 , &
math_Mandel3333to66 , &
math_Mandel66to3333 , &
math_Voigt66to3333 , &
math_qRand , &
math_qMul , &
math_qDot , &
math_qConj , &
math_qInv , &
math_qRot , &
math_RtoEuler , &
math_RtoQ , &
math_EulerToR , &
math_EulerToQ , &
2013-06-06 00:40:37 +05:30
math_EulerAxisAngleToR , &
math_axisAngleToR , &
math_EulerAxisAngleToQ , &
math_axisAngleToQ , &
2013-01-31 21:58:08 +05:30
math_qToRodrig , &
math_qToEuler , &
2013-06-06 00:40:37 +05:30
math_qToEulerAxisAngle , &
2013-01-31 21:58:08 +05:30
math_qToAxisAngle , &
math_qToR , &
math_EulerMisorientation , &
math_sampleRandomOri , &
math_sampleGaussOri , &
math_sampleFiberOri , &
math_sampleGaussVar , &
math_symmetricEulers , &
2016-02-27 00:51:47 +05:30
math_eigenvectorBasisSym33 , &
math_eigenvectorBasisSym , &
2016-02-26 20:06:24 +05:30
math_eigenValuesVectorsSym33 , &
math_eigenValuesVectorsSym , &
2016-01-13 03:09:31 +05:30
math_rotationalPart33 , &
2016-02-02 17:53:45 +05:30
math_invariantsSym33 , &
2016-01-12 16:30:23 +05:30
math_eigenvaluesSym33 , &
2014-08-21 14:03:55 +05:30
math_factorial , &
math_binomial , &
math_multinomial , &
2013-01-31 21:58:08 +05:30
math_volTetrahedron , &
2013-04-09 23:37:30 +05:30
math_areaTriangle , &
2013-01-31 21:58:08 +05:30
math_rotate_forward33 , &
math_rotate_backward33 , &
2016-07-27 18:02:03 +05:30
math_rotate_forward3333 , &
math_limit
2013-01-31 21:58:08 +05:30
private :: &
halton , &
halton_memory , &
halton_ndim_set , &
2017-08-12 09:33:40 +05:30
halton_seed_set
2013-01-31 21:58:08 +05:30
2012-03-09 01:55:28 +05:30
contains
2012-08-25 17:16:36 +05:30
2012-07-19 00:09:59 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief initialization of random seed generator
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
subroutine math_init
2007-03-28 12:51:47 +05:30
2012-07-19 00:09:59 +05:30
use , intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
2011-12-01 17:31:13 +05:30
use prec , only : tol_math_check
2014-10-10 01:53:06 +05:30
use numerics , only : &
fixedSeed
2013-02-25 22:04:59 +05:30
use IO , only : IO_error , IO_timeStamp
2012-08-25 17:16:36 +05:30
2007-03-28 12:51:47 +05:30
implicit none
2011-11-04 15:59:35 +05:30
integer ( pInt ) :: i
2010-05-06 19:37:21 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: R , R2
2013-06-06 00:40:37 +05:30
real ( pReal ) , dimension ( 3 ) :: Eulers , v
2011-11-04 15:59:35 +05:30
real ( pReal ) , dimension ( 4 ) :: q , q2 , axisangle , randTest
2011-12-01 17:31:13 +05:30
! the following variables are system dependend and shound NOT be pInt
2012-08-25 17:16:36 +05:30
integer :: randSize ! gfortran requires a variable length to compile
2012-07-19 00:09:59 +05:30
integer , dimension ( : ) , allocatable :: randInit ! if recalculations of former randomness (with given seed) is necessary
! comment the first random_seed call out, set randSize to 1, and use ifort
2012-01-20 02:08:52 +05:30
character ( len = 64 ) :: error_msg
2012-08-25 17:16:36 +05:30
2016-07-09 01:42:02 +05:30
write ( 6 , '(/,a)' ) ' <<<+- math init -+>>>'
write ( 6 , '(a15,a)' ) ' Current time: ' , IO_timeStamp ( )
2012-02-01 00:48:55 +05:30
#include "compilation_info.f90"
2012-08-25 17:16:36 +05:30
2011-11-04 15:59:35 +05:30
call random_seed ( size = randSize )
2012-08-25 17:16:36 +05:30
if ( allocated ( randInit ) ) deallocate ( randInit )
2011-11-04 15:59:35 +05:30
allocate ( randInit ( randSize ) )
2012-03-09 01:55:28 +05:30
if ( fixedSeed > 0_pInt ) then
2012-08-30 20:28:30 +05:30
randInit ( 1 : randSize ) = int ( fixedSeed ) ! fixedSeed is of type pInt, randInit not
2009-08-27 21:00:40 +05:30
call random_seed ( put = randInit )
else
call random_seed ( )
2012-08-30 20:28:30 +05:30
call random_seed ( get = randInit )
randInit ( 2 : randSize ) = randInit ( 1 )
call random_seed ( put = randInit )
2009-08-27 21:00:40 +05:30
endif
2011-12-01 17:31:13 +05:30
do i = 1_pInt , 4_pInt
2011-11-04 15:59:35 +05:30
call random_number ( randTest ( i ) )
enddo
2016-07-14 18:59:53 +05:30
write ( 6 , '(a,I2)' ) ' size of random seed: ' , randSize
2016-07-09 01:42:02 +05:30
do i = 1 , randSize
2016-07-14 18:59:53 +05:30
write ( 6 , '(a,I2,I14)' ) ' value of random seed: ' , i , randInit ( i )
2016-07-09 01:42:02 +05:30
enddo
write ( 6 , '(a,4(/,26x,f17.14),/)' ) ' start of random sequence: ' , randTest
2012-08-25 17:16:36 +05:30
2012-08-30 20:28:30 +05:30
call random_seed ( put = randInit )
2011-11-04 15:59:35 +05:30
2012-02-14 19:13:36 +05:30
call halton_seed_set ( int ( randInit ( 1 ) , pInt ) )
2011-12-01 17:31:13 +05:30
call halton_ndim_set ( 3_pInt )
2010-05-04 18:33:35 +05:30
2010-05-06 19:37:21 +05:30
! --- check rotation dictionary ---
2013-06-06 00:40:37 +05:30
q = math_qRand ( ) ! random quaternion
2010-05-06 19:37:21 +05:30
! +++ q -> a -> q +++
2013-06-06 00:40:37 +05:30
axisangle = math_qToAxisAngle ( q )
q2 = math_axisAngleToQ ( axisangle ( 1 : 3 ) , axisangle ( 4 ) )
2011-12-01 17:31:13 +05:30
if ( any ( abs ( q - q2 ) > tol_math_check ) . and . &
2012-01-20 02:08:52 +05:30
any ( abs ( - q - q2 ) > tol_math_check ) ) then
2012-02-01 00:48:55 +05:30
write ( error_msg , '(a,e14.6)' ) 'maximum deviation ' , min ( maxval ( abs ( q - q2 ) ) , maxval ( abs ( - q - q2 ) ) )
2012-02-13 23:11:27 +05:30
call IO_error ( 401_pInt , ext_msg = error_msg )
2012-08-25 17:16:36 +05:30
endif
2010-05-06 19:37:21 +05:30
! +++ q -> R -> q +++
2013-06-06 00:40:37 +05:30
R = math_qToR ( q )
q2 = math_RtoQ ( R )
2011-12-01 17:31:13 +05:30
if ( any ( abs ( q - q2 ) > tol_math_check ) . and . &
2012-01-20 02:08:52 +05:30
any ( abs ( - q - q2 ) > tol_math_check ) ) then
2012-02-01 00:48:55 +05:30
write ( error_msg , '(a,e14.6)' ) 'maximum deviation ' , min ( maxval ( abs ( q - q2 ) ) , maxval ( abs ( - q - q2 ) ) )
2012-02-13 23:11:27 +05:30
call IO_error ( 402_pInt , ext_msg = error_msg )
2012-08-25 17:16:36 +05:30
endif
2010-05-06 19:37:21 +05:30
! +++ q -> euler -> q +++
2013-06-06 00:40:37 +05:30
Eulers = math_qToEuler ( q )
2013-01-31 21:58:08 +05:30
q2 = math_EulerToQ ( Eulers )
2011-12-01 17:31:13 +05:30
if ( any ( abs ( q - q2 ) > tol_math_check ) . and . &
2012-01-20 02:08:52 +05:30
any ( abs ( - q - q2 ) > tol_math_check ) ) then
2012-02-01 00:48:55 +05:30
write ( error_msg , '(a,e14.6)' ) 'maximum deviation ' , min ( maxval ( abs ( q - q2 ) ) , maxval ( abs ( - q - q2 ) ) )
2012-02-13 23:11:27 +05:30
call IO_error ( 403_pInt , ext_msg = error_msg )
2012-08-25 17:16:36 +05:30
endif
2012-01-20 02:08:52 +05:30
2010-05-06 19:37:21 +05:30
! +++ R -> euler -> R +++
2013-06-06 00:40:37 +05:30
Eulers = math_RtoEuler ( R )
2010-05-06 19:37:21 +05:30
R2 = math_EulerToR ( Eulers )
2012-01-20 02:08:52 +05:30
if ( any ( abs ( R - R2 ) > tol_math_check ) ) then
2012-02-01 00:48:55 +05:30
write ( error_msg , '(a,e14.6)' ) 'maximum deviation ' , maxval ( abs ( R - R2 ) )
2012-02-13 23:11:27 +05:30
call IO_error ( 404_pInt , ext_msg = error_msg )
2012-08-25 17:16:36 +05:30
endif
2013-06-06 00:40:37 +05:30
! +++ check rotation sense of q and R +++
q = math_qRand ( ) ! random quaternion
call halton ( 3_pInt , v ) ! random vector
R = math_qToR ( q )
if ( any ( abs ( math_mul33x3 ( R , v ) - math_qRot ( q , v ) ) > tol_math_check ) ) then
2013-06-06 14:55:27 +05:30
write ( 6 , '(a,4(f8.3,1x))' ) 'q' , q
2013-06-06 00:40:37 +05:30
call IO_error ( 409_pInt )
endif
2012-03-09 01:55:28 +05:30
end subroutine math_init
2007-03-28 12:51:47 +05:30
2009-01-26 18:28:58 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Quicksort algorithm for two-dimensional integer arrays
2007-04-03 13:47:58 +05:30
! Sorting is done with respect to array(1,:)
! and keeps array(2:N,:) linked to it.
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2013-02-01 12:56:21 +05:30
recursive subroutine math_qsort ( a , istart , iend )
2007-04-03 13:47:58 +05:30
implicit none
2012-02-09 21:28:15 +05:30
integer ( pInt ) , dimension ( : , : ) , intent ( inout ) :: a
integer ( pInt ) , intent ( in ) :: istart , iend
integer ( pInt ) :: ipivot
2007-04-03 13:47:58 +05:30
if ( istart < iend ) then
2017-08-12 09:33:40 +05:30
ipivot = qsort_partition ( a , istart , iend )
2013-02-01 12:56:21 +05:30
call math_qsort ( a , istart , ipivot - 1_pInt )
call math_qsort ( a , ipivot + 1_pInt , iend )
2007-04-03 13:47:58 +05:30
endif
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2017-08-12 09:33:40 +05:30
contains
2007-04-03 13:47:58 +05:30
2017-08-12 09:33:40 +05:30
!-------------------------------------------------------------------------------------------------
!> @brief Partitioning required for quicksort
!-------------------------------------------------------------------------------------------------
integer ( pInt ) function qsort_partition ( a , istart , iend )
implicit none
integer ( pInt ) , dimension ( : , : ) , intent ( inout ) :: a
integer ( pInt ) , intent ( in ) :: istart , iend
integer ( pInt ) :: i , j , k , tmp
do
! find the first element on the right side less than or equal to the pivot point
do j = iend , istart , - 1_pInt
if ( a ( 1 , j ) < = a ( 1 , istart ) ) exit
2007-04-03 13:47:58 +05:30
enddo
2017-08-12 09:33:40 +05:30
! find the first element on the left side greater than the pivot point
do i = istart , iend
if ( a ( 1 , i ) > a ( 1 , istart ) ) exit
2007-04-03 13:47:58 +05:30
enddo
2017-08-12 09:33:40 +05:30
if ( i < j ) then ! if the indexes do not cross, exchange values
do k = 1_pInt , int ( size ( a , 1_pInt ) , pInt )
tmp = a ( k , i )
a ( k , i ) = a ( k , j )
a ( k , j ) = tmp
enddo
else ! if they do cross, exchange left value with pivot and return with the partition index
do k = 1_pInt , int ( size ( a , 1_pInt ) , pInt )
tmp = a ( k , istart )
a ( k , istart ) = a ( k , j )
a ( k , j ) = tmp
enddo
qsort_partition = j
return
endif
enddo
end function qsort_partition
2007-04-03 13:47:58 +05:30
2017-08-12 09:33:40 +05:30
end subroutine math_qsort
2007-04-03 13:47:58 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief range of integers starting at one
!--------------------------------------------------------------------------------------------------
pure function math_range ( N )
2009-03-04 17:18:54 +05:30
implicit none
2013-03-05 20:05:26 +05:30
integer ( pInt ) , intent ( in ) :: N !< length of range
2011-12-01 17:31:13 +05:30
integer ( pInt ) :: i
2009-03-04 17:18:54 +05:30
integer ( pInt ) , dimension ( N ) :: math_range
2013-12-28 01:33:28 +05:30
math_range = [ ( i , i = 1_pInt , N ) ]
2009-03-04 17:18:54 +05:30
2012-03-09 01:55:28 +05:30
end function math_range
2009-03-04 17:18:54 +05:30
2011-12-01 17:31:13 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief second rank identity tensor of specified dimension
!--------------------------------------------------------------------------------------------------
pure function math_identity2nd ( dimen )
2007-03-29 21:02:52 +05:30
implicit none
2013-03-05 20:05:26 +05:30
integer ( pInt ) , intent ( in ) :: dimen !< tensor dimension
2011-12-01 17:31:13 +05:30
integer ( pInt ) :: i
2007-04-11 15:34:22 +05:30
real ( pReal ) , dimension ( dimen , dimen ) :: math_identity2nd
2007-03-29 21:02:52 +05:30
2012-08-25 17:16:36 +05:30
math_identity2nd = 0.0_pReal
forall ( i = 1_pInt : dimen ) math_identity2nd ( i , i ) = 1.0_pReal
2007-03-29 21:02:52 +05:30
2012-03-09 01:55:28 +05:30
end function math_identity2nd
2007-03-29 21:02:52 +05:30
2013-01-31 21:58:08 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief symmetric fourth rank identity tensor of specified dimension
! from http://en.wikipedia.org/wiki/Tensor_derivative_(continuum_mechanics)#Derivative_of_a_second-order_tensor_with_respect_to_itself
!--------------------------------------------------------------------------------------------------
pure function math_identity4th ( dimen )
implicit none
2013-03-05 20:05:26 +05:30
integer ( pInt ) , intent ( in ) :: dimen !< tensor dimension
2013-01-31 21:58:08 +05:30
integer ( pInt ) :: i , j , k , l
real ( pReal ) , dimension ( dimen , dimen , dimen , dimen ) :: math_identity4th
forall ( i = 1_pInt : dimen , j = 1_pInt : dimen , k = 1_pInt : dimen , l = 1_pInt : dimen ) math_identity4th ( i , j , k , l ) = &
0.5_pReal * ( math_I3 ( i , k ) * math_I3 ( j , l ) + math_I3 ( i , l ) * math_I3 ( j , k ) )
end function math_identity4th
2011-12-01 17:31:13 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief permutation tensor e_ijk used for computing cross product of two tensors
2008-03-26 19:05:01 +05:30
! e_ijk = 1 if even permutation of ijk
! e_ijk = -1 if odd permutation of ijk
! e_ijk = 0 otherwise
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2016-01-10 19:04:26 +05:30
real ( pReal ) pure function math_civita ( i , j , k )
2008-03-26 19:05:01 +05:30
implicit none
2009-12-14 16:32:10 +05:30
integer ( pInt ) , intent ( in ) :: i , j , k
2008-03-26 19:05:01 +05:30
2009-07-31 17:32:20 +05:30
math_civita = 0.0_pReal
2011-12-01 17:31:13 +05:30
if ( ( ( i == 1_pInt ) . and . ( j == 2_pInt ) . and . ( k == 3_pInt ) ) . or . &
( ( i == 2_pInt ) . and . ( j == 3_pInt ) . and . ( k == 1_pInt ) ) . or . &
( ( i == 3_pInt ) . and . ( j == 1_pInt ) . and . ( k == 2_pInt ) ) ) math_civita = 1.0_pReal
if ( ( ( i == 1_pInt ) . and . ( j == 3_pInt ) . and . ( k == 2_pInt ) ) . or . &
( ( i == 2_pInt ) . and . ( j == 1_pInt ) . and . ( k == 3_pInt ) ) . or . &
( ( i == 3_pInt ) . and . ( j == 2_pInt ) . and . ( k == 1_pInt ) ) ) math_civita = - 1.0_pReal
2008-03-27 17:24:34 +05:30
2012-03-09 01:55:28 +05:30
end function math_civita
2008-03-27 17:24:34 +05:30
2011-12-01 17:31:13 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief kronecker delta function d_ij
2008-03-27 17:24:34 +05:30
! d_ij = 1 if i = j
! d_ij = 0 otherwise
2015-05-05 12:07:59 +05:30
! inspired by http://fortraninacworld.blogspot.de/2012/12/ternary-operator.html
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2013-01-31 21:58:08 +05:30
real ( pReal ) pure function math_delta ( i , j )
2008-03-27 17:24:34 +05:30
implicit none
2009-12-14 16:32:10 +05:30
integer ( pInt ) , intent ( in ) :: i , j
2008-03-27 17:24:34 +05:30
2015-05-05 12:07:59 +05:30
math_delta = merge ( 0.0_pReal , 1.0_pReal , i / = j )
2008-03-27 17:24:34 +05:30
2012-03-09 01:55:28 +05:30
end function math_delta
2007-03-29 21:02:52 +05:30
2011-12-01 17:31:13 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2016-01-10 19:04:26 +05:30
!> @brief cross product a x b
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2016-01-10 19:04:26 +05:30
pure function math_crossproduct ( A , B )
2009-01-20 00:40:58 +05:30
implicit none
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: A , B
2016-01-10 19:04:26 +05:30
real ( pReal ) , dimension ( 3 ) :: math_crossproduct
2009-01-20 00:40:58 +05:30
2016-01-10 19:04:26 +05:30
math_crossproduct = [ A ( 2 ) * B ( 3 ) - A ( 3 ) * B ( 2 ) , &
A ( 3 ) * B ( 1 ) - A ( 1 ) * B ( 3 ) , &
A ( 1 ) * B ( 2 ) - A ( 2 ) * B ( 1 ) ]
2009-01-20 00:40:58 +05:30
2016-01-10 19:04:26 +05:30
end function math_crossproduct
2009-01-20 00:40:58 +05:30
2009-03-05 20:07:59 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2016-02-26 20:06:24 +05:30
!> @brief tensor product A \otimes B of arbitrary sized vectors A and B
!--------------------------------------------------------------------------------------------------
pure function math_tensorproduct ( A , B )
implicit none
real ( pReal ) , dimension ( : ) , intent ( in ) :: A , B
real ( pReal ) , dimension ( size ( A , 1 ) , size ( B , 1 ) ) :: math_tensorproduct
integer ( pInt ) :: i , j
forall ( i = 1_pInt : size ( A , 1 ) , j = 1_pInt : size ( B , 1 ) ) math_tensorproduct ( i , j ) = A ( i ) * B ( j )
end function math_tensorproduct
!--------------------------------------------------------------------------------------------------
!> @brief tensor product A \otimes B of leght-3 vectors A and B
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2016-01-10 19:04:26 +05:30
pure function math_tensorproduct33 ( A , B )
2009-03-17 20:43:17 +05:30
implicit none
2016-01-10 19:04:26 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: math_tensorproduct33
2013-01-31 21:58:08 +05:30
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: A , B
2011-12-01 17:31:13 +05:30
integer ( pInt ) :: i , j
2012-08-25 17:16:36 +05:30
2016-01-10 19:04:26 +05:30
forall ( i = 1_pInt : 3_pInt , j = 1_pInt : 3_pInt ) math_tensorproduct33 ( i , j ) = A ( i ) * B ( j )
2009-03-17 20:43:17 +05:30
2016-01-10 19:04:26 +05:30
end function math_tensorproduct33
2009-03-17 20:43:17 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief matrix multiplication 3x3 = 1
!--------------------------------------------------------------------------------------------------
2013-01-31 21:58:08 +05:30
real ( pReal ) pure function math_mul3x3 ( A , B )
2009-03-05 20:07:59 +05:30
implicit none
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: A , B
2013-01-31 21:58:08 +05:30
math_mul3x3 = sum ( A * B )
2009-03-05 20:07:59 +05:30
2012-03-09 01:55:28 +05:30
end function math_mul3x3
2009-03-05 20:07:59 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief matrix multiplication 6x6 = 1
!--------------------------------------------------------------------------------------------------
2013-01-31 21:58:08 +05:30
real ( pReal ) pure function math_mul6x6 ( A , B )
2009-03-05 20:07:59 +05:30
implicit none
real ( pReal ) , dimension ( 6 ) , intent ( in ) :: A , B
2013-01-31 21:58:08 +05:30
math_mul6x6 = sum ( A * B )
2009-03-05 20:07:59 +05:30
2012-03-09 01:55:28 +05:30
end function math_mul6x6
2009-01-20 00:40:58 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2013-08-07 22:50:05 +05:30
!> @brief matrix multiplication 33xx33 = 1 (double contraction --> ij * ij)
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2013-01-31 21:58:08 +05:30
real ( pReal ) pure function math_mul33xx33 ( A , B )
2010-09-30 15:02:49 +05:30
implicit none
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: A , B
2013-01-31 21:58:08 +05:30
integer ( pInt ) :: i , j
2010-09-30 15:02:49 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: C
2011-12-01 17:31:13 +05:30
forall ( i = 1_pInt : 3_pInt , j = 1_pInt : 3_pInt ) C ( i , j ) = A ( i , j ) * B ( i , j )
2010-09-30 15:02:49 +05:30
math_mul33xx33 = sum ( C )
2012-03-09 01:55:28 +05:30
end function math_mul33xx33
2010-09-30 15:02:49 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief matrix multiplication 3333x33 = 33 (double contraction --> ijkl *kl = ij)
!--------------------------------------------------------------------------------------------------
pure function math_mul3333xx33 ( A , B )
2010-10-13 21:34:44 +05:30
implicit none
2013-01-31 21:58:08 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: math_mul3333xx33
2010-10-13 21:34:44 +05:30
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) , intent ( in ) :: A
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: B
2013-01-31 21:58:08 +05:30
integer ( pInt ) :: i , j
2010-10-13 21:34:44 +05:30
2013-01-31 21:58:08 +05:30
forall ( i = 1_pInt : 3_pInt , j = 1_pInt : 3_pInt ) &
2012-02-13 19:38:07 +05:30
math_mul3333xx33 ( i , j ) = sum ( A ( i , j , 1 : 3 , 1 : 3 ) * B ( 1 : 3 , 1 : 3 ) )
2012-03-09 01:55:28 +05:30
end function math_mul3333xx33
2010-10-13 21:34:44 +05:30
2010-09-30 15:02:49 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief matrix multiplication 3333x3333 = 3333 (ijkl *klmn = ijmn)
!--------------------------------------------------------------------------------------------------
pure function math_mul3333xx3333 ( A , B )
2012-02-23 01:41:09 +05:30
implicit none
integer ( pInt ) :: i , j , k , l
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) , intent ( in ) :: A
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) , intent ( in ) :: B
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: math_mul3333xx3333
2016-01-10 19:04:26 +05:30
forall ( i = 1_pInt : 3_pInt , j = 1_pInt : 3_pInt , k = 1_pInt : 3_pInt , l = 1_pInt : 3_pInt ) &
math_mul3333xx3333 ( i , j , k , l ) = sum ( A ( i , j , 1 : 3 , 1 : 3 ) * B ( 1 : 3 , 1 : 3 , k , l ) )
2012-02-23 01:41:09 +05:30
2012-03-09 01:55:28 +05:30
end function math_mul3333xx3333
2012-02-23 01:41:09 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief matrix multiplication 33x33 = 33
!--------------------------------------------------------------------------------------------------
pure function math_mul33x33 ( A , B )
2009-01-20 00:40:58 +05:30
implicit none
real ( pReal ) , dimension ( 3 , 3 ) :: math_mul33x33
2013-01-31 21:58:08 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: A , B
integer ( pInt ) :: i , j
2009-01-20 00:40:58 +05:30
2016-01-10 19:04:26 +05:30
forall ( i = 1_pInt : 3_pInt , j = 1_pInt : 3_pInt ) &
math_mul33x33 ( i , j ) = A ( i , 1 ) * B ( 1 , j ) + A ( i , 2 ) * B ( 2 , j ) + A ( i , 3 ) * B ( 3 , j )
2009-01-20 00:40:58 +05:30
2012-03-09 01:55:28 +05:30
end function math_mul33x33
2009-01-20 00:40:58 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief matrix multiplication 66x66 = 66
!--------------------------------------------------------------------------------------------------
pure function math_mul66x66 ( A , B )
2008-07-09 01:08:22 +05:30
implicit none
2009-01-20 00:40:58 +05:30
real ( pReal ) , dimension ( 6 , 6 ) :: math_mul66x66
2013-01-31 21:58:08 +05:30
real ( pReal ) , dimension ( 6 , 6 ) , intent ( in ) :: A , B
integer ( pInt ) :: i , j
2008-07-09 01:08:22 +05:30
2011-12-01 17:31:13 +05:30
forall ( i = 1_pInt : 6_pInt , j = 1_pInt : 6_pInt ) math_mul66x66 ( i , j ) = &
2008-07-09 01:08:22 +05:30
A ( i , 1 ) * B ( 1 , j ) + A ( i , 2 ) * B ( 2 , j ) + A ( i , 3 ) * B ( 3 , j ) + &
A ( i , 4 ) * B ( 4 , j ) + A ( i , 5 ) * B ( 5 , j ) + A ( i , 6 ) * B ( 6 , j )
2012-03-09 01:55:28 +05:30
end function math_mul66x66
2008-07-09 01:08:22 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief matrix multiplication 99x99 = 99
!--------------------------------------------------------------------------------------------------
pure function math_mul99x99 ( A , B )
2009-08-11 22:01:57 +05:30
2012-03-30 01:24:31 +05:30
implicit none
2009-08-11 22:01:57 +05:30
real ( pReal ) , dimension ( 9 , 9 ) :: math_mul99x99
2013-01-31 21:58:08 +05:30
real ( pReal ) , dimension ( 9 , 9 ) , intent ( in ) :: A , B
integer ( pInt ) i , j
2009-08-11 22:01:57 +05:30
2011-12-01 17:31:13 +05:30
forall ( i = 1_pInt : 9_pInt , j = 1_pInt : 9_pInt ) math_mul99x99 ( i , j ) = &
2009-08-11 22:01:57 +05:30
A ( i , 1 ) * B ( 1 , j ) + A ( i , 2 ) * B ( 2 , j ) + A ( i , 3 ) * B ( 3 , j ) + &
A ( i , 4 ) * B ( 4 , j ) + A ( i , 5 ) * B ( 5 , j ) + A ( i , 6 ) * B ( 6 , j ) + &
A ( i , 7 ) * B ( 7 , j ) + A ( i , 8 ) * B ( 8 , j ) + A ( i , 9 ) * B ( 9 , j )
2012-03-09 01:55:28 +05:30
end function math_mul99x99
2009-01-20 00:40:58 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief matrix multiplication 33x3 = 3
!--------------------------------------------------------------------------------------------------
pure function math_mul33x3 ( A , B )
2009-08-11 22:01:57 +05:30
implicit none
2013-01-31 21:58:08 +05:30
real ( pReal ) , dimension ( 3 ) :: math_mul33x3
2009-08-11 22:01:57 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: A
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: B
2013-01-31 21:58:08 +05:30
integer ( pInt ) :: i
2009-08-11 22:01:57 +05:30
2012-01-25 20:01:21 +05:30
forall ( i = 1_pInt : 3_pInt ) math_mul33x3 ( i ) = sum ( A ( i , 1 : 3 ) * B )
2009-08-11 22:01:57 +05:30
2012-03-09 01:55:28 +05:30
end function math_mul33x3
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief matrix multiplication complex(33) x real(3) = complex(3)
!--------------------------------------------------------------------------------------------------
pure function math_mul33x3_complex ( A , B )
2010-09-22 17:34:43 +05:30
implicit none
2013-01-31 21:58:08 +05:30
complex ( pReal ) , dimension ( 3 ) :: math_mul33x3_complex
2010-09-22 17:34:43 +05:30
complex ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: A
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: B
2013-01-31 21:58:08 +05:30
integer ( pInt ) :: i
2010-09-22 17:34:43 +05:30
2012-02-10 16:54:53 +05:30
forall ( i = 1_pInt : 3_pInt ) math_mul33x3_complex ( i ) = sum ( A ( i , 1 : 3 ) * cmplx ( B , 0.0_pReal , pReal ) )
2010-09-22 17:34:43 +05:30
2012-03-09 01:55:28 +05:30
end function math_mul33x3_complex
2009-08-11 22:01:57 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief matrix multiplication 66x6 = 6
!--------------------------------------------------------------------------------------------------
pure function math_mul66x6 ( A , B )
2009-01-20 00:40:58 +05:30
implicit none
2013-01-31 21:58:08 +05:30
real ( pReal ) , dimension ( 6 ) :: math_mul66x6
2009-01-20 00:40:58 +05:30
real ( pReal ) , dimension ( 6 , 6 ) , intent ( in ) :: A
real ( pReal ) , dimension ( 6 ) , intent ( in ) :: B
2013-01-31 21:58:08 +05:30
integer ( pInt ) :: i
2009-01-20 00:40:58 +05:30
2011-12-01 17:31:13 +05:30
forall ( i = 1_pInt : 6_pInt ) math_mul66x6 ( i ) = &
2009-01-20 00:40:58 +05:30
A ( i , 1 ) * B ( 1 ) + A ( i , 2 ) * B ( 2 ) + A ( i , 3 ) * B ( 3 ) + &
A ( i , 4 ) * B ( 4 ) + A ( i , 5 ) * B ( 5 ) + A ( i , 6 ) * B ( 6 )
2012-03-09 01:55:28 +05:30
end function math_mul66x6
2010-05-06 19:37:21 +05:30
2012-08-25 17:16:36 +05:30
2012-10-12 23:24:20 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief 3x3 matrix exponential up to series approximation order n (default 5)
!--------------------------------------------------------------------------------------------------
pure function math_exp33 ( A , n )
implicit none
integer ( pInt ) :: i , order
integer ( pInt ) , intent ( in ) , optional :: n
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: A
real ( pReal ) , dimension ( 3 , 3 ) :: B , math_exp33
real ( pReal ) :: invfac
2016-01-10 19:04:26 +05:30
order = merge ( n , 5_pInt , present ( n ) )
2012-10-12 23:24:20 +05:30
2016-01-10 19:04:26 +05:30
B = math_I3 ! init
2013-05-08 21:22:29 +05:30
invfac = 1.0_pReal ! 0!
math_exp33 = B ! A^0 = eye2
2012-10-12 23:24:20 +05:30
do i = 1_pInt , n
2016-02-26 20:06:24 +05:30
invfac = invfac / real ( i , pReal ) ! invfac = 1/i!
2012-10-12 23:24:20 +05:30
B = math_mul33x33 ( B , A )
2013-05-08 21:22:29 +05:30
math_exp33 = math_exp33 + invfac * B ! exp = SUM (A^i)/i!
2012-10-12 23:24:20 +05:30
enddo
end function math_exp33
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief transposition of a 33 matrix
!--------------------------------------------------------------------------------------------------
2012-01-26 19:20:00 +05:30
pure function math_transpose33 ( A )
2008-07-09 01:08:22 +05:30
2013-01-31 21:58:08 +05:30
implicit none
2012-01-26 19:20:00 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: math_transpose33
2013-01-31 21:58:08 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: A
2011-12-01 17:31:13 +05:30
integer ( pInt ) :: i , j
2012-08-25 17:16:36 +05:30
2012-01-26 19:20:00 +05:30
forall ( i = 1_pInt : 3_pInt , j = 1_pInt : 3_pInt ) math_transpose33 ( i , j ) = A ( j , i )
2008-07-09 01:08:22 +05:30
2012-03-09 01:55:28 +05:30
end function math_transpose33
2007-03-29 21:02:52 +05:30
2009-03-31 13:01:38 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Cramer inversion of 33 matrix (function)
2009-03-31 13:01:38 +05:30
! direct Cramer inversion of matrix A.
! returns all zeroes if not possible, i.e. if det close to zero
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
pure function math_inv33 ( A )
2016-05-27 15:16:34 +05:30
use prec , only : &
2016-06-30 14:41:35 +05:30
dNeq0
2009-03-31 13:01:38 +05:30
implicit none
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: A
2011-12-01 17:31:13 +05:30
real ( pReal ) :: DetA
2012-01-26 19:20:00 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: math_inv33
2012-08-25 17:16:36 +05:30
2015-05-05 12:07:59 +05:30
math_inv33 ( 1 , 1 ) = A ( 2 , 2 ) * A ( 3 , 3 ) - A ( 2 , 3 ) * A ( 3 , 2 )
math_inv33 ( 2 , 1 ) = - A ( 2 , 1 ) * A ( 3 , 3 ) + A ( 2 , 3 ) * A ( 3 , 1 )
math_inv33 ( 3 , 1 ) = A ( 2 , 1 ) * A ( 3 , 2 ) - A ( 2 , 2 ) * A ( 3 , 1 )
2009-03-31 13:01:38 +05:30
2015-05-05 12:07:59 +05:30
DetA = A ( 1 , 1 ) * math_inv33 ( 1 , 1 ) + A ( 1 , 2 ) * math_inv33 ( 2 , 1 ) + A ( 1 , 3 ) * math_inv33 ( 3 , 1 )
2009-03-31 13:01:38 +05:30
2016-06-30 14:41:35 +05:30
if ( dNeq0 ( DetA ) ) then
2015-05-05 12:07:59 +05:30
math_inv33 ( 1 , 2 ) = - A ( 1 , 2 ) * A ( 3 , 3 ) + A ( 1 , 3 ) * A ( 3 , 2 )
math_inv33 ( 2 , 2 ) = A ( 1 , 1 ) * A ( 3 , 3 ) - A ( 1 , 3 ) * A ( 3 , 1 )
math_inv33 ( 3 , 2 ) = - A ( 1 , 1 ) * A ( 3 , 2 ) + A ( 1 , 2 ) * A ( 3 , 1 )
2009-03-31 13:01:38 +05:30
2015-05-05 12:07:59 +05:30
math_inv33 ( 1 , 3 ) = A ( 1 , 2 ) * A ( 2 , 3 ) - A ( 1 , 3 ) * A ( 2 , 2 )
math_inv33 ( 2 , 3 ) = - A ( 1 , 1 ) * A ( 2 , 3 ) + A ( 1 , 3 ) * A ( 2 , 1 )
math_inv33 ( 3 , 3 ) = A ( 1 , 1 ) * A ( 2 , 2 ) - A ( 1 , 2 ) * A ( 2 , 1 )
math_inv33 = math_inv33 / DetA
else
2016-01-12 16:30:23 +05:30
math_inv33 = 0.0_pReal
2009-03-31 13:01:38 +05:30
endif
2012-03-09 01:55:28 +05:30
end function math_inv33
2009-03-31 13:01:38 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Cramer inversion of 33 matrix (subroutine)
2015-05-05 12:07:59 +05:30
! direct Cramer inversion of matrix A.
! also returns determinant
! returns error if not possible, i.e. if det close to zero
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
pure subroutine math_invert33 ( A , InvA , DetA , error )
2016-05-29 14:15:03 +05:30
use prec , only : &
2016-06-30 14:41:35 +05:30
dEq0
2008-02-15 18:12:27 +05:30
implicit none
logical , intent ( out ) :: error
2007-04-11 15:34:22 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: A
real ( pReal ) , dimension ( 3 , 3 ) , intent ( out ) :: InvA
real ( pReal ) , intent ( out ) :: DetA
2008-02-15 18:12:27 +05:30
2015-05-05 12:07:59 +05:30
InvA ( 1 , 1 ) = A ( 2 , 2 ) * A ( 3 , 3 ) - A ( 2 , 3 ) * A ( 3 , 2 )
InvA ( 2 , 1 ) = - A ( 2 , 1 ) * A ( 3 , 3 ) + A ( 2 , 3 ) * A ( 3 , 1 )
InvA ( 3 , 1 ) = A ( 2 , 1 ) * A ( 3 , 2 ) - A ( 2 , 2 ) * A ( 3 , 1 )
DetA = A ( 1 , 1 ) * InvA ( 1 , 1 ) + A ( 1 , 2 ) * InvA ( 2 , 1 ) + A ( 1 , 3 ) * InvA ( 3 , 1 )
2008-02-15 18:12:27 +05:30
2016-06-30 14:41:35 +05:30
if ( dEq0 ( DetA ) ) then
2016-02-26 20:06:24 +05:30
InvA = 0.0_pReal
2008-02-15 18:12:27 +05:30
error = . true .
2007-04-11 15:34:22 +05:30
else
2015-05-05 12:07:59 +05:30
InvA ( 1 , 2 ) = - A ( 1 , 2 ) * A ( 3 , 3 ) + A ( 1 , 3 ) * A ( 3 , 2 )
InvA ( 2 , 2 ) = A ( 1 , 1 ) * A ( 3 , 3 ) - A ( 1 , 3 ) * A ( 3 , 1 )
InvA ( 3 , 2 ) = - A ( 1 , 1 ) * A ( 3 , 2 ) + A ( 1 , 2 ) * A ( 3 , 1 )
2008-02-15 18:12:27 +05:30
2015-05-05 12:07:59 +05:30
InvA ( 1 , 3 ) = A ( 1 , 2 ) * A ( 2 , 3 ) - A ( 1 , 3 ) * A ( 2 , 2 )
InvA ( 2 , 3 ) = - A ( 1 , 1 ) * A ( 2 , 3 ) + A ( 1 , 3 ) * A ( 2 , 1 )
InvA ( 3 , 3 ) = A ( 1 , 1 ) * A ( 2 , 2 ) - A ( 1 , 2 ) * A ( 2 , 1 )
2012-08-25 17:16:36 +05:30
2015-05-05 12:07:59 +05:30
InvA = InvA / DetA
2008-02-15 18:12:27 +05:30
error = . false .
2007-04-11 15:34:22 +05:30
endif
2012-03-09 01:55:28 +05:30
end subroutine math_invert33
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Inversion of symmetriced 3x3x3x3 tensor.
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
function math_invSym3333 ( A )
2014-02-28 18:58:27 +05:30
use IO , only : &
IO_error
2012-08-25 17:16:36 +05:30
2012-03-09 01:55:28 +05:30
implicit none
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: math_invSym3333
2012-08-25 17:16:36 +05:30
2012-03-09 01:55:28 +05:30
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) , intent ( in ) :: A
2012-08-31 01:56:28 +05:30
integer ( pInt ) :: ierr
2012-03-09 01:55:28 +05:30
integer ( pInt ) , dimension ( 6 ) :: ipiv6
real ( pReal ) , dimension ( 6 , 6 ) :: temp66_Real
real ( pReal ) , dimension ( 6 ) :: work6
2016-03-21 03:50:58 +05:30
external :: &
dgetrf , &
dgetri
2012-08-25 17:16:36 +05:30
2012-03-09 01:55:28 +05:30
temp66_real = math_Mandel3333to66 ( A )
2012-08-31 01:56:28 +05:30
call dgetrf ( 6 , 6 , temp66_real , 6 , ipiv6 , ierr )
call dgetri ( 6 , temp66_real , 6 , ipiv6 , work6 , 6 , ierr )
if ( ierr == 0_pInt ) then
2012-03-09 01:55:28 +05:30
math_invSym3333 = math_Mandel66to3333 ( temp66_real )
2012-08-25 17:16:36 +05:30
else
2012-03-09 01:55:28 +05:30
call IO_error ( 400_pInt , ext_msg = 'math_invSym3333' )
endif
2008-02-15 18:12:27 +05:30
2012-03-09 01:55:28 +05:30
end function math_invSym3333
2007-03-29 21:02:52 +05:30
2008-02-15 18:12:27 +05:30
2012-08-27 13:34:47 +05:30
!--------------------------------------------------------------------------------------------------
2012-08-28 22:29:45 +05:30
!> @brief invert matrix of arbitrary dimension
2012-08-27 13:34:47 +05:30
!--------------------------------------------------------------------------------------------------
2012-08-28 22:29:45 +05:30
subroutine math_invert ( myDim , A , InvA , error )
2008-02-15 18:12:27 +05:30
implicit none
2012-08-28 22:29:45 +05:30
integer ( pInt ) , intent ( in ) :: myDim
real ( pReal ) , dimension ( myDim , myDim ) , intent ( in ) :: A
2012-08-24 18:57:55 +05:30
2012-08-28 22:29:45 +05:30
integer ( pInt ) :: ierr
integer ( pInt ) , dimension ( myDim ) :: ipiv
real ( pReal ) , dimension ( myDim ) :: work
real ( pReal ) , dimension ( myDim , myDim ) , intent ( out ) :: invA
2011-12-01 17:31:13 +05:30
logical , intent ( out ) :: error
2016-04-26 13:06:27 +05:30
external :: &
dgetrf , &
dgetri
2012-08-28 22:29:45 +05:30
invA = A
call dgetrf ( myDim , myDim , invA , myDim , ipiv , ierr )
call dgetri ( myDim , InvA , myDim , ipiv , work , myDim , ierr )
2015-05-05 12:07:59 +05:30
error = merge ( . true . , . false . , ierr / = 0_pInt ) ! http://fortraninacworld.blogspot.de/2012/12/ternary-operator.html
2012-08-28 22:29:45 +05:30
end subroutine math_invert
2008-02-15 18:12:27 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief symmetrize a 33 matrix
!--------------------------------------------------------------------------------------------------
2016-01-10 19:04:26 +05:30
pure function math_symmetric33 ( m )
2008-02-15 18:12:27 +05:30
implicit none
2012-01-26 19:20:00 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: math_symmetric33
2011-12-01 17:31:13 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m
2008-02-15 18:12:27 +05:30
2016-01-10 19:04:26 +05:30
math_symmetric33 = 0.5_pReal * ( m + transpose ( m ) )
2008-02-15 18:12:27 +05:30
2012-03-09 01:55:28 +05:30
end function math_symmetric33
2008-02-15 18:12:27 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief symmetrize a 66 matrix
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
pure function math_symmetric66 ( m )
2008-02-15 18:12:27 +05:30
implicit none
2012-01-26 19:20:00 +05:30
real ( pReal ) , dimension ( 6 , 6 ) :: math_symmetric66
2016-01-10 19:04:26 +05:30
real ( pReal ) , dimension ( 6 , 6 ) , intent ( in ) :: m
2012-08-25 17:16:36 +05:30
2016-01-10 19:04:26 +05:30
math_symmetric66 = 0.5_pReal * ( m + transpose ( m ) )
2012-01-26 19:20:00 +05:30
2012-03-09 01:55:28 +05:30
end function math_symmetric66
2008-02-15 18:12:27 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief skew part of a 33 matrix
!--------------------------------------------------------------------------------------------------
2012-02-09 21:28:15 +05:30
pure function math_skew33 ( m )
2012-01-25 16:00:39 +05:30
implicit none
2012-01-26 19:20:00 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: math_skew33
2012-01-25 16:00:39 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m
2012-08-25 17:16:36 +05:30
2016-01-10 19:04:26 +05:30
math_skew33 = m - math_symmetric33 ( m )
2012-01-26 19:20:00 +05:30
2012-03-09 01:55:28 +05:30
end function math_skew33
2008-02-15 18:12:27 +05:30
2016-01-09 00:27:37 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief hydrostatic part of a 33 matrix
!--------------------------------------------------------------------------------------------------
pure function math_spherical33 ( m )
implicit none
real ( pReal ) , dimension ( 3 , 3 ) :: math_spherical33
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m
2016-01-10 19:04:26 +05:30
math_spherical33 = math_I3 * math_trace33 ( m ) / 3.0_pReal
2016-01-09 00:27:37 +05:30
end function math_spherical33
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief deviatoric part of a 33 matrix
!--------------------------------------------------------------------------------------------------
2012-02-09 21:28:15 +05:30
pure function math_deviatoric33 ( m )
2012-01-26 19:20:00 +05:30
2012-02-09 21:28:15 +05:30
implicit none
real ( pReal ) , dimension ( 3 , 3 ) :: math_deviatoric33
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m
2012-01-26 19:20:00 +05:30
2016-01-10 19:04:26 +05:30
math_deviatoric33 = m - math_spherical33 ( m )
2012-01-26 19:20:00 +05:30
2012-03-09 01:55:28 +05:30
end function math_deviatoric33
2012-01-26 19:20:00 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2015-05-05 12:07:59 +05:30
!> @brief equivalent scalar quantity of a full symmetric strain tensor
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
pure function math_equivStrain33 ( m )
2010-03-24 18:50:12 +05:30
implicit none
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m
2015-05-05 12:07:59 +05:30
real ( pReal ) , dimension ( 3 ) :: e , s
real ( pReal ) :: math_equivStrain33
real ( pReal ) , parameter :: TWOTHIRD = 2.0_pReal / 3.0_pReal
2010-03-24 18:50:12 +05:30
2015-05-05 12:07:59 +05:30
e = [ 2.0_pReal * m ( 1 , 1 ) - m ( 2 , 2 ) - m ( 3 , 3 ) , &
2.0_pReal * m ( 2 , 2 ) - m ( 3 , 3 ) - m ( 1 , 1 ) , &
2.0_pReal * m ( 3 , 3 ) - m ( 1 , 1 ) - m ( 2 , 2 ) ] / 3.0_pReal
s = [ m ( 1 , 2 ) , m ( 2 , 3 ) , m ( 1 , 3 ) ] * 2.0_pReal
2010-03-24 18:50:12 +05:30
2015-05-05 12:07:59 +05:30
math_equivStrain33 = TWOTHIRD * ( 1.50_pReal * ( sum ( e ** 2.0_pReal ) ) + &
0.75_pReal * ( sum ( s ** 2.0_pReal ) ) ) ** ( 0.5_pReal )
2010-03-24 18:50:12 +05:30
2012-03-09 01:55:28 +05:30
end function math_equivStrain33
2010-03-24 18:50:12 +05:30
2015-05-05 12:07:59 +05:30
2015-01-13 15:13:05 +05:30
!--------------------------------------------------------------------------------------------------
2015-05-05 12:07:59 +05:30
!> @brief von Mises equivalent of a full symmetric stress tensor
2015-01-13 15:13:05 +05:30
!--------------------------------------------------------------------------------------------------
pure function math_equivStress33 ( m )
implicit none
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m
real ( pReal ) :: math_equivStress33
2015-05-05 12:07:59 +05:30
math_equivStress33 = ( ( ( m ( 1 , 1 ) - m ( 2 , 2 ) ) ** 2.0_pReal + &
( m ( 2 , 2 ) - m ( 3 , 3 ) ) ** 2.0_pReal + &
( m ( 3 , 3 ) - m ( 1 , 1 ) ) ** 2.0_pReal + &
6.0_pReal * ( m ( 1 , 2 ) ** 2.0_pReal + &
m ( 2 , 3 ) ** 2.0_pReal + &
m ( 1 , 3 ) ** 2.0_pReal &
) &
) ** 0.5_pReal &
) / sqrt ( 2.0_pReal )
2015-01-13 15:13:05 +05:30
end function math_equivStress33
2015-05-05 12:07:59 +05:30
2012-10-12 23:24:20 +05:30
!--------------------------------------------------------------------------------------------------
2014-10-14 09:00:59 +05:30
!> @brief trace of a 33 matrix
2012-10-12 23:24:20 +05:30
!--------------------------------------------------------------------------------------------------
2013-01-31 21:58:08 +05:30
real ( pReal ) pure function math_trace33 ( m )
2012-10-12 23:24:20 +05:30
implicit none
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m
math_trace33 = m ( 1 , 1 ) + m ( 2 , 2 ) + m ( 3 , 3 )
end function math_trace33
2015-05-05 12:07:59 +05:30
2012-08-27 13:34:47 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief determinant of a 33 matrix
!--------------------------------------------------------------------------------------------------
2013-01-31 21:58:08 +05:30
real ( pReal ) pure function math_det33 ( m )
2007-03-20 19:25:22 +05:30
implicit none
2009-01-16 20:57:13 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m
2013-01-31 21:58:08 +05:30
2016-01-10 19:04:26 +05:30
math_det33 = m ( 1 , 1 ) * ( m ( 2 , 2 ) * m ( 3 , 3 ) - m ( 2 , 3 ) * m ( 3 , 2 ) ) &
- m ( 1 , 2 ) * ( m ( 2 , 1 ) * m ( 3 , 3 ) - m ( 2 , 3 ) * m ( 3 , 1 ) ) &
+ m ( 1 , 3 ) * ( m ( 2 , 1 ) * m ( 3 , 2 ) - m ( 2 , 2 ) * m ( 3 , 1 ) )
2007-03-20 19:25:22 +05:30
2012-03-09 01:55:28 +05:30
end function math_det33
2007-03-21 15:50:25 +05:30
2012-08-25 17:16:36 +05:30
2016-02-02 17:53:45 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief determinant of a symmetric 33 matrix
!--------------------------------------------------------------------------------------------------
real ( pReal ) pure function math_detSym33 ( m )
implicit none
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m
math_detSym33 = - ( m ( 1 , 1 ) * m ( 2 , 3 ) ** 2_pInt + m ( 2 , 2 ) * m ( 1 , 3 ) ** 2_pInt + m ( 3 , 3 ) * m ( 1 , 2 ) ** 2_pInt ) &
2016-02-02 23:29:04 +05:30
+ m ( 1 , 1 ) * m ( 2 , 2 ) * m ( 3 , 3 ) - 2.0_pReal * m ( 1 , 2 ) * m ( 1 , 3 ) * m ( 2 , 3 )
2016-02-02 17:53:45 +05:30
end function math_detSym33
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief convert 33 matrix into vector 9
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
pure function math_Plain33to9 ( m33 )
2008-02-15 18:12:27 +05:30
implicit none
real ( pReal ) , dimension ( 9 ) :: math_Plain33to9
2013-01-31 21:58:08 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m33
2011-12-01 17:31:13 +05:30
integer ( pInt ) :: i
2012-08-25 17:16:36 +05:30
2011-12-01 17:31:13 +05:30
forall ( i = 1_pInt : 9_pInt ) math_Plain33to9 ( i ) = m33 ( mapPlain ( 1 , i ) , mapPlain ( 2 , i ) )
2008-02-15 18:12:27 +05:30
2012-03-09 01:55:28 +05:30
end function math_Plain33to9
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief convert Plain 9 back to 33 matrix
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
pure function math_Plain9to33 ( v9 )
2008-02-15 18:12:27 +05:30
implicit none
real ( pReal ) , dimension ( 3 , 3 ) :: math_Plain9to33
2013-01-31 21:58:08 +05:30
real ( pReal ) , dimension ( 9 ) , intent ( in ) :: v9
2011-12-01 17:31:13 +05:30
integer ( pInt ) :: i
2012-08-25 17:16:36 +05:30
2011-12-01 17:31:13 +05:30
forall ( i = 1_pInt : 9_pInt ) math_Plain9to33 ( mapPlain ( 1 , i ) , mapPlain ( 2 , i ) ) = v9 ( i )
2008-02-15 18:12:27 +05:30
2012-03-09 01:55:28 +05:30
end function math_Plain9to33
2008-02-15 18:12:27 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief convert symmetric 33 matrix into Mandel vector 6
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
pure function math_Mandel33to6 ( m33 )
2007-03-28 12:51:47 +05:30
implicit none
real ( pReal ) , dimension ( 6 ) :: math_Mandel33to6
2013-01-31 21:58:08 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m33
2011-12-01 17:31:13 +05:30
integer ( pInt ) :: i
2012-08-25 17:16:36 +05:30
2011-12-01 17:31:13 +05:30
forall ( i = 1_pInt : 6_pInt ) math_Mandel33to6 ( i ) = nrmMandel ( i ) * m33 ( mapMandel ( 1 , i ) , mapMandel ( 2 , i ) )
2007-03-28 12:51:47 +05:30
2012-03-09 01:55:28 +05:30
end function math_Mandel33to6
2007-03-28 12:51:47 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief convert Mandel 6 back to symmetric 33 matrix
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
pure function math_Mandel6to33 ( v6 )
2007-03-28 12:51:47 +05:30
implicit none
2009-01-16 20:57:13 +05:30
real ( pReal ) , dimension ( 6 ) , intent ( in ) :: v6
2007-03-28 12:51:47 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: math_Mandel6to33
2011-12-01 17:31:13 +05:30
integer ( pInt ) :: i
2012-08-25 17:16:36 +05:30
2011-12-01 17:31:13 +05:30
forall ( i = 1_pInt : 6_pInt )
2007-03-28 13:50:50 +05:30
math_Mandel6to33 ( mapMandel ( 1 , i ) , mapMandel ( 2 , i ) ) = invnrmMandel ( i ) * v6 ( i )
math_Mandel6to33 ( mapMandel ( 2 , i ) , mapMandel ( 1 , i ) ) = invnrmMandel ( i ) * v6 ( i )
2007-03-28 12:51:47 +05:30
end forall
2012-03-09 01:55:28 +05:30
end function math_Mandel6to33
2007-03-28 12:51:47 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief convert 3333 tensor into plain matrix 99
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
pure function math_Plain3333to99 ( m3333 )
2008-02-15 18:12:27 +05:30
implicit none
2009-01-16 20:57:13 +05:30
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) , intent ( in ) :: m3333
2008-02-15 18:12:27 +05:30
real ( pReal ) , dimension ( 9 , 9 ) :: math_Plain3333to99
2011-12-01 17:31:13 +05:30
integer ( pInt ) :: i , j
2012-08-25 17:16:36 +05:30
2011-12-01 17:31:13 +05:30
forall ( i = 1_pInt : 9_pInt , j = 1_pInt : 9_pInt ) math_Plain3333to99 ( i , j ) = &
2008-02-15 18:12:27 +05:30
m3333 ( mapPlain ( 1 , i ) , mapPlain ( 2 , i ) , mapPlain ( 1 , j ) , mapPlain ( 2 , j ) )
2012-03-09 01:55:28 +05:30
end function math_Plain3333to99
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief plain matrix 99 into 3333 tensor
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
pure function math_Plain99to3333 ( m99 )
2010-09-22 17:34:43 +05:30
implicit none
real ( pReal ) , dimension ( 9 , 9 ) , intent ( in ) :: m99
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: math_Plain99to3333
2011-12-01 17:31:13 +05:30
integer ( pInt ) :: i , j
2012-08-25 17:16:36 +05:30
2011-12-01 17:31:13 +05:30
forall ( i = 1_pInt : 9_pInt , j = 1_pInt : 9_pInt ) math_Plain99to3333 ( mapPlain ( 1 , i ) , mapPlain ( 2 , i ) , &
2016-06-30 14:41:35 +05:30
mapPlain ( 1 , j ) , mapPlain ( 2 , j ) ) = m99 ( i , j )
2008-02-15 18:12:27 +05:30
2012-03-09 01:55:28 +05:30
end function math_Plain99to3333
2008-02-15 18:12:27 +05:30
2011-07-29 21:27:39 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief convert Mandel matrix 66 into Plain matrix 66
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
pure function math_Mandel66toPlain66 ( m66 )
2011-07-29 21:27:39 +05:30
implicit none
real ( pReal ) , dimension ( 6 , 6 ) , intent ( in ) :: m66
real ( pReal ) , dimension ( 6 , 6 ) :: math_Mandel66toPlain66
2011-12-01 17:31:13 +05:30
integer ( pInt ) :: i , j
2012-08-25 17:16:36 +05:30
2011-12-01 17:31:13 +05:30
forall ( i = 1_pInt : 6_pInt , j = 1_pInt : 6_pInt ) &
2011-07-29 21:27:39 +05:30
math_Mandel66toPlain66 ( i , j ) = invnrmMandel ( i ) * invnrmMandel ( j ) * m66 ( i , j )
2013-01-31 21:58:08 +05:30
end function math_Mandel66toPlain66
2011-07-29 21:27:39 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief convert Plain matrix 66 into Mandel matrix 66
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
pure function math_Plain66toMandel66 ( m66 )
2011-07-29 21:27:39 +05:30
implicit none
real ( pReal ) , dimension ( 6 , 6 ) , intent ( in ) :: m66
real ( pReal ) , dimension ( 6 , 6 ) :: math_Plain66toMandel66
2015-05-05 12:07:59 +05:30
integer ( pInt ) :: i , j
2012-08-25 17:16:36 +05:30
2011-12-01 17:31:13 +05:30
forall ( i = 1_pInt : 6_pInt , j = 1_pInt : 6_pInt ) &
2011-07-29 21:27:39 +05:30
math_Plain66toMandel66 ( i , j ) = nrmMandel ( i ) * nrmMandel ( j ) * m66 ( i , j )
2013-01-31 21:58:08 +05:30
end function math_Plain66toMandel66
2011-07-29 21:27:39 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief convert symmetric 3333 tensor into Mandel matrix 66
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
pure function math_Mandel3333to66 ( m3333 )
2007-03-28 12:51:47 +05:30
implicit none
2009-01-16 20:57:13 +05:30
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) , intent ( in ) :: m3333
2007-03-28 12:51:47 +05:30
real ( pReal ) , dimension ( 6 , 6 ) :: math_Mandel3333to66
2011-12-01 17:31:13 +05:30
integer ( pInt ) :: i , j
2012-08-25 17:16:36 +05:30
2011-12-01 17:31:13 +05:30
forall ( i = 1_pInt : 6_pInt , j = 1_pInt : 6_pInt ) math_Mandel3333to66 ( i , j ) = &
2007-03-28 13:50:50 +05:30
nrmMandel ( i ) * nrmMandel ( j ) * m3333 ( mapMandel ( 1 , i ) , mapMandel ( 2 , i ) , mapMandel ( 1 , j ) , mapMandel ( 2 , j ) )
2007-03-28 12:51:47 +05:30
2012-03-09 01:55:28 +05:30
end function math_Mandel3333to66
2010-05-06 19:37:21 +05:30
2011-12-01 17:31:13 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief convert Mandel matrix 66 back to symmetric 3333 tensor
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
pure function math_Mandel66to3333 ( m66 )
2008-02-15 18:12:27 +05:30
implicit none
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: math_Mandel66to3333
2013-01-31 21:58:08 +05:30
real ( pReal ) , dimension ( 6 , 6 ) , intent ( in ) :: m66
2011-12-01 17:31:13 +05:30
integer ( pInt ) :: i , j
2012-08-25 17:16:36 +05:30
forall ( i = 1_pInt : 6_pInt , j = 1_pInt : 6_pInt )
2013-01-31 21:58:08 +05:30
math_Mandel66to3333 ( mapMandel ( 1 , i ) , mapMandel ( 2 , i ) , mapMandel ( 1 , j ) , mapMandel ( 2 , j ) ) = &
invnrmMandel ( i ) * invnrmMandel ( j ) * m66 ( i , j )
math_Mandel66to3333 ( mapMandel ( 2 , i ) , mapMandel ( 1 , i ) , mapMandel ( 1 , j ) , mapMandel ( 2 , j ) ) = &
invnrmMandel ( i ) * invnrmMandel ( j ) * m66 ( i , j )
math_Mandel66to3333 ( mapMandel ( 1 , i ) , mapMandel ( 2 , i ) , mapMandel ( 2 , j ) , mapMandel ( 1 , j ) ) = &
invnrmMandel ( i ) * invnrmMandel ( j ) * m66 ( i , j )
math_Mandel66to3333 ( mapMandel ( 2 , i ) , mapMandel ( 1 , i ) , mapMandel ( 2 , j ) , mapMandel ( 1 , j ) ) = &
invnrmMandel ( i ) * invnrmMandel ( j ) * m66 ( i , j )
2008-02-15 18:12:27 +05:30
end forall
2012-03-09 01:55:28 +05:30
end function math_Mandel66to3333
2008-02-15 18:12:27 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief convert Voigt matrix 66 back to symmetric 3333 tensor
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
pure function math_Voigt66to3333 ( m66 )
2008-02-15 18:12:27 +05:30
implicit none
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: math_Voigt66to3333
2013-01-31 21:58:08 +05:30
real ( pReal ) , dimension ( 6 , 6 ) , intent ( in ) :: m66
2011-12-01 17:31:13 +05:30
integer ( pInt ) :: i , j
2012-08-25 17:16:36 +05:30
forall ( i = 1_pInt : 6_pInt , j = 1_pInt : 6_pInt )
2013-05-08 21:40:21 +05:30
math_Voigt66to3333 ( mapVoigt ( 1 , i ) , mapVoigt ( 2 , i ) , mapVoigt ( 1 , j ) , mapVoigt ( 2 , j ) ) = &
invnrmVoigt ( i ) * invnrmVoigt ( j ) * m66 ( i , j )
math_Voigt66to3333 ( mapVoigt ( 2 , i ) , mapVoigt ( 1 , i ) , mapVoigt ( 1 , j ) , mapVoigt ( 2 , j ) ) = &
invnrmVoigt ( i ) * invnrmVoigt ( j ) * m66 ( i , j )
math_Voigt66to3333 ( mapVoigt ( 1 , i ) , mapVoigt ( 2 , i ) , mapVoigt ( 2 , j ) , mapVoigt ( 1 , j ) ) = &
invnrmVoigt ( i ) * invnrmVoigt ( j ) * m66 ( i , j )
math_Voigt66to3333 ( mapVoigt ( 2 , i ) , mapVoigt ( 1 , i ) , mapVoigt ( 2 , j ) , mapVoigt ( 1 , j ) ) = &
invnrmVoigt ( i ) * invnrmVoigt ( j ) * m66 ( i , j )
2008-02-15 18:12:27 +05:30
end forall
2012-03-09 01:55:28 +05:30
end function math_Voigt66to3333
2008-02-15 18:12:27 +05:30
2013-01-31 21:58:08 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief random quaternion
2017-04-03 02:47:22 +05:30
! http://math.stackexchange.com/questions/131336/uniform-random-quaternion-in-a-restricted-angle-range
! K. Shoemake. Uniform random rotations. In D. Kirk, editor, Graphics Gems III, pages 124-132.
! Academic, New York, 1992.
2013-01-31 21:58:08 +05:30
!--------------------------------------------------------------------------------------------------
function math_qRand ( )
implicit none
real ( pReal ) , dimension ( 4 ) :: math_qRand
real ( pReal ) , dimension ( 3 ) :: rnd
call halton ( 3_pInt , rnd )
2016-07-14 18:59:53 +05:30
math_qRand = [ cos ( 2.0_pReal * PI * rnd ( 1 ) ) * sqrt ( rnd ( 3 ) ) , &
sin ( 2.0_pReal * PI * rnd ( 2 ) ) * sqrt ( 1.0_pReal - rnd ( 3 ) ) , &
cos ( 2.0_pReal * PI * rnd ( 2 ) ) * sqrt ( 1.0_pReal - rnd ( 3 ) ) , &
2016-06-30 14:41:35 +05:30
sin ( 2.0_pReal * PI * rnd ( 1 ) ) * sqrt ( rnd ( 3 ) ) ]
2013-01-31 21:58:08 +05:30
end function math_qRand
!--------------------------------------------------------------------------------------------------
!> @brief quaternion multiplication q1xq2 = q12
!--------------------------------------------------------------------------------------------------
pure function math_qMul ( A , B )
implicit none
real ( pReal ) , dimension ( 4 ) :: math_qMul
real ( pReal ) , dimension ( 4 ) , intent ( in ) :: A , B
2016-01-10 19:04:26 +05:30
math_qMul = [ A ( 1 ) * B ( 1 ) - A ( 2 ) * B ( 2 ) - A ( 3 ) * B ( 3 ) - A ( 4 ) * B ( 4 ) , &
A ( 1 ) * B ( 2 ) + A ( 2 ) * B ( 1 ) + A ( 3 ) * B ( 4 ) - A ( 4 ) * B ( 3 ) , &
A ( 1 ) * B ( 3 ) - A ( 2 ) * B ( 4 ) + A ( 3 ) * B ( 1 ) + A ( 4 ) * B ( 2 ) , &
A ( 1 ) * B ( 4 ) + A ( 2 ) * B ( 3 ) - A ( 3 ) * B ( 2 ) + A ( 4 ) * B ( 1 ) ]
2013-01-31 21:58:08 +05:30
end function math_qMul
!--------------------------------------------------------------------------------------------------
!> @brief quaternion dotproduct
!--------------------------------------------------------------------------------------------------
real ( pReal ) pure function math_qDot ( A , B )
implicit none
real ( pReal ) , dimension ( 4 ) , intent ( in ) :: A , B
math_qDot = sum ( A * B )
end function math_qDot
!--------------------------------------------------------------------------------------------------
!> @brief quaternion conjugation
!--------------------------------------------------------------------------------------------------
pure function math_qConj ( Q )
implicit none
real ( pReal ) , dimension ( 4 ) :: math_qConj
real ( pReal ) , dimension ( 4 ) , intent ( in ) :: Q
2016-01-10 19:04:26 +05:30
math_qConj = [ Q ( 1 ) , - Q ( 2 : 4 ) ]
2013-01-31 21:58:08 +05:30
end function math_qConj
!--------------------------------------------------------------------------------------------------
!> @brief quaternion norm
!--------------------------------------------------------------------------------------------------
real ( pReal ) pure function math_qNorm ( Q )
implicit none
real ( pReal ) , dimension ( 4 ) , intent ( in ) :: Q
2016-01-10 19:04:26 +05:30
math_qNorm = norm2 ( Q )
2013-01-31 21:58:08 +05:30
end function math_qNorm
!--------------------------------------------------------------------------------------------------
!> @brief quaternion inversion
!--------------------------------------------------------------------------------------------------
pure function math_qInv ( Q )
2016-05-29 14:15:03 +05:30
use prec , only : &
2016-06-30 14:41:35 +05:30
dNeq0
2013-01-31 21:58:08 +05:30
implicit none
real ( pReal ) , dimension ( 4 ) , intent ( in ) :: Q
real ( pReal ) , dimension ( 4 ) :: math_qInv
real ( pReal ) :: squareNorm
math_qInv = 0.0_pReal
squareNorm = math_qDot ( Q , Q )
2016-06-30 14:41:35 +05:30
if ( dNeq0 ( squareNorm ) ) math_qInv = math_qConj ( Q ) / squareNorm
2013-01-31 21:58:08 +05:30
end function math_qInv
!--------------------------------------------------------------------------------------------------
!> @brief action of a quaternion on a vector (rotate vector v with Q)
!--------------------------------------------------------------------------------------------------
pure function math_qRot ( Q , v )
implicit none
real ( pReal ) , dimension ( 4 ) , intent ( in ) :: Q
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: v
real ( pReal ) , dimension ( 3 ) :: math_qRot
real ( pReal ) , dimension ( 4 , 4 ) :: T
integer ( pInt ) :: i , j
do i = 1_pInt , 4_pInt
do j = 1_pInt , i
T ( i , j ) = Q ( i ) * Q ( j )
enddo
enddo
2014-05-15 14:22:16 +05:30
math_qRot = [ - v ( 1 ) * ( T ( 3 , 3 ) + T ( 4 , 4 ) ) + v ( 2 ) * ( T ( 3 , 2 ) - T ( 4 , 1 ) ) + v ( 3 ) * ( T ( 4 , 2 ) + T ( 3 , 1 ) ) , &
v ( 1 ) * ( T ( 3 , 2 ) + T ( 4 , 1 ) ) - v ( 2 ) * ( T ( 2 , 2 ) + T ( 4 , 4 ) ) + v ( 3 ) * ( T ( 4 , 3 ) - T ( 2 , 1 ) ) , &
v ( 1 ) * ( T ( 4 , 2 ) - T ( 3 , 1 ) ) + v ( 2 ) * ( T ( 4 , 3 ) + T ( 2 , 1 ) ) - v ( 3 ) * ( T ( 2 , 2 ) + T ( 3 , 3 ) ) ]
2013-01-31 21:58:08 +05:30
math_qRot = 2.0_pReal * math_qRot + v
end function math_qRot
2013-06-06 00:40:37 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Euler angles (in radians) from rotation matrix
2013-06-06 00:40:37 +05:30
!> @details rotation matrix is meant to represent a PASSIVE rotation,
!> composed of INTRINSIC rotations around the axes of the
!> rotating reference frame
!> (see http://en.wikipedia.org/wiki/Euler_angles for definitions)
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
pure function math_RtoEuler ( R )
2007-03-20 19:25:22 +05:30
implicit none
2009-01-16 20:57:13 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: R
2007-03-29 21:02:52 +05:30
real ( pReal ) , dimension ( 3 ) :: math_RtoEuler
2016-01-10 19:04:26 +05:30
real ( pReal ) :: sqhkl , squvw , sqhk
2007-03-29 21:02:52 +05:30
sqhkl = sqrt ( R ( 1 , 3 ) * R ( 1 , 3 ) + R ( 2 , 3 ) * R ( 2 , 3 ) + R ( 3 , 3 ) * R ( 3 , 3 ) )
squvw = sqrt ( R ( 1 , 1 ) * R ( 1 , 1 ) + R ( 2 , 1 ) * R ( 2 , 1 ) + R ( 3 , 1 ) * R ( 3 , 1 ) )
2016-01-10 19:04:26 +05:30
sqhk = sqrt ( R ( 1 , 3 ) * R ( 1 , 3 ) + R ( 2 , 3 ) * R ( 2 , 3 ) )
2012-08-25 17:16:36 +05:30
2016-01-10 19:04:26 +05:30
! calculate PHI
math_RtoEuler ( 2 ) = acos ( math_limit ( R ( 3 , 3 ) / sqhkl , - 1.0_pReal , 1.0_pReal ) )
2007-03-21 15:50:25 +05:30
2014-01-15 14:02:43 +05:30
if ( ( math_RtoEuler ( 2 ) < 1.0e-8_pReal ) . or . ( pi - math_RtoEuler ( 2 ) < 1.0e-8_pReal ) ) then
2016-01-10 19:04:26 +05:30
math_RtoEuler ( 3 ) = 0.0_pReal
math_RtoEuler ( 1 ) = acos ( math_limit ( R ( 1 , 1 ) / squvw , - 1.0_pReal , 1.0_pReal ) )
if ( R ( 2 , 1 ) > 0.0_pReal ) math_RtoEuler ( 1 ) = 2.0_pReal * pi - math_RtoEuler ( 1 )
2007-03-29 21:02:52 +05:30
else
2016-01-10 19:04:26 +05:30
math_RtoEuler ( 3 ) = acos ( math_limit ( R ( 2 , 3 ) / sqhk , - 1.0_pReal , 1.0_pReal ) )
if ( R ( 1 , 3 ) < 0.0 ) math_RtoEuler ( 3 ) = 2.0_pReal * pi - math_RtoEuler ( 3 )
math_RtoEuler ( 1 ) = acos ( math_limit ( - R ( 3 , 2 ) / sin ( math_RtoEuler ( 2 ) ) , - 1.0_pReal , 1.0_pReal ) )
if ( R ( 3 , 1 ) < 0.0 ) math_RtoEuler ( 1 ) = 2.0_pReal * pi - math_RtoEuler ( 1 )
2007-03-29 21:02:52 +05:30
end if
2012-08-25 17:16:36 +05:30
2012-03-09 01:55:28 +05:30
end function math_RtoEuler
2010-05-06 19:37:21 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2013-06-06 00:40:37 +05:30
!> @brief converts a rotation matrix into a quaternion (w+ix+jy+kz)
!> @details math adopted from http://arxiv.org/pdf/math/0701759v1.pdf
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2013-01-31 21:58:08 +05:30
pure function math_RtoQ ( R )
2010-05-06 19:37:21 +05:30
implicit none
2013-06-06 00:40:37 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: R
2013-01-31 21:58:08 +05:30
real ( pReal ) , dimension ( 4 ) :: absQ , math_RtoQ
2011-12-01 17:31:13 +05:30
real ( pReal ) :: max_absQ
2012-02-14 19:13:36 +05:30
integer , dimension ( 1 ) :: largest !no pInt, maxloc returns integer default
2010-05-26 21:22:54 +05:30
2013-01-31 21:58:08 +05:30
math_RtoQ = 0.0_pReal
2011-03-03 16:17:07 +05:30
2016-01-10 19:04:26 +05:30
absQ = [ + R ( 1 , 1 ) + R ( 2 , 2 ) + R ( 3 , 3 ) , &
+ R ( 1 , 1 ) - R ( 2 , 2 ) - R ( 3 , 3 ) , &
- R ( 1 , 1 ) + R ( 2 , 2 ) - R ( 3 , 3 ) , &
- R ( 1 , 1 ) - R ( 2 , 2 ) + R ( 3 , 3 ) ] + 1.0_pReal
2011-03-03 16:17:07 +05:30
2013-06-06 00:40:37 +05:30
largest = maxloc ( absQ )
2011-03-03 16:17:07 +05:30
2016-01-10 19:04:26 +05:30
largestComponent : select case ( largest ( 1 ) )
case ( 1 ) largestComponent
2011-12-01 17:31:13 +05:30
!1----------------------------------
2013-06-06 00:40:37 +05:30
math_RtoQ ( 2 ) = R ( 3 , 2 ) - R ( 2 , 3 )
math_RtoQ ( 3 ) = R ( 1 , 3 ) - R ( 3 , 1 )
math_RtoQ ( 4 ) = R ( 2 , 1 ) - R ( 1 , 2 )
2016-01-10 19:04:26 +05:30
case ( 2 ) largestComponent
2013-06-06 00:40:37 +05:30
math_RtoQ ( 1 ) = R ( 3 , 2 ) - R ( 2 , 3 )
2011-12-01 17:31:13 +05:30
!2----------------------------------
2013-06-06 00:40:37 +05:30
math_RtoQ ( 3 ) = R ( 2 , 1 ) + R ( 1 , 2 )
math_RtoQ ( 4 ) = R ( 1 , 3 ) + R ( 3 , 1 )
2016-01-10 19:04:26 +05:30
case ( 3 ) largestComponent
2013-06-06 00:40:37 +05:30
math_RtoQ ( 1 ) = R ( 1 , 3 ) - R ( 3 , 1 )
math_RtoQ ( 2 ) = R ( 2 , 1 ) + R ( 1 , 2 )
2011-12-01 17:31:13 +05:30
!3----------------------------------
2013-06-06 00:40:37 +05:30
math_RtoQ ( 4 ) = R ( 3 , 2 ) + R ( 2 , 3 )
2016-01-10 19:04:26 +05:30
case ( 4 ) largestComponent
2013-06-06 00:40:37 +05:30
math_RtoQ ( 1 ) = R ( 2 , 1 ) - R ( 1 , 2 )
math_RtoQ ( 2 ) = R ( 1 , 3 ) + R ( 3 , 1 )
math_RtoQ ( 3 ) = R ( 2 , 3 ) + R ( 3 , 2 )
2011-12-01 17:31:13 +05:30
!4----------------------------------
2016-01-10 19:04:26 +05:30
end select largestComponent
2010-05-26 21:22:54 +05:30
2013-06-06 00:40:37 +05:30
max_absQ = 0.5_pReal * sqrt ( absQ ( largest ( 1 ) ) )
math_RtoQ = math_RtoQ * 0.25_pReal / max_absQ
2013-01-31 21:58:08 +05:30
math_RtoQ ( largest ( 1 ) ) = max_absQ
2012-08-25 17:16:36 +05:30
2013-01-31 21:58:08 +05:30
end function math_RtoQ
2007-03-21 15:50:25 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief rotation matrix from Euler angles (in radians)
2013-06-06 00:40:37 +05:30
!> @details rotation matrix is meant to represent a PASSIVE rotation,
!> @details composed of INTRINSIC rotations around the axes of the
!> @details rotating reference frame
!> @details (see http://en.wikipedia.org/wiki/Euler_angles for definitions)
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
pure function math_EulerToR ( Euler )
2010-03-18 17:53:17 +05:30
implicit none
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: Euler
real ( pReal ) , dimension ( 3 , 3 ) :: math_EulerToR
real ( pReal ) c1 , c , c2 , s1 , s , s2
2011-02-25 14:55:53 +05:30
C1 = cos ( Euler ( 1 ) )
2015-06-19 12:59:01 +05:30
C = cos ( Euler ( 2 ) )
2011-02-25 14:55:53 +05:30
C2 = cos ( Euler ( 3 ) )
S1 = sin ( Euler ( 1 ) )
2015-06-19 12:59:01 +05:30
S = sin ( Euler ( 2 ) )
2011-02-25 14:55:53 +05:30
S2 = sin ( Euler ( 3 ) )
2011-12-01 17:31:13 +05:30
2010-03-18 17:53:17 +05:30
math_EulerToR ( 1 , 1 ) = C1 * C2 - S1 * S2 * C
2013-06-06 00:40:37 +05:30
math_EulerToR ( 1 , 2 ) = - C1 * S2 - S1 * C2 * C
math_EulerToR ( 1 , 3 ) = S1 * S
math_EulerToR ( 2 , 1 ) = S1 * C2 + C1 * S2 * C
2010-03-18 17:53:17 +05:30
math_EulerToR ( 2 , 2 ) = - S1 * S2 + C1 * C2 * C
2013-06-06 00:40:37 +05:30
math_EulerToR ( 2 , 3 ) = - C1 * S
math_EulerToR ( 3 , 1 ) = S2 * S
math_EulerToR ( 3 , 2 ) = C2 * S
2010-03-18 17:53:17 +05:30
math_EulerToR ( 3 , 3 ) = C
2013-06-06 00:40:37 +05:30
math_EulerToR = transpose ( math_EulerToR ) ! convert to passive rotation
2012-08-25 17:16:36 +05:30
2012-03-09 01:55:28 +05:30
end function math_EulerToR
2010-03-18 17:53:17 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief quaternion (w+ix+jy+kz) from 3-1-3 Euler angles (in radians)
2013-06-06 00:40:37 +05:30
!> @details quaternion is meant to represent a PASSIVE rotation,
!> @details composed of INTRINSIC rotations around the axes of the
!> @details rotating reference frame
!> @details (see http://en.wikipedia.org/wiki/Euler_angles for definitions)
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2013-01-31 21:58:08 +05:30
pure function math_EulerToQ ( eulerangles )
2010-03-18 17:53:17 +05:30
implicit none
2010-05-06 19:37:21 +05:30
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: eulerangles
2013-01-31 21:58:08 +05:30
real ( pReal ) , dimension ( 4 ) :: math_EulerToQ
2010-05-06 19:37:21 +05:30
real ( pReal ) , dimension ( 3 ) :: halfangles
2011-12-01 17:31:13 +05:30
real ( pReal ) :: c , s
2012-08-25 17:16:36 +05:30
2010-05-06 19:37:21 +05:30
halfangles = 0.5_pReal * eulerangles
2012-08-25 17:16:36 +05:30
2011-02-25 14:55:53 +05:30
c = cos ( halfangles ( 2 ) )
s = sin ( halfangles ( 2 ) )
2012-08-25 17:16:36 +05:30
2016-01-10 19:04:26 +05:30
math_EulerToQ = [ cos ( halfangles ( 1 ) + halfangles ( 3 ) ) * c , &
cos ( halfangles ( 1 ) - halfangles ( 3 ) ) * s , &
sin ( halfangles ( 1 ) - halfangles ( 3 ) ) * s , &
sin ( halfangles ( 1 ) + halfangles ( 3 ) ) * c ]
2016-09-21 09:58:24 +05:30
math_EulerToQ = math_qConj ( math_EulerToQ ) ! convert to passive rotation
2012-08-25 17:16:36 +05:30
2013-01-31 21:58:08 +05:30
end function math_EulerToQ
2010-03-18 17:53:17 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief rotation matrix from axis and angle (in radians)
2013-06-06 00:40:37 +05:30
!> @details rotation matrix is meant to represent a ACTIVE rotation
!> @details (see http://en.wikipedia.org/wiki/Euler_angles for definitions)
!> @details formula for active rotation taken from http://mathworld.wolfram.com/RodriguesRotationFormula.html
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2013-06-06 00:40:37 +05:30
pure function math_axisAngleToR ( axis , omega )
2010-03-18 17:53:17 +05:30
implicit none
2013-06-06 00:40:37 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: math_axisAngleToR
2010-05-06 19:37:21 +05:30
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: axis
real ( pReal ) , intent ( in ) :: omega
real ( pReal ) , dimension ( 3 ) :: axisNrm
2011-12-01 17:31:13 +05:30
real ( pReal ) :: norm , s , c , c1
2010-05-06 19:37:21 +05:30
2016-06-30 14:41:35 +05:30
norm = norm2 ( axis )
2011-12-01 17:31:13 +05:30
if ( norm > 1.0e-8_pReal ) then ! non-zero rotation
2012-11-23 22:32:50 +05:30
axisNrm = axis / norm ! normalize axis to be sure
2010-05-06 19:37:21 +05:30
2011-02-25 14:55:53 +05:30
s = sin ( omega )
c = cos ( omega )
2010-05-06 19:37:21 +05:30
c1 = 1.0_pReal - c
2012-08-25 17:16:36 +05:30
2013-06-06 00:40:37 +05:30
math_axisAngleToR ( 1 , 1 ) = c + c1 * axisNrm ( 1 ) ** 2.0_pReal
math_axisAngleToR ( 1 , 2 ) = - s * axisNrm ( 3 ) + c1 * axisNrm ( 1 ) * axisNrm ( 2 )
math_axisAngleToR ( 1 , 3 ) = s * axisNrm ( 2 ) + c1 * axisNrm ( 1 ) * axisNrm ( 3 )
math_axisAngleToR ( 2 , 1 ) = s * axisNrm ( 3 ) + c1 * axisNrm ( 2 ) * axisNrm ( 1 )
math_axisAngleToR ( 2 , 2 ) = c + c1 * axisNrm ( 2 ) ** 2.0_pReal
math_axisAngleToR ( 2 , 3 ) = - s * axisNrm ( 1 ) + c1 * axisNrm ( 2 ) * axisNrm ( 3 )
math_axisAngleToR ( 3 , 1 ) = - s * axisNrm ( 2 ) + c1 * axisNrm ( 3 ) * axisNrm ( 1 )
math_axisAngleToR ( 3 , 2 ) = s * axisNrm ( 1 ) + c1 * axisNrm ( 3 ) * axisNrm ( 2 )
math_axisAngleToR ( 3 , 3 ) = c + c1 * axisNrm ( 3 ) ** 2.0_pReal
2010-05-06 19:37:21 +05:30
else
2013-06-06 00:40:37 +05:30
math_axisAngleToR = math_I3
2010-05-06 19:37:21 +05:30
endif
2010-03-18 17:53:17 +05:30
2013-06-06 00:40:37 +05:30
end function math_axisAngleToR
!--------------------------------------------------------------------------------------------------
!> @brief rotation matrix from axis and angle (in radians)
!> @details rotation matrix is meant to represent a PASSIVE rotation
!> @details (see http://en.wikipedia.org/wiki/Euler_angles for definitions)
!--------------------------------------------------------------------------------------------------
pure function math_EulerAxisAngleToR ( axis , omega )
implicit none
real ( pReal ) , dimension ( 3 , 3 ) :: math_EulerAxisAngleToR
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: axis
real ( pReal ) , intent ( in ) :: omega
2010-05-06 19:37:21 +05:30
2016-09-21 09:58:24 +05:30
math_EulerAxisAngleToR = transpose ( math_axisAngleToR ( axis , omega ) ) ! convert to passive rotation
2013-06-06 00:40:37 +05:30
end function math_EulerAxisAngleToR
!--------------------------------------------------------------------------------------------------
!> @brief quaternion (w+ix+jy+kz) from Euler axis and angle (in radians)
!> @details quaternion is meant to represent a PASSIVE rotation
!> @details (see http://en.wikipedia.org/wiki/Euler_angles for definitions)
!> @details formula for active rotation taken from
!> @details http://en.wikipedia.org/wiki/Rotation_representation_%28mathematics%29#Rodrigues_parameters
!--------------------------------------------------------------------------------------------------
pure function math_EulerAxisAngleToQ ( axis , omega )
implicit none
real ( pReal ) , dimension ( 4 ) :: math_EulerAxisAngleToQ
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: axis
real ( pReal ) , intent ( in ) :: omega
2016-09-21 09:58:24 +05:30
math_EulerAxisAngleToQ = math_qConj ( math_axisAngleToQ ( axis , omega ) ) ! convert to passive rotation
2013-06-06 00:40:37 +05:30
end function math_EulerAxisAngleToQ
2010-05-06 19:37:21 +05:30
2012-08-27 13:34:47 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief quaternion (w+ix+jy+kz) from axis and angle (in radians)
2013-06-06 00:40:37 +05:30
!> @details quaternion is meant to represent an ACTIVE rotation
!> @details (see http://en.wikipedia.org/wiki/Euler_angles for definitions)
!> @details formula for active rotation taken from
!> @details http://en.wikipedia.org/wiki/Rotation_representation_%28mathematics%29#Rodrigues_parameters
2012-08-27 13:34:47 +05:30
!--------------------------------------------------------------------------------------------------
2013-06-06 00:40:37 +05:30
pure function math_axisAngleToQ ( axis , omega )
2010-05-06 19:37:21 +05:30
implicit none
2013-06-06 00:40:37 +05:30
real ( pReal ) , dimension ( 4 ) :: math_axisAngleToQ
2010-05-06 19:37:21 +05:30
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: axis
real ( pReal ) , intent ( in ) :: omega
real ( pReal ) , dimension ( 3 ) :: axisNrm
2016-01-10 19:04:26 +05:30
real ( pReal ) :: norm
2010-05-06 19:37:21 +05:30
2011-02-25 14:55:53 +05:30
norm = sqrt ( math_mul3x3 ( axis , axis ) )
2016-01-10 19:04:26 +05:30
rotation : if ( norm > 1.0e-8_pReal ) then
2016-09-21 09:58:24 +05:30
axisNrm = axis / norm ! normalize axis to be sure
2016-01-10 19:04:26 +05:30
math_axisAngleToQ = [ cos ( 0.5_pReal * omega ) , sin ( 0.5_pReal * omega ) * axisNrm ( 1 : 3 ) ]
else rotation
math_axisAngleToQ = [ 1.0_pReal , 0.0_pReal , 0.0_pReal , 0.0_pReal ]
endif rotation
2010-05-06 19:37:21 +05:30
2013-06-06 00:40:37 +05:30
end function math_axisAngleToQ
2010-03-18 17:53:17 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief orientation matrix from quaternion (w+ix+jy+kz)
2013-06-06 00:40:37 +05:30
!> @details taken from http://arxiv.org/pdf/math/0701759v1.pdf
2013-09-19 01:24:39 +05:30
!> @details see also http://en.wikipedia.org/wiki/Rotation_formalisms_in_three_dimensions
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2013-06-06 00:40:37 +05:30
pure function math_qToR ( q )
2010-03-18 17:53:17 +05:30
implicit none
2013-06-06 00:40:37 +05:30
real ( pReal ) , dimension ( 4 ) , intent ( in ) :: q
2013-01-31 21:58:08 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: math_qToR , T , S
2011-12-01 17:31:13 +05:30
integer ( pInt ) :: i , j
2012-08-25 17:16:36 +05:30
2011-12-01 17:31:13 +05:30
forall ( i = 1_pInt : 3_pInt , j = 1_pInt : 3_pInt ) &
2013-06-06 00:40:37 +05:30
T ( i , j ) = q ( i + 1_pInt ) * q ( j + 1_pInt )
2016-01-10 19:04:26 +05:30
2013-06-06 00:40:37 +05:30
S = reshape ( [ 0.0_pReal , - q ( 4 ) , q ( 3 ) , &
q ( 4 ) , 0.0_pReal , - q ( 2 ) , &
2016-06-30 14:41:35 +05:30
- q ( 3 ) , q ( 2 ) , 0.0_pReal ] , [ 3 , 3 ] ) ! notation is transposed
2010-05-06 19:37:21 +05:30
2013-06-06 00:40:37 +05:30
math_qToR = ( 2.0_pReal * q ( 1 ) * q ( 1 ) - 1.0_pReal ) * math_I3 &
+ 2.0_pReal * T - 2.0_pReal * q ( 1 ) * S
2011-12-01 17:31:13 +05:30
2013-01-31 21:58:08 +05:30
end function math_qToR
2010-03-18 17:53:17 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief 3-1-3 Euler angles (in radians) from quaternion (w+ix+jy+kz)
2013-06-06 00:40:37 +05:30
!> @details quaternion is meant to represent a PASSIVE rotation,
!> @details composed of INTRINSIC rotations around the axes of the
!> @details rotating reference frame
!> @details (see http://en.wikipedia.org/wiki/Euler_angles for definitions)
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2013-06-06 00:40:37 +05:30
pure function math_qToEuler ( qPassive )
2010-03-18 17:53:17 +05:30
implicit none
2013-06-06 00:40:37 +05:30
real ( pReal ) , dimension ( 4 ) , intent ( in ) :: qPassive
real ( pReal ) , dimension ( 4 ) :: q
2013-01-31 21:58:08 +05:30
real ( pReal ) , dimension ( 3 ) :: math_qToEuler
2010-05-06 19:37:21 +05:30
2013-06-06 00:40:37 +05:30
q = math_qConj ( qPassive ) ! convert to active rotation, since formulas are defined for active rotations
2016-06-30 14:41:35 +05:30
math_qToEuler ( 2 ) = acos ( 1.0_pReal - 2.0_pReal * ( q ( 2 ) ** 2 + q ( 3 ) ** 2 ) )
2010-05-06 19:37:21 +05:30
2013-06-06 00:40:37 +05:30
if ( abs ( math_qToEuler ( 2 ) ) < 1.0e-6_pReal ) then
2016-01-10 19:04:26 +05:30
math_qToEuler ( 1 ) = sign ( 2.0_pReal * acos ( math_limit ( q ( 1 ) , - 1.0_pReal , 1.0_pReal ) ) , q ( 4 ) )
2013-01-31 21:58:08 +05:30
math_qToEuler ( 3 ) = 0.0_pReal
2010-05-26 21:22:54 +05:30
else
2016-06-30 14:41:35 +05:30
math_qToEuler ( 1 ) = atan2 ( + q ( 1 ) * q ( 3 ) + q ( 2 ) * q ( 4 ) , q ( 1 ) * q ( 2 ) - q ( 3 ) * q ( 4 ) )
2013-06-06 00:40:37 +05:30
math_qToEuler ( 3 ) = atan2 ( - q ( 1 ) * q ( 3 ) + q ( 2 ) * q ( 4 ) , q ( 1 ) * q ( 2 ) + q ( 3 ) * q ( 4 ) )
2010-05-26 21:22:54 +05:30
endif
2010-03-19 21:41:53 +05:30
2016-06-30 14:41:35 +05:30
math_qToEuler = merge ( math_qToEuler + [ 2.0_pReal * PI , PI , 2.0_pReal * PI ] , & ! ensure correct range
2016-01-10 19:04:26 +05:30
math_qToEuler , math_qToEuler < 0.0_pReal )
2010-03-19 21:41:53 +05:30
2013-01-31 21:58:08 +05:30
end function math_qToEuler
2010-03-18 17:53:17 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief axis-angle (x, y, z, ang in radians) from quaternion (w+ix+jy+kz)
2013-06-06 00:40:37 +05:30
!> @details quaternion is meant to represent an ACTIVE rotation
!> @details (see http://en.wikipedia.org/wiki/Euler_angles for definitions)
!> @details formula for active rotation taken from
!> @details http://en.wikipedia.org/wiki/Rotation_representation_%28mathematics%29#Rodrigues_parameters
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2013-01-31 21:58:08 +05:30
pure function math_qToAxisAngle ( Q )
2010-04-12 16:37:25 +05:30
implicit none
real ( pReal ) , dimension ( 4 ) , intent ( in ) :: Q
2011-12-01 17:31:13 +05:30
real ( pReal ) :: halfAngle , sinHalfAngle
2013-01-31 21:58:08 +05:30
real ( pReal ) , dimension ( 4 ) :: math_qToAxisAngle
2010-04-12 16:37:25 +05:30
2016-03-23 19:42:35 +05:30
halfAngle = acos ( math_limit ( Q ( 1 ) , - 1.0_pReal , 1.0_pReal ) )
2011-02-25 14:55:53 +05:30
sinHalfAngle = sin ( halfAngle )
2012-08-25 17:16:36 +05:30
2016-03-23 19:42:35 +05:30
smallRotation : if ( sinHalfAngle < = 1.0e-4_pReal ) then
2013-01-31 21:58:08 +05:30
math_qToAxisAngle = 0.0_pReal
2016-03-23 19:42:35 +05:30
else smallRotation
2016-01-10 19:04:26 +05:30
math_qToAxisAngle = [ Q ( 2 : 4 ) / sinHalfAngle , halfAngle * 2.0_pReal ]
2016-03-23 19:42:35 +05:30
endif smallRotation
2011-08-01 15:41:32 +05:30
2013-01-31 21:58:08 +05:30
end function math_qToAxisAngle
2010-05-06 19:37:21 +05:30
2010-04-12 16:37:25 +05:30
2013-06-06 00:40:37 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Euler axis-angle (x, y, z, ang in radians) from quaternion (w+ix+jy+kz)
!> @details quaternion is meant to represent a PASSIVE rotation
!> @details (see http://en.wikipedia.org/wiki/Euler_angles for definitions)
!--------------------------------------------------------------------------------------------------
pure function math_qToEulerAxisAngle ( qPassive )
implicit none
real ( pReal ) , dimension ( 4 ) , intent ( in ) :: qPassive
real ( pReal ) , dimension ( 4 ) :: q
real ( pReal ) , dimension ( 4 ) :: math_qToEulerAxisAngle
q = math_qConj ( qPassive ) ! convert to active rotation
math_qToEulerAxisAngle = math_qToAxisAngle ( q )
end function math_qToEulerAxisAngle
2012-08-27 13:34:47 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Rodrigues vector (x, y, z) from unit quaternion (w+ix+jy+kz)
!--------------------------------------------------------------------------------------------------
2013-01-31 21:58:08 +05:30
pure function math_qToRodrig ( Q )
2017-05-04 04:02:44 +05:30
use , intrinsic :: &
IEEE_arithmetic
2014-02-28 18:58:27 +05:30
use prec , only : &
2015-03-29 18:24:13 +05:30
tol_math_check
2010-04-28 22:49:58 +05:30
implicit none
real ( pReal ) , dimension ( 4 ) , intent ( in ) :: Q
2013-01-31 21:58:08 +05:30
real ( pReal ) , dimension ( 3 ) :: math_qToRodrig
2010-04-28 22:49:58 +05:30
2017-05-04 04:02:44 +05:30
math_qToRodrig = merge ( Q ( 2 : 4 ) / Q ( 1 ) , IEEE_value ( 1.0_pReal , IEEE_quiet_NaN ) , abs ( Q ( 1 ) ) > tol_math_check ) ! NaN for 180 deg since Rodrig is unbound
2010-04-28 22:49:58 +05:30
2013-01-31 21:58:08 +05:30
end function math_qToRodrig
2007-03-20 19:25:22 +05:30
2007-03-29 21:02:52 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief misorientation angle between two sets of Euler angles
!--------------------------------------------------------------------------------------------------
2013-01-31 21:58:08 +05:30
real ( pReal ) pure function math_EulerMisorientation ( EulerA , EulerB )
2007-03-29 21:02:52 +05:30
2007-03-20 19:25:22 +05:30
implicit none
2009-01-16 20:57:13 +05:30
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: EulerA , EulerB
2007-03-29 21:02:52 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: r
2013-01-31 21:58:08 +05:30
real ( pReal ) :: tr
2007-03-29 21:02:52 +05:30
2008-07-23 18:19:40 +05:30
r = math_mul33x33 ( math_EulerToR ( EulerB ) , transpose ( math_EulerToR ( EulerA ) ) )
2008-07-09 01:08:22 +05:30
2016-01-10 19:04:26 +05:30
tr = ( math_trace33 ( r ) - 1.0_pReal ) * 0.4999999_pReal
math_EulerMisorientation = abs ( 0.5_pReal * PI - asin ( tr ) )
2007-03-20 19:25:22 +05:30
2012-03-09 01:55:28 +05:30
end function math_EulerMisorientation
2010-05-06 19:37:21 +05:30
2007-03-20 19:25:22 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief draw a random sample from Euler space
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
function math_sampleRandomOri ( )
2007-03-21 15:50:25 +05:30
2007-03-20 19:25:22 +05:30
implicit none
2007-03-29 21:02:52 +05:30
real ( pReal ) , dimension ( 3 ) :: math_sampleRandomOri , rnd
2007-03-21 15:50:25 +05:30
2011-12-01 17:31:13 +05:30
call halton ( 3_pInt , rnd )
2016-01-10 19:04:26 +05:30
math_sampleRandomOri = [ rnd ( 1 ) * 2.0_pReal * PI , &
acos ( 2.0_pReal * rnd ( 2 ) - 1.0_pReal ) , &
rnd ( 3 ) * 2.0_pReal * PI ]
2007-03-20 19:25:22 +05:30
2012-03-09 01:55:28 +05:30
end function math_sampleRandomOri
2007-03-20 19:25:22 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief draw a random sample from Gauss component with noise (in radians) half-width
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
function math_sampleGaussOri ( center , noise )
2015-03-29 18:24:13 +05:30
use prec , only : &
tol_math_check
2007-03-21 15:50:25 +05:30
2007-03-20 19:25:22 +05:30
implicit none
2015-03-29 18:24:13 +05:30
real ( pReal ) , intent ( in ) :: noise
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: center
real ( pReal ) :: cosScatter , scatter
real ( pReal ) , dimension ( 3 ) :: math_sampleGaussOri , disturb
2014-02-28 18:58:27 +05:30
real ( pReal ) , dimension ( 3 ) , parameter :: ORIGIN = [ 0.0_pReal , 0.0_pReal , 0.0_pReal ]
2007-03-29 21:02:52 +05:30
real ( pReal ) , dimension ( 5 ) :: rnd
2015-03-29 18:24:13 +05:30
integer ( pInt ) :: i
2008-02-15 18:12:27 +05:30
2015-03-29 18:24:13 +05:30
if ( abs ( noise ) < tol_math_check ) then
2014-02-28 18:58:27 +05:30
math_sampleGaussOri = center
return
endif
2008-02-15 18:12:27 +05:30
2007-03-29 21:02:52 +05:30
! Helming uses different distribution with Bessel functions
! therefore the gauss scatter width has to be scaled differently
scatter = 0.95_pReal * noise
cosScatter = cos ( scatter )
do
2011-12-01 17:31:13 +05:30
call halton ( 5_pInt , rnd )
2012-08-25 17:16:36 +05:30
forall ( i = 1_pInt : 3_pInt ) rnd ( i ) = 2.0_pReal * rnd ( i ) - 1.0_pReal ! expand 1:3 to range [-1,+1]
2016-01-10 19:04:26 +05:30
disturb = [ scatter * rnd ( 1 ) , & ! phi1
sign ( 1.0_pReal , rnd ( 2 ) ) * acos ( cosScatter + ( 1.0_pReal - cosScatter ) * rnd ( 4 ) ) , & ! Phi
scatter * rnd ( 2 ) ] ! phi2
2014-02-28 18:58:27 +05:30
if ( rnd ( 5 ) < = exp ( - 1.0_pReal * ( math_EulerMisorientation ( ORIGIN , disturb ) / scatter ) ** 2_pReal ) ) exit
2008-07-09 01:08:22 +05:30
enddo
2008-07-23 18:19:40 +05:30
math_sampleGaussOri = math_RtoEuler ( math_mul33x33 ( math_EulerToR ( disturb ) , math_EulerToR ( center ) ) )
2012-08-25 17:16:36 +05:30
2012-03-09 01:55:28 +05:30
end function math_sampleGaussOri
2012-08-25 17:16:36 +05:30
2007-03-20 19:25:22 +05:30
2012-08-24 18:57:55 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief draw a random sample from Fiber component with noise (in radians)
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
function math_sampleFiberOri ( alpha , beta , noise )
2015-03-29 18:24:13 +05:30
use prec , only : &
tol_math_check
2007-03-21 15:50:25 +05:30
2007-03-20 19:25:22 +05:30
implicit none
2007-03-29 21:02:52 +05:30
real ( pReal ) , dimension ( 3 ) :: math_sampleFiberOri , fiberInC , fiberInS , axis
2012-08-24 18:57:55 +05:30
real ( pReal ) , dimension ( 2 ) , intent ( in ) :: alpha , beta
real ( pReal ) , dimension ( 6 ) :: rnd
2007-03-29 21:02:52 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: oRot , fRot , pRot
2011-12-01 17:31:13 +05:30
real ( pReal ) :: noise , scatter , cos2Scatter , angle
2015-05-05 12:07:59 +05:30
integer ( pInt ) , dimension ( 2 , 3 ) , parameter :: ROTMAP = reshape ( [ 2_pInt , 3_pInt , &
3_pInt , 1_pInt , &
1_pInt , 2_pInt ] , [ 2 , 3 ] )
2011-12-01 17:31:13 +05:30
integer ( pInt ) :: i
2007-03-21 15:50:25 +05:30
2007-03-29 21:02:52 +05:30
! Helming uses different distribution with Bessel functions
! therefore the gauss scatter width has to be scaled differently
scatter = 0.95_pReal * noise
cos2Scatter = cos ( 2.0_pReal * scatter )
! fiber axis in crystal coordinate system
2016-01-10 19:04:26 +05:30
fiberInC = [ sin ( alpha ( 1 ) ) * cos ( alpha ( 2 ) ) , &
sin ( alpha ( 1 ) ) * sin ( alpha ( 2 ) ) , &
cos ( alpha ( 1 ) ) ]
2007-03-29 21:02:52 +05:30
! fiber axis in sample coordinate system
2016-01-10 19:04:26 +05:30
fiberInS = [ sin ( beta ( 1 ) ) * cos ( beta ( 2 ) ) , &
sin ( beta ( 1 ) ) * sin ( beta ( 2 ) ) , &
cos ( beta ( 1 ) ) ]
2007-03-29 21:02:52 +05:30
! ---# rotation matrix from sample to crystal system #---
2011-02-25 14:55:53 +05:30
angle = - acos ( dot_product ( fiberInC , fiberInS ) )
2015-03-29 18:24:13 +05:30
if ( abs ( angle ) > tol_math_check ) then
2007-03-29 21:02:52 +05:30
! rotation axis between sample and crystal system (cross product)
2015-05-05 12:07:59 +05:30
forall ( i = 1_pInt : 3_pInt ) axis ( i ) = fiberInC ( ROTMAP ( 1 , i ) ) * fiberInS ( ROTMAP ( 2 , i ) ) - fiberInC ( ROTMAP ( 2 , i ) ) * fiberInS ( ROTMAP ( 1 , i ) )
2016-01-10 19:04:26 +05:30
oRot = math_EulerAxisAngleToR ( math_crossproduct ( fiberInC , fiberInS ) , angle )
2007-03-29 21:02:52 +05:30
else
oRot = math_I3
end if
! ---# rotation matrix about fiber axis (random angle) #---
2012-08-24 18:57:55 +05:30
do
call halton ( 6_pInt , rnd )
2013-06-06 00:40:37 +05:30
fRot = math_EulerAxisAngleToR ( fiberInS , rnd ( 1 ) * 2.0_pReal * pi )
2007-03-29 21:02:52 +05:30
! ---# rotation about random axis perpend to fiber #---
2012-08-24 18:57:55 +05:30
! random axis pependicular to fiber axis
axis ( 1 : 2 ) = rnd ( 2 : 3 )
2015-03-29 18:24:13 +05:30
if ( abs ( fiberInS ( 3 ) ) > tol_math_check ) then
2007-03-29 21:02:52 +05:30
axis ( 3 ) = - ( axis ( 1 ) * fiberInS ( 1 ) + axis ( 2 ) * fiberInS ( 2 ) ) / fiberInS ( 3 )
2015-03-29 18:24:13 +05:30
else if ( abs ( fiberInS ( 2 ) ) > tol_math_check ) then
2007-03-29 21:02:52 +05:30
axis ( 3 ) = axis ( 2 )
axis ( 2 ) = - ( axis ( 1 ) * fiberInS ( 1 ) + axis ( 3 ) * fiberInS ( 3 ) ) / fiberInS ( 2 )
2015-03-29 18:24:13 +05:30
else if ( abs ( fiberInS ( 1 ) ) > tol_math_check ) then
2007-03-29 21:02:52 +05:30
axis ( 3 ) = axis ( 1 )
axis ( 1 ) = - ( axis ( 2 ) * fiberInS ( 2 ) + axis ( 3 ) * fiberInS ( 3 ) ) / fiberInS ( 1 )
2012-08-24 18:57:55 +05:30
end if
2007-03-29 21:02:52 +05:30
2012-08-24 18:57:55 +05:30
! scattered rotation angle
if ( noise > 0.0_pReal ) then
angle = acos ( cos2Scatter + ( 1.0_pReal - cos2Scatter ) * rnd ( 4 ) )
if ( rnd ( 5 ) < = exp ( - 1.0_pReal * ( angle / scatter ) ** 2.0_pReal ) ) exit
else
angle = 0.0_pReal
exit
end if
2007-03-29 21:02:52 +05:30
enddo
2012-08-24 18:57:55 +05:30
if ( rnd ( 6 ) < = 0.5 ) angle = - angle
2013-06-06 00:40:37 +05:30
pRot = math_EulerAxisAngleToR ( axis , angle )
2007-03-29 21:02:52 +05:30
! ---# apply the three rotations #---
2012-08-24 18:57:55 +05:30
math_sampleFiberOri = math_RtoEuler ( math_mul33x33 ( pRot , math_mul33x33 ( fRot , oRot ) ) )
2008-07-09 01:08:22 +05:30
2012-03-09 01:55:28 +05:30
end function math_sampleFiberOri
2008-02-15 18:12:27 +05:30
2013-01-31 21:58:08 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief draw a random sample from Gauss variable
!--------------------------------------------------------------------------------------------------
real ( pReal ) function math_sampleGaussVar ( meanvalue , stddev , width )
2015-03-29 18:24:13 +05:30
use prec , only : &
tol_math_check
2013-01-31 21:58:08 +05:30
implicit none
real ( pReal ) , intent ( in ) :: meanvalue , & ! meanvalue of gauss distribution
stddev ! standard deviation of gauss distribution
real ( pReal ) , intent ( in ) , optional :: width ! width of considered values as multiples of standard deviation
real ( pReal ) , dimension ( 2 ) :: rnd ! random numbers
real ( pReal ) :: scatter , & ! normalized scatter around meanvalue
myWidth
2015-03-29 18:24:13 +05:30
if ( abs ( stddev ) < tol_math_check ) then
2013-01-31 21:58:08 +05:30
math_sampleGaussVar = meanvalue
return
endif
2015-05-05 12:07:59 +05:30
myWidth = merge ( width , 3.0_pReal , present ( width ) ) ! use +-3*sigma as default value for scatter if not given
2013-01-31 21:58:08 +05:30
do
call halton ( 2_pInt , rnd )
scatter = myWidth * ( 2.0_pReal * rnd ( 1 ) - 1.0_pReal )
2016-01-10 19:04:26 +05:30
if ( rnd ( 2 ) < = exp ( - 0.5_pReal * scatter ** 2.0_pReal ) ) exit ! test if scattered value is drawn
2013-01-31 21:58:08 +05:30
enddo
math_sampleGaussVar = scatter * stddev
end function math_sampleGaussVar
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2016-09-21 09:58:24 +05:30
!> @brief symmetrically equivalent Euler angles for given sample symmetry
!> @detail 1 (equivalent to != 2 and !=4):triclinic, 2:monoclinic, 4:orthotropic
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
pure function math_symmetricEulers ( sym , Euler )
2008-02-15 18:12:27 +05:30
implicit none
2016-09-21 09:58:24 +05:30
integer ( pInt ) , intent ( in ) :: sym !< symmetry Class
2008-02-15 18:12:27 +05:30
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: Euler
real ( pReal ) , dimension ( 3 , 3 ) :: math_symmetricEulers
2016-09-21 09:58:24 +05:30
math_symmetricEulers = transpose ( reshape ( [ PI + Euler ( 1 ) , PI - Euler ( 1 ) , 2.0_pReal * PI - Euler ( 1 ) , &
Euler ( 2 ) , PI - Euler ( 2 ) , PI - Euler ( 2 ) , &
Euler ( 3 ) , PI + Euler ( 3 ) , PI + Euler ( 3 ) ] , [ 3 , 3 ] ) ) ! transpose is needed to have symbolic notation instead of column-major
2008-02-15 18:12:27 +05:30
2016-09-21 09:58:24 +05:30
math_symmetricEulers = modulo ( math_symmetricEulers , 2.0_pReal * pi )
2008-02-15 18:12:27 +05:30
select case ( sym )
2016-09-21 09:58:24 +05:30
case ( 4_pInt ) ! orthotropic: all done
2008-02-15 18:12:27 +05:30
2016-09-21 09:58:24 +05:30
case ( 2_pInt ) ! monoclinic: return only first
2011-12-01 17:31:13 +05:30
math_symmetricEulers ( 1 : 3 , 2 : 3 ) = 0.0_pReal
2008-02-15 18:12:27 +05:30
2016-09-21 09:58:24 +05:30
case default ! triclinic: return blank
2008-02-15 18:12:27 +05:30
math_symmetricEulers = 0.0_pReal
end select
2012-03-09 01:55:28 +05:30
end function math_symmetricEulers
2008-02-15 18:12:27 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2016-01-31 16:55:26 +05:30
!> @brief eigenvalues and eigenvectors of symmetric matrix m
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2016-02-26 20:06:24 +05:30
subroutine math_eigenValuesVectorsSym ( m , values , vectors , error )
2011-09-14 18:56:00 +05:30
2012-08-25 17:16:36 +05:30
implicit none
2016-01-12 16:30:23 +05:30
real ( pReal ) , dimension ( : , : ) , intent ( in ) :: m
real ( pReal ) , dimension ( size ( m , 1 ) ) , intent ( out ) :: values
real ( pReal ) , dimension ( size ( m , 1 ) , size ( m , 1 ) ) , intent ( out ) :: vectors
2011-09-14 18:56:00 +05:30
logical , intent ( out ) :: error
2011-12-01 17:31:13 +05:30
integer ( pInt ) :: info
2016-01-12 16:30:23 +05:30
real ( pReal ) , dimension ( ( 64 + 2 ) * size ( m , 1 ) ) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f
2016-03-21 03:50:58 +05:30
external :: &
dsyev
2012-08-25 17:16:36 +05:30
2016-02-26 20:06:24 +05:30
vectors = m ! copy matrix to input (doubles as output) array
2016-01-12 16:30:23 +05:30
call dsyev ( 'V' , 'U' , size ( m , 1 ) , vectors , size ( m , 1 ) , values , work , ( 64 + 2 ) * size ( m , 1 ) , info )
2011-09-14 18:56:00 +05:30
error = ( info == 0_pInt )
2012-08-25 17:16:36 +05:30
2016-02-26 20:06:24 +05:30
end subroutine math_eigenValuesVectorsSym
2016-01-31 16:55:26 +05:30
!--------------------------------------------------------------------------------------------------
2016-02-27 00:51:47 +05:30
!> @brief eigenvalues and eigenvectors of symmetric 33 matrix m using an analytical expression
!> and the general LAPACK powered version for arbritrary sized matrices as fallback
!> @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)
2016-01-31 16:55:26 +05:30
!--------------------------------------------------------------------------------------------------
2016-02-26 20:06:24 +05:30
subroutine math_eigenValuesVectorsSym33 ( m , values , vectors )
2016-02-27 00:51:47 +05:30
2016-01-31 16:55:26 +05:30
implicit none
2016-02-27 00:51:47 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m
real ( pReal ) , dimension ( 3 ) , intent ( out ) :: values
real ( pReal ) , dimension ( 3 , 3 ) , intent ( out ) :: vectors
real ( pReal ) :: T , U , norm , threshold
2016-01-31 16:55:26 +05:30
logical :: error
2016-02-27 00:51:47 +05:30
values = math_eigenvaluesSym33 ( m )
vectors ( 1 : 3 , 2 ) = [ m ( 1 , 2 ) * m ( 2 , 3 ) - m ( 1 , 3 ) * m ( 2 , 2 ) , &
m ( 1 , 3 ) * m ( 1 , 2 ) - m ( 2 , 3 ) * m ( 1 , 1 ) , &
m ( 1 , 2 ) ** 2_pInt ]
2016-01-31 16:55:26 +05:30
2016-02-27 00:51:47 +05:30
T = maxval ( abs ( values ) )
U = max ( T , T ** 2_pInt )
threshold = sqrt ( 5.68e-14_pReal * U ** 2_pInt )
! Calculate first eigenvector by the formula v[0] = (m - lambda[0]).e1 x (m - lambda[0]).e2
vectors ( 1 : 3 , 1 ) = [ vectors ( 1 , 2 ) + m ( 1 , 3 ) * values ( 1 ) , &
vectors ( 2 , 2 ) + m ( 2 , 3 ) * values ( 1 ) , &
( m ( 1 , 1 ) - values ( 1 ) ) * ( m ( 2 , 2 ) - values ( 1 ) ) - vectors ( 3 , 2 ) ]
norm = norm2 ( vectors ( 1 : 3 , 1 ) )
fallback1 : if ( norm < threshold ) then
call math_eigenValuesVectorsSym ( m , values , vectors , error )
return
endif fallback1
vectors ( 1 : 3 , 1 ) = vectors ( 1 : 3 , 1 ) / norm
! Calculate second eigenvector by the formula v[1] = (m - lambda[1]).e1 x (m - lambda[1]).e2
vectors ( 1 : 3 , 2 ) = [ vectors ( 1 , 2 ) + m ( 1 , 3 ) * values ( 2 ) , &
vectors ( 2 , 2 ) + m ( 2 , 3 ) * values ( 2 ) , &
( m ( 1 , 1 ) - values ( 2 ) ) * ( m ( 2 , 2 ) - values ( 2 ) ) - vectors ( 3 , 2 ) ]
norm = norm2 ( vectors ( 1 : 3 , 2 ) )
fallback2 : if ( norm < threshold ) then
call math_eigenValuesVectorsSym ( m , values , vectors , error )
return
endif fallback2
vectors ( 1 : 3 , 2 ) = vectors ( 1 : 3 , 2 ) / norm
! Calculate third eigenvector according to v[2] = v[0] x v[1]
vectors ( 1 : 3 , 3 ) = math_crossproduct ( vectors ( 1 : 3 , 1 ) , vectors ( 1 : 3 , 2 ) )
2016-02-26 20:06:24 +05:30
end subroutine math_eigenValuesVectorsSym33
2016-01-31 16:55:26 +05:30
2016-02-26 20:06:24 +05:30
!--------------------------------------------------------------------------------------------------
2016-02-27 00:51:47 +05:30
!> @brief eigenvector basis of symmetric matrix m
2016-02-26 20:06:24 +05:30
!--------------------------------------------------------------------------------------------------
2016-02-27 00:51:47 +05:30
function math_eigenvectorBasisSym ( m )
2016-01-31 16:55:26 +05:30
2016-02-26 20:06:24 +05:30
implicit none
real ( pReal ) , dimension ( : , : ) , intent ( in ) :: m
real ( pReal ) , dimension ( size ( m , 1 ) ) :: values
real ( pReal ) , dimension ( size ( m , 1 ) , size ( m , 1 ) ) :: vectors
2016-02-27 00:51:47 +05:30
real ( pReal ) , dimension ( size ( m , 1 ) , size ( m , 1 ) ) :: math_eigenvectorBasisSym
2016-02-26 20:06:24 +05:30
logical :: error
integer ( pInt ) :: i
2016-01-31 16:55:26 +05:30
2016-02-27 00:51:47 +05:30
math_eigenvectorBasisSym = 0.0_pReal
2016-02-26 20:06:24 +05:30
call math_eigenValuesVectorsSym ( m , values , vectors , error )
if ( error ) return
2016-01-31 16:55:26 +05:30
2016-02-26 20:06:24 +05:30
do i = 1_pInt , size ( m , 1 )
2016-02-27 00:51:47 +05:30
math_eigenvectorBasisSym = math_eigenvectorBasisSym &
+ sqrt ( values ( i ) ) * math_tensorproduct ( vectors ( : , i ) , vectors ( : , i ) )
2016-02-26 20:06:24 +05:30
enddo
2016-01-31 16:55:26 +05:30
2016-02-27 00:51:47 +05:30
end function math_eigenvectorBasisSym
2016-01-31 16:55:26 +05:30
2016-02-26 20:06:24 +05:30
!--------------------------------------------------------------------------------------------------
2016-02-27 00:51:47 +05:30
!> @brief eigenvector basis of symmetric 33 matrix m
2016-02-26 20:06:24 +05:30
!--------------------------------------------------------------------------------------------------
2016-02-27 00:51:47 +05:30
function math_eigenvectorBasisSym33 ( m )
2016-02-26 20:06:24 +05:30
implicit none
2016-02-27 00:51:47 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: math_eigenvectorBasisSym33
real ( pReal ) , dimension ( 3 ) :: invariants , values
2016-02-26 21:05:55 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m
2016-02-27 00:07:10 +05:30
real ( pReal ) :: P , Q , rho , phi
2016-02-26 21:05:55 +05:30
real ( pReal ) , parameter :: TOL = 1.e-14_pReal
2016-02-27 00:51:47 +05:30
real ( pReal ) , dimension ( 3 , 3 , 3 ) :: N , EB
2016-02-26 20:06:24 +05:30
2016-02-26 21:05:55 +05:30
invariants = math_invariantsSym33 ( m )
2016-02-27 00:51:47 +05:30
EB = 0.0_pReal
2016-02-26 21:05:55 +05:30
2016-02-27 00:51:47 +05:30
P = invariants ( 2 ) - invariants ( 1 ) ** 2.0_pReal / 3.0_pReal
Q = - 2.0_pReal / 2 7.0_pReal * invariants ( 1 ) ** 3.0_pReal + product ( invariants ( 1 : 2 ) ) / 3.0_pReal - invariants ( 3 )
2016-02-26 21:05:55 +05:30
2016-02-27 00:51:47 +05:30
threeSimilarEigenvalues : if ( all ( abs ( [ P , Q ] ) < TOL ) ) then
2016-02-26 21:05:55 +05:30
values = invariants ( 1 ) / 3.0_pReal
2016-02-27 00:51:47 +05:30
! this is not really correct, but at least the basis is correct
EB ( 1 , 1 , 1 ) = 1.0_pReal
EB ( 2 , 2 , 2 ) = 1.0_pReal
EB ( 3 , 3 , 3 ) = 1.0_pReal
else threeSimilarEigenvalues
2016-02-26 21:05:55 +05:30
rho = sqrt ( - 3.0_pReal * P ** 3.0_pReal ) / 9.0_pReal
phi = acos ( math_limit ( - Q / rho * 0.5_pReal , - 1.0_pReal , 1.0_pReal ) )
2016-02-26 22:21:34 +05:30
values = 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 ) &
] + invariants ( 1 ) / 3.0_pReal
2016-02-27 00:51:47 +05:30
N ( 1 : 3 , 1 : 3 , 1 ) = m - values ( 1 ) * math_I3
N ( 1 : 3 , 1 : 3 , 2 ) = m - values ( 2 ) * math_I3
N ( 1 : 3 , 1 : 3 , 3 ) = m - values ( 3 ) * math_I3
twoSimilarEigenvalues : if ( abs ( values ( 1 ) - values ( 2 ) ) < TOL ) then
EB ( 1 : 3 , 1 : 3 , 3 ) = math_mul33x33 ( N ( 1 : 3 , 1 : 3 , 1 ) , N ( 1 : 3 , 1 : 3 , 2 ) ) / &
( ( values ( 3 ) - values ( 1 ) ) * ( values ( 3 ) - values ( 2 ) ) )
EB ( 1 : 3 , 1 : 3 , 1 ) = math_I3 - EB ( 1 : 3 , 1 : 3 , 3 )
elseif ( abs ( values ( 2 ) - values ( 3 ) ) < TOL ) then twoSimilarEigenvalues
EB ( 1 : 3 , 1 : 3 , 1 ) = math_mul33x33 ( N ( 1 : 3 , 1 : 3 , 2 ) , N ( 1 : 3 , 1 : 3 , 3 ) ) / &
( ( values ( 1 ) - values ( 2 ) ) * ( values ( 1 ) - values ( 3 ) ) )
EB ( 1 : 3 , 1 : 3 , 2 ) = math_I3 - EB ( 1 : 3 , 1 : 3 , 1 )
elseif ( abs ( values ( 3 ) - values ( 1 ) ) < TOL ) then twoSimilarEigenvalues
EB ( 1 : 3 , 1 : 3 , 2 ) = math_mul33x33 ( N ( 1 : 3 , 1 : 3 , 1 ) , N ( 1 : 3 , 1 : 3 , 3 ) ) / &
( ( values ( 2 ) - values ( 1 ) ) * ( values ( 2 ) - values ( 3 ) ) )
EB ( 1 : 3 , 1 : 3 , 1 ) = math_I3 - EB ( 1 : 3 , 1 : 3 , 2 )
else twoSimilarEigenvalues
EB ( 1 : 3 , 1 : 3 , 1 ) = math_mul33x33 ( N ( 1 : 3 , 1 : 3 , 2 ) , N ( 1 : 3 , 1 : 3 , 3 ) ) / &
( ( values ( 1 ) - values ( 2 ) ) * ( values ( 1 ) - values ( 3 ) ) )
EB ( 1 : 3 , 1 : 3 , 2 ) = math_mul33x33 ( N ( 1 : 3 , 1 : 3 , 1 ) , N ( 1 : 3 , 1 : 3 , 3 ) ) / &
( ( values ( 2 ) - values ( 1 ) ) * ( values ( 2 ) - values ( 3 ) ) )
EB ( 1 : 3 , 1 : 3 , 3 ) = math_mul33x33 ( N ( 1 : 3 , 1 : 3 , 1 ) , N ( 1 : 3 , 1 : 3 , 2 ) ) / &
( ( values ( 3 ) - values ( 1 ) ) * ( values ( 3 ) - values ( 2 ) ) )
endif twoSimilarEigenvalues
endif threeSimilarEigenvalues
math_eigenvectorBasisSym33 = sqrt ( values ( 1 ) ) * EB ( 1 : 3 , 1 : 3 , 1 ) &
+ sqrt ( values ( 2 ) ) * EB ( 1 : 3 , 1 : 3 , 2 ) &
+ sqrt ( values ( 3 ) ) * EB ( 1 : 3 , 1 : 3 , 3 )
end function math_eigenvectorBasisSym33
2011-08-01 15:41:32 +05:30
2013-01-31 21:58:08 +05:30
!--------------------------------------------------------------------------------------------------
2016-02-26 20:06:24 +05:30
!> @brief rotational part from polar decomposition of 33 tensor m
2013-01-31 21:58:08 +05:30
!--------------------------------------------------------------------------------------------------
2016-01-13 03:09:31 +05:30
function math_rotationalPart33 ( m )
2016-05-29 14:15:03 +05:30
use prec , only : &
2016-06-30 14:41:35 +05:30
dEq0
2016-01-13 03:09:31 +05:30
use IO , only : &
IO_warning
2013-01-31 21:58:08 +05:30
implicit none
2016-01-13 03:09:31 +05:30
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: m
real ( pReal ) , dimension ( 3 , 3 ) :: math_rotationalPart33
2016-02-26 20:06:24 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: U , Uinv
2016-01-10 19:04:26 +05:30
2016-02-27 00:51:47 +05:30
U = math_eigenvectorBasisSym33 ( math_mul33x33 ( transpose ( m ) , m ) )
2016-01-10 19:04:26 +05:30
Uinv = math_inv33 ( U )
2016-02-26 20:06:24 +05:30
2016-06-30 14:41:35 +05:30
inversionFailed : if ( all ( dEq0 ( Uinv ) ) ) then
2016-01-13 03:09:31 +05:30
math_rotationalPart33 = math_I3
call IO_warning ( 650_pInt )
2016-05-29 14:15:03 +05:30
else inversionFailed
2016-01-13 03:09:31 +05:30
math_rotationalPart33 = math_mul33x33 ( m , Uinv )
2016-05-29 14:15:03 +05:30
endif inversionFailed
2013-01-31 21:58:08 +05:30
2016-01-13 03:09:31 +05:30
end function math_rotationalPart33
2007-03-20 19:25:22 +05:30
2011-12-01 17:31:13 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2016-02-02 02:07:27 +05:30
!> @brief Eigenvalues of symmetric matrix m
2017-05-04 04:02:44 +05:30
! will return NaN on error
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2016-02-02 02:07:27 +05:30
function math_eigenvaluesSym ( m )
2017-05-04 04:02:44 +05:30
use , intrinsic :: &
IEEE_arithmetic
2016-01-10 19:04:26 +05:30
2011-08-26 19:36:37 +05:30
implicit none
2016-01-12 16:30:23 +05:30
real ( pReal ) , dimension ( : , : ) , intent ( in ) :: m
2016-02-02 02:07:27 +05:30
real ( pReal ) , dimension ( size ( m , 1 ) ) :: math_eigenvaluesSym
2016-01-12 16:30:23 +05:30
real ( pReal ) , dimension ( size ( m , 1 ) , size ( m , 1 ) ) :: vectors
2016-01-10 19:04:26 +05:30
integer ( pInt ) :: info
2016-01-12 16:30:23 +05:30
real ( pReal ) , dimension ( ( 64 + 2 ) * size ( m , 1 ) ) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f
2016-03-21 03:50:58 +05:30
external :: &
dsyev
2012-08-25 17:16:36 +05:30
2016-01-12 16:30:23 +05:30
vectors = m ! copy matrix to input (doubles as output) array
2016-02-02 02:07:27 +05:30
call dsyev ( 'N' , 'U' , size ( m , 1 ) , vectors , size ( m , 1 ) , math_eigenvaluesSym , work , ( 64 + 2 ) * size ( m , 1 ) , info )
2017-05-04 04:02:44 +05:30
if ( info / = 0_pInt ) math_eigenvaluesSym = IEEE_value ( 1.0_pReal , IEEE_quiet_NaN )
2016-01-10 19:04:26 +05:30
2016-02-02 02:07:27 +05:30
end function math_eigenvaluesSym
2007-03-20 19:25:22 +05:30
2016-02-02 13:39:07 +05:30
2016-02-01 12:18:42 +05:30
!--------------------------------------------------------------------------------------------------
2016-02-02 23:29:04 +05:30
!> @brief eigenvalues of symmetric 33 matrix m using an analytical expression
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @details similar to http://arxiv.org/abs/physics/0610206 (DSYEVC3)
!> but apparently more stable solution and has general LAPACK powered version for arbritrary sized
!> matrices as fallback
2016-02-01 12:18:42 +05:30
!--------------------------------------------------------------------------------------------------
2016-02-02 02:07:27 +05:30
function math_eigenvaluesSym33 ( m )
implicit none
2016-02-02 13:39:07 +05:30
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: m
2016-02-02 23:29:04 +05:30
real ( pReal ) , dimension ( 3 ) :: math_eigenvaluesSym33 , invariants
real ( pReal ) :: P , Q , rho , phi
2016-02-02 13:39:07 +05:30
real ( pReal ) , parameter :: TOL = 1.e-14_pReal
2016-02-02 02:07:27 +05:30
2016-02-02 23:29:04 +05:30
invariants = math_invariantsSym33 ( m ) ! invariants are coefficients in characteristic polynomial apart for the sign of c0 and c2 in http://arxiv.org/abs/physics/0610206
2016-02-02 17:53:45 +05:30
2016-02-02 23:29:04 +05:30
P = invariants ( 2 ) - invariants ( 1 ) ** 2.0_pReal / 3.0_pReal ! different from http://arxiv.org/abs/physics/0610206 (this formulation was in DAMASK)
Q = - 2.0_pReal / 2 7.0_pReal * invariants ( 1 ) ** 3.0_pReal + product ( invariants ( 1 : 2 ) ) / 3.0_pReal - invariants ( 3 ) ! different from http://arxiv.org/abs/physics/0610206 (this formulation was in DAMASK)
2016-02-02 02:07:27 +05:30
2016-02-27 00:51:47 +05:30
if ( all ( abs ( [ P , Q ] ) < TOL ) ) then
2016-02-02 23:29:04 +05:30
math_eigenvaluesSym33 = math_eigenvaluesSym ( m )
2016-02-02 13:39:07 +05:30
else
rho = sqrt ( - 3.0_pReal * P ** 3.0_pReal ) / 9.0_pReal
2016-02-02 23:29:04 +05:30
phi = acos ( math_limit ( - Q / rho * 0.5_pReal , - 1.0_pReal , 1.0_pReal ) )
2016-02-02 13:39:07 +05:30
math_eigenvaluesSym33 = 2.0_pReal * rho ** ( 1.0_pReal / 3.0_pReal ) * &
[ cos ( phi / 3.0_pReal ) , &
2016-02-02 23:29:04 +05:30
cos ( ( phi + 2.0_pReal * PI ) / 3.0_pReal ) , &
cos ( ( phi + 4.0_pReal * PI ) / 3.0_pReal ) &
] + invariants ( 1 ) / 3.0_pReal
2016-02-02 13:39:07 +05:30
endif
2016-02-02 02:07:27 +05:30
end function math_eigenvaluesSym33
2011-12-01 17:31:13 +05:30
2016-02-02 17:53:45 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2016-02-02 17:53:45 +05:30
!> @brief invariants of symmetrix 33 matrix m
2016-02-02 14:14:51 +05:30
!--------------------------------------------------------------------------------------------------
pure function math_invariantsSym33 ( m )
implicit none
2016-02-02 17:53:45 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m
2016-02-02 14:14:51 +05:30
real ( pReal ) , dimension ( 3 ) :: math_invariantsSym33
math_invariantsSym33 ( 1 ) = math_trace33 ( m )
2016-02-02 17:53:45 +05:30
math_invariantsSym33 ( 2 ) = m ( 1 , 1 ) * m ( 2 , 2 ) + m ( 1 , 1 ) * m ( 3 , 3 ) + m ( 2 , 2 ) * m ( 3 , 3 ) &
- ( m ( 1 , 2 ) ** 2 + m ( 1 , 3 ) ** 2 + m ( 2 , 3 ) ** 2 )
math_invariantsSym33 ( 3 ) = math_detSym33 ( m )
2016-02-02 14:14:51 +05:30
end function math_invariantsSym33
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2013-03-07 01:04:30 +05:30
!> @brief computes the next element in the Halton sequence.
!> @author John Burkardt
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
subroutine halton ( ndim , r )
2013-03-07 01:04:30 +05:30
2007-03-20 19:25:22 +05:30
implicit none
2013-03-07 01:04:30 +05:30
integer ( pInt ) , intent ( in ) :: ndim !< dimension of the element
real ( pReal ) , intent ( out ) , dimension ( ndim ) :: r !< next element of the current Halton sequence
2011-12-01 17:31:13 +05:30
integer ( pInt ) , dimension ( ndim ) :: base
integer ( pInt ) :: seed
integer ( pInt ) , dimension ( 1 ) :: value_halton
2007-03-20 19:25:22 +05:30
2011-12-01 17:31:13 +05:30
call halton_memory ( 'GET' , 'SEED' , 1_pInt , value_halton )
seed = value_halton ( 1 )
2007-03-20 19:25:22 +05:30
2011-12-01 17:31:13 +05:30
call halton_memory ( 'GET' , 'BASE' , ndim , base )
2007-03-20 19:25:22 +05:30
2011-12-01 17:31:13 +05:30
call i_to_halton ( seed , base , ndim , r )
2007-03-20 19:25:22 +05:30
2011-12-01 17:31:13 +05:30
value_halton ( 1 ) = 1_pInt
call halton_memory ( 'INC' , 'SEED' , 1_pInt , value_halton )
2007-03-21 15:50:25 +05:30
2017-08-12 09:33:40 +05:30
!--------------------------------------------------------------------------------------------------
contains
!-------------------------------------------------------------------------------------------------
!> @brief computes an element of a Halton sequence.
!> @details Only the absolute value of SEED is considered. SEED = 0 is allowed, and returns R = 0.
!> @details Halton Bases should be distinct prime numbers. This routine only checks that each base
!> @details is greater than 1.
!> @details Reference:
!> @details J.H. Halton: On the efficiency of certain quasi-random sequences of points in evaluating
!> @details multi-dimensional integrals, Numerische Mathematik, Volume 2, pages 84-90, 1960.
!> @author John Burkardt
!-------------------------------------------------------------------------------------------------
subroutine i_to_halton ( seed , base , ndim , r )
use IO , only : &
IO_error
implicit none
integer ( pInt ) , intent ( in ) :: &
ndim , & !< dimension of the sequence
seed !< index of the desired element
integer ( pInt ) , intent ( in ) , dimension ( ndim ) :: base !< Halton bases
real ( pReal ) , intent ( out ) , dimension ( ndim ) :: r !< the SEED-th element of the Halton sequence for the given bases
real ( pReal ) , dimension ( ndim ) :: base_inv
integer ( pInt ) , dimension ( ndim ) :: &
digit , &
seed2
seed2 = abs ( seed )
r = 0.0_pReal
if ( any ( base ( 1 : ndim ) < = 1_pInt ) ) call IO_error ( error_ID = 405_pInt )
base_inv ( 1 : ndim ) = 1.0_pReal / real ( base ( 1 : ndim ) , pReal )
do while ( any ( seed2 ( 1 : ndim ) / = 0_pInt ) )
digit ( 1 : ndim ) = mod ( seed2 ( 1 : ndim ) , base ( 1 : ndim ) )
r ( 1 : ndim ) = r ( 1 : ndim ) + real ( digit ( 1 : ndim ) , pReal ) * base_inv ( 1 : ndim )
base_inv ( 1 : ndim ) = base_inv ( 1 : ndim ) / real ( base ( 1 : ndim ) , pReal )
seed2 ( 1 : ndim ) = seed2 ( 1 : ndim ) / base ( 1 : ndim )
enddo
end subroutine i_to_halton
2012-03-09 01:55:28 +05:30
end subroutine halton
2007-03-21 15:50:25 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2013-03-07 01:04:30 +05:30
!> @brief sets or returns quantities associated with the Halton sequence.
!> @details If action_halton is 'SET' and action_halton is 'BASE', then NDIM is input, and
!> @details is the number of entries in value_halton to be put into BASE.
!> @details If action_halton is 'SET', then on input, value_halton contains values to be assigned
!> @details to the internal variable.
!> @details If action_halton is 'GET', then on output, value_halton contains the values of
!> @details the specified internal variable.
!> @details If action_halton is 'INC', then on input, value_halton contains the increment to
!> @details be added to the specified internal variable.
!> @author John Burkardt
2012-08-27 13:34:47 +05:30
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
subroutine halton_memory ( action_halton , name_halton , ndim , value_halton )
2017-08-12 09:33:40 +05:30
use IO , only : &
IO_lc
2007-03-21 15:50:25 +05:30
2014-02-28 18:58:27 +05:30
implicit none
2013-03-07 01:04:30 +05:30
character ( len = * ) , intent ( in ) :: &
action_halton , & !< desired action: GET the value of a particular quantity, SET the value of a particular quantity, INC the value of a particular quantity (only for SEED)
name_halton !< name of the quantity: BASE: Halton base(s), NDIM: spatial dimension, SEED: current Halton seed
2011-12-01 17:31:13 +05:30
integer ( pInt ) , dimension ( * ) , intent ( inout ) :: value_halton
integer ( pInt ) , allocatable , save , dimension ( : ) :: base
2007-03-20 19:25:22 +05:30
logical , save :: first_call = . true .
2013-03-07 01:04:30 +05:30
integer ( pInt ) , intent ( in ) :: ndim !< dimension of the quantity
2011-12-01 17:31:13 +05:30
integer ( pInt ) , save :: ndim_save = 0_pInt , seed = 1_pInt
2017-08-12 09:33:40 +05:30
integer ( pInt ) :: i
2011-12-01 17:31:13 +05:30
if ( first_call ) then
ndim_save = 1_pInt
allocate ( base ( ndim_save ) )
base ( 1 ) = 2_pInt
2007-03-20 19:25:22 +05:30
first_call = . false .
2011-12-01 17:31:13 +05:30
endif
2013-03-07 01:04:30 +05:30
!--------------------------------------------------------------------------------------------------
! Set
2017-08-12 09:33:40 +05:30
actionHalton : if ( IO_lc ( action_halton ( 1 : 1 ) ) == 's' ) then
2007-03-20 19:25:22 +05:30
2017-08-12 09:33:40 +05:30
nameSet : if ( IO_lc ( name_halton ( 1 : 1 ) ) == 'b' ) then
if ( ndim_save / = ndim ) ndim_save = ndim
base = value_halton ( 1 : ndim )
elseif ( IO_lc ( name_halton ( 1 : 1 ) ) == 'n' ) then nameSet
2011-12-01 17:31:13 +05:30
if ( ndim_save / = value_halton ( 1 ) ) then
ndim_save = value_halton ( 1 )
2017-08-12 09:33:40 +05:30
base = [ ( prime ( i ) , i = 1_pInt , ndim_save ) ]
2007-03-20 19:25:22 +05:30
else
2011-12-01 17:31:13 +05:30
ndim_save = value_halton ( 1 )
endif
2017-08-12 09:33:40 +05:30
elseif ( IO_lc ( name_halton ( 1 : 1 ) ) == 's' ) then nameSet
2011-12-01 17:31:13 +05:30
seed = value_halton ( 1 )
2017-08-12 09:33:40 +05:30
endif nameSet
2013-03-07 01:04:30 +05:30
!--------------------------------------------------------------------------------------------------
! Get
2017-08-12 09:33:40 +05:30
elseif ( IO_lc ( action_halton ( 1 : 1 ) ) == 'g' ) then actionHalton
nameGet : if ( IO_lc ( name_halton ( 1 : 1 ) ) == 'b' ) then
2011-12-01 17:31:13 +05:30
if ( ndim / = ndim_save ) then
2017-08-12 09:33:40 +05:30
ndim_save = ndim
base = [ ( prime ( i ) , i = 1_pInt , ndim_save ) ]
2011-12-01 17:31:13 +05:30
endif
value_halton ( 1 : ndim_save ) = base ( 1 : ndim_save )
2017-08-12 09:33:40 +05:30
elseif ( IO_lc ( name_halton ( 1 : 1 ) ) == 'n' ) then nameGet
2011-12-01 17:31:13 +05:30
value_halton ( 1 ) = ndim_save
2017-08-12 09:33:40 +05:30
elseif ( IO_lc ( name_halton ( 1 : 1 ) ) == 's' ) then nameGet
2011-12-01 17:31:13 +05:30
value_halton ( 1 ) = seed
2017-08-12 09:33:40 +05:30
endif nameGet
2013-03-07 01:04:30 +05:30
!--------------------------------------------------------------------------------------------------
! Increment
2017-08-12 09:33:40 +05:30
elseif ( IO_lc ( action_halton ( 1 : 1 ) ) == 'i' ) then actionHalton
if ( IO_lc ( name_halton ( 1 : 1 ) ) == 's' ) seed = seed + value_halton ( 1 )
endif actionHalton
2011-12-01 17:31:13 +05:30
2017-08-12 10:05:44 +05:30
!--------------------------------------------------------------------------------------------------
contains
!--------------------------------------------------------------------------------------------------
!> @brief returns any of the first 1500 prime numbers.
!> @details n = 0 is legal, returning PRIME = 1.
!> @details Reference:
!> @details Milton Abramowitz and Irene Stegun: Handbook of Mathematical Functions,
!> @details US Department of Commerce, 1964, pages 870-873.
!> @details Daniel Zwillinger: CRC Standard Mathematical Tables and Formulae,
!> @details 30th Edition, CRC Press, 1996, pages 95-98.
!> @author John Burkardt
!--------------------------------------------------------------------------------------------------
integer ( pInt ) function prime ( n )
use IO , only : &
IO_error
implicit none
integer ( pInt ) , intent ( in ) :: n !< index of the desired prime number
integer ( pInt ) , dimension ( 0 : 1500 ) , parameter :: &
npvec = int ( [ &
1 , &
2 , 3 , 5 , 7 , 11 , 13 , 17 , 19 , 23 , 29 , &
31 , 37 , 41 , 43 , 47 , 53 , 59 , 61 , 67 , 71 , &
73 , 79 , 83 , 89 , 97 , 101 , 103 , 107 , 109 , 113 , &
127 , 131 , 137 , 139 , 149 , 151 , 157 , 163 , 167 , 173 , &
179 , 181 , 191 , 193 , 197 , 199 , 211 , 223 , 227 , 229 , &
233 , 239 , 241 , 251 , 257 , 263 , 269 , 271 , 277 , 281 , &
283 , 293 , 307 , 311 , 313 , 317 , 331 , 337 , 347 , 349 , &
353 , 359 , 367 , 373 , 379 , 383 , 389 , 397 , 401 , 409 , &
419 , 421 , 431 , 433 , 439 , 443 , 449 , 457 , 461 , 463 , &
467 , 479 , 487 , 491 , 499 , 503 , 509 , 521 , 523 , 541 , &
! 101:200
547 , 557 , 563 , 569 , 571 , 577 , 587 , 593 , 599 , 601 , &
607 , 613 , 617 , 619 , 631 , 641 , 643 , 647 , 653 , 659 , &
661 , 673 , 677 , 683 , 691 , 701 , 709 , 719 , 727 , 733 , &
739 , 743 , 751 , 757 , 761 , 769 , 773 , 787 , 797 , 809 , &
811 , 821 , 823 , 827 , 829 , 839 , 853 , 857 , 859 , 863 , &
877 , 881 , 883 , 887 , 907 , 911 , 919 , 929 , 937 , 941 , &
947 , 953 , 967 , 971 , 977 , 983 , 991 , 997 , 1009 , 1013 , &
1019 , 1021 , 1031 , 1033 , 1039 , 1049 , 1051 , 1061 , 1063 , 1069 , &
1087 , 1091 , 1093 , 1097 , 1103 , 1109 , 1117 , 1123 , 1129 , 1151 , &
1153 , 1163 , 1171 , 1181 , 1187 , 1193 , 1201 , 1213 , 1217 , 1223 , &
! 201:300
1229 , 1231 , 1237 , 1249 , 1259 , 1277 , 1279 , 1283 , 1289 , 1291 , &
1297 , 1301 , 1303 , 1307 , 1319 , 1321 , 1327 , 1361 , 1367 , 1373 , &
1381 , 1399 , 1409 , 1423 , 1427 , 1429 , 1433 , 1439 , 1447 , 1451 , &
1453 , 1459 , 1471 , 1481 , 1483 , 1487 , 1489 , 1493 , 1499 , 1511 , &
1523 , 1531 , 1543 , 1549 , 1553 , 1559 , 1567 , 1571 , 1579 , 1583 , &
1597 , 1601 , 1607 , 1609 , 1613 , 1619 , 1621 , 1627 , 1637 , 1657 , &
1663 , 1667 , 1669 , 1693 , 1697 , 1699 , 1709 , 1721 , 1723 , 1733 , &
1741 , 1747 , 1753 , 1759 , 1777 , 1783 , 1787 , 1789 , 1801 , 1811 , &
1823 , 1831 , 1847 , 1861 , 1867 , 1871 , 1873 , 1877 , 1879 , 1889 , &
1901 , 1907 , 1913 , 1931 , 1933 , 1949 , 1951 , 1973 , 1979 , 1987 , &
! 301:400
1993 , 1997 , 1999 , 2003 , 2011 , 2017 , 2027 , 2029 , 2039 , 2053 , &
2063 , 2069 , 2081 , 2083 , 2087 , 2089 , 2099 , 2111 , 2113 , 2129 , &
2131 , 2137 , 2141 , 2143 , 2153 , 2161 , 2179 , 2203 , 2207 , 2213 , &
2221 , 2237 , 2239 , 2243 , 2251 , 2267 , 2269 , 2273 , 2281 , 2287 , &
2293 , 2297 , 2309 , 2311 , 2333 , 2339 , 2341 , 2347 , 2351 , 2357 , &
2371 , 2377 , 2381 , 2383 , 2389 , 2393 , 2399 , 2411 , 2417 , 2423 , &
2437 , 2441 , 2447 , 2459 , 2467 , 2473 , 2477 , 2503 , 2521 , 2531 , &
2539 , 2543 , 2549 , 2551 , 2557 , 2579 , 2591 , 2593 , 2609 , 2617 , &
2621 , 2633 , 2647 , 2657 , 2659 , 2663 , 2671 , 2677 , 2683 , 2687 , &
2689 , 2693 , 2699 , 2707 , 2711 , 2713 , 2719 , 2729 , 2731 , 2741 , &
! 401:500
2749 , 2753 , 2767 , 2777 , 2789 , 2791 , 2797 , 2801 , 2803 , 2819 , &
2833 , 2837 , 2843 , 2851 , 2857 , 2861 , 2879 , 2887 , 2897 , 2903 , &
2909 , 2917 , 2927 , 2939 , 2953 , 2957 , 2963 , 2969 , 2971 , 2999 , &
3001 , 3011 , 3019 , 3023 , 3037 , 3041 , 3049 , 3061 , 3067 , 3079 , &
3083 , 3089 , 3109 , 3119 , 3121 , 3137 , 3163 , 3167 , 3169 , 3181 , &
3187 , 3191 , 3203 , 3209 , 3217 , 3221 , 3229 , 3251 , 3253 , 3257 , &
3259 , 3271 , 3299 , 3301 , 3307 , 3313 , 3319 , 3323 , 3329 , 3331 , &
3343 , 3347 , 3359 , 3361 , 3371 , 3373 , 3389 , 3391 , 3407 , 3413 , &
3433 , 3449 , 3457 , 3461 , 3463 , 3467 , 3469 , 3491 , 3499 , 3511 , &
3517 , 3527 , 3529 , 3533 , 3539 , 3541 , 3547 , 3557 , 3559 , 3571 , &
! 501:600
3581 , 3583 , 3593 , 3607 , 3613 , 3617 , 3623 , 3631 , 3637 , 3643 , &
3659 , 3671 , 3673 , 3677 , 3691 , 3697 , 3701 , 3709 , 3719 , 3727 , &
3733 , 3739 , 3761 , 3767 , 3769 , 3779 , 3793 , 3797 , 3803 , 3821 , &
3823 , 3833 , 3847 , 3851 , 3853 , 3863 , 3877 , 3881 , 3889 , 3907 , &
3911 , 3917 , 3919 , 3923 , 3929 , 3931 , 3943 , 3947 , 3967 , 3989 , &
4001 , 4003 , 4007 , 4013 , 4019 , 4021 , 4027 , 4049 , 4051 , 4057 , &
4073 , 4079 , 4091 , 4093 , 4099 , 4111 , 4127 , 4129 , 4133 , 4139 , &
4153 , 4157 , 4159 , 4177 , 4201 , 4211 , 4217 , 4219 , 4229 , 4231 , &
4241 , 4243 , 4253 , 4259 , 4261 , 4271 , 4273 , 4283 , 4289 , 4297 , &
4327 , 4337 , 4339 , 4349 , 4357 , 4363 , 4373 , 4391 , 4397 , 4409 , &
! 601:700
4421 , 4423 , 4441 , 4447 , 4451 , 4457 , 4463 , 4481 , 4483 , 4493 , &
4507 , 4513 , 4517 , 4519 , 4523 , 4547 , 4549 , 4561 , 4567 , 4583 , &
4591 , 4597 , 4603 , 4621 , 4637 , 4639 , 4643 , 4649 , 4651 , 4657 , &
4663 , 4673 , 4679 , 4691 , 4703 , 4721 , 4723 , 4729 , 4733 , 4751 , &
4759 , 4783 , 4787 , 4789 , 4793 , 4799 , 4801 , 4813 , 4817 , 4831 , &
4861 , 4871 , 4877 , 4889 , 4903 , 4909 , 4919 , 4931 , 4933 , 4937 , &
4943 , 4951 , 4957 , 4967 , 4969 , 4973 , 4987 , 4993 , 4999 , 5003 , &
5009 , 5011 , 5021 , 5023 , 5039 , 5051 , 5059 , 5077 , 5081 , 5087 , &
5099 , 5101 , 5107 , 5113 , 5119 , 5147 , 5153 , 5167 , 5171 , 5179 , &
5189 , 5197 , 5209 , 5227 , 5231 , 5233 , 5237 , 5261 , 5273 , 5279 , &
! 701:800
5281 , 5297 , 5303 , 5309 , 5323 , 5333 , 5347 , 5351 , 5381 , 5387 , &
5393 , 5399 , 5407 , 5413 , 5417 , 5419 , 5431 , 5437 , 5441 , 5443 , &
5449 , 5471 , 5477 , 5479 , 5483 , 5501 , 5503 , 5507 , 5519 , 5521 , &
5527 , 5531 , 5557 , 5563 , 5569 , 5573 , 5581 , 5591 , 5623 , 5639 , &
5641 , 5647 , 5651 , 5653 , 5657 , 5659 , 5669 , 5683 , 5689 , 5693 , &
5701 , 5711 , 5717 , 5737 , 5741 , 5743 , 5749 , 5779 , 5783 , 5791 , &
5801 , 5807 , 5813 , 5821 , 5827 , 5839 , 5843 , 5849 , 5851 , 5857 , &
5861 , 5867 , 5869 , 5879 , 5881 , 5897 , 5903 , 5923 , 5927 , 5939 , &
5953 , 5981 , 5987 , 6007 , 6011 , 6029 , 6037 , 6043 , 6047 , 6053 , &
6067 , 6073 , 6079 , 6089 , 6091 , 6101 , 6113 , 6121 , 6131 , 6133 , &
! 801:900
6143 , 6151 , 6163 , 6173 , 6197 , 6199 , 6203 , 6211 , 6217 , 6221 , &
6229 , 6247 , 6257 , 6263 , 6269 , 6271 , 6277 , 6287 , 6299 , 6301 , &
6311 , 6317 , 6323 , 6329 , 6337 , 6343 , 6353 , 6359 , 6361 , 6367 , &
6373 , 6379 , 6389 , 6397 , 6421 , 6427 , 6449 , 6451 , 6469 , 6473 , &
6481 , 6491 , 6521 , 6529 , 6547 , 6551 , 6553 , 6563 , 6569 , 6571 , &
6577 , 6581 , 6599 , 6607 , 6619 , 6637 , 6653 , 6659 , 6661 , 6673 , &
6679 , 6689 , 6691 , 6701 , 6703 , 6709 , 6719 , 6733 , 6737 , 6761 , &
6763 , 6779 , 6781 , 6791 , 6793 , 6803 , 6823 , 6827 , 6829 , 6833 , &
6841 , 6857 , 6863 , 6869 , 6871 , 6883 , 6899 , 6907 , 6911 , 6917 , &
6947 , 6949 , 6959 , 6961 , 6967 , 6971 , 6977 , 6983 , 6991 , 6997 , &
! 901:1000
7001 , 7013 , 7019 , 7027 , 7039 , 7043 , 7057 , 7069 , 7079 , 7103 , &
7109 , 7121 , 7127 , 7129 , 7151 , 7159 , 7177 , 7187 , 7193 , 7207 , &
7211 , 7213 , 7219 , 7229 , 7237 , 7243 , 7247 , 7253 , 7283 , 7297 , &
7307 , 7309 , 7321 , 7331 , 7333 , 7349 , 7351 , 7369 , 7393 , 7411 , &
7417 , 7433 , 7451 , 7457 , 7459 , 7477 , 7481 , 7487 , 7489 , 7499 , &
7507 , 7517 , 7523 , 7529 , 7537 , 7541 , 7547 , 7549 , 7559 , 7561 , &
7573 , 7577 , 7583 , 7589 , 7591 , 7603 , 7607 , 7621 , 7639 , 7643 , &
7649 , 7669 , 7673 , 7681 , 7687 , 7691 , 7699 , 7703 , 7717 , 7723 , &
7727 , 7741 , 7753 , 7757 , 7759 , 7789 , 7793 , 7817 , 7823 , 7829 , &
7841 , 7853 , 7867 , 7873 , 7877 , 7879 , 7883 , 7901 , 7907 , 7919 , &
! 1001:1100
7927 , 7933 , 7937 , 7949 , 7951 , 7963 , 7993 , 8009 , 8011 , 8017 , &
8039 , 8053 , 8059 , 8069 , 8081 , 8087 , 8089 , 8093 , 8101 , 8111 , &
8117 , 8123 , 8147 , 8161 , 8167 , 8171 , 8179 , 8191 , 8209 , 8219 , &
8221 , 8231 , 8233 , 8237 , 8243 , 8263 , 8269 , 8273 , 8287 , 8291 , &
8293 , 8297 , 8311 , 8317 , 8329 , 8353 , 8363 , 8369 , 8377 , 8387 , &
8389 , 8419 , 8423 , 8429 , 8431 , 8443 , 8447 , 8461 , 8467 , 8501 , &
8513 , 8521 , 8527 , 8537 , 8539 , 8543 , 8563 , 8573 , 8581 , 8597 , &
8599 , 8609 , 8623 , 8627 , 8629 , 8641 , 8647 , 8663 , 8669 , 8677 , &
8681 , 8689 , 8693 , 8699 , 8707 , 8713 , 8719 , 8731 , 8737 , 8741 , &
8747 , 8753 , 8761 , 8779 , 8783 , 8803 , 8807 , 8819 , 8821 , 8831 , &
! 1101:1200
8837 , 8839 , 8849 , 8861 , 8863 , 8867 , 8887 , 8893 , 8923 , 8929 , &
8933 , 8941 , 8951 , 8963 , 8969 , 8971 , 8999 , 9001 , 9007 , 9011 , &
9013 , 9029 , 9041 , 9043 , 9049 , 9059 , 9067 , 9091 , 9103 , 9109 , &
9127 , 9133 , 9137 , 9151 , 9157 , 9161 , 9173 , 9181 , 9187 , 9199 , &
9203 , 9209 , 9221 , 9227 , 9239 , 9241 , 9257 , 9277 , 9281 , 9283 , &
9293 , 9311 , 9319 , 9323 , 9337 , 9341 , 9343 , 9349 , 9371 , 9377 , &
9391 , 9397 , 9403 , 9413 , 9419 , 9421 , 9431 , 9433 , 9437 , 9439 , &
9461 , 9463 , 9467 , 9473 , 9479 , 9491 , 9497 , 9511 , 9521 , 9533 , &
9539 , 9547 , 9551 , 9587 , 9601 , 9613 , 9619 , 9623 , 9629 , 9631 , &
9643 , 9649 , 9661 , 9677 , 9679 , 9689 , 9697 , 9719 , 9721 , 9733 , &
! 1201:1300
9739 , 9743 , 9749 , 9767 , 9769 , 9781 , 9787 , 9791 , 9803 , 9811 , &
9817 , 9829 , 9833 , 9839 , 9851 , 9857 , 9859 , 9871 , 9883 , 9887 , &
9901 , 9907 , 9923 , 9929 , 9931 , 9941 , 9949 , 9967 , 9973 , 10007 , &
10009 , 10037 , 10039 , 10061 , 10067 , 10069 , 10079 , 10091 , 10093 , 10099 , &
10103 , 10111 , 10133 , 10139 , 10141 , 10151 , 10159 , 10163 , 10169 , 10177 , &
10181 , 10193 , 10211 , 10223 , 10243 , 10247 , 10253 , 10259 , 10267 , 10271 , &
10273 , 10289 , 10301 , 10303 , 10313 , 10321 , 10331 , 10333 , 10337 , 10343 , &
10357 , 10369 , 10391 , 10399 , 10427 , 10429 , 10433 , 10453 , 10457 , 10459 , &
10463 , 10477 , 10487 , 10499 , 10501 , 10513 , 10529 , 10531 , 10559 , 10567 , &
10589 , 10597 , 10601 , 10607 , 10613 , 10627 , 10631 , 10639 , 10651 , 10657 , &
! 1301:1400
10663 , 10667 , 10687 , 10691 , 10709 , 10711 , 10723 , 10729 , 10733 , 10739 , &
10753 , 10771 , 10781 , 10789 , 10799 , 10831 , 10837 , 10847 , 10853 , 10859 , &
10861 , 10867 , 10883 , 10889 , 10891 , 10903 , 10909 , 19037 , 10939 , 10949 , &
10957 , 10973 , 10979 , 10987 , 10993 , 11003 , 11027 , 11047 , 11057 , 11059 , &
11069 , 11071 , 11083 , 11087 , 11093 , 11113 , 11117 , 11119 , 11131 , 11149 , &
11159 , 11161 , 11171 , 11173 , 11177 , 11197 , 11213 , 11239 , 11243 , 11251 , &
11257 , 11261 , 11273 , 11279 , 11287 , 11299 , 11311 , 11317 , 11321 , 11329 , &
11351 , 11353 , 11369 , 11383 , 11393 , 11399 , 11411 , 11423 , 11437 , 11443 , &
11447 , 11467 , 11471 , 11483 , 11489 , 11491 , 11497 , 11503 , 11519 , 11527 , &
11549 , 11551 , 11579 , 11587 , 11593 , 11597 , 11617 , 11621 , 11633 , 11657 , &
! 1401:1500
11677 , 11681 , 11689 , 11699 , 11701 , 11717 , 11719 , 11731 , 11743 , 11777 , &
11779 , 11783 , 11789 , 11801 , 11807 , 11813 , 11821 , 11827 , 11831 , 11833 , &
11839 , 11863 , 11867 , 11887 , 11897 , 11903 , 11909 , 11923 , 11927 , 11933 , &
11939 , 11941 , 11953 , 11959 , 11969 , 11971 , 11981 , 11987 , 12007 , 12011 , &
12037 , 12041 , 12043 , 12049 , 12071 , 12073 , 12097 , 12101 , 12107 , 12109 , &
12113 , 12119 , 12143 , 12149 , 12157 , 12161 , 12163 , 12197 , 12203 , 12211 , &
12227 , 12239 , 12241 , 12251 , 12253 , 12263 , 12269 , 12277 , 12281 , 12289 , &
12301 , 12323 , 12329 , 12343 , 12347 , 12373 , 12377 , 12379 , 12391 , 12401 , &
12409 , 12413 , 12421 , 12433 , 12437 , 12451 , 12457 , 12473 , 12479 , 12487 , &
12491 , 12497 , 12503 , 12511 , 12517 , 12527 , 12539 , 12541 , 12547 , 12553 ] , pInt )
if ( n < size ( npvec ) ) then
prime = npvec ( n )
else
call IO_error ( error_ID = 406_pInt )
end if
end function prime
2017-08-12 09:33:40 +05:30
end subroutine halton_memory
!--------------------------------------------------------------------------------------------------
!> @brief sets the dimension for a Halton sequence
!> @author John Burkardt
!--------------------------------------------------------------------------------------------------
2017-08-12 10:05:44 +05:30
subroutine halton_ndim_set ( ndim )
2017-08-12 09:33:40 +05:30
implicit none
integer ( pInt ) , intent ( in ) :: ndim !< dimension of the Halton vectors
integer ( pInt ) :: value_halton ( 1 )
value_halton ( 1 ) = ndim
call halton_memory ( 'SET' , 'NDIM' , 1_pInt , value_halton )
end subroutine halton_ndim_set
!--------------------------------------------------------------------------------------------------
!> @brief sets the seed for the Halton sequence.
!> @details Calling HALTON repeatedly returns the elements of the Halton sequence in order,
!> @details starting with element number 1.
!> @details An internal counter, called SEED, keeps track of the next element to return. Each time
!> @details is computed, and then SEED is incremented by 1.
!> @details To restart the Halton sequence, it is only necessary to reset SEED to 1. It might also
!> @details be desirable to reset SEED to some other value. This routine allows the user to specify
!> @details any value of SEED.
!> @details The default value of SEED is 1, which restarts the Halton sequence.
!> @author John Burkardt
!--------------------------------------------------------------------------------------------------
subroutine halton_seed_set ( seed )
implicit none
integer ( pInt ) , parameter :: NDIM = 1_pInt
integer ( pInt ) , intent ( in ) :: seed !< seed for the Halton sequence.
integer ( pInt ) :: value_halton ( ndim )
value_halton ( 1 ) = seed
call halton_memory ( 'SET' , 'SEED' , NDIM , value_halton )
end subroutine halton_seed_set
2014-08-21 14:03:55 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief factorial
!--------------------------------------------------------------------------------------------------
integer ( pInt ) pure function math_factorial ( n )
implicit none
integer ( pInt ) , intent ( in ) :: n
integer ( pInt ) :: i
2014-09-16 15:20:37 +05:30
math_factorial = product ( [ ( i , i = 1 , n ) ] )
2014-08-21 14:03:55 +05:30
end function math_factorial
!--------------------------------------------------------------------------------------------------
!> @brief binomial coefficient
!--------------------------------------------------------------------------------------------------
integer ( pInt ) pure function math_binomial ( n , k )
implicit none
integer ( pInt ) , intent ( in ) :: n , k
integer ( pInt ) :: i , j
j = min ( k , n - k )
2014-09-16 15:20:37 +05:30
math_binomial = product ( [ ( i , i = n , n - j + 1 , - 1 ) ] ) / math_factorial ( j )
2014-08-21 14:03:55 +05:30
end function math_binomial
!--------------------------------------------------------------------------------------------------
!> @brief multinomial coefficient
!--------------------------------------------------------------------------------------------------
integer ( pInt ) pure function math_multinomial ( alpha )
implicit none
2015-03-29 18:24:13 +05:30
integer ( pInt ) , intent ( in ) , dimension ( : ) :: alpha
2014-08-21 14:03:55 +05:30
integer ( pInt ) :: i
2015-03-29 18:24:13 +05:30
math_multinomial = 1_pInt
2014-08-21 14:03:55 +05:30
do i = 1 , size ( alpha )
math_multinomial = math_multinomial * math_binomial ( sum ( alpha ( 1 : i ) ) , alpha ( i ) )
enddo
end function math_multinomial
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief volume of tetrahedron given by four vertices
!--------------------------------------------------------------------------------------------------
2013-01-31 21:58:08 +05:30
real ( pReal ) pure function math_volTetrahedron ( v1 , v2 , v3 , v4 )
2009-01-20 00:40:58 +05:30
implicit none
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: v1 , v2 , v3 , v4
real ( pReal ) , dimension ( 3 , 3 ) :: m
2007-03-20 19:25:22 +05:30
2011-12-01 17:31:13 +05:30
m ( 1 : 3 , 1 ) = v1 - v2
m ( 1 : 3 , 2 ) = v2 - v3
m ( 1 : 3 , 3 ) = v3 - v4
2009-01-20 00:40:58 +05:30
2012-08-25 17:16:36 +05:30
math_volTetrahedron = math_det33 ( m ) / 6.0_pReal
2009-01-20 00:40:58 +05:30
2012-03-09 01:55:28 +05:30
end function math_volTetrahedron
2009-01-20 00:40:58 +05:30
2011-12-01 17:31:13 +05:30
2013-04-09 23:37:30 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief area of triangle given by three vertices
!--------------------------------------------------------------------------------------------------
real ( pReal ) pure function math_areaTriangle ( v1 , v2 , v3 )
implicit none
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: v1 , v2 , v3
2016-01-10 19:04:26 +05:30
math_areaTriangle = 0.5_pReal * norm2 ( math_crossproduct ( v1 - v2 , v1 - v3 ) )
2013-04-09 23:37:30 +05:30
end function math_areaTriangle
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief rotate 33 tensor forward
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
pure function math_rotate_forward33 ( tensor , rot_tensor )
2011-10-24 23:56:34 +05:30
implicit none
2012-01-26 19:20:00 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: math_rotate_forward33
2011-10-24 23:56:34 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: tensor , rot_tensor
2012-08-25 17:16:36 +05:30
2012-01-26 19:20:00 +05:30
math_rotate_forward33 = math_mul33x33 ( rot_tensor , &
math_mul33x33 ( tensor , math_transpose33 ( rot_tensor ) ) )
2012-08-25 17:16:36 +05:30
2012-03-09 01:55:28 +05:30
end function math_rotate_forward33
2011-10-24 23:56:34 +05:30
2011-12-01 17:31:13 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief rotate 33 tensor backward
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
pure function math_rotate_backward33 ( tensor , rot_tensor )
2011-10-24 23:56:34 +05:30
implicit none
2012-01-26 19:20:00 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: math_rotate_backward33
2011-10-24 23:56:34 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: tensor , rot_tensor
2012-08-25 17:16:36 +05:30
2012-01-26 19:20:00 +05:30
math_rotate_backward33 = math_mul33x33 ( math_transpose33 ( rot_tensor ) , &
2011-10-24 23:56:34 +05:30
math_mul33x33 ( tensor , rot_tensor ) )
2012-08-25 17:16:36 +05:30
2012-03-09 01:55:28 +05:30
end function math_rotate_backward33
2011-10-24 23:56:34 +05:30
2011-12-01 17:31:13 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief rotate 3333 tensor C'_ijkl=g_im*g_jn*g_ko*g_lp*C_mnop
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
pure function math_rotate_forward3333 ( tensor , rot_tensor )
2011-10-24 23:56:34 +05:30
implicit none
2012-01-26 19:20:00 +05:30
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: math_rotate_forward3333
2011-10-24 23:56:34 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: rot_tensor
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) , intent ( in ) :: tensor
2011-10-25 19:08:24 +05:30
integer ( pInt ) :: i , j , k , l , m , n , o , p
2012-08-25 17:16:36 +05:30
2012-01-26 19:20:00 +05:30
math_rotate_forward3333 = 0.0_pReal
2011-10-25 19:08:24 +05:30
2011-12-01 17:31:13 +05:30
do i = 1_pInt , 3_pInt ; do j = 1_pInt , 3_pInt ; do k = 1_pInt , 3_pInt ; do l = 1_pInt , 3_pInt
do m = 1_pInt , 3_pInt ; do n = 1_pInt , 3_pInt ; do o = 1_pInt , 3_pInt ; do p = 1_pInt , 3_pInt
2014-01-31 19:58:20 +05:30
math_rotate_forward3333 ( i , j , k , l ) = math_rotate_forward3333 ( i , j , k , l ) &
+ rot_tensor ( m , i ) * rot_tensor ( n , j ) &
* rot_tensor ( o , k ) * rot_tensor ( p , l ) * tensor ( m , n , o , p )
2011-10-25 19:08:24 +05:30
enddo ; enddo ; enddo ; enddo ; enddo ; enddo ; enddo ; enddo
2012-08-25 17:16:36 +05:30
2012-03-09 01:55:28 +05:30
end function math_rotate_forward3333
2009-01-20 00:40:58 +05:30
2011-12-01 17:31:13 +05:30
2016-01-10 19:04:26 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief limits a scalar value to a certain range (either one or two sided)
! Will return NaN if left > right
!--------------------------------------------------------------------------------------------------
real ( pReal ) pure function math_limit ( a , left , right )
2017-05-04 04:02:44 +05:30
use , intrinsic :: &
IEEE_arithmetic
2016-01-10 19:04:26 +05:30
implicit none
real ( pReal ) , intent ( in ) :: a
real ( pReal ) , intent ( in ) , optional :: left , right
math_limit = min ( &
max ( merge ( left , - huge ( a ) , present ( left ) ) , a ) , &
merge ( right , huge ( a ) , present ( right ) ) &
)
2017-05-04 04:02:44 +05:30
if ( present ( left ) . and . present ( right ) ) &
math_limit = merge ( IEEE_value ( 1.0_pReal , IEEE_quiet_NaN ) , math_limit , left > right )
2016-01-10 19:04:26 +05:30
end function math_limit
2016-02-26 21:05:55 +05:30
end module math