📄 function.f
字号:
***********************************************************************
* making B matrix
* This subroutine must be performed after subroutine 'makeaa'
***********************************************************************
subroutine bmatrx(ng,mnod,md,nnod,bl,ai,dadx,dldx,bmat)
implicit real*8(a-h,o-z)
dimension ai(6,6), bl(6,mnod)
dimension dadx(6,6,2), dldx(6,mnod,ng)
dimension bmat(3,mnod*ng)
dimension daidx(6,6,2)
dimension p(6), dpdx(6,2)
do 100 i= 1,nnod*2
do 110 j= 1,3
bmat(j,i)= 0.0d0
110 continue
100 continue
do 120 i= 1,md
do 130 j= 1,md
do 140 k= 1,2
daidx(i,j,k)= 0.0d0
140 continue
130 continue
120 continue
call pvec(md,p,dpdx)
do 240 i= 1,2
do 200 m= 1,md
do 210 j= 1,md
do 220 k= 1,md
do 230 l= 1,md
daidx(j,m,i)= daidx(j,m,i) -
$ ai(j,k)*dadx(k,l,i)*ai(l,m)
230 continue
220 continue
210 continue
200 continue
240 continue
do 300 l= 1,nnod
do 310 k= 1,md
do 320 j= 1,md
bmat(1,(l-1)*2+1)= bmat(1,(l-1)*2+1) +
$ dpdx(j,1)*ai(j,k)*bl(k,l) +
$ p(j)*(daidx(j,k,1)*bl(k,l)+ai(j,k)*dldx(k,l,1))
bmat(2,(l-1)*2+2)= bmat(2,(l-1)*2+2) +
$ dpdx(j,2)*ai(j,k)*bl(k,l) +
$ p(j)*(daidx(j,k,2)*bl(k,l)+ai(j,k)*dldx(k,l,2))
320 continue
310 continue
bmat(3,(l-1)*2+1)= bmat(2,(l-1)*2+2)
bmat(3,(l-1)*2+2)= bmat(1,(l-1)*2+1)
300 continue
return
end
*
***********************************************************************
* making interpolated function phi
* This subroutine must be performed after subroutine 'makeaa'
***********************************************************************
subroutine itfunc(ng,mnod,md,nnod,bl,ai,phi)
implicit real*8(a-h,o-z)
dimension ai(6,3), bl(6,mnod)
dimension phi(2,mnod*ng)
dimension p(6),dpdx(6,2)
do 100 i= 1,nnod*ng
do 110 j= 1,ng
phi(j,i)= 0.0d0
110 continue
100 continue
call pvec(md,p,dpdx)
do 200 k= 1,nnod
do 210 j= 1,md
do 220 i= 1,md
phi(1,ng*k-1)= phi(1,ng*k-1) + p(i)*ai(i,j)*bl(j,k)
220 continue
210 continue
200 continue
do 230 k= 1,nnod
phi(2,ng*k)= phi(1,ng*k-1)
230 continue
return
end
*
***********************************************************************
* making matrix A, L, A-1, dA/dX, dL/dX
***********************************************************************
subroutine makeaa(ng,mnod,md,nnod,xi,dm,cm,
$ kk,aa,bl,ai,dadx,dldx)
implicit real*8(a-h,o-z)
dimension xi(ng,mnod)
dimension aa(6,6), ai(6,6), bl(6,mnod)
dimension dadx(6,6,2), dldx(6,mnod,ng)
dimension xxi(2), pi(6), dw(2)
do 100 i= 1,md
do 110 j= 1,md
aa(i,j)= 0.0d0
dadx(i,j,1)= 0.0d0
dadx(i,j,2)= 0.0d0
110 continue
100 continue
do 130 i= 1,nnod
do 140 j= 1,md
bl(j,i)= 0.0d0
dldx(j,i,1)= 0.0d0
dldx(j,i,2)= 0.0d0
140 continue
130 continue
do 200 ip= 1,nnod
xxi(1)= xi(1,ip)
xxi(2)= xi(2,ip)
call pvecaa(md,xxi,pi)
wait= waitf(xxi,dm,cm,kk)
call dwait(xxi,dm,cm,kk,dw)
do 300 i= 1,md
do 310 j= 1,md
pipj= pi(i)*pi(j)
aa(i,j)= aa(i,j) + pipj*wait
dadx(i,j,1)= dadx(i,j,1) + pipj*dw(1)
dadx(i,j,2)= dadx(i,j,2) + pipj*dw(2)
310 continue
300 continue
do 350 i= 1,md
bl(i,ip)= wait*pi(i)
350 continue
do 400 i= 1,2
do 410 j= 1,md
dldx(j,ip,i)= dw(i)*pi(j)
410 continue
400 continue
200 continue
call inverse(md,aa,ai)
return
end
*
**********************************************************************
* Make D-Matrix
**********************************************************************
subroutine dmatrx(ntype,props,nstrs,dmat)
implicit real*8(a-h,o-z)
dimension dmat(4,4), props(7)
character emsg*70
do 10 j= 1,4
do 15 i= 1,4
dmat(i,j)= 0.0d0
15 continue
10 continue
eel= props(1)
por= props(2)
if(ntype.eq.1) then
nstrs= 3
coeff= eel/(1.0d0-por*por)
dmat(1,1)= coeff
dmat(1,2)= coeff*por
dmat(2,1)= dmat(1,2)
dmat(2,2)= coeff
dmat(3,3)= coeff*(1.0d0-por)/2.0d0
elseif(ntype.eq.2) then
nstrs= 3
coeff= eel*(1.0d0-por)/(1.0d0+por)/(1.0d0-2.0d0*por)
dmat(1,1)= coeff
dmat(1,2)= coeff*por/(1.0d0-por)
dmat(2,1)= dmat(1,2)
dmat(2,2)= coeff
dmat(3,3)= eel/(1.0d0+por)/2.0d0
elseif(ntype.eq.3) then
nstrs= 4
coeff=eel*(1.0d0-por)/(1.0d0+por)/(1.0d0-2.0d0*por)
dmat(1,1)= coeff
dmat(1,2)= coeff*por/(1.0d0-por)
dmat(1,3)= dmat(1,2)
dmat(2,1)= dmat(1,2)
dmat(2,2)= coeff
dmat(2,3)= dmat(1,2)
dmat(3,1)= dmat(1,2)
dmat(3,2)= dmat(1,2)
dmat(3,3)= coeff
dmat(4,4)= coeff*(1.0d0-2.0d0*por)/2.0d0/(1.0d0-por)
else
write(7777,*) 'Error Message from subroutine dmat 1 !'
write(7777,*) 'Irregular Analysis Type !
$ The ntype must be 1-4.'
stop
endif
return
end
*
***********************************************************************
* making inverse matrix of a 3*3 matrix
***********************************************************************
subroutine inv33(a,ai)
implicit real*8(a-h,o-z)
dimension a(6,6), ai(6,6)
det = a(1,2)*a(2,3)*a(3,1)+a(1,3)*a(2,1)*a(3,2)
$ +a(1,1)*a(2,2)*a(3,3)-a(1,3)*a(2,2)*a(3,1)
$ -a(1,1)*a(2,3)*a(3,2)-a(1,2)*a(2,1)*a(3,3)
if (det.eq.0) then
write(7777,*) 'Error Message from subroutine inv33 1 ! '
write(7777,*) 'determinant A is zero!!!'
write(7777,600) ((a(i,j),j=1,3),i=1,3)
600 format(3e15.7)
stop
endif
ai(1,1)=(a(2,2)*a(3,3)-a(2,3)*a(3,2))/det
ai(1,2)=(a(1,3)*a(3,2)-a(1,2)*a(3,3))/det
ai(1,3)=(a(1,2)*a(2,3)-a(1,3)*a(2,2))/det
ai(2,1)=(a(2,3)*a(3,1)-a(2,1)*a(3,3))/det
ai(2,2)=(a(1,1)*a(3,3)-a(1,3)*a(3,1))/det
ai(2,3)=(a(1,3)*a(2,1)-a(1,1)*a(2,3))/det
ai(3,1)=(a(2,1)*a(3,2)-a(2,2)*a(3,1))/det
ai(3,2)=(a(1,2)*a(3,1)-a(1,1)*a(3,2))/det
ai(3,3)=(a(1,1)*a(2,2)-a(1,2)*a(2,1))/det
return
end
*
***********************************************************************
* weight function
***********************************************************************
double precision function waitf(xxi,dm,cm,kk)
implicit real*8(a-h,o-z)
dimension xxi(2)
di=sqrt(xxi(1)**2.0d0+xxi(2)**2.0d0)
if(di.gt.dm) then
waitf= 0.0d0
else
waitf = exp(-(di/cm)**(2.0d0*kk))-exp(-(dm/cm)**(2.0d0*kk))
waitf = waitf/(1-exp(-(dm/cm)**(2.0d0*kk)))
endif
return
end
*
***********************************************************************
* partial differencials of the weight function
***********************************************************************
subroutine dwait(xxi,dm,cm,kk,dw)
implicit real*8(a-h,o-z)
dimension xxi(2), dw(2)
di=sqrt(xxi(1)**2.0d0+xxi(2)**2.0d0)
if(di.lt.dm) then
coe1= -((di/cm)**(2*kk))
coe2= -((dm/cm)**(2*kk))
dwdd2k= -exp(coe1)*(cm**(-2*kk))/(1.0d0-exp(coe2))
do 110 i= 1,2
dw(i)= - 2.0d0*kk*(di**(2*kk-2))*xxi(i)*dwdd2k
110 continue
else
dw(1)= 0.0d0
dw(2)= 0.0d0
endif
return
end
*
*********************************************************************
* 0 clear of matrix (real*8)
*********************************************************************
subroutine zerod(npsize,dmatrx)
real*8 dmatrx(npsize)
do 10 i=1,npsize
dmatrx(i)= 0.0d0
10 continue
return
end
*
*********************************************************************
* 0 clear of matrix (integer)
*********************************************************************
subroutine zeroi(npsize,imatrx)
integer imatrx(npsize)
do 10 i=1,npsize
imatrx(i)= 0
10 continue
return
end
*
***********************************************************************
* making inverse matrix of a matrix
***********************************************************************
subroutine inverse(md,aa,ai)
implicit real*8(a-h,o-z)
dimension aa(6,6), ai(6,6)
if(md.eq.3) then
call inv33(aa,ai)
elseif(md.eq.6) then
call invn(md,ai,aa)
else
write(7777,*) 'Error Message from subroutine inverse 1 ! '
write(7777,*) 'md is not an available number!!!'
stop
endif
return
end
*
***********************************************************************
* set p program
***********************************************************************
subroutine pvec(md,p,dpdx)
implicit real*8(a-h,o-z)
dimension p(6), dpdx(6,2)
if(md.eq.3) then
p(1)= 1.0d0
p(2)= 0.0d0
p(3)= 0.0d0
dpdx(1,1)= 0.0d0
dpdx(2,1)= 1.0d0
dpdx(3,1)= 0.0d0
dpdx(1,2)= 0.0d0
dpdx(2,2)= 0.0d0
dpdx(3,2)= 1.0d0
elseif(md.eq.6) then
p(1)= 1.0d0
p(2)= 0.0d0
p(3)= 0.0d0
p(4)= 0.0d0
p(5)= 0.0d0
p(6)= 0.0d0
dpdx(1,1)= 0.0d0
dpdx(2,1)= 1.0d0
dpdx(3,1)= 0.0d0
dpdx(4,1)= 0.0d0
dpdx(5,1)= 0.0d0
dpdx(6,1)= 0.0d0
dpdx(1,2)= 0.0d0
dpdx(2,2)= 0.0d0
dpdx(3,2)= 1.0d0
dpdx(4,2)= 0.0d0
dpdx(5,2)= 0.0d0
dpdx(6,2)= 0.0d0
else
write(7777,*) 'Error Message from subroutine pvec 1 ! '
write(7777,*) 'md is not an available number!!!'
stop
endif
return
end
*
***********************************************************************
* set p program for assemble A matrix
***********************************************************************
subroutine pvecaa(md,xxi,pi)
implicit real*8(a-h,o-z)
dimension xxi(2)
dimension pi(6)
if(md.eq.3) then
pi(1)= 1.0d0
pi(2)= xxi(1)
pi(3)= xxi(2)
elseif(md.eq.6) then
pi(1)= 1.0d0
pi(2)= xxi(1)
pi(3)= xxi(2)
pi(4)= xxi(1) * xxi(1)
pi(5)= xxi(1) * xxi(2)
pi(6)= xxi(2) * xxi(2)
else
write(7777,*) 'Error Message from subroutine pvecaa 1 ! '
write(7777,*) 'md is not an available number!!!'
stop
endif
return
end
*
***********************************************************************
* making inverse matrix of a 6*6 matrix
***********************************************************************
subroutine invn(n,a,b)
implicit real*8(a-h,o-z)
parameter(ipx=6)
parameter(eps=0.0d0)
DIMENSION A(IPX,IPX),B(IPX,IPX),IP(50),IQ(50)
DO 20 I=1,N
IP(I)=0
IQ(I)=0
DO 20 J=1,N
A(I,J)=B(I,J)
20 CONTINUE
DO 25 K=1,N
PMAX=0.0d0
DO 30 I=1,N
IF(IQ(I).EQ.0)THEN
AA=ABS(A(I,K))
IF(AA.GT.PMAX)THEN
PMAX=AA
L=I
END IF
END IF
30 CONTINUE
IF(PMAX.LE.EPS)THEN
write(7777,*) 'Error Message from subroutine invn 1 ! '
write(7777,*) 'matrix inversion error !!!'
stop
END IF
IP(K)=L
IQ(L)=K
ALK=A(L,K)
DO 40 J=1,N
IF(J.NE.K)THEN
A(L,J)=A(L,J)/ALK
W=A(L,J)
DO 50 I=1,N
IF(I.NE.L) THEN
A(I,J)=A(I,J)-W*A(I,K)
END IF
50 CONTINUE
END IF
40 CONTINUE
A(L,K)=1.0d0/A(L,K)
W=A(L,K)
DO 60 I=1,N
IF(I.NE.L)THEN
A(I,K)=-A(I,K)*W
END IF
60 CONTINUE
25 CONTINUE
DO 70 K=1,N
L=IP(K)
IF(L.NE.K)THEN
DO 80 J=1,N
W=A(K,J)
A(K,J)=A(L,J)
A(L,J)=W
80 CONTINUE
M=IQ(K)
DO 90 I=1,N
W=A(I,K)
A(I,K)=A(I,M)
A(I,M)=W
90 CONTINUE
IP(M)=L
IP(K)=K
IQ(L)=M
IQ(K)=K
END IF
70 CONTINUE
RETURN
END
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -