Multiple regions - The Boundary Element Method with Programming - page 300

REAL(KIND=8), INTENT(OUT) :: KBE(:,:) ! Stiffness matrix
REAL(KIND=8), INTENT(OUT) :: A(:,:) ! u due to u c =1
REAL(KIND=8), INTENT(OUT) :: TC( :) ! t due to u c =0
INTEGER, ALLOCATABLE :: Ldeste(:,:)! Element destinations
REAL(KIND=8), ALLOCATABLE :: dUe(:,:),dTe(:,:),Diag(:,:)
REAL(KIND=8), ALLOCATABLE :: Lhs(:,:)
REAL(KIND=8), ALLOCATABLE :: Rhs(:),RhsM(:,:) ! right hand sides
REAL(KIND=8), ALLOCATABLE :: u1(:),u2(:,:) ! results
REAL, ALLOCATABLE :: Elcor(:,:)
REAL :: v3(3),v1(3),v2(3)
INTEGER :: Nodes,Dof,k,l,nel
INTEGER :: n,m,Ndofs,Pos,i,j,nd
Nodes= UBOUND (xP,2) ! total number of nodes of region
Ndofs= Nodes*Ndof ! Total degrees of freedom of region
ALLOCATE (Ldeste(maxe,Ndofe)) ! Elem. destination vector
!------------------------------------------
! Determine Element destination vector
!---------------------------------------------
Elements:&
DO Nel=1,Maxe
k=0
DO n=1,Nodel
DO m=1,Ndof
k=k+1
IF (Ndof > 1) THEN
Ldeste(Nel,k)= ((Incie(Nel,n)-1)*Ndof + m)
ELSE
Ldeste(Nel,k)= Incie(Nel,n)
END IF
END DO
END DO
END DO &
Elements
ALLOCATE (dTe(Ndofs,Ndofe),dUe(Ndofs,Ndofe))
ALLOCATE (Diag(Ndofs,Ndof))
ALLOCATE (Lhs(Ndofs,Ndofs),Rhs(Ndofs),RhsM(Ndofs,Ndofs))
ALLOCATE (u1(Ndofs),u2(Ndofs,Ndofs))
ALLOCATE (Elcor(Cdim,Nodel))
!---------------------------------------------------------------
! Compute and assemble element coefficient matrices
!---------------------------------------------------------------
Lhs= 0
Rhs= 0
Elements_1:&
DO Nel=1,Maxe
Elcor(:,:)= xP(:,Incie(nel,:))! gather element coords
IF (Cdim == 2) THEN
IF (Ndof == 1) THEN
CALL Integ2P(Elcor,Incie(nel,:),Nodel,Nodes,xP,C1,dUe,dTe)
ELSE
CALL Integ2E(Elcor,Incie(nel,:),Nodel,Nodes&