program haowu integer N0,BCT(1000),N,i,ians double precision xb(1000),yb(1000),xm(1000),ym(1000), & nx(1000),ny(1000),lg(1000),BCV(1000),BBB, & phi(1000),dphi(1000),pint,dl,xi,eta,pi,tau,F print*,'Enter number of elements:' read*,N tau=0.25d0 c Putting N well spaced out points in anticlockwise direction on the circle pi=4d0*datan(1d0) dl=2d0*pi/dfloat(N) do 10 i=1,N xb(i)=dcos(dfloat(i-1)*dl) yb(i)=dsin(dfloat(i-1)*dl) 10 continue xb(N+1)=xb(1) yb(N+1)=yb(1) c Working out the collocation points, 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 do 30 i=1,N BCT(i)=1 BCV(i)=2d0*((1d0-tau)*F(xb(i),yb(i))+tau*F(xb(i+1),yb(i+1))) & -8d0*xm(i)*ym(i) BCV(N+i)=2d0*((1d0-tau)*F(xb(i+1),yb(i+1))+tau*F(xb(i),yb(i))) & -8d0*xm(N+i)*ym(N+i) 30 continue call DLELAP1(N,tau,xm,ym,xb,yb,nx,ny,lg,BCT,BCV,phi,dphi) c Post processing c Estimation of B BBB=0d0 do 40 i=1,N BBB=BBB+0.50*(F(xb(i),yb(i))+F(xb(i+1),yb(i+1))) & -0.5d0*(phi(i)+phi(i+N)) & +(phi(i+N)-phi(i)) & /dsqrt((xm(N+i)-xm(i))**2d0+(ym(N+i)-ym(i))**2d0) 40 continue BBB=BBB/dfloat(N) 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) dl=xi**2d0-eta**2d0 write(*,60)pint+BBB,dl,xi,eta 60 format('Numerical and exact values are:', & F14.6,' and',F14.6,' respectively at (',F14.6,',',F14.6,')') print*,'To continue with another point enter 1:' read*,ians if (ians.eq.1) goto 50 end function F(x,y) implicit none double precision F,x,y F=x**2d0-y**2d0+4d0*x*y return end subroutine DLELAP1(N,tau,xm,ym,xb,yb,nx,ny,lg, & BCT,BCV,phi,dphi) c BVP governed by 2D Laplace's equation c Determining unknowns on the boundary by discontinuous linear elements c N: number of elements c tau*L is the distance of the two selected points from endpoints of each element c xm,ym: real arrays containing the two selected points on each of the elements c xb,yb: real arrays containing starting points of elements c nx,ny: real arrays containing normal vectors to elements c lg: real arrays containing lengths of elements c BCT: integer arrays containing type of boundary conditions c BCT(k)=0 if phi is specified on k-th element c BCV: real arrays containing values specified on the two points on each element c phi: real arrays containing values of phi on the two points on each element (returned) c dphi: real arrays containing values of normal derivative of phi on the two points of each element (returned) integer m,k,N,BCT(1000) double precision xm(1000),ym(1000),xb(1000),yb(1000), & nx(1000),ny(1000),lg(1000),BCV(1000),A(1000,1000), & B(1000),pi,PF1,PF2,delkm,delnkm,phi(1000),dphi(1000),F1,F2, & Z(1000),PF3,PF4,F3,F4,tau pi=4d0*datan(1d0) do 10 m=1,2*N B(m)=0d0 do 5 k=1,N call DPF(xm(m),ym(m),xb(k),yb(k),nx(k),ny(k),lg(k), & PF1,PF2,PF3,PF4) F1=PF1/pi F2=PF2/pi F3=PF3/pi F4=PF4/pi if (k.eq.m) then delkm=1d0 else delkm=0d0 endif if (k.eq.(m-N)) then delnkm=1d0 else delnkm=0d0 endif if (BCT(k).eq.0) then A(m,k)=((1d0-tau)*F1-F3/lg(k))/(2d0*tau-1d0) A(m,N+k)=(-tau*F1+F3/lg(k))/(2d0*tau-1d0) B(m)=B(m)-((-(1d0-tau)*F2+F4/lg(k)) & /(2d0*tau-1d0)-0.5d0*delkm)*BCV(k) & -((tau*F2-F4/lg(k)) & /(2d0*tau-1d0)-0.5d0*delnkm)*BCV(N+k) else A(m,k)=(-(1d0-tau)*F2+F4/lg(k))/(2d0*tau-1d0)-0.5d0*delkm A(m,N+k)=(tau*F2-F4/lg(k))/(2d0*tau-1d0)-0.5d0*delnkm B(m)=B(m)-(((1d0-tau)*F1-F3/lg(k))*BCV(k) & +(-tau*F1+F3/lg(k))*BCV(N+k))/(2d0*tau-1d0) endif 5 continue 10 continue call solver(A,B,2*N,1,Z) do 15 m=1,N if (BCT(m).eq.0) then phi(m)=BCV(m) phi(N+m)=BCV(N+m) dphi(m)=Z(m) dphi(N+m)=Z(N+m) else phi(m)=Z(m) phi(N+m)=Z(N+m) dphi(m)=BCV(m) dphi(N+m)=BCV(N+m) endif 15 continue return end subroutine DLELAP2(N,tau,xi,eta,xb,yb,nx,ny,lg,phi,dphi,pint) c BVP governed by 2D Laplace's equation c Discontinuous linear element calculation of solution at interior point (xi,eta). c N: number of elements c tau*L: distance of two selected points from the endpoints of the boundary c xb,yb: real arrays containing points on the boundary c nx,ny: real arrays containing normal vectors to elements c lg: real arrays containing the lengths of elements c phi: real arrays containing values of phi on two selected points of each element c dphi: real arrays containing values of normal derivative of phi on two selected points of each element c pint: real value of phi at the point (xi,eta) integer N,i double precision xi,eta,xb(1000),yb(1000),nx(1000),ny(1000), & lg(1000),phi(1000),dphi(1000),pint,sum,pi,PF1,PF2,PF3,PF4, & tau pi=4d0*datan(1d0) sum=0d0 do 10 i=1,N call DPF(xi,eta,xb(i),yb(i),nx(i),ny(i),lg(i), & PF1,PF2,PF3,PF4) sum=sum+(phi(i)*(-(1d0-tau)*lg(i)*PF2+PF4) & +phi(N+i)*(tau*lg(i)*PF2-PF4) & -dphi(i)*(-(1d0-tau)*lg(i)*PF1+PF3) & -dphi(N+i)*(tau*lg(i)*PF1-PF3))/lg(i) 10 continue pint=sum/(pi*(2d0*tau-1d0)) return end subroutine DPF(xi,eta,xk,yk,nkx,nky,L,PF1,PF2,PF3,PF4) double precision xi,eta,xk,yk,nkx,nky,L,PF1,PF2, & PF3,PF4,A,B,E,D,BA,EA A=L**2d0 B=2d0*L*(-nky*(xk-xi)+nkx*(yk-eta)) E=(xk-xi)**2d0+(yk-eta)**2d0 D=dsqrt(dabs(4d0*A*E-B**2d0)) BA=B/A EA=E/A if (D.lt.0.0000000001d0) then PF1=0.5d0*L*(dlog(L) & +(1d0+0.5d0*BA)*dlog(dabs(1d0+0.5d0*BA)) & -0.5d0*BA*dlog(dabs(0.5d0*BA))-1d0) PF2=0d0 PF4=0d0 else PF1=0.25d0*L*(2d0*(dlog(L)-1d0)-0.5d0*BA*dlog(dabs(EA)) & +(1d0+0.5d0*BA)*dlog(dabs(1d0+BA+EA)) & +(D/A)*(datan((2d0*A+B)/D)-datan(B/D))) PF2=L*(nkx*(xk-xi)+nky*(yk-eta))/D & *(datan((2d0*A+B)/D)-datan(B/D)) PF4=-0.5d0*BA*L*PF2 & +(0.25d0*L*L/A)*(nkx*(xk-xi)+nky*(yk-eta)) & *(dlog(dabs(A+B+E))-dlog(dabs(E))) endif PF3=-0.5d0*BA*L*PF1+(L*L/(8d0*A)) & *((A+B+E)*(dlog(dabs(A+B+E))-1d0)-E*(dlog(dabs(E))-1d0)) return end subroutine SOLVER(A,B,N,lud,Z) c Solution of AX=B (N unknowns) by LU decomposition and then backward substitutions. c Solution is returned in X with A and B unmodified. c If an LU decomposition is needed on A, set lud to any integer value except 0 c when calling the subroutine. Otherwise, if the decomposition is not needed, c give lud the value 0 (zero) on calling. c Adapted from the subroutines DGEFA and DGESL in LINPACK originally written c by Cleve Moler, University of New Mexico, Argonne National Lab. c Supporting subroutines DAXPY and DSCAL and function IDAMAX c are also available in LINPACK. integer lda,N,ipvt(1000),info,lud,IDAMAX, & j,k,kp1,l,nm1,kb double precision A(1000,1000),B(1000),Z(1000),t,AMD(1000,1000) common /ludcmp/ipvt,AMD nm1=N-1 do 5 i=1,N Z(i)=B(i) 5 continue if (lud.eq.0) goto 99 do 6 i=1,N do 6 j=1,N AMD(i,j)=A(i,j) 6 continue c This part for decomposing A is taken from DGEFA. info=0 if (nm1.lt.1) go to 70 do 60 k=1,nm1 kp1=k+1 l=IDAMAX(N-k+1,AMD(k,k),1)+k-1 ipvt(k)=l if (AMD(l,k).eq.0.0d0) goto 40 if (l.eq.k) goto 10 t=AMD(l,k) AMD(l,k)=AMD(k,k) AMD(k,k)=t 10 continue t=-1.0d0/AMD(k,k) call DSCAL(N-k,t,AMD(k+1,k),1) do 30 j=kp1,N t=AMD(l,j) if (l.eq.k) go to 20 AMD(l,j)=AMD(k,j) AMD(k,j)=t 20 continue call DAXPY(N-k,t,AMD(k+1,k),1,AMD(k+1,j),1) 30 continue go to 50 40 continue info=k 50 continue 60 continue 70 continue ipvt(N)=N if (AMD(N,N).eq.0.0d0) info=N if (info.ne.0) pause 'Division by zero in SOLVER!' 99 continue c This part for finding solution of AX=B is taken from DGESL. if (nm1.lt.1) goto 130 do 120 k=1,nm1 l=ipvt(k) t=Z(l) if (l.eq.k) goto 110 Z(l)=Z(k) Z(k)=t 110 continue call DAXPY(N-k,t,AMD(k+1,k),1,Z(k+1),1) 120 continue 130 continue do 140 kb=1,N k=N+1-kb Z(k) = Z(k)/AMD(k,k) t=-Z(k) call DAXPY(k-1,t,AMD(1,k),1,Z(1),1) 140 continue return end subroutine DAXPY(N,da,dx,incx,dy,incy) c LINPACK subroutine originally written by Jack Dongarra double precision dx(1000),dy(1000),da integer i,incx,incy,ix,iy,m,mp1,N if(N.le.0) return if (da .eq. 0.0d0) return if(incx.eq.1.and.incy.eq.1) goto 20 ix=1 iy=1 if(incx.lt.0) ix=(-N+1)*incx+1 if(incy.lt.0) iy=(-N+1)*incy+1 do 10 i=1,N dy(iy)=dy(iy)+da*dx(ix) ix=ix+incx iy=iy+incy 10 continue return 20 m=mod(N,4) if( m.eq. 0 ) go to 40 do 30 i=1,m dy(i)=dy(i)+da*dx(i) 30 continue if(N.lt.4) return 40 mp1=m+1 do 50 i=mp1,N,4 dy(i)=dy(i)+da*dx(i) dy(i+1)=dy(i+1)+da*dx(i+1) dy(i+2)=dy(i+2)+da*dx(i+2) dy(i+3)=dy(i+3)+da*dx(i+3) 50 continue return end subroutine DSCAL(N,da,dx,incx) c LINPACK subroutine originally written by Jack Dongarra double precision da,dx(1000) integer i,incx,m,mp1,N,nincx if(N.le.0.or.incx.le.0) return if(incx.eq.1) goto 20 nincx = N*incx do 10 i=1,nincx,incx dx(i)=da*dx(i) 10 continue return 20 m=mod(N,5) if(m.eq.0) goto 40 do 30 i=1,m dx(i) = da*dx(i) 30 continue if(N.lt.5) return 40 mp1 = m + 1 do 50 i=mp1,N,5 dx(i)=da*dx(i) dx(i+1)=da*dx(i+1) dx(i+2)=da*dx(i+2) dx(i+3)=da*dx(i+3) dx(i+4)=da*dx(i+4) 50 continue return end function IDAMAX(N,dx,incx) c LINPACK function originally written by Jack Dongarra double precision dx(1000),dmax integer i,incx,ix,N,IDAMAX IDAMAX = 0 if(N.lt.1.or.incx.le.0 ) return IDAMAX = 1 if(N.eq.1)return if(incx.eq.1) goto 20 ix = 1 dmax = dabs(dx(1)) ix = ix + incx do 10 i=2,N if(dabs(dx(ix)).le.dmax) goto 5 IDAMAX=i dmax=dabs(dx(ix)) 5 ix=ix+incx 10 continue return 20 dmax=dabs(dx(1)) do 30 i=2,N if(dabs(dx(i)).le.dmax) goto 30 IDAMAX=i dmax=dabs(dx(i)) 30 continue return end c Copyright 2006-2007 WT Ang, A Beginner's Course in Boundary Element Methods