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

Colloc_points: DO i=1,Ncol
Rmin= Min_dist1 (Elcor,xP(:,i),Nodel,inci,ELeng,Eleng,ldim)
RonL= Rmin/Eleng ! R/L
!-----------------------------------------------------
! Integration off-diagonal coeff. -> normal Gauss Quadrature
!---------------------------------------------------------
Mi= Ngaus (RonL,1,Rlim) ! Number of Gauss points for o(1/r)
NDIVS= 1
RJacB=1.0
IF (Mi == 5) THEN ! Determine number of subdiv. In [
IF (RonL > 0.0) NDIVS= INT(RLim(2)/RonL) + 1
IF (NDIVS > MAXDIVS) MAXDIVS= NDIVS
RJacB= 1.0/NDIVS
Mi=4
END IF
Call Gauss_coor(Glcor,Wi,Mi) ! Assign coords/Weights
Xsi1=-1
Subdivisions: DO NDIV=1,NDIVS
Xsi2= Xsi1 + 2.0/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= UK(dxr,r,E,ny,Cdim) ; TP= TK(dxr,r,Vnorm,ny,Cdim)
Node_points: DO n=1,Nodel
Direction_P: DO j=1,2
IF (Isym == 0)THEN
iD= 2*(i-1) + j
ELSE
iD= Ndest(i,j) ! line number in array
END IF
IF (id == 0) CYCLE
Direction_Q: DO k= 1,2
nD= 2*(n-1) + k ! column number in array
IF (Dist(Elcor(:,n),xP(:,i),cdim) > epsi) THEN
dUe(iD,nD)= dUe(iD,nD) + Ni(n)*UP(j,k)*Jac*Wi(m)*RJacB
dTe(iD,nD)= dTe(iD,nD) + Ni(n)*TP(j,k)*Jac*Wi(m)*RJacB
ELSE
dUe(iD,nD)= dUe(iD,nD) + &
Ni(n)*C*dxr(j)*dxr(k)*Jac*Wi(m)*RJacB
END IF
END DO Direction_Q
END DO Direction_P
END DO Node_points
END DO Gauss_points
Xsi1= Xsi2
END DO Subdivisions