Dealing with corners and changing geometry - The Boundary Element Method with Programming - page 333

! -> not singular -> normal gauss quadrature
!-------------------------------------------------
Element_nodes:DO n=1,Nodel
i=Inci(n) ! Collocation node
Shape_function: DO ns=1, Nodel
! Node at element -> shape function in integral term
IF(n /= ns) CYCLE
Region_Loop: DO nreg=1, 2
Mi= 8
Call Gauss_coor(Glcor,Wi,Mi)
Gauss_points: DO m=1,Mi
SELECT CASE (n)
CASE (1) ! Node1
IF (nreg==1) THEN ! right
xsi= (1.-d2)/2. + Glcor(m) * (1.+d2)/2.
dxdxb= (1.+d2)/2.
ELSE IF (nreg==2) THEN ! left
xsi= (-1.-d2)/2. + Glcor(m) * (1.-d2)/2
dxdxb= (1.-d2)/2
END IF
CASE (2) ! Node2
IF (nreg==1) THEN ! right
xsi= (1.+d1)/2. + Glcor(m)*(1.-d1)/2
dxdxb= (1.-d1)/2
ELSE IF (nreg==2) THEN ! left
xsi= (-1.+d1)/2. + Glcor(m)*(1.+d1)/2
dxdxb= (1.+d1)/2
END IF
CASE (3) ! Node3
IF (nreg==1) THEN ! right
xsi= 0.5 + Glcor(m) * 0.5
dxdxb= 0.5
ELSE IF (nreg==2) THEN ! left
xsi= -0.5 + Glcor(m) * 0.5
dxdxb= 0.5
END IF
CASE DEFAULT
END SELECT
CALL Normal_Jac(Vnorm,Jac,xsi,eta,ldim,nodel,Inci,elcor)
CALL Serendip_func(Ni,xsi,eta,ldim,nodel,Inci)
CALL Cartesian(GCcor,Ni,ldim,elcor)
r= Dist(GCcor,xP(:,i),cdim) ! Dist. P,Q
dxr= (GCcor-xP(:,i))/r ! rx/r , ry/r
CALL ShapefunctionDisc(Ni,xsi,eta,ldim,nodel,Inci)
Direction1:DO j=1,2
IF (Isym == 0) THEN
iD= 2*(i-1) + j
ELSE
iD= Ndest(i,j)
END IF
IF (id == 0) CYCLE