An efficient solution algorithm for boundary element equations JOSI~ L. ORTIZ and C. V. G I R I J A V A L L A B H A N Texas Tech UniverstO', Lubbock, Texas 79409, USA
Boundary element techniques result in the solution of a linear system of equations of the type HU = GQ + B, which can be transformed into a system of equations of the type AXF. The coefficient matrix A requires the storage of a full matrix on the computer. This storage requirement, of the order of n*n memory positions (n = number of equations), for a very large n is often considered negative for the boundary element method. Here, two algorithms are presented where the memory requirements to solve the system are only n*(n1)/2 and n'n/4 respectively. The algorithms do not necessitate any external storage devices nor do they increase the computational efforts. Key Words: Boundary Element Method, Two efficient solution algorithms, Less computer storage, Computer code. INTRODUCTION The boundary element method, like other d~screte numerical methods such as finite element and finite &fferences methods, results m a set of linear simultaneous equauons In other numerical methods, the global system of equauons ts represented by one coefficient matrix, which ~s symmetric, often has banded characteristics when properly numbered, and is stored in a long half banded array. On the other hand, in the boundary element method, there are two full coefficient matrices, H and G. Even though the matrix H is symmetric, matrix G Js unsymmetnc, and since they are full matrices, the use of partiuoning and storing on auxiliary devices can be extremely complex and difficult to program. Earlier computer codes 1 ~ used 2*n*n storage positions for H and G matrices, where n is the total number of unknowns m the model. This excessive storage requirement for very large systems has been pointed out negatively while dehberating the capabilmes of the boundary element method. Computer codes 2"3 ~ have been made more efficient, where only a single full matrix A is employed. However. for very large systems, the storage requirement is relatwely high. The authors have developed two soluuon algorithms that reqmre n ' I n  1)/2 and n'n/4 storage locations, respectwely, for handhng the coefficient
Paper accepted Februar3 1991 Discussioncloses Januar3 1992
matrices. The first algorithm is similar to a work done by M. C. Au 5. These algorithms are developed in such a way that they neither necessitate any external storage devices nor increase the overall computational effort. It is also shown that the algorithms do not decrease the accuracy of the solution vector.
ALGORITHM I The computer program developed here is similar to a work done by M. C. Au 5 and uses the Gauss Method. A set of subroutines (see Listing 1 in the Appendix) is presented here to illustrate the capability of this algorithm. The equation that has to be solved is of the form HU=GQ+B
(1)
where U and Q are vectors representing the field variable and its normal derivatives on the boundary. B is a vector resulting from the domain integral and H and G are the coefficient matrices. Later, this system is converted to the form AX = F
(2)
where X is the vector of all unknown quantities. This is done with the use of the known prescribed boundary conditions and a K O D E vector which indicates, for each equation, if the field or its derivative is prescribed on the elements on the boundary. It should be observed that the ith rows of the H and G matnces in equation (1) are formed due to a collocation of the fundamental solution at some point i; and that the resulting coefficients are completely independent of the rest of the rows. This characteristic makes it possible to form any row at any time. Furthermore, it has been shown 2"3 that there is no need to form the entire H and G matrices before transforming them into the A matrix; it is possible to transform any row of H and G into the corresponding row of A. Explmting this idea, the algorithm, in one main loop, forms the ith row of the A matrix (see GaussALG1 subroutine m the Appendix). Instead of forming and storing the ~th row in an array A, the row is formed and stored m a working array V of dimension equal to n. If z= 1, the array V is normalized making V(I)= 1 (diagonal
C 1991 ElsevierScience Pubhshers Ltd
Adt. Eng Sqftware, 1991. Vol. 13, No, 4
197
,In effictent solutton algortthm for boundary element equattons
3rd eqn.
.
.
.
01 001 000 000
.
1 01
0 0 11
0 0 0
Ith eqn v, vz
ELEMENTSNOT YETCALCULATED
J L Orttz and C I" Gu tla l'al/ahhwz
voiv ,.
. v, I
I Vo
ALL ZE~ROS ~ 1
ELEMENTSNOTYET CALCULATED a. "r'~o ROWS TRANSFORMED 3 ro ROW READY FOR GAUSS ELIMINATION
b. (I1)ROW$ TRANSFORMED Ith ROW READY FOR GAUSS ELIMINATION
C. MAXIMUM STORAGE
REQUIREMENTFORA( )
Fig. 1. Algortthm I
coefficient) a n d o n l y the next n  1 t e r m s are s t o r e d m a r r a y A. If i = 2, the a r r a y V ns t r a n s f o r m e d by G a u s s e l i m i n a t m n (see T r a n s f R o w l s u b r o u t i n e m the A p p e n d i x ) . a n d n o r m a h z e d by m a k i n g V ( 2 ) = 1, and finally o n l y the next n  2 t e r m s are sent to a r r a y A for storage. W h e n d e a l i n g w i t h the ith row, the a l g o r t t h m t r a n s f o r m s the array, n o r m a l i z e s by m a k i n g V(0 = 1. a n d o n l y ~ t coefficients are sent to a r r a y A, after V has been t r a n s f o r m e d (see S t o r e R o w l s u b r o u t i n e nn the A p p e n d i x ) T h e c o r r e s p o n d i n g e l e m e n t of the force v e c t o r is also t r a n s f o r m e d a c c o r d i n g l y wnthm this loop. F i g u r e 1 deptcts h o w this m a r e l o o p w o r k s a n d s h o w s t h a t there ~s n o n e e d to h a v e the rest of the e q u a t t o n s of m a t r i x A m m e m o r y to factortze a p a r t i c u l a r row. Yh~s results in a r e q u i r e d s t o r a g e of o n l y n * ( n  1)/2 m e m o r y p o s m o n s for the a r r a y A. T a b l e 1 s h o w s the s t o r a g e r e q u i r e m e n t s a n d T a b l e 2 s h o w s the t o t a l n u m b e r of operatmns.
Table 1 Storage requtrement Jor A matrzx alter the nth row transJormatzon Storage reqmrement for A matrix Row Algorithm I 1
n
Algornthm II
I
l*(n
1)
2 3
2*n  3 3*n  6
2*(n  2) 3*(n  31
t
t*n  t*(z + 1)/2
t*(n  I)
n'In
(n
n
1
n
1)/2

11"1
0
Tahle 2 ,Vumhero/operatums m Algortthm 1 Gau~ Method bl Baek~ubstLtuuon
al Factonzatlon and normahzatnon Factorlzatlon
Normahzatmn Row
Row Row of matrnx
Load vector
1
m=~=0
re=s=0
d=Inl)
d=p
2 3 4 5
m = ~ = ( n  I) m=~=lnI)+tn2t m=~=(nI)+fn2l+ln
m=s=t*p rn=~=2*p m=~=3*p
d=(n21 d=tn3)
d=p d=p
m=~=In2)*p m=~=(nI)*p
d=l d=0
nI n
m=s=(nI)+(n2)+ m=s=(nI)+(n2)+
Sum
m = ~ = ~ t"
. J
3)
+2+2 +3+2+1
~
In
m=s=
I
=",V ~ t' +(n  I)*n*p] ,
N()TF
n = number of equa/nons p = number o[ load vet.It)r,,
i ", =
~,um
d =
dlvl~,l(Hl:
oz =
198
~, ' . , u b t r L l c t l o n ' . , ",
mu.ltlpln~.,ttkms
Adt' Enq SoJtware, 1991, Vol 13, No 4
Row matrnx
I)*n*p


(n
d=
d=p d=p l)*n

Load ,rector
Vector
d=n*p
n
HZ =
~, ~ 0
I'll =
~, =
t*p
n  2
m
', =
~*p
n  3
m = s =
1l

I
=
3*p
2
m = ~ = ( n  2)*p
I
m=s=lnl)*p
Sum
In D1 =
1 }*n*p
An eO~clent .~olutton algorlthnl ~or boundary element equation.s" J. L. Ortiz and C V. Gtrila Vallahhan Tahh'
3 ,~Utllhet o/ opet a11o11~ II1 A I,k,ot 111111111 Gau~Jordem M(,thod Factorlzatlon
Normalization
Updating
upperIdentity matrix
Ro~ Rov~ o f m a t r i x
Load vector
Row matrix
l
re=s=0
m=
d=(n
l)
2 3 4 5
i i 1 = s = 1 " 1 1 1  11 m = s = 2 " 1 1 7  21 m = s = 3 " 0 1  3)
ii1 = s = I * p m = s = 2*p nl = s = 3 * p
d =(nd =In
2) 3)
hi=s=111_) m=s=(11li*l
117=s=(n2)*p
d=l
m=s=(nI)*p
d=O ( n  1)*n d =  
I11 i1

~=0
.1 Sum
m=~=
111 l)*n*p t*(nr~
~
nl=~
i=1
[
Total sum
NOTE
n
d
=
~
( n  I)* n "1
number
p = number
+ n*p
1
+
[
m
Rows of matrix
Load vector
d=p
m = s =0
m = s =0
d = p d = p
m = ~ = I * ( n  21 m = s = 2 * I n  3) m = s = 3"(n4)
m = s = l*p m = s = 2*p m = s = 3*p
d=p d=p
m=s=ln2)*l m=,,=O
m=s=ln2)*p m=s=(nI)*p
d=n*p
re=s=
.2
"
n1
n2
i=1
i=l
~
( n  I)*n*p t*(nll)
m=s
1=1
= s= ~" t'Int)+ ~ i * ( n   i  I ) + l n   l l * n * p

1
= sums/subtractions d = divisions m = multlphcanons
of equanons of load vectors
The second main loop is the calculation of the solution vector which is accomplished by the standard backsubstitution. Usual checking of the zero on the diagonal can also be accomplished without any difficulty. The storage requirement for this algorithm, compared to the usual storage required in textbook codes is approximately 50 %, without any additional computational effort.
algorithm I. In the same way, the row Js stored in a working array V which is transformed (see TransfRow2 subroutine in the Appendix) and normalized. On the other hand, before sending the transformed row to the array A, a GaussJordan elimination is performed in the rows already stores in A (see UpdateIdentity subroutine in the Appendix). Figure 21a) depicts the situation before row 4 has been transformed. Observe that, at this moment, the required storage for A is only 3 * ( n  3) positions. After the fourth equation is transformed and normalized, a top left identity matrix 4*4 is updated. Note that at this time only ( n  4) coefficients are to be sent to A and only 4*(n  4) coefficients are needed to be stored (see StoreRow2 subroutine). The force vector is also transformed respectively inside this main loop. Notice that there is no need to store the rows below the row being transformed and, as shown in Table 1, the storage requirement increases from the beginning to the middle of the matrix and then it gets reduced. The maximum storage requirement for the coefficients related
A L G O R I T H M II Algorithm II is very similar to the Algorithm I. Here, along with the Gauss elimination, a GaussJordan transformation is applied to the A matrix as its coefficients are developed row after row. (See Listing 2 in the Appendix.) Table 1 shows the storage requirements and Table 3 shows the total number of operations. The algorithm, in this case, consists of one main loop where it forms the ith row of matrix A as in the case of
II
0 o I STORED AS ONE ~0 1 0 DIMENSIONAL 0 0 1 ARRAY 4th eqn V~ i/,2 V~ V4VIi . . . .
Vector
JV n
ELEMENTS NOT YET CALCULATED
11
10 01
0t0 0i0
000 000
0 0 0 0 0 0 00
1 0 0 1 0 0 00
0 0 0 STORED 0 0 0 AS ONE 1 0 IO DIMENSIONAL 01:0 ARRAY
00
00
00:1
ith eqn. V I V~, V~ . . . . . . . . .
//////+
//L////~.,'+
0
Vi
ALL ZEROS
ELEMENTS NOT YET CALCULATED
/ / / / / / +(n.a) / / / / / / 41n4) ////// I 1 I I
1
(n.1)1 1
b 01) ROWS TRANSFORMED ith ROW READY FOR GAUSS JORDAN ELIMINATION
FZ~ 2
b 01)ROWS TRANSFORMED ith ROW READY FOR GAUSS JORDAN ELIMINATION
c. MAXIMUM STORAGE REQUIREMENT FOR A ( . ) =r~/4 If n Is even ,,(r~l)/41f n Is odd
A/gorlthm II
Adt Eng Software, 1991, Vol 13, No 4
199
An effictent solutton algortthm jor boundary element equattons J L. Ortlz and C ~'. Gwtla ~'allahhan to the matrix A is the maximum of z * ( n  0 , which ~s n'n~4 if n is even and (n*n  I)/4 if n is odd. When n is large, the maximum storage requirement for the matrxx A is slightly over 25 % of the one used'in conventional boundary element programs 2"3. In algorithm II, the transformed force vector automatlcally becomes the unknown vector X, m other words. there ,s no backsubstitution as required m Algorithm I. A minor drawback in Algorithm II is that the storage m A must be rearranged before an equation ~s added. However, the code presented here uses a check ff a new equation can be added before a rearrangement. If not, the entire array is rearranged. As shown below, the total number of matrix operations in both Algorithms I and II is the same. From Table 3 we notice that the total number of operations for the GaussJordan method involves the following summations: .1
n2
l*(ni)+ ~ i=l
CONCLUSIONS It is apparent that both of the algorithms presented here are efficient and do not require any additional computational effort as compared to techniques published earlier Algorithm II is much more efficient than the Algorithm I While algorithm I can be generalized to use external storage when deahng with a ~ery large set of equat,ons. the second algorithm should be the worst choice to program when usmg external storage
REFERENCES 1 Brebbla. C A The Boundar~ Element '~lethod/or Enqlneets Pentech Press, London. 1978 2 Banerjee, P K and Butterfield. R Bounda~ Element ~lerhod~ m Eno,neerm O Sczence. McGrawH~ll Book Co. New York. 1981 3 Brebbla. C A. Telles, J C F and Wrobet L C Boundar~ Element Teehmques. SprmgerVerlag. Berhn. 1984 4 Brebbta. C A and Dommguez, J Boundary Elements. ~m Introductory Course. Computational Mechamcs Pubhcatlons and McGrawHall Book Co. New York. 1989 5 Au. M C Boundary element dLrect ehmmatton procedure m multiregion problems Proceedings o! the Boundar~ Element oth Internattonal Conference, 1984
nI
i*[nz1)=
1=1
tz
~ l=l
where the equality can easdy be demonstrated. Thus, by comparing Tables 2 and 3 we see that the total number of operat,ons in the solution process ~s the same.
APPENDIX The following listings present general purpose basic codes for boundary element programs using Algorithms I and II. The user should update the main program, add parameters m the subroutines calls and write additional subroutines depending on the type of problem involved. Three points ( ... ) are used to indicate the additional parameters which must be supplied by the reader. The code is written in QuickBASIC language.
Ltstmg 1
' Main
program
n X[ n
using
Algorithm
I
= number of equations = force vector
]
A[ V[
n*(n1)/2 n ] ip[ n ]
tol
] = = =
array to array to auxiliary
= tolerance
store the A matrix store a row of A matrix vector to swap rows
used
to
check
singularity
T h e f o l l o w i n g s u b r o u t i n e s h o u l d be p r e p a r e d d e p e n d i n g on t h e t y p e of p r o b l e m d e a l i n g w i t h : sub FormRow( i,V(1),xi, ... ) s t a t i c sub InputData( ... ) s t a t i c sub OutputResults( ... ) s t a t i c sub InteriorPoints(...) ) static defdbl ah,oz :defint in dim X(100),A(4950),V(100),ip(100) common
const
shared
n, t o l
ok=l,notOk=0
' Tolerance
tol
200
= 0.000001
Adv Enq. Sojtware, 1991, Vol 13, No 4
' Global
variables
An e~ctent ,solutton algortthm /or boundal:v element equations J L. Orttz and C. V Gt,'t/a Vallahhan ' Read
data call
' Form
InputData(
...
)
and solve system of equations call GaussALGl(X(),V(),ip(),A(),
...
' Compute potential in i n t e r i o r p o i n t s • call InteriorPoints( ... ) call OutputResults( ... )
)
print
results
end sub
StoreRowl(k,V(1),A(1),xi,X(1)
' Stores
transformed
kth
row
)
static
of A matrix
' D o n o t s t o r e t h e l a s t e q u a t i o n b e c a u s e it if k < n t h e n ii = ( n  l ) * ( k  l )  (kl) * ( k  2 ) / 2 + 1 for i = k+l to n A(ii) = V(i) : ii = ii + 1 next i e n d if ' Store
the
X(k)
=
end
sub
sub
coefficient
of
the
force
is
a one.
vector
xi
TransfRowl(
• Transforms
kth
k•V(1),A(1),xi,X(1),status row
of A matrix.
Gauss
) static
elimination
• R o w I is n o t t r a n s f o r m e d if k > l t h e n ii = 1 for i = 1 to k1 c
=
v(i)
for
end
j=i+l to n V ( j ) = V ( j )  c * A (ii) next j xi = x i  c * X (i) next i if
• Check pivot less than zero c = V (k) if a b s ( c ) < t o l t h e n status = notOk : exit else status = Ok e n d if • Normalize row making if k < n t h e n for i=k+l to n V(i) =V(i)/c next i end if
: ii=ii+l
sub
diagonal
coefficient
I.
Adt En9 Software, 1991, Vol 13, No 4
201
An ephctent solutton algortthm /or boundary element equations J L Ortt and C I . G,'tla I allahhan
' Normalize xi=xi/c end
corresponding
force
vector
value
sub
sub GaussALGI(X(1),V(1),ip(1),A(1), ' Solves
A
* U = X.
' ip() u s e d t o s w a p for i=l to n ip(i) = i next i
U is rows
returned when
...
) static
in a r r a y
pivot
is
close
X. to
zero
' M a i n a l g o r i t h m . S t a t u s = O k if p i v o t n o t c l o s e to z e r o f o r k=l t o n j =k+ 1 do call FormRow ( ip(k),V(),xi, ... ) call TransfRowl ( k,V(),A(),xi,X(),status ) if s t a t u s = O k t h e n call StoreRowl (k,V(),A(),xi,X() ) else if j = n + 1 then goto ErrorHandling else jk=ip(k) : ip(k)=ip(j) : ip(j)=jk : j=j+l e n d if e n d if loop while status = notOk next k ' Back substitution ii = ( n  1 ) * n / 2 for b=l to n1 k=n~ : kl=k+l f o r j=n t o kl s t e p i X(k)=X(k)A(ii)*X(j) : ii = ii next j next exit
 1
sub
' E r r o r h a n d l i n g in c a s e o f s i n g u l a r i t y ErrorHandling: print "Singularity in r o w = " ; k ; '' Program end end
202
sub
Adv Enq So/tware, 1991, Vol 13, No 4
aborted"
An e~hctent solution algorithm/'or houndarv element equatton.s J L. Ortiz and C. V. Gtrqa Vallahhan Ltsting 2 ' Main
program
' n • X[ n ] • A[ n'n/4 ' ip[n] ' V[ n ] • tol • iwb • ihb •nm
using
= = = = = = = = =
]
Algorithm
II
number of equations force vector array to store A matrix. auxiliary vector to swap array to store a row of tolerance used to check width of block storing height of block storing maximum size of array A
• The following subroutine ' depending on the type of
rows the A singularity A matrix A matrix
should be prepared problem dealing with:
f
'
sub sub sub sub
• •
•
FormRow( i,V(1),xi, ... ) s t a t i c InputData( ... ) s t a t i c OutputResults( ... ) s t a t i c InteriorPoints(...) ) static
#
defdbl ah,oz :defint in dim X(100) ,A(2500), ip(100) ,V(100) c o m m o n s h a r e d n, t o l , i w b , i h b , n m const
data call
• Form
variables
ok=l, notOk=0
' Tolerance and maximum to1 = 0.000001 nm = 2500 • Read
• Global
size
InputData(
...
of
array
A
)
and solve system of equations call GaussJorALG2(X(),V(),ip(),A(),
• Compute potential in interior points• call InteriorPoints( ... ) call OutputResults( ... )
... print
)
results
end
sub
StoreRow2(k,Vll),A(1),xi,X(I)
' Stores
the
kth
row
of
the
)
static
A matrix
• Rearrange storage if no space for new equation if k>ihb and k
)
if
no
space.
Program
aborted"
Adv En9 Software, 1991, Vol 13, No 4
203
An effictent solutton algortthm /or boundary element equatton~
f o r i  2 t o k  i ii = ( i  l ) * i w w + l f o r j=l t o n  k A(jj) = A(ii) next j next i e n d if
J L Ortlz amt C I Glrtla I'alh~hhan
: jj
=
(il)*iwb+l
: ii
=
ii
• Store the new equation except if k < n t h e n ii = ( k  1 ) * i w b + l f o r i = n t o k + l s t e p I A(ii)=V(i) : ii = i i + l next i
+ 1
last
: jj
=
jj
+
1
one
end if ' Stores
the
X(k)
=
end
sub
value
of the
force
vector
xi
• ~ Z ~ ~  k ~ B ~ B ~ D I R U S m l l l l S m S ~ l l n l S J ~ U m ~
sub TransfRow2( • Transforms
kth
k,V(1),A(1),xi,X(1),status row
of A matrix.
Gauss
) static
elimination
• R o w 1 is if k > l for c = for
not transformed then i = 1 to kI V(i) : ii = ( i  l ) * i w b +I j = n t o k s t e p i V(j)=V(j)c*A(ii) : ii=ii+l next j xi = x i  c * X ( i ) next i e n d if
' Check pivot less than zero c=V(k) if a b s ( c ) < t o l then status = notOk : exit else status = Ok end if • Normalize row making if k
204
force
diagonal
vector
sub
Adv Enq So/tware, 1991, Vol 13, No 4
sub
coefficient
i.
An e~cient solunon a~ortthm /or boun~ry element equattons J. L Orttz and C. V. Gtroa ~llabhan
sub OpdateIdentity
(k,V(1)
• Updates top left identity ' GaussJordan method
,A(1) , x i , X ( 1 ) matrix
) static
in A m a t r i x .
' If it is r o w 1 we do n o t h i n g if k>l t h e n f o r i = 1 to k1 jj = ( i  l ) * i w b + 1 : c = A(jj+nk) for j=n to k+l s t e p 1 A(jj)=A(jj)c*V(j) : jj=jj+l next j X(i) = X ( i )  c * x i next i e n d if end
sub
sub G a u s s J o r A L G 2 ( X ( 1 ) , V ( 1 ) , i p ( 1 ) , A ( 1 ) , • Solves
A
* U 
X.
• ip() u s e d t o s w a p f o r i1 t o n ip(i) = i next i
U is
rows
returned
when
' Main algorithm. for k1 to n j = k + 1 do call call
StatusOk
FormRow TransfRow2
in
pivot
• Initial dimensions of b l o c k iwb  n1 i h b = int( n m / i w b )
for
... array
is
if p i v o t
X.
close
storing
not
) static
to
zero
A matrix
close
( ip(k),V(),xi, ... ) ( k,V(),A(),xi,X()•status
status ~ Ok then call UpdateIdentity (k,V(),A(),xi,X() call StoreRow2 (k,V(),A(),xi,X() else if j  n + 1 t h e n goto ErrorHandling else jk=ip(k) : ip(k)=ip(j) : ip(j)=jk
to
zero
)
if
e n d if e n d if loop while next k exit
status
: j=j+l
= notOk
sub
• E r r o r h a n d l i n g of s i n g u l a r i t y ErrorHandling: print "Singularity in r o w = " ; k ; " end end
) )
Program
aborted."
sub
Adt Eng, Software. 1991, Vol 13, No 4
205