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

END IF
IF (id == 0) CYCLE
Direction_Q: DO jj=1,Ndof
Node_points: DO n=1,Nodel
nD= Ndof*(n-1) + jj ! column number in array
dUe(iD,nD)= dUe(iD,nD) + Ni(n)*UP(ii,jj)*Jac*Weit
dTe(iD,nD)= dTe(iD,nD) + Ni(n)*TP(ii,jj)*Jac*Weit
END DO Node_points
END DO Direction_Q
END DO Direction_P
END DO Gauss_points_eta
END DO Gauss_points_xsi
Eta1= Eta2
END DO Subdivisions_eta
Xsi1= Xsi2
END DO Subdivisions_xsi
END DO Colloc_points
!---------------------------------------------------------
! Part 1 : Pi is one of the element nodes
!---------------------------------------------------------
Colloc_points1: DO i=1,Ncol
lnod= 0
DO n= 1,Nodel ! Determine which local node is Pi
IF (Inci(n) .EQ. i) THEN
lnod=n
END IF
END DO
IF (lnod .EQ. 0) CYCLE ! None -> next Pi
Ntri= 2
IF (lnod > 4) Ntri=3 ! Number of triangles
Triangles: DO ntr=1,Ntri
CALL Tri_RL(RLx,RLe,Elengx,Elenge,lnod,ntr)
Mi= Ngaus(RLx,2,Rlim)
IF(Mi == 5) Mi=4 ! Triangles are not sub-divided
Call Gauss_coor(Glcorx,Wix,Mi)
Ki= Ngaus(RLe,2,Rlim)
IF(Ki == 5) Ki=4
Call Gauss_coor(Glcore,Wie,Ki)
Gauss_points_xsi1: DO m=1,Mi
xsib= Glcorx(m)
Gauss_points_eta1: DO k=1,Ki
etab= Glcore(k)
Weit= Wix(m)*Wie(k)
CALL Trans_Tri(ntr,lnod,xsib,etab,xsi,eta,Jacb)
CALL Serendip_func(Ni,xsi,eta,ldim,nodel,Inci)
Call Normal_Jac(Vnorm,Jac,xsi,eta,ldim,nodel,Inci,elcor)
Jac= Jac*Jacb
CALL Cartesian(GCcor,Ni,ldim,elcor)
r= Dist(GCcor,xPi(:,i),Cdim)
dxr= (GCcor-xPi(:,i))/r
IF (Ndof .EQ. 1) THEN