efg2d.f95

用无网格EFG方法求解悬臂梁程序

Subroutine GaussEqSolver_Sym(n,ma,a,b,ep,kwji)
!------------------------------------------------------------------
! Solve sysmmetric linear equation ax=b by using Gauss elimination.
! If kwji=1, no solution;if kwji=0,has solution
! Input--n,ma,a(ma,n),b(n),ep,
! Output--b,kwji
!------------------------------------------------------------------
implicit real*8 (a-h,o-z)
dimension a(ma,n),b(n),m(n+1)
do 10 i=1,n
10 m(i)=i
do 120 k=1,n
p=0.0
do 20 i=k,n
do 20 j=k,n
if(dabs(a(i,j)).gt.dabs(p)) then
p=a(i,j)
io=i
jo=j
endif
20 continue
if(dabs(p)-ep) 30,30,35
30 kwji=1
return
35 continue
if(jo.eq.k) go to 45
do 40 i=1,n
t=a(i,jo)
a(i,jo)=a(i,k)
a(i,k)=t
40 continue
j=m(k)
m(k)=m(jo)
m(jo)=j
45 if(io.eq.k) go to 55
do 50 j=k,n
t=a(io,j)
a(io,j)=a(k,j)
a(k,j)=t
50 continue
t=b(io)
b(io)=b(k)
b(k)=t
55 p=1./p
in=n-1
if(k.eq.n) go to 65
do 60 j=k,in
60 a(k,j+1)=a(k,j+1)*p
65 b(k)=b(k)*p
if(k.eq.n) go to 120
do 80 i=k,in
do 70 j=k,in
70 a(i+1,j+1)=a(i+1,j+1)-a(i+1,k)*a(k,j+1)
80 b(i+1)=b(i+1)-a(i+1,k)*b(k)
120 continue
do 130 i1=2,n
i=n+1-i1
do 130 j=i,in
130 b(i)=b(i)-a(i,j+1)*b(j+1)
do 140 k=1,n
i=m(k)
140 a(1,i)=b(k)
do 150 k=1,n
150 b(k)=a(1,k)
kwji=0
return
END




SUBROUTINE Weight_W1(dif,nv,ds,w,nx,ndex,numnode)
!------------------------------------------------------------------
! Cubic spline weight function
! input--dif,nv,ds,nx,ndex,numnode
! output--w
! from 1 to 10 column of w denotes w,dwdx,dwdy,dwdxx,dwdxy,dwdyy
!------------------------------------------------------------------
implicit real*8 (a-h,o-z)
dimension dif(nx,ndex),nv(numnode),ds(nx,numnode),w(ndex,10)
ep=1.0e-20
do 10 i=1,ndex
nn=nv(i)
difx=dif(1,i)
dify=dif(2,i)
if(dabs(difx).le.ep) then
drdx=0.
else
drdx=(difx/dabs(difx))/ds(1,nn)
end if
if (dabs(dify).le.ep) then
drdy=0.
else
drdy=(dify/dabs(dify))/ds(2,nn)
end if
rx=dabs(dif(1,i))/ds(1,nn)
ry=dabs(dif(2,i))/ds(2,nn)
if(rx.gt.0.5) then
wx=(4./3.)-4.*rx+4.*rx*rx-(4./3.)*rx*rx*rx
dwxdx=(-4.+8.*rx-4.*rx*rx)*drdx
dwxdxx=(8.-8.*rx)*drdx*drdx
dwxdxxx=(-8.)*drdx*drdx*drdx
else if(rx.le.0.5) then
wx=(2./3.)-4.*rx*rx+4.*rx*rx*rx
dwxdx=(-8.*rx+12.*rx*rx)*drdx
dwxdxx=(-8.+24.*rx)*drdx*drdx
dwxdxxx=(24.)*drdx*drdx*drdx
end if
if(ry.gt.0.5) then
wy=(4./3.)-4.*ry+4.*ry*ry-(4./3.)*ry*ry*ry
dwydy=(-4.+8.*ry-4.*ry*ry)*drdy
dwydyy=(8.-8.*ry)*drdy*drdy
dwydyyy=(-8.)*drdy*drdy*drdy
else if(ry.le.0.5) then
wy=(2./3.)-4.*ry*ry+4.*ry*ry*ry
dwydy=(-8.*ry+12.*ry*ry)*drdy
dwydyy=(-8.+24.*ry)*drdy*drdy
dwydyyy=(24.)*drdy*drdy*drdy
end if
w(i,1)=wx*wy
w(i,2)=wy*dwxdx
w(i,3)=wx*dwydy
w(i,4)=wy*dwxdxx
w(i,5)=dwxdx*dwydy
w(i,6)=wx*dwydyy
w(i,7)=wy*dwxdxxx
w(i,8)=dwxdxx*dwydy
w(i,9)=dwxdx*dwydyy
w(i,10)=wx*dwydyyy
10 continue
return
end





SUBROUTINE Compute_Basis(gpos,gp,nx,mm)
!------------------------------------------------------------------
! Compute basis functions and their derivatives
! Input-gpos,nx,mm
! Output-gp
! From 1 to 10 columns of gp: p,dpdx,dpdy,dpdxx,dpdxy,dpdyy,
! dpdxxx,dpdxxy,dpdxyy,dpdyyy
!------------------------------------------------------------------
implicit real*8 (a-h,o-z)
dimension gpos(nx),gp(10,mm)
do i=1,mm
do j=1,10
gp(i,j)=0.0
enddo
enddo
gp(1,1)=1.0
gp(1,2)=gpos(1)
gp(1,3)=gpos(2)
gp(1,4)=gpos(1)*gpos(1)
gp(1,5)=gpos(1)*gpos(2)
gp(1,6)=gpos(2)*gpos(2)
gp(2,2)=1.0
gp(2,4)=2.0*gpos(1)
gp(2,5)=gpos(2)
gp(3,3)=1.0
gp(3,5)=gpos(1)
gp(3,6)=2.0*gpos(2)
gp(4,4)=2.0
gp(5,5)=1.0
gp(6,6)=2.0
return
end



SUBROUTINE Compute_AB(gpos,x,nv,ds,a,b,nx,numnode,ndex,mm)
!------------------------------------------------------------------
! Compute A matrix and B matrix and their derivatives
! input--gpos,x,nv,dm,nx,numnode,ndex,mm
! output--a,b
! From 1 to 10 of the third dimension of a denotes
! a,dax,day,daxx,daxy,dayy, dadxxx,dadxxy,dadxyy,dadyyy
! From 1 to 10 of the third dimension of b denotes
! b,dbx,dby,dbxx,dbxy,dbyy,dbdxxx,dbdxxy,dbdxyy,dbdyyy
!------------------------------------------------------------------
implicit real*8 (a-h,o-z)
dimension gpos(nx),x(nx,numnode),nv(numnode),ds(nx,numnode)
dimension a(mm,mm,10),b(mm,ndex,10)
dimension xv(nx,ndex),dif(nx,ndex),w(ndex,10),p(6,ndex),pp(mm,mm)
do i=1,ndex
nn=nv(i)
xv(1,i)=x(1,nn)
xv(2,i)=x(2,nn)
p(1,i)=1.0
p(2,i)=xv(1,i)
p(3,i)=xv(2,i)
p(4,i)=xv(1,i)*xv(1,i)
p(5,i)=xv(1,i)*xv(2,i)
p(6,i)=xv(2,i)*xv(2,i)
dif(1,i)=gpos(1)-xv(1,i)
dif(2,i)=gpos(2)-xv(2,i)
enddo
call Weight_W1(dif,nv,ds,w,nx,ndex,numnode)
! ************* Compute b and its derivatives
do 20 ii=1,mm
do 20 jj=1,ndex
do 20 kk=1,10
b(ii,jj,kk)=p(ii,jj)*w(jj,kk)
20 continue
! ************* Compute a and its derivatives
do 25 ie=1,mm
do 25 je=1,mm
do 25 ke=1,10
a(ie,je,ke)=0.
25 continue
do 30 iii=1,ndex
do 40 ik=1,mm
do 40 jk=1,mm
pp(ik,jk)=p(ik,iii)*p(jk,iii)
40 continue
do 50 ikk=1,mm
do 50 jkk=1,mm
do 50 kkk=1,10
a(ikk,jkk,kkk)=a(ikk,jkk,kkk)+w(iii,kkk)*pp(ikk,jkk)
50 continue
30 continue
return
end



SUBROUTINE MLS_ShapeFunc_2D(gpos,x,nv,ds,phi,nx,numnode,ndex,mm)
!------------------------------------------------------------------
! Compute MLS shape functions and their derivatives
! Input--gpos,x,nv,ds,nx,numnode,ndex,mm
! Output--phi
! From 1 to 10 of the two dimension of phi denotes
! phi,dphix,dphiy,dphixx,dphixy,dphiyy
! dphidxxx,dphidxxy, dphidxyy, dphidyyy
!------------------------------------------------------------------
implicit real*8 (a-h,o-z)
dimension gpos(nx),x(nx,numnode),nv(numnode)
dimension ds(nx,numnode),xv(nx,ndex)
dimension gp(10,mm),gam(mm,10),a(mm,mm,10)
dimension b(mm,ndex,10),c(mm),aa(mm,mm),phi(10,ndex)
do i1=1,mm
do j1=1,10
gp(j1,i1)=0.0
enddo
enddo
call Compute_Basis(gpos,gp,nx,mm)
call Compute_AB(gpos,x,nv,ds,a,b,nx,numnode,ndex,mm)
ep=1.0e-20
do 10 in=1,mm
c(in)=gp(1,in)
10 continue
do 20 i1=1,mm
do 20 j1=1,mm
aa(i1,j1)=a(i1,j1,1)
20 continue
do i1=1,mm
do j1=1,10
gam(i1,j1)=0.0
enddo
enddo
! ************* Compute gam
call GaussEqSolver_Sym(mm,mm,aa,c,ep,kwji)
21 format(1x,' gam kwji=',i2)
do 25 k1=1,mm
gam(k1,1)=c(k1)
25 continue
! ************* Compute dgamdx
do 30 in=1,mm
c(in)=0.
do 30 jn=1,mm
c(in)=c(in)+a(in,jn,2)*gam(jn,1)
30 continue
do 35 kn=1,mm
c(kn)=gp(2,kn)-c(kn)
35 continue
do 40 i1=1,mm
do 40 j1=1,mm
aa(i1,j1)=a(i1,j1,1)
40 continue
call GaussEqSolver_Sym(mm,mm,aa,c,ep,kwji)
do 45 k1=1,mm
gam(k1,2)=c(k1)
45 continue
! ************* Compute dgamdy
do 50 in=1,mm
c(in)=0.
do 50 jn=1,mm
c(in)=c(in)+a(in,jn,3)*gam(jn,1)
50 continue
do 55 kn=1,mm
c(kn)=gp(3,kn)-c(kn)
55 continue
do 60 i1=1,mm
do 60 j1=1,mm
aa(i1,j1)=a(i1,j1,1)
60 continue
call GaussEqSolver_Sym(mm,mm,aa,c,ep,kwji)
do 65 k1=1,mm
gam(k1,3)=c(k1)
65 continue
! ************* Compute dgamdxx
do 70 in=1,mm
c(in)=0.
do 70 jn=1,mm
c(in)=c(in)+a(in,jn,4)*gam(jn,1)+2.0*a(in,jn,2)*gam(jn,2)
70 continue
do 75 kn=1,mm
c(kn)=gp(4,kn)-c(kn)
75 continue
do 80 i1=1,mm
do 80 j1=1,mm
aa(i1,j1)=a(i1,j1,1)
80 continue
call GaussEqSolver_Sym(mm,mm,aa,c,ep,kwji)
do 85 k1=1,mm
gam(k1,4)=c(k1)
85 continue
! ************* Compute dgamdxy
do 90 in=1,mm
c(in)=0.
do 90 jn=1,mm
c(in)=c(in)+a(in,jn,5)*gam(jn,1)+a(in,jn,2)*gam(jn,3)+ &
a(in,jn,3)*gam(jn,2)
90 continue
do 95 kn=1,mm
c(kn)=gp(5,kn)-c(kn)
95 continue
do 100 i1=1,mm
do 100 j1=1,mm
aa(i1,j1)=a(i1,j1,1)
100 continue
call GaussEqSolver_Sym(mm,mm,aa,c,ep,kwji)
do 105 k1=1,mm
gam(k1,5)=c(k1)
105 continue
! ************* Compute dgamdyy
do 110 in=1,mm
c(in)=0.
do 110 jn=1,mm
c(in)=c(in)+a(in,jn,6)*gam(jn,1)+2.0*a(in,jn,3)*gam(jn,3)
110 continue
do 115 kn=1,mm
c(kn)=gp(6,kn)-c(kn)
115 continue
do 120 i1=1,mm
do 120 j1=1,mm
aa(i1,j1)=a(i1,j1,1)
120 continue
call GaussEqSolver_Sym(mm,mm,aa,c,ep,kwji)
do 125 k1=1,mm
gam(k1,6)=c(k1)
125 continue
! ************* Compute dgamdxxx
do in=1,mm
c(in)=0.
do jn=1,mm
c(in)=c(in)+a(in,jn,7)*gam(jn,1)+3*a(in,jn,4)*gam(jn,2)+ &
3*a(in,jn,2)*gam(jn,4)
enddo
enddo
do kn=1,mm
c(kn)=gp(7,kn)-c(kn)
enddo
do i1=1,mm
do j1=1,mm
aa(i1,j1)=a(i1,j1,1)
enddo
enddo
call GaussEqSolver_Sym(mm,mm,aa,c,ep,kwji)
do k1=1,mm
gam(k1,7)=c(k1)
enddo
! ************* Compute dgamdxxy
do in=1,mm
c(in)=0.
do jn=1,mm
c(in)=c(in)+a(in,jn,8)*gam(jn,1)+ &
a(in,jn,4)*gam(jn,3)+2*a(in,jn,5)*gam(jn,2)+ &
2*a(in,jn,2)*gam(jn,5)+a(in,jn,3)*gam(jn,4)
enddo
enddo
do kn=1,mm
c(kn)=gp(8,kn)-c(kn)
enddo
do i1=1,mm
do j1=1,mm
aa(i1,j1)=a(i1,j1,1)
enddo
enddo
call GaussEqSolver_Sym(mm,mm,aa,c,ep,kwji)
do k1=1,mm
gam(k1,8)=c(k1)
enddo
! ************* Compute dgamdxyy
do in=1,mm
c(in)=0.
do jn=1,mm
c(in)=c(in)+a(in,jn,9)*gam(jn,1)+ &
a(in,jn,6)*gam(jn,2)+2*a(in,jn,5)*gam(jn,3)+ &
2*a(in,jn,3)*gam(jn,5)+a(in,jn,2)*gam(jn,6)
enddo
enddo
do kn=1,mm
c(kn)=gp(9,kn)-c(kn)
enddo
do i1=1,mm
do j1=1,mm
aa(i1,j1)=a(i1,j1,1)
enddo
enddo
call GaussEqSolver_Sym(mm,mm,aa,c,ep,kwji)
do k1=1,mm
gam(k1,9)=c(k1)
enddo
! ************* Compute dgamdyyy
do in=1,mm
c(in)=0.
do jn=1,mm
c(in)=c(in)+a(in,jn,10)*gam(jn,1)+ &
3*a(in,jn,6)*gam(jn,3)+3*a(in,jn,3)*gam(jn,6)
enddo
enddo
do kn=1,mm
c(kn)=gp(10,kn)-c(kn)
enddo
do i1=1,mm
do j1=1,mm
aa(i1,j1)=a(i1,j1,1)
enddo
enddo
call GaussEqSolver_Sym(mm,mm,aa,c,ep,kwji)
do k1=1,mm
gam(k1,10)=c(k1)
enddo
!! ************* Compute Phi and their derivatives
do 130 iph=1,ndex
do iiii=1,10
phi(iiii,iph)=0.0
enddo
do 130 jph=1,mm
phi(1,iph)=phi(1,iph)+gam(jph,1)*b(jph,iph,1)
phi(2,iph)=phi(2,iph)+gam(jph,2)*b(jph,iph,1)+ &
gam(jph,1)*b(jph,iph,2)
phi(3,iph)=phi(3,iph)+gam(jph,3)*b(jph,iph,1)+ &
gam(jph,1)*b(jph,iph,3)
phi(4,iph)=phi(4,iph)+gam(jph,4)*b(jph,iph,1)+ &
2.0*gam(jph,2)*b(jph,iph,2)+gam(jph,1)*b(jph,iph,4)
phi(5,iph)=phi(5,iph)+gam(jph,5)*b(jph,iph,1)+ &
gam(jph,2)*b(jph,iph,3)+gam(jph,3)*b(jph,iph,2)+ &
gam(jph,1)*b(jph,iph,5)
phi(6,iph)=phi(6,iph)+gam(jph,6)*b(jph,iph,1)+ &
2.0*gam(jph,3)*b(jph,iph,3)+gam(jph,1)*b(jph,iph,6)
phi(7,iph)=phi(7,iph)+gam(jph,7)*b(jph,iph,1)+ &
3.0*gam(jph,4)*b(jph,iph,2)+3*gam(jph,2)*b(jph,iph,4)+ &
gam(jph,1)*b(jph,iph,7)
phi(8,iph)=phi(8,iph)+gam(jph,8)*b(jph,iph,1)+ &
2.0*gam(jph,5)*b(jph,iph,2)+2*gam(jph,2)*b(jph,iph,5)+ &
gam(jph,1)*b(jph,iph,8)+gam(jph,4)*b(jph,iph,3)+ &
gam(jph,3)*b(jph,iph,4)
phi(9,iph)=phi(9,iph)+gam(jph,9)*b(jph,iph,1)+ &
2.0*gam(jph,5)*b(jph,iph,3)+2*gam(jph,3)*b(jph,iph,5)+ &
gam(jph,1)*b(jph,iph,9)+gam(jph,6)*b(jph,iph,2)+ &
gam(jph,2)*b(jph,iph,6)
phi(10,iph)=phi(10,iph)+gam(jph,10)*b(jph,iph,1)+ &
3.0*gam(jph,6)*b(jph,iph,3)+3*gam(jph,3)*b(jph,iph,6)+ &
gam(jph,1)*b(jph,iph,10)
130 continue
return
end





!----------------------------------------------------------------------------
! main program--2D FORTRAN 90 CODE-MFree global weak-form methods
! Using square support domain and square background cells
! input file -- input.dat
! output file -- result.dat
! include file -- parameter.h, variable.h
!----------------------------------------------------------------------------
implicit real*8 (a-h,o-z)
include 'parameter.h'
include 'variables.h'
open(ninput,file='Input175_55.dat',status='unknown')
open(noutput,file='result.dat',status='unknown')
! ************* Input data
call input(x,numd,nx,numnode,ndivx,ndivy,ndivxq,ndivyq,&
nconn2,nquado,pAlf,Dmat,ALFs,numcell,numq,noCell,ncn,xc,&
npEBCnum,npEBC,pEBC,npNBCnum,npNBC,pNBC)
numgauss=nquado*nquado !total number of Gauss points in a cell
! ************* Determine sizes of influence domains -- uniform nodal spacing
xspace=xlength/ndivx
yspace=ylength/ndivy
do i=1,numnode
ds(1,i)=alfs*xspace
ds(2,i)=alfs*yspace
enddo
! ************* Coefficients of Gauss points,Weights and Jacobian for a cell
call GaussCoefficient(nquado,gauss)
do ik=1,ng
do jk=1,numgauss
gs(ik,jk)=0
enddo
enddo
do ik=1,2*numd
force(ik)=0.
do jk=1,2*numd
ak(ik,jk)=0.
enddo
enddo
! ************* Loop for background cells
do 10 ibk=1,numcell
write(*,*)'Cell No.=',ibk
! ************* Set Gauss points for this cell
call CellGaussPoints(ibk,numcell,nquado,numq,numgauss, &
xc,noCell,gauss,gs)
! ************* Loop over Gauss points to assemble discrete equations
do 20 ie=1,numgauss