-
Dave Whipp authoredDave Whipp authored
ps.f90 6.38 KiB
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ROUTINE Nov. 2006 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine PS (cc11,cc22,cc33,cc12,cc23,cc31,l1,l2,l3, &
n11,n12,n13, &
n21,n22,n23, &
n31,n32,n33)
implicit none
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! calculates principal directions and eigenvalues of symmetrical
! tensor cc11,cc22,cc33,cc12,cc23,cc31
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
double precision cc11,cc22,cc33,cc12,cc23,cc31
double precision l1,l2,l3
double precision n11,n12,n13,n21,n22,n23,n31,n32,n33
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
double precision a1,a2,a3,q
double precision r, theta, a13, xnorm
double precision q1,q2,q3
double precision x1,x2,x3
double precision pi3,pi23,pi43
double precision c11,c22,c33,c12,c23,c13
double precision mi(3),mimax
integer im
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
pi3=3.141592654d0/3.d0
pi23=pi3*2.d0
pi43=pi3*4.d0
c11=cc11
c22=cc22
c33=cc33
c12=cc12
c23=cc23
c13=cc31
! find the stress tensor eigenvalues (principal stresses)
a3=-c11*c22*c33-c12*c23*c13-c12*c23*c13 &
+c13*c22*c13+c12*c12*c33+c11*c23*c23a2=c11*c22+c22*c33+c33*c11-c13*c13-c23*c23-c12*c12
a1=-c11-c22-c33
q=(a1*a1-3.d0*a2)/9.d0
r=(2.d0*a1*a1*a1-9.d0*a1*a2+27.d0*a3)/54.d0
q3=q*q*q
! check for multiple eigenvalues
if (q3-r*r.le.0.d0) then
l1=c11
l2=c22
l3=c33
else
theta=acos(r/sqrt(q3))/3.d0
q2=sqrt(q)
a13=a1/3.d0
l1=-2.d0*q2*cos(theta)-a13
l2=-2.d0*q2*cos(theta+pi23)-a13
l3=-2.d0*q2*cos(theta+pi43)-a13
endif
! first eigenvector
c11=c11-l1
c22=c22-l1
c33=c33-l1
mi(1)=c22*c33-c23*c23
mi(2)=c33*c11-c13*c13
mi(3)=c11*c22-c12*c12
im=2
if (abs(mi(1)).gt.abs(mi(2))) im=1
if (abs(mi(im)).lt.abs(mi(3))) im=3
mimax=mi(im)
if (im.eq.1) then
x2=(c13*c23-c12*c33)/mimax
x3=(c12*c23-c13*c22)/mimax
xnorm=sqrt(1.d0+x2*x2+x3*x3)
n11=1.d0/xnorm
n12=x2/xnorm
n13=x3/xnorm
elseif (im.eq.2) then
x3=(c12*c13-c23*c11)/mimax
x1=(c13*c23-c12*c33)/mimax
xnorm=sqrt(1.d0+x3*x3+x1*x1)
n11=x1/xnorm
n12=1.d0/xnorm
n13=x3/xnorm
elseif (im.eq.3) then
x1=(c23*c12-c13*c22)/mimax
x2=(c13*c12-c23*c11)/mimax
xnorm=sqrt(1.d0+x1*x1+x2*x2)
n11=x1/xnorm
n12=x2/xnorm
n13=1.d0/xnorm
endif
! second eigenvector
c11=c11+l1-l2
c22=c22+l1-l2
c33=c33+l1-l2 mi(1)=c22*c33-c23*c23
mi(2)=c33*c11-c13*c13
mi(3)=c11*c22-c12*c12
im=2
if (abs(mi(1)).gt.abs(mi(2))) im=1
if (abs(mi(im)).lt.abs(mi(3))) im=3
mimax=mi(im)
if (im.eq.1) then
x2=(c13*c23-c12*c33)/mimax
x3=(c12*c23-c13*c22)/mimax
xnorm=sqrt(1.d0+x2*x2+x3*x3)
n21=1.d0/xnorm
n22=x2/xnorm
n23=x3/xnorm
elseif (im.eq.2) then
x3=(c12*c13-c23*c11)/mimax
x1=(c13*c23-c12*c33)/mimax
xnorm=sqrt(1.d0+x3*x3+x1*x1)
n21=x1/xnorm
n22=1.d0/xnorm
n23=x3/xnorm
elseif (im.eq.3) then
x1=(c23*c12-c13*c22)/mimax
x2=(c13*c12-c23*c11)/mimax
xnorm=sqrt(1.d0+x1*x1+x2*x2)
n21=x1/xnorm
n22=x2/xnorm
n23=1.d0/xnorm
endif
! third eigenvector
c11=c11+l2-l3
c22=c22+l2-l3
c33=c33+l2-l3
mi(1)=c22*c33-c23*c23
mi(2)=c33*c11-c13*c13
mi(3)=c11*c22-c12*c12
im=2
if (abs(mi(1)).gt.abs(mi(2))) im=1
if (abs(mi(im)).lt.abs(mi(3))) im=3
mimax=mi(im)
if (im.eq.1) then
x2=(c13*c23-c12*c33)/mimax
x3=(c12*c23-c13*c22)/mimax
xnorm=sqrt(1.+x2*x2+x3*x3)
n31=1.d0/xnorm
n32=x2/xnorm
n33=x3/xnorm
elseif (im.eq.2) then
x3=(c12*c13-c23*c11)/mimax
x1=(c13*c23-c12*c33)/mimax
xnorm=sqrt(1.d0+x3*x3+x1*x1)
n31=x1/xnorm
n32=1.d0/xnorm
n33=x3/xnorm
elseif (im.eq.3) then
x1=(c23*c12-c13*c22)/mimax
x2=(c13*c12-c23*c11)/mimax
xnorm=sqrt(1.d0+x1*x1+x2*x2)
n31=x1/xnorm
n32=x2/xnorm
n33=1.d0/xnorm
endif
return
end
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------