aselimov/cacmb
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Working continuum metric calculation code

Browse source
This commit is contained in:
Alex Selimov 2020-10-20 01:03:31 -04:00
parent 95e2ad0b4d
commit 6e08517697
8 changed files with 531 additions and 10 deletions
Download patch file Download diff file Expand all files Collapse all files

120
src/functions.f90
View file

@ -271,4 +271,124 @@ END FUNCTION StrDnCase
norm_dis(1:3) = (rk - rl)
norm_dis(4) = norm2(rk-rl)
end function
pure function matinv3(A) result(B)
!! Performs a direct calculation of the inverse of a 3×3 matrix.
real(kind=dp), intent(in) :: A(3,3) !! Matrix
real(kind=dp) :: B(3,3) !! Inverse matrix
real(kind=dp) :: detinv
if(abs(A(1,1)*A(2,2)*A(3,3) - A(1,1)*A(2,3)*A(3,2)&
- A(1,2)*A(2,1)*A(3,3) + A(1,2)*A(2,3)*A(3,1)&
+ A(1,3)*A(2,1)*A(3,2) - A(1,3)*A(2,2)*A(3,1)) < lim_zero) then
B(:,:) = 0
return
else
! Calculate the inverse determinant of the matrix
detinv = 1/(A(1,1)*A(2,2)*A(3,3) - A(1,1)*A(2,3)*A(3,2)&
- A(1,2)*A(2,1)*A(3,3) + A(1,2)*A(2,3)*A(3,1)&
+ A(1,3)*A(2,1)*A(3,2) - A(1,3)*A(2,2)*A(3,1))
! Calculate the inverse of the matrix
B(1,1) = +detinv * (A(2,2)*A(3,3) - A(2,3)*A(3,2))
B(2,1) = -detinv * (A(2,1)*A(3,3) - A(2,3)*A(3,1))
B(3,1) = +detinv * (A(2,1)*A(3,2) - A(2,2)*A(3,1))
B(1,2) = -detinv * (A(1,2)*A(3,3) - A(1,3)*A(3,2))
B(2,2) = +detinv * (A(1,1)*A(3,3) - A(1,3)*A(3,1))
B(3,2) = -detinv * (A(1,1)*A(3,2) - A(1,2)*A(3,1))
B(1,3) = +detinv * (A(1,2)*A(2,3) - A(1,3)*A(2,2))
B(2,3) = -detinv * (A(1,1)*A(2,3) - A(1,3)*A(2,1))
B(3,3) = +detinv * (A(1,1)*A(2,2) - A(1,2)*A(2,1))
end if
end function
pure function transpose3(A) result(B)
!!Transposes matrix A
real(kind=dp), intent(in) :: A(3,3)
real(kind=dp) :: B(3,3)
integer :: i, j
forall(i =1:3,j=1:3) B(i,j) = A(j,i)
end function transpose3
pure function sqrt3(A) result(B)
!This calculates the square of matrix A fulfilling the equation B*B = A according to
!the algorithm by Franca1989
real(kind=dp), intent(in) :: A(3,3)
real(kind=dp) :: B(3,3)
real(kind=dp) :: Ione, Itwo, Ithree, l, k, phi, Asq(3,3), lambda, Bone, Btwo, Bthree, p
!Step 1 is calculating the invariants of C
Ione = A(1,1) + A(2,2) + A(3,3)
Asq = matmul(A,A)
Itwo = 0.5_dp *(Ione*Ione - (Asq(1,1) + Asq(2,2) + Asq(3,3)))
Ithree = (A(1,1)*A(2,2)*A(3,3) - A(1,1)*A(2,3)*A(3,2)&
- A(1,2)*A(2,1)*A(3,3) + A(1,2)*A(2,3)*A(3,1)&
+ A(1,3)*A(2,1)*A(3,2) - A(1,3)*A(2,2)*A(3,1))
if (Ithree < 0) then
B(:,:)=0.0_dp
return
end if
!Check for an isotropic matrix
k = Ione*Ione - 3*Itwo
if (k < lim_zero) then
lambda = sqrt(Ione/3.0_dp)
B = lambda*identity_mat(3)
else
l = Ione*(Ione*Ione - 9.0_dp/2.0_dp * Itwo) + 27.0_dp/2.0_dp * Ithree
p = l/(k**(1.5_dp))
if (p > 1.0 ) then
B(:,:) = 0.0_dp
return
end if
if ((p< -1).or.(p>1)) then
B(:,:) = 0.0_dp
return
end if
phi = acos(p)
lambda = sqrt(1.0_dp/3.0_dp * (Ione + 2*sqrt(k)*cos(phi/3)))
!Now calculate invariantes of B
Bthree = sqrt(Ithree)
if((-lambda*lambda + Ione + 2*Ithree/lambda) < 0) then
B(:,:) = 0.0_dp
return
end if
Bone = lambda + sqrt(- lambda*lambda + Ione + 2*Ithree/lambda)
Btwo = (Bone*Bone - Ione)/2.0_dp
!Now calculate B
if(abs(Bone*Btwo -Bthree) < lim_zero) then
B(:,:) = 0.0_dp
return
end if
B = 1/(Bone*Btwo - Bthree) *(Bone*Bthree*identity_mat(3) + (Bone*Bone - Btwo)*A - matmul(A,A))
end if
end function sqrt3
pure function permutation(i,j,k) result(e)
!Calculates the permutation tensor
integer, intent(in) :: i,j,k
integer :: e
if ( ((i==1).and.(j==2).and.(k==3)).or. &
((i==2).and.(j==3).and.(k==1)).or. &
((i==3).and.(j==1).and.(k==2))) then
e=1
else if( ((i==2).and.(j==1).and.(k==3)).or. &
((i==1).and.(j==3).and.(k==2)).or. &
((i==3).and.(j==2).and.(k==1))) then
e=-1
else
e=0
end if
end function permutation
end module functions