Merge remote-tracking branch 'origin/development' into order4-polynomial
This commit is contained in:
Martin Diehl
commit
5b8e5591ed
|
|
@ -1966,8 +1966,8 @@ end function buildCoordinateSystem
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Helper function to define transformation systems
|
||||
! Needed to calculate Schmid matrix and rotated stiffness matrices.
|
||||
! @details: use c/a for cF -> cI transformation
|
||||
! use a_cX for cF -> hP transformation
|
||||
! @details: use c/a for cF -> hP transformation
|
||||
! use a_cX for cF -> cI transformation
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_cF,a_cI)
|
||||
|
||||
|
|
|
|||
45
src/math.f90
45
src/math.f90
|
|
@ -127,8 +127,10 @@ pure recursive subroutine math_sort(a, istart, iend, sortDim)
|
|||
|
||||
integer, dimension(:,:), intent(inout) :: a
|
||||
integer, intent(in),optional :: istart,iend, sortDim
|
||||
|
||||
integer :: ipivot,s,e,d
|
||||
|
||||
|
||||
if (present(istart)) then
|
||||
s = istart
|
||||
else
|
||||
|
|
@ -164,9 +166,11 @@ pure recursive subroutine math_sort(a, istart, iend, sortDim)
|
|||
integer, dimension(:,:), intent(inout) :: a
|
||||
integer, intent(out) :: p ! Pivot element
|
||||
integer, intent(in) :: istart,iend,sort
|
||||
integer, dimension(size(a,1)) :: tmp
|
||||
|
||||
integer, dimension(size(a,1)) :: tmp
|
||||
integer :: i,j
|
||||
|
||||
|
||||
do
|
||||
! find the first element on the right side less than or equal to the pivot point
|
||||
do j = iend, istart, -1
|
||||
|
|
@ -204,8 +208,10 @@ pure function math_expand(what,how)
|
|||
real(pReal), dimension(:), intent(in) :: what
|
||||
integer, dimension(:), intent(in) :: how
|
||||
real(pReal), dimension(sum(how)) :: math_expand
|
||||
|
||||
integer :: i
|
||||
|
||||
|
||||
if (sum(how) == 0) return
|
||||
|
||||
do i = 1, size(how)
|
||||
|
|
@ -221,9 +227,11 @@ end function math_expand
|
|||
pure function math_range(N)
|
||||
|
||||
integer, intent(in) :: N !< length of range
|
||||
integer :: i
|
||||
integer, dimension(N) :: math_range
|
||||
|
||||
integer :: i
|
||||
|
||||
|
||||
math_range = [(i,i=1,N)]
|
||||
|
||||
end function math_range
|
||||
|
|
@ -235,9 +243,11 @@ end function math_range
|
|||
pure function math_eye(d)
|
||||
|
||||
integer, intent(in) :: d !< tensor dimension
|
||||
integer :: i
|
||||
real(pReal), dimension(d,d) :: math_eye
|
||||
|
||||
integer :: i
|
||||
|
||||
|
||||
math_eye = 0.0_pReal
|
||||
do i=1,d
|
||||
math_eye(i,i) = 1.0_pReal
|
||||
|
|
@ -302,6 +312,7 @@ real(pReal) pure function math_delta(i,j)
|
|||
|
||||
integer, intent (in) :: i,j
|
||||
|
||||
|
||||
math_delta = merge(0.0_pReal, 1.0_pReal, i /= j)
|
||||
|
||||
end function math_delta
|
||||
|
|
@ -315,6 +326,7 @@ pure function math_cross(A,B)
|
|||
real(pReal), dimension(3), intent(in) :: A,B
|
||||
real(pReal), dimension(3) :: math_cross
|
||||
|
||||
|
||||
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) ]
|
||||
|
|
@ -329,6 +341,7 @@ pure function math_outer(A,B)
|
|||
|
||||
real(pReal), dimension(:), intent(in) :: A,B
|
||||
real(pReal), dimension(size(A,1),size(B,1)) :: math_outer
|
||||
|
||||
integer :: i,j
|
||||
|
||||
|
||||
|
|
@ -351,6 +364,7 @@ real(pReal) pure function math_inner(A,B)
|
|||
real(pReal), dimension(:), intent(in) :: A
|
||||
real(pReal), dimension(size(A,1)), intent(in) :: B
|
||||
|
||||
|
||||
math_inner = sum(A*B)
|
||||
|
||||
end function math_inner
|
||||
|
|
@ -363,6 +377,7 @@ real(pReal) pure function math_tensordot(A,B)
|
|||
|
||||
real(pReal), dimension(3,3), intent(in) :: A,B
|
||||
|
||||
|
||||
math_tensordot = sum(A*B)
|
||||
|
||||
end function math_tensordot
|
||||
|
|
@ -376,6 +391,7 @@ pure function math_mul3333xx33(A,B)
|
|||
real(pReal), dimension(3,3,3,3), intent(in) :: A
|
||||
real(pReal), dimension(3,3), intent(in) :: B
|
||||
real(pReal), dimension(3,3) :: math_mul3333xx33
|
||||
|
||||
integer :: i,j
|
||||
|
||||
|
||||
|
|
@ -395,11 +411,12 @@ end function math_mul3333xx33
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
pure function math_mul3333xx3333(A,B)
|
||||
|
||||
integer :: 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
|
||||
|
||||
integer :: i,j,k,l
|
||||
|
||||
|
||||
#ifndef __INTEL_COMPILER
|
||||
do concurrent(i=1:3, j=1:3, k=1:3, l=1:3)
|
||||
|
|
@ -424,6 +441,7 @@ pure function math_exp33(A,n)
|
|||
real(pReal) :: invFac
|
||||
integer :: n_,i
|
||||
|
||||
|
||||
if (present(n)) then
|
||||
n_ = n
|
||||
else
|
||||
|
|
@ -456,6 +474,7 @@ pure function math_inv33(A)
|
|||
real(pReal) :: DetA
|
||||
logical :: error
|
||||
|
||||
|
||||
call math_invert33(math_inv33,DetA,error,A)
|
||||
if (error) math_inv33 = 0.0_pReal
|
||||
|
||||
|
|
@ -474,6 +493,7 @@ pure subroutine math_invert33(InvA, DetA, error, A)
|
|||
logical, intent(out) :: error
|
||||
real(pReal), dimension(3,3), intent(in) :: A
|
||||
|
||||
|
||||
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)
|
||||
|
|
@ -504,7 +524,7 @@ end subroutine math_invert33
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
pure function math_invSym3333(A)
|
||||
|
||||
real(pReal),dimension(3,3,3,3) :: math_invSym3333
|
||||
real(pReal),dimension(3,3,3,3) :: math_invSym3333
|
||||
|
||||
real(pReal),dimension(3,3,3,3),intent(in) :: A
|
||||
|
||||
|
|
@ -513,6 +533,7 @@ pure function math_invSym3333(A)
|
|||
real(pReal), dimension(6*6) :: work
|
||||
integer :: ierr_i, ierr_f
|
||||
|
||||
|
||||
temp66 = math_sym3333to66(A)
|
||||
call dgetrf(6,6,temp66,6,ipiv6,ierr_i)
|
||||
call dgetri(6,temp66,6,ipiv6,work,size(work,1),ierr_f)
|
||||
|
|
@ -538,6 +559,7 @@ pure subroutine math_invert(InvA, error, A)
|
|||
real(pReal), dimension(size(A,1)**2) :: work
|
||||
integer :: ierr
|
||||
|
||||
|
||||
invA = A
|
||||
call dgetrf(size(A,1),size(A,1),invA,size(A,1),ipiv,ierr)
|
||||
error = (ierr /= 0)
|
||||
|
|
@ -555,6 +577,7 @@ pure function math_symmetric33(m)
|
|||
real(pReal), dimension(3,3) :: math_symmetric33
|
||||
real(pReal), dimension(3,3), intent(in) :: m
|
||||
|
||||
|
||||
math_symmetric33 = 0.5_pReal * (m + transpose(m))
|
||||
|
||||
end function math_symmetric33
|
||||
|
|
@ -568,6 +591,7 @@ pure function math_skew33(m)
|
|||
real(pReal), dimension(3,3) :: math_skew33
|
||||
real(pReal), dimension(3,3), intent(in) :: m
|
||||
|
||||
|
||||
math_skew33 = m - math_symmetric33(m)
|
||||
|
||||
end function math_skew33
|
||||
|
|
@ -581,6 +605,7 @@ pure function math_spherical33(m)
|
|||
real(pReal), dimension(3,3) :: math_spherical33
|
||||
real(pReal), dimension(3,3), intent(in) :: m
|
||||
|
||||
|
||||
math_spherical33 = math_I3 * math_trace33(m)/3.0_pReal
|
||||
|
||||
end function math_spherical33
|
||||
|
|
@ -594,6 +619,7 @@ pure function math_deviatoric33(m)
|
|||
real(pReal), dimension(3,3) :: math_deviatoric33
|
||||
real(pReal), dimension(3,3), intent(in) :: m
|
||||
|
||||
|
||||
math_deviatoric33 = m - math_spherical33(m)
|
||||
|
||||
end function math_deviatoric33
|
||||
|
|
@ -606,6 +632,7 @@ real(pReal) pure function math_trace33(m)
|
|||
|
||||
real(pReal), dimension(3,3), intent(in) :: m
|
||||
|
||||
|
||||
math_trace33 = m(1,1) + m(2,2) + m(3,3)
|
||||
|
||||
end function math_trace33
|
||||
|
|
@ -618,6 +645,7 @@ real(pReal) pure function math_det33(m)
|
|||
|
||||
real(pReal), dimension(3,3), intent(in) :: m
|
||||
|
||||
|
||||
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))
|
||||
|
|
@ -632,6 +660,7 @@ real(pReal) pure function math_detSym33(m)
|
|||
|
||||
real(pReal), dimension(3,3), intent(in) :: m
|
||||
|
||||
|
||||
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)
|
||||
|
||||
|
|
@ -761,6 +790,7 @@ pure function math_99to3333(m99)
|
|||
|
||||
integer :: i,j
|
||||
|
||||
|
||||
#ifndef __INTEL_COMPILER
|
||||
do concurrent(i=1:9, j=1:9)
|
||||
math_99to3333(MAPPLAIN(1,i),MAPPLAIN(2,i),MAPPLAIN(1,j),MAPPLAIN(2,j)) = m99(i,j)
|
||||
|
|
@ -1012,6 +1042,7 @@ pure subroutine math_eigh33(w,v,m)
|
|||
real(pReal) :: T, U, norm, threshold
|
||||
logical :: error
|
||||
|
||||
|
||||
w = math_eigvalsh33(m)
|
||||
|
||||
v(1:3,2) = [ m(1, 2) * m(2, 3) - m(1, 3) * m(2, 2), &
|
||||
|
|
@ -1066,6 +1097,7 @@ pure function math_rotationalPart(F) result(R)
|
|||
I_F ! first two invariants of F
|
||||
real(pReal) :: x,Phi
|
||||
|
||||
|
||||
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)))]
|
||||
|
|
@ -1105,6 +1137,7 @@ pure function math_eigvalsh(m)
|
|||
integer :: ierr
|
||||
real(pReal), dimension(size(m,1)**2) :: work
|
||||
|
||||
|
||||
m_= m ! copy matrix to input (will be destroyed)
|
||||
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)
|
||||
|
|
@ -1126,6 +1159,7 @@ pure function math_eigvalsh33(m)
|
|||
real(pReal) :: P, Q, rho, phi
|
||||
real(pReal), parameter :: TOL=1.e-14_pReal
|
||||
|
||||
|
||||
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
|
||||
|
||||
P = I(2)-I(1)**2/3.0_pReal ! different from http://arxiv.org/abs/physics/0610206 (this formulation was in DAMASK)
|
||||
|
|
@ -1157,6 +1191,7 @@ pure function math_invariantsSym33(m)
|
|||
real(pReal), dimension(3,3), intent(in) :: m
|
||||
real(pReal), dimension(3) :: math_invariantsSym33
|
||||
|
||||
|
||||
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)
|
||||
|
|
|
|||
Loading…
Reference in New Issue