# An efficient solution algorithm for boundary element equations

## An efficient solution algorithm for boundary element equations

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, ...

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 AX--F. 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*(n-1)/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 i-th 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 i-th 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

3-rd eqn.

.

.

.

01 001 000 000

.

1 01

0 0 11

0 0 0

I-th 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. (I-1)ROW\$ TRANSFORMED I-th 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 i-th 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 n-th 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

1

m=~=0

re=s=0

d=In-l)

d=p

2 3 4 5

m = ~ = ( n - I) m=~=ln-I)+tn-2t m=~=(n-I)+fn-2l+ln

m=s=t*p rn=~=2*p m=~=3*p

d=(n-21 d=tn-3)

d=p d=p

m=~=In-2)*p m=~=(n-I)*p

d=l d=0

n-I n

m=s=(n-I)+(n-2)+ m=s=(n-I)+(n-2)+

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

--

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=ln-l)*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

upper-Identity matrix

Ro~ Rov~ o f m a t r i x

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=(11-li*l

117=s=(n--2)*p

d=l

m=s=(n-I)*p

d=O ( n - 1)*n d = - -

I1-1 i1

-

~=0

.-1 Sum

m=~=

111- l)*n*p t*(n--r~

~

nl=~

i=1

[

Total sum

NOTE

n

d-

=

~

( n - I-)-* n "1

number

p = number

+ n*p

1

+

[

m

Rows of matrix

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"(n-4)

m = s = l*p m = s = 2*p m = s = 3*p

d=p d=p

m=s=ln-2)*l m=,,=O

m=s=ln-2)*p m=s=(n-I)*p

d=n*p

re=s=

.-2

"

n-1

n-2

i=1

i=l

~

( n - I)*n*p t*(n-l-l)

m=s

1=1

= s= ~" t'In--t)+ ~- i * ( n - - i - I ) + l n - - l l * n * p

-

1

= sums/subtractions d = divisions m = multlphcanons

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 Gauss-Jordan 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 Gauss-Jordan 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 i-th row of matrix A as in the case of

II

0 o I STORED AS ONE ~0 1 0 DIMENSIONAL 0 0 1 ARRAY 4-th 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

i-th 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 0-1) ROWS TRANSFORMED i-th ROW READY FOR GAUSS JORDAN ELIMINATION

FZ~ 2

b 0-1)ROWS TRANSFORMED i-th 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 back-substitution 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 Gauss-Jordan method involves the following summations: .-1

n-2

l*(n--i)+ ~ 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. McGraw-H~ll Book Co. New York. 1981 3 Brebbla. C A. Telles, J C F and Wrobet L C Boundar~ Element Teehmques. Sprmger-Verlag. Berhn. 1984 4 Brebbta. C A and Dommguez, J Boundary Elements. ~m Introductory Course. Computational Mechamcs Pubhcatlons and McGraw-Hall Book Co. New York. 1989 5 Au. M C Boundary element dLrect ehmmatton procedure m multi-region problems Proceedings o! the Boundar~ Element oth Internattonal Conference, 1984

n-I

i*[n-z-1)=

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*(n-1)/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 a-h,o-z :defint i-n 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

k-th

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 ) - (k-l) * ( 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

k-th

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 k-1 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 n-1 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 a-h,o-z :defint i-n 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

k-th

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

=

(i-l)*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

k-th

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 k-I 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 ' Gauss-Jordan 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 k-1 jj = ( i - l ) * i w b + 1 : c = A(jj+n-k) 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 i-1 t o n ip(i) = i next i

U is

rows

returned

when

' Main algorithm. for k-1 to n j = k + 1 do call call

Status-Ok

FormRow TransfRow2

in

pivot

• Initial dimensions of b l o c k iwb - n-1 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