Postprocessing - The Boundary Element Method with Programming - page 248

Element_loop1: &
DO NEL= 1,MAXE
Symmetry_loop1:&
DO nsy=1,Nsym
Inci= Incie(nel,:)
Elcor= xp(:,Inci(:))
IF (ldim == 2) THEN
ELengx= Dist((Elcor(:,3)+Elcor(:,2))/2.&
,(Elcor(:,4)+Elcor(:,1))/2.,Cdim) ! Lxsi
ELenge= Dist((Elcor(:,2)+Elcor(:,1))/2.&
,(Elcor(:,3)+Elcor(:,4))/2.,Cdim) ! Leta
ELSE
Call Elength(Elengx,Elcor,nodel,ldim)
END IF
Ldest= 1
Fac= 1.0
Fac_nod=1.0
El_ue(:,:)=El_u(Nel,:,:)
El_te(:,:)=El_t(Nel,:,:)
IF (Isym > 0) THEN
DO Nod=1,Nodel
dofa= (nod-1)*Ndof+1
dofe= dofa+Ndof-1
El_trac(dofa:dofe)= El_te(Nod,:)
EL_disp(dofa:dofe)= El_ue(Nod,:)
END DO
CALL Mirror(Isym,nsy,Nodes,Elcor,Fac,Inci&
,Ldest,El_trac,EL_disp &
,nodel,ndof,Cdim)
DO Nod=1,Nodel
dofa= (nod-1)*Ndof+1
dofe= dofa+Ndof-1
El_te(Nod,:)= El_trac(dofa:dofe)
El_ue(Nod,:)= El_disp(dofa:dofe)
Fac_nod(Nod,:)= Fac(dofa:dofe)
END DO
END IF
Rmin= Min_dist(Elcor,xPnt,Nodel,ldim,Inci)
Mi= Ngaus(Rmin/Elengx,Cdim-1,Rlim)
NDIVSX= 1 ; NDIVSE= 1
RJacB=1.0
RonL= Rmin/Elengx
IF (Mi == 5) THEN ! subdivision required
IF (RonL > 0.0) NDIVSX= INT(RLim(2)/RonL) + 1
IF (NDIVSX > MAXDIVS) MAXDIVS= NDIVSX
Mi=4
END IF
CALL Gauss_coor(Glcorx,Wix,Mi) ! Coords/Wghts x dir
Ki= 1 ; Wie(1)= 1.0 ; Glcore(1)= 0.0
IF (Cdim == 3) THEN
Ki= Ngaus(Rmin/Elenge,Cdim-1,Rlim) !