efg2d.f95

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

SUBROUTINE SolverBand(ak,fp,neq,nmat)
!------------------------------------------------------------------
! Sloving linear equations; it calls BandSolver & GaussSolver
! input-ak,fp,neq,nmat
! output--fp
!------------------------------------------------------------------
implicit real*8 (a-h,o-z)
dimension ak(nmat,nEq),fp(nmat)
real*8, allocatable :: tp(:,:)
real*8,allocatable :: stfp(:,:)
allocate (tp(1:neq,1:nmat))
allocate (stfp(1:neq,1:neq))
ep=1.d-10
do i=1,nEq
do j=1,nEq
stfp(i,j)=0.
tp(i,j)=0.
enddo
enddo
do i=1,nEq
do j=1,nEq
stfp(i,j)=ak(i,j)
enddo
enddo
ni=nEq
Lp=0 ! half band width
do 20 i=1,ni
do j=ni,i,-1
if(stfp(i,j).ne.0.) then ! stfp[,] stifness matrix
if(abs(j-i).gt.Lp) Lp=abs(j-i)
go to 21
endif
enddo
21 continue
do j=1,i
if(stfp(i,j).ne.0.) then
if(abs(j-i).gt.Lp) Lp=abs(j-i)
go to 20
endif
enddo
20 continue
ilp=2*lp+1 ! band width
nm=nEq
if(ilp.lt.nEq) then
call BandSolver(stfp,fp,tp,nm,lp,ilp,nmat) ! solver for band matrix
else
call GaussSolver(nEq,nmat,ak,fp,ep,kkkk) ! standard solver
endif
deallocate (tp)
deallocate (stfp)
END
SUBROUTINE BANDSOLVER(A,F,B,N,L,IL,nmat)
!------------------------------------------------------------------
! Slover for banded linear equations
!------------------------------------------------------------------
implicit real*8 (a-h,o-z)
DIMENSION A(N,N),F(N)
DIMENSION B(N,nmat),d(n,1)
M=1
LP1=L+1
DO I=1,N
DO K=1,IL
B(I,K)=0.
IF(I.LE.LP1) B(I,K)=A(I,K)
IF(I.GT.LP1.AND.I.LE.(N-L)) B(I,K)=A(I,I+K-LP1)
IF(I.GT.(N-L).AND.(I+K-LP1).LE.N) B(I,K)=A(I,I+K-LP1)
ENDDO
ENDDO
DO I=1,N
D(I,1)=F(I)
ENDDO
IT=1
IF (IL.NE.2*L+1) THEN
IT=-1
WRITE(*,*)'***FAIL***'
RETURN
END IF
LS=L+1
DO 100 K=1,N-1
P=0.0
DO I=K,LS
IF (ABS(B(I,1)).GT.P) THEN
P=ABS(B(I,1))
IS=I
END IF
ENDDO
IF (P+1.0.EQ.1.0) THEN
IT=0
WRITE(*,*)'***FAIL***'
RETURN
END IF
DO J=1,M
T=D(K,J)
D(K,J)=D(IS,J)
D(IS,J)=T
ENDDO
DO J=1,IL
T=B(K,J)
B(K,J)=B(IS,J)
B(IS,J)=T
ENDDO
DO J=1,M
D(K,J)=D(K,J)/B(K,1)
ENDDO
DO J=2,IL
B(K,J)=B(K,J)/B(K,1)
ENDDO
DO I=K+1,LS
T=B(I,1)
DO J=1,M
D(I,J)=D(I,J)-T*D(K,J)
ENDDO
DO J=2,IL
B(I,J-1)=B(I,J)-T*B(K,J)
ENDDO
B(I,IL)=0.0
ENDDO
IF (LS.NE.N) LS=LS+1
100 CONTINUE
IF (ABS(B(N,1))+1.0.EQ.1.0) THEN
IT=0
WRITE(*,*)'***FAIL***'
RETURN
END IF
DO J=1,M
D(N,J)=D(N,J)/B(N,1)
ENDDO
JS=2
DO 150 I=N-1,1,-1
DO K=1,M
DO J=2,JS
D(I,K)=D(I,K)-B(I,J)*D(I+J-1,K)
ENDDO
ENDDO
IF (JS.NE.IL) JS=JS+1
150 CONTINUE
if(it.le.0) write(*,*) "BandSolver failed"
DO I=1,N
F(I)=D(I,1)
ENDDO
RETURN
END
SUBROUTINE GaussSolver(n,mk,a,b,ep,kwji)
!------------------------------------------------------------------
! Stnadard Gauss elimination slover for linear equations that are
! not suitably solved by BandSolver.
!------------------------------------------------------------------
implicit real*8 (a-h,o-z)
dimension a(mk,mk),b(mk)
integer, allocatable :: m(:)
allocate (m(2*mk))
ep=1.0e-10
kwji=0
do i=1,n
m(i)=i
enddo
do 20 k=1,n
p=0.0
do 30 i=k,n
do 30 j=k,n
if(abs(a(i,j)).le.abs(p)) goto 30
p=a(i,j)
io=i
jo=j
30 continue
if(abs(p)-ep) 200,200,300
200 kwji=1
return
300 if(jo.eq.k) goto 400
do i=1,n
t=a(i,jo)
a(i,jo)=a(i,k)
a(i,k)=t
enddo
j=m(k)
m(k)=m(jo)
m(jo)=j
400 if(io.eq.k) goto 500
do j=k,n
t=a(io,j)
a(io,j)=a(k,j)
a(k,j)=t
enddo
t=b(io)
b(io)=b(k)
b(k)=t
500 p=1.0/p
in=n-1
if(k.eq.n) goto 600
do j=k,in
a(k,j+1)=a(k,j+1)*p
enddo
600 b(k)=b(k)*p
if(k.eq.n) goto 20
do i=k,in
do j=k,in
a(i+1,j+1)=a(i+1,j+1)-a(i+1,k)*a(k,j+1)
enddo
b(i+1)=b(I+1)-a(i+1,k)*b(k)
enddo
20 continue
do i1=2,n
i=n+1-i1
do j=i,in
b(i)=b(i)-a(i,j+1)*b(j+1)
enddo
enddo
do k=1,n
i=m(k)
a(1,i)=b(k)
enddo
do k=1,n
b(k)=a(1,k)
enddo
kwji=0
deallocate (m)
return
END



SUBROUTINE GetDisplacement(x,ds,u2,disp,alfs,nx,numnode)
!----------------------------------------------------------------------------
! This subroutine to get the final displacements from
! displacement parameters using the MFree interpolation;
! input--numnode: total number of field nodes;
! alfs: coefficent of support support
! x(nx,numnode): coordinates of all field nodes;
! u2(2*numnode): displacement parameters;
! input and output-- disp: final displacements.
!---------------------------------------------------------------------------
implicit real*8 (a-h,o-z)
COMMON/rpim/ALFC,DC,Q,nRBF
common/basis/mbasis
dimension x(nx,numnode),ds(nx,numnode),gpos(nx),u2(nx,numnode)
dimension disp(2*numnode)
dimension ph(10,numnode), nv(numnode)
write(2,*)'Displacements of field nodes'
nn=2*numnode
do i=1,nn
disp(i)=0.
enddo
ind=0
do 50 id=1,numnode
ind=ind+1
gpos(1)= x(1,id)
gpos(2)=x(2,id)
ndex=0
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)
nc1=2*ind-1
nc2=2*ind
do kk=1,ndex
m=nv(kk)
disp(nc1)=disp(nc1)+ph(1,kk)*u2(1,m)
disp(nc2)=disp(nc2)+ph(1,kk)*u2(2,m)
enddo
50 continue
do ii=1,numnode
write(2,52) ii,disp(2*ii-1),disp(2*ii)
enddo
52 format(1x,i5,2e20.5)
RETURN
END




SUBROUTINE GetStress(x,noCell,ds,Dmat,u2,alfs,nx,numnode,numgauss,&
xc,gauss,nquado,ng,numq,numcell, ENORM,Stressnode)
!----------------------------------------------------------------------------
! This subroutine to get the stress and energy error;
! input--numnode: total number of field nodes;
! numcell: number of cells;
! numq: total number of points for cells;
! alfs: coefficent of support support;
! x(nx,numnode): coordinates of all field nodes;
! xc(nx,numcell): coordinates of all points for cells;
! u2(2*numnode): displacement parameters;
! ds(nx,numnode): sizes of influence domain;
! Dmat(3,3): material matrix;
! nquado: number of Gauss points in a cell;
! gauss(nx,nquado): coefficients of Gauss points;
! numgauss: total number of Gauss points in all cells;
! output-- Enorm: energy error;
! Stressnode:stress for field nodes;
! compute out--Stress: stress for Gauss points;
! stressex, strne: exact stresse for beam problem.
!---------------------------------------------------------------------------
implicit real*8 (a-h,o-z)
common/para/xlength,ylength,p,young,anu,aimo
COMMON/rpim/ALFC,DC,Q,nRBF
common/basis/mbasis
dimension noCell(4,numcell),ds(nx,numnode),x(nx,numnode),u(2*numnode)
dimension xc(nx,numnode),gauss(nx,nquado)
dimension Dmat(3,3),u2(nx,numnode)
dimension Stressnode(3,numnode),strne(3,numnode)
dimension stress(3),stressex(3),err(3),Dinv(3,3),der(3)
integer, allocatable :: nv(:)
integer, allocatable :: ne(:)
real*8, allocatable :: ph(:,:)
real*8,allocatable :: gpos(:)
real*8, allocatable :: gs(:,:)
real*8,allocatable :: bmat(:,:)
allocate ( nv(1:numnode) )
allocate ( ne(1:2*numnode) )
allocate ( ph(1:10,1:3*numnode) )
allocate ( gpos(1:nx) )
allocate ( gs(1:ng,1:numgauss) )
allocate ( bmat(1:3,1:2*numnode) )
close(37)
open(37, file='midstr.dat')
do id=1,3
do jd=1,3
Dinv(id,jd)=Dmat(id,jd)
enddo
enddo
invd=3
ep=1.d-10
call GetINVASY(invd,invd,Dinv,EP)
do iu=1,numnode
ju=2*iu-1
ku=2*iu
u(ju)=u2(1,iu)
u(ku)=u2(2,iu)
enddo
enorm=0.
!****************Compute energy error
do 100 ibk=1,numcell
ind=0
call CellGaussPoints(ibk,numcell,nquado,numq,numgauss,&
xc,noCell,gauss,gs)
do 200 is=1,numgauss
do i=1,3
stress(i)=0.
stressex(i)=0.
enddo
ind=ind+1
gpos(1)= gs(1,is)
gpos(2)=gs(2,is)
weight=gs(3,is)
ajac=gs(4,is)
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)
nb=2*ndex
do in=1,nb
ne(in)=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
do ib=1,3
do jb=1,nb
Bmat(ib,jb)=0.
enddo
enddo
do inn=1,ndex
j=2*inn-1
k=2*inn
Bmat(1,j)=ph(2,inn)
Bmat(1,k)=0.
Bmat(2,j)=0.
Bmat(2,k)=ph(3,inn)
Bmat(3,j)=ph(3,inn)
Bmat(3,k)=ph(2,inn)
enddo
do ii=1,3
do kk=1,3
do mm=1,nb
mn=ne(mm)
stress(ii)=stress(ii)+&
Dmat(ii,kk)*Bmat(kk,mm)*u(mn)
enddo
enddo
enddo
!****************Exact stress for beam problem
stressex(1)=(1./aimo)*p*(xlength-gpos(1))*gpos(2)
stressex(2)=0.
stressex(3)=-0.5*(p/aimo)*(0.25*ylength*ylength-gpos(2)*gpos(2))
do ier=1,3
err(ier)=stress(ier)-stressex(ier)
enddo
do jer=1,3
der(jer)=0.
do ker=1,3
der(jer)=der(jer)+Dinv(jer,ker)*err(ker)
enddo
enddo
err2=0.
do mer=1,3
err2=err2+weight*ajac*(0.5*der(mer)*err(mer))
enddo
enorm=enorm+err2
200 continue
100 continue
!****************Compute nodal stresses
write(2,*)'Stress of field nodes'
do is=1,numnode
gpos(1)= x(1,is)
gpos(2)=x(2,is)
do ii=1,3
Stressnode(ii,is)=0.
enddo
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)
nb=2*ndex
do in=1,nb
ne(in)=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
do ib=1,3
do jb=1,nb
Bmat(ib,jb)=0.
enddo
enddo
do inn=1,ndex
j=2*inn-1
k=2*inn
Bmat(1,j)=ph(2,inn)
Bmat(1,k)=0.
Bmat(2,j)=0.
Bmat(2,k)=ph(3,inn)
Bmat(3,j)=ph(3,inn)
Bmat(3,k)=ph(2,inn)
enddo
do ii=1,3
do kk=1,3
do mm=1,nb
mn=ne(mm)
Stressnode(ii,is)=Stressnode(ii,is)+&
Dmat(ii,kk)*Bmat(kk,mm)*u(mn)
enddo
enddo
enddo
strne(1,is)=(1./aimo)*p*(xlength-gpos(1))*gpos(2)
strne(2,is)=0.
strne(3,is)=-0.5*(p/aimo)*(0.25*ylength*ylength-gpos(2)*gpos(2))
write(2,220) is,Stressnode(1,is),Stressnode(2,is),Stressnode(3,is)
if(abs(gpos(1)-24).le.1.d-5) then
write(37,240) is,gpos(2),Stressnode(1,is),Stressnode(2,is), &
Stressnode(3,is),strne(1,is),strne(2,is),strne(3,is)
endif
enddo
enorm=dsqrt(enorm)
write(2,230) enorm
220 format(1x,i5,3e20.5)
230 format(1x,'Energy Error=',e20.10)
240 format(1x,i5,f8.3,6e15.5)
deallocate ( nv)
deallocate ( ne)
deallocate ( ph)
deallocate ( gpos)
deallocate ( gs )
deallocate ( bmat)
RETURN
END




SUBROUTINE GetInvasy(N,MA,A,EPS)
!----------------------------------------------------------------------------
! This subroutine to get INVARSION OF A(N,N) USING THE GAUSS-JODON METHOD.
! MATRIX A MUST BE DEFINITE BUT MAY BE ASYMMETRIC.
! input--N: dimension of A;
! MA: max number of rows of A;
! EPS: tolerance;
! Input and output--A[N,N]: the matrix in the input and invarsion in output;
!---------------------------------------------------------------------------
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION A(MA,N)
DO 10 K=1,N
C=A(K,K)
IF(DABS(C).LE.EPS)pause
C=1.0/C
A(K,K)=1.0
DO J=1,N
A(K,J)=A(K,J)*C
ENDDO
DO 10 I=1,N
IF(I.EQ.K)GOTO 10
C=A(I,K)
A(I,K)=0.0
DO J=1,N
A(I,J)=A(I,J)-A(K,J)*C
ENDDO
10 CONTINUE
RETURN
END