program CVBEM c CVBEM for Example 1.2 in "A Beginner's Course in Boundary Element Methods". c Revised 14 Sept 2010: Function AG integer N0,BCT(1000),N,i,ians,p double precision xb(1000),yb(1000),xm(1000),ym(1000), & nx(1000),ny(1000),lg(1000),BCV(1000),ov2p,pi,qq,rr, & dl,xi,eta,AA(1000,1000),BB(1000),gam,thet,THETA, & UV(1000) double complex zb(1000),zm(1000),nxy(1000),FZ,zint pi=4d0*datan(1d0) ov2p=0.5d0/pi print*,'Enter integer N0 (<42):' read*,N0 N=12*N0 c Putting N well spaced out points in anticlockwise direction on the quarter annular boundary do 10 i=1,8*N0 dl=pi/dfloat(16*N0) xb(i+N0)=2d0*dcos(dfloat(i-1)*dl) yb(i+N0)=2d0*dsin(dfloat(i-1)*dl) if (i.le.N0) then dl=1d0/dfloat(N0) xb(i)=1d0+dfloat(i-1)*dl yb(i)=0d0 xb(i+9*N0)=0d0 yb(i+9*N0)=2d0-dfloat(i-1)*dl endif if (i.le.(2*N0)) then dl=pi/dfloat(4*N0) xb(i+10*N0)=dsin(dfloat(i-1)*dl) yb(i+10*N0)=dcos(dfloat(i-1)*dl) endif 10 continue xb(N+1)=xb(1) yb(N+1)=yb(1) zb(N+1)=dcmplx(xb(N+1),yb(N+1)) c Working out the midpoints, lengths and normal vectors of the elements do 20 i=1,N zb(i)=dcmplx(xb(i),yb(i)) xm(i)=0.5d0*(xb(i)+xb(i+1)) ym(i)=0.5d0*(yb(i)+yb(i+1)) zm(i)=dcmplx(xm(i),ym(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) nxy(i)=dcmplx(nx(i),ny(i)) 20 continue c Specifying the boundary conditions do 30 i=1,N if ((i.le.N0).or.((i.gt.(9*N0)).and.(i.le.(10*N0)))) then BCT(i)=1 BCV(i)=0d0 else if ((i.gt.N0).and.(i.le.(9*N0))) then BCT(i)=0 BCV(i)=3d0*dcos(4d0*datan(ym(i)/xm(i))) else BCT(i)=0 BCV(i)=dcos(4d0*datan(ym(i)/xm(i))) endif 30 continue c The CVBEM part: setting up system of linear algebraic equations do 35 p=1,2*N-1 do 35 i=1,2*N-1 AA(p,i)=0d0 35 continue c The first N-1 equations (from discretisation of the first Cauchy integral formula) c First N unknowns are the constants denoted by u and the next N-1 unknowns are v. do 50 p=1,N-1 BB(p)=0d0 do 40 i=1,N gam=dlog(cdabs((zb(i+1)-zm(p))/(zb(i)-zm(p)))) thet=THETA(zb(i),zb(i+1),zm(p)) AA(p,i)=-ov2p*gam if (i.ne.N) then AA(p,i+N)=ov2p*thet if (p.eq.i) then AA(p,i+N)=AA(p,i+N)-1d0 endif endif 40 continue 50 continue c The next N equations (from the boundary conditions) do 70 p=1,N BB(N-1+p)=BCV(p) if (BCT(p).eq.0) then AA(N-1+p,p)=1d0 else do 60 i=1,N call QIR(zb(i),zb(i+1),zm(p),nxy(p),qq,rr) AA(N-1+p,i)=qq if (i.ne.N) AA(N-1+p,N+i)=-rr 60 continue endif 70 continue c Solving the system of linear algebraic equations call solver(AA,BB,2*N-1,1,UV) UV(2*N)=0d0 c Post processing 150 print*,'Enter the coordinates xi and eta of interior point' read*,xi,eta zint=dcmplx(xi,eta) FZ=dcmplx(0d0,0d0) do 155 p=1,N gam=dlog(cdabs((zb(p+1)-zint)/(zb(p)-zint))) thet=THETA(zb(p),zb(p+1),zint) FZ=FZ+(ov2p/dcmplx(0d0,1d0)) & *dcmplx(gam,thet)*dcmplx(UV(p),UV(N+p)) 155 continue dl=(xi**2d0+eta**2d0)**2d0 dl=((16d0/85d0)*(dl-1d0/dl) & -(16d0/255d0)*(dl/16d0-16d0/dl)) & *dcos(4d0*datan(eta/xi)) write(*,160) dreal(FZ),dl,xi,eta 160 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 150 end function THETA(z1,z2,zm) implicit none double precision THETA,AG,pi,a1,a2,dif double complex z1,z2,zm pi=4d0*datan(1d0) a1=AG(dreal(z1-zm),dimag(z1-zm)) a2=AG(dreal(z2-zm),dimag(z2-zm)) dif=a2-a1 if (dabs(dif-pi).lt.0.000001d0) then THETA=pi goto 231 endif if (dabs(dif+pi).lt.0.0000001d0) then THETA=pi goto 231 endif if (dabs(dif).le.pi) then THETA=dif else if (dif.gt.(pi)) then THETA=dif-2d0*pi else if (dif.lt.(-pi)) then THETA=dif+2d0*pi endif 231 continue return end function AG(x,y) implicit none double precision x,y,AG,pi pi=4d0*datan(1d0) if (x.gt.0d0) then AG=datan(y/x) else if (y.gt.0d0) then AG=pi-datan(dabs(y/x)) else if (y.lt.0d0) then AG=datan(dabs(y/x))-pi else AG=pi endif return end subroutine QIR(z1,z2,zm,nz,q,r) implicit none double precision q,r,pi double complex z1,z2,zm,nz,qr pi=4d0*datan(1d0) qr=(nz/(pi*dcmplx(0d0,1d0))) & *(-1d0/(z2-zm)+1d0/(z1-zm)) q=dreal(qr) r=dimag(qr) 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 2010 WT Ang