2012-07-19 00:09:59 +05:30
!--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @author Christoph Kords, Max-Planck-Institut für Eisenforschung GmbH
2016-01-12 16:30:23 +05:30
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
2017-09-20 03:09:19 +05:30
!> @brief Mathematical library, including random number generation and tensor representations
2012-07-19 00:09:59 +05:30
!--------------------------------------------------------------------------------------------------
2012-08-25 17:16:36 +05:30
module math
2019-05-11 15:21:57 +05:30
use prec
2019-06-10 00:50:38 +05:30
use IO
2020-08-15 19:32:10 +05:30
use config
2020-06-17 16:54:31 +05:30
use YAML_types
2020-04-10 16:22:27 +05:30
use LAPACK_interface
2020-01-29 18:31:14 +05:30
2019-05-11 15:21:57 +05:30
implicit none
2019-05-17 10:54:01 +05:30
public
#if __INTEL_COMPILER >= 1900
! do not make use associated entities available to other modules
private :: &
2020-11-28 11:23:06 +05:30
IO , &
config
2019-05-17 10:54:01 +05:30
#endif
real ( pReal ) , parameter :: PI = acos ( - 1.0_pReal ) !< ratio of a circle's circumference to its diameter
real ( pReal ) , parameter :: INDEG = 18 0.0_pReal / PI !< conversion from radian into degree
real ( pReal ) , parameter :: INRAD = PI / 18 0.0_pReal !< conversion from degree into radian
complex ( pReal ) , parameter :: TWOPIIMG = cmplx ( 0.0_pReal , 2.0_pReal * PI ) !< Re(0.0), Im(2xPi)
real ( pReal ) , dimension ( 3 , 3 ) , parameter :: &
2020-01-29 18:31:14 +05:30
math_I3 = reshape ( [ &
2019-05-11 15:21:57 +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 &
2020-03-29 22:47:24 +05:30
] , shape ( math_I3 ) ) !< 3x3 Identity
2019-05-11 15:21:57 +05:30
2020-03-29 22:47:24 +05:30
real ( pReal ) , dimension ( * ) , parameter , private :: &
NRMMANDEL = [ 1.0_pReal , 1.0_pReal , 1.0_pReal , sqrt ( 2.0_pReal ) , sqrt ( 2.0_pReal ) , sqrt ( 2.0_pReal ) ] !< forward weighting for Mandel notation
2019-05-11 15:21:57 +05:30
2020-03-29 22:47:24 +05:30
real ( pReal ) , dimension ( * ) , parameter , private :: &
2020-02-04 02:09:00 +05:30
INVNRMMANDEL = 1.0_pReal / NRMMANDEL !< backward weighting for Mandel notation
2019-05-11 15:21:57 +05:30
integer , dimension ( 2 , 6 ) , parameter , private :: &
2020-01-29 18:31:14 +05:30
MAPNYE = reshape ( [ &
2019-05-11 15:21:57 +05:30
1 , 1 , &
2 , 2 , &
3 , 3 , &
1 , 2 , &
2 , 3 , &
1 , 3 &
2020-03-29 22:47:24 +05:30
] , shape ( MAPNYE ) ) !< arrangement in Nye notation.
2019-05-11 15:21:57 +05:30
integer , dimension ( 2 , 6 ) , parameter , private :: &
2020-01-29 18:31:14 +05:30
MAPVOIGT = reshape ( [ &
2019-05-11 15:21:57 +05:30
1 , 1 , &
2 , 2 , &
3 , 3 , &
2 , 3 , &
1 , 3 , &
1 , 2 &
2020-03-29 22:47:24 +05:30
] , shape ( MAPVOIGT ) ) !< arrangement in Voigt notation
2019-05-11 15:21:57 +05:30
integer , dimension ( 2 , 9 ) , parameter , private :: &
2020-01-29 18:31:14 +05:30
MAPPLAIN = reshape ( [ &
2019-05-11 15:21:57 +05:30
1 , 1 , &
1 , 2 , &
1 , 3 , &
2 , 1 , &
2 , 2 , &
2 , 3 , &
3 , 1 , &
3 , 2 , &
3 , 3 &
2020-03-29 22:47:24 +05:30
] , shape ( MAPPLAIN ) ) !< arrangement in Plain notation
2020-03-01 14:17:20 +05:30
2019-01-14 11:57:18 +05:30
!---------------------------------------------------------------------------------------------------
2013-01-31 21:58:08 +05:30
private :: &
2020-05-16 20:35:03 +05:30
selfTest
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
!--------------------------------------------------------------------------------------------------
2019-05-11 15:21:57 +05:30
!> @brief initialization of random seed generator and internal checks
2012-07-19 00:09:59 +05:30
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
subroutine math_init
2012-08-25 17:16:36 +05:30
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 4 ) :: randTest
2020-06-17 20:49:21 +05:30
integer :: &
2020-06-29 20:35:11 +05:30
randSize , &
randomSeed !< fixed seeding for pseudo-random number generator, Default 0: use random seed
2019-05-11 15:21:57 +05:30
integer , dimension ( : ) , allocatable :: randInit
2020-06-17 16:54:31 +05:30
class ( tNode ) , pointer :: &
num_generic
2020-09-22 16:39:12 +05:30
print '(/,a)' , ' <<<+- math init -+>>>' ; flush ( IO_STDOUT )
2019-05-11 15:21:57 +05:30
2020-09-13 14:09:17 +05:30
num_generic = > config_numerics % get ( 'generic' , defaultVal = emptyDict )
2020-06-29 20:35:11 +05:30
randomSeed = num_generic % get_asInt ( 'random_seed' , defaultVal = 0 )
2020-06-17 16:54:31 +05:30
2019-05-11 15:21:57 +05:30
call random_seed ( size = randSize )
allocate ( randInit ( randSize ) )
2020-06-29 20:35:11 +05:30
if ( randomSeed > 0 ) then
randInit = randomSeed
2019-05-11 15:21:57 +05:30
else
call random_seed ( )
call random_seed ( get = randInit )
randInit ( 2 : randSize ) = randInit ( 1 )
endif
2020-03-17 00:09:18 +05:30
call random_seed ( put = randInit )
call random_number ( randTest )
2019-05-11 15:21:57 +05:30
2020-09-17 22:58:41 +05:30
print '(a,i2)' , ' size of random seed: ' , randSize
print '(a,i0)' , ' value of random seed: ' , randInit ( 1 )
print '(a,4(/,26x,f17.14),/)' , ' start of random sequence: ' , randTest
2019-05-11 15:21:57 +05:30
call random_seed ( put = randInit )
2020-02-04 02:09:00 +05:30
2020-05-16 20:35:03 +05:30
call selfTest
2017-09-15 00:55:22 +05:30
end subroutine math_init
2019-05-17 10:54:01 +05:30
2009-01-26 18:28:58 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-07-16 18:56:00 +05:30
!> @brief Sorting of two-dimensional integer arrays
!> @details Based on quicksort.
! Sorting is done with respect to array(sortDim,:) and keeps array(/=sortDim,:) linked to it.
! Default: sortDim=1
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-07-16 18:56:00 +05:30
pure recursive subroutine math_sort ( a , istart , iend , sortDim )
2007-04-03 13:47:58 +05:30
2019-05-11 15:21:57 +05:30
integer , dimension ( : , : ) , intent ( inout ) :: a
integer , intent ( in ) , optional :: istart , iend , sortDim
integer :: ipivot , s , e , d
2020-01-29 18:31:14 +05:30
2019-05-11 15:21:57 +05:30
if ( present ( istart ) ) then
s = istart
else
s = lbound ( a , 2 )
endif
2020-01-29 18:31:14 +05:30
2019-05-11 15:21:57 +05:30
if ( present ( iend ) ) then
e = iend
else
e = ubound ( a , 2 )
endif
2020-03-01 14:17:20 +05:30
2019-06-16 00:12:16 +05:30
if ( present ( sortDim ) ) then
2019-05-11 15:21:57 +05:30
d = sortDim
else
d = 1
endif
2020-01-29 18:31:14 +05:30
2019-05-11 15:21:57 +05:30
if ( s < e ) then
2020-07-16 18:56:00 +05:30
call qsort_partition ( a , ipivot , s , e , d )
2019-05-11 15:21:57 +05:30
call math_sort ( a , s , ipivot - 1 , d )
call math_sort ( a , ipivot + 1 , e , d )
endif
2012-08-25 17:16:36 +05:30
2007-04-03 13:47:58 +05:30
2019-05-11 15:21:57 +05:30
contains
!-------------------------------------------------------------------------------------------------
!> @brief Partitioning required for quicksort
!-------------------------------------------------------------------------------------------------
2020-07-16 18:56:00 +05:30
pure subroutine qsort_partition ( a , p , istart , iend , sort )
2020-03-01 14:17:20 +05:30
2019-05-11 15:21:57 +05:30
integer , dimension ( : , : ) , intent ( inout ) :: a
2020-07-16 18:56:00 +05:30
integer , intent ( out ) :: p ! Pivot element
2019-05-11 15:21:57 +05:30
integer , intent ( in ) :: istart , iend , sort
integer , dimension ( size ( a , 1 ) ) :: tmp
integer :: i , j
2020-03-01 14:17:20 +05:30
2019-05-11 15:21:57 +05:30
do
! find the first element on the right side less than or equal to the pivot point
do j = iend , istart , - 1
if ( a ( sort , j ) < = a ( sort , istart ) ) exit
enddo
! find the first element on the left side greater than the pivot point
do i = istart , iend
if ( a ( sort , i ) > a ( sort , istart ) ) exit
enddo
2020-02-04 04:18:09 +05:30
cross : if ( i > = j ) then ! exchange left value with pivot and return with the partition index
2019-05-11 15:21:57 +05:30
tmp = a ( : , istart )
a ( : , istart ) = a ( : , j )
a ( : , j ) = tmp
2020-07-16 18:56:00 +05:30
p = j
2019-05-11 15:21:57 +05:30
return
2020-02-04 04:18:09 +05:30
else cross ! exchange values
2019-05-11 15:21:57 +05:30
tmp = a ( : , i )
a ( : , i ) = a ( : , j )
a ( : , j ) = tmp
endif cross
enddo
2020-03-01 14:17:20 +05:30
2020-07-16 18:56:00 +05:30
end subroutine qsort_partition
2007-04-03 13:47:58 +05:30
2019-05-10 18:33:54 +05:30
end subroutine math_sort
2007-04-03 13:47:58 +05:30
2012-08-25 17:16:36 +05:30
2017-09-15 00:55:22 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief vector expansion
2017-09-20 03:09:19 +05:30
!> @details takes a set of numbers (a,b,c,...) and corresponding multiples (x,y,z,...)
2020-03-01 14:17:20 +05:30
!> to return a vector of x times a, y times b, z times c, ...
2017-09-15 00:55:22 +05:30
!--------------------------------------------------------------------------------------------------
pure function math_expand ( what , how )
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( : ) , intent ( in ) :: what
2020-02-04 04:18:09 +05:30
integer , dimension ( : ) , intent ( in ) :: how
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( sum ( how ) ) :: math_expand
integer :: i
2020-03-01 14:17:20 +05:30
2019-05-11 15:21:57 +05:30
if ( sum ( how ) == 0 ) return
2020-03-01 14:17:20 +05:30
2019-05-11 15:21:57 +05:30
do i = 1 , size ( how )
math_expand ( sum ( how ( 1 : i - 1 ) ) + 1 : sum ( how ( 1 : i ) ) ) = what ( mod ( i - 1 , size ( what ) ) + 1 )
enddo
2017-09-15 00:55:22 +05:30
end function math_expand
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
2019-05-11 15:21:57 +05:30
integer , intent ( in ) :: N !< length of range
integer :: i
integer , dimension ( N ) :: math_range
2009-03-04 17:18:54 +05:30
2019-05-11 15:21:57 +05:30
math_range = [ ( i , i = 1 , 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
!--------------------------------------------------------------------------------------------------
2020-09-13 14:48:09 +05:30
pure function math_eye ( d )
2007-03-29 21:02:52 +05:30
2020-01-29 18:31:14 +05:30
integer , intent ( in ) :: d !< tensor dimension
2019-05-11 15:21:57 +05:30
integer :: i
2020-09-13 14:48:09 +05:30
real ( pReal ) , dimension ( d , d ) :: math_eye
2007-03-29 21:02:52 +05:30
2020-09-13 14:48:09 +05:30
math_eye = 0.0_pReal
2020-01-29 18:31:14 +05:30
do i = 1 , d
2020-09-13 14:48:09 +05:30
math_eye ( i , i ) = 1.0_pReal
2019-05-11 15:21:57 +05:30
enddo
2007-03-29 21:02:52 +05:30
2020-09-13 14:48:09 +05:30
end function math_eye
2007-03-29 21:02:52 +05:30
2019-05-17 10:19:25 +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
!--------------------------------------------------------------------------------------------------
2020-01-29 18:31:14 +05:30
pure function math_identity4th ( d )
2013-01-31 21:58:08 +05:30
2020-01-29 18:31:14 +05:30
integer , intent ( in ) :: d !< tensor dimension
2019-05-11 15:21:57 +05:30
integer :: i , j , k , l
2020-01-29 18:31:14 +05:30
real ( pReal ) , dimension ( d , d , d , d ) :: math_identity4th
real ( pReal ) , dimension ( d , d ) :: identity2nd
2013-01-31 21:58:08 +05:30
2020-09-13 14:48:09 +05:30
identity2nd = math_eye ( d )
2020-01-29 18:31:14 +05:30
do i = 1 , d ; do j = 1 , d ; do k = 1 , d ; do l = 1 , d
2020-02-04 04:18:09 +05:30
math_identity4th ( i , j , k , l ) = 0.5_pReal &
* ( identity2nd ( i , k ) * identity2nd ( j , l ) + identity2nd ( i , l ) * identity2nd ( j , k ) )
2020-01-29 18:31:14 +05:30
enddo ; enddo ; enddo ; enddo
2013-01-31 21:58:08 +05:30
end function math_identity4th
2011-12-01 17:31:13 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-02 03:30:06 +05:30
!> @brief permutation tensor e_ijk
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
!--------------------------------------------------------------------------------------------------
2020-03-02 03:30:06 +05:30
real ( pReal ) pure function math_LeviCivita ( i , j , k )
2008-03-26 19:05:01 +05:30
2019-05-11 15:21:57 +05:30
integer , intent ( in ) :: i , j , k
2008-03-26 19:05:01 +05:30
2020-12-29 02:04:14 +05:30
integer :: o
if ( any ( [ ( all ( cshift ( [ i , j , k ] , o ) == [ 1 , 2 , 3 ] ) , o = 0 , 2 ) ] ) ) then
2020-03-02 03:30:06 +05:30
math_LeviCivita = + 1.0_pReal
2020-12-29 02:04:14 +05:30
elseif ( any ( [ ( all ( cshift ( [ i , j , k ] , o ) == [ 3 , 2 , 1 ] ) , o = 0 , 2 ) ] ) ) then
2020-03-02 03:30:06 +05:30
math_LeviCivita = - 1.0_pReal
else
math_LeviCivita = 0.0_pReal
endif
2008-03-27 17:24:34 +05:30
2020-03-02 03:30:06 +05:30
end function math_LeviCivita
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
2019-05-11 15:21:57 +05:30
integer , intent ( in ) :: i , j
2008-03-27 17:24:34 +05:30
2019-05-11 15:21:57 +05:30
math_delta = merge ( 0.0_pReal , 1.0_pReal , i / = j )
2008-03-27 17:24:34 +05:30
2012-03-09 01:55:28 +05:30
end function math_delta
2007-03-29 21:02:52 +05:30
2011-12-01 17:31:13 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2016-01-10 19:04:26 +05:30
!> @brief cross product a x b
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-03-08 02:52:49 +05:30
pure function math_cross ( A , B )
2009-01-20 00:40:58 +05:30
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: A , B
real ( pReal ) , dimension ( 3 ) :: math_cross
2009-01-20 00:40:58 +05:30
2019-05-11 15:21:57 +05:30
math_cross = [ A ( 2 ) * B ( 3 ) - A ( 3 ) * B ( 2 ) , &
A ( 3 ) * B ( 1 ) - A ( 1 ) * B ( 3 ) , &
A ( 1 ) * B ( 2 ) - A ( 2 ) * B ( 1 ) ]
2009-01-20 00:40:58 +05:30
2019-03-08 02:52:49 +05:30
end function math_cross
2009-01-20 00:40:58 +05:30
2009-03-05 20:07:59 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-02 20:00:39 +05:30
!> @brief outer product of arbitrary sized vectors (A ⊗ B / i,j)
2016-02-26 20:06:24 +05:30
!--------------------------------------------------------------------------------------------------
2019-03-08 02:52:49 +05:30
pure function math_outer ( A , B )
2016-02-26 20:06:24 +05:30
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( : ) , intent ( in ) :: A , B
real ( pReal ) , dimension ( size ( A , 1 ) , size ( B , 1 ) ) :: math_outer
integer :: i , j
2016-02-26 20:06:24 +05:30
2020-01-29 18:31:14 +05:30
do i = 1 , size ( A , 1 ) ; do j = 1 , size ( B , 1 )
math_outer ( i , j ) = A ( i ) * B ( j )
enddo ; enddo
2009-03-05 20:07:59 +05:30
2019-03-08 02:52:49 +05:30
end function math_outer
2009-03-05 20:07:59 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-02 20:00:39 +05:30
!> @brief inner product of arbitrary sized vectors (A · B / i,i)
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-03-08 02:52:49 +05:30
real ( pReal ) pure function math_inner ( A , B )
2009-03-05 20:07:59 +05:30
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( : ) , intent ( in ) :: A
real ( pReal ) , dimension ( size ( A , 1 ) ) , intent ( in ) :: B
2009-03-05 20:07:59 +05:30
2019-05-11 15:21:57 +05:30
math_inner = sum ( A * B )
2009-03-05 20:07:59 +05:30
2019-03-08 02:52:49 +05:30
end function math_inner
2009-01-20 00:40:58 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-02 20:00:39 +05:30
!> @brief double contraction of 3x3 matrices (A : B / ij,ij)
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-01 14:17:20 +05:30
real ( pReal ) pure function math_tensordot ( A , B )
2010-09-30 15:02:49 +05:30
2020-01-29 18:31:14 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: A , B
2010-09-30 15:02:49 +05:30
2020-03-01 14:17:20 +05:30
math_tensordot = sum ( A * B )
2010-09-30 15:02:49 +05:30
2020-03-01 14:17:20 +05:30
end function math_tensordot
2010-09-30 15:02:49 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-02 20:00:39 +05:30
!> @brief matrix double contraction 3333x33 = 33 (ijkl,kl)
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
pure function math_mul3333xx33 ( A , B )
2010-10-13 21:34:44 +05:30
2020-01-29 18:31:14 +05:30
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) , intent ( in ) :: A
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: B
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: math_mul3333xx33
2020-03-01 14:17:20 +05:30
integer :: i , j
2010-10-13 21:34:44 +05:30
2020-01-29 18:31:14 +05:30
do i = 1 , 3 ; do j = 1 , 3
math_mul3333xx33 ( i , j ) = sum ( A ( i , j , 1 : 3 , 1 : 3 ) * B ( 1 : 3 , 1 : 3 ) )
enddo ; enddo
2019-01-14 11:57:18 +05:30
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
!--------------------------------------------------------------------------------------------------
2020-03-02 20:00:39 +05:30
!> @brief matrix multiplication 3333x3333 = 3333 (ijkl,klmn)
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
pure function math_mul3333xx3333 ( A , B )
2012-02-23 01:41:09 +05:30
2019-05-11 15:21:57 +05:30
integer :: i , j , k , l
2020-01-29 18:31:14 +05:30
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) , intent ( in ) :: A
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) , intent ( in ) :: B
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: math_mul3333xx3333
2012-02-23 01:41:09 +05:30
2020-01-29 18:31:14 +05:30
do i = 1 , 3 ; do j = 1 , 3 ; do k = 1 , 3 ; do l = 1 , 3
2019-05-11 15:21:57 +05:30
math_mul3333xx3333 ( i , j , k , l ) = sum ( A ( i , j , 1 : 3 , 1 : 3 ) * B ( 1 : 3 , 1 : 3 , k , l ) )
2020-01-29 18:31:14 +05:30
enddo ; enddo ; enddo ; enddo
2012-02-23 01:41:09 +05:30
2012-03-09 01:55:28 +05:30
end function math_mul3333xx3333
2012-02-23 01:41:09 +05:30
2012-08-25 17:16:36 +05:30
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 )
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: A
2020-03-01 14:11:42 +05:30
integer , intent ( in ) , optional :: n
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: B , math_exp33
2020-03-01 14:17:20 +05:30
2019-05-11 15:21:57 +05:30
real ( pReal ) :: invFac
2020-01-29 18:31:14 +05:30
integer :: n_ , i
2012-10-12 23:24:20 +05:30
2019-09-20 02:50:02 +05:30
if ( present ( n ) ) then
2020-01-10 06:41:19 +05:30
n_ = n
2019-09-20 02:50:02 +05:30
else
2020-01-10 06:41:19 +05:30
n_ = 5
2019-09-20 02:50:02 +05:30
endif
2020-01-10 06:41:19 +05:30
invFac = 1.0_pReal ! 0!
B = math_I3
math_exp33 = math_I3 ! A^0 = I
do i = 1 , n_
2020-01-03 01:19:02 +05:30
invFac = invFac / real ( i , pReal ) ! invfac = 1/(i!)
2019-05-11 15:21:57 +05:30
B = matmul ( B , A )
2020-01-03 01:19:02 +05:30
math_exp33 = math_exp33 + invFac * B ! exp = SUM (A^i)/(i!)
2019-05-11 15:21:57 +05:30
enddo
2012-10-12 23:24:20 +05:30
end function math_exp33
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-15 18:30:35 +05:30
!> @brief Cramer inversion of 3x3 matrix (function)
2020-03-01 14:17:20 +05:30
!> @details Direct Cramer inversion of matrix A. Returns all zeroes if not possible, i.e.
2019-01-14 11:57:18 +05:30
! if determinant is close to zero
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
pure function math_inv33 ( A )
2009-03-31 13:01:38 +05:30
2020-03-01 14:11:42 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: A
real ( pReal ) , dimension ( 3 , 3 ) :: math_inv33
2020-03-01 14:17:20 +05:30
2020-03-01 14:11:42 +05:30
real ( pReal ) :: DetA
logical :: error
2012-08-25 17:16:36 +05:30
2019-09-21 07:03:12 +05:30
call math_invert33 ( math_inv33 , DetA , error , A )
if ( error ) math_inv33 = 0.0_pReal
2009-03-31 13:01:38 +05:30
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
!--------------------------------------------------------------------------------------------------
2020-03-15 18:30:35 +05:30
!> @brief Cramer inversion of 3x3 matrix (subroutine)
2019-01-14 11:57:18 +05:30
!> @details Direct Cramer inversion of matrix A. Also returns determinant
! Returns an error if not possible, i.e. if determinant is close to zero
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-09-21 06:58:46 +05:30
pure subroutine math_invert33 ( InvA , DetA , error , A )
2008-02-15 18:12:27 +05:30
2020-03-01 14:11:42 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( out ) :: InvA
real ( pReal ) , intent ( out ) :: DetA
logical , intent ( out ) :: error
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: A
2008-02-15 18:12:27 +05:30
2019-05-11 15:21:57 +05:30
InvA ( 1 , 1 ) = A ( 2 , 2 ) * A ( 3 , 3 ) - A ( 2 , 3 ) * A ( 3 , 2 )
InvA ( 2 , 1 ) = - A ( 2 , 1 ) * A ( 3 , 3 ) + A ( 2 , 3 ) * A ( 3 , 1 )
InvA ( 3 , 1 ) = A ( 2 , 1 ) * A ( 3 , 2 ) - A ( 2 , 2 ) * A ( 3 , 1 )
2015-05-05 12:07:59 +05:30
2019-05-11 15:21:57 +05:30
DetA = A ( 1 , 1 ) * InvA ( 1 , 1 ) + A ( 1 , 2 ) * InvA ( 2 , 1 ) + A ( 1 , 3 ) * InvA ( 3 , 1 )
2008-02-15 18:12:27 +05:30
2019-05-11 15:21:57 +05:30
if ( dEq0 ( DetA ) ) then
InvA = 0.0_pReal
error = . true .
else
InvA ( 1 , 2 ) = - A ( 1 , 2 ) * A ( 3 , 3 ) + A ( 1 , 3 ) * A ( 3 , 2 )
InvA ( 2 , 2 ) = A ( 1 , 1 ) * A ( 3 , 3 ) - A ( 1 , 3 ) * A ( 3 , 1 )
InvA ( 3 , 2 ) = - A ( 1 , 1 ) * A ( 3 , 2 ) + A ( 1 , 2 ) * A ( 3 , 1 )
2008-02-15 18:12:27 +05:30
2019-05-11 15:21:57 +05:30
InvA ( 1 , 3 ) = A ( 1 , 2 ) * A ( 2 , 3 ) - A ( 1 , 3 ) * A ( 2 , 2 )
InvA ( 2 , 3 ) = - A ( 1 , 1 ) * A ( 2 , 3 ) + A ( 1 , 3 ) * A ( 2 , 1 )
InvA ( 3 , 3 ) = A ( 1 , 1 ) * A ( 2 , 2 ) - A ( 1 , 2 ) * A ( 2 , 1 )
2012-08-25 17:16:36 +05:30
2019-05-11 15:21:57 +05:30
InvA = InvA / DetA
error = . false .
endif
2007-04-11 15:34:22 +05:30
2012-03-09 01:55:28 +05:30
end subroutine math_invert33
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-15 18:30:35 +05:30
!> @brief Inversion of symmetriced 3x3x3x3 matrix
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
function math_invSym3333 ( A )
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: math_invSym3333
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) , intent ( in ) :: A
2020-04-11 03:24:38 +05:30
integer , dimension ( 6 ) :: ipiv6
real ( pReal ) , dimension ( 6 , 6 ) :: temp66
real ( pReal ) , dimension ( 6 * 6 ) :: work
integer :: ierr_i , ierr_f
2019-05-11 15:21:57 +05:30
2019-09-22 10:22:58 +05:30
temp66 = math_sym3333to66 ( A )
2020-03-29 23:34:24 +05:30
call dgetrf ( 6 , 6 , temp66 , 6 , ipiv6 , ierr_i )
call dgetri ( 6 , temp66 , 6 , ipiv6 , work , size ( work , 1 ) , ierr_f )
if ( ierr_i / = 0 . or . ierr_f / = 0 ) then
2020-11-11 20:27:05 +05:30
error stop 'matrix inversion error'
2019-09-22 10:22:58 +05:30
else
math_invSym3333 = math_66toSym3333 ( temp66 )
2019-05-11 15:21:57 +05:30
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
2019-01-13 13:17:01 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-14 11:57:18 +05:30
!> @brief invert quadratic matrix of arbitrary dimension
2019-01-13 13:17:01 +05:30
!--------------------------------------------------------------------------------------------------
2019-09-21 06:50:33 +05:30
subroutine math_invert ( InvA , error , A )
2019-01-13 13:17:01 +05:30
2019-09-21 07:14:23 +05:30
real ( pReal ) , dimension ( : , : ) , intent ( in ) :: A
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( size ( A , 1 ) , size ( A , 1 ) ) , intent ( out ) :: invA
2019-09-21 07:14:23 +05:30
logical , intent ( out ) :: error
2008-02-15 18:12:27 +05:30
2020-04-11 03:24:38 +05:30
integer , dimension ( size ( A , 1 ) ) :: ipiv
real ( pReal ) , dimension ( size ( A , 1 ) ** 2 ) :: work
integer :: ierr
2019-09-21 07:14:23 +05:30
2020-03-01 14:17:20 +05:30
invA = A
2019-09-21 07:14:23 +05:30
call dgetrf ( size ( A , 1 ) , size ( A , 1 ) , invA , size ( A , 1 ) , ipiv , ierr )
2019-09-22 10:22:58 +05:30
error = ( ierr / = 0 )
call dgetri ( size ( A , 1 ) , InvA , size ( A , 1 ) , ipiv , work , size ( work , 1 ) , ierr )
error = error . or . ( ierr / = 0 )
2020-03-01 14:17:20 +05:30
2019-09-21 07:14:23 +05:30
end subroutine math_invert
2008-02-15 18:12:27 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-15 18:30:35 +05:30
!> @brief symmetrize a 3x3 matrix
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2016-01-10 19:04:26 +05:30
pure function math_symmetric33 ( m )
2008-02-15 18:12:27 +05:30
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: math_symmetric33
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m
2019-09-22 09:55:55 +05:30
2019-05-11 15:21:57 +05:30
math_symmetric33 = 0.5_pReal * ( m + transpose ( m ) )
2008-02-15 18:12:27 +05:30
2012-03-09 01:55:28 +05:30
end function math_symmetric33
2008-02-15 18:12:27 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-15 18:30:35 +05:30
!> @brief symmetrize a 6x6 matrix
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
pure function math_symmetric66 ( m )
2008-02-15 18:12:27 +05:30
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 6 , 6 ) :: math_symmetric66
real ( pReal ) , dimension ( 6 , 6 ) , intent ( in ) :: m
2012-08-25 17:16:36 +05:30
2019-05-11 15:21:57 +05:30
math_symmetric66 = 0.5_pReal * ( m + transpose ( m ) )
2012-01-26 19:20:00 +05:30
2012-03-09 01:55:28 +05:30
end function math_symmetric66
2008-02-15 18:12:27 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-15 18:30:35 +05:30
!> @brief skew part of a 3x3 matrix
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2012-02-09 21:28:15 +05:30
pure function math_skew33 ( m )
2012-01-25 16:00:39 +05:30
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: math_skew33
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m
2012-08-25 17:16:36 +05:30
2019-05-11 15:21:57 +05:30
math_skew33 = m - math_symmetric33 ( m )
2012-01-26 19:20:00 +05:30
2012-03-09 01:55:28 +05:30
end function math_skew33
2008-02-15 18:12:27 +05:30
2019-05-17 10:19:25 +05:30
2016-01-09 00:27:37 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-15 18:30:35 +05:30
!> @brief hydrostatic part of a 3x3 matrix
2016-01-09 00:27:37 +05:30
!--------------------------------------------------------------------------------------------------
pure function math_spherical33 ( m )
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: math_spherical33
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m
2016-01-09 00:27:37 +05:30
2019-05-11 15:21:57 +05:30
math_spherical33 = math_I3 * math_trace33 ( m ) / 3.0_pReal
2016-01-09 00:27:37 +05:30
end function math_spherical33
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-15 18:30:35 +05:30
!> @brief deviatoric part of a 3x3 matrix
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2012-02-09 21:28:15 +05:30
pure function math_deviatoric33 ( m )
2012-01-26 19:20:00 +05:30
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: math_deviatoric33
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m
2012-01-26 19:20:00 +05:30
2019-05-11 15:21:57 +05:30
math_deviatoric33 = m - math_spherical33 ( m )
2012-01-26 19:20:00 +05:30
2012-03-09 01:55:28 +05:30
end function math_deviatoric33
2012-01-26 19:20:00 +05:30
2012-10-12 23:24:20 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-15 18:30:35 +05:30
!> @brief trace of a 3x3 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
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m
2012-10-12 23:24:20 +05:30
2019-05-11 15:21:57 +05:30
math_trace33 = m ( 1 , 1 ) + m ( 2 , 2 ) + m ( 3 , 3 )
2012-10-12 23:24:20 +05:30
end function math_trace33
2015-05-05 12:07:59 +05:30
2012-08-27 13:34:47 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-15 18:30:35 +05:30
!> @brief determinant of a 3x3 matrix
2012-08-27 13:34:47 +05:30
!--------------------------------------------------------------------------------------------------
2013-01-31 21:58:08 +05:30
real ( pReal ) pure function math_det33 ( m )
2007-03-20 19:25:22 +05:30
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m
2020-03-01 14:17:20 +05:30
2019-05-11 15:21:57 +05:30
math_det33 = m ( 1 , 1 ) * ( m ( 2 , 2 ) * m ( 3 , 3 ) - m ( 2 , 3 ) * m ( 3 , 2 ) ) &
- m ( 1 , 2 ) * ( m ( 2 , 1 ) * m ( 3 , 3 ) - m ( 2 , 3 ) * m ( 3 , 1 ) ) &
+ m ( 1 , 3 ) * ( m ( 2 , 1 ) * m ( 3 , 2 ) - m ( 2 , 2 ) * m ( 3 , 1 ) )
2007-03-20 19:25:22 +05:30
2012-03-09 01:55:28 +05:30
end function math_det33
2007-03-21 15:50:25 +05:30
2012-08-25 17:16:36 +05:30
2016-02-02 17:53:45 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-15 18:30:35 +05:30
!> @brief determinant of a symmetric 3x3 matrix
2016-02-02 17:53:45 +05:30
!--------------------------------------------------------------------------------------------------
real ( pReal ) pure function math_detSym33 ( m )
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m
2020-03-01 14:17:20 +05:30
2019-03-13 11:16:59 +05:30
math_detSym33 = - ( m ( 1 , 1 ) * m ( 2 , 3 ) ** 2 + m ( 2 , 2 ) * m ( 1 , 3 ) ** 2 + m ( 3 , 3 ) * m ( 1 , 2 ) ** 2 ) &
+ m ( 1 , 1 ) * m ( 2 , 2 ) * m ( 3 , 3 ) + 2.0_pReal * m ( 1 , 2 ) * m ( 1 , 3 ) * m ( 2 , 3 )
2016-02-02 17:53:45 +05:30
end function math_detSym33
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-15 18:30:35 +05:30
!> @brief convert 3x3 matrix into vector 9
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-14 11:57:18 +05:30
pure function math_33to9 ( m33 )
2008-02-15 18:12:27 +05:30
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 9 ) :: math_33to9
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m33
2020-03-01 14:17:20 +05:30
2019-05-11 15:21:57 +05:30
integer :: i
2020-03-01 14:17:20 +05:30
2019-05-11 15:21:57 +05:30
do i = 1 , 9
2020-01-29 18:31:14 +05:30
math_33to9 ( i ) = m33 ( MAPPLAIN ( 1 , i ) , MAPPLAIN ( 2 , i ) )
2019-05-11 15:21:57 +05:30
enddo
2008-02-15 18:12:27 +05:30
2019-01-14 11:57:18 +05:30
end function math_33to9
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-15 18:30:35 +05:30
!> @brief convert 9 vector into 3x3 matrix
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-14 11:57:18 +05:30
pure function math_9to33 ( v9 )
2008-02-15 18:12:27 +05:30
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: math_9to33
real ( pReal ) , dimension ( 9 ) , intent ( in ) :: v9
2020-03-01 14:17:20 +05:30
2019-05-11 15:21:57 +05:30
integer :: i
2012-08-25 17:16:36 +05:30
2019-05-11 15:21:57 +05:30
do i = 1 , 9
2020-01-29 18:31:14 +05:30
math_9to33 ( MAPPLAIN ( 1 , i ) , MAPPLAIN ( 2 , i ) ) = v9 ( i )
2019-05-11 15:21:57 +05:30
enddo
2008-02-15 18:12:27 +05:30
2019-01-14 11:57:18 +05:30
end function math_9to33
2008-02-15 18:12:27 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-15 18:30:35 +05:30
!> @brief convert symmetric 3x3 matrix into 6 vector
2019-01-14 11:57:18 +05:30
!> @details Weighted conversion (default) rearranges according to Nye and weights shear
! components according to Mandel. Advisable for matrix operations.
! Unweighted conversion only changes order according to Nye
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-14 21:06:08 +05:30
pure function math_sym33to6 ( m33 , weighted )
2007-03-28 12:51:47 +05:30
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 6 ) :: math_sym33to6
2020-03-15 17:39:27 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m33 !< symmetric 3x3 matrix (no internal check)
logical , optional , intent ( in ) :: weighted !< weight according to Mandel (.true. by default)
2020-03-01 14:17:20 +05:30
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 6 ) :: w
integer :: i
2020-03-01 14:17:20 +05:30
2019-05-11 15:21:57 +05:30
if ( present ( weighted ) ) then
2020-01-29 18:31:14 +05:30
w = merge ( NRMMANDEL , 1.0_pReal , weighted )
2019-05-11 15:21:57 +05:30
else
2020-01-29 18:31:14 +05:30
w = NRMMANDEL
2019-05-11 15:21:57 +05:30
endif
2012-08-25 17:16:36 +05:30
2019-05-11 15:21:57 +05:30
do i = 1 , 6
2020-01-29 18:31:14 +05:30
math_sym33to6 ( i ) = w ( i ) * m33 ( MAPNYE ( 1 , i ) , MAPNYE ( 2 , i ) )
2020-03-01 14:17:20 +05:30
enddo
2007-03-28 12:51:47 +05:30
2019-01-14 21:06:08 +05:30
end function math_sym33to6
2007-03-28 12:51:47 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-15 18:30:35 +05:30
!> @brief convert 6 vector into symmetric 3x3 matrix
2019-01-14 11:57:18 +05:30
!> @details Weighted conversion (default) rearranges according to Nye and weights shear
! components according to Mandel. Advisable for matrix operations.
! Unweighted conversion only changes order according to Nye
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-14 21:06:08 +05:30
pure function math_6toSym33 ( v6 , weighted )
2007-03-28 12:51:47 +05:30
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: math_6toSym33
2020-03-15 17:39:27 +05:30
real ( pReal ) , dimension ( 6 ) , intent ( in ) :: v6 !< 6 vector
logical , optional , intent ( in ) :: weighted !< weight according to Mandel (.true. by default)
2019-01-14 11:57:18 +05:30
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 6 ) :: w
integer :: i
2020-03-01 14:17:20 +05:30
2019-05-11 15:21:57 +05:30
if ( present ( weighted ) ) then
2020-01-29 18:31:14 +05:30
w = merge ( INVNRMMANDEL , 1.0_pReal , weighted )
2019-05-11 15:21:57 +05:30
else
2020-01-29 18:31:14 +05:30
w = INVNRMMANDEL
2019-05-11 15:21:57 +05:30
endif
2012-08-25 17:16:36 +05:30
2019-05-11 15:21:57 +05:30
do i = 1 , 6
2020-01-29 18:31:14 +05:30
math_6toSym33 ( MAPNYE ( 1 , i ) , MAPNYE ( 2 , i ) ) = w ( i ) * v6 ( i )
math_6toSym33 ( MAPNYE ( 2 , i ) , MAPNYE ( 1 , i ) ) = w ( i ) * v6 ( i )
2019-05-11 15:21:57 +05:30
enddo
2007-03-28 12:51:47 +05:30
2019-01-14 21:06:08 +05:30
end function math_6toSym33
2007-03-28 12:51:47 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-15 18:30:35 +05:30
!> @brief convert 3x3x3x3 matrix into 9x9 matrix
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-14 11:57:18 +05:30
pure function math_3333to99 ( m3333 )
2008-02-15 18:12:27 +05:30
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 9 , 9 ) :: math_3333to99
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) , intent ( in ) :: m3333
2008-02-15 18:12:27 +05:30
2019-05-11 15:21:57 +05:30
integer :: i , j
2012-08-25 17:16:36 +05:30
2019-05-11 15:21:57 +05:30
do i = 1 , 9 ; do j = 1 , 9
2020-01-29 18:31:14 +05:30
math_3333to99 ( i , j ) = m3333 ( MAPPLAIN ( 1 , i ) , MAPPLAIN ( 2 , i ) , MAPPLAIN ( 1 , j ) , MAPPLAIN ( 2 , j ) )
2019-05-11 15:21:57 +05:30
enddo ; enddo
2008-02-15 18:12:27 +05:30
2019-01-14 11:57:18 +05:30
end function math_3333to99
2008-02-15 18:12:27 +05:30
2011-07-29 21:27:39 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-15 18:30:35 +05:30
!> @brief convert 9x9 matrix into 3x3x3x3 matrix
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-14 11:57:18 +05:30
pure function math_99to3333 ( m99 )
2011-07-29 21:27:39 +05:30
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: math_99to3333
real ( pReal ) , dimension ( 9 , 9 ) , intent ( in ) :: m99
2011-07-29 21:27:39 +05:30
2019-05-11 15:21:57 +05:30
integer :: i , j
2012-08-25 17:16:36 +05:30
2019-05-11 15:21:57 +05:30
do i = 1 , 9 ; do j = 1 , 9
2020-01-29 18:31:14 +05:30
math_99to3333 ( MAPPLAIN ( 1 , i ) , MAPPLAIN ( 2 , i ) , MAPPLAIN ( 1 , j ) , MAPPLAIN ( 2 , j ) ) = m99 ( i , j )
2019-05-11 15:21:57 +05:30
enddo ; enddo
2011-07-29 21:27:39 +05:30
2019-01-14 11:57:18 +05:30
end function math_99to3333
2011-07-29 21:27:39 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-15 18:30:35 +05:30
!> @brief convert symmetric 3x3x3x3 matrix into 6x6 matrix
2019-01-14 11:57:18 +05:30
!> @details Weighted conversion (default) rearranges according to Nye and weights shear
! components according to Mandel. Advisable for matrix operations.
2020-03-15 17:39:27 +05:30
! Unweighted conversion only rearranges order according to Nye
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-14 21:06:08 +05:30
pure function math_sym3333to66 ( m3333 , weighted )
2007-03-28 12:51:47 +05:30
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 6 , 6 ) :: math_sym3333to66
2020-03-15 17:39:27 +05:30
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) , intent ( in ) :: m3333 !< symmetric 3x3x3x3 matrix (no internal check)
logical , optional , intent ( in ) :: weighted !< weight according to Mandel (.true. by default)
2019-01-14 11:57:18 +05:30
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 6 ) :: w
integer :: i , j
2020-03-01 14:17:20 +05:30
2019-05-11 15:21:57 +05:30
if ( present ( weighted ) ) then
2020-01-29 18:31:14 +05:30
w = merge ( NRMMANDEL , 1.0_pReal , weighted )
2019-05-11 15:21:57 +05:30
else
2020-01-29 18:31:14 +05:30
w = NRMMANDEL
2019-05-11 15:21:57 +05:30
endif
2012-08-25 17:16:36 +05:30
2019-05-11 15:21:57 +05:30
do i = 1 , 6 ; do j = 1 , 6
2020-01-29 18:31:14 +05:30
math_sym3333to66 ( i , j ) = w ( i ) * w ( j ) * m3333 ( MAPNYE ( 1 , i ) , MAPNYE ( 2 , i ) , MAPNYE ( 1 , j ) , MAPNYE ( 2 , j ) )
2019-05-11 15:21:57 +05:30
enddo ; enddo
2007-03-28 12:51:47 +05:30
2019-01-14 21:06:08 +05:30
end function math_sym3333to66
2010-05-06 19:37:21 +05:30
2011-12-01 17:31:13 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-15 18:30:35 +05:30
!> @brief convert 66 matrix into symmetric 3x3x3x3 matrix
2019-01-14 11:57:18 +05:30
!> @details Weighted conversion (default) rearranges according to Nye and weights shear
! components according to Mandel. Advisable for matrix operations.
2020-03-15 17:39:27 +05:30
! Unweighted conversion only rearranges order according to Nye
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-14 21:06:08 +05:30
pure function math_66toSym3333 ( m66 , weighted )
2008-02-15 18:12:27 +05:30
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: math_66toSym3333
2020-03-15 17:39:27 +05:30
real ( pReal ) , dimension ( 6 , 6 ) , intent ( in ) :: m66 !< 6x6 matrix
logical , optional , intent ( in ) :: weighted !< weight according to Mandel (.true. by default)
2020-03-01 14:17:20 +05:30
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 6 ) :: w
integer :: i , j
2020-03-01 14:17:20 +05:30
2019-05-11 15:21:57 +05:30
if ( present ( weighted ) ) then
2020-01-29 18:31:14 +05:30
w = merge ( INVNRMMANDEL , 1.0_pReal , weighted )
2019-05-11 15:21:57 +05:30
else
2020-01-29 18:31:14 +05:30
w = INVNRMMANDEL
2019-05-11 15:21:57 +05:30
endif
do i = 1 , 6 ; do j = 1 , 6
2020-01-29 18:31:14 +05:30
math_66toSym3333 ( MAPNYE ( 1 , i ) , MAPNYE ( 2 , i ) , MAPNYE ( 1 , j ) , MAPNYE ( 2 , j ) ) = w ( i ) * w ( j ) * m66 ( i , j )
math_66toSym3333 ( MAPNYE ( 2 , i ) , MAPNYE ( 1 , i ) , MAPNYE ( 1 , j ) , MAPNYE ( 2 , j ) ) = w ( i ) * w ( j ) * m66 ( i , j )
math_66toSym3333 ( MAPNYE ( 1 , i ) , MAPNYE ( 2 , i ) , MAPNYE ( 2 , j ) , MAPNYE ( 1 , j ) ) = w ( i ) * w ( j ) * m66 ( i , j )
math_66toSym3333 ( MAPNYE ( 2 , i ) , MAPNYE ( 1 , i ) , MAPNYE ( 2 , j ) , MAPNYE ( 1 , j ) ) = w ( i ) * w ( j ) * m66 ( i , j )
2019-05-11 15:21:57 +05:30
enddo ; enddo
2008-02-15 18:12:27 +05:30
2019-01-14 21:06:08 +05:30
end function math_66toSym3333
2008-02-15 18:12:27 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-15 18:30:35 +05:30
!> @brief convert 66 Voigt matrix into symmetric 3x3x3x3 matrix
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
pure function math_Voigt66to3333 ( m66 )
2008-02-15 18:12:27 +05:30
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: math_Voigt66to3333
2020-03-15 17:39:27 +05:30
real ( pReal ) , dimension ( 6 , 6 ) , intent ( in ) :: m66 !< 6x6 matrix
2019-05-11 15:21:57 +05:30
integer :: i , j
2012-08-25 17:16:36 +05:30
2019-05-11 15:21:57 +05:30
do i = 1 , 6 ; do j = 1 , 6
2020-01-29 18:31:14 +05:30
math_Voigt66to3333 ( MAPVOIGT ( 1 , i ) , MAPVOIGT ( 2 , i ) , MAPVOIGT ( 1 , j ) , MAPVOIGT ( 2 , j ) ) = m66 ( i , j )
math_Voigt66to3333 ( MAPVOIGT ( 2 , i ) , MAPVOIGT ( 1 , i ) , MAPVOIGT ( 1 , j ) , MAPVOIGT ( 2 , j ) ) = m66 ( i , j )
math_Voigt66to3333 ( MAPVOIGT ( 1 , i ) , MAPVOIGT ( 2 , i ) , MAPVOIGT ( 2 , j ) , MAPVOIGT ( 1 , j ) ) = m66 ( i , j )
math_Voigt66to3333 ( MAPVOIGT ( 2 , i ) , MAPVOIGT ( 1 , i ) , MAPVOIGT ( 2 , j ) , MAPVOIGT ( 1 , j ) ) = m66 ( i , j )
2019-05-11 15:21:57 +05:30
enddo ; enddo
2008-02-15 18:12:27 +05:30
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 draw a random sample from Gauss variable
!--------------------------------------------------------------------------------------------------
real ( pReal ) function math_sampleGaussVar ( meanvalue , stddev , width )
2020-03-15 17:39:27 +05:30
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
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 2 ) :: rnd ! random numbers
real ( pReal ) :: scatter , & ! normalized scatter around meanvalue
2020-01-10 06:41:19 +05:30
width_
2019-05-11 15:21:57 +05:30
if ( abs ( stddev ) < tol_math_check ) then
math_sampleGaussVar = meanvalue
else
2020-01-10 06:41:19 +05:30
if ( present ( width ) ) then
width_ = width
else
width_ = 3.0_pReal ! use +-3*sigma as default scatter
endif
2019-05-11 15:21:57 +05:30
do
call random_number ( rnd )
2020-01-10 06:41:19 +05:30
scatter = width_ * ( 2.0_pReal * rnd ( 1 ) - 1.0_pReal )
2019-05-11 15:21:57 +05:30
if ( rnd ( 2 ) < = exp ( - 0.5_pReal * scatter ** 2.0_pReal ) ) exit ! test if scattered value is drawn
enddo
math_sampleGaussVar = scatter * stddev
endif
2013-01-31 21:58:08 +05:30
end function math_sampleGaussVar
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-15 18:30:35 +05:30
!> @brief eigenvalues and eigenvectors of symmetric matrix
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-07-25 03:10:42 +05:30
subroutine math_eigh ( w , v , error , m )
2011-09-14 18:56:00 +05:30
2020-03-15 17:39:27 +05:30
real ( pReal ) , dimension ( : , : ) , intent ( in ) :: m !< quadratic matrix to compute eigenvectors and values of
real ( pReal ) , dimension ( size ( m , 1 ) ) , intent ( out ) :: w !< eigenvalues
real ( pReal ) , dimension ( size ( m , 1 ) , size ( m , 1 ) ) , intent ( out ) :: v !< eigenvectors
2020-07-25 03:10:42 +05:30
logical , intent ( out ) :: error
2019-09-22 10:22:58 +05:30
integer :: ierr
2020-04-11 03:24:38 +05:30
real ( pReal ) , dimension ( size ( m , 1 ) ** 2 ) :: work
2012-08-25 17:16:36 +05:30
2020-03-09 18:47:05 +05:30
v = m ! copy matrix to input (doubles as output) array
call dsyev ( 'V' , 'U' , size ( m , 1 ) , v , size ( m , 1 ) , w , work , size ( work , 1 ) , ierr )
2019-09-22 10:22:58 +05:30
error = ( ierr / = 0 )
2012-08-25 17:16:36 +05:30
2020-03-09 18:47:05 +05:30
end subroutine math_eigh
2016-01-31 16:55:26 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-15 18:30:35 +05:30
!> @brief eigenvalues and eigenvectors of symmetric 3x3 matrix using an analytical expression
2016-02-27 00:51:47 +05:30
!> and the general LAPACK powered version for arbritrary sized matrices as fallback
2019-01-27 13:09:11 +05:30
!> @author Joachim Kopp, Max-Planck-Institut für Kernphysik, Heidelberg (Copyright (C) 2006)
2016-02-27 00:51:47 +05:30
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @details See http://arxiv.org/abs/physics/0610206 (DSYEVH3)
2016-01-31 16:55:26 +05:30
!--------------------------------------------------------------------------------------------------
2020-07-25 03:10:42 +05:30
subroutine math_eigh33 ( w , v , m )
2020-03-01 14:17:20 +05:30
2020-03-15 17:39:27 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m !< 3x3 matrix to compute eigenvectors and values of
real ( pReal ) , dimension ( 3 ) , intent ( out ) :: w !< eigenvalues
real ( pReal ) , dimension ( 3 , 3 ) , intent ( out ) :: v !< eigenvectors
2019-05-11 15:21:57 +05:30
real ( pReal ) :: T , U , norm , threshold
logical :: error
2020-03-01 14:17:20 +05:30
2020-03-09 18:47:05 +05:30
w = math_eigvalsh33 ( m )
2020-03-01 14:17:20 +05:30
2020-03-09 18:47:05 +05:30
v ( 1 : 3 , 2 ) = [ m ( 1 , 2 ) * m ( 2 , 3 ) - m ( 1 , 3 ) * m ( 2 , 2 ) , &
m ( 1 , 3 ) * m ( 1 , 2 ) - m ( 2 , 3 ) * m ( 1 , 1 ) , &
m ( 1 , 2 ) ** 2 ]
2020-03-01 14:17:20 +05:30
2020-03-09 18:47:05 +05:30
T = maxval ( abs ( w ) )
2019-05-11 15:21:57 +05:30
U = max ( T , T ** 2 )
threshold = sqrt ( 5.68e-14_pReal * U ** 2 )
2020-03-01 14:17:20 +05:30
2020-03-09 18:47:05 +05:30
v ( 1 : 3 , 1 ) = [ v ( 1 , 2 ) + m ( 1 , 3 ) * w ( 1 ) , &
v ( 2 , 2 ) + m ( 2 , 3 ) * w ( 1 ) , &
( m ( 1 , 1 ) - w ( 1 ) ) * ( m ( 2 , 2 ) - w ( 1 ) ) - v ( 3 , 2 ) ]
norm = norm2 ( v ( 1 : 3 , 1 ) )
2019-05-11 15:21:57 +05:30
fallback1 : if ( norm < threshold ) then
2020-07-25 03:10:42 +05:30
call math_eigh ( w , v , error , m )
2020-03-09 18:47:05 +05:30
else fallback1
v ( 1 : 3 , 1 ) = v ( 1 : 3 , 1 ) / norm
v ( 1 : 3 , 2 ) = [ v ( 1 , 2 ) + m ( 1 , 3 ) * w ( 2 ) , &
v ( 2 , 2 ) + m ( 2 , 3 ) * w ( 2 ) , &
( m ( 1 , 1 ) - w ( 2 ) ) * ( m ( 2 , 2 ) - w ( 2 ) ) - v ( 3 , 2 ) ]
norm = norm2 ( v ( 1 : 3 , 2 ) )
fallback2 : if ( norm < threshold ) then
2020-07-25 03:10:42 +05:30
call math_eigh ( w , v , error , m )
2020-03-09 18:47:05 +05:30
else fallback2
v ( 1 : 3 , 2 ) = v ( 1 : 3 , 2 ) / norm
v ( 1 : 3 , 3 ) = math_cross ( v ( 1 : 3 , 1 ) , v ( 1 : 3 , 2 ) )
endif fallback2
2019-05-11 15:21:57 +05:30
endif fallback1
2020-03-01 14:17:20 +05:30
2020-03-09 18:47:05 +05:30
end subroutine math_eigh33
2016-01-31 16:55:26 +05:30
2016-02-26 20:06:24 +05:30
2013-01-31 21:58:08 +05:30
!--------------------------------------------------------------------------------------------------
2020-07-25 02:38:05 +05:30
!> @brief Calculate rotational part of a deformation gradient
!> @details https://www.jstor.org/stable/43637254
!! https://www.jstor.org/stable/43637372
!! https://doi.org/10.1023/A:1007407802076
!--------------------------------------------------------------------------------------------------
pure function math_rotationalPart ( F ) result ( R )
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: &
F ! deformation gradient
real ( pReal ) , dimension ( 3 , 3 ) :: &
C , & ! right Cauchy-Green tensor
R ! rotational part
real ( pReal ) , dimension ( 3 ) :: &
lambda , & ! principal stretches
I_C , & ! invariants of C
I_U ! invariants of U
real ( pReal ) , dimension ( 2 ) :: &
I_F ! first two invariants of F
real ( pReal ) :: x , Phi
2013-01-31 21:58:08 +05:30
2020-07-25 02:38:05 +05:30
C = matmul ( transpose ( F ) , F )
I_C = math_invariantsSym33 ( C )
I_F = [ math_trace33 ( F ) , 0.5 * ( math_trace33 ( F ) ** 2 - math_trace33 ( matmul ( F , F ) ) ) ]
2016-01-10 19:04:26 +05:30
2020-07-25 02:38:05 +05:30
x = math_clip ( I_C ( 1 ) ** 2 - 3.0_pReal * I_C ( 2 ) , 0.0_pReal ) ** ( 3.0_pReal / 2.0_pReal )
if ( dNeq0 ( x ) ) then
Phi = acos ( math_clip ( ( I_C ( 1 ) ** 3 - 4.5_pReal * I_C ( 1 ) * I_C ( 2 ) + 1 3.5_pReal * I_C ( 3 ) ) / x , - 1.0_pReal , 1.0_pReal ) )
lambda = I_C ( 1 ) + ( 2.0_pReal * sqrt ( math_clip ( I_C ( 1 ) ** 2 - 3.0_pReal * I_C ( 2 ) , 0.0_pReal ) ) ) &
* cos ( ( Phi - 2.0_pReal * PI * [ 1.0_pReal , 2.0_pReal , 3.0_pReal ] ) / 3.0_pReal )
lambda = sqrt ( math_clip ( lambda , 0.0_pReal ) / 3.0_pReal )
else
lambda = sqrt ( I_C ( 1 ) / 3.0_pReal )
endif
2016-02-26 20:06:24 +05:30
2020-07-25 02:38:05 +05:30
I_U = [ sum ( lambda ) , lambda ( 1 ) * lambda ( 2 ) + lambda ( 2 ) * lambda ( 3 ) + lambda ( 3 ) * lambda ( 1 ) , product ( lambda ) ]
2013-01-31 21:58:08 +05:30
2020-07-25 02:38:05 +05:30
R = I_U ( 1 ) * I_F ( 2 ) * math_I3 &
+ ( I_U ( 1 ) ** 2 - I_U ( 2 ) ) * F &
- I_U ( 1 ) * I_F ( 1 ) * transpose ( F ) &
+ I_U ( 1 ) * transpose ( matmul ( F , F ) ) &
- matmul ( F , C )
R = R / ( I_U ( 1 ) * I_U ( 2 ) - I_U ( 3 ) )
2020-03-15 20:41:28 +05:30
2020-03-15 18:51:11 +05:30
end function math_rotationalPart
2007-03-20 19:25:22 +05:30
2011-12-01 17:31:13 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-15 18:30:35 +05:30
!> @brief Eigenvalues of symmetric matrix
2017-05-04 04:02:44 +05:30
! will return NaN on error
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-09 18:47:05 +05:30
function math_eigvalsh ( m )
2020-03-01 14:17:20 +05:30
2020-03-15 17:39:27 +05:30
real ( pReal ) , dimension ( : , : ) , intent ( in ) :: m !< symmetric matrix to compute eigenvalues of
2020-03-09 18:47:05 +05:30
real ( pReal ) , dimension ( size ( m , 1 ) ) :: math_eigvalsh
2020-03-15 17:39:27 +05:30
2019-09-22 10:22:58 +05:30
real ( pReal ) , dimension ( size ( m , 1 ) , size ( m , 1 ) ) :: m_
integer :: ierr
2020-04-11 03:24:38 +05:30
real ( pReal ) , dimension ( size ( m , 1 ) ** 2 ) :: work
2019-05-11 15:21:57 +05:30
2019-09-22 10:22:58 +05:30
m_ = m ! copy matrix to input (will be destroyed)
2020-03-09 18:47:05 +05:30
call dsyev ( 'N' , 'U' , size ( m , 1 ) , m_ , size ( m , 1 ) , math_eigvalsh , work , size ( work , 1 ) , ierr )
if ( ierr / = 0 ) math_eigvalsh = IEEE_value ( 1.0_pReal , IEEE_quiet_NaN )
2016-01-10 19:04:26 +05:30
2020-03-09 18:47:05 +05:30
end function math_eigvalsh
2007-03-20 19:25:22 +05:30
2016-02-02 13:39:07 +05:30
2016-02-01 12:18:42 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-15 18:30:35 +05:30
!> @brief eigenvalues of symmetric 3x3 matrix using an analytical expression
2016-02-02 23:29:04 +05:30
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @details similar to http://arxiv.org/abs/physics/0610206 (DSYEVC3)
!> but apparently more stable solution and has general LAPACK powered version for arbritrary sized
!> matrices as fallback
2016-02-01 12:18:42 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-09 18:47:05 +05:30
function math_eigvalsh33 ( m )
2016-02-02 02:07:27 +05:30
2020-03-15 17:39:27 +05:30
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: m !< 3x3 symmetric matrix to compute eigenvalues of
2020-03-15 18:30:35 +05:30
real ( pReal ) , dimension ( 3 ) :: math_eigvalsh33 , I
2019-05-11 15:21:57 +05:30
real ( pReal ) :: P , Q , rho , phi
real ( pReal ) , parameter :: TOL = 1.e-14_pReal
2016-02-02 02:07:27 +05:30
2020-03-15 18:30:35 +05:30
I = math_invariantsSym33 ( m ) ! invariants are coefficients in characteristic polynomial apart for the sign of c0 and c2 in http://arxiv.org/abs/physics/0610206
2016-02-02 17:53:45 +05:30
2020-03-15 18:30:35 +05:30
P = I ( 2 ) - I ( 1 ) ** 2.0_pReal / 3.0_pReal ! different from http://arxiv.org/abs/physics/0610206 (this formulation was in DAMASK)
Q = product ( I ( 1 : 2 ) ) / 3.0_pReal &
- 2.0_pReal / 2 7.0_pReal * I ( 1 ) ** 3.0_pReal &
- I ( 3 ) ! different from http://arxiv.org/abs/physics/0610206 (this formulation was in DAMASK)
2016-02-02 02:07:27 +05:30
2019-05-11 15:21:57 +05:30
if ( all ( abs ( [ P , Q ] ) < TOL ) ) then
2020-03-09 18:47:05 +05:30
math_eigvalsh33 = math_eigvalsh ( m )
2019-05-11 15:21:57 +05:30
else
rho = sqrt ( - 3.0_pReal * P ** 3.0_pReal ) / 9.0_pReal
phi = acos ( math_clip ( - Q / rho * 0.5_pReal , - 1.0_pReal , 1.0_pReal ) )
2020-03-09 18:47:05 +05:30
math_eigvalsh33 = 2.0_pReal * rho ** ( 1.0_pReal / 3.0_pReal ) * &
2020-07-25 03:10:42 +05:30
[ cos ( phi / 3.0_pReal ) , &
2020-03-09 18:47:05 +05:30
cos ( ( phi + 2.0_pReal * PI ) / 3.0_pReal ) , &
cos ( ( phi + 4.0_pReal * PI ) / 3.0_pReal ) &
] &
2020-03-15 18:30:35 +05:30
+ I ( 1 ) / 3.0_pReal
2019-05-11 15:21:57 +05:30
endif
2017-09-06 19:50:24 +05:30
2020-03-09 18:47:05 +05:30
end function math_eigvalsh33
2011-12-01 17:31:13 +05:30
2016-02-02 17:53:45 +05:30
2012-08-25 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-15 18:30:35 +05:30
!> @brief invariants of symmetrix 3x3 matrix
2016-02-02 14:14:51 +05:30
!--------------------------------------------------------------------------------------------------
pure function math_invariantsSym33 ( m )
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: m
real ( pReal ) , dimension ( 3 ) :: math_invariantsSym33
2016-02-02 14:14:51 +05:30
2019-05-11 15:21:57 +05:30
math_invariantsSym33 ( 1 ) = math_trace33 ( m )
math_invariantsSym33 ( 2 ) = m ( 1 , 1 ) * m ( 2 , 2 ) + m ( 1 , 1 ) * m ( 3 , 3 ) + m ( 2 , 2 ) * m ( 3 , 3 ) &
- ( m ( 1 , 2 ) ** 2 + m ( 1 , 3 ) ** 2 + m ( 2 , 3 ) ** 2 )
math_invariantsSym33 ( 3 ) = math_detSym33 ( m )
2016-02-02 14:14:51 +05:30
end function math_invariantsSym33
2014-08-21 14:03:55 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief factorial
!--------------------------------------------------------------------------------------------------
2019-03-13 11:16:59 +05:30
integer pure function math_factorial ( n )
2014-08-21 14:03:55 +05:30
2019-05-11 15:21:57 +05:30
integer , intent ( in ) :: n
2020-01-29 18:31:14 +05:30
2019-09-22 19:11:55 +05:30
math_factorial = product ( math_range ( n ) )
2014-08-21 14:03:55 +05:30
end function math_factorial
!--------------------------------------------------------------------------------------------------
!> @brief binomial coefficient
!--------------------------------------------------------------------------------------------------
2019-03-13 11:16:59 +05:30
integer pure function math_binomial ( n , k )
2014-08-21 14:03:55 +05:30
2019-05-11 15:21:57 +05:30
integer , intent ( in ) :: n , k
2019-09-23 10:11:00 +05:30
integer :: i , k_ , n_
2020-03-01 14:17:20 +05:30
2019-09-23 10:11:00 +05:30
k_ = min ( k , n - k )
n_ = n
math_binomial = merge ( 1 , 0 , k_ > - 1 ) ! handling special cases k < 0 or k > n
do i = 1 , k_
math_binomial = ( math_binomial * n_ ) / i
n_ = n_ - 1
enddo
2014-08-21 14:03:55 +05:30
end function math_binomial
!--------------------------------------------------------------------------------------------------
!> @brief multinomial coefficient
!--------------------------------------------------------------------------------------------------
2019-03-13 11:16:59 +05:30
integer pure function math_multinomial ( alpha )
2014-08-21 14:03:55 +05:30
2019-05-11 15:21:57 +05:30
integer , intent ( in ) , dimension ( : ) :: alpha
integer :: i
2020-01-29 18:31:14 +05:30
2019-05-11 15:21:57 +05:30
math_multinomial = 1
do i = 1 , size ( alpha )
math_multinomial = math_multinomial * math_binomial ( sum ( alpha ( 1 : i ) ) , alpha ( i ) )
2020-03-01 14:17:20 +05:30
enddo
2014-08-21 14:03:55 +05:30
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
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: v1 , v2 , v3 , v4
real ( pReal ) , dimension ( 3 , 3 ) :: m
2007-03-20 19:25:22 +05:30
2019-05-11 15:21:57 +05:30
m ( 1 : 3 , 1 ) = v1 - v2
2019-10-17 00:26:31 +05:30
m ( 1 : 3 , 2 ) = v1 - v3
m ( 1 : 3 , 3 ) = v1 - v4
2009-01-20 00:40:58 +05:30
2019-10-17 00:26:31 +05:30
math_volTetrahedron = abs ( 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 )
2019-05-11 15:21:57 +05:30
real ( pReal ) , dimension ( 3 ) , intent ( in ) :: v1 , v2 , v3
2013-04-09 23:37:30 +05:30
2019-05-11 15:21:57 +05:30
math_areaTriangle = 0.5_pReal * norm2 ( math_cross ( v1 - v2 , v1 - v3 ) )
2013-04-09 23:37:30 +05:30
end function math_areaTriangle
2016-01-10 19:04:26 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief limits a scalar value to a certain range (either one or two sided)
! Will return NaN if left > right
!--------------------------------------------------------------------------------------------------
2018-12-05 00:05:29 +05:30
real ( pReal ) pure elemental function math_clip ( a , left , right )
2016-01-10 19:04:26 +05:30
2019-05-11 15:21:57 +05:30
real ( pReal ) , intent ( in ) :: a
real ( pReal ) , intent ( in ) , optional :: left , right
2016-01-10 19:04:26 +05:30
2019-05-11 15:21:57 +05:30
math_clip = a
2019-09-21 05:09:33 +05:30
if ( present ( left ) ) math_clip = max ( left , math_clip )
if ( present ( right ) ) math_clip = min ( right , math_clip )
2019-05-11 15:21:57 +05:30
if ( present ( left ) . and . present ( right ) ) &
math_clip = merge ( IEEE_value ( 1.0_pReal , IEEE_quiet_NaN ) , math_clip , left > right )
2016-01-10 19:04:26 +05:30
2018-07-28 05:01:02 +05:30
end function math_clip
2018-02-14 17:33:50 +05:30
2019-09-22 09:55:55 +05:30
!--------------------------------------------------------------------------------------------------
2020-09-13 15:39:32 +05:30
!> @brief Check correctness of some math functions.
2019-09-22 09:55:55 +05:30
!--------------------------------------------------------------------------------------------------
2020-05-16 20:35:03 +05:30
subroutine selfTest
2020-01-29 18:31:14 +05:30
2019-09-22 09:55:55 +05:30
integer , dimension ( 2 , 4 ) :: &
sort_in_ = reshape ( [ + 1 , + 5 , + 5 , + 6 , - 1 , - 1 , + 3 , - 2 ] , [ 2 , 4 ] )
integer , dimension ( 2 , 4 ) , parameter :: &
sort_out_ = reshape ( [ - 1 , - 1 , + 1 , + 5 , + 5 , + 6 , + 3 , - 2 ] , [ 2 , 4 ] )
2020-01-29 18:31:14 +05:30
2019-09-23 18:59:32 +05:30
integer , dimension ( 5 ) :: range_out_ = [ 1 , 2 , 3 , 4 , 5 ]
2020-03-02 03:30:06 +05:30
integer , dimension ( 3 ) :: ijk
2019-09-22 09:55:55 +05:30
real ( pReal ) :: det
2019-10-17 00:43:26 +05:30
real ( pReal ) , dimension ( 3 ) :: v3_1 , v3_2 , v3_3 , v3_4
2019-09-22 09:55:55 +05:30
real ( pReal ) , dimension ( 6 ) :: v6
real ( pReal ) , dimension ( 9 ) :: v9
real ( pReal ) , dimension ( 3 , 3 ) :: t33 , t33_2
real ( pReal ) , dimension ( 6 , 6 ) :: t66
real ( pReal ) , dimension ( 9 , 9 ) :: t99 , t99_2
2020-01-29 18:31:14 +05:30
real ( pReal ) , dimension ( : , : ) , &
allocatable :: txx , txx_2
real ( pReal ) :: r
integer :: d
2019-09-22 09:55:55 +05:30
logical :: e
2020-01-29 18:31:14 +05:30
2019-09-22 09:55:55 +05:30
if ( any ( abs ( [ 1.0_pReal , 2.0_pReal , 2.0_pReal , 3.0_pReal , 3.0_pReal , 3.0_pReal ] - &
math_expand ( [ 1.0_pReal , 2.0_pReal , 3.0_pReal ] , [ 1 , 2 , 3 , 0 ] ) ) > tol_math_check ) ) &
2020-09-13 15:39:32 +05:30
error stop 'math_expand [1,2,3] by [1,2,3,0] => [1,2,2,3,3,3]'
2019-09-22 09:55:55 +05:30
if ( any ( abs ( [ 1.0_pReal , 2.0_pReal , 2.0_pReal ] - &
math_expand ( [ 1.0_pReal , 2.0_pReal , 3.0_pReal ] , [ 1 , 2 ] ) ) > tol_math_check ) ) &
2020-09-13 15:39:32 +05:30
error stop 'math_expand [1,2,3] by [1,2] => [1,2,2]'
2019-09-22 09:55:55 +05:30
if ( any ( abs ( [ 1.0_pReal , 2.0_pReal , 2.0_pReal , 1.0_pReal , 1.0_pReal , 1.0_pReal ] - &
math_expand ( [ 1.0_pReal , 2.0_pReal ] , [ 1 , 2 , 3 ] ) ) > tol_math_check ) ) &
2020-09-13 15:39:32 +05:30
error stop 'math_expand [1,2] by [1,2,3] => [1,2,2,1,1,1]'
2019-09-22 09:55:55 +05:30
call math_sort ( sort_in_ , 1 , 3 , 2 )
if ( any ( sort_in_ / = sort_out_ ) ) &
2020-09-13 15:39:32 +05:30
error stop 'math_sort'
2020-01-29 18:31:14 +05:30
2019-09-22 09:55:55 +05:30
if ( any ( math_range ( 5 ) / = range_out_ ) ) &
2020-09-13 15:39:32 +05:30
error stop 'math_range'
2020-01-29 18:31:14 +05:30
2020-09-13 15:39:32 +05:30
if ( any ( dNeq ( math_exp33 ( math_I3 , 0 ) , math_I3 ) ) ) &
error stop 'math_exp33(math_I3,1)'
2020-11-12 01:27:17 +05:30
if ( any ( dNeq ( math_exp33 ( math_I3 , 128 ) , exp ( 1.0_pReal ) * math_I3 ) ) ) &
error stop 'math_exp33(math_I3,128)'
2020-01-29 18:31:14 +05:30
2019-09-22 09:55:55 +05:30
call random_number ( v9 )
if ( any ( dNeq ( math_33to9 ( math_9to33 ( v9 ) ) , v9 ) ) ) &
2020-09-13 15:39:32 +05:30
error stop 'math_33to9/math_9to33'
2020-03-01 14:17:20 +05:30
2019-09-22 09:55:55 +05:30
call random_number ( t99 )
if ( any ( dNeq ( math_3333to99 ( math_99to3333 ( t99 ) ) , t99 ) ) ) &
2020-09-13 15:39:32 +05:30
error stop 'math_3333to99/math_99to3333'
2020-01-29 18:31:14 +05:30
2019-09-22 09:55:55 +05:30
call random_number ( v6 )
if ( any ( dNeq ( math_sym33to6 ( math_6toSym33 ( v6 ) ) , v6 ) ) ) &
2020-09-13 15:39:32 +05:30
error stop 'math_sym33to6/math_6toSym33'
2020-01-29 18:31:14 +05:30
2019-09-22 09:55:55 +05:30
call random_number ( t66 )
if ( any ( dNeq ( math_sym3333to66 ( math_66toSym3333 ( t66 ) ) , t66 ) ) ) &
2020-09-13 15:39:32 +05:30
error stop 'math_sym3333to66/math_66toSym3333'
2020-01-29 18:31:14 +05:30
2019-09-22 09:55:55 +05:30
call random_number ( v6 )
if ( any ( dNeq0 ( math_6toSym33 ( v6 ) - math_symmetric33 ( math_6toSym33 ( v6 ) ) ) ) ) &
2020-09-13 15:39:32 +05:30
error stop 'math_symmetric33'
2020-01-29 18:31:14 +05:30
2019-10-17 00:43:26 +05:30
call random_number ( v3_1 )
call random_number ( v3_2 )
call random_number ( v3_3 )
call random_number ( v3_4 )
2020-01-29 18:31:14 +05:30
2019-10-17 00:43:26 +05:30
if ( dNeq ( abs ( dot_product ( math_cross ( v3_1 - v3_4 , v3_2 - v3_4 ) , v3_3 - v3_4 ) ) / 6.0 , &
math_volTetrahedron ( v3_1 , v3_2 , v3_3 , v3_4 ) , tol = 1.0e-12_pReal ) ) &
2020-09-13 15:39:32 +05:30
error stop 'math_volTetrahedron'
2020-01-29 18:31:14 +05:30
2019-09-22 09:55:55 +05:30
call random_number ( t33 )
if ( dNeq ( math_det33 ( math_symmetric33 ( t33 ) ) , math_detSym33 ( math_symmetric33 ( t33 ) ) , tol = 1.0e-12_pReal ) ) &
2020-09-13 15:39:32 +05:30
error stop 'math_det33/math_detSym33'
2020-03-01 14:17:20 +05:30
2020-09-13 14:48:09 +05:30
if ( any ( dNeq0 ( math_eye ( 3 ) , math_inv33 ( math_I3 ) ) ) ) &
2020-09-13 15:39:32 +05:30
error stop 'math_inv33(math_I3)'
2020-01-29 18:31:14 +05:30
2019-09-22 19:11:55 +05:30
do while ( abs ( math_det33 ( t33 ) ) < 1.0e-9_pReal )
call random_number ( t33 )
enddo
2020-09-13 14:48:09 +05:30
if ( any ( dNeq0 ( matmul ( t33 , math_inv33 ( t33 ) ) - math_eye ( 3 ) , tol = 1.0e-9_pReal ) ) ) &
2020-09-13 15:39:32 +05:30
error stop 'math_inv33'
2020-01-29 18:31:14 +05:30
2019-09-22 19:11:55 +05:30
call math_invert33 ( t33_2 , det , e , t33 )
2020-09-13 14:48:09 +05:30
if ( any ( dNeq0 ( matmul ( t33 , t33_2 ) - math_eye ( 3 ) , tol = 1.0e-9_pReal ) ) . or . e ) &
2020-09-13 15:39:32 +05:30
error stop 'math_invert33: T:T^-1 != I'
2019-09-22 19:11:55 +05:30
if ( dNeq ( det , math_det33 ( t33 ) , tol = 1.0e-12_pReal ) ) &
2020-09-13 15:39:32 +05:30
error stop 'math_invert33 (determinant)'
2020-01-29 18:31:14 +05:30
2019-09-22 19:11:55 +05:30
call math_invert ( t33_2 , e , t33 )
2020-09-13 14:48:09 +05:30
if ( any ( dNeq0 ( matmul ( t33 , t33_2 ) - math_eye ( 3 ) , tol = 1.0e-9_pReal ) ) . or . e ) &
2020-09-13 15:39:32 +05:30
error stop 'math_invert t33'
2020-01-29 18:31:14 +05:30
2020-03-10 12:24:43 +05:30
do while ( math_det33 ( t33 ) < 1.0e-2_pReal ) ! O(det(F)) = 1
call random_number ( t33 )
enddo
2020-03-15 18:51:11 +05:30
t33_2 = math_rotationalPart ( transpose ( t33 ) )
t33 = math_rotationalPart ( t33 )
if ( any ( dNeq0 ( matmul ( t33_2 , t33 ) - math_I3 , tol = 1.0e-10_pReal ) ) ) &
2020-09-13 15:39:32 +05:30
error stop 'math_rotationalPart'
2019-09-22 19:11:55 +05:30
2020-01-29 18:31:14 +05:30
call random_number ( r )
d = int ( r * 5.0_pReal ) + 1
2020-09-13 14:48:09 +05:30
txx = math_eye ( d )
2020-01-29 18:31:14 +05:30
allocate ( txx_2 ( d , d ) )
call math_invert ( txx_2 , e , txx )
if ( any ( dNeq0 ( txx_2 , txx ) . or . e ) ) &
2020-09-13 15:39:32 +05:30
error stop 'math_invert(txx)/math_eye'
2019-09-22 09:55:55 +05:30
2020-03-01 14:17:20 +05:30
call math_invert ( t99_2 , e , t99 ) ! not sure how likely it is that we get a singular matrix
2020-09-13 14:48:09 +05:30
if ( any ( dNeq0 ( matmul ( t99_2 , t99 ) - math_eye ( 9 ) , tol = 1.0e-9_pReal ) ) . or . e ) &
2020-09-13 15:39:32 +05:30
error stop 'math_invert(t99)'
2019-09-22 09:55:55 +05:30
2019-09-23 19:00:21 +05:30
if ( any ( dNeq ( math_clip ( [ 4.0_pReal , 9.0_pReal ] , 5.0_pReal , 6.5_pReal ) , [ 5.0_pReal , 6.5_pReal ] ) ) ) &
2020-09-13 15:39:32 +05:30
error stop 'math_clip'
2019-09-22 19:11:55 +05:30
2019-09-23 10:11:00 +05:30
if ( math_factorial ( 10 ) / = 3628800 ) &
2020-09-13 15:39:32 +05:30
error stop 'math_factorial'
2019-09-23 10:11:00 +05:30
if ( math_binomial ( 49 , 6 ) / = 13983816 ) &
2020-09-13 15:39:32 +05:30
error stop 'math_binomial'
2019-09-23 10:11:00 +05:30
2020-03-02 03:30:06 +05:30
ijk = cshift ( [ 1 , 2 , 3 ] , int ( r * 1.0e2_pReal ) )
if ( dNeq ( math_LeviCivita ( ijk ( 1 ) , ijk ( 2 ) , ijk ( 3 ) ) , + 1.0_pReal ) ) &
2020-09-13 15:39:32 +05:30
error stop 'math_LeviCivita(even)'
2020-03-02 03:30:06 +05:30
ijk = cshift ( [ 3 , 2 , 1 ] , int ( r * 2.0e2_pReal ) )
if ( dNeq ( math_LeviCivita ( ijk ( 1 ) , ijk ( 2 ) , ijk ( 3 ) ) , - 1.0_pReal ) ) &
2020-09-13 15:39:32 +05:30
error stop 'math_LeviCivita(odd)'
2020-03-02 03:30:06 +05:30
ijk = cshift ( [ 2 , 2 , 1 ] , int ( r * 2.0e2_pReal ) )
2020-09-13 15:39:32 +05:30
if ( dNeq0 ( math_LeviCivita ( ijk ( 1 ) , ijk ( 2 ) , ijk ( 3 ) ) ) ) &
error stop 'math_LeviCivita'
2020-03-02 03:30:06 +05:30
2020-05-16 20:35:03 +05:30
end subroutine selfTest
2019-09-22 09:55:55 +05:30
2016-02-26 21:05:55 +05:30
end module math