2012-07-19 00:09:59 +05:30
!--------------------------------------------------------------------------------------------------
2013-03-07 01:04:30 +05:30
! $Id$
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
!> @brief Mathematical library, including random number generation and tensor represenations
!--------------------------------------------------------------------------------------------------
2012-08-25 17:16:36 +05:30
module math
2012-01-13 21:48:16 +05:30
use , intrinsic :: iso_c_binding
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
2013-03-07 01:04:30 +05:30
real ( pReal ) , parameter , public :: PI = 3.14159265358979323846264338327950288419716939937510_pReal !< ratio of a circle's circumference to its diameter
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
complex ( pReal ) , parameter , public :: TWOPIIMG = ( 0.0_pReal , 2.0_pReal ) * PI !< 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
2012-06-15 21:40:21 +05:30
#ifdef Spectral
2012-02-01 00:48:55 +05:30
include 'fftw3.f03'
2012-06-15 21:40:21 +05:30
#endif
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 , &
math_vectorproduct , &
math_tensorproduct , &
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 , &
math_deviatoric33 , &
math_equivStrain33 , &
2015-01-13 15:13:05 +05:30
math_equivStress33 , &
2013-01-31 21:58:08 +05:30
math_trace33 , &
2014-10-14 06:40:42 +05:30
math_j3_33 , &
2013-01-31 21:58:08 +05:30
math_det33 , &
math_norm33 , &
math_norm3 , &
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_qNorm , &
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 , &
math_spectralDecompositionSym33 , &
math_spectralDecomposition , &
math_pDecomposition , &
math_hi , &
math_eigenvalues33 , &
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 , &
math_rotate_forward3333
2012-08-27 13:34:47 +05:30
#ifdef Spectral
2013-01-31 21:58:08 +05:30
public :: &
2013-06-06 14:55:27 +05:30
fftw_set_timelimit , &
fftw_plan_dft_3d , &
fftw_plan_many_dft_r2c , &
fftw_plan_many_dft_c2r , &
fftw_plan_with_nthreads , &
fftw_init_threads , &
fftw_alloc_complex , &
fftw_execute_dft , &
fftw_execute_dft_r2c , &
fftw_execute_dft_c2r , &
fftw_destroy_plan , &
2013-01-31 21:58:08 +05:30
math_curlFFT , &
2013-03-07 13:28:42 +05:30
math_gradFFT , &
2013-01-31 21:58:08 +05:30
math_divergenceFFT , &
math_divergenceFDM , &
math_tensorAvg , &
math_logstrainSpat , &
math_logstrainMat , &
math_cauchy , &
2013-02-19 20:26:26 +05:30
math_periodicNearestNeighbor , &
2013-06-07 12:57:12 +05:30
math_nearestNeighbor , &
2013-02-19 20:26:26 +05:30
math_periodicNearestNeighborDistances
2013-02-11 15:14:17 +05:30
#endif
2013-01-31 21:58:08 +05:30
private :: &
math_partition , &
halton , &
halton_memory , &
halton_ndim_set , &
halton_seed_set , &
i_to_halton , &
prime
2013-02-11 15:14:17 +05:30
external :: &
dsyev , &
dgetrf , &
dgetri
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 : &
worldrank , &
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
2014-10-10 21:51:10 +05:30
mainProcess : if ( worldrank == 0 ) then
write ( 6 , '(/,a)' ) ' <<<+- math init -+>>>'
write ( 6 , '(a)' ) ' $Id$'
write ( 6 , '(a15,a)' ) ' Current time: ' , IO_timeStamp ( )
2012-02-01 00:48:55 +05:30
#include "compilation_info.f90"
2014-10-10 21:51:10 +05:30
endif mainProcess
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
2014-10-13 14:14:49 +05:30
mainProcess2 : if ( worldrank == 0 ) then
write ( 6 , * ) 'size of random seed: ' , randSize
do i = 1 , randSize
write ( 6 , * ) 'value of random seed: ' , i , randInit ( i )
enddo
write ( 6 , '(a,4(/,26x,f17.14),/)' ) ' start of random sequence: ' , randTest
endif mainProcess2
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
ipivot = math_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
2013-02-01 12:56:21 +05:30
end subroutine math_qsort
2007-04-03 13:47:58 +05:30
2011-12-01 17:31:13 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Partitioning required for quicksort
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
integer ( pInt ) function math_partition ( 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 ) :: d , i , j , k , x , tmp
2007-04-03 13:47:58 +05:30
2012-02-16 00:28:38 +05:30
d = int ( size ( a , 1_pInt ) , pInt ) ! number of linked data
2008-02-15 18:12:27 +05:30
! set the starting and ending points, and the pivot point
i = istart
2007-04-04 14:19:48 +05:30
j = iend
2007-04-03 13:47:58 +05:30
x = a ( 1 , istart )
do
! find the first element on the right side less than or equal to the pivot point
2011-12-01 17:31:13 +05:30
do j = j , istart , - 1_pInt
2007-04-03 13:47:58 +05:30
if ( a ( 1 , j ) < = x ) exit
enddo
! find the first element on the left side greater than the pivot point
do i = i , iend
if ( a ( 1 , i ) > x ) exit
enddo
2011-12-01 17:31:13 +05:30
if ( i < j ) then ! if the indexes do not cross, exchange values
do k = 1_pInt , d
2007-04-03 13:47:58 +05:30
tmp = a ( k , i )
a ( k , i ) = a ( k , j )
a ( k , j ) = tmp
enddo
2007-12-14 19:06:04 +05:30
else ! if they do cross, exchange left value with pivot and return with the partition index
2011-12-01 17:31:13 +05:30
do k = 1_pInt , d
2007-04-03 13:47:58 +05:30
tmp = a ( k , istart )
a ( k , istart ) = a ( k , j )
a ( k , j ) = tmp
enddo
2007-04-04 14:19:48 +05:30
math_partition = j
2007-04-03 13:47:58 +05:30
return
endif
enddo
2012-03-09 01:55:28 +05:30
end function math_partition
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
!--------------------------------------------------------------------------------------------------
2013-01-31 21:58:08 +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
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
2013-01-31 21:58:08 +05:30
if ( i / = j ) then
math_delta = 0.0_pReal
else
math_delta = 1.0_pReal
endif
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
!--------------------------------------------------------------------------------------------------
!> @brief vector product a x b
!--------------------------------------------------------------------------------------------------
pure function math_vectorproduct ( A , B )
2009-01-20 00:40:58 +05:30
implicit none
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: A , B
real ( pReal ) , dimension ( 3 ) :: math_vectorproduct
math_vectorproduct ( 1 ) = A ( 2 ) * B ( 3 ) - A ( 3 ) * B ( 2 )
math_vectorproduct ( 2 ) = A ( 3 ) * B ( 1 ) - A ( 1 ) * B ( 3 )
math_vectorproduct ( 3 ) = A ( 1 ) * B ( 2 ) - A ( 2 ) * B ( 1 )
2012-03-09 01:55:28 +05:30
end function math_vectorproduct
2009-01-20 00:40:58 +05:30
2009-03-05 20:07:59 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief tensor product a \otimes b
!--------------------------------------------------------------------------------------------------
pure function math_tensorproduct ( A , B )
2009-03-17 20:43:17 +05:30
implicit none
real ( pReal ) , dimension ( 3 , 3 ) :: math_tensorproduct
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
2011-12-01 17:31:13 +05:30
forall ( i = 1_pInt : 3_pInt , j = 1_pInt : 3_pInt ) math_tensorproduct ( i , j ) = A ( i ) * B ( j )
2009-03-17 20:43:17 +05:30
2012-03-09 01:55:28 +05:30
end function math_tensorproduct
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
do i = 1_pInt , 3_pInt
do j = 1_pInt , 3_pInt
do k = 1_pInt , 3_pInt
do 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 ) )
enddo ; enddo ; enddo ; enddo
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
2011-12-01 17:31:13 +05:30
forall ( i = 1_pInt : 3_pInt , j = 1_pInt : 3_pInt ) math_mul33x33 ( i , j ) = &
2009-01-20 00:40:58 +05:30
A ( i , 1 ) * B ( 1 , j ) + A ( i , 2 ) * B ( 2 , j ) + A ( i , 3 ) * B ( 3 , j )
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
order = 5
if ( present ( n ) ) order = n
2013-05-08 21:22:29 +05:30
B = math_identity2nd ( 3 ) ! init
invfac = 1.0_pReal ! 0!
math_exp33 = B ! A^0 = eye2
2012-10-12 23:24:20 +05:30
do i = 1_pInt , n
2013-05-08 21:22:29 +05:30
invfac = invfac / real ( i ) ! 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 )
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
2012-01-26 19:20:00 +05:30
math_inv33 = 0.0_pReal
2009-03-31 13:01:38 +05:30
2011-12-01 17:31:13 +05:30
DetA = A ( 1 , 1 ) * ( A ( 2 , 2 ) * A ( 3 , 3 ) - A ( 2 , 3 ) * A ( 3 , 2 ) ) &
- A ( 1 , 2 ) * ( A ( 2 , 1 ) * A ( 3 , 3 ) - A ( 2 , 3 ) * A ( 3 , 1 ) ) &
+ A ( 1 , 3 ) * ( A ( 2 , 1 ) * A ( 3 , 2 ) - A ( 2 , 2 ) * A ( 3 , 1 ) )
2009-03-31 13:01:38 +05:30
2015-04-11 13:55:23 +05:30
if ( abs ( DetA ) > tiny ( DetA ) ) then
2012-01-26 19:20:00 +05:30
math_inv33 ( 1 , 1 ) = ( A ( 2 , 2 ) * A ( 3 , 3 ) - A ( 2 , 3 ) * A ( 3 , 2 ) ) / DetA
math_inv33 ( 2 , 1 ) = ( - A ( 2 , 1 ) * A ( 3 , 3 ) + A ( 2 , 3 ) * A ( 3 , 1 ) ) / DetA
math_inv33 ( 3 , 1 ) = ( A ( 2 , 1 ) * A ( 3 , 2 ) - A ( 2 , 2 ) * A ( 3 , 1 ) ) / DetA
2009-03-31 13:01:38 +05:30
2012-01-26 19:20:00 +05:30
math_inv33 ( 1 , 2 ) = ( - A ( 1 , 2 ) * A ( 3 , 3 ) + A ( 1 , 3 ) * A ( 3 , 2 ) ) / DetA
math_inv33 ( 2 , 2 ) = ( A ( 1 , 1 ) * A ( 3 , 3 ) - A ( 1 , 3 ) * A ( 3 , 1 ) ) / DetA
math_inv33 ( 3 , 2 ) = ( - A ( 1 , 1 ) * A ( 3 , 2 ) + A ( 1 , 2 ) * A ( 3 , 1 ) ) / DetA
2009-03-31 13:01:38 +05:30
2012-01-26 19:20:00 +05:30
math_inv33 ( 1 , 3 ) = ( A ( 1 , 2 ) * A ( 2 , 3 ) - A ( 1 , 3 ) * A ( 2 , 2 ) ) / DetA
math_inv33 ( 2 , 3 ) = ( - A ( 1 , 1 ) * A ( 2 , 3 ) + A ( 1 , 3 ) * A ( 2 , 1 ) ) / DetA
math_inv33 ( 3 , 3 ) = ( A ( 1 , 1 ) * A ( 2 , 2 ) - A ( 1 , 2 ) * A ( 2 , 1 ) ) / DetA
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)
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
pure subroutine math_invert33 ( A , InvA , DetA , error )
2008-02-15 18:12:27 +05:30
2012-01-26 19:20:00 +05:30
! Bestimmung der Determinanten und Inversen einer 33-Matrix
2008-02-15 18:12:27 +05:30
! A = Matrix A
! InvA = Inverse of A
! DetA = Determinant of A
2007-04-11 15:34:22 +05:30
! error = logical
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
2011-12-01 17:31:13 +05:30
DetA = A ( 1 , 1 ) * ( A ( 2 , 2 ) * A ( 3 , 3 ) - A ( 2 , 3 ) * A ( 3 , 2 ) ) &
- A ( 1 , 2 ) * ( A ( 2 , 1 ) * A ( 3 , 3 ) - A ( 2 , 3 ) * A ( 3 , 1 ) ) &
+ A ( 1 , 3 ) * ( A ( 2 , 1 ) * A ( 3 , 2 ) - A ( 2 , 2 ) * A ( 3 , 1 ) )
2008-02-15 18:12:27 +05:30
2015-04-11 13:55:23 +05:30
if ( abs ( DetA ) < = tiny ( DetA ) ) then
2008-02-15 18:12:27 +05:30
error = . true .
2007-04-11 15:34:22 +05:30
else
2011-12-01 17:31:13 +05:30
InvA ( 1 , 1 ) = ( A ( 2 , 2 ) * A ( 3 , 3 ) - A ( 2 , 3 ) * A ( 3 , 2 ) ) / DetA
InvA ( 2 , 1 ) = ( - A ( 2 , 1 ) * A ( 3 , 3 ) + A ( 2 , 3 ) * A ( 3 , 1 ) ) / DetA
InvA ( 3 , 1 ) = ( A ( 2 , 1 ) * A ( 3 , 2 ) - A ( 2 , 2 ) * A ( 3 , 1 ) ) / DetA
2008-02-15 18:12:27 +05:30
2011-12-01 17:31:13 +05:30
InvA ( 1 , 2 ) = ( - A ( 1 , 2 ) * A ( 3 , 3 ) + A ( 1 , 3 ) * A ( 3 , 2 ) ) / DetA
InvA ( 2 , 2 ) = ( A ( 1 , 1 ) * A ( 3 , 3 ) - A ( 1 , 3 ) * A ( 3 , 1 ) ) / DetA
InvA ( 3 , 2 ) = ( - A ( 1 , 1 ) * A ( 3 , 2 ) + A ( 1 , 2 ) * A ( 3 , 1 ) ) / DetA
2008-02-15 18:12:27 +05:30
2011-12-01 17:31:13 +05:30
InvA ( 1 , 3 ) = ( A ( 1 , 2 ) * A ( 2 , 3 ) - A ( 1 , 3 ) * A ( 2 , 2 ) ) / DetA
InvA ( 2 , 3 ) = ( - A ( 1 , 1 ) * A ( 2 , 3 ) + A ( 1 , 3 ) * A ( 2 , 1 ) ) / DetA
InvA ( 3 , 3 ) = ( A ( 1 , 1 ) * A ( 2 , 2 ) - A ( 1 , 2 ) * A ( 2 , 1 ) ) / DetA
2012-08-25 17:16:36 +05:30
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
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
#if(FLOAT==8)
call dgetrf ( 6 , 6 , temp66_real , 6 , ipiv6 , ierr )
call dgetri ( 6 , temp66_real , 6 , ipiv6 , work6 , 6 , ierr )
#elif(FLOAT==4)
call sgetrf ( 6 , 6 , temp66_real , 6 , ipiv6 , ierr )
call sgetri ( 6 , temp66_real , 6 , ipiv6 , work6 , 6 , ierr )
#endif
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
2012-08-28 22:29:45 +05:30
invA = A
#if(FLOAT==8)
call dgetrf ( myDim , myDim , invA , myDim , ipiv , ierr )
call dgetri ( myDim , InvA , myDim , ipiv , work , myDim , ierr )
#elif(FLOAT==4)
call sgetrf ( myDim , myDim , invA , myDim , ipiv , ierr )
call sgetri ( myDim , InvA , myDim , ipiv , work , myDim , ierr )
#endif
if ( ierr == 0_pInt ) then
error = . false .
else
error = . true .
endif
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
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
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
integer ( pInt ) :: i , j
2008-02-15 18:12:27 +05:30
2012-01-26 19:20:00 +05:30
forall ( i = 1_pInt : 3_pInt , j = 1_pInt : 3_pInt ) math_symmetric33 ( i , j ) = 0.5_pReal * ( m ( i , j ) + m ( j , i ) )
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
2011-12-01 17:31:13 +05:30
integer ( pInt ) :: i , j
2009-01-20 00:40:58 +05:30
real ( pReal ) , dimension ( 6 , 6 ) , intent ( in ) :: m
2012-01-26 19:20:00 +05:30
real ( pReal ) , dimension ( 6 , 6 ) :: math_symmetric66
2012-08-25 17:16:36 +05:30
2012-01-26 19:20:00 +05:30
forall ( i = 1_pInt : 6_pInt , j = 1_pInt : 6_pInt ) math_symmetric66 ( i , j ) = 0.5_pReal * ( m ( i , j ) + m ( j , i ) )
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
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_skew33 ( i , j ) = m ( i , j ) - 0.5_pReal * ( m ( i , j ) + m ( j , i ) )
2012-03-09 01:55:28 +05:30
end function math_skew33
2008-02-15 18:12:27 +05:30
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
integer ( pInt ) :: i
real ( pReal ) :: hydrostatic
2012-01-26 19:20:00 +05:30
2012-02-09 21:28:15 +05:30
math_deviatoric33 = m
2012-10-12 23:24:20 +05:30
hydrostatic = ( m ( 1 , 1 ) + m ( 2 , 2 ) + m ( 3 , 3 ) ) / 3.0_pReal
2012-02-09 21:28:15 +05:30
forall ( i = 1_pInt : 3_pInt ) math_deviatoric33 ( i , i ) = m ( i , i ) - hydrostatic
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
!--------------------------------------------------------------------------------------------------
!> @brief equivalent scalar quantity of a full strain tensor
!--------------------------------------------------------------------------------------------------
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
2011-12-01 17:31:13 +05:30
real ( pReal ) :: math_equivStrain33 , e11 , e22 , e33 , s12 , s23 , s31
2010-03-24 18:50:12 +05:30
e11 = ( 2.0_pReal * m ( 1 , 1 ) - m ( 2 , 2 ) - m ( 3 , 3 ) ) / 3.0_pReal
e22 = ( 2.0_pReal * m ( 2 , 2 ) - m ( 3 , 3 ) - m ( 1 , 1 ) ) / 3.0_pReal
e33 = ( 2.0_pReal * m ( 3 , 3 ) - m ( 1 , 1 ) - m ( 2 , 2 ) ) / 3.0_pReal
s12 = 2.0_pReal * m ( 1 , 2 )
s23 = 2.0_pReal * m ( 2 , 3 )
s31 = 2.0_pReal * m ( 3 , 1 )
math_equivStrain33 = 2.0_pReal * ( 1.50_pReal * ( e11 ** 2.0_pReal + e22 ** 2.0_pReal + e33 ** 2.0_pReal ) + &
0.75_pReal * ( s12 ** 2.0_pReal + s23 ** 2.0_pReal + s31 ** 2.0_pReal ) ) ** ( 0.5_pReal ) / 3.0_pReal
2012-03-09 01:55:28 +05:30
end function math_equivStrain33
2010-03-24 18:50:12 +05:30
2015-01-13 15:13:05 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief von Mises equivalent of a full stress tensor
!--------------------------------------------------------------------------------------------------
pure function math_equivStress33 ( m )
implicit none
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m
real ( pReal ) :: math_equivStress33
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 * ( m ( 1 , 2 ) ** 2.0_pReal + m ( 3 , 1 ) ** 2.0_pReal + m ( 2 , 3 ) ** 2.0_pReal ) ) ** ( 0.5_pReal ) ) / sqrt ( 2.0_pReal )
end function math_equivStress33
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
2014-10-14 06:40:42 +05:30
!--------------------------------------------------------------------------------------------------
2014-10-14 09:00:59 +05:30
!> @brief invariant 3 of a 33 matrix
2014-10-14 06:40:42 +05:30
!--------------------------------------------------------------------------------------------------
real ( pReal ) pure function math_j3_33 ( m )
implicit none
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m
2014-10-14 09:00:59 +05:30
math_j3_33 = sum ( m ** 3.0_pReal ) / 3.0_pReal
2014-10-14 06:40:42 +05:30
end function math_j3_33
2012-10-12 23:24:20 +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
2007-03-21 15:50:25 +05:30
2012-01-26 19:20:00 +05:30
math_det33 = m ( 1 , 1 ) * ( m ( 2 , 2 ) * m ( 3 , 3 ) - m ( 2 , 3 ) * m ( 3 , 2 ) ) &
2013-01-31 21:58:08 +05:30
- 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
!--------------------------------------------------------------------------------------------------
!> @brief norm of a 33 matrix
!--------------------------------------------------------------------------------------------------
2013-01-31 21:58:08 +05:30
real ( pReal ) pure function math_norm33 ( m )
2011-07-29 21:27:39 +05:30
implicit none
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m
2013-01-31 21:58:08 +05:30
2011-07-29 21:27:39 +05:30
math_norm33 = sqrt ( sum ( m ** 2.0_pReal ) )
2012-03-09 01:55:28 +05:30
end function
2011-07-29 21:27:39 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2012-10-12 23:24:20 +05:30
!> @brief euclidian norm of a 3 vector
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2013-01-31 21:58:08 +05:30
real ( pReal ) pure function math_norm3 ( v )
2007-03-21 15:50:25 +05:30
2009-08-11 22:01:57 +05:30
implicit none
2010-03-18 17:53:17 +05:30
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: v
2009-08-11 22:01:57 +05:30
2013-01-08 16:39:20 +05:30
math_norm3 = sqrt ( sum ( v ** 2.0_pReal ) )
2012-08-25 17:16:36 +05:30
2012-03-09 01:55:28 +05:30
end function math_norm3
2009-08-11 22:01:57 +05:30
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 ) , &
2010-09-22 17:34:43 +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
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
!--------------------------------------------------------------------------------------------------
function math_qRand ( )
implicit none
real ( pReal ) , dimension ( 4 ) :: math_qRand
real ( pReal ) , dimension ( 3 ) :: rnd
call halton ( 3_pInt , rnd )
2013-05-08 21:40:21 +05:30
math_qRand ( 1 ) = cos ( 2.0_pReal * PI * rnd ( 1 ) ) * sqrt ( rnd ( 3 ) )
math_qRand ( 2 ) = sin ( 2.0_pReal * PI * rnd ( 2 ) ) * sqrt ( 1.0_pReal - rnd ( 3 ) )
math_qRand ( 3 ) = cos ( 2.0_pReal * PI * rnd ( 2 ) ) * sqrt ( 1.0_pReal - rnd ( 3 ) )
math_qRand ( 4 ) = 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
math_qMul ( 1 ) = A ( 1 ) * B ( 1 ) - A ( 2 ) * B ( 2 ) - A ( 3 ) * B ( 3 ) - A ( 4 ) * B ( 4 )
math_qMul ( 2 ) = A ( 1 ) * B ( 2 ) + A ( 2 ) * B ( 1 ) + A ( 3 ) * B ( 4 ) - A ( 4 ) * B ( 3 )
math_qMul ( 3 ) = A ( 1 ) * B ( 3 ) - A ( 2 ) * B ( 4 ) + A ( 3 ) * B ( 1 ) + A ( 4 ) * B ( 2 )
math_qMul ( 4 ) = A ( 1 ) * B ( 4 ) + A ( 2 ) * B ( 3 ) - A ( 3 ) * B ( 2 ) + A ( 4 ) * B ( 1 )
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
math_qConj ( 1 ) = Q ( 1 )
math_qConj ( 2 : 4 ) = - Q ( 2 : 4 )
end function math_qConj
!--------------------------------------------------------------------------------------------------
!> @brief quaternion norm
!--------------------------------------------------------------------------------------------------
real ( pReal ) pure function math_qNorm ( Q )
implicit none
real ( pReal ) , dimension ( 4 ) , intent ( in ) :: Q
math_qNorm = sqrt ( max ( 0.0_pReal , sum ( Q * Q ) ) )
end function math_qNorm
!--------------------------------------------------------------------------------------------------
!> @brief quaternion inversion
!--------------------------------------------------------------------------------------------------
pure function math_qInv ( Q )
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 )
2015-04-11 13:55:23 +05:30
if ( abs ( squareNorm ) > tiny ( squareNorm ) ) &
2013-01-31 21:58:08 +05:30
math_qInv = math_qConj ( Q ) / squareNorm
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
2012-03-09 01:55:28 +05:30
real ( pReal ) :: sqhkl , squvw , sqhk , myVal
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 ) )
sqhk = sqrt ( R ( 1 , 3 ) * R ( 1 , 3 ) + R ( 2 , 3 ) * R ( 2 , 3 ) )
! calculate PHI
2012-03-09 01:55:28 +05:30
myVal = R ( 3 , 3 ) / sqhkl
2012-08-25 17:16:36 +05:30
2014-05-15 14:22:16 +05:30
myVal = min ( myVal , 1.0_pReal )
myVal = max ( myVal , - 1.0_pReal )
2012-08-25 17:16:36 +05:30
2012-03-09 01:55:28 +05:30
math_RtoEuler ( 2 ) = acos ( myVal )
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
2007-03-29 21:02:52 +05:30
! calculate phi2
math_RtoEuler ( 3 ) = 0.0_pReal
! calculate phi1
2012-03-09 01:55:28 +05:30
myVal = R ( 1 , 1 ) / squvw
2014-05-15 14:22:16 +05:30
myVal = min ( myVal , 1.0_pReal )
myVal = max ( myVal , - 1.0_pReal )
2012-08-25 17:16:36 +05:30
2012-03-09 01:55:28 +05:30
math_RtoEuler ( 1 ) = acos ( myVal )
2007-03-29 21:02:52 +05:30
if ( R ( 2 , 1 ) > 0.0_pReal ) math_RtoEuler ( 1 ) = 2.0_pReal * pi - math_RtoEuler ( 1 )
else
! calculate phi2
2012-03-09 01:55:28 +05:30
myVal = R ( 2 , 3 ) / sqhk
2014-05-15 14:22:16 +05:30
myVal = min ( myVal , 1.0_pReal )
myVal = max ( myVal , - 1.0_pReal )
2012-08-25 17:16:36 +05:30
2012-03-09 01:55:28 +05:30
math_RtoEuler ( 3 ) = acos ( myVal )
2007-03-29 21:02:52 +05:30
if ( R ( 1 , 3 ) < 0.0 ) math_RtoEuler ( 3 ) = 2.0_pReal * pi - math_RtoEuler ( 3 )
! calculate phi1
2012-03-09 01:55:28 +05:30
myVal = - R ( 3 , 2 ) / sin ( math_RtoEuler ( 2 ) )
2014-05-15 14:22:16 +05:30
myVal = min ( myVal , 1.0_pReal )
myVal = max ( myVal , - 1.0_pReal )
2012-08-25 17:16:36 +05:30
2012-03-09 01:55:28 +05:30
math_RtoEuler ( 1 ) = acos ( myVal )
2007-03-29 21:02:52 +05:30
if ( R ( 3 , 1 ) < 0.0 ) math_RtoEuler ( 1 ) = 2.0_pReal * pi - math_RtoEuler ( 1 )
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
2013-06-06 00:40:37 +05:30
absQ ( 1 ) = 1.0_pReal + R ( 1 , 1 ) + R ( 2 , 2 ) + R ( 3 , 3 )
absQ ( 2 ) = 1.0_pReal + R ( 1 , 1 ) - R ( 2 , 2 ) - R ( 3 , 3 )
absQ ( 3 ) = 1.0_pReal - R ( 1 , 1 ) + R ( 2 , 2 ) - R ( 3 , 3 )
absQ ( 4 ) = 1.0_pReal - R ( 1 , 1 ) - R ( 2 , 2 ) + R ( 3 , 3 )
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
2010-05-26 21:22:54 +05:30
select case ( largest ( 1 ) )
2013-06-06 00:40:37 +05:30
case ( 1 )
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 )
case ( 2 )
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 )
case ( 3 )
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 )
case ( 4 )
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----------------------------------
2010-05-26 21:22:54 +05:30
end select
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 ) )
C = cos ( Euler ( 2 ) )
C2 = cos ( Euler ( 3 ) )
S1 = sin ( Euler ( 1 ) )
S = sin ( Euler ( 2 ) )
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
2013-01-31 21:58:08 +05:30
math_EulerToQ ( 1 ) = cos ( halfangles ( 1 ) + halfangles ( 3 ) ) * c
math_EulerToQ ( 2 ) = cos ( halfangles ( 1 ) - halfangles ( 3 ) ) * s
math_EulerToQ ( 3 ) = sin ( halfangles ( 1 ) - halfangles ( 3 ) ) * s
math_EulerToQ ( 4 ) = sin ( halfangles ( 1 ) + halfangles ( 3 ) ) * c
2013-06-06 00:40:37 +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
2011-02-25 14:55:53 +05:30
norm = sqrt ( math_mul3x3 ( axis , 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
2013-06-06 00:40:37 +05:30
math_EulerAxisAngleToR = transpose ( math_axisAngleToR ( axis , omega ) ) ! convert to passive rotation
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
math_EulerAxisAngleToQ = math_qConj ( math_axisAngleToQ ( axis , omega ) ) ! convert to passive rotation
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
2011-12-01 17:31:13 +05:30
real ( pReal ) :: s , c , norm
2010-05-06 19:37:21 +05:30
2011-02-25 14:55:53 +05:30
norm = sqrt ( math_mul3x3 ( axis , axis ) )
2010-05-06 19:37:21 +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
s = sin ( 0.5_pReal * omega )
c = cos ( 0.5_pReal * omega )
2013-06-06 00:40:37 +05:30
math_axisAngleToQ ( 1 ) = c
math_axisAngleToQ ( 2 : 4 ) = s * axisNrm ( 1 : 3 )
2010-05-06 19:37:21 +05:30
else
2013-06-06 00:40:37 +05:30
math_axisAngleToQ = [ 1.0_pReal , 0.0_pReal , 0.0_pReal , 0.0_pReal ] ! no rotation
2010-05-06 19:37:21 +05:30
endif
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 )
S = reshape ( [ 0.0_pReal , - q ( 4 ) , q ( 3 ) , &
q ( 4 ) , 0.0_pReal , - q ( 2 ) , &
2013-09-19 01:24:39 +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
2011-12-01 17:31:13 +05:30
real ( pReal ) :: acos_arg
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
math_qToEuler ( 2 ) = acos ( 1.0_pReal - 2.0_pReal * ( q ( 2 ) * q ( 2 ) + q ( 3 ) * q ( 3 ) ) )
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
acos_arg = q ( 1 )
if ( acos_arg > 1.0_pReal ) acos_arg = 1.0_pReal
if ( acos_arg < - 1.0_pReal ) acos_arg = - 1.0_pReal
math_qToEuler ( 1 ) = sign ( 2.0_pReal * acos ( acos_arg ) , 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
2013-06-06 00:40:37 +05:30
math_qToEuler ( 1 ) = atan2 ( q ( 1 ) * q ( 3 ) + q ( 2 ) * q ( 4 ) , q ( 1 ) * q ( 2 ) - q ( 3 ) * q ( 4 ) )
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
2013-06-06 00:40:37 +05:30
if ( math_qToEuler ( 1 ) < 0.0_pReal ) &
math_qToEuler ( 1 ) = math_qToEuler ( 1 ) + 2.0_pReal * pi
2013-01-31 21:58:08 +05:30
if ( math_qToEuler ( 2 ) < 0.0_pReal ) &
math_qToEuler ( 2 ) = math_qToEuler ( 2 ) + pi
2013-06-06 00:40:37 +05:30
if ( math_qToEuler ( 3 ) < 0.0_pReal ) &
math_qToEuler ( 3 ) = math_qToEuler ( 3 ) + 2.0_pReal * pi
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
2011-02-25 14:55:53 +05:30
halfAngle = acos ( max ( - 1.0_pReal , min ( 1.0_pReal , Q ( 1 ) ) ) ) ! limit to [-1,1] --> 0 to 180 deg
sinHalfAngle = sin ( halfAngle )
2012-08-25 17:16:36 +05:30
2010-05-06 19:37:21 +05:30
if ( sinHalfAngle < = 1.0e-4_pReal ) then ! very small rotation angle?
2013-01-31 21:58:08 +05:30
math_qToAxisAngle = 0.0_pReal
2010-04-29 15:31:09 +05:30
else
2013-01-31 21:58:08 +05:30
math_qToAxisAngle ( 1 : 3 ) = Q ( 2 : 4 ) / sinHalfAngle
math_qToAxisAngle ( 4 ) = halfAngle * 2.0_pReal
2010-04-29 15:31:09 +05:30
endif
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 )
2014-02-28 18:58:27 +05:30
use prec , only : &
2015-03-29 18:24:13 +05:30
DAMASK_NaN , &
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
2015-03-29 18:24:13 +05:30
if ( abs ( Q ( 1 ) ) > tol_math_check ) then ! unless rotation by 180 deg
2013-01-31 21:58:08 +05:30
math_qToRodrig = Q ( 2 : 4 ) / Q ( 1 )
2010-04-28 22:49:58 +05:30
else
2015-03-29 18:24:13 +05:30
math_qToRodrig = DAMASK_NaN ! NaN since Rodrig is unbound for 180 deg...
2010-04-28 22:49:58 +05:30
endif
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
2007-03-29 21:02:52 +05:30
tr = ( r ( 1 , 1 ) + r ( 2 , 2 ) + r ( 3 , 3 ) - 1.0_pReal ) * 0.4999999_pReal
2010-04-28 22:49:58 +05:30
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 )
2007-03-29 21:02:52 +05:30
math_sampleRandomOri ( 1 ) = rnd ( 1 ) * 2.0_pReal * pi
math_sampleRandomOri ( 2 ) = acos ( 2.0_pReal * rnd ( 2 ) - 1.0_pReal )
math_sampleRandomOri ( 3 ) = 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]
2007-03-29 21:02:52 +05:30
disturb ( 1 ) = scatter * rnd ( 1 ) ! phi1
disturb ( 2 ) = sign ( 1.0_pReal , rnd ( 2 ) ) * acos ( cosScatter + ( 1.0_pReal - cosScatter ) * rnd ( 4 ) ) ! Phi
disturb ( 3 ) = scatter * rnd ( 2 ) ! phi2
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
2014-02-28 18:58:27 +05:30
integer ( pInt ) , dimension ( 2 , 3 ) , parameter :: rotMap = reshape ( [ 2_pInt , 3_pInt , &
2011-12-01 17:31:13 +05:30
3_pInt , 1_pInt , &
2014-02-28 18:58:27 +05:30
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
fiberInC ( 1 ) = sin ( alpha ( 1 ) ) * cos ( alpha ( 2 ) )
fiberInC ( 2 ) = sin ( alpha ( 1 ) ) * sin ( alpha ( 2 ) )
fiberInC ( 3 ) = cos ( alpha ( 1 ) )
! fiber axis in sample coordinate system
fiberInS ( 1 ) = sin ( beta ( 1 ) ) * cos ( beta ( 2 ) )
fiberInS ( 2 ) = sin ( beta ( 1 ) ) * sin ( beta ( 2 ) )
fiberInS ( 3 ) = cos ( beta ( 1 ) )
! ---# 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)
2012-02-14 19:13:36 +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 ) )
2013-06-06 00:40:37 +05:30
oRot = math_EulerAxisAngleToR ( math_vectorproduct ( 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
if ( present ( width ) ) then
myWidth = width
else
myWidth = 3.0_pReal ! use +-3*sigma as default value for scatter
endif
do
call halton ( 2_pInt , rnd )
scatter = myWidth * ( 2.0_pReal * rnd ( 1 ) - 1.0_pReal )
if ( rnd ( 2 ) < = exp ( - 0.5_pReal * scatter ** 2.0_pReal ) ) & ! test if scattered value is drawn
exit
enddo
math_sampleGaussVar = scatter * stddev
end function math_sampleGaussVar
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2013-02-19 20:26:26 +05:30
!> @brief symmetric Euler angles for given symmetry 1: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
2009-03-04 17:18:54 +05:30
integer ( pInt ) , intent ( in ) :: sym
2008-02-15 18:12:27 +05:30
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: Euler
real ( pReal ) , dimension ( 3 , 3 ) :: math_symmetricEulers
2011-12-01 17:31:13 +05:30
integer ( pInt ) :: i , j
2012-08-25 17:16:36 +05:30
2008-02-15 18:12:27 +05:30
math_symmetricEulers ( 1 , 1 ) = pi + Euler ( 1 )
math_symmetricEulers ( 2 , 1 ) = Euler ( 2 )
math_symmetricEulers ( 3 , 1 ) = Euler ( 3 )
math_symmetricEulers ( 1 , 2 ) = pi - Euler ( 1 )
math_symmetricEulers ( 2 , 2 ) = pi - Euler ( 2 )
math_symmetricEulers ( 3 , 2 ) = pi + Euler ( 3 )
math_symmetricEulers ( 1 , 3 ) = 2.0_pReal * pi - Euler ( 1 )
math_symmetricEulers ( 2 , 3 ) = pi - Euler ( 2 )
math_symmetricEulers ( 3 , 3 ) = pi + Euler ( 3 )
2011-12-01 17:31:13 +05:30
forall ( i = 1_pInt : 3_pInt , j = 1_pInt : 3_pInt ) math_symmetricEulers ( j , i ) = modulo ( math_symmetricEulers ( j , i ) , 2.0_pReal * pi )
2008-02-15 18:12:27 +05:30
select case ( sym )
2011-12-01 17:31:13 +05:30
case ( 4_pInt ) ! all done
2008-02-15 18:12:27 +05:30
2011-12-01 17:31:13 +05:30
case ( 2_pInt ) ! return only first
math_symmetricEulers ( 1 : 3 , 2 : 3 ) = 0.0_pReal
2008-02-15 18:12:27 +05:30
case default ! return blank
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
!--------------------------------------------------------------------------------------------------
!> @brief not yet done
!--------------------------------------------------------------------------------------------------
2012-01-26 19:20:00 +05:30
subroutine math_spectralDecompositionSym33 ( M , values , vectors , error )
2011-09-14 18:56:00 +05:30
2012-08-25 17:16:36 +05:30
implicit none
2011-09-14 18:56:00 +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
logical , intent ( out ) :: error
2011-12-01 17:31:13 +05:30
integer ( pInt ) :: info
2011-09-14 18:56:00 +05:30
real ( pReal ) , dimension ( ( 64 + 2 ) * 3 ) :: work ! block size of 64 taken from http://www.netlib.org/lapack/double/dsyev.f
2012-08-25 17:16:36 +05:30
2011-09-14 18:56:00 +05:30
vectors = M ! copy matrix to input (doubles as output) array
2012-08-28 22:29:45 +05:30
#if(FLOAT==8)
call dsyev ( 'V' , 'U' , 3 , vectors , 3 , values , work , ( 64 + 2 ) * 3 , info )
#elif(FLOAT==4)
call ssyev ( 'V' , 'U' , 3 , vectors , 3 , values , work , ( 64 + 2 ) * 3 , info )
#endif
2011-09-14 18:56:00 +05:30
error = ( info == 0_pInt )
2012-08-25 17:16:36 +05:30
2011-09-14 18:56:00 +05:30
end subroutine
2011-02-04 21:11:32 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief EIGENWERTE UND EIGENWERTBASIS DER SYMMETRISCHEN 3X3 MATRIX M
!--------------------------------------------------------------------------------------------------
2013-01-31 21:58:08 +05:30
pure subroutine math_spectralDecomposition ( M , EW1 , EW2 , EW3 , EB1 , EB2 , EB3 )
2007-03-20 19:25:22 +05:30
implicit none
2011-12-01 17:31:13 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: M
real ( pReal ) , dimension ( 3 , 3 ) , intent ( out ) :: EB1 , EB2 , EB3
real ( pReal ) , intent ( out ) :: EW1 , EW2 , EW3
real ( pReal ) HI1M , HI2M , HI3M , R , S , T , P , Q , RHO , PHI , Y1 , Y2 , Y3 , D1 , D2 , D3
real ( pReal ) , parameter :: TOL = 1.e-14_pReal
real ( pReal ) , dimension ( 3 , 3 ) :: M1 , M2 , M3
real ( pReal ) C1 , C2 , C3 , arg
2013-01-31 21:58:08 +05:30
call math_hi ( M , HI1M , HI2M , HI3M )
2007-03-20 19:25:22 +05:30
R = - HI1M
S = HI2M
T = - HI3M
P = S - R ** 2.0_pReal / 3.0_pReal
Q = 2.0_pReal / 2 7.0_pReal * R ** 3.0_pReal - R * S / 3.0_pReal + T
EB1 = 0.0_pReal
EB2 = 0.0_pReal
EB3 = 0.0_pReal
2013-01-31 21:58:08 +05:30
if ( ( ABS ( P ) < TOL ) . AND . ( ABS ( Q ) < TOL ) ) then
2007-03-21 15:50:25 +05:30
! DREI GLEICHE EIGENWERTE
2007-03-20 19:25:22 +05:30
EW1 = HI1M / 3.0_pReal
EW2 = EW1
EW3 = EW1
2007-03-21 15:50:25 +05:30
! this is not really correct, but this way U is calculated
! correctly in PDECOMPOSITION (correct is EB?=I)
2007-03-20 19:25:22 +05:30
EB1 ( 1 , 1 ) = 1.0_pReal
EB2 ( 2 , 2 ) = 1.0_pReal
EB3 ( 3 , 3 ) = 1.0_pReal
2013-01-31 21:58:08 +05:30
else
2011-02-25 14:55:53 +05:30
RHO = sqrt ( - 3.0_pReal * P ** 3.0_pReal ) / 9.0_pReal
2007-03-20 19:25:22 +05:30
arg = - Q / RHO / 2.0_pReal
2013-01-31 21:58:08 +05:30
if ( arg > 1.0_pReal ) arg = 1.0_pReal
if ( arg < - 1.0_pReal ) arg = - 1.0_pReal
2011-02-25 14:55:53 +05:30
PHI = acos ( arg )
2011-12-01 17:31:13 +05:30
Y1 = 2.0_pReal * RHO ** ( 1.0_pReal / 3.0_pReal ) * cos ( PHI / 3.0_pReal )
Y2 = 2.0_pReal * RHO ** ( 1.0_pReal / 3.0_pReal ) * cos ( PHI / 3.0_pReal + 2.0_pReal / 3.0_pReal * PI )
Y3 = 2.0_pReal * RHO ** ( 1.0_pReal / 3.0_pReal ) * cos ( PHI / 3.0_pReal + 4.0_pReal / 3.0_pReal * PI )
2007-03-20 19:25:22 +05:30
EW1 = Y1 - R / 3.0_pReal
EW2 = Y2 - R / 3.0_pReal
EW3 = Y3 - R / 3.0_pReal
C1 = ABS ( EW1 - EW2 )
2012-08-25 17:16:36 +05:30
C2 = ABS ( EW2 - EW3 )
2007-03-20 19:25:22 +05:30
C3 = ABS ( EW3 - EW1 )
2013-01-31 21:58:08 +05:30
if ( C1 < TOL ) then
2007-03-21 15:50:25 +05:30
! EW1 is equal to EW2
2013-01-31 21:58:08 +05:30
D3 = 1.0_pReal / ( EW3 - EW1 ) / ( EW3 - EW2 )
M1 = M - EW1 * math_I3
M2 = M - EW2 * math_I3
EB3 = math_mul33x33 ( M1 , M2 ) * D3
2008-07-09 01:08:22 +05:30
2013-01-31 21:58:08 +05:30
EB1 = math_I3 - EB3
2007-03-21 15:50:25 +05:30
! both EB2 and EW2 are set to zero so that they do not
! contribute to U in PDECOMPOSITION
2013-01-31 21:58:08 +05:30
EW2 = 0.0_pReal
elseif ( C2 < TOL ) then
2007-03-21 15:50:25 +05:30
! EW2 is equal to EW3
2013-01-31 21:58:08 +05:30
D1 = 1.0_pReal / ( EW1 - EW2 ) / ( EW1 - EW3 )
M2 = M - math_I3 * EW2
M3 = M - math_I3 * EW3
EB1 = math_mul33x33 ( M2 , M3 ) * D1
EB2 = math_I3 - EB1
2007-03-21 15:50:25 +05:30
! both EB3 and EW3 are set to zero so that they do not
! contribute to U in PDECOMPOSITION
2013-01-31 21:58:08 +05:30
EW3 = 0.0_pReal
elseif ( C3 < TOL ) then
2007-03-21 15:50:25 +05:30
! EW1 is equal to EW3
2013-01-31 21:58:08 +05:30
D2 = 1.0_pReal / ( EW2 - EW1 ) / ( EW2 - EW3 )
M1 = M - math_I3 * EW1
M3 = M - math_I3 * EW3
EB2 = math_mul33x33 ( M1 , M3 ) * D2
EB1 = math_I3 - EB2
2007-03-21 15:50:25 +05:30
! both EB3 and EW3 are set to zero so that they do not
! contribute to U in PDECOMPOSITION
2013-01-31 21:58:08 +05:30
EW3 = 0.0_pReal
else
2007-03-21 15:50:25 +05:30
! all three eigenvectors are different
2013-01-31 21:58:08 +05:30
D1 = 1.0_pReal / ( EW1 - EW2 ) / ( EW1 - EW3 )
D2 = 1.0_pReal / ( EW2 - EW1 ) / ( EW2 - EW3 )
D3 = 1.0_pReal / ( EW3 - EW1 ) / ( EW3 - EW2 )
M1 = M - EW1 * math_I3
M2 = M - EW2 * math_I3
M3 = M - EW3 * math_I3
EB1 = math_mul33x33 ( M2 , M3 ) * D1
EB2 = math_mul33x33 ( M1 , M3 ) * D2
EB3 = math_mul33x33 ( M1 , M2 ) * D3
endif
endif
2008-07-09 01:08:22 +05:30
2013-01-31 21:58:08 +05:30
end subroutine math_spectralDecomposition
2011-08-01 15:41:32 +05:30
2013-01-31 21:58:08 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief FE = R.U
!--------------------------------------------------------------------------------------------------
pure subroutine math_pDecomposition ( FE , U , R , error )
implicit none
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: FE
real ( pReal ) , intent ( out ) , dimension ( 3 , 3 ) :: R , U
logical , intent ( out ) :: error
real ( pReal ) , dimension ( 3 , 3 ) :: CE , EB1 , EB2 , EB3 , UI
real ( pReal ) :: EW1 , EW2 , EW3 , det
error = . false .
ce = math_mul33x33 ( math_transpose33 ( FE ) , FE )
CALL math_spectralDecomposition ( CE , EW1 , EW2 , EW3 , EB1 , EB2 , EB3 )
U = sqrt ( EW1 ) * EB1 + sqrt ( EW2 ) * EB2 + sqrt ( EW3 ) * EB3
call math_invert33 ( U , UI , det , error )
if ( . not . error ) R = math_mul33x33 ( FE , UI )
end subroutine math_pDecomposition
2007-03-20 19:25:22 +05:30
2011-12-01 17:31:13 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Eigenvalues of symmetric 3X3 matrix M
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
function math_eigenvalues33 ( M )
2011-08-26 19:36:37 +05:30
implicit none
2011-12-01 17:31:13 +05:30
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: M
real ( pReal ) , dimension ( 3 , 3 ) :: EB1 = 0.0_pReal , EB2 = 0.0_pReal , EB3 = 0.0_pReal
2012-01-26 19:20:00 +05:30
real ( pReal ) , dimension ( 3 ) :: math_eigenvalues33
2011-12-01 17:31:13 +05:30
real ( pReal ) :: HI1M , HI2M , HI3M , R , S , T , P , Q , RHO , PHI , Y1 , Y2 , Y3 , arg
real ( pReal ) , parameter :: TOL = 1.e-14_pReal
2011-08-26 19:36:37 +05:30
CALL math_hi ( M , HI1M , HI2M , HI3M )
R = - HI1M
S = HI2M
T = - HI3M
P = S - R ** 2.0_pReal / 3.0_pReal
Q = 2.0_pReal / 2 7.0_pReal * R ** 3.0_pReal - R * S / 3.0_pReal + T
2012-08-25 17:16:36 +05:30
2011-08-26 19:36:37 +05:30
if ( ( abs ( P ) < TOL ) . and . ( abs ( Q ) < TOL ) ) THEN
! three equivalent eigenvalues
2012-01-26 19:20:00 +05:30
math_eigenvalues33 ( 1 ) = HI1M / 3.0_pReal
math_eigenvalues33 ( 2 ) = math_eigenvalues33 ( 1 )
math_eigenvalues33 ( 3 ) = math_eigenvalues33 ( 1 )
2011-08-26 19:36:37 +05:30
! this is not really correct, but this way U is calculated
! correctly in PDECOMPOSITION (correct is EB?=I)
EB1 ( 1 , 1 ) = 1.0_pReal
EB2 ( 2 , 2 ) = 1.0_pReal
EB3 ( 3 , 3 ) = 1.0_pReal
else
RHO = sqrt ( - 3.0_pReal * P ** 3.0_pReal ) / 9.0_pReal
arg = - Q / RHO / 2.0_pReal
2011-12-01 17:31:13 +05:30
if ( arg . GT . 1.0_pReal ) arg = 1.0_pReal
if ( arg . LT . - 1.0_pReal ) arg = - 1.0_pReal
2011-08-26 19:36:37 +05:30
PHI = acos ( arg )
Y1 = 2 * RHO ** ( 1.0_pReal / 3.0_pReal ) * cos ( PHI / 3.0_pReal )
Y2 = 2 * RHO ** ( 1.0_pReal / 3.0_pReal ) * cos ( PHI / 3.0_pReal + 2.0_pReal / 3.0_pReal * PI )
Y3 = 2 * RHO ** ( 1.0_pReal / 3.0_pReal ) * cos ( PHI / 3.0_pReal + 4.0_pReal / 3.0_pReal * PI )
2012-01-26 19:20:00 +05:30
math_eigenvalues33 ( 1 ) = Y1 - R / 3.0_pReal
math_eigenvalues33 ( 2 ) = Y2 - R / 3.0_pReal
math_eigenvalues33 ( 3 ) = Y3 - R / 3.0_pReal
2011-08-26 19:36:37 +05:30
endif
2012-03-09 01:55:28 +05:30
end function math_eigenvalues33
2007-03-20 19:25:22 +05:30
2011-12-01 17:31:13 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief HAUPTINVARIANTEN HI1M, HI2M, HI3M DER 3X3 MATRIX M
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
pure subroutine math_hi ( M , HI1M , HI2M , HI3M )
2012-08-25 17:16:36 +05:30
2007-03-20 19:25:22 +05:30
implicit none
2012-08-25 17:16:36 +05:30
real ( pReal ) , intent ( in ) :: M ( 3 , 3 )
real ( pReal ) , intent ( out ) :: HI1M , HI2M , HI3M
2007-03-21 15:50:25 +05:30
HI1M = M ( 1 , 1 ) + M ( 2 , 2 ) + M ( 3 , 3 )
2011-12-01 17:31:13 +05:30
HI2M = HI1M ** 2.0_pReal / 2.0_pReal - ( M ( 1 , 1 ) ** 2.0_pReal + M ( 2 , 2 ) ** 2.0_pReal + M ( 3 , 3 ) ** 2.0_pReal ) &
/ 2.0_pReal - M ( 1 , 2 ) * M ( 2 , 1 ) - M ( 1 , 3 ) * M ( 3 , 1 ) - M ( 2 , 3 ) * M ( 3 , 2 )
2012-01-26 19:20:00 +05:30
HI3M = math_det33 ( M )
2007-03-21 15:50:25 +05:30
2012-03-09 01:55:28 +05:30
end subroutine math_hi
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 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
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 )
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 ) :: i
integer ( pInt ) , save :: ndim_save = 0_pInt , seed = 1_pInt
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
2011-12-01 17:31:13 +05:30
if ( action_halton ( 1 : 1 ) == 'S' . or . action_halton ( 1 : 1 ) == 's' ) then
2007-03-20 19:25:22 +05:30
2011-12-01 17:31:13 +05:30
if ( name_halton ( 1 : 1 ) == 'B' . or . name_halton ( 1 : 1 ) == 'b' ) then
2007-03-20 19:25:22 +05:30
2011-12-01 17:31:13 +05:30
if ( ndim_save / = ndim ) then
deallocate ( base )
2010-05-06 19:37:21 +05:30
ndim_save = ndim
2011-12-01 17:31:13 +05:30
allocate ( base ( ndim_save ) )
endif
2007-03-20 19:25:22 +05:30
2011-12-01 17:31:13 +05:30
base ( 1 : ndim ) = value_halton ( 1 : ndim )
2007-03-20 19:25:22 +05:30
2011-12-01 17:31:13 +05:30
elseif ( name_halton ( 1 : 1 ) == 'N' . or . name_halton ( 1 : 1 ) == 'n' ) then
2007-03-20 19:25:22 +05:30
2011-12-01 17:31:13 +05:30
if ( ndim_save / = value_halton ( 1 ) ) then
deallocate ( base )
ndim_save = value_halton ( 1 )
allocate ( base ( ndim_save ) )
do i = 1_pInt , ndim_save
base ( i ) = prime ( i )
2010-05-06 19:37:21 +05:30
enddo
2007-03-20 19:25:22 +05:30
else
2011-12-01 17:31:13 +05:30
ndim_save = value_halton ( 1 )
endif
elseif ( name_halton ( 1 : 1 ) == 'S' . or . name_halton ( 1 : 1 ) == 's' ) then
seed = value_halton ( 1 )
endif
2013-03-07 01:04:30 +05:30
!--------------------------------------------------------------------------------------------------
! Get
2011-12-01 17:31:13 +05:30
elseif ( action_halton ( 1 : 1 ) == 'G' . or . action_halton ( 1 : 1 ) == 'g' ) then
if ( name_halton ( 1 : 1 ) == 'B' . or . name_halton ( 1 : 1 ) == 'b' ) then
if ( ndim / = ndim_save ) then
deallocate ( base )
2007-03-20 19:25:22 +05:30
ndim_save = ndim
2011-12-01 17:31:13 +05:30
allocate ( base ( ndim_save ) )
do i = 1_pInt , ndim_save
2007-03-20 19:25:22 +05:30
base ( i ) = prime ( i )
2009-06-29 20:59:07 +05:30
enddo
2011-12-01 17:31:13 +05:30
endif
value_halton ( 1 : ndim_save ) = base ( 1 : ndim_save )
elseif ( name_halton ( 1 : 1 ) == 'N' . or . name_halton ( 1 : 1 ) == 'n' ) then
value_halton ( 1 ) = ndim_save
elseif ( name_halton ( 1 : 1 ) == 'S' . or . name_halton ( 1 : 1 ) == 's' ) then
value_halton ( 1 ) = seed
endif
2013-03-07 01:04:30 +05:30
!--------------------------------------------------------------------------------------------------
! Increment
2011-12-01 17:31:13 +05:30
elseif ( action_halton ( 1 : 1 ) == 'I' . or . action_halton ( 1 : 1 ) == 'i' ) then
if ( name_halton ( 1 : 1 ) == 'S' . or . name_halton ( 1 : 1 ) == 's' ) then
seed = seed + value_halton ( 1 )
2007-03-20 19:25:22 +05:30
end if
2011-12-01 17:31:13 +05:30
endif
2007-03-21 15:50:25 +05:30
2012-03-09 01:55:28 +05:30
end subroutine halton_memory
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 the dimension for a Halton sequence
!> @author John Burkardt
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
subroutine halton_ndim_set ( ndim )
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
integer ( pInt ) , intent ( in ) :: ndim !< dimension of the Halton vectors
2011-12-01 17:31:13 +05:30
integer ( pInt ) :: value_halton ( 1 )
2007-03-20 19:25:22 +05:30
2011-12-01 17:31:13 +05:30
value_halton ( 1 ) = ndim
call halton_memory ( 'SET' , 'NDIM' , 1_pInt , value_halton )
2007-03-21 15:50:25 +05:30
2012-03-09 01:55:28 +05:30
end subroutine halton_ndim_set
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 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 )
2007-03-20 19:25:22 +05:30
implicit none
2007-03-21 15:50:25 +05:30
2011-12-01 17:31:13 +05:30
integer ( pInt ) , parameter :: ndim = 1_pInt
2013-03-07 01:04:30 +05:30
integer ( pInt ) , intent ( in ) :: seed !< seed for the Halton sequence.
2011-12-01 17:31:13 +05:30
integer ( pInt ) :: value_halton ( ndim )
2007-03-20 19:25:22 +05:30
2011-12-01 17:31:13 +05:30
value_halton ( 1 ) = seed
call halton_memory ( 'SET' , 'SEED' , ndim , value_halton )
2007-03-21 15:50:25 +05:30
2012-03-09 01:55:28 +05:30
end subroutine halton_seed_set
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 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
2012-08-27 13:34:47 +05:30
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
subroutine i_to_halton ( seed , base , ndim , r )
2013-03-07 01:04:30 +05:30
use IO , only : &
IO_error
2011-12-01 17:31:13 +05:30
implicit none
2013-03-07 01:04:30 +05:30
integer ( pInt ) , intent ( in ) :: ndim !< dimension of the sequence
integer ( pInt ) , intent ( in ) , dimension ( ndim ) :: base !< Halton bases
2011-12-01 17:31:13 +05:30
real ( pReal ) , dimension ( ndim ) :: base_inv
integer ( pInt ) , dimension ( ndim ) :: digit
2013-03-07 01:04:30 +05:30
real ( pReal ) , dimension ( ndim ) , intent ( out ) :: r !< the SEED-th element of the Halton sequence for the given bases
integer ( pInt ) , intent ( in ) :: seed !< index of the desired element
2011-12-01 17:31:13 +05:30
integer ( pInt ) , dimension ( ndim ) :: seed2
2007-03-21 15:50:25 +05:30
2011-12-01 17:31:13 +05:30
seed2 ( 1 : ndim ) = abs ( seed )
2007-03-20 19:25:22 +05:30
r ( 1 : ndim ) = 0.0_pReal
2012-02-13 23:11:27 +05:30
if ( any ( base ( 1 : ndim ) < = 1_pInt ) ) call IO_error ( error_ID = 405_pInt )
2007-03-20 19:25:22 +05:30
2011-12-01 17:31:13 +05:30
base_inv ( 1 : ndim ) = 1.0_pReal / real ( base ( 1 : ndim ) , pReal )
2007-03-20 19:25:22 +05:30
2011-12-01 17:31:13 +05:30
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 )
2007-03-20 19:25:22 +05:30
seed2 ( 1 : ndim ) = seed2 ( 1 : ndim ) / base ( 1 : ndim )
2009-06-29 20:59:07 +05:30
enddo
2007-03-20 19:25:22 +05:30
2012-03-09 01:55:28 +05:30
end subroutine i_to_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 returns any of the first 1500 prime numbers.
!> @details n <= 0 returns 1500, the index of the largest prime (12553) available.
!> @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
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2013-01-31 21:58:08 +05:30
integer ( pInt ) function prime ( n )
2013-03-07 01:04:30 +05:30
use IO , only : &
IO_error
2007-03-21 15:50:25 +05:30
2013-01-31 21:58:08 +05:30
implicit none
2013-03-07 01:04:30 +05:30
integer ( pInt ) , intent ( in ) :: n !< index of the desired prime number
2012-02-14 19:13:36 +05:30
integer ( pInt ) , parameter :: prime_max = 1500_pInt
2011-12-01 17:31:13 +05:30
integer ( pInt ) , save :: icall = 0_pInt
integer ( pInt ) , save , dimension ( prime_max ) :: npvec
2007-03-21 15:50:25 +05:30
2011-12-01 17:31:13 +05:30
if ( icall == 0_pInt ) then
icall = 1_pInt
2012-08-25 17:16:36 +05:30
npvec = [ &
2011-12-01 17:31:13 +05:30
2_pInt , 3_pInt , 5_pInt , 7_pInt , 11_pInt , 13_pInt , 17_pInt , 19_pInt , 23_pInt , 29_pInt , &
31_pInt , 37_pInt , 41_pInt , 43_pInt , 47_pInt , 53_pInt , 59_pInt , 61_pInt , 67_pInt , 71_pInt , &
73_pInt , 79_pInt , 83_pInt , 89_pInt , 97_pInt , 101_pInt , 103_pInt , 107_pInt , 109_pInt , 113_pInt , &
127_pInt , 131_pInt , 137_pInt , 139_pInt , 149_pInt , 151_pInt , 157_pInt , 163_pInt , 167_pInt , 173_pInt , &
179_pInt , 181_pInt , 191_pInt , 193_pInt , 197_pInt , 199_pInt , 211_pInt , 223_pInt , 227_pInt , 229_pInt , &
233_pInt , 239_pInt , 241_pInt , 251_pInt , 257_pInt , 263_pInt , 269_pInt , 271_pInt , 277_pInt , 281_pInt , &
283_pInt , 293_pInt , 307_pInt , 311_pInt , 313_pInt , 317_pInt , 331_pInt , 337_pInt , 347_pInt , 349_pInt , &
353_pInt , 359_pInt , 367_pInt , 373_pInt , 379_pInt , 383_pInt , 389_pInt , 397_pInt , 401_pInt , 409_pInt , &
419_pInt , 421_pInt , 431_pInt , 433_pInt , 439_pInt , 443_pInt , 449_pInt , 457_pInt , 461_pInt , 463_pInt , &
2012-08-25 17:16:36 +05:30
467_pInt , 479_pInt , 487_pInt , 491_pInt , 499_pInt , 503_pInt , 509_pInt , 521_pInt , 523_pInt , 541_pInt , &
! 101:200
2011-12-01 17:31:13 +05:30
547_pInt , 557_pInt , 563_pInt , 569_pInt , 571_pInt , 577_pInt , 587_pInt , 593_pInt , 599_pInt , 601_pInt , &
607_pInt , 613_pInt , 617_pInt , 619_pInt , 631_pInt , 641_pInt , 643_pInt , 647_pInt , 653_pInt , 659_pInt , &
661_pInt , 673_pInt , 677_pInt , 683_pInt , 691_pInt , 701_pInt , 709_pInt , 719_pInt , 727_pInt , 733_pInt , &
739_pInt , 743_pInt , 751_pInt , 757_pInt , 761_pInt , 769_pInt , 773_pInt , 787_pInt , 797_pInt , 809_pInt , &
811_pInt , 821_pInt , 823_pInt , 827_pInt , 829_pInt , 839_pInt , 853_pInt , 857_pInt , 859_pInt , 863_pInt , &
877_pInt , 881_pInt , 883_pInt , 887_pInt , 907_pInt , 911_pInt , 919_pInt , 929_pInt , 937_pInt , 941_pInt , &
947_pInt , 953_pInt , 967_pInt , 971_pInt , 977_pInt , 983_pInt , 991_pInt , 997_pInt , 1009_pInt , 1013_pInt , &
1019_pInt , 1021_pInt , 1031_pInt , 1033_pInt , 1039_pInt , 1049_pInt , 1051_pInt , 1061_pInt , 1063_pInt , 1069_pInt , &
1087_pInt , 1091_pInt , 1093_pInt , 1097_pInt , 1103_pInt , 1109_pInt , 1117_pInt , 1123_pInt , 1129_pInt , 1151_pInt , &
2012-08-25 17:16:36 +05:30
1153_pInt , 1163_pInt , 1171_pInt , 1181_pInt , 1187_pInt , 1193_pInt , 1201_pInt , 1213_pInt , 1217_pInt , 1223_pInt , &
! 201:300
2011-12-01 17:31:13 +05:30
1229_pInt , 1231_pInt , 1237_pInt , 1249_pInt , 1259_pInt , 1277_pInt , 1279_pInt , 1283_pInt , 1289_pInt , 1291_pInt , &
1297_pInt , 1301_pInt , 1303_pInt , 1307_pInt , 1319_pInt , 1321_pInt , 1327_pInt , 1361_pInt , 1367_pInt , 1373_pInt , &
1381_pInt , 1399_pInt , 1409_pInt , 1423_pInt , 1427_pInt , 1429_pInt , 1433_pInt , 1439_pInt , 1447_pInt , 1451_pInt , &
1453_pInt , 1459_pInt , 1471_pInt , 1481_pInt , 1483_pInt , 1487_pInt , 1489_pInt , 1493_pInt , 1499_pInt , 1511_pInt , &
1523_pInt , 1531_pInt , 1543_pInt , 1549_pInt , 1553_pInt , 1559_pInt , 1567_pInt , 1571_pInt , 1579_pInt , 1583_pInt , &
1597_pInt , 1601_pInt , 1607_pInt , 1609_pInt , 1613_pInt , 1619_pInt , 1621_pInt , 1627_pInt , 1637_pInt , 1657_pInt , &
1663_pInt , 1667_pInt , 1669_pInt , 1693_pInt , 1697_pInt , 1699_pInt , 1709_pInt , 1721_pInt , 1723_pInt , 1733_pInt , &
1741_pInt , 1747_pInt , 1753_pInt , 1759_pInt , 1777_pInt , 1783_pInt , 1787_pInt , 1789_pInt , 1801_pInt , 1811_pInt , &
1823_pInt , 1831_pInt , 1847_pInt , 1861_pInt , 1867_pInt , 1871_pInt , 1873_pInt , 1877_pInt , 1879_pInt , 1889_pInt , &
2012-08-25 17:16:36 +05:30
1901_pInt , 1907_pInt , 1913_pInt , 1931_pInt , 1933_pInt , 1949_pInt , 1951_pInt , 1973_pInt , 1979_pInt , 1987_pInt , &
! 301:400
2011-12-01 17:31:13 +05:30
1993_pInt , 1997_pInt , 1999_pInt , 2003_pInt , 2011_pInt , 2017_pInt , 2027_pInt , 2029_pInt , 2039_pInt , 2053_pInt , &
2063_pInt , 2069_pInt , 2081_pInt , 2083_pInt , 2087_pInt , 2089_pInt , 2099_pInt , 2111_pInt , 2113_pInt , 2129_pInt , &
2131_pInt , 2137_pInt , 2141_pInt , 2143_pInt , 2153_pInt , 2161_pInt , 2179_pInt , 2203_pInt , 2207_pInt , 2213_pInt , &
2221_pInt , 2237_pInt , 2239_pInt , 2243_pInt , 2251_pInt , 2267_pInt , 2269_pInt , 2273_pInt , 2281_pInt , 2287_pInt , &
2293_pInt , 2297_pInt , 2309_pInt , 2311_pInt , 2333_pInt , 2339_pInt , 2341_pInt , 2347_pInt , 2351_pInt , 2357_pInt , &
2371_pInt , 2377_pInt , 2381_pInt , 2383_pInt , 2389_pInt , 2393_pInt , 2399_pInt , 2411_pInt , 2417_pInt , 2423_pInt , &
2437_pInt , 2441_pInt , 2447_pInt , 2459_pInt , 2467_pInt , 2473_pInt , 2477_pInt , 2503_pInt , 2521_pInt , 2531_pInt , &
2539_pInt , 2543_pInt , 2549_pInt , 2551_pInt , 2557_pInt , 2579_pInt , 2591_pInt , 2593_pInt , 2609_pInt , 2617_pInt , &
2621_pInt , 2633_pInt , 2647_pInt , 2657_pInt , 2659_pInt , 2663_pInt , 2671_pInt , 2677_pInt , 2683_pInt , 2687_pInt , &
2012-08-25 17:16:36 +05:30
2689_pInt , 2693_pInt , 2699_pInt , 2707_pInt , 2711_pInt , 2713_pInt , 2719_pInt , 2729_pInt , 2731_pInt , 2741_pInt , &
! 401:500
2011-12-01 17:31:13 +05:30
2749_pInt , 2753_pInt , 2767_pInt , 2777_pInt , 2789_pInt , 2791_pInt , 2797_pInt , 2801_pInt , 2803_pInt , 2819_pInt , &
2833_pInt , 2837_pInt , 2843_pInt , 2851_pInt , 2857_pInt , 2861_pInt , 2879_pInt , 2887_pInt , 2897_pInt , 2903_pInt , &
2909_pInt , 2917_pInt , 2927_pInt , 2939_pInt , 2953_pInt , 2957_pInt , 2963_pInt , 2969_pInt , 2971_pInt , 2999_pInt , &
3001_pInt , 3011_pInt , 3019_pInt , 3023_pInt , 3037_pInt , 3041_pInt , 3049_pInt , 3061_pInt , 3067_pInt , 3079_pInt , &
3083_pInt , 3089_pInt , 3109_pInt , 3119_pInt , 3121_pInt , 3137_pInt , 3163_pInt , 3167_pInt , 3169_pInt , 3181_pInt , &
3187_pInt , 3191_pInt , 3203_pInt , 3209_pInt , 3217_pInt , 3221_pInt , 3229_pInt , 3251_pInt , 3253_pInt , 3257_pInt , &
3259_pInt , 3271_pInt , 3299_pInt , 3301_pInt , 3307_pInt , 3313_pInt , 3319_pInt , 3323_pInt , 3329_pInt , 3331_pInt , &
3343_pInt , 3347_pInt , 3359_pInt , 3361_pInt , 3371_pInt , 3373_pInt , 3389_pInt , 3391_pInt , 3407_pInt , 3413_pInt , &
3433_pInt , 3449_pInt , 3457_pInt , 3461_pInt , 3463_pInt , 3467_pInt , 3469_pInt , 3491_pInt , 3499_pInt , 3511_pInt , &
2012-08-25 17:16:36 +05:30
3517_pInt , 3527_pInt , 3529_pInt , 3533_pInt , 3539_pInt , 3541_pInt , 3547_pInt , 3557_pInt , 3559_pInt , 3571_pInt , &
! 501:600
2011-12-01 17:31:13 +05:30
3581_pInt , 3583_pInt , 3593_pInt , 3607_pInt , 3613_pInt , 3617_pInt , 3623_pInt , 3631_pInt , 3637_pInt , 3643_pInt , &
3659_pInt , 3671_pInt , 3673_pInt , 3677_pInt , 3691_pInt , 3697_pInt , 3701_pInt , 3709_pInt , 3719_pInt , 3727_pInt , &
3733_pInt , 3739_pInt , 3761_pInt , 3767_pInt , 3769_pInt , 3779_pInt , 3793_pInt , 3797_pInt , 3803_pInt , 3821_pInt , &
3823_pInt , 3833_pInt , 3847_pInt , 3851_pInt , 3853_pInt , 3863_pInt , 3877_pInt , 3881_pInt , 3889_pInt , 3907_pInt , &
3911_pInt , 3917_pInt , 3919_pInt , 3923_pInt , 3929_pInt , 3931_pInt , 3943_pInt , 3947_pInt , 3967_pInt , 3989_pInt , &
4001_pInt , 4003_pInt , 4007_pInt , 4013_pInt , 4019_pInt , 4021_pInt , 4027_pInt , 4049_pInt , 4051_pInt , 4057_pInt , &
4073_pInt , 4079_pInt , 4091_pInt , 4093_pInt , 4099_pInt , 4111_pInt , 4127_pInt , 4129_pInt , 4133_pInt , 4139_pInt , &
4153_pInt , 4157_pInt , 4159_pInt , 4177_pInt , 4201_pInt , 4211_pInt , 4217_pInt , 4219_pInt , 4229_pInt , 4231_pInt , &
4241_pInt , 4243_pInt , 4253_pInt , 4259_pInt , 4261_pInt , 4271_pInt , 4273_pInt , 4283_pInt , 4289_pInt , 4297_pInt , &
2012-08-25 17:16:36 +05:30
4327_pInt , 4337_pInt , 4339_pInt , 4349_pInt , 4357_pInt , 4363_pInt , 4373_pInt , 4391_pInt , 4397_pInt , 4409_pInt , &
! 601:700
2011-12-01 17:31:13 +05:30
4421_pInt , 4423_pInt , 4441_pInt , 4447_pInt , 4451_pInt , 4457_pInt , 4463_pInt , 4481_pInt , 4483_pInt , 4493_pInt , &
4507_pInt , 4513_pInt , 4517_pInt , 4519_pInt , 4523_pInt , 4547_pInt , 4549_pInt , 4561_pInt , 4567_pInt , 4583_pInt , &
4591_pInt , 4597_pInt , 4603_pInt , 4621_pInt , 4637_pInt , 4639_pInt , 4643_pInt , 4649_pInt , 4651_pInt , 4657_pInt , &
4663_pInt , 4673_pInt , 4679_pInt , 4691_pInt , 4703_pInt , 4721_pInt , 4723_pInt , 4729_pInt , 4733_pInt , 4751_pInt , &
4759_pInt , 4783_pInt , 4787_pInt , 4789_pInt , 4793_pInt , 4799_pInt , 4801_pInt , 4813_pInt , 4817_pInt , 4831_pInt , &
4861_pInt , 4871_pInt , 4877_pInt , 4889_pInt , 4903_pInt , 4909_pInt , 4919_pInt , 4931_pInt , 4933_pInt , 4937_pInt , &
4943_pInt , 4951_pInt , 4957_pInt , 4967_pInt , 4969_pInt , 4973_pInt , 4987_pInt , 4993_pInt , 4999_pInt , 5003_pInt , &
5009_pInt , 5011_pInt , 5021_pInt , 5023_pInt , 5039_pInt , 5051_pInt , 5059_pInt , 5077_pInt , 5081_pInt , 5087_pInt , &
5099_pInt , 5101_pInt , 5107_pInt , 5113_pInt , 5119_pInt , 5147_pInt , 5153_pInt , 5167_pInt , 5171_pInt , 5179_pInt , &
2012-08-25 17:16:36 +05:30
5189_pInt , 5197_pInt , 5209_pInt , 5227_pInt , 5231_pInt , 5233_pInt , 5237_pInt , 5261_pInt , 5273_pInt , 5279_pInt , &
! 701:800
2011-12-01 17:31:13 +05:30
5281_pInt , 5297_pInt , 5303_pInt , 5309_pInt , 5323_pInt , 5333_pInt , 5347_pInt , 5351_pInt , 5381_pInt , 5387_pInt , &
5393_pInt , 5399_pInt , 5407_pInt , 5413_pInt , 5417_pInt , 5419_pInt , 5431_pInt , 5437_pInt , 5441_pInt , 5443_pInt , &
5449_pInt , 5471_pInt , 5477_pInt , 5479_pInt , 5483_pInt , 5501_pInt , 5503_pInt , 5507_pInt , 5519_pInt , 5521_pInt , &
5527_pInt , 5531_pInt , 5557_pInt , 5563_pInt , 5569_pInt , 5573_pInt , 5581_pInt , 5591_pInt , 5623_pInt , 5639_pInt , &
5641_pInt , 5647_pInt , 5651_pInt , 5653_pInt , 5657_pInt , 5659_pInt , 5669_pInt , 5683_pInt , 5689_pInt , 5693_pInt , &
5701_pInt , 5711_pInt , 5717_pInt , 5737_pInt , 5741_pInt , 5743_pInt , 5749_pInt , 5779_pInt , 5783_pInt , 5791_pInt , &
5801_pInt , 5807_pInt , 5813_pInt , 5821_pInt , 5827_pInt , 5839_pInt , 5843_pInt , 5849_pInt , 5851_pInt , 5857_pInt , &
5861_pInt , 5867_pInt , 5869_pInt , 5879_pInt , 5881_pInt , 5897_pInt , 5903_pInt , 5923_pInt , 5927_pInt , 5939_pInt , &
5953_pInt , 5981_pInt , 5987_pInt , 6007_pInt , 6011_pInt , 6029_pInt , 6037_pInt , 6043_pInt , 6047_pInt , 6053_pInt , &
2012-08-25 17:16:36 +05:30
6067_pInt , 6073_pInt , 6079_pInt , 6089_pInt , 6091_pInt , 6101_pInt , 6113_pInt , 6121_pInt , 6131_pInt , 6133_pInt , &
! 801:900
2011-12-01 17:31:13 +05:30
6143_pInt , 6151_pInt , 6163_pInt , 6173_pInt , 6197_pInt , 6199_pInt , 6203_pInt , 6211_pInt , 6217_pInt , 6221_pInt , &
6229_pInt , 6247_pInt , 6257_pInt , 6263_pInt , 6269_pInt , 6271_pInt , 6277_pInt , 6287_pInt , 6299_pInt , 6301_pInt , &
6311_pInt , 6317_pInt , 6323_pInt , 6329_pInt , 6337_pInt , 6343_pInt , 6353_pInt , 6359_pInt , 6361_pInt , 6367_pInt , &
6373_pInt , 6379_pInt , 6389_pInt , 6397_pInt , 6421_pInt , 6427_pInt , 6449_pInt , 6451_pInt , 6469_pInt , 6473_pInt , &
6481_pInt , 6491_pInt , 6521_pInt , 6529_pInt , 6547_pInt , 6551_pInt , 6553_pInt , 6563_pInt , 6569_pInt , 6571_pInt , &
6577_pInt , 6581_pInt , 6599_pInt , 6607_pInt , 6619_pInt , 6637_pInt , 6653_pInt , 6659_pInt , 6661_pInt , 6673_pInt , &
6679_pInt , 6689_pInt , 6691_pInt , 6701_pInt , 6703_pInt , 6709_pInt , 6719_pInt , 6733_pInt , 6737_pInt , 6761_pInt , &
6763_pInt , 6779_pInt , 6781_pInt , 6791_pInt , 6793_pInt , 6803_pInt , 6823_pInt , 6827_pInt , 6829_pInt , 6833_pInt , &
6841_pInt , 6857_pInt , 6863_pInt , 6869_pInt , 6871_pInt , 6883_pInt , 6899_pInt , 6907_pInt , 6911_pInt , 6917_pInt , &
2012-08-25 17:16:36 +05:30
6947_pInt , 6949_pInt , 6959_pInt , 6961_pInt , 6967_pInt , 6971_pInt , 6977_pInt , 6983_pInt , 6991_pInt , 6997_pInt , &
! 901:1000
2011-12-01 17:31:13 +05:30
7001_pInt , 7013_pInt , 7019_pInt , 7027_pInt , 7039_pInt , 7043_pInt , 7057_pInt , 7069_pInt , 7079_pInt , 7103_pInt , &
7109_pInt , 7121_pInt , 7127_pInt , 7129_pInt , 7151_pInt , 7159_pInt , 7177_pInt , 7187_pInt , 7193_pInt , 7207_pInt , &
7211_pInt , 7213_pInt , 7219_pInt , 7229_pInt , 7237_pInt , 7243_pInt , 7247_pInt , 7253_pInt , 7283_pInt , 7297_pInt , &
7307_pInt , 7309_pInt , 7321_pInt , 7331_pInt , 7333_pInt , 7349_pInt , 7351_pInt , 7369_pInt , 7393_pInt , 7411_pInt , &
7417_pInt , 7433_pInt , 7451_pInt , 7457_pInt , 7459_pInt , 7477_pInt , 7481_pInt , 7487_pInt , 7489_pInt , 7499_pInt , &
7507_pInt , 7517_pInt , 7523_pInt , 7529_pInt , 7537_pInt , 7541_pInt , 7547_pInt , 7549_pInt , 7559_pInt , 7561_pInt , &
7573_pInt , 7577_pInt , 7583_pInt , 7589_pInt , 7591_pInt , 7603_pInt , 7607_pInt , 7621_pInt , 7639_pInt , 7643_pInt , &
7649_pInt , 7669_pInt , 7673_pInt , 7681_pInt , 7687_pInt , 7691_pInt , 7699_pInt , 7703_pInt , 7717_pInt , 7723_pInt , &
7727_pInt , 7741_pInt , 7753_pInt , 7757_pInt , 7759_pInt , 7789_pInt , 7793_pInt , 7817_pInt , 7823_pInt , 7829_pInt , &
2012-08-25 17:16:36 +05:30
7841_pInt , 7853_pInt , 7867_pInt , 7873_pInt , 7877_pInt , 7879_pInt , 7883_pInt , 7901_pInt , 7907_pInt , 7919_pInt , &
! 1001:1100
2011-12-01 17:31:13 +05:30
7927_pInt , 7933_pInt , 7937_pInt , 7949_pInt , 7951_pInt , 7963_pInt , 7993_pInt , 8009_pInt , 8011_pInt , 8017_pInt , &
8039_pInt , 8053_pInt , 8059_pInt , 8069_pInt , 8081_pInt , 8087_pInt , 8089_pInt , 8093_pInt , 8101_pInt , 8111_pInt , &
8117_pInt , 8123_pInt , 8147_pInt , 8161_pInt , 8167_pInt , 8171_pInt , 8179_pInt , 8191_pInt , 8209_pInt , 8219_pInt , &
8221_pInt , 8231_pInt , 8233_pInt , 8237_pInt , 8243_pInt , 8263_pInt , 8269_pInt , 8273_pInt , 8287_pInt , 8291_pInt , &
8293_pInt , 8297_pInt , 8311_pInt , 8317_pInt , 8329_pInt , 8353_pInt , 8363_pInt , 8369_pInt , 8377_pInt , 8387_pInt , &
8389_pInt , 8419_pInt , 8423_pInt , 8429_pInt , 8431_pInt , 8443_pInt , 8447_pInt , 8461_pInt , 8467_pInt , 8501_pInt , &
8513_pInt , 8521_pInt , 8527_pInt , 8537_pInt , 8539_pInt , 8543_pInt , 8563_pInt , 8573_pInt , 8581_pInt , 8597_pInt , &
8599_pInt , 8609_pInt , 8623_pInt , 8627_pInt , 8629_pInt , 8641_pInt , 8647_pInt , 8663_pInt , 8669_pInt , 8677_pInt , &
8681_pInt , 8689_pInt , 8693_pInt , 8699_pInt , 8707_pInt , 8713_pInt , 8719_pInt , 8731_pInt , 8737_pInt , 8741_pInt , &
2012-08-25 17:16:36 +05:30
8747_pInt , 8753_pInt , 8761_pInt , 8779_pInt , 8783_pInt , 8803_pInt , 8807_pInt , 8819_pInt , 8821_pInt , 8831_pInt , &
! 1101:1200
2011-12-01 17:31:13 +05:30
8837_pInt , 8839_pInt , 8849_pInt , 8861_pInt , 8863_pInt , 8867_pInt , 8887_pInt , 8893_pInt , 8923_pInt , 8929_pInt , &
8933_pInt , 8941_pInt , 8951_pInt , 8963_pInt , 8969_pInt , 8971_pInt , 8999_pInt , 9001_pInt , 9007_pInt , 9011_pInt , &
9013_pInt , 9029_pInt , 9041_pInt , 9043_pInt , 9049_pInt , 9059_pInt , 9067_pInt , 9091_pInt , 9103_pInt , 9109_pInt , &
9127_pInt , 9133_pInt , 9137_pInt , 9151_pInt , 9157_pInt , 9161_pInt , 9173_pInt , 9181_pInt , 9187_pInt , 9199_pInt , &
9203_pInt , 9209_pInt , 9221_pInt , 9227_pInt , 9239_pInt , 9241_pInt , 9257_pInt , 9277_pInt , 9281_pInt , 9283_pInt , &
9293_pInt , 9311_pInt , 9319_pInt , 9323_pInt , 9337_pInt , 9341_pInt , 9343_pInt , 9349_pInt , 9371_pInt , 9377_pInt , &
9391_pInt , 9397_pInt , 9403_pInt , 9413_pInt , 9419_pInt , 9421_pInt , 9431_pInt , 9433_pInt , 9437_pInt , 9439_pInt , &
9461_pInt , 9463_pInt , 9467_pInt , 9473_pInt , 9479_pInt , 9491_pInt , 9497_pInt , 9511_pInt , 9521_pInt , 9533_pInt , &
9539_pInt , 9547_pInt , 9551_pInt , 9587_pInt , 9601_pInt , 9613_pInt , 9619_pInt , 9623_pInt , 9629_pInt , 9631_pInt , &
2012-08-25 17:16:36 +05:30
9643_pInt , 9649_pInt , 9661_pInt , 9677_pInt , 9679_pInt , 9689_pInt , 9697_pInt , 9719_pInt , 9721_pInt , 9733_pInt , &
! 1201:1300
2011-12-01 17:31:13 +05:30
9739_pInt , 9743_pInt , 9749_pInt , 9767_pInt , 9769_pInt , 9781_pInt , 9787_pInt , 9791_pInt , 9803_pInt , 9811_pInt , &
9817_pInt , 9829_pInt , 9833_pInt , 9839_pInt , 9851_pInt , 9857_pInt , 9859_pInt , 9871_pInt , 9883_pInt , 9887_pInt , &
9901_pInt , 9907_pInt , 9923_pInt , 9929_pInt , 9931_pInt , 9941_pInt , 9949_pInt , 9967_pInt , 9973_pInt , 10007_pInt , &
10009_pInt , 10037_pInt , 10039_pInt , 10061_pInt , 10067_pInt , 10069_pInt , 10079_pInt , 10091_pInt , 10093_pInt , 10099_pInt , &
10103_pInt , 10111_pInt , 10133_pInt , 10139_pInt , 10141_pInt , 10151_pInt , 10159_pInt , 10163_pInt , 10169_pInt , 10177_pInt , &
10181_pInt , 10193_pInt , 10211_pInt , 10223_pInt , 10243_pInt , 10247_pInt , 10253_pInt , 10259_pInt , 10267_pInt , 10271_pInt , &
10273_pInt , 10289_pInt , 10301_pInt , 10303_pInt , 10313_pInt , 10321_pInt , 10331_pInt , 10333_pInt , 10337_pInt , 10343_pInt , &
10357_pInt , 10369_pInt , 10391_pInt , 10399_pInt , 10427_pInt , 10429_pInt , 10433_pInt , 10453_pInt , 10457_pInt , 10459_pInt , &
10463_pInt , 10477_pInt , 10487_pInt , 10499_pInt , 10501_pInt , 10513_pInt , 10529_pInt , 10531_pInt , 10559_pInt , 10567_pInt , &
2012-08-25 17:16:36 +05:30
10589_pInt , 10597_pInt , 10601_pInt , 10607_pInt , 10613_pInt , 10627_pInt , 10631_pInt , 10639_pInt , 10651_pInt , 10657_pInt , &
! 1301:1400
2011-12-01 17:31:13 +05:30
10663_pInt , 10667_pInt , 10687_pInt , 10691_pInt , 10709_pInt , 10711_pInt , 10723_pInt , 10729_pInt , 10733_pInt , 10739_pInt , &
10753_pInt , 10771_pInt , 10781_pInt , 10789_pInt , 10799_pInt , 10831_pInt , 10837_pInt , 10847_pInt , 10853_pInt , 10859_pInt , &
10861_pInt , 10867_pInt , 10883_pInt , 10889_pInt , 10891_pInt , 10903_pInt , 10909_pInt , 19037_pInt , 10939_pInt , 10949_pInt , &
10957_pInt , 10973_pInt , 10979_pInt , 10987_pInt , 10993_pInt , 11003_pInt , 11027_pInt , 11047_pInt , 11057_pInt , 11059_pInt , &
11069_pInt , 11071_pInt , 11083_pInt , 11087_pInt , 11093_pInt , 11113_pInt , 11117_pInt , 11119_pInt , 11131_pInt , 11149_pInt , &
11159_pInt , 11161_pInt , 11171_pInt , 11173_pInt , 11177_pInt , 11197_pInt , 11213_pInt , 11239_pInt , 11243_pInt , 11251_pInt , &
11257_pInt , 11261_pInt , 11273_pInt , 11279_pInt , 11287_pInt , 11299_pInt , 11311_pInt , 11317_pInt , 11321_pInt , 11329_pInt , &
11351_pInt , 11353_pInt , 11369_pInt , 11383_pInt , 11393_pInt , 11399_pInt , 11411_pInt , 11423_pInt , 11437_pInt , 11443_pInt , &
11447_pInt , 11467_pInt , 11471_pInt , 11483_pInt , 11489_pInt , 11491_pInt , 11497_pInt , 11503_pInt , 11519_pInt , 11527_pInt , &
2012-08-25 17:16:36 +05:30
11549_pInt , 11551_pInt , 11579_pInt , 11587_pInt , 11593_pInt , 11597_pInt , 11617_pInt , 11621_pInt , 11633_pInt , 11657_pInt , &
! 1401:1500
2011-12-01 17:31:13 +05:30
11677_pInt , 11681_pInt , 11689_pInt , 11699_pInt , 11701_pInt , 11717_pInt , 11719_pInt , 11731_pInt , 11743_pInt , 11777_pInt , &
11779_pInt , 11783_pInt , 11789_pInt , 11801_pInt , 11807_pInt , 11813_pInt , 11821_pInt , 11827_pInt , 11831_pInt , 11833_pInt , &
11839_pInt , 11863_pInt , 11867_pInt , 11887_pInt , 11897_pInt , 11903_pInt , 11909_pInt , 11923_pInt , 11927_pInt , 11933_pInt , &
11939_pInt , 11941_pInt , 11953_pInt , 11959_pInt , 11969_pInt , 11971_pInt , 11981_pInt , 11987_pInt , 12007_pInt , 12011_pInt , &
12037_pInt , 12041_pInt , 12043_pInt , 12049_pInt , 12071_pInt , 12073_pInt , 12097_pInt , 12101_pInt , 12107_pInt , 12109_pInt , &
12113_pInt , 12119_pInt , 12143_pInt , 12149_pInt , 12157_pInt , 12161_pInt , 12163_pInt , 12197_pInt , 12203_pInt , 12211_pInt , &
12227_pInt , 12239_pInt , 12241_pInt , 12251_pInt , 12253_pInt , 12263_pInt , 12269_pInt , 12277_pInt , 12281_pInt , 12289_pInt , &
12301_pInt , 12323_pInt , 12329_pInt , 12343_pInt , 12347_pInt , 12373_pInt , 12377_pInt , 12379_pInt , 12391_pInt , 12401_pInt , &
12409_pInt , 12413_pInt , 12421_pInt , 12433_pInt , 12437_pInt , 12451_pInt , 12457_pInt , 12473_pInt , 12479_pInt , 12487_pInt , &
2012-08-25 17:16:36 +05:30
12491_pInt , 12497_pInt , 12503_pInt , 12511_pInt , 12517_pInt , 12527_pInt , 12539_pInt , 12541_pInt , 12547_pInt , 12553_pInt ]
2011-12-01 17:31:13 +05:30
endif
2007-03-20 19:25:22 +05:30
2013-03-07 01:04:30 +05:30
if ( n < 0_pInt ) then
2007-03-20 19:25:22 +05:30
prime = prime_max
2011-12-01 17:31:13 +05:30
else if ( n == 0_pInt ) then
prime = 1_pInt
else if ( n < = prime_max ) then
2007-03-20 19:25:22 +05:30
prime = npvec ( n )
2012-02-13 19:38:07 +05:30
else
2013-02-11 15:14:17 +05:30
prime = - 1_pInt
2012-02-13 23:11:27 +05:30
call IO_error ( error_ID = 406_pInt )
2007-03-20 19:25:22 +05:30
end if
2013-02-11 15:14:17 +05:30
2012-03-09 01:55:28 +05:30
end function prime
2007-03-20 19:25:22 +05:30
2011-12-01 17:31:13 +05:30
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
math_areaTriangle = 0.5_pReal * math_norm3 ( math_vectorproduct ( v1 - v2 , v1 - v3 ) )
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
2012-06-15 21:40:21 +05:30
#ifdef Spectral
2012-08-27 13:34:47 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculates curl field using differentation in Fourier space
2013-11-22 16:01:52 +05:30
!> @todo enable odd resolution
2012-08-27 13:34:47 +05:30
!--------------------------------------------------------------------------------------------------
function math_curlFFT ( geomdim , field )
2013-03-07 01:04:30 +05:30
use IO , only : &
IO_error
use numerics , only : &
fftw_timelimit , &
fftw_planner_flag
use debug , only : &
debug_math , &
debug_level , &
debug_levelBasic
2012-03-09 01:55:28 +05:30
2011-12-01 17:31:13 +05:30
implicit none
2013-03-07 13:28:42 +05:30
real ( pReal ) , dimension ( : , : , : , : , : ) , intent ( in ) :: field !< field of data, first three dimensions are resolution, 4th is 1 or 3 (vector or tensor), 5th is 3
real ( pReal ) , dimension ( size ( field , 1 ) , size ( field , 2 ) , size ( field , 3 ) , size ( field , 4 ) , size ( field , 5 ) ) :: &
math_curlFFT
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: geomdim !< physical length dimension in three directions
real ( pReal ) , dimension ( : , : , : , : , : ) , pointer :: field_real
complex ( pReal ) , dimension ( : , : , : , : , : ) , pointer :: field_fourier
real ( pReal ) , dimension ( : , : , : , : , : ) , pointer :: curl_real
complex ( pReal ) , dimension ( : , : , : , : , : ) , pointer :: curl_fourier
integer ( pInt ) , dimension ( 3 ) :: &
k_s , &
res
complex ( pReal ) , dimension ( 3 ) :: &
xi
type ( C_PTR ) :: fftw_forth , fftw_back
type ( C_PTR ) :: field_fftw , curl_fftw
integer ( pInt ) :: i , j , k , l , res1_red , vec_tens
real ( pReal ) :: wgt
2013-03-07 01:04:30 +05:30
2012-08-27 13:34:47 +05:30
res = [ size ( field , 1 ) , size ( field , 2 ) , size ( field , 3 ) ]
vec_tens = size ( field , 4 )
2012-12-14 23:00:22 +05:30
2012-07-05 15:24:50 +05:30
if ( iand ( debug_level ( debug_math ) , debug_levelBasic ) / = 0_pInt ) then
2013-01-31 21:58:08 +05:30
if ( vec_tens == 1_pInt ) write ( 6 , '(a)' ) 'Calculating curl of vector field'
if ( vec_tens == 3_pInt ) write ( 6 , '(a)' ) 'Calculating curl of tensor field'
write ( 6 , '(a,3(e12.5))' ) ' Dimension: ' , geomdim
write ( 6 , '(a,3(i5))' ) ' Resolution:' , res
2012-01-13 21:48:16 +05:30
endif
2012-12-14 23:00:22 +05:30
2013-03-07 13:28:42 +05:30
!--------------------------------------------------------------------------------------------------
! sanity checks
2012-12-14 23:00:22 +05:30
if ( vec_tens / = 1_pInt . and . vec_tens / = 3_pInt ) &
call IO_error ( 0_pInt , ext_msg = 'Invalid data type in math_curlFFT' )
if ( ( mod ( res ( 3 ) , 2_pInt ) / = 0_pInt . and . res ( 3 ) / = 1_pInt ) . or . &
mod ( res ( 2 ) , 2_pInt ) / = 0_pInt . or . &
mod ( res ( 1 ) , 2_pInt ) / = 0_pInt ) &
call IO_error ( 0_pInt , ext_msg = 'Resolution in math_curlFFT' )
if ( pReal / = C_DOUBLE . or . pInt / = C_INT ) &
call IO_error ( 0_pInt , ext_msg = 'Fortran to C in math_curlFFT' )
2013-03-07 13:28:42 +05:30
wgt = 1.0_pReal / real ( product ( res ) , pReal )
res1_red = res ( 1 ) / 2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c)
2012-01-25 16:00:39 +05:30
2013-03-07 13:28:42 +05:30
!--------------------------------------------------------------------------------------------------
! allocation and FFTW initialization
2012-01-13 21:48:16 +05:30
call fftw_set_timelimit ( fftw_timelimit )
field_fftw = fftw_alloc_complex ( int ( res1_red * res ( 2 ) * res ( 3 ) * vec_tens * 3_pInt , C_SIZE_T ) ) !C_SIZE_T is of type integer(8)
2012-02-09 21:28:15 +05:30
call c_f_pointer ( field_fftw , field_real , [ res ( 1 ) + 2_pInt , res ( 2 ) , res ( 3 ) , vec_tens , 3_pInt ] )
call c_f_pointer ( field_fftw , field_fourier , [ res1_red , res ( 2 ) , res ( 3 ) , vec_tens , 3_pInt ] )
2013-03-07 13:28:42 +05:30
curl_fftw = fftw_alloc_complex ( int ( res1_red * res ( 2 ) * res ( 3 ) * vec_tens * 3_pInt , C_SIZE_T ) ) !C_SIZE_T is of type integer(8)
2012-02-09 21:28:15 +05:30
call c_f_pointer ( curl_fftw , curl_real , [ res ( 1 ) + 2_pInt , res ( 2 ) , res ( 3 ) , vec_tens , 3_pInt ] )
call c_f_pointer ( curl_fftw , curl_fourier , [ res1_red , res ( 2 ) , res ( 3 ) , vec_tens , 3_pInt ] )
2012-01-13 21:48:16 +05:30
2013-03-07 13:28:42 +05:30
fftw_forth = fftw_plan_many_dft_r2c ( 3_pInt , [ res ( 3 ) , res ( 2 ) , res ( 1 ) ] , vec_tens * 3_pInt , & ! dimensions, length in each dimension in reversed order, total # of transforms
field_real , [ res ( 3 ) , res ( 2 ) , res ( 1 ) + 2_pInt ] , & ! input data, physical length in each dimension in reversed order
1_pInt , res ( 3 ) * res ( 2 ) * ( res ( 1 ) + 2_pInt ) , & ! striding, product of physical lenght in the 3 dimensions
field_fourier , [ res ( 3 ) , res ( 2 ) , res1_red ] , & ! output data, physical length in each dimension in reversed order
1_pInt , res ( 3 ) * res ( 2 ) * res1_red , fftw_planner_flag ) ! striding, product of physical lenght in the 3 dimensions, planner mode
2013-03-07 01:04:30 +05:30
fftw_back = fftw_plan_many_dft_c2r ( 3_pInt , [ res ( 3 ) , res ( 2 ) , res ( 1 ) ] , vec_tens * 3_pInt , &
curl_fourier , [ res ( 3 ) , res ( 2 ) , res1_red ] , &
2012-02-09 21:28:15 +05:30
1_pInt , res ( 3 ) * res ( 2 ) * res1_red , &
2013-03-07 01:04:30 +05:30
curl_real , [ res ( 3 ) , res ( 2 ) , res ( 1 ) + 2_pInt ] , &
2012-02-09 21:28:15 +05:30
1_pInt , res ( 3 ) * res ( 2 ) * ( res ( 1 ) + 2_pInt ) , fftw_planner_flag )
2011-12-01 17:31:13 +05:30
2013-03-07 13:28:42 +05:30
field_real ( 1 : res ( 1 ) , 1 : res ( 2 ) , 1 : res ( 3 ) , 1 : vec_tens , 1 : 3 ) = field ! field_real is overwritten during plan creation and is larger (padding)
2011-12-01 17:31:13 +05:30
2013-03-07 13:28:42 +05:30
!--------------------------------------------------------------------------------------------------
! FFT
2012-02-09 21:28:15 +05:30
call fftw_execute_dft_r2c ( fftw_forth , field_real , field_fourier )
2012-08-25 17:16:36 +05:30
2013-03-07 13:28:42 +05:30
!--------------------------------------------------------------------------------------------------
2013-03-27 17:58:55 +05:30
! remove highest frequency in each direction, in third direction only if not 2D
2013-03-07 13:28:42 +05:30
field_fourier ( res ( 1 ) / 2_pInt + 1_pInt , 1 : res ( 2 ) , 1 : res ( 3 ) , &
1 : vec_tens , 1 : 3 ) = cmplx ( 0.0_pReal , 0.0_pReal , pReal )
field_fourier ( 1 : res1_red , res ( 2 ) / 2_pInt + 1_pInt , 1 : res ( 3 ) , &
1 : vec_tens , 1 : 3 ) = cmplx ( 0.0_pReal , 0.0_pReal , pReal )
2012-02-09 21:28:15 +05:30
if ( res ( 3 ) > 1_pInt ) &
2013-03-07 13:28:42 +05:30
field_fourier ( 1 : res1_red , 1 : res ( 2 ) , res ( 3 ) / 2_pInt + 1_pInt , &
1 : vec_tens , 1 : 3 ) = cmplx ( 0.0_pReal , 0.0_pReal , pReal )
2012-08-25 17:16:36 +05:30
2013-03-07 13:28:42 +05:30
!--------------------------------------------------------------------------------------------------
! differentiation in Fourier space
2013-03-07 01:04:30 +05:30
do k = 1_pInt , res ( 3 )
2012-01-20 02:08:52 +05:30
k_s ( 3 ) = k - 1_pInt
if ( k > res ( 3 ) / 2_pInt + 1_pInt ) k_s ( 3 ) = k_s ( 3 ) - res ( 3 )
2013-03-07 01:04:30 +05:30
do j = 1_pInt , res ( 2 )
k_s ( 2 ) = j - 1_pInt
if ( j > res ( 2 ) / 2_pInt + 1_pInt ) k_s ( 2 ) = k_s ( 2 ) - res ( 2 )
do i = 1_pInt , res1_red
k_s ( 1 ) = i - 1_pInt
xi = cmplx ( real ( k_s , pReal ) / geomdim , 0.0_pReal )
do l = 1_pInt , vec_tens
curl_fourier ( i , j , k , l , 1 ) = ( field_fourier ( i , j , k , l , 3 ) * xi ( 2 ) &
- field_fourier ( i , j , k , l , 2 ) * xi ( 3 ) ) * TWOPIIMG
curl_fourier ( i , j , k , l , 2 ) = ( - field_fourier ( i , j , k , l , 3 ) * xi ( 1 ) &
+ field_fourier ( i , j , k , l , 1 ) * xi ( 3 ) ) * TWOPIIMG
curl_fourier ( i , j , k , l , 3 ) = ( field_fourier ( i , j , k , l , 2 ) * xi ( 1 ) &
- field_fourier ( i , j , k , l , 1 ) * xi ( 2 ) ) * TWOPIIMG
enddo
2011-12-01 17:31:13 +05:30
enddo ; enddo ; enddo
2013-03-07 13:28:42 +05:30
!--------------------------------------------------------------------------------------------------
! iFFT
2012-02-09 21:28:15 +05:30
call fftw_execute_dft_c2r ( fftw_back , curl_fourier , curl_real )
2013-03-07 13:28:42 +05:30
math_curlFFT = curl_real ( 1 : res ( 1 ) , 1 : res ( 2 ) , 1 : res ( 3 ) , 1 : vec_tens , 1 : 3 ) * wgt ! copy to output and weight
if ( vec_tens == 3_pInt ) &
2013-03-27 17:58:55 +05:30
forall ( k = 1_pInt : res ( 3 ) , j = 1_pInt : res ( 2 ) , i = 1_pInt : res ( 1 ) ) &
2013-03-07 13:28:42 +05:30
math_curlFFT ( i , j , k , 1 : 3 , 1 : 3 ) = math_transpose33 ( math_curlFFT ( i , j , k , 1 : 3 , 1 : 3 ) ) ! results are stored transposed
2013-03-07 01:04:30 +05:30
2012-08-27 13:34:47 +05:30
call fftw_destroy_plan ( fftw_forth )
call fftw_destroy_plan ( fftw_back )
call fftw_free ( field_fftw )
call fftw_free ( curl_fftw )
2011-12-01 17:31:13 +05:30
2012-08-27 13:34:47 +05:30
end function math_curlFFT
2013-03-07 13:28:42 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculates gradient field using differentation in Fourier space
2013-11-22 16:01:52 +05:30
!> @todo enable odd resolution
2013-03-07 13:28:42 +05:30
!--------------------------------------------------------------------------------------------------
function math_gradFFT ( geomdim , field )
use IO , only : &
IO_error
use numerics , only : &
fftw_timelimit , &
fftw_planner_flag
use debug , only : &
debug_math , &
debug_level , &
debug_levelBasic
2011-12-01 17:31:13 +05:30
2013-03-07 13:28:42 +05:30
implicit none
real ( pReal ) , dimension ( : , : , : , : ) , intent ( in ) :: field !< field of data, first three dimensions are resolution, 4th is 1 or 3 (scalar or vector)
real ( pReal ) , dimension ( size ( field , 1 ) , size ( field , 2 ) , size ( field , 3 ) , 3 , size ( field , 4 ) ) :: &
math_gradFFT
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: geomdim !< physical length dimension in three directions
real ( pReal ) , dimension ( : , : , : , : ) , pointer :: field_real
complex ( pReal ) , dimension ( : , : , : , : ) , pointer :: field_fourier
real ( pReal ) , dimension ( : , : , : , : , : ) , pointer :: grad_real
complex ( pReal ) , dimension ( : , : , : , : , : ) , pointer :: grad_fourier
integer ( pInt ) , dimension ( 3 ) :: &
k_s , &
res
complex ( pReal ) , dimension ( 3 ) :: xi
type ( C_PTR ) :: fftw_forth , fftw_back
type ( C_PTR ) :: field_fftw , grad_fftw
integer ( pInt ) :: i , j , k , l , res1_red , vec_tens
real ( pReal ) :: wgt
res = [ size ( field , 1 ) , size ( field , 2 ) , size ( field , 3 ) ]
vec_tens = size ( field , 4 )
if ( iand ( debug_level ( debug_math ) , debug_levelBasic ) / = 0_pInt ) then
if ( vec_tens == 1_pInt ) write ( 6 , '(a)' ) 'Calculating gradient of scalar field'
if ( vec_tens == 3_pInt ) write ( 6 , '(a)' ) 'Calculating gradeint of vector field'
write ( 6 , '(a,3(e12.5))' ) ' Dimension: ' , geomdim
write ( 6 , '(a,3(i5))' ) ' Resolution:' , res
endif
!--------------------------------------------------------------------------------------------------
! sanity checks
if ( vec_tens / = 1_pInt . and . vec_tens / = 3_pInt ) &
call IO_error ( 0_pInt , ext_msg = 'Invalid data type in math_gradFFT' )
if ( ( mod ( res ( 3 ) , 2_pInt ) / = 0_pInt . and . res ( 3 ) / = 1_pInt ) . or . &
mod ( res ( 2 ) , 2_pInt ) / = 0_pInt . or . &
mod ( res ( 1 ) , 2_pInt ) / = 0_pInt ) &
call IO_error ( 0_pInt , ext_msg = 'Resolution in math_gradFFT' )
if ( pReal / = C_DOUBLE . or . pInt / = C_INT ) &
call IO_error ( 0_pInt , ext_msg = 'Fortran to C in math_gradFFT' )
wgt = 1.0_pReal / real ( product ( res ) , pReal )
res1_red = res ( 1 ) / 2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c)
!--------------------------------------------------------------------------------------------------
! allocation and FFTW initialization
call fftw_set_timelimit ( fftw_timelimit )
field_fftw = fftw_alloc_complex ( int ( res1_red * res ( 2 ) * res ( 3 ) * vec_tens , C_SIZE_T ) ) ! C_SIZE_T is of type integer(8)
call c_f_pointer ( field_fftw , field_real , [ res ( 1 ) + 2_pInt , res ( 2 ) , res ( 3 ) , vec_tens ] )
call c_f_pointer ( field_fftw , field_fourier , [ res1_red , res ( 2 ) , res ( 3 ) , vec_tens ] )
grad_fftw = fftw_alloc_complex ( int ( res1_red * res ( 2 ) * res ( 3 ) * vec_tens * 3_pInt , C_SIZE_T ) ) ! C_SIZE_T is of type integer(8)
call c_f_pointer ( grad_fftw , grad_real , [ res ( 1 ) + 2_pInt , res ( 2 ) , res ( 3 ) , 3_pInt , vec_tens ] )
call c_f_pointer ( grad_fftw , grad_fourier , [ res1_red , res ( 2 ) , res ( 3 ) , 3_pInt , vec_tens ] )
fftw_forth = fftw_plan_many_dft_r2c ( 3_pInt , [ res ( 3 ) , res ( 2 ) , res ( 1 ) ] , vec_tens * 3_pInt , & ! dimensions, length in each dimension in reversed order, total # of transforms
field_real , [ res ( 3 ) , res ( 2 ) , res ( 1 ) + 2_pInt ] , & ! input data, physical length in each dimension in reversed order
1_pInt , res ( 3 ) * res ( 2 ) * ( res ( 1 ) + 2_pInt ) , & ! striding, product of physical lenght in the 3 dimensions
field_fourier , [ res ( 3 ) , res ( 2 ) , res1_red ] , & ! output data, physical length in each dimension in reversed order
1_pInt , res ( 3 ) * res ( 2 ) * res1_red , fftw_planner_flag ) ! striding, product of physical lenght in the 3 dimensions, planner mode
fftw_back = fftw_plan_many_dft_c2r ( 3_pInt , [ res ( 3 ) , res ( 2 ) , res ( 1 ) ] , vec_tens * 3_pInt , &
grad_fourier , [ res ( 3 ) , res ( 2 ) , res1_red ] , &
1_pInt , res ( 3 ) * res ( 2 ) * res1_red , &
grad_real , [ res ( 3 ) , res ( 2 ) , res ( 1 ) + 2_pInt ] , &
1_pInt , res ( 3 ) * res ( 2 ) * ( res ( 1 ) + 2_pInt ) , fftw_planner_flag )
field_real ( 1 : res ( 1 ) , 1 : res ( 2 ) , 1 : res ( 3 ) , 1 : vec_tens ) = field ! field_real is overwritten during plan creation and is larger (padding)
!--------------------------------------------------------------------------------------------------
! FFT
call fftw_execute_dft_r2c ( fftw_forth , field_real , field_fourier )
!--------------------------------------------------------------------------------------------------
2013-03-27 17:58:55 +05:30
! remove highest frequency in each direction, in third direction only if not 2D
2013-03-07 13:28:42 +05:30
field_fourier ( res ( 1 ) / 2_pInt + 1_pInt , 1 : res ( 2 ) , 1 : res ( 3 ) , &
1 : vec_tens ) = cmplx ( 0.0_pReal , 0.0_pReal , pReal )
field_fourier ( 1 : res1_red , res ( 2 ) / 2_pInt + 1_pInt , 1 : res ( 3 ) , &
1 : vec_tens ) = cmplx ( 0.0_pReal , 0.0_pReal , pReal )
if ( res ( 3 ) > 1_pInt ) &
field_fourier ( 1 : res1_red , 1 : res ( 2 ) , res ( 3 ) / 2_pInt + 1_pInt , &
1 : vec_tens ) = cmplx ( 0.0_pReal , 0.0_pReal , pReal )
!--------------------------------------------------------------------------------------------------
! differentiation in Fourier space
do k = 1_pInt , res ( 3 )
k_s ( 3 ) = k - 1_pInt
if ( k > res ( 3 ) / 2_pInt + 1_pInt ) k_s ( 3 ) = k_s ( 3 ) - res ( 3 )
do j = 1_pInt , res ( 2 )
k_s ( 2 ) = j - 1_pInt
if ( j > res ( 2 ) / 2_pInt + 1_pInt ) k_s ( 2 ) = k_s ( 2 ) - res ( 2 )
do i = 1_pInt , res1_red
k_s ( 1 ) = i - 1_pInt
xi = cmplx ( real ( k_s , pReal ) / geomdim , 0.0_pReal )
do l = 1_pInt , vec_tens
grad_fourier ( i , j , k , 1 , l ) = field_fourier ( i , j , k , l ) * xi ( 1 ) * TWOPIIMG
grad_fourier ( i , j , k , 2 , l ) = field_fourier ( i , j , k , l ) * xi ( 2 ) * TWOPIIMG
grad_fourier ( i , j , k , 3 , l ) = field_fourier ( i , j , k , l ) * xi ( 3 ) * TWOPIIMG
enddo
enddo ; enddo ; enddo
!--------------------------------------------------------------------------------------------------
! iFFT
call fftw_execute_dft_c2r ( fftw_back , grad_fourier , grad_real )
math_gradFFT = grad_real ( 1 : res ( 1 ) , 1 : res ( 2 ) , 1 : res ( 3 ) , 1 : 3 , 1 : vec_tens ) * wgt ! copy to output and weight
if ( vec_tens == 3_pInt ) &
forall ( k = 1_pInt : res ( 3 ) , j = 1_pInt : res ( 2 ) , i = 1_pInt : res ( 1 ) ) &
math_gradFFT ( i , j , k , 1 : 3 , 1 : 3 ) = math_transpose33 ( math_gradFFT ( i , j , k , 1 : 3 , 1 : 3 ) ) ! results are stored transposed
call fftw_destroy_plan ( fftw_forth )
call fftw_destroy_plan ( fftw_back )
call fftw_free ( field_fftw )
call fftw_free ( grad_fftw )
end function math_gradFFT
2013-03-07 01:04:30 +05:30
2012-08-27 13:34:47 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculates divergence field using integration in Fourier space
2013-11-22 16:01:52 +05:30
!> @todo enable odd resolution
2012-08-27 13:34:47 +05:30
!--------------------------------------------------------------------------------------------------
function math_divergenceFFT ( geomdim , field )
2013-03-07 01:04:30 +05:30
use IO , only : &
IO_error
use numerics , only : &
fftw_timelimit , &
fftw_planner_flag
use debug , only : &
debug_math , &
debug_level , &
debug_levelBasic
2012-03-09 01:55:28 +05:30
2011-12-01 17:31:13 +05:30
implicit none
2013-03-07 13:28:42 +05:30
real ( pReal ) , dimension ( : , : , : , : , : ) , intent ( in ) :: field !< field of data, first three dimensions are resolution, 4th is 1 or 3 (vector or tensor), 5th is 3
real ( pReal ) , dimension ( size ( field , 1 ) , size ( field , 2 ) , size ( field , 3 ) , size ( field , 4 ) ) :: &
math_divergenceFFT
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: geomdim !< physical length dimension in three directions
2012-01-13 21:48:16 +05:30
real ( pReal ) , dimension ( : , : , : , : , : ) , pointer :: field_real
2012-02-09 21:28:15 +05:30
complex ( pReal ) , dimension ( : , : , : , : , : ) , pointer :: field_fourier
2012-01-13 21:48:16 +05:30
real ( pReal ) , dimension ( : , : , : , : ) , pointer :: divergence_real
2012-02-09 21:28:15 +05:30
complex ( pReal ) , dimension ( : , : , : , : ) , pointer :: divergence_fourier
2013-03-07 13:28:42 +05:30
integer ( pInt ) , dimension ( 3 ) :: &
k_s , &
res
complex ( pReal ) , dimension ( 3 ) :: &
xi
type ( C_PTR ) :: fftw_forth , fftw_back
type ( C_PTR ) :: field_fftw , divergence_fftw
integer ( pInt ) :: i , j , k , l , res1_red , vec_tens
real ( pReal ) :: wgt
2013-03-07 01:04:30 +05:30
2012-08-27 13:34:47 +05:30
res = [ size ( field , 1 ) , size ( field , 2 ) , size ( field , 3 ) ]
vec_tens = size ( field , 4 )
2011-12-01 17:31:13 +05:30
2012-07-05 15:24:50 +05:30
if ( iand ( debug_level ( debug_math ) , debug_levelBasic ) / = 0_pInt ) then
2013-01-31 21:58:08 +05:30
if ( vec_tens == 1_pInt ) write ( 6 , '(a)' ) 'Calculating FFT divergence of vector field'
if ( vec_tens == 3_pInt ) write ( 6 , '(a)' ) 'Calculating FFT divergence of tensor field'
write ( 6 , '(a,3(e12.5))' ) ' Dimension: ' , geomdim
write ( 6 , '(a,3(i5))' ) ' Resolution:' , res
2012-01-13 21:48:16 +05:30
endif
2012-12-14 23:00:22 +05:30
2013-03-07 13:28:42 +05:30
!--------------------------------------------------------------------------------------------------
! sanity checks
2012-12-14 23:00:22 +05:30
if ( vec_tens / = 1_pInt . and . vec_tens / = 3_pInt ) &
call IO_error ( 0_pInt , ext_msg = 'Invalid data type in math_divergenceFFT' )
if ( ( mod ( res ( 3 ) , 2_pInt ) / = 0_pInt . and . res ( 3 ) / = 1_pInt ) . or . &
mod ( res ( 2 ) , 2_pInt ) / = 0_pInt . or . &
mod ( res ( 1 ) , 2_pInt ) / = 0_pInt ) &
call IO_error ( 0_pInt , ext_msg = 'Resolution in math_divergenceFFT' )
if ( pReal / = C_DOUBLE . or . pInt / = C_INT ) &
call IO_error ( 0_pInt , ext_msg = 'Fortran to C in math_divergenceFFT' )
2013-03-07 01:04:30 +05:30
2012-02-13 23:11:27 +05:30
res1_red = res ( 1 ) / 2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c)
2013-03-07 13:28:42 +05:30
wgt = 1.0_pReal / real ( product ( res ) , pReal )
2012-01-25 16:00:39 +05:30
2013-03-07 13:28:42 +05:30
!--------------------------------------------------------------------------------------------------
! allocation and FFTW initialization
2012-01-13 21:48:16 +05:30
call fftw_set_timelimit ( fftw_timelimit )
2013-03-07 13:28:42 +05:30
field_fftw = fftw_alloc_complex ( int ( res1_red * res ( 2 ) * res ( 3 ) * vec_tens * 3_pInt , C_SIZE_T ) ) ! C_SIZE_T is of type integer(8)
2012-02-09 21:28:15 +05:30
call c_f_pointer ( field_fftw , field_real , [ res ( 1 ) + 2_pInt , res ( 2 ) , res ( 3 ) , vec_tens , 3_pInt ] )
call c_f_pointer ( field_fftw , field_fourier , [ res1_red , res ( 2 ) , res ( 3 ) , vec_tens , 3_pInt ] )
2013-03-07 13:28:42 +05:30
divergence_fftw = fftw_alloc_complex ( int ( res1_red * res ( 2 ) * res ( 3 ) * vec_tens , C_SIZE_T ) ) ! C_SIZE_T is of type integer(8)
2012-01-13 21:48:16 +05:30
call c_f_pointer ( divergence_fftw , divergence_real , [ res ( 1 ) + 2_pInt , res ( 2 ) , res ( 3 ) , vec_tens ] )
2012-02-09 21:28:15 +05:30
call c_f_pointer ( divergence_fftw , divergence_fourier , [ res1_red , res ( 2 ) , res ( 3 ) , vec_tens ] )
2012-01-13 21:48:16 +05:30
2013-03-07 13:28:42 +05:30
fftw_forth = fftw_plan_many_dft_r2c ( 3_pInt , [ res ( 3 ) , res ( 2 ) , res ( 1 ) ] , vec_tens * 3_pInt , & ! dimensions, length in each dimension in reversed order, total # of transforms
field_real , [ res ( 3 ) , res ( 2 ) , res ( 1 ) + 2_pInt ] , & ! input data, physical length in each dimension in reversed order
1_pInt , res ( 3 ) * res ( 2 ) * ( res ( 1 ) + 2_pInt ) , & ! striding, product of physical lenght in the 3 dimensions
field_fourier , [ res ( 3 ) , res ( 2 ) , res1_red ] , & ! output data, physical length in each dimension in reversed order
1_pInt , res ( 3 ) * res ( 2 ) * res1_red , fftw_planner_flag ) ! striding, product of physical lenght in the 3 dimensions, planner mode
2013-03-07 01:04:30 +05:30
fftw_back = fftw_plan_many_dft_c2r ( 3_pInt , [ res ( 3 ) , res ( 2 ) , res ( 1 ) ] , vec_tens , &
divergence_fourier , [ res ( 3 ) , res ( 2 ) , res1_red ] , &
2012-02-09 21:28:15 +05:30
1_pInt , res ( 3 ) * res ( 2 ) * res1_red , &
2013-03-07 01:04:30 +05:30
divergence_real , [ res ( 3 ) , res ( 2 ) , res ( 1 ) + 2_pInt ] , &
2013-03-07 13:28:42 +05:30
1_pInt , res ( 3 ) * res ( 2 ) * ( res ( 1 ) + 2_pInt ) , fftw_planner_flag )
2013-03-07 01:04:30 +05:30
2013-03-07 13:28:42 +05:30
field_real ( 1 : res ( 1 ) , 1 : res ( 2 ) , 1 : res ( 3 ) , 1 : vec_tens , 1 : 3 ) = field ! field_real is overwritten during plan creation and is larger (padding)
2012-08-25 17:16:36 +05:30
2013-03-07 13:28:42 +05:30
!--------------------------------------------------------------------------------------------------
! FFT
2012-02-09 21:28:15 +05:30
call fftw_execute_dft_r2c ( fftw_forth , field_real , field_fourier )
2012-08-25 17:16:36 +05:30
2013-03-07 13:28:42 +05:30
!--------------------------------------------------------------------------------------------------
2013-03-27 17:58:55 +05:30
! remove highest frequency in each direction, in third direction only if not 2D
2013-03-07 13:28:42 +05:30
field_fourier ( res ( 1 ) / 2_pInt + 1_pInt , 1 : res ( 2 ) , 1 : res ( 3 ) , &
1 : vec_tens , 1 : 3 ) = cmplx ( 0.0_pReal , 0.0_pReal , pReal )
field_fourier ( 1 : res1_red , res ( 2 ) / 2_pInt + 1_pInt , 1 : res ( 3 ) , &
1 : vec_tens , 1 : 3 ) = cmplx ( 0.0_pReal , 0.0_pReal , pReal )
2012-02-09 21:28:15 +05:30
if ( res ( 3 ) > 1_pInt ) &
2013-03-07 13:28:42 +05:30
field_fourier ( 1 : res1_red , 1 : res ( 2 ) , res ( 3 ) / 2_pInt + 1_pInt , &
1 : vec_tens , 1 : 3 ) = cmplx ( 0.0_pReal , 0.0_pReal , pReal )
2012-08-25 17:16:36 +05:30
2013-03-07 13:28:42 +05:30
!--------------------------------------------------------------------------------------------------
! differentiation in Fourier space
2013-03-07 01:04:30 +05:30
do k = 1_pInt , res ( 3 )
k_s ( 3 ) = k - 1_pInt
if ( k > res ( 3 ) / 2_pInt + 1_pInt ) k_s ( 3 ) = k_s ( 3 ) - res ( 3 )
do j = 1_pInt , res ( 2 )
k_s ( 2 ) = j - 1_pInt
if ( j > res ( 2 ) / 2_pInt + 1_pInt ) k_s ( 2 ) = k_s ( 2 ) - res ( 2 )
do i = 1_pInt , res1_red
k_s ( 1 ) = i - 1_pInt
xi = cmplx ( real ( k_s , pReal ) / geomdim , 0.0_pReal )
do l = 1_pInt , vec_tens
divergence_fourier ( i , j , k , l ) = sum ( field_fourier ( i , j , k , l , 1 : 3 ) * xi ) * TWOPIIMG
enddo
2011-12-01 17:31:13 +05:30
enddo ; enddo ; enddo
2013-03-07 01:04:30 +05:30
2013-03-07 13:28:42 +05:30
!--------------------------------------------------------------------------------------------------
! iFFT
2012-02-09 21:28:15 +05:30
call fftw_execute_dft_c2r ( fftw_back , divergence_fourier , divergence_real )
2013-03-07 13:28:42 +05:30
math_divergenceFFT = divergence_real ( 1 : res ( 1 ) , 1 : res ( 2 ) , 1 : res ( 3 ) , 1 : vec_tens ) * wgt ! copy to output and weight
2011-12-01 17:31:13 +05:30
2012-08-27 13:34:47 +05:30
call fftw_destroy_plan ( fftw_forth )
call fftw_destroy_plan ( fftw_back )
call fftw_free ( field_fftw )
call fftw_free ( divergence_fftw )
2011-12-01 17:31:13 +05:30
2012-08-27 13:34:47 +05:30
end function math_divergenceFFT
2012-08-24 18:57:55 +05:30
2011-12-01 17:31:13 +05:30
2012-08-27 13:34:47 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculates divergence field using FDM with variable accuracy
!--------------------------------------------------------------------------------------------------
function math_divergenceFDM ( geomdim , order , field )
2013-03-07 01:04:30 +05:30
use IO , only : &
IO_error
use debug , only : &
debug_math , &
debug_level , &
debug_levelBasic
2011-12-01 17:31:13 +05:30
implicit none
2013-03-07 13:28:42 +05:30
real ( pReal ) , dimension ( : , : , : , : , : ) , intent ( in ) :: field !< field of data, first three dimensions are resolution, 4th is 1 or 3 (vector or tensor), 5th is 3
real ( pReal ) , dimension ( size ( field , 1 ) , size ( field , 2 ) , size ( field , 3 ) , size ( field , 4 ) ) :: &
math_divergenceFDM
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: geomdim !< physical length dimension in three directions
integer ( pInt ) , intent ( in ) :: order !< order of Finite Differences
real ( pReal ) , dimension ( 4 , 4 ) , parameter :: FDcoefficient = reshape ( [ &
1.0_pReal / 2.0_pReal , 0.0_pReal , 0.0_pReal , 0.0_pReal , & ! from http://en.wikipedia.org/wiki/Finite_difference_coefficients
2.0_pReal / 3.0_pReal , - 1.0_pReal / 1 2.0_pReal , 0.0_pReal , 0.0_pReal , &
3.0_pReal / 4.0_pReal , - 3.0_pReal / 2 0.0_pReal , 1.0_pReal / 6 0.0_pReal , 0.0_pReal , &
4.0_pReal / 5.0_pReal , - 1.0_pReal / 5.0_pReal , 4.0_pReal / 10 5.0_pReal , - 1.0_pReal / 28 0.0_pReal ] , [ 4 , 4 ] )
integer ( pInt ) , dimension ( 6 , 3 ) :: coordinates
integer ( pInt ) , dimension ( 3 ) :: res
integer ( pInt ) :: i , j , k , m , l , vec_tens
2012-08-27 13:34:47 +05:30
res = [ size ( field , 1 ) , size ( field , 2 ) , size ( field , 3 ) ]
vec_tens = size ( field , 4 )
2011-12-01 17:31:13 +05:30
2012-07-05 15:24:50 +05:30
if ( iand ( debug_level ( debug_math ) , debug_levelBasic ) / = 0_pInt ) then
2013-01-31 21:58:08 +05:30
if ( vec_tens == 1_pInt ) write ( 6 , '(a)' ) 'Calculating FDM divergence of vector field'
if ( vec_tens == 3_pInt ) write ( 6 , '(a)' ) 'Calculating FDM divergence of tensor field'
write ( 6 , '(a,3(e12.5))' ) ' Dimension: ' , geomdim
write ( 6 , '(a,3(i5))' ) ' Resolution:' , res
2012-01-13 21:48:16 +05:30
endif
2011-12-01 17:31:13 +05:30
2013-03-07 13:28:42 +05:30
!--------------------------------------------------------------------------------------------------
! sanity checks
if ( vec_tens / = 1_pInt . and . vec_tens / = 3_pInt ) &
call IO_error ( 0_pInt , ext_msg = 'Invalid data type in math_divergenceFDM' )
2012-12-14 23:00:22 +05:30
2013-03-07 13:28:42 +05:30
!--------------------------------------------------------------------------------------------------
! differentiation in real space
2012-08-27 13:34:47 +05:30
math_divergenceFDM = 0.0_pReal
2011-12-01 17:31:13 +05:30
do k = 0_pInt , res ( 3 ) - 1_pInt ; do j = 0_pInt , res ( 2 ) - 1_pInt ; do i = 0_pInt , res ( 1 ) - 1_pInt
2012-08-27 13:34:47 +05:30
do m = 1_pInt , order + 1_pInt
2013-03-07 13:28:42 +05:30
coordinates ( 1 , 1 : 3 ) = periodic_location ( periodic_index ( [ i + m , j , k ] , res ) , res )
coordinates ( 2 , 1 : 3 ) = periodic_location ( periodic_index ( [ i - m , j , k ] , res ) , res )
coordinates ( 3 , 1 : 3 ) = periodic_location ( periodic_index ( [ i , j + m , k ] , res ) , res )
coordinates ( 4 , 1 : 3 ) = periodic_location ( periodic_index ( [ i , j - m , k ] , res ) , res )
coordinates ( 5 , 1 : 3 ) = periodic_location ( periodic_index ( [ i , j , k + m ] , res ) , res )
coordinates ( 6 , 1 : 3 ) = periodic_location ( periodic_index ( [ i , j , k - m ] , res ) , res )
2013-01-31 21:58:08 +05:30
coordinates = coordinates + 1_pInt
2011-12-01 17:31:13 +05:30
do l = 1_pInt , vec_tens
2012-08-27 13:34:47 +05:30
math_divergenceFDM ( i + 1_pInt , j + 1_pInt , k + 1_pInt , l ) = math_divergenceFDM ( i + 1_pInt , j + 1_pInt , k + 1_pInt , l ) &
2012-08-27 21:13:17 +05:30
+ FDcoefficient ( m , order + 1_pInt ) * &
2011-12-01 17:31:13 +05:30
( ( field ( coordinates ( 1 , 1 ) , coordinates ( 1 , 2 ) , coordinates ( 1 , 3 ) , l , 1 ) - &
field ( coordinates ( 2 , 1 ) , coordinates ( 2 , 2 ) , coordinates ( 2 , 3 ) , l , 1 ) ) * real ( res ( 1 ) , pReal ) / geomdim ( 1 ) + &
( field ( coordinates ( 3 , 1 ) , coordinates ( 3 , 2 ) , coordinates ( 3 , 3 ) , l , 2 ) - &
field ( coordinates ( 4 , 1 ) , coordinates ( 4 , 2 ) , coordinates ( 4 , 3 ) , l , 2 ) ) * real ( res ( 2 ) , pReal ) / geomdim ( 2 ) + &
( field ( coordinates ( 5 , 1 ) , coordinates ( 5 , 2 ) , coordinates ( 5 , 3 ) , l , 3 ) - &
field ( coordinates ( 6 , 1 ) , coordinates ( 6 , 2 ) , coordinates ( 6 , 3 ) , l , 3 ) ) * real ( res ( 3 ) , pReal ) / geomdim ( 3 ) )
enddo
enddo
enddo ; enddo ; enddo
2013-03-07 01:04:30 +05:30
2013-01-31 21:58:08 +05:30
contains
!--------------------------------------------------------------------------------------------------
2013-11-22 16:01:52 +05:30
!> @brief small helper functions for indexing. CAREFUL: index and location runs from
2013-01-31 21:58:08 +05:30
! 0 to N-1 (python style)
!--------------------------------------------------------------------------------------------------
pure function periodic_location ( idx , res )
implicit none
integer ( pInt ) , intent ( in ) :: idx
integer ( pInt ) , intent ( in ) , dimension ( 3 ) :: res
integer ( pInt ) , dimension ( 3 ) :: periodic_location
periodic_location = [ modulo ( idx / res ( 3 ) / res ( 2 ) , res ( 1 ) ) , &
modulo ( idx / res ( 3 ) , res ( 2 ) ) , &
modulo ( idx , res ( 3 ) ) ]
end function periodic_location
!--------------------------------------------------------------------------------------------------
2013-11-22 16:01:52 +05:30
!> @brief small helper functions for indexing CAREFUL: index and location runs from
2013-01-31 21:58:08 +05:30
! 0 to N-1 (python style)
!--------------------------------------------------------------------------------------------------
integer ( pInt ) pure function periodic_index ( location , res )
2013-03-07 01:04:30 +05:30
2013-01-31 21:58:08 +05:30
implicit none
integer ( pInt ) , intent ( in ) , dimension ( 3 ) :: res , location
periodic_index = modulo ( location ( 3 ) , res ( 3 ) ) + &
( modulo ( location ( 2 ) , res ( 2 ) ) ) * res ( 3 ) + &
( modulo ( location ( 1 ) , res ( 1 ) ) ) * res ( 3 ) * res ( 2 )
2011-12-01 17:31:13 +05:30
2013-01-31 21:58:08 +05:30
end function periodic_index
2012-08-27 13:34:47 +05:30
end function math_divergenceFDM
!--------------------------------------------------------------------------------------------------
2013-03-27 15:44:44 +05:30
!> @brief Obtain the nearest neighbor from periodic domainSet at points in querySet
2012-08-27 13:34:47 +05:30
!--------------------------------------------------------------------------------------------------
2013-01-31 21:58:08 +05:30
function math_periodicNearestNeighbor ( geomdim , Favg , querySet , domainSet )
2012-08-27 13:34:47 +05:30
use kdtree2_module
use IO , only : &
IO_error
2013-03-07 01:04:30 +05:30
2012-08-27 13:34:47 +05:30
implicit none
2013-01-31 21:58:08 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: Favg
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: geomdim
real ( pReal ) , dimension ( : , : ) , intent ( in ) :: querySet
real ( pReal ) , dimension ( : , : ) , intent ( in ) :: domainSet
integer ( pInt ) , dimension ( size ( querySet , 2 ) ) :: math_periodicNearestNeighbor
2013-03-07 01:04:30 +05:30
real ( pReal ) , dimension ( size ( domainSet , 1 ) , ( 3_pInt ** size ( domainSet , 1 ) ) * size ( domainSet , 2 ) ) :: &
domainSetLarge
2013-01-31 21:58:08 +05:30
integer ( pInt ) :: i , j , l , m , n , spatialDim
2012-08-27 13:34:47 +05:30
type ( kdtree2 ) , pointer :: tree
type ( kdtree2_result ) , dimension ( 1 ) :: Results
2013-01-31 21:58:08 +05:30
if ( size ( querySet , 1 ) / = size ( domainSet , 1 ) ) call IO_error ( 407_pInt , ext_msg = 'query set' )
spatialDim = size ( querySet , 1 )
2012-08-24 18:57:55 +05:30
2012-08-27 13:34:47 +05:30
i = 0_pInt
if ( spatialDim == 2_pInt ) then
2013-01-31 21:58:08 +05:30
do j = 1_pInt , size ( domainSet , 2 )
2012-08-27 13:34:47 +05:30
do l = - 1_pInt , 1_pInt ; do m = - 1_pInt , 1_pInt
i = i + 1_pInt
domainSetLarge ( 1 : 2 , i ) = domainSet ( 1 : 2 , j ) + matmul ( Favg ( 1 : 2 , 1 : 2 ) , real ( [ l , m ] , pReal ) * geomdim ( 1 : 2 ) )
enddo ; enddo
enddo
else
2013-01-31 21:58:08 +05:30
do j = 1_pInt , size ( domainSet , 2 )
2012-08-27 13:34:47 +05:30
do l = - 1_pInt , 1_pInt ; do m = - 1_pInt , 1_pInt ; do n = - 1_pInt , 1_pInt
i = i + 1_pInt
domainSetLarge ( 1 : 3 , i ) = domainSet ( 1 : 3 , j ) + math_mul33x3 ( Favg , real ( [ l , m , n ] , pReal ) * geomdim )
enddo ; enddo ; enddo
enddo
endif
tree = > kdtree2_create ( domainSetLarge , sort = . true . , rearrange = . true . )
2013-01-31 21:58:08 +05:30
do j = 1_pInt , size ( querySet , 2 )
2012-08-27 13:34:47 +05:30
call kdtree2_n_nearest ( tp = tree , qv = querySet ( 1 : spatialDim , j ) , nn = 1_pInt , results = Results )
2013-01-31 21:58:08 +05:30
math_periodicNearestNeighbor ( j ) = Results ( 1 ) % idx
2012-08-27 13:34:47 +05:30
enddo
2013-01-31 21:58:08 +05:30
math_periodicNearestNeighbor = math_periodicNearestNeighbor - 1_pInt ! let them run from 0 to domainPoints -1
end function math_periodicNearestNeighbor
2013-02-19 20:26:26 +05:30
2013-06-07 12:57:12 +05:30
2013-03-27 15:44:44 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Obtain the nearest neighbor from domainSet at points in querySet
!--------------------------------------------------------------------------------------------------
function math_nearestNeighbor ( querySet , domainSet )
use kdtree2_module
use IO , only : &
IO_error
implicit none
real ( pReal ) , dimension ( : , : ) , intent ( in ) :: querySet
real ( pReal ) , dimension ( : , : ) , intent ( in ) :: domainSet
integer ( pInt ) , dimension ( size ( querySet , 2 ) ) :: math_nearestNeighbor
2013-03-27 17:58:55 +05:30
integer ( pInt ) :: j , spatialDim
2013-03-27 15:44:44 +05:30
type ( kdtree2 ) , pointer :: tree
type ( kdtree2_result ) , dimension ( 1 ) :: Results
if ( size ( querySet , 1 ) / = size ( domainSet , 1 ) ) call IO_error ( 407_pInt , ext_msg = 'query set' )
spatialDim = size ( querySet , 1 )
tree = > kdtree2_create ( domainSet , sort = . true . , rearrange = . true . )
do j = 1_pInt , size ( querySet , 2 )
call kdtree2_n_nearest ( tp = tree , qv = querySet ( 1 : spatialDim , j ) , nn = 1_pInt , results = Results )
math_nearestNeighbor ( j ) = Results ( 1 ) % idx
enddo
math_nearestNeighbor = math_nearestNeighbor - 1_pInt ! let them run from 0 to domainPoints -1
end function math_nearestNeighbor
2013-02-19 20:26:26 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Obtain the distances to the next N nearest neighbors from domainSet at points in querySet
!--------------------------------------------------------------------------------------------------
function math_periodicNearestNeighborDistances ( geomdim , Favg , querySet , domainSet , Ndist ) result ( distances )
use kdtree2_module
use IO , only : &
IO_error
2014-02-28 18:58:27 +05:30
2013-02-19 20:26:26 +05:30
implicit none
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: geomdim
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: Favg
integer ( pInt ) , intent ( in ) :: Ndist
real ( pReal ) , dimension ( : , : ) , intent ( in ) :: querySet
real ( pReal ) , dimension ( : , : ) , intent ( in ) :: domainSet
! output variable
real ( pReal ) , dimension ( Ndist , size ( querySet , 2 ) ) :: distances
real ( pReal ) , dimension ( size ( domainSet , 1 ) , ( 3_pInt ** size ( domainSet , 1 ) ) * size ( domainSet , 2 ) ) &
:: domainSetLarge
integer ( pInt ) :: i , j , l , m , n , spatialDim
type ( kdtree2 ) , pointer :: tree
type ( kdtree2_result ) , dimension ( : ) , allocatable :: Results
allocate ( Results ( Ndist ) )
if ( size ( querySet , 1 ) / = size ( domainSet , 1 ) ) call IO_error ( 407_pInt , ext_msg = 'query set' )
spatialDim = size ( querySet , 1 )
i = 0_pInt
if ( spatialDim == 2_pInt ) then
do j = 1_pInt , size ( domainSet , 2 )
do l = - 1_pInt , 1_pInt ; do m = - 1_pInt , 1_pInt
i = i + 1_pInt
domainSetLarge ( 1 : 2 , i ) = domainSet ( 1 : 2 , j ) + matmul ( Favg ( 1 : 2 , 1 : 2 ) , real ( [ l , m ] , pReal ) * geomdim ( 1 : 2 ) )
enddo ; enddo
enddo
else
do j = 1_pInt , size ( domainSet , 2 )
do l = - 1_pInt , 1_pInt ; do m = - 1_pInt , 1_pInt ; do n = - 1_pInt , 1_pInt
i = i + 1_pInt
domainSetLarge ( 1 : 3 , i ) = domainSet ( 1 : 3 , j ) + math_mul33x3 ( Favg , real ( [ l , m , n ] , pReal ) * geomdim )
enddo ; enddo ; enddo
enddo
endif
tree = > kdtree2_create ( domainSetLarge , sort = . true . , rearrange = . true . )
do j = 1_pInt , size ( querySet , 2 )
call kdtree2_n_nearest ( tp = tree , qv = querySet ( 1 : spatialDim , j ) , nn = Ndist , results = Results )
distances ( 1 : Ndist , j ) = sqrt ( Results ( 1 : Ndist ) % dis )
enddo
deallocate ( Results )
end function math_periodicNearestNeighborDistances
2012-08-27 13:34:47 +05:30
#endif
!--------------------------------------------------------------------------------------------------
!> @brief calculate average of tensor field
!--------------------------------------------------------------------------------------------------
function math_tensorAvg ( field )
2013-03-07 01:04:30 +05:30
2011-12-01 17:31:13 +05:30
implicit none
2012-08-27 13:34:47 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: math_tensorAvg
2013-03-07 01:04:30 +05:30
real ( pReal ) , intent ( in ) , dimension ( : , : , : , : , : ) :: field
2012-08-27 13:34:47 +05:30
real ( pReal ) :: wgt
2011-12-01 17:31:13 +05:30
2013-01-31 21:58:08 +05:30
wgt = 1.0_pReal / real ( size ( field , 3 ) * size ( field , 4 ) * size ( field , 5 ) , pReal )
math_tensorAvg = sum ( sum ( sum ( field , dim = 5 ) , dim = 4 ) , dim = 3 ) * wgt
2011-12-01 17:31:13 +05:30
2012-08-27 13:34:47 +05:30
end function math_tensorAvg
2011-12-01 17:31:13 +05:30
2013-01-31 21:58:08 +05:30
2012-08-27 13:34:47 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculate logarithmic strain in spatial configuration for given F field
!--------------------------------------------------------------------------------------------------
function math_logstrainSpat ( F )
2013-03-07 01:04:30 +05:30
2011-12-01 17:31:13 +05:30
implicit none
2012-08-27 13:34:47 +05:30
real ( pReal ) , intent ( in ) , dimension ( : , : , : , : , : ) :: F
2013-01-31 21:58:08 +05:30
real ( pReal ) , dimension ( 3 , 3 , size ( F , 3 ) , size ( F , 4 ) , size ( F , 5 ) ) :: math_logstrainSpat
integer ( pInt ) , dimension ( 3 ) :: res
real ( pReal ) , dimension ( 3 , 3 ) :: temp33_Real , temp33_Real2
real ( pReal ) , dimension ( 3 , 3 , 3 ) :: evbasis
real ( pReal ) , dimension ( 3 ) :: eigenvalue
2011-12-01 17:31:13 +05:30
integer ( pInt ) :: i , j , k
logical :: errmatinv
2013-01-31 21:58:08 +05:30
res = [ size ( F , 3 ) , size ( F , 4 ) , size ( F , 5 ) ]
2011-12-01 17:31:13 +05:30
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res ( 1 )
2013-01-31 21:58:08 +05:30
call math_pDecomposition ( F ( 1 : 3 , 1 : 3 , i , j , k ) , temp33_Real2 , temp33_Real , errmatinv ) !store R in temp33_Real
2012-01-26 19:20:00 +05:30
temp33_Real2 = math_inv33 ( temp33_Real )
2013-01-31 21:58:08 +05:30
temp33_Real = math_mul33x33 ( F ( 1 : 3 , 1 : 3 , i , j , k ) , temp33_Real2 ) ! v = F o inv(R), store in temp33_Real2
call math_spectralDecomposition ( temp33_Real , eigenvalue ( 1 ) , eigenvalue ( 2 ) , eigenvalue ( 3 ) , &
evbasis ( 1 : 3 , 1 : 3 , 1 ) , evbasis ( 1 : 3 , 1 : 3 , 2 ) , evbasis ( 1 : 3 , 1 : 3 , 3 ) )
2011-12-01 17:31:13 +05:30
eigenvalue = log ( sqrt ( eigenvalue ) )
2013-01-31 21:58:08 +05:30
math_logstrainSpat ( 1 : 3 , 1 : 3 , i , j , k ) = eigenvalue ( 1 ) * evbasis ( 1 : 3 , 1 : 3 , 1 ) + &
eigenvalue ( 2 ) * evbasis ( 1 : 3 , 1 : 3 , 2 ) + &
eigenvalue ( 3 ) * evbasis ( 1 : 3 , 1 : 3 , 3 )
2011-12-01 17:31:13 +05:30
enddo ; enddo ; enddo
2012-08-27 13:34:47 +05:30
end function math_logstrainSpat
2011-12-01 17:31:13 +05:30
2013-01-31 21:58:08 +05:30
2012-08-27 13:34:47 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculate logarithmic strain in material configuration for given F field
!--------------------------------------------------------------------------------------------------
function math_logstrainMat ( F )
2013-03-07 01:04:30 +05:30
2011-12-01 17:31:13 +05:30
implicit none
2012-08-27 13:34:47 +05:30
real ( pReal ) , intent ( in ) , dimension ( : , : , : , : , : ) :: F
2013-01-31 21:58:08 +05:30
real ( pReal ) , dimension ( 3 , 3 , size ( F , 3 ) , size ( F , 4 ) , size ( F , 5 ) ) :: math_logstrainMat
integer ( pInt ) , dimension ( 3 ) :: res
real ( pReal ) , dimension ( 3 , 3 ) :: temp33_Real , temp33_Real2
real ( pReal ) , dimension ( 3 , 3 , 3 ) :: evbasis
real ( pReal ) , dimension ( 3 ) :: eigenvalue
2011-12-01 17:31:13 +05:30
integer ( pInt ) :: i , j , k
logical :: errmatinv
2013-01-31 21:58:08 +05:30
res = [ size ( F , 3 ) , size ( F , 4 ) , size ( F , 5 ) ]
2011-12-01 17:31:13 +05:30
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res ( 1 )
2013-01-31 21:58:08 +05:30
call math_pDecomposition ( F ( 1 : 3 , 1 : 3 , i , j , k ) , temp33_Real , temp33_Real2 , errmatinv ) !store U in temp33_Real
call math_spectralDecomposition ( temp33_Real , eigenvalue ( 1 ) , eigenvalue ( 2 ) , eigenvalue ( 3 ) , &
evbasis ( 1 : 3 , 1 : 3 , 1 ) , evbasis ( 1 : 3 , 1 : 3 , 2 ) , evbasis ( 1 : 3 , 1 : 3 , 3 ) )
2011-12-01 17:31:13 +05:30
eigenvalue = log ( sqrt ( eigenvalue ) )
2013-01-31 21:58:08 +05:30
math_logstrainMat ( 1 : 3 , 1 : 3 , i , j , k ) = eigenvalue ( 1 ) * evbasis ( 1 : 3 , 1 : 3 , 1 ) + &
eigenvalue ( 2 ) * evbasis ( 1 : 3 , 1 : 3 , 2 ) + &
eigenvalue ( 3 ) * evbasis ( 1 : 3 , 1 : 3 , 3 )
2011-12-01 17:31:13 +05:30
enddo ; enddo ; enddo
2012-08-27 13:34:47 +05:30
end function math_logstrainMat
2011-12-01 17:31:13 +05:30
2013-01-31 21:58:08 +05:30
2012-08-27 13:34:47 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculate cauchy stress for given PK1 stress and F field
!--------------------------------------------------------------------------------------------------
function math_cauchy ( F , P )
2013-03-07 01:04:30 +05:30
2011-12-01 17:31:13 +05:30
implicit none
2012-08-27 13:34:47 +05:30
real ( pReal ) , intent ( in ) , dimension ( : , : , : , : , : ) :: F
real ( pReal ) , intent ( in ) , dimension ( : , : , : , : , : ) :: P
2013-01-31 21:58:08 +05:30
real ( pReal ) , dimension ( 3 , 3 , size ( F , 3 ) , size ( F , 4 ) , size ( F , 5 ) ) :: math_cauchy
integer ( pInt ) , dimension ( 3 ) :: res
2011-12-01 17:31:13 +05:30
real ( pReal ) :: jacobi
integer ( pInt ) :: i , j , k
2013-01-31 21:58:08 +05:30
res = [ size ( F , 3 ) , size ( F , 4 ) , size ( F , 5 ) ]
2011-12-01 17:31:13 +05:30
do k = 1_pInt , res ( 3 ) ; do j = 1_pInt , res ( 2 ) ; do i = 1_pInt , res ( 1 )
2013-01-31 21:58:08 +05:30
jacobi = math_det33 ( F ( 1 : 3 , 1 : 3 , i , j , k ) )
math_cauchy ( 1 : 3 , 1 : 3 , i , j , k ) = matmul ( P ( 1 : 3 , 1 : 3 , i , j , k ) , transpose ( F ( 1 : 3 , 1 : 3 , i , j , k ) ) ) / jacobi
2011-12-01 17:31:13 +05:30
enddo ; enddo ; enddo
2012-08-27 13:34:47 +05:30
end function math_cauchy
2012-01-13 21:48:16 +05:30
2012-03-09 01:55:28 +05:30
end module math