efg2d.f95

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

gpos(1)=gs(1,ie) ! Gauss point x
gpos(2)=gs(2,ie) ! Gauss point y
weight=gs(3,ie) ! weight coefficent
ajac=gs(4,ie) ! Jacobian
! ************* Determine the support domain of Gauss point
ndex=0
call SupportDomain(numnode,nx,gpos,x,ds,ndex,nv)
do ik=1,3*ndex
do jk=1,10
ph(jk,ik)=0.
enddo
enddo
! ************* Construct RPIM shape functions for a Gauss point
call MLS_ShapeFunc_2D(gpos,x,nv,ds,phi,nx,numnode,ndex,mm)
do ik=1,2*ndex
ne(ik)=0
enddo
do ine=1,ndex
n1=2*ine-1
n2=2*ine
ne(n1)=2*nv(ine)-1
ne(n2)=2*nv(ine)
enddo
mbdb=4*ndex*ndex
do kbdb=1,mbdb
GSPk(kbdb)=0.
enddo
! ************* Compute the stiffness matrix for a Gauss point
call PointStiffnessMatrix(ndex,weight,ajac,ph,Dmat,GSPk)
nb=2*ndex
do ikk=1,nb
do jkk=1,nb
m1=ne(ikk)
m2=ne(jkk)
nbdb=(jkk-1)*nb+ikk
ak(m1,m2)=ak(m1,m2)+GSPk(nbdb)
enddo
enddo
20 continue ! end of loop for Gauss points
! ************* Implement natural BC
in=0
jn=0
nn=noCell(3,ibk)
if(xc(1,nn).eq.xlength) in=nn
nn= noCell(4,ibk)
if(xc(1,nn).eq.xlength) jn=nn
if((in.ne.0).and.(jn.ne.0)) then
call naturalBC_distributed(numnode,numq,in,jn, &
alfs,x,xc,ds,gauss,nquado,force)
endif
10 continue ! end of loop for cells
! ************* Boundary conditions: essential BC
write(*,*)' Boundary conditions....'
nak=2*numd
call EssentialBC(numnode,pAlf,alfs,x,ds,ak,force,npEBCnum,npEBC,pEBC)
! ************* Boundary conditions: concentrated natural BC
call NaturalBC_concentrated(x,nx,numnode,force,ds,alfs,npNBCnum,npNBC,pNBC)
nak=2*numd
b=1.d-10
! ************* Solve equation to get the solutions
write(*,*)' Solving....'
call SolverBand(ak,force,2*numnode,2*numd)
nnn=2*numd
do ik=1,nx
do jk=1,numnode
u2(ik,jk)=0.
enddo
enddo
do ik=1,numnode
jk=2*ik-1
u2(1,ik)=force(jk)
u2(2,ik)=force(jk+1)
enddo
! ************* Get the final displacement
call GetDisplacement(x,ds,u2,disp,alfs,nx,numnode)
! ************* Get stress
call GetStress(x,noCell,ds,Dmat,u2,alfs,nx,numnode,numgauss,&
xc,gauss,nquado,ng,numq,numcell, ENORM,Stressnode)
STOP
END




SUBROUTINE 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)
!------------------------------------------------------------------
! Input data from outside
! Output-all variables are output
!------------------------------------------------------------------
implicit real*8 (a-h,o-z)
common/para/xlength,ylength,p,young,anu,aimo
COMMON/rpim/ALFC,DC,Q,nRBF
common /basis/mbasis
CHARACTER*40 NAM
dimension npEBC(3,100),pEBC(2,100),npNBC(3,100),pNBC(2,100)
dimension x(nx,numd),Dmat(3,3),noCell(4,ncn),xc(nx,numd)
read(4,10)nam
read(4,*)xlength,ylength,young,anu,p
read(4,10)nam
read(4,*)numnode,nconn2
read(4,10)nam
read(4,*)ndivx,ndivy
read(4,10)nam
read(4,*)numq,numcell
read(4,10)nam
read(4,*)ndivxq,ndivyq
read(4,10)nam
read(4,*)nquado,pAlf
read(4,10)nam
read(4,*)ALFs
numgauss=nquado*nquado
read(4,10)nam
do i=1,numnode
read(4,*)j,x(1,i),x(2,i)
enddo
read(4,10)nam
do i=1,numq
read(4,*)j,xc(1,i),xc(2,i)
enddo
read(4,10)nam
do j=1,numcell
read(4,*)i,noCell(1,j),noCell(2,j),noCell(3,j),noCell(4,j)
enddo
read(4,10)nam
read(4,*)npEBCnum
read(4,10)nam
do i=1,npEBCnum
read(4,*)npEBC(1,i),npEBC(2,i),npEBC(3,i),pEBC(1,i),pEBC(2,i)
enddo
read(4,10)nam
read(4,*)npNBCnum
read(4,10)nam
do i=1,npNBCnum
read(4,*)npNBC(1,i),npNBC(2,i),npNBC(3,i),pNBC(1,i),pNBC(2,i)
enddo
read(4,10)nam
READ(4,*)nRBF, alfc,dc, q
read(4,10)nam
READ(4,*)mbasis
! ************* Compute material matrix D[] for the plane stress
you=young/(1.-anu*anu)
aimo=(1./12.)*ylength**3
Dmat(1,1)=you
Dmat(1,2)=anu*you
Dmat(1,3)=0.
Dmat(2,1)=anu*you
Dmat(2,2)=you
Dmat(2,3)=0.
Dmat(3,1)=0.
Dmat(3,2)=0.
Dmat(3,3)=0.5*(1.-anu)*you
10 format(a40)
RETURN
END



SUBROUTINE GaussCoefficient(k,v)
!----------------------------------------------------------------------------
! This subroutine returns a matrix with Gauss points and their weights
! input--k: k -- number of Gauss points;
! output--v(2,k): weight matrix of k Gauss points
!---------------------------------------------------------------------------
implicit real*8 (a-h,o-z)
dimension v(2,k)
SELECT CASE (k)
Case (2)
v(1,1)=-.57735
v(1,2)=-v(1,1)
v(2,1)=1.00000
v(2,2)=v(2,1)
Case (3)
v(1,1)=-.77459
v(1,2)=-.00000
v(1,3)=-v(1,1)
v(2,1)=.55555
v(2,2)=.88888
v(2,3)=v(2,1)
Case (4)
v(1,1)=-.86113
v(1,2)=-.33998
v(1,3)=-v(1,2)
v(1,4)=-v(1,1)
v(2,1)=.34785
v(2,2)=.65214
v(2,3)=v(2,2)
v(2,4)=v(2,1)
Case (6)
v(1,1)=-.93246
v(1,2)=-.66120
v(1,3)=-.23861
v(1,4)=-v(1,3)
v(1,5)=-v(1,2)
v(1,6)=-v(1,1)
v(2,1)=.17132
v(2,2)=.36076
v(2,3)=.46791
v(2,4)=v(2,3)
v(2,5)=v(2,2)
v(2,6)=v(2,1)
Case (8)
v(1,1)=-.96028
v(1,2)=-.79666
v(1,3)=-.52553
v(1,4)=-.18343
v(1,5)=-v(1,4)
v(1,6)=-v(1,3)
v(1,7)=-v(1,2)
v(1,8)=-v(1,1)
v(2,1)=.10122
v(2,2)=.22238
v(2,3)=.31370
v(2,4)=.36268
v(2,5)=v(2,4)
v(2,6)=v(2,3)
v(2,7)=v(2,2)
v(2,8)=v(2,1)
end select
RETURN
END


SUBROUTINE CellGaussPoints(ibk,numcell,k,numq,numgauss,xc,noCell,gauss,gs)
!----------------------------------------------------------------------------
! This subroutine to set up Gauss points,Jacobian and weights for a cell
! input--ibk: the No. of the consider cell;
! numq: number of points for background cells;
! numcell: number of background cells;
! numgauss: number of Gauss points in a cell;
! k: number of Gauss points used, numgauss=k*k for 2D cell;
! xc(nx,numq): coordinates of points for background cells;
! noCell(ng,numcell): No. of points to form this cell;
! gauss(2,k): coefficients of Gauss points;
! nx,ng: parameters are defined in file parameter.h.
! output--gs(ng,numgauss): coordinate of the Gauss points, weight and Jacobian
!---------------------------------------------------------------------------
implicit real*8 (a-h,o-z)
include 'parameter.h'
dimension xc(nx,numq),noCell(ng,numcell),gauss(nx,k),gs(ng,numgauss)
dimension psiJ(ng),etaJ(ng),xe(ng),ye(ng),aN(ng),aNJpsi(ng),aNJeta(ng)
index=0
psiJ(1)=-1.
psiJ(2)=1.
psiJ(3)=1.
psiJ(4)=-1.
etaJ(1)=-1.
etaJ(2)=-1.
etaJ(3)=1.
etaJ(4)=1.
l=k
ie=ibk
do j=1,ng
je=noCell(j,ie)
xe(j)=xc(1,je)
ye(j)=xc(2,je)
enddo
do 10 i=1,l
do 10 j=1,l
index=index+1
eta=gauss(1,i)
psi=gauss(1,j)
do ik=1,ng
aN(ik)=.25*(1.+psi*psiJ(ik))*(1.+eta*etaJ(ik))
aNJpsi(ik)=.25*psiJ(ik)*(1.+eta*etaJ(ik))
aNJeta(ik)=.25*etaJ(ik)*(1.+psi*psiJ(ik))
enddo
xpsi=0.
ypsi=0.
xeta=0.
yeta=0.
do jk=1,ng
xpsi=xpsi+aNJpsi(jk)*xe(jk)
ypsi=ypsi+aNJpsi(jk)*ye(jk)
xeta=xeta+aNJeta(jk)*xe(jk)
yeta=yeta+aNJeta(jk)*ye(jk)
enddo
ajcob=xpsi*yeta-xeta*ypsi
xq=0.
yq=0.
do kk=1,ng
xq=xq+aN(kk)*xe(kk)
yq=yq+aN(kk)*ye(kk)
enddo
gs(1,index)=xq
gs(2,index)=yq
gs(3,index)=gauss(2,i)*gauss(2,j)
gs(4,index)=ajcob
10 continue
RETURN
END



SUBROUTINE SupportDomain(numnode,nx,gpos,x,ds,ndex,nv)
!----------------------------------------------------------------------------
! This subroutine to determines nodes in the support domain of a Gauss point
! input--numnode: total number of field nodes;
! nx=2: for 2D problem;
! x(nx,numnode): coordinates of all field nodes;
! numgauss: number of Gauss points in a cell;
! gpos(2): x and y coordinate of a Gauss point;
! ds(nx,numnode): sizes of support domain;
! input and output-- ndex: when input ndex=0;
! when return ndex is the number of nodes in the support domain
! output--nv(ndex): No. of field nodes in the support domain
!---------------------------------------------------------------------------
implicit real*8 (a-h,o-z)
dimension gpos(nx),x(nx,numnode),ds(nx,numnode),nv(numnode)
eps=1.e-16
ndex=0
do ik=1,numnode
nv(ik)=0
enddo
do ik=1,numnode
dx=ds(1,ik)-dabs(gpos(1)-x(1,ik))
dy=ds(2,ik)-dabs(gpos(2)-x(2,ik))
if((dx.ge.eps).and.(dy.ge.eps)) then
ndex=ndex+1
nv(ndex)=ik
end if
enddo
RETURN
END




SUBROUTINE PointStiffnessMatrix(ndex,weight,ajac,ph,Dmat,GSPk)
!----------------------------------------------------------------------------
! This subroutine to calculate sparse stiff matrix
! input--ndex: the number of nodes in the support domain;
! weight: weight of Gauss quadrature;
! ajac: Jacobian;
! dphix: first dirivetive of x of shape function;
! dphiy: first dirivetive of y of shape function;
! Dmat(3,3): the matrix of strain-stress;
! output--GSPk(2ndex,2ndex): sub-stiffness matrix of the Gauss point
!---------------------------------------------------------------------------
implicit real*8 (a-h,o-z)
dimension ph(10,ndex),Dmat(3,3),GSPk(2*ndex,2*ndex)
dimension bmat(3,2*ndex),dphix(ndex),dphiy(ndex)
nb=2*ndex
do i=1,ndex
dphix(i)=ph(2,i)
dphiy(i)=ph(3,i)
enddo
do ib=1,3
do jb=1,nb
Bmat(ib,jb)=0.
enddo
enddo
do in=1,ndex
j=2*in-1
k=2*in
Bmat(1,j)=dphix(in)
Bmat(1,k)=0.
Bmat(2,j)=0.
Bmat(2,k)=dphiy(in)
Bmat(3,j)=dphiy(in)
Bmat(3,k)=dphix(in)
enddo
do ii=1,nb
do jj=1,nb
GSPk(ii,jj)=0.
enddo
enddo
do ii=1,nb
do jj=1,nb
do kk=1,3
do mm=1,3
GSPk(ii,jj)=GSPk(ii,jj)+weight*ajac*Bmat(kk,ii)* &
Dmat(kk,mm)*Bmat(mm,jj)
enddo
enddo
enddo
enddo
RETURN
END




SUBROUTINE EssentialBC(numnode,pAlf,alfs,x,ds,ak,af,npEBCnum,npEBC,pEBC)
!----------------------------------------------------------------------------
! This subroutine to enforce point essential bc's using penalty method;
! input--numnode: total number of field nodes;
! pAlf: penalty Fac; npEBCnum: number of e. b.c points
! alfs: coefficient of support domain
! x(nx,numnode): coordinates of all field nodes;
! input and output-- ak[]: stiffness matrix;
! af{}:force vector.
!---------------------------------------------------------------------------
implicit real*8 (a-h,o-z)
include 'parameter.h'
COMMON/rpim/ALFC,DC,Q,nRBF
common/basis/mbasis
dimension npEBC(3,100),pEBC(2,100)
dimension x(nx,numnode),ds(2,numnode),ak(2*numd,2*(numnode)),af(2*numnode)
dimension nv(numnode),ph(10,numnode), x2(2)
maxak=0.
do iebc=1,2*numnode
if(abs(ak(iebc,iebc)).gt.maxak) maxak=abs(ak(iebc,iebc))
enddo
do 10 iebc=1,npEBCnum
ie=npEBC(1,iebc)
x2(1)=x(1,ie)
x2(2)=x(2,ie)
ndex=0
! call support(x2,x,ds,nv(1),numnode,nx,ndex)
call SupportDomain(numnode,nx,x2,x,ds,ndex,nv)
do ik=1,ndex
do jk=1,10
ph(jk,ik)=0.
enddo
enddo
call MLS_ShapeFunc_2D(gpos,x,nv,ds,phi,nx,numnode,ndex,mm)
do iee=1,ndex
ine=nv(iee)
do ii=1,ndex
jne=nv(ii)
if(npEBC(2,iebc).eq.1) then
ak((ine*2-1),(jne*2-1))=ak((ine*2-1),(jne*2-1))-pAlf*maxak* &
ph(1,iee)*ph(1,ii)
endif
if(npEBC(3,iebc).eq.1) then
ak((ine*2),(jne*2))=ak((ine*2),(jne*2))-pAlf*maxak* &
ph(1,iee)*ph(1,ii)
endif
enddo
if(npEBC(2,iebc).eq.1) then
uu=pEBC(1,iebc)
af(ine*2-1)=af(ine*2-1)-pAlf*uu*maxak*ph(1,iee)
endif
if(npEBC(3,iebc).eq.1) then
uu=pEBC(2,iebc)
af(ine*2)=af(ine*2)-pAlf*uu*maxak*ph(1,iee)
endif
enddo
10 continue
RETURN
END





SUBROUTINE NaturalBC_concentrated(x,nx,numnode,af,ds,alfs,npNBCnum,npNBC,pNBC)
implicit real*8 (a-h,o-z)
dimension npNBC(3,100),pNBC(2,100)
COMMON/rpim/ALFC,DC,Q,nRBF
common/basis/mbasis
dimension af(2*numnode),x(nx,numnode),ds(nx,numnode)
dimension ph(10,numnode),gpos(2),nv(numnode)
do 10 iebc=1,npNBCnum
ie=npNBC(1,iebc)
gpos(1)=x(1,ie)
gpos(2)=x(2,ie)
ndex=0
call SupportDomain(numnode,nx,gpos,x,ds,ndex,nv)
do kph=1,3*ndex
do ik=1,10
ph(ik,kph)=0.
enddo
enddo
call MLS_ShapeFunc_2D(gpos,x,nv,ds,phi,nx,numnode,ndex,mm)
do iee=1,ndex
ie=nv(iee)
uu=pNBC(1,iebc)
af(ie*2-1)=af(ie*2-1)+ph(1,iee)*uu
uu=pNBC(2,iebc)
af(ie*2)=af(ie*2)+ph(1,iee)*uu
enddo
10 continue
RETURN
END





SUBROUTINE naturalBC_distributed(numnode,numq,in,jn,alfs,x,xc,ds, &
gauss,nquado,force)
!----------------------------------------------------------------------------
! This subroutine to enforce point natural bc's;
! input-numnode, numq, in,jn,alfs,x,xc,ds,gauss, nquado.
! input and output-- force{}:force vector.
!---------------------------------------------------------------------------
implicit real*8 (a-h,o-z)
include 'parameter.h'
common/para/xlength,ylength,p,young,anu,aimo
COMMON/rpim/ALFC,DC,Q,nRBF
COMMON /basis/mbasis
dimension xc(nx,numq),gauss(2,nquado),force(2*numnode),x(nx,numnode)
dimension ph(10,numnode),gpos(2),nv(numnode),ds(nx,numnode)
ax=0.5*(xc(1,in)-xc(1,jn))
ay=0.5*(xc(2,in)-xc(2,jn))
bx=0.5*(xc(1,in)+xc(1,jn))
by=0.5*(xc(2,in)+xc(2,jn))
do il=1,nquado
gpos(1)=ax*gauss(1,il)+bx
gpos(2)=ay*gauss(1,il)+by
weight=gauss(2,il)
ajac=0.5*sqrt((xc(1,in)-xc(1,jn))**2+(xc(2,in)-xc(2,jn))**2)
aimo=(1./12.)*ylength**3
ty=(-1000./(2.*aimo))*(ylength*ylength/4.-gpos(2)*gpos(2))
call SupportDomain(numnode,nx,gpos,x,ds,ndex,nv)
do kph=1,ndex
do ik=1,10
ph(ik,kph)=0.
enddo
enddo
call MLS_ShapeFunc_2D(gpos,x,nv,ds,phi,nx,numnode,ndex,mm)
do ie=1,ndex
nn=nv(ie)
force(2*nn)=force(2*nn)+weight*ajac*ph(1,ie)*ty
enddo
enddo
END