# 12th GAMM-Seminar Kiel

#### Information about 12th GAMM-Seminar Kiel Mathematisches Seminar Bereich 2 Christian-Albrechts-Universität zu Kiel

The GAMM Committee Efficient numerical methods for partial differential equations in cooperation with the Christian-Albrechts-Universität Kiel organizes the

# Boundary Elements: Implementation and Analysis of Advanced Algorithms

### January 19th to 21st, 1996.

#### Programme for Thursday, January 18th, 1996:

Those arriving in Kiel on Thursday are invited to meet about 19.00 h in a restaurant and pub called 'DÜPPEL 88' (located at Düppelstr. 88, 24105 Kiel).

### Friday, January 19th, 1996

 9.00         Opening
9.15 - 9.55  R. Schneider (Darmstadt):
Multi-Scale Methods for Boundary Integral Equations

Coffee-break

10.30 - 10.55 B.H. Kleemann (Berlin):
Wavelet method for a logarithmic singular integral equation
arising in scattering
11.00 - 11.30 W.S.Hall, R.A.McKenzie (Middlesbrough, U.K.):
A Study of Multi-Wavelet Transforms for a Non-Conforming
Boundary Element Formulation
11.35 - 12.00 A. Rathsfeld (Berlin):
A wavelet algorithm for the BEM corresponding to the fixed
boundary value problem of geodesy

Lunch

14.30 - 15.00 O. Steinbach (Stuttgart):
Fast Solvers for Boundary Element Methods:
Parallelization and Preconditioning
15.05 - 15.15 R. Lehmann (Karlsruhe), R. Klees (Delft, The Netherlands):
Parallel Setup for Galerkin equation system in a
BEM-solution for a geodetic BVP
15.20 - 15.50 A. Greenbaum, A. Mayo, V. Sonnad (Austin, U.S.A.):
Rapid, Parallel Evaluation of Integrals in Potential
Theory on General Three Dimensional Regions

Coffee-break

16.15 - 16.45 M. Kuhn (Linz, Austria):
Parallel Solution of DD-BEM equations using local
Multigrid Preconditioners
16.50 - 17.15 S.A. Funken, E.P. Stephan (Hannover):
Fast Solvers for Adaptive FEM--BEM Coupling
17.20 - 17.50 K.Türke, E. Schnack (Karlsruhe):
A Two Grid Method for Coupling FEM and BEM in Elasticity

->[GAMM-Homepage] ->[Saturday] ->[Sunday]

### Saturday, January 20th, 1996


9.00 -  9.30 S.A. Sauter (Kiel):
Fast Algorithms for Galerkin BEM:
Panel-Clustering and Cubature
9.30 - 10.00 Ch. Lage (Kiel):
Object-oriented design aspects for BEM

Coffee-break

10.30 - 10.55 K. Giebermann (Karlsruhe):
On the Panel Clustering Method for the Helmholtz equation
11.00 - 11.30 C. Gaspar (Gyor, Hungary):
Multigrid and Multipole Techniques in the
Boundary Integral Equation Method
11.35 - 12.05 Y. Yamada (Kyoto, Japan), K. Hayami (Tokyo, Japan):
A multipole boundary element method for two dimensional
elastostatics

Lunch

14.45 - 15.15 B. Faermann (Kiel):
Local a-posteriori error estimators for the discretization
of boundary integral equations
15.20 - 15.40 R. Hochmuth (Berlin):
A-posteriori Estimates for Boundary Elements

Coffee-break

16.20 - 16.50 J.O. Nygaard, J. Grue, H.P. Langtangen,  K. Mørken (Oslo, Norway):
On adaptive spline and wavelet methods
for an integral formulation of inviscid flow
16.55 - 17.35 W.L. Wendland (Stuttgart):
An Extraction and Window Technique for BEM

17.45 - ...   Poster session, Short presentations will be given by:

-  Ke Chen (Liverpool):
Preconditioning Boundary Element Equations
-  T. Finck (Chemnitz):
Spline approximation methods for a class of singular
integral equations over plane domains over plane domains
-  J.-P. Mayer (Kiel):
Solution to Geodetic Problems by Using the Double-Layer
Potential
-  S. Szikrai (Karlsruhe):
A parallel, iterative Neumann solver for
three dimensional elastostatic problems

18.30 Reception

->[GAMM-Homepage] ->[Friday] ->[Sunday]

### Sunday, January 21st, 1996



9.00 -  9.30 Ch. Schwab (Zürich):
Quadrature error analysis for hp-Galerkin BEM on polyhedra
9.35 - 10.00 H. Schippers, F.P. Grooteman (Amsterdam, The Netherlands):
A symmetrical boundary element formulation through cabin walls

Coffee-break

10.20 - 10.50 H. Andrä (Karlsruhe):
A Galerkin-Type Boundary Element Implementation for
3D Elasticity Problems by Using a Computer Algebra System
10.50 - 11.20 Ken Hayami (Tokyo, Japan):
In search of an optimum variable transformation for
nearly singular integrals in BEM

End of 12th GAMM-Seminar Kiel

->[GAMM-Homepage] ->[Friday] ->[Saturday]

# Boundary Elements: Implementation and Analysis of Advanced Algorithms

### A Galerkin-Type Boundary Element Implementation for 3D Elasticity Problems by Using a Computer Algebra System

H. Andrä
Institut für Technische Mechanik/Festigkeitslehre
Universität Karlsruhe (T.H.), Kaiserstr. 12, D-76128 Karlsruhe
Email: heiko@imfsun2.mach.uni-karlsruhe.de

A Galerkin approximation of both strongly and hypersingular boundary integral equation (BIE) is considered for the solution of a mixed boundary value problem in 3D elasticity leading to a symmetric system of linear equations. The evaluation of Cauchy principal values (v.p.) and finite parts (p.f.) of double integrals is one of the most difficult parts within the implementation of such boundary element methods (BEMs). A new integration method, which is strictly derived for the cases of coincident elements as well as edge-adjacent and vertex-adjacent elements, leads to explicitly given regular integrand functions which can be integrated by the standard Gauss-Legendre and Gauss-Jacobi quadrature rules. Problems of a wide range of integral kernels on curved surfaces can be treated by this integration method. We give estimates of the quadrature errors of the singular four-dimensional integrals.
Because the effort for the analytical evaluation of the integrands of the element stiffness matrices increases with the complexity of the kernel functions and with the polynomial order of the shape functions, the computer algebra system Maple V [MAPLE is a registered trademark of Waterloo Maple Software] was employed for this task. A Maple V program greatly simplifies the implementation of this new integration method especially for matrix-valued kernel functions. The computer algebra program produces C or F77 source code for the computation of the integrand functions which can be integrated by standard quadrature rules.
Numerical aspects are discussed.

### Preconditioning Boundary Element Equations

Ke CHEN
The University of Liverpool, England
Email: k.chen@liv.ac.uk

Boundary element equations are linear systems of equations Ax=b with full and dense matrices A. Although wavelet bases may give rise to nearly sparse matrices, we here consider the commonly used boundary element methods based on collocation with polynomial functions or on quadrature.

Many preconditioning techniques have been proposed in the literature for full linear systems, though some are designed with heuristic assumptions or on sparity consideration only. We investigate a class of such preconditioners that may be applied for a fast solution of singular and weakly-singular boundary integral equations. Our results, both theoretical and numerical, show that those preconditioners that contain the essential singularity of a singular operator or the dominant part of a weakly-singular operator are more efficient and reliable. Examples of an exterior Neumann problem and some integral equations with Cauchy singularity are illustrated.

### Local a-posteriori error estimators for the discretization of boundary integral equations

Birgit Faermann
Lehrstuhl Praktische Mathematik, Mathematisches Seminar
Christian-Albrechts-Universität Kiel
Hermann-Rodewald-Str. 3, 24098 Kiel

Adaptive methods are already well established in the finite element field. In boundary element methods, on the other hand, only a few theoretical and numerical results concerning adaptive procedures exist.
The refinement process in adaptive methods is generally controlled by local a-posteriori error estimators. We have shown that some special local a-posteriori error estimators of Babuska and Rheinboldt --which were developed for finite element methods-- are also applicable to boundary element methods. They are applicable to the Galerkin method, if the boundary integral operator has an order alpha in (-1/2,infty) and they are also applicable to the nodal collocation method with odd degree splines.
The Babuska-Rheinboldt error estimators are bounded by local weighted Sobolev-Slobodeckij norms of the residual. These computable local values of the residual are used to control the mesh refinement process.

### Spline approximation methods for a class of singular integral equations over plane domains over plane domains

Tilo Finck
TU Chemnitz-Zwickau, Fakultät für Mathematik, D-09107 Chemnitz

We consider operator equations $(aI+bS_{G})u=f$ in the $L^2(G)$--space, where

  (S_Gu)(x)=-\frac{1}{\pi}\int\limits_G \frac{u(y)}{(x-y)^2}dy

is a Cauchy--type operator over a plane domain $G$ and the coefficients $a$ and $b$ are continuous functions on the closure $\bar{G}$. In his book "'The Method of singular integral equations"', A.D. Dzuraev pointed out the connection between the singular integral operator and certain boundary value problems. We present a uniform Banach algebra approach to the stability problem for a bulk of special approximation methods (including Galerkin and Collocation methods). This approach is based on non--commutative Gelfand theories (local principles). We therefore interprete the stability problem as invertibility problem in suitable constructed Banach algebras. We prove that those approximation method is stable if and only if a certain class of discrete convolution operators and Toeplitz operators on a half plane is invertible. We are able to solve the invertibility problem for the corresponding Toeplitz operators, which are generated by non--continuous functions, up to now for some simple spline functions only. In this case we formulate the necessary and sufficient stability conditions in terms of the coefficients $a$ and $b$ of the considered operator.

### Fast Solvers for Adaptive FEM--BEM Coupling

S.A. Funken & E.P. Stephan
Institut für Angewandte Mathematik
Universität Hannover, Welfengarten 1, 30167 Hannover, GERMANY

We present two nearly optimal preconditioned iterative methods to solve indefinite linear systems of equations arising from h-adaptive procedures for the symmetric coupling nof Finite Elements and Boundary Elements. These solvers are nearly optimal in the sense, that their convergence rate grows only logarithmically with the number of unknowns
They are based either on the conjugate residual method with block diagonal preconditioning, where no Schur Complement construction is required, or on an inner--outer iteration of Axelsson and Vassilevski. Both methods use multilevel additive Schwarz method of seperate positive semi--definite and negative definite parts of the coupled operator.
The efficiency of the solvers is shown by numerical experiments yielding fast convergence.

### Multigrid and Multipole Techniques in the Boundary Integral Equation Method

Csaba Gaspar
senior research associate
Szechenyi Istvan College, Dept. of Mathematics
P.O.Box 701, H-9026 Gyor, Hungary

The computational cost of the traditional version of the Boundary Integral Equation Method is O(N^3) where N is the number of boundary nodes. This is due to the relatively bad properties of the boundary element matrices, since they are generally neither self-adjoint nor sparse. In addition to it, if the original differential equation is supplied with mixed boundary condition, the corresponding boundary integral equation is not of the second kind, so that the traditional well-known iterative techniques can hardly be applied. We present a method which shows some similarities to the multiplicative Schwarz method defined along the boundary. In each iteration step, the boundary integral equations of certain pure Dirichlet and Neumann subproblems are to be solved (even if the original boundary condition is of mixed type), which allows the use of standard multigrid tools. We give a theoretical foundation of the method in Sobolev spaces. Convergence theorems as well as numerical examples are presented, which show that the number of the necessary arithmetic operations can be reduced from O(N^3) to O(N^2). Moreover, we also derive a multipole-based technique to evaluate the appearing boundary integrals in an economic way: the overall computational cost can thus be reduced further to O(NlogN) only.

### On the Panel Clustering Method for the Helmholtz equation

Klaus Giebermann
Universität Karlsruhe, Institut f. Praktische Mathematik
Englerstr. 2, 76128 Karlsruhe

We present an implementation of the panel clustering method for the Helmholtz equation, which occurs for example in acoustic scattering theory. Using the ansatz of Brakhage and Werner, i.e., representing the solution as a combination of acoustic single and double layer potential leads to a Fredholm integral equation of the second kind. This approach avoids problem with wavenumbers close to eigenfrequencys of the Laplace operator. We show the influence of the coupling ansatz to the linear system of equations resulting from disretization procedure, i.e., from the collocation or Galerkin method.
We also describe a method for the fast computation of a wavelet based approximation of the collocation or Galerkin matrix using the panel clustering method.

### In search of an optimum variable transformation for nearly singular integrals in BEM

Ken Hayami
Department of Mathematical Engineering and Information Physics,
University of Tokyo, Japan,
e--mail: hayami@simplex.t.u-tokyo.ac.jp

The author proposed a variable transformation method for the accurate calculation of nearly singular integrals over curved surface elements, which occur in the three dimensional boundary element method when the source point is very near the element [1].
Taking, as an example, a curved quadrilateral element S described by x(eta_1,eta_2), the method first finds the closest point x(overline{eta}_1,overline{eta}_2) on S to the source point x_s. Then, it calculates the point tilde{x}_s = tilde{x}(overline{eta}_1,overline{eta}_2) on the bilinear element tilde{S} defined by the four corner nodes of S. Next, each triangle tilde{triangle}_j, (j=1,...,4) formed by tilde{x}_s and two adjacent corner nodes of the element S, is linearly mapped to the corresponding triangle in the parametric space (eta_1,eta_2). Then, polar coordinates (rho,theta) centred at tilde{x}_s are introduced in each triangle tilde{triangle}_j. Next, a radial variable transformation such as R(rho)=log(rho+d) is introduced in order to weaken the near singularity of the integrand when the source point is near the element, where d is the distance between x_s and S. Also, an angular variable transformation t(theta) is employed to weaken the angular near singularity when the triangle tilde{triangle}_j is thin. Finally, numerical integration is performed in the transformed variables R and t.
In the talk, we will first review the method and its theoretical error estimates, and then introduce some new developments including automatic numerical integration using the trapezium rule and an attempt to further optimize the radial variable transformation for the Gauss-Legendre rule.

[1] K. Hayami, A Projection Transformation Method for Nearly Singular Surface Boundary Element Integrals, Lecture Notes in Engineering, Vol.73, Springer-Verlag, 1992.

### A multipole boundary element method for two dimensional elastostatics

Software Development Center,
Ken Hayami
Department of Mathematical Engineering and Information Physics,
University of Tokyo, Japan
e--mail: hayami@simplex.t.u-tokyo.ac.jp

We will present a multipole boundary element method (MBEM) for two dimensional elastostatics. Unlike the biharmonic equation formulation by Greenbaum et al. [1], we give a direct formulation in terms of displacement and traction variables, which seems more convenient for general boundary conditions and applications.
MBEM formulations are given for the Dirichlet, Neumann and mixed boundary value problems, including the evaluation of the internal stress. At present, the implemented algorithm requires O(N log N) computational work and memory, where N is the number of boundary elements. Theoretical error estimates for the MBEM are also derived.
Numerical examples for Dirichlet, Neumann and mixed boundary problems are given. The GCR and Bi-CGSTAB were used as the iterative solvers for the MBEM. The MBEM was compared with the conventional BEM using the LU decomposition or the same iterative solvers. For the Dirichlet problem, the MBEM with GCR was the most efficient in CPU time for N>300. However, for the Neumann and mixed boundary value problems, the MBEM could not compete in CPU-time, since more iterations were required for convergence. Simple preconditioners such as diagonal scaling and ILU decomposition only gave minor improvements. The MBEM was still the most efficient in terms of elapsed time, which is governed by the memory requirement.

[1] Greenbaum, A., Greengard, L. and Mayo A., On the numerical solution of the biharmonic equation in the plane, Physica D, Vol.60, pp.216-225, 1992.

### A-posteriori Estimates for Boundary Elements

Reinhard Hochmuth
I.Mathematisches Institut, FU-Berlin,
Arnimallee 3, 14195 Berlin

The understanding of a-posteriori error estimates for boundary integral methods appears to be by far less developed than for differential equations. The main reason is that unlike e.g. partial differential operators integral operators possess only pseudolocal properties. The objective of this lecture is to consider a-posteriori estimates in the context of multiscale bases oriented methods. We discuss some basic concepts and ideas from this point of view and relate them to existing techniques in more conventional settings. Our local error estimates arise in a fairly unified fashion essentialy as coefficients of corresponding multiscale expansions. In this way we obtain residual based a-posteriori estimates which are reliable and efficient. A concrete adaptive strategy based on these estimates can be proved to converge. In principle, the results apply to a wide class of elliptic problems covering also operators of negative order.

### Wavelet method for a logarithmic singular integral equation arising in scattering

Bernd H. Kleemann
Berlin Institute of Optics (BIFO)
Rudower Chaussee 5, D-12484 Berlin

3D scattering at cylindrical objects leads e.g. to the exterior Dirichlet problem for the Helmholtz equation in 2D. The problem is transformed into an equivalent boundary integral equation with logarithmic kernel function. For this equation a fully discrete collocation method with biorthogonal wavelets is presented. The arising stiffness matrix is known to be sparse in the wavelet basis. The number of nonzero entries then is of the order O(N logN) with N the number of unknowns. Using an a-priori compression criterion only these entries are necessary to calculate if the stiffness matrix is assembled directly in the wavelet basis. From this a computational amount of the same order does not directly follow because of the necessity of highly accurate quadrature of the arising integrals if the support of the wavelet is near to the support of the test functional. Therefore great attention has to be put to an efficient quadrature within the process of direct matrix assemblation. The rules are to be adapted to the logarithmic singular kernel and to the wavelet structure. The arising sparse linear system is solved iteratively by the Krylov-subspace method GMRES. Because of the logarithmic singularity the condition number grows with N. Therefore a special preconditioner is applied to get a constant iteration count. Results achieved for the accuracy, convergence, compression and time complexity of the method for various scattering objects and wavenumbers are discussed and compared with theoretical ones from the wavelet theory for pseudodifferential equations.

### Object-oriented design aspects for BEM

Christian Lage,
Lehrstuhl Praktische Mathematik, Mathematisches Seminar
Christian-Albrechts-Universität Kiel
Hermann-Rodewald-Str. 3, D-24098 Kiel

Software design for the BEM has to encounter several tasks each with a great variety of parameters. Treating these parameters, e.g. discretization schemes or cubature techniques, too restrictive means developing software adequate only for a few applications, i.e., the reusability of the software or even parts of the software is very limited. To circumvent this nuisance one has to isolate the essential concepts of the BEM as well as their interactions to specify an appropriate and flexible design.

We present an object-oriented approach considering these design rules. The identification of key abstractions like geometry, spaces, dualforms, operators and functions is discussed as well as the extensibility of the developed system to advanced methods, especially the panel clustering algorithm.

### Parallel Setup for Galerkin equation system in a BEM-solution for a geodetic BVP

Rüdiger Lehmann,
University of Karlsruhe, Geodetic Institute, GERMANY
Roland Klees,
Delft University of Technology, Faculty of Geodetic Engineering,
THE NETHERLANDS

A particular problem in setting up large BEM equation systems is the huge number of offdiagonal elements, which in Galerkin method are mainly regular integrals over pairs of boundary elements. The numerical effort of cubature depends strongly on the distance between these elements. Therefore, the polynomial degree of exactness of the cubature formula may be choosen based on an estimate of this distance. In a straightforeward implementation on a parallel computer, this may lead to severe imbalances of workload, which drastically reduces the speedup. To overcome this difficulty, dynamic load balancing schemes are applied. Also, different distributions of boundary elements in the parallel computing environment must be considered, which is equivalent to a static balance of workload. It is shown how this works in the solution of the classical oblique BVP of potential theory (linearized version of the fixed gravimetric BVP) on a IBM 9076 SP/2.

### Solution to Geodetic Problems by Using the Double-Layer Potential

J.-P. Mayer
Mathematisches Seminar, Bereich II,
Christian-Albrechts-Universität Kiel,
Hermann-Rodewald-Str. 3, 24098 Kiel

One problem of physical geodesy is the determination of the shape of the Earth and its gravity field from astrogeodetic and gravimetric data. The mathematical formulation of this problem results in a free boundary value problem (the so-called Molodensky problem) or a mixed, i.e. partly free, boundary value problem (the altimetry-gravimetry problem) for exterior domains. We apply the integral equation method using the double-layer ansatz. Let $D_\varphi^+ (f)$ be the extension onto the boundary of the double-layer potential from the exterior domain. In the same manner we define $\frac{\partial}{\partial l} D_\varphi^+ ( f)$ as the extension of the oblique derivative of the double-layer potential. Here $\varphi$ is a Lipschitz surface and $f$ a density. It is proved that $\frac{\partial}{\partial l} D_\varphi^+$ is invertible in $W_2^{\frac{1}{2}} / \rset$, provided the direction l(x) is not tangential to the surface. For given $G,H: \SS ^2 \longrightarrow \rset$ we consider the problem of seeking a surface $\varphi$ as the solution of

    (*) D_\varphi^+ \left(  \left( \frac{\partial}{\partial l}
D_\varphi^+  \right)^{-1} (H) \right) \ = \ G \qquad .

Knowing that this equation has one and only one solution $\varphi$, we get existence and uniqueness of the linearized Molodensky problem, too. In order to prove local existence and uniqueness of (*), we apply the implicit function theorem. In this light we evaluate the Fr\'{e}chet derivative of (*) w.r.t. surfaces at the sphere. If the arising operator is invertible we can apply the implicit function theorem to get local existence and uniqueness. In this context local means the closeness of the Earth's shape to the sphere. It is planed to use formulation (*) also for constructing an iterative procedure to solve the problem.

### A Study of Multi-Wavelet Transforms for a Non-Conforming Boundary Element Formulation

W.S.Hall, R.A.McKenzie
School of Computing and Mathematics, University of Teesside,
Middlesbrough, Cleveland, TS1 3BA, U.K

We have constructed a non-conforming Boundary Element formulation for 2D Potential problems. The geometry is approximated by quadratic shape functions and the potential and flux are approximated by Legendre polynomials up to degree two. The purpose of which is to produce a linear system of equations which are appropriate for the the application of the Discrete Multi-Wavelet Transform. The solutions of known problems are compared with the approximations produced by this method once truncation has been applied to the transformed matrix. Plots of accuracy against the sparsity are produced showing the trade off between speed of obtaining the solution and the accuracy. The Schultz method is used to solve the system of equations. This method is efficient when the sparsity of the system matrix is $O(n)$. Higher order approximations are thus considered to keep the matrices sparse enough so the Schultz method remains efficient. The benefits of this wavelet based method are compared to more conventional solution procedures. Extensions and generalisations of the work will be discussed.

### Parallel Solution of DD-BEM equations using local Multigrid Preconditioners

Michael Kuhn
Institut f. Mathematik, Johannes Kepler Universität Linz
Altenberger Str. 69, A-4040 Linz, Austria

The Domain Decomposition Method (DD) is a powerful tool for establishing weak formulations and for constructing the corresponding parallel solvers. Although the method allows the coupling of different discretization techniques, i.e. Boundary Element Methods (BEM) and Finite Element Methods, as it is desired in various applications, we restrict ourselves to pure BEM formulations.

In the talk, we give a brief introduction to the preprocessing tools which perform an automatic decomposition of the domain of interest into $p$~subdomains, where $p$ is the number of processors to be used. Furthermore, the parallel algorithm and, in particular, the preconditioners being involved will be discussed in detail. The latter are required for the Schur--complement system and for the single layer potential operator and both are based on local multigrid methods.

Numerical examples, including potential and linear elasticity problems, which demonstrate the high efficiency of the algorithm will be presented.

### On adaptive spline and wavelet methods for an integral formulation of inviscid flow

Jens Olav Nygaard, John Grue, Hans Petter Langtangen, Knut Mørken
University of Oslo, Norway

Traditionally, integral equations have been solved numerically by boundary element methods where the solution is expanded in a sum of lower order piecewise polynomial basis functions along the boundary. It is well known that higher order expansions lead to a faster convergence rate, and hence the possibility of a significant reduction of the computational work in order to reach a particular level of accuracy.
In this paper, higher order boundary element methods based on B-splines and wavelets are described and compared. By numerical examples, it is shown that the process of choosing a good knot vector for the spline case corresponds to the compression of the wavelet matrix. The wavelet method therefore eliminates the possibly tricky procedure of choosing knot vectors for the splines. On the other hand, there are limits to how much the matrix of the wavelet method can be compressed.
The two-dimensional problem is a natural starting point for a thorough comparison, because it involves all the principal mathematical and numerical difficulties. The ultimate aim is to develop better methods for the three-dimensional case.

### A wavelet algorithm for the BEM corresponding to the fixed boundary value problem of geodesy

Andreas Rathsfeld
Weierstraß-Institut für Angewandte Analysis und Stochastik,
Mohrenstr.39, Berlin 10117

One of the fundamental problems of geodesy consists in the determination of the potential field of gravity from the gravity data measured on the surface of the earth. With the help of linearization and boundary element techniques the problem can be reduced to the numerical solution of the singular integral equation

\begin{eqnarray}  \label{1}
-2\pi\cdot\cos [n(x),l(x)]\cdot f(x) &&
\nonumber \vspace*{5mm}  \nonumber
+ p.v.\int_\Gamma
\frac{\cos [l(x),y-x]}{|y-x|^2}\cdot f(y)\cdot d\Gamma (y)&=&4\pi\cdot
\delta g(x).
\end{eqnarray}


Here n(x) stands for the normal at the point x of the surface of the earth Gamma, l is the direction of the reference gravity and the right-hand side delta g is a certain function of the reference gravity and the measured gravity. The solution f is the unknown density function of the single layer representation for the unknown difference gravity potential. Equation ([?]) can be solved numerically by e.g. a piecewise bilinear collocation method. For an accurate approximate solution a high resolution of the geometry (part of the surface of the earth) is necessary. This leads to large and dense linear systems. To save storage and computation time we have implemented a wavelet algorithm for the solution of the collocation system. We will present and discuss numerical experiments with this algorithm.

### Fast Algorithms for Galerkin BEM: Panel-Clustering and Cubature

Stefan A. Sauter,
Lehrstuhl Praktische Mathematik, Mathematisches Seminar
Christian-Albrechts-Universität Kiel
Hermann-Rodewald-Str. 3, D-24098 Kiel

In our talk we will present efficient quadrature strategies for the approximation of the matrix elements which can be applied to any kind of kernel functions, surface approximation techniques, and shape functions. Only the subroutines which evaluates the kernel functions and the shape functions have to be modified. Furthermore this integration technique is relatively easy to implement since it can be tested isolated on simple examples. In the second part of the talk we will formulate the Panel-Clustering method for the Galerkin BEM in a way which is easy to implement. Only the fundamental solution which does not depend on possibly complicated surfaces has to be expanded. The expansions of the kernel functions which are derivatives of the fundamental solution can be derived very fast from that expansion.

### A symmetrical boundary element formulation through cabin walls

H. Schippers, F.P. Grooteman
National Aerospace Laboratory NLR
P.O. Box 90502, 1006 BM AMSTERDAM, The Netherlands

Investigations on the acoustic field around and inside commercial aircraft are motivated by the urge to reduce cabin noise inside the aircraft. The exterior acoustic field causes vibrations of the cabin wall, which contri- bute to the interior noise. The cabin noise inside the aircraft can be reduced by putting porous acoustically absorbing material (insulation) just on the inside of the skin panel of the cabin wall.
In this paper a symmetrical boundary element is described for the trans- mission of sound through a panel of a cabin wall. The geometrical model consists of a double panel configuration including a cavity partly filled with porous material and with air. The lower panel is a stiffened skin panel, while the upper panel is an un-stiffened trim panel. The trim panel is backed by a semi-infinite room modelling the interior of the aircraft. The displacement of the skin panel is assumed to be given by the exterior acoustic field at the outer side of the aircraft. The sound transmission problem requires the calculation of the acoustic pressure in the air and in the porous material as well as the fluid-structure interaction with the un-stiffened trim panel. The acoustic pressure in the air and porous material is modelled by the Helmholtz equation (in the case of glass wool insulation with a complex wavenumber, which models the damping). The porous material is essentially assumed to be described by its porosity, resistivi- ty, effective density and effective speed of sound, all of which are frequency dependent.
The acoustic pressure in the vibrating medium (air and porous material) is modelled by a boundary integral ansatz. The vibrating panels and the vibrating medium are coupled by the boundary condition for flexible walls. Classical boundary integral formulations result in a non-symmetric boundary element matrix, while the discretization of the Helmholtz equation using finite elements would yield a symmetrical stiffness matrix. In the present paper a symmetrical boundary integral formulation is described for the acoustic pressure inside and outside the cabin wall. The formulation maps the normal derivative of the pressure along the boundary on the pressure itself. In the case of a bounded domain (e.g. the volume covered by the porous medium or the volume covered by the air inside the cabin wall) the symmetrical boundary integral formulation has been taken from $[1]$. The symmetrical boundary element matrix follows from applying Galerkin's method. The present formulation is advantageous from the point of view that the boundary element matrix has the same mapping properties as the finite element matrix of the weak formulation of the Helmholtz equation. These properties are not preserved by the classical boundary element formulation, which yields an asymmetric system matrix. The basis and test functions are given by the classical piecewise linear functions on flat non-overlapping triangular elements.
The boundary elements are put on the panels of the cabin wall and on the interface between the air and the porous material. The distance between the panels and the interface is small, which requires an accurate evaluation of the influence coefficients of the boundary element matrix. Cubature formulas for the nearly singular surface integrals will be discussed.
The computations involve the discretization of a hyper-singular integral operator, a weakly singular integral operator and a regular operator. The hyper-singular integral operator is regularized by integration of parts. Once this operator has been regularized, the computation of the coeffi- cients of the corresponding boundary element matrix requires the same number of kernel evaluations as the computation of the coefficients of the weakly singular operator. Therefore, the computational costs of the present symmetrical boundary element formulation are comparable with the costs of classical boundary element formulations in acoustics, which are based on the direct method.

[1] SCHMIDT, G.: Boundary element discretization of Poincare-Steklov operators, Numer. Math., Vol. 69, 1994, pp. 83-101.

### Multi-Scale Methods for Boundary Integral Equations

Reinhold Schneider
Fachbereich Mathematik

Boundary element methods offer an appropriate tool for the numerical solution of certain boundary value problems in engineering. A major drawback of BEM is the fact that the arising system matrices are densely populated, which is limiting the discretization of realistic 3D problems with complex geometry. Like panel-clustering and multi-pole expansion, biorthogonal wavelet bases remedy this situation by approximating the discrete scheme in an efficient way. Multi-scale methods achieve this by approximating the system matrix relative to a biorthogonal wavelet basis by a sparse matrix. We propose a boundary element discretization for static or low frequency problems which is based on parametric surface representation. The surfaces are assumed to be piecewise analytic and parametrized over quadrilateral surface (macro-)elements. Trial functions are supposed to be globally continuous and piecewise bilinear in each parameter domain. We construct a biorthogonal wavelet basis with desired vanishing moments on each surface (macro-)element. In addition, this basis satisfy a stability condition so that the Galerkin wavelet method can be immediately preconditioned. The proposed basis have all desired properties achieving an asymptotically optimal balance between accuracy and efficiency. This means that we directly compute a sparse approximation of the system matrix causing an additional error proportional the optimal error bound of the Galerkin discretization. This can be done such that the total number of nonzero matrix coefficients increases only linearly with the total number of unknowns $N$. In order to compute the nonzero matrix coefficients directly, we propose an adaptive quadrature method, with the desired accuracy requiring totally O (N) floating point operations.
An outlook over essential implementation features and particularities of wavelet basis approach including pre-processing for surface and mesh generation, hierarchical data structure, a priory compression strategy, adaptive feedback matrix computation or oracle, adaptive boundary element methods and application to post-processing for field computation close and far from the boundary should be appended.

### Quadrature error analysis for hp-Galerkin BEM on polyhedra

Christoph Schwab
Seminar fuer Angewandte Mathematik, ETH Zürich,
Rämistrasse 101, CH-8092 Zürich

We analyze the quadrature error for a hp-Galerkin BEM on curved polyhedra. The method is based on subspaces of high degree polynomials combined with geometric mesh refinement towards the edges of the domain. Under suitable regularity assumptions on the solution generally believed to be valid (and proved in certain special cases by ELSCHNER and STEPHAN and coworkers) the method converges exponentially PROVIDED the entries in the stiffness matrix are evaluated exactly.
We propose a numerical quadrature scheme for the fully discrete evaluation of the singular and near singular integrals in the stiffness matrix. We pay particular attention to the numerical evaluation of the high aspect ratio elements arising in the vicinity of the edges of the polyhedron.
We show that using certain regularizing coordinate transformations combined with proper geometric subdivisions of the domains of integration all entries in the stiffness matrix can be computed to the required exponential accuracy with work that grows only algebraically in terms of the number of degrees of freedom.
We present numerical computations which support the theoretical predictions.

### Rapid, Parallel Evaluation of Integrals in Potential Theory on General Three Dimensional Regions

Ann Greenbaum, Anita Mayo, Vijay Sonnad

We present a new, high order accurate method for the rapid, parallel evaluation of certain integrals in potential theory on general three-dimensional regions. The kernels of the integrals are a fundamental solution, or a linear combination of the derivatives of a fundamental solution of some second order linear elliptic differential equation. What is different and important about these methods is that they avoid the the use of any standard quadrature formula or the exact evaluation of any integral. Instead, they rely on rapid methods of solving the differential equation which the kernel satisfies. In fact, the number of operations needed to evaluate the integrals is essentially equal to the number of operations needed to solve the differential equation on a uniform cubic grid.
In particular, one can evaluate integrals whose kernels are the Green's function for Poisson's equation by using Fourier methods on a cubic grid, or a fast Poisson solver. Furthermore, the method requires no evaluation of the kernel. Before applying the Poisson solver one only needs to compute special correction terms to the right hand side of the Poisson equation at mesh points near the boundary of the irregular region, and these correction terms can be obtained by local calculations. (We note that it is not possible to to evaluate these integrals by straightforward use of Poisson solvers since they have discontinuities in their second derivatives across the region of integration.) The method we present can also be thought of as a way of solving elliptic differential equations whose solutions and gradients are continuous, but which have discontinuities in their second derivatives across some irregular boundary.
In the talk we present results of numerical experiments, including some on the IBM SP2 parallel computer.

### Fast Solvers for Boundary Element Methods: Parallelization and Preconditioning

Olaf Steinbach
Mathematisches Institut A
Universität Stuttgart

The discretization of mixed boundary value problems by Galerkin or collocation boundary element methods leads to a system of linear equations with a dense matrix, whose spectral condition number depends from the mesh size h.
To construct efficient iterative solvers the use of good preconditioners is necessary. We will give a general approach based on the Calderon projector related to the problem. However, the construction of fast solvers for mixed boundary value problems depends strongly on the boundary integral formulation, which is in general not unique.
Due to the effort of integration of singular integral operators a parallelization of these methods seems to be more suitable. Here we consider algebraic block decompositions of the matrix as well as domain decomposition methods. In both cases we have to discuss related preconditioners.
Some numerical results concerning the potential equation and problems in linear elasticity will be given.

### A parallel, iterative Neumann solver for three dimensional elastostatic problems

Szabolcs Szikrai
Institute of Solid Mechanics, University of Karlsruhe

In the FE/BE coupling method by Schnack et al. ([1],[3]) we need a numerical technique to map the Neumann data (traction field) into the Dirichlet data (displacement field) on the coupling boundary. The basis of the method is the Calderon projector with the strong singular boundary integral equation of the linear elasticity. To discretize the problem, we use the collocation method. As the linear system of algebraic equations is very sensitive to small perturbations (e.g. roundoff error) in the matrix and in the right--hand side, we must resort to special techniques. Türke [3] developed an iterative algorithm to solve the linear system based on the Neumann series. In the lecture we present the efficient parallelized version of this iterative method. The algorithm is well suited for MIMD parallel architectures. We examine the convergence properties and compare some message passing communication structures on several test problems.

[1] SCHNACK E.: "Stress Analysis with a Combination of HSM and BEM" Proceedings of the MAFELAP 1984 Conference on "The Mathematics of Finite Elements and Applications", Uxbridge, 1-4 May 1984. Edited by J. R. Whiteman, Academic Press 1985, pp. 273-281.

[2] KARAOSMANOGLU N.: "Kopplung von Randelement- und Finite Element Verfahren für dreidimensionale elastische Strukturen" Ph.D. thesis at the Institute of Solid Mechanics, Karlsruhe University (1989).

[3] TÜRKE K.: "Eine Zweigitter-Methode zur Kopplung von FEM und BEM" Ph.D. thesis at the Institute of Solid Mechanics, Karlsruhe University (1995).

### A Two Grid Method for Coupling FEM and BEM in Elasticity

K. Türke, E.Schnack
Institute of Solid Mechanics, Karlsruhe University

A coupling algorithm of FEM and BEM for solving mixed boundary value problems in elas\-tostatics is presented. A domain decomposition of the bounded domain Omega leads to a basic substructure Omega_1 with a skeleton H, where the FEM is applied, and several macroelements Omega_i (i=2,...,p), where the BEM is used on a fine grid h. As fundamental equations the well known energy bilinear form for the whole domain Omega and, additionally, a second bilinear form on the macroelements Omega_i as a coupling condition is used. This way a symmetric and nonconforming -- that means different, independent grids H and h -- coupling algorithm can be obtained. This can be realized by applying the Poincaré--Steklov operator on the macroelement surfaces in the strong singular form and by avoiding the use of the hypersingular boundary integral equation.
The construction of a robust and reliable numerical algorithm depends on the adaptive control of symmetry and definiteness of the coupling matrix. Therefore we use an iterative method for solving the boundary integral equation by expanding the {\em Calderon} projector in a Neumann series. The proof of convergence for this expansion is based on the fundamental work of Kupradze [1].
Numerical results in 2D and 3D will show the preciseness and efficiency of the method.

[1] Kupradze V.D., T.V. Burchuladze, T.G. Gegelia, M.O. Bacheleishvilii: Three dimensional problems of elasticity and thermoelasticity. North--Holland, Amsterdam (Nauka, Moscow, 1976).

### An Extraction and Window Technique for BEM

W.L. Wendland
Universität Stuttgart

The boundary integral equation for the missing Cauchy datum can be used to compute arbitrary derivatives of the original solution with a boot--strapping algorithm. A related algorithm leads to a window technique which can be used for local error indicators as well as for overlapping Schwarz methods. The lecture presents joint work with C. Schwab and H. Schulz.

### Matrix free iterative solution strategies for 3d industrial BEM applications

A. Yeremin
Institute of Numerical Mathematics
of the Russian Academy of Sciences

The principal bottleneck of many 3D industrial boundary integral applications is related to necessity to store using O(n^2) memory locations and to solve by a direct solver using O(n^3) arithmetic operations large dense linear systems. We suggest a so-called algebraic matrix free iterative solution strategy with an O(n^alpha) alpha < 2, serial arithmetic complexity which requires only O(n^beta), beta < 2, memory locations. The main idea is to replace the original linear system AX = B by the so-called compressed linear system tilde{A} tilde{X} = B, where only O(n^beta), beta < 2, memory locations is required to store tilde{A} and only O(n^alpha), alpha < 2, multiplications is required to compute tilde{A} Y. At the same time we quarantee that ||tilde{X} - X||/||X|| < epsilon for any prescribed epsilon if ||A - tilde{A}||_F < = epsilon_A ||A||_F, where tilde{A} is a block low rank matrix approximation. This strategy allows to solve on existing supercomputers very large complex linear systems of sizes up to several hundreds thousands. For example, the complex dense linear system of size n = 80802 with 722 multiple right hand sides originated from the standard Metallic NASA Almond CEM benchmark can be solved on 1 CPU of CRAY C90 in 25h 20min by the matrix free iterative solver LRA_CDENSE using only 147 MWords of the main memory and 21.2 GBytes of the disk space. The incore direct solver LAPACK will require more than 104 GBytes and more than 390h (assuming its peak performance on CRAY C90). The real dense unsymmetric linear system of size n = 154089 with a single right hand side originated from the industrial CFD benchmark can be solved on 1 CPU of CRAY C90 in 7h by the matrix free iterative solver LRA_RDENSE using only 103 Mwords of the main memory and 13 GBytes of the disk space. The incore direct solver LAPACK will require more than 192 Gbytes and more than 685h.