Boundary Element Methods–Numerical Implementation - The Boundary Element Method with Programming - page 145

Call Gauss_coor(Glcor,Wi,Mi) ! Assign coords/Weights
Xsi1=-1
Subregions: &
DO NDIV=1,NDIVS
IF (NDIVS > 1) THEN
Xsi2= Xsi1+2/NDIVS
Gauss_points:&
DO m=1,Mi
xsi= Glcor(m)
IF (NDIVS > 1) Xsi= 0.5*(Xsi1+Xsi2)+xsi/NDIVS
CALL Serendip_func(Ni,xsi,eta,ldim,Nodel,Inci)
Call Normal_Jac(Vnorm,Jac,xsi,eta,ldim,Nodel,Inci,elcor)
CALL Cartesian(GCcor,Ni,ldim,elcor) ! Coords of Gauss pt
r= Dist(GCcor,xP(:,i),cdim) ! Dist. P,Q
dxr= (GCcor-xP(:,i))/r ! rx/r , ry/r
UP= U(r,k,cdim) ; TP= T(r,dxr,Vnorm,cdim) ! Kernels
Node_points: &
DO n=1,Nodel
IF (Dist(Elcor(:,n),xP(:,i),cdim) < epsi) EXIT ! P i is n
dUe(i,n)= dUe(i,n) + Ni(n)*UP*Jac*Wi(m)*RJACB
dTe(i,n)= dTe(i,n) + Ni(n)*TP*Jac*Wi(m)*RJACB
END DO &
Node_points
END DO &
Gauss_points
END DO &
Subregions
END DO &
Colloc_points
!------------------------------
! Diagonal terms of dUe
!------------------------------
c1= 1/(2.0*pi*k)
Colloc_points1: &
DO i=1,Ncol
Node_points1: &
DO n=1,Nodel
IF (Dist(Elcor(:,n),xP(:,i),cdim) > Epsi) CYCLE ! Pi not n
Nreg=1
IF (n == 3) nreg= 2
!-------------------------------------------
! Integration of logarithmic term
!--------------------------------------------
Subregions: &
DO nr=1,Nreg
Mi= 4
Call Gauss_Laguerre_coor(Glcor,Wi,Mi)
Gauss_points1: &
DO m=1,Mi
SELECT CASE (n)
CASE (1)