program EXPRD2 C Example for 2D potential problem with periodic conditions solved with C discontinuous linear elements integer N1,N2,BCT(1000),N,i,ians double precision xb(1000),yb(1000),xm(1000),ym(1000), & nx(1000),ny(1000),lg(1000),BCV(1000), & phi(1000),dphi(1000),pint,dl,xi,eta,pi,tau N1=10 N2=30 N=2*N1+2*N2 tau=.22d0 C Putting N well spaced out points in anticlockwise direction on the square boundary pi=4d0*datan(1d0) dl=0.5d0/dfloat(N1) do 10 i=1,N1 xb(i)=0d0 yb(i)=0.50d0-dfloat(i-1)*dl xb(N1+N2+i)=2d0*pi yb(N1+N2+i)=dfloat(i-1)*dl 10 continue dl=2d0*pi/dfloat(N2) do 12 i=1,N2 xb(N1+i)=dfloat(i-1)*dl yb(N1+i)=0d0 xb(N1+N2+N1+i)=2d0*pi-dfloat(i-1)*dl yb(N1+N2+N1+i)=0.5d0 12 continue xb(N+1)=xb(1) yb(N+1)=yb(1) C Working out the midpoints, lengths and normal vectors of the elements do 20 i=1,N xm(i)=xb(i)+tau*(xb(i+1)-xb(i)) ym(i)=yb(i)+tau*(yb(i+1)-yb(i)) xm(N+i)=xb(i)+(1d0-tau)*(xb(i+1)-xb(i)) ym(N+i)=yb(i)+(1d0-tau)*(yb(i+1)-yb(i)) lg(i)=dsqrt((xb(i+1)-xb(i))**2d0+(yb(i+1)-yb(i))**2d0) nx(i)=(yb(i+1)-yb(i))/lg(i) ny(i)=(xb(i)-xb(i+1))/lg(i) 20 continue C Specifying the boundary conditions on the sides of the square do 30 i=1,N if ((i.gt.N1).and.(i.le.(N1+N2))) then BCT(i)=0 BCV(i)=dsin(xm(i)) BCV(N+i)=dsin(xm(N+i)) else if (i.gt.(2*N1+N2).and.i.le.N) then BCT(i)=1 BCV(i)=dexp(0.5d0)*dsin(xm(i)) BCV(N+i)=dexp(0.5d0)*dsin(xm(N+i)) else BCT(i)=2 endif 30 continue call DLELAPPER(N,N1,N2,tau,xm,ym,xb,yb,nx,ny,lg,BCT,BCV,phi,dphi) 50 print*,'Enter coordinates xi and eta of an interior point:' read*,xi,eta call DLELAP2(N,tau,xi,eta,xb,yb,nx,ny,lg,phi,dphi,pint) write(*,60)pint,dexp(eta)*dsin(xi) 60 format('Numerical and exact values are:', & F14.6,' and',F14.6,' respectively') print*,'To continue with another point enter 1:' read*,ians if (ians.eq.1) goto 50 end