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

of Gauss points. For simplicity, both are integrated using the number of Gauss points
required for the higher order singularity. Indeed, the subroutine presented has not been
programmed very efficiently but, for the purpose of this topic, simplicity was the
paramount factor. Additional improvements in efficiency can, for example, be made by
carefully examining if the operations in the DO loops actually depend on the DO loop
variable. If they do not, then that operation should be taken outside of the corresponding
DO loop. Substantial savings can be made here for a program that involves up to seven
implied DO loops and which has to be executed for all boundary elements.
The second part of the SUBROUTINE deals with the case where P i is one of the
nodes of the element which we integrate over. To deal with the singularity of the
integrand the element has to be subdivided into 2 or 3 triangles, as explained previously.
Since there are a lot of implied DO loops involved, we show a structure chart of this part
of the program in Figure 6.18.
A subdivision of the integration region has been implemented, but in order to improve
clarity of the structure chart is not shown there. The subdivision of integration involves
two more DO loops.
SUBROUTINE
Integ3(Elcor,Inci,Nodel,Ncol,xPi,Ndof,E,ny,ko,dUe,dTe,Ndest,Isym)
!--------------------------------------------------
! Computes [dT]e and [dU]e for 3-D problems
!-------------------------------------------------
IMPLICIT NONE
REAL, INTENT(IN) :: Elcor(:,:) ! Element coordinates
INTEGER, INTENT(IN) :: Ndest(:,:) ! Node destination vector
INTEGER, INTENT(IN) :: Inci(:) ! Element Incidences
INTEGER, INTENT(IN) :: Nodel ! No. of Element Nodes
INTEGER , INTENT(IN):: Ncol ! Number of points Pi
REAL , INTENT(IN) : : xPi(:,:) ! coll. points coords.
INTEGER , INTENT(IN):: Ndof ! Number DoF /node (1 or 3)
INTEGER , INTENT(IN):: Isym
REAL , INTENT(IN) :: E,ny ! Elastic constants
REAL , INTENT(IN) :: ko
REAL(KIND=8) , INTENT(OUT) :: dUe(:,:),dTe(:,:)
REAL :: Elengx,Elenge,Rmin,RLx,RLe,Glcorx(8),Wix(8),Glcore(8)&
,Wie(8),Weit,r
REAL :: Ni(Nodel),Vnorm(3),GCcor(3),dxr(3),Jac,Jacb,xsi,eta,xsib&
,etab,Rlim(2)
REAL :: Xsi1,Xsi2,Eta1,Eta2,RJacB,RonL
REAL :: UP(Ndof,Ndof),TP(Ndof,Ndof) ! for storing kernels
INTEGER :: i,m,n,k,ii,jj,ntr,Mi,Ki,id,nd,lnod,Ntri,NDIVX&
,NDIVSX,NDIVE,NDIVSE,MAXDIVS
INTEGER :: ldim= 2 ! Element dimension
INTEGER :: Cdim= 3 ! Cartesian dimension
ELengx=&
Dist((Elcor(:,3)+Elcor(:,2))/2.,(Elcor(:,4)+Elcor(:,1))/2.,Cdim)
ELenge=&
Dist((Elcor(:,2)+Elcor(:,1))/2.,(Elcor(:,3)+Elcor(:,4))/2.,Cdim)
dUe= 0.0 ; dTe= 0.0 ! Clear arrays for summation