HILBERT — a MATLAB implementation of adaptive 2D-BEM

We report on the Matlab program package HILBERT. It provides an easily-accessible implementation of lowest order adaptive Galerkin boundary element methods for the numerical solution of the Poisson equation in 2D. The library was designed to serve several purposes: The stable implementation of the integral operators may be used in research code. The framework of Matlab ensures usability in lectures on boundary element methods or scientific computing. Finally, we emphasize the use of adaptivity as general concept and for boundary element methods in particular. In this work, we summarize recent analytical results on adaptivity in the context of BEM and illustrate the use of HILBERT. Various benchmarks are performed to empirically analyze the performance of the proposed adaptive algorithms and to compare adaptive and uniform mesh-refinements. In particular, we do not only focus on mathematical convergence behavior but also on the usage of critical system resources such as memory consumption and computational time. In any case, the superiority of the proposed adaptive approach is empirically supported.


Mathematics Subject Classifications (2010) 65N38 • 65Y20 • 65N50 1 Introduction
In many applications, the (Galerkin) boundary element method (BEM) has established itself as a possible strategy for the numerical solution of a certain class of partial differential equations.One of the strengths of BEM is its potentially high order of convergence: For smooth analytical solutions, the error of the lowest-order Galerkin BEM behaves like O(h 3/2 ) with respect to the natural energy norm, where h denotes the global mesh-size of a partition of the boundary.Moreover, by means of the representation formula, one can obtain high-order pointwise approximations of the solution of the underlying partial differential equation.However, due to generic singularities of concentrations and/or fluxes on the boundary, these high convergence rates are not observed with uniform meshes in practice.
The MATLAB program package HILBERT has been designed to make adaptive BEM (ABEM) more accessible to a broader audience.It provides stable implementations of the discrete boundary integral operators corresponding to the Laplace operator in 2D.Various error estimators, including the weighted-residual error estimator from [11,12], (h − h/2)-type error estimators from [5,16,21], and two-level estimators from [25,27,31,32] are provided.In addition, mesh-refining algorithms, functions for the discretization of given boundary data, visualization routines, and many other functions necessary for the implementation and use of ABEM for the Poisson problem are included.
Several benchmarks demonstrate and compare the behavior of adaptive and uniform BEM in presence of different types of singularities.In our experiments, we observe that our adaptive algorithm resolves the singularities of both, analytical solution and given data.Throughout, it shows the optimal order of convergence.Moreover, the reachable accuracy-usually limited by system resources such as time and memory-is increased drastically by use of adaptivity.

Model problem
As model problem, we consider the Poisson equation on some bounded Lipschitz domain ⊂ R 2 with polygonal boundary = ∂ .By scaling of the domain , we assume diam( ) < 1 to ensure ellipticity of the single-layer integral operator V , see Section 4 below.Each solution u ∈ H 1 ( ) can explicitly be written in the form where g := u| ∈ H 1/2 ( ) := v ∈ L 2 ( ) | there is an extension v ∈ H 1 ( ) such that v| = v is the trace of u and φ := ∂ n u ∈ H −1/2 ( ) := H 1/2 ( ) * is the normal derivative of u on .For this fact and all other statements and details on the spaces and integral operators involved, the reader is referred to the literature, e.g., the monographs [30,34,35,38].For x ∈ , the involved linear integral operators read where n y denotes the outer unit vector of at some point y ∈ and where d (•) denotes integration along the boundary.Put differently, the solution u of the Poisson problem ( 1) is known as soon as its Cauchy data (u| , ∂ n u) are known on the entire boundary .
If one considers the trace and the normal derivative of u, the so-called representation formula (2) becomes the Calderón system It involves six linear integral operators: The single-layer integral operator V , the double-layer integral operator K with adjoint operator K , the hypersingular integral operator W , which is the negative normal derivative of K, see (39) below, as well as the trace N 0 and the normal derivative N 1 of the Newtonian potential N .We stress that N 0 , V , and K formally coincide with the operators N , V , and K from (3)-( 5), but are now evaluated for x ∈ instead of x ∈ .Moreover, the jump ± 1 2 in (6) holds only almost everywhere on , more precisely, in all x ∈ but corners, where 1 2 becomes 1 − σ (x) in (6a) and σ (x) in (6b) with σ (x) := α 2π and α being the interior angle at x.While this is important for collocation BEM, one can use the value 1   2   for Galerkin BEM throughout.
For direct BEM, the Poisson equation with given boundary data is equivalently stated in terms of the Calderón system (6), which leads to a boundary integral equation (BIE) on .In order to solve this BIE, we discretize a variational formulation by a Galerkin method.To that end, let E 1 , . . ., E N be a partition of the boundary into closed affine line segments.We then use piecewise constant functions to discretize the fluxes ∂ n u and piecewise linear and globally continuous functions for the concentrations u| .In a post-processing step, the computed Cauchy data are then plugged into the representation formula (2) to obtain an approximation of the solution u of the Poisson problem (1).

Installation & basic structure of HILBERT
HILBERT is available from Netlib (http://www.netlib.org/numeralgo/)as the na38 package.It is distributed as a single file hilberttools.zip.To install the library, download and unpack the zip archive, start MATLAB, change to the root folder hilberttools/ and type make at the MATLAB prompt.You may test the installation by running some examples from the demos/ or paper/ folder.
The MATLAB functions provided by HILBERT are contained in the root folder hilberttools/, while the subfolders demos/laplace/ and demos/poisson/ contain demo implementations of adaptive algorithms for different model problems.The subfolder paper/ contains functions to reproduce the figures of this paper.The implementation of the integral operators is done in C via the MATLAB MEX-Interface.The C sources are found in source/, see also Section 3.3 below.
The file documentation.pdf in the root folder hilberttools/ contains a detailed description of HILBERT and its usage.The documentation reports on the aspects of the implementation and provides mathematical background and references.It may therefore serve as users and developers handbook.
HILBERT should compile without problems using MATLAB versions ≥ 7.4 and compatible ANSI-C compilers.HILBERT is known to compile on Mac OS X and Linux using the GNU Compiler Collection (gcc) and on WINDOWS using Microsoft's Visual C Compiler.Furthermore, it is known to compile on WINDOWS using MATLAB's LCC compiler, if multithreading is disabled.Further information on the installation is also provided by README.txtwhich is found in the HILBERT root folder hilberttools/.

Outline of the paper
First, Section 2 deals with an abstract formulation of the proposed adaptive algorithm.Preliminaries of practical relevance, such as mesh administration with HILBERT and comments on the stable implementation of the discrete boundary integral operators, are given in Section 3. In the following Sections 4-5, we give possible choices of error indicators for the boundary integral formulation of the Dirichlet and the Neumann problem.Finally, Section 6 deals with a mixed boundary value problem and non-homogeneous volume forces f = 0.Sections 4-6 are supported each by its own set of numerical experiments.A discussion and a summary of the functionality of HILBERT finally conclude the paper in Section 7.
Throughout, we use the notation A B to indicate that expression B gives an upper bound A ≤ C • B up to some generic constant C that is independent of relevant unknown data.The notation A B is used to express A B A.

Adaptive algorithm
Using a lowest-order ansatz, BEM is known to converge optimally with order O(h 3/2 ) with respect to the energy norm if the data and unknown solution are sufficiently smooth, e.g., have piecewise regularity φ ∈ H 1 , g ∈ H 2 , f ∈ H 1 [34].However, fluxes and concentrations on the boundary usually lack the regularity necessary to observe the optimal order with uniform meshes.We therefore aim at providing a fully adaptive mesh-refining algorithm that leads to more efficient use of available system resources.

Formulation of adaptive algorithm
Let H be a Hilbert space with corresponding norm ||| • |||, which is made more precise in the Sections 4-6 below.For discrete subspaces X corresponding to successively generated meshes E of , let U ∈ X denote the computed discrete approximation of the unknown solution u ∈ H. Furthermore, we assume that we can compute a posteriori error bounds ρ = E∈E ρ (E) 2 1/2 such that the local contributions ρ (E) measure, at least heuristically, the error between u and U locally on the element E ∈ E .Under these assumptions and given an initial mesh E 0 as well as an adaptivity parameter 0 < θ < 1, the usual adaptive algorithm reads as follows: (iv) refine (at least) elements E ∈ M and obtain E +1 (v) update counter → + 1 and goto (1)

(h − h/2)-type error estimators
Let E denote the uniform refinement of E which is obtained by splitting all elements E ∈ E at their midpoint.The associated discrete space is denoted by X .Under the Céa-type quasi-optimality which is met in the applications at hand due to the uniform stability of the Galerkin BEM, the Galerkin solution U corresponding to the space X is a better approximation to the unknown solution u ∈ H than U .In [21], we first analyzed the simple in the context of BEM.It always provides a lower bound for the Galerkin error Moreover, an upper bound |||u − U ||| η (11) follows from the so-called saturation assumption with some uniform constant q ∈ (0, 1).I.e., uniform mesh-refinement leads to a uniform improvement of the discretization error.From a heuristic point of view, the saturation assumption states that the BEM error exhibits asymptotics O(h α ) for some α > 0, see the discussion in [21].Note that the saturation assumption (12) is even equivalent to the upper bound (11) in case of symmetric problems, where ||| • ||| stems from an energy scalar product.Moreover, we stress that (12) can, for instance, be proved for finite element model problems and up to data approximation terms [14], see also [20,Section 2.3].For BEM, however, (12) still remains mathematically open although observed throughout academic experiments.In [3, Prop.A.2], a weaker form of ( 12) is verified for Symm's integral equation of Section 4 if E is obtained from k ≥ 2 uniform refinements of E : There exists a constant k ∈ N ≥2 which depends only on and the given data, such that the saturation assumption (12) holds up to higher-order terms which involve the smooth parts of the solution u only.
Due to the non-local nature of the integral operators involved, the energy norm ||| • |||, which is equivalent to a fractional Sobolev norm • H ±1/2 ( ) (or a sum of those for mixed boundary value problems), cannot be easily written as a sum of local contributions, as is, e.g., the case for the L 2 -norm . Based on certain localization techniques, the works [8,9,15,16,21], for instance, give suitable error estimators μ 2 = E∈E μ (E) 2 which satisfy η μ , but whose contributions μ (E), at least heuristically, measure the local error.Moreover, based on recent results from [5], a compound error estimator ρ which additionally controls data oscillations, is used for all problem settings discussed in this paper.

State of the art
In practice, Algorithm 1 yields a sequence of discrete solutions U ∈ X which converge towards the exact solution u ∈ H in the energy norm.Usually, one empirically even observes a quasi-optimal decay of the energy error |||u − U ||| with respect to the number of elements #E .However, since Algorithm 1 does not enforce that the meshsize tends to zero everywhere, even plain convergence lim U = u is a mathematical issue.In this section, we aim at a short overview on the state of the art of convergence of adaptive BEM.We stress that the numerical analysis is still in its infancy in this respect, and all results are so far tailored to the energy norm.
For Symm's integral equation and lowest-order discretization, the first convergence result on adaptive BEM was [10], where Algorithm 1 is extended and the (h − h/2)-error estimator of [21] and the weighted-residual error of [11,12] are combined to ensure the saturation assumption (12) algorithmically.In [20], we prove that under the saturation assumption (12), Algorithm 1 driven by the (h − h/2)-error estimator μ yields contraction of a quasi-error quantity μ 2 in each step of the adaptive loop.
In [6], we introduced the estimator reduction principle to prove that, independently of the saturation assumption (12), Algorithm 1 yields estimator convergence lim μ = 0.In [5], we consider (h − h/2)-error estimation for a mixed boundary value problem with data approximation.Using the estimator reduction principle, we prove that Algorithm 1, driven by a compound error estimator ρ which additionally controls data oscillations, guarantees lim ρ = 0. Here, given Dirichlet data are discretized by nodal interpolation, whereas given Neumann data are discretized by the L 2 -projection onto piecewise constants.In [23], we extend the concept of fullydiscrete adaptive BEM to 3D and the discretization of Dirichlet data by means of the L 2 -projection onto the continuous and piecewise affine functions.Dealing with discrete integral operators only, i.e. matrices, will be a starting point to combine adaptive BEM with fast methods for the discrete integral operators, see e.g.[33] for a recent overview on fast methods.This will, however, be a challenging topic for future research.
Very recently, it has been proved that also weighted-residual error estimators yield contraction of an appropriate quasi-error in each step of the adaptive loop.Moreover, for sufficiently small 0 < θ < 1, i.e. strongly adapted meshes, it is proved that adaptive BEM leads to quasi-optimal decay of the error with respect to certain approximation classes, see [19] for Symm's integral equation of Section 4 and [36] for the hypersingular integral equation of Section 5, both for the respective lowestorder BEM.Furthermore, for technical reasons, the analysis of [36] is restricted to C 1,1 -boundaries which formally excludes polygonal geometries.In [3], the approximation class for the lowest-order BEM for Symm's integral equation in 2D is characterized in terms of the energy error |||u − U ||| only.In particular, we derive that each possible convergence rate 0 ≤ s ≤ 3/2 is in fact achieved by Algorithm 1 driven by the weighted-residual error estimator.
In the work [4] the heart of the analysis of [19,36], namely certain novel inverse-type estimates for the boundary integral operators involved, are generalized to arbitrary polynomial degree and polygonal (resp.polyhedral in 3D) boundaries.Therefore, convergence and quasi-optimality of adaptive BEM driven by weightedresidual error estimators also hold for general polynomial ansatz spaces.Since the boundary integral formulation of the mixed boundary value problem (see Section 6) involves a non-symmetric operator, a quasi-optimality result remains open in this context.
Combining the techniques of [4,19,36] with the data approximation of [5,23], convergence of a fully discrete adaptive BEM follows, where Algorithm 1 is driven by the sum of a weighted-residual error estimator and appropriate data oscillation terms.Moreover, even quasi-optimal rates can be guaranteed for polynomial ansatz of arbitrary, but fixed polynomial order [17,18].

Benchmarks
To study the performance of our adaptive approach compared to a uniform meshrefining strategy empirically, we propose two benchmarks.The first example covers the Laplace equation, i.e. f = 0, and is constructed in such a way that Dirichlet and Neumann data, both, have singularities.In particular, they lack smoothness properties sufficiently to reveal optimal convergence with uniform mesh refinement.The second example covers the case of nontrivial volume forces f = 0 and is designed such that our adaptive algorithm needs to resolve singularities of Cauchy and volume data.
Below, the results of all numerical experiments are visualized with three figures.Since we prescribe the analytical solution, a reliable error bound can be computed and is shown for reference along with the error estimators and data oscillations.In the first figure, we plot all quantities over the number of boundary elements #E .We recall that the optimal rate of convergence of lowest-order BEM is O #E −3/2 , since h ∼ (#E ) −1 for uniform meshes.Second, we plot the quantities over the computational time.Since an adaptively generated solution depends on the entire history of solutions, whereas this is not the case for uniform meshes, the time consumption is measured differently for the uniform and the adaptive approach.We define the computational time as follows: -For uniform mesh-refinement, t ,unif is the time elapsed for uniform refinements of the initial mesh E 0 , the assembly of the Galerkin data on E , and the computation of the Galerkin solution U .
-For ≥ 0, we set t ,adap = t −1,adap plus the time elapsed for the assembly of the Galerkin data on E , the computation of the Galerkin solution U , the computation of the error indicators, and the adaptive refinement of the mesh E to obtain E +1 .
Finally, in a third figure, we plot the quantities over the memory consumption which is understood as follows: -For uniform mesh-refinement, we count the memory which is occupied by the data structure for the mesh, the discrete integral operators, and the solution vector.-For the adaptive version, we count the memory which is occupied by the data structure for the mesh, the integral operators, the solution vector, the error estimators, and the data oscillations.
With these three different figures, we empirically evaluate the quality of a performed computation with respect to both, mathematical order of convergence and computational effectivity.Throughout, the numerical experiments were conducted on a common Linux PC (Ubuntu 12.04, 64 Bit) with 32 GB of RAM and Intel i7-3930K processor with 6 cores running at 3.2 GHz.We used MATLAB R2012a (release 7.14.0.739).MATLAB scripts which reproduce the Figs.2-6 from the experiments below, are part of the HILBERT library and found in the path paper/.

Geometry
We choose to be a rotated L-shaped domain with diam( ) < 1 as shown in Fig. 1.The initial boundary mesh E 0 consist of #E 0 = 8 elements.In the case of nonvanishing volume forces, we discretize the given data f ∈ L 2 ( ).To that end, we use a regular triangulation T of that satisfies T | = E and approximate f by its
An analytical solution u of the Laplace equation where z denotes the uppermost corner of .This choice of u effects that the Dirichlet data g has a singularity at the uppermost corner, whereas the Neumann data φ has an additional singularity at the reentrant corner, as in general arises for the reentrant corner with angle 3π/2.

Example 2
Let z = (0.14, 0.14) ∈ denote a point inside the domain.We then prescribe the analytical solution The volume force f = − u has a singularity at the point z, and the normal derivative A triangulation or mesh of is a finite set E = {E 1 , . . ., E N } such that the elements E j ∈ E are closed affine line segments and their intersection has vanishing measure, i.e., it holds that E j = [a j , b j ] := conv{a j , b j } for a j , b j ∈ as well as Furthermore, if is partitioned into D and N , one usually assumes that this partition is resolved by E , i.e., E j ∈ E satisfies either E j ⊆ D or E j ⊆ N .With K = {z 1 , . . ., z N } the set of all nodes of the triangulation E , it holds that #E = #K for the closed boundary .

Data structure
The set of nodes K = {z 1 , . . ., z N } of the triangulation E = {E 1 , . . ., E N } is represented by an N × 2 array coordinates.The j -th row of coordinates stores the coordinates of the j -th node z j = (x j , y j ) ∈ R 2 as coordinates(j ,:) = [ x j y j ].
If is not split into several parts, the triangulation E is represented by an N × 2 array elements.The i-th boundary element where the nodes are given in counterclockwise order, i.e., the parametrization of the boundary element E i ⊂ is mathematically positive.This assumption allows to compute the exterior unit normal vector of on each element E i ∈ E in a unique way.
If is split into Dirichlet boundary D and Neumann boundary N , the triangulation E is represented by an N D × 2 array dirichlet and an N N × 2 array neumann which describe the elements E i ⊆ D and E i ⊆ N , resp., as before.Then, elements = [dirichlet;neumann] with N = N D + N N .

Mesh refinement
For boundary element meshes, HILBERT includes an optimal local mesh-refinement function refineBoundaryMesh, where the sense of optimality is explained in the subsequent Section 3.2.1 below.Let marked be an (M × 1)-column vector containing the indices of marked elements, i.e. the set M in step (3) If an element E i ∈ E is not refined, one has E i = e j ∈ E +1 , where the link between these indices is given by If the optional parameter marked is omitted, the uniformly refined mesh provides an accordingly refined mesh.The parameters marked dirichlet and marked neumann are optional and may again be omitted to obtain uniform refinement.

Discretization of the domain
In case of non-homogeneuous volume data, the evaluation of the Newtonian potential requires the discretization of .To this end, we restrict to a regular triangulation T of into compact non-degenerate triangles T ∈ T .For the data structure, we use a standard list format as is also done in [1,22].The set of vertices V = {v 1 , . . ., v m } of the triangulation T = {T 1 , . . ., T n } is represented by an m × 2 array vertices.The j -th row stores the coordinates of the j -th vertex v j = (x j , y j ) ∈ R 2 as vertices(j ,:) = [ x j y j ].
The p-th triangle T p = [z i , z j , z k ] with vertices z i , z j , z k ∈ V is stored as where the nodes are given in counterclockwise order, i.e., the parametrization of the element's boundary T p is mathematically positive.
HILBERT contains functions for the local mesh-refinement of the volume triangulation T and for the extraction of the induced boundary partition E = T | .Finally, by use of the same data structure as in [1,22], it is easily possible to link HILBERT to existing MATLAB finite element codes to realize, e.g., the (adaptive) coupling of FEM and BEM.

Boundedness of K-mesh constant
Many estimates in numerical analysis depend on local quantities of the mesh, e.g., on an upper bound of the K-mesh constant which is the maximal ratio of the element widths of neighboring elements.To avoid the blow-up of the K-mesh constant, the mesh-refinement algorithm implemented in refineBoundaryMesh guarantees by one bisection of all elements in a certain superset R ⊇ M .To that end, the function refineBoundaryMesh implements a mesh-refinement strategy analyzed in [3].
Note that #M ≤ #R = #E +1 − #E , whereas the converse estimate #E +1 − #E #M cannot hold under the additional constraint (15).However, as is proven in [3, Theorem 3], our mesh-refinement algorithm ensures that the overall number of refined elements does not exceed the overall number of marked elements arbitrarily.More precisely, refineBoundaryMesh ensures The constant hidden in the symbol depends only on the initial mesh E 0 .The validity of such an estimate is crucial for proving quasi-optimal convergence rates for residual-based ABEM in [19,36] as well as [3,17,18].

Discrete function spaces
Let P p (E ) be the space of all E -piecewise polynomials of degree p ∈ N 0 with respect to the arc-length.Note that functions ∈ P p (E ) are, in general, not continuous, but have jumps at the nodes of E .In particular, P 0 (E ) denotes the space of all E -piecewise constant functions.If χ j ∈ P 0 (E ) denotes the characteristic function of E j ∈ E , the set {χ 1 , . . ., χ N } is the basis of P 0 (E ) which is used throughout our implementation.
One particular example for a function in P 0 (E ) is the local mesh-width h ∈ P 0 (E ) which is defined E -elementwise by Let S 1 (E ) := P 1 (E ) ∩ C( ) denote the set of all continuous and (with respect to the arc-length) E -piecewise affine functions.For each node z j ∈ K of E , let ζ j ∈ S 1 (E ) be the nodal "hat function" associated with the node z j ∈ K , i.e., ζ j (z k ) = δ jk with Kronecker's delta.Then, the set {ζ 1 , . . ., ζ N } is the basis of S 1 (E ), which is used throughout our implementation.
In the following, we only consider the lowest-order BEM, where the spaces P 0 (E ) and S 1 (E ) are used to discretize fluxes and concentrations, respectively.

Discrete integral operators
The Calderón system (6) essentially relies on four linear integral operators V , K, W , and N 0 .The remaining two operators K and N 1 can then be expressed in terms of the other four, see (56) below for N 1 .Following our lowest-order ansatz, HILBERT provides a C-implementation of integral operators for discrete fluxes and concentrations: -the single-layer integral operator matrix V for P 0 (E ) ansatz and test functions is given by -the double-layer integral operator matrix K for S 1 (E ) ansatz and P 0 (E ) test functions is given by -the hypersingular integral operator matrix W for S 1 (E ) ansatz and test functions is given by -Given a regular triangulation T of the domain , the trace of the Newtonian potential N for P 0 (T ) ansatz and P 0 (E ) test functions reads Antiderivatives, i.e. explicit formulas, for the integrals of V and K are taken from [26].With similar techniques, we developed antiderivatives for the computation of N. The computation of W is implemented by use of Maue's formula [28] W u, v = V u , v , where (•) denotes the arc-length derivative.However, for adaptive meshes, the analytical computation of the integrals leads to instabilities due to cancellation effects when using the closed formulas for the (continuous) antiderivatives, e.g.b a f dt = f (b) − f (a) for a ≈ b.For two boundary elements E i and E j , we change the order of integration in case of diam(E j ) < diam(E i ) to ensure that the outer integration is performed over the the smaller boundary element.Given some fixed parameter η > 0, we call two boundary elements E i , E j admissible if they satisfy In this case, the outer integration is replaced by Gaussian quadrature.In particular, this guarantees that the matrices V and W related to the symmetric operators V and W are symmetric as well.For admissible elements E i , E j , there holds exponential convergence of the semi-analytic quadrature used with respect to the order of the Gaussian quadrature rule chosen, see [29] or [24] in the context of hierarchical matrices.
In HILBERT, the discrete integral operators are provided by the following MEXfunctions: The optional admissibility parameter eta may be omitted, a sane default value is used in that case.The assembly of the Newtonian potential N requires a volume mesh for the quadrature over the domain .To that end, we use the standard format as described in [1] and buildN uses mesh-manipulation routines as described in [22].The assembly of the discrete boundary integral operators is significantly time consuming.Therefore, HILBERT uses POSIX-threads, a simple parallelization library, to increase efficiency on multi-core systems.In our experiments, runtime scales properly with number of cores, see Fig. 2.

Evaluation of integral operators
Additionally to the discrete integral operator matrices, we provide routines for the evaluation of the integral operators V , K, K , W , and N introduced in Section 1.1.
The approximate evaluation of N 1 can be performed by postprocessing of the other evaluation functions, see (56) below.The operators V , K, N can be evaluated in the whole plane R 2 , whereas the operators K , W , and N 1 can be evaluated on .
In HILBERT, the evaluation routines of the integral operators are provided by the following MEX-functions: -V x = evaluateV(coordinates,elements,phih,x[,eta]) -K x = evaluateK(coordinates,elements,gh,x[,eta]) -N x = evaluateN(vertices,volumes,fh,x) Note that evaluateK returns the value of Kg h (x) for x ∈ R 2 \ and Kg h (x) for x ∈ .For an (M × 2) matrix x containing the evaluation points, these function calls return an (M × 1) vector of evaluations of the associated integral operators.As usual, the boundary mesh E is described by coordinates and elements, whereas vertices and volumes describe the volume mesh T .The optional admissibility parameter eta may be omitted, a sane default value is used in that case.The vectors phih, gh, and fh represent the data ∈ P 0 (E ), G ∈ S 1 (E ), and F ∈ P 0 (T ), respectively.The (M × 2) matrix nx describes the outer normal unit vector of , i.e., if x(k,:) corresponds to a point x ∈ , then nx(k,:) is n x .For the implementation, we write the boundary integrals as sum over the elements E j , for j = 1, . . ., N. For an evaluation point x and an element E j , the evaluation of the corresponding integral is carried out by quadrature in case of admissibility, diam(E j ) ≤ η distance(x, E j ), and by analytical formulas otherwise.The parameter η can be specified by providing eta to the functions.

Galerkin discretization
To discretize (22), we first replace the continuous Dirichlet data g ∈ H 1 ( ) by its nodal interpolant Alternatively, g can be discretized by means of the L 2 -projection onto S 1 (E ), i.e., with Second, we replace the function space H −1/2 ( ) in ( 22) by the finite-dimensional space X = P 0 (E ).Since the discrete space 23) is also a scalar product on P 0 (E ).Consequently, there is a unique Galerkin solution Let x ∈ R N denote the coefficient vector of the ansatz x j χ j (27) and let g ∈ R N be defined by g j := G (z j ) for all z j ∈ K .With the matrices V, K ∈ R N ×N defined in Section 3.3.2and the mass matrix M ∈ R N ×N defined by the Galerkin formulation ( 26) is equivalently stated in terms of the linear system We stress that V is symmetric and positive definite since it stems from a scalar product, i.e., V ij = χ i , χ j V and hence y • gy = , V > 0 with = j y j χ j = 0 for y ∈ R N \ {0}.

(h − h/2)-error estimation and refinement indicators
Instead of discretizing the correct variational form (22), in fact, we discretize with perturbed right-hand side, where we use the approximation G ≈ g.It is an analytical observation that the error between the exact solution φ ∈ H −1/2 ( ) of ( 22) and the exact solution φ ∈ H −1/2 ( ) of the perturbed formulation ( 29) is controlled by where (•) denotes the arc-length derivative, cf.[5,Section 4] for nodal interpolation and [23,Section 3] for the L 2 -projection, respectively.
As already mentioned, the non-locality of the integral operators leads to difficulties measuring the local contribution of a function to its norm.In [21,Theorem 3.4], we prove that Here, : L 2 ( ) → P 0 (E ) is the L 2 -orthogonal projection onto the space of piecewise constants, which is just the piecewise integral mean The estimator μ D, is stated in a weighted L 2 -norm and may thus be used to steer the local mesh-refinement.
If we plot η and μ D, over the number of elements, from the equivalence (34) of estimators, one can predict that the corresponding curves, for a sequence of arbitrarily refined meshes, are parallel, cf.[15,21].
The equivalence (34) as well as the error control (33) lead to the choice of ρ (E j ) 2 := μ D, (E j ) 2 + osc D, (E j ) 2 as error indicator to steer Algorithm 1.Based on results from [5,6], one can prove that this choice of ρ and the Dörfler marking criterion (7) guarantee lim ρ = 0, see [5,Theorem 5.4].Therefore, if the saturation assumption ( 12) holds (at least for the meshes generated), we obtain convergence of to φ.

Implementation of adaptive algorithm
The MATLAB script adaptiveSymmMuTilde.m which is contained in demos/laplace/ demonstrates the use of HILBERT.The adaptive loop merely contains eight command lines.The discretization of the boundary data is performed by the functions computeOscDirichlet and computeOscDirichletL2, which besides computing the data oscillations, also return the coefficient vector of the discretized data G from ( 24) or (25), respectively.The function buildSymmRHS computes the right-hand side vector of ( 28) which corresponds to (K + 1/2)G .
The linear system is solved with the MATLAB backslash-operator.The function computeEstSlpMuTilde provides the local quantities μ D, (E) 2 .The marking criterion ( 7) is implemented in markElements.

Numerical experiment
We perform Example 1 from Section 2.4.2 as benchmark and compare the results obtained by the proposed adaptive algorithm with those obtained from a uniform approach.The Dirichlet data g = u| is given, and the missing information φ = ∂ n u is computed by solving Symm's integral equation.Figure 3 shows errors, error estimators, and data oscillations with respect to the number of boundary elements, the computational time, and the memory consumption, where we use nodal interpolation (24) of the Dirichlet data.Since the Galerkin error in the energy norm can hardly be computed, we use the additional regularity φ ∈ L 2 ( ) and argue as in [5,Section 7] to prove the upper bound We stress that the inequalities hold without appealing to the saturation assumption (12) and may thus serve to check (12) numerically.The singularities of the analytical solution lead to a reduced order of convergence O #E −2/3 when using uniform meshes.The adaptive algorithm, however, resolves the singularities of the solution φ as well as of the given data g and reveals the optimal convergence behavior.The more interesting information for practitioners is effectivity with respect to computational resources.Figure 3 clearly shows that, in our experiment, the overhead introduced by the adaptive algorithm is soon overcome due to higher order of convergence.Using the L 2 -projection (25) to discretize the Dirichlet data, we observe the same results (not displayed).
As mentioned in the introduction, a particular strength of the boundary element method is that the representation formula (2) admits the definition of an approximate solution which, under certain assumptions, exhibits a high order pointwise convergence towards the solution u of (18), see e.g.[35,Section 12.1]:More precisely, for smooth solution φ = ∂ n u, lowest-order boundary elements ∈ P 0 (E ), and exact Dirichlet data g = G , theory allows for for uniform meshes E with mesh-size h .For nodal interpolation (24) of the Dirichlet data, this order is reduced to while L 2 -projection of the Dirichlet data (25) onto S 1 (E ) preserves third-order convergence, cf.[35].
Although, the current theory of adaptive BEM is tailored to the energy norm, it might be of interest how the adaptive algorithm behaves with respect to the pointwise error.To that end, we compare uniform and adaptive mesh-refinement for both, nodal interpolation (24) and L 2 -projection (25) of the Dirichlet data, in Fig. 4. We plot the energy errors as well as the point error for an example point x = (0.2, 0.13) ∈ .Whereas the discretization of the Dirichlet data does essentially not influence the convergence order of the energy errors, the influence on the point error is visible.Moreover, the proposed adaptive algorithm recovers the optimal orders even for the Fig. 3 Example 1 computed with given Dirichlet and unknown Neumann data, see Section 2.4 for the experimental setup.For uniform mesh-refinement, the singularity of φ leads to a reduced order of convergence O #E −2/3 , whereas the adaptive strategy recovers the optimal order of convergence . Moreover, the adaptive scheme is also superior with respect to computational time and memory consumption We assume that is connected, i.e. has no holes.Note that the PDE formulation (36) implies a side constraint on φ, namely The second equation of the Calderón system (6) yields the hypersingular integral equation with the hypersingular integral operator and K the adjoint double layer integral operator defined by Kφ, ψ = φ, K ψ .We stress that ( 38) is an equivalent formulation of (36).Recall that W : H 1/2+s ( ) → H −1/2+s ( ) and K : H −1/2+s ( ) → H −1/2+s ( ) are linear and continuous operators for all −1/2 ≤ s ≤ 1/2.Moreover, if the constant functions are factored out, W : ( ), where the subscript abbreviates the constraint φ , 1 = 0. We will, however, assume additional regularity φ ∈ L 2 ( ).The exact solution g ∈ H 1/2 ( ) of the integral formulation (38) is the Dirichlet data u| of the solution u ∈ H 1 ( ) of (36).Moreover, the additional regularity φ ∈ L 2 ( ) of the data provides additional regularity g ∈ H 1 ( ) of the trace.
Due to W c = 0 for all constant functions c ∈ R, the solutions of ( 36) and ( 38) are only unique up to additive constants.To fix the additive constant, we consider the bilinear form which leads to the following variational form: One can prove that • , • W +S from (40) defines a scalar product such that the induced norm ||| • ||| W +S is an equivalent norm on H 1/2 ( ).Consequently, (41) has a unique solution g which depends continuously on the Neumann data φ.Moreover, the solution automatically satisfies g ∈ H 1/2 * ( ), which follows from the constraint (37) on φ, and therefore also solves (37).

Galerkin discretization
To discretize (41), we first replace the Neumann data φ ∈ L 2 ( ) by its L 2 -projection ∈ P 0 (E ), According to this definition and φ ∈ H ( ).Second, we replace the function space H = H 1/2 ( ) in (41) by the finite-dimensional space X = S 1 (E ).Since S 1 (E ) is a subspace of H 1/2 ( ), there exists a unique Galerkin solution G ∈ S 1 (E ) of the discretized problem As in the continuous case, the discrete solution G automatically satisfies G d = 0, i.e., G ∈ S 1 (E ) ∩ H 1/2 * ( ).Recall the definition of p ∈ R N from (42) and let x ∈ R N denote the coefficient vector of the ansatz We recall the definition of the mass matrix M ∈ R N ×N from Section 4 as well as the definition of the matrices K, W ∈ R N ×N from Section 3.3.2.With the rank-1 stabilization matrix S∈ R N ×N defined by the Galerkin system (43) is equivalently stated in terms of the linear system Finally, we stress that the matrix W + S from ( 45) is symmetric and positive definite since it stems from a scalar product, i.e., (W + S) ij = ζ i , ζ j W +S and hence y• (W + S)y = V , V W +S > 0 with V = j y j ζ j = 0 for y ∈ R N \ {0}.

(h − h/2)-error estimation and refinement indicators
Instead of discretizing the correct variational form (41), we discretize with perturbed right-hand side, where we use the approximation ≈ φ.Analytically, the error between the exact solution g ∈ H 1/2 ( ) of (41) and the exact solution g ∈ H 1/2 ( ) of the perturbed formulation (46) is controlled by the data oscillations [5, Section 4], As for Symm's integral equation, the error estimator η N, = ||| G − G ||| W from ( 9) controls the discretization error by where the upper bound relies on the saturation assumption for the non-perturbed problem.This again leads to according to the triangle inequality and (47).
We stress that we are now dealing with different norms than in the case of the Dirichlet problem.In [16], it is proven that The nodal interpolation operator I : C( ) → S 1 (E ) is defined as in (24).We need to estimate both, the Galerkin error and the data-discretization error which leads to the choice of ρ (E j ) 2 := μ N, (E j ) 2 + osc N, (E j ) 2 as error indicator to steer Algorithm 1.As discussed in Section 4.4, and similar to [5, Section 5], we prove that our adaptive algorithm guarantees lim ρ = 0 for the hypersingular integral equation.Therefore, if the saturation assumption holds (at least for the meshes E generated), we obtain convergence of G to g.

Implementation of adaptive algorithm
The MATLAB script adaptiveHypMuTilde.m which is contained in demos/laplace/ demonstrates the use of HILBERT.The adaptive loop merely contains 9 command lines.The discretization of the boundary data is performed by the function computeOscNeumann which besides computing the data oscillations also returns the coefficient vector of the discretized data.The stabilization matrix S is computed by the function buildHypsingStabilization. The function buildHypsingRHS computes the right-hand side (1/2 − K ) in (45).The linear system is solved with the MATLAB backslash operator.The function computeEstHypsingMuTilde provides the local quantities μ N, (E) 2 .

Numerical experiment
We perform Example 1 from Section 2.4.2 as benchmark and compare the results obtained by the proposed adaptive algorithm with those from a uniform approach.The Neumann data φ = ∂ n u is given, and the missing information g = u| is computed by solving the hypersingular integral equation.Since the Galerkin error in the energy norm can hardly be computed, we use the additional regularity g ∈ H 1 ( ) and argue as in in the case of the Dirichlet problem to obtain the upper bound We stress that the inequalities hold without appealing to the saturation assumption (12) and may thus serve to check (12) numerically.Figure 5 shows error, error estimators, and data oscillations with respect to the number of boundary elements, the computational time, and the memory consumption.The singularities of the Cauchy data lead to a reduced order of convergence O #E −2/3 when using uniform meshes.Again, the adaptive For uniform mesh-refinement, the singularities of the Cauchy data lead to a reduced order of convergence O #E −2/3 , whereas the adaptive strategy recovers the optimal order of convergence O #E −3/2 .Moreover, the adaptive scheme is also superior with respect to computational time and memory consumption algorithm reveals the optimal convergence behavior.Moreover, the adaptive approach in our experiment is clearly superior to a uniform strategy in the sense that available computational resources are used more efficiently.

A mixed boundary value problem for the Poisson equation
In this section, we use the techniques introduced in Sections 4-5 to solve a mixed boundary value problem.Additionally, we consider non-vanishing volume forces: where N is assumed to be connected.For the equivalent integral formulation of (51), we choose (and fix) arbitrary extensions g D ∈ H the unknown data g N ∈ H 1/2 ( N ) and φ D ∈ H −1/2 ( D ) then satisfy the system of integral equations The definition which is equivalent to the usual product norm, i.e., the bilinear form on the right-hand side is continuous and elliptic, but now lacks symmetry.
In order to provide black-box schemes to solve (54), we do not only discretize the given boundary data, but also the volume forces.To that end, let T denote a triangulation of the domain with T | = E .We then replace N 0 f by N 0 F with F = π f ∈ P 0 (T ) the L 2 -orthogonal projection of f onto the space of piecewise constants over T .To compute an approximation of N 1 f , we use the well-known identity see [35,Lemma 6.20], which follows from the idempotency of the Calderón system (6) for f = 0.The discrete scheme now reads as follows: First, we solve Symm's integral equation to obtain an approximation ≈ V −1 N 0 f .In the next step, we seek a solution where The algorithm for solving a mixed problem with volume force may be sketched as follows: -Input: volume mesh T , boundary mesh -Solve the linear system For details, the reader is referred to [5, Section 6].
. For a posteriori error estimation of the Galerkin error and the data oscillations on the boundary, we may therefore use the estimators introduced above for the Dirichlet and Neumann problem.
To include the discretization error introduced by solving (58), we use the local error estimator is the solution of (57) with respect to the fine mesh E .To steer an adaptive mesh-refining algorithm, it is therefore natural to use the combined error estimator where osc 2 V , := h (F − f ) 2 L 2 denote the data oscillations of the volume forces.The estimator μ , (E) and osc , (E) are either μ D, (E), osc D, (E) or μ N, (E), osc N, (E) in the case of E ⊆ D or E ⊆ N , respectively.The complexity of the implementation of the adaptive algorithm for the mixed boundary value problem is significantly higher than in the simpler cases discussed earlier.
The example file adaptiveMixedVolMuTilde.m, which is contained in the Fig. 6 Example 2 computed with mixed boundary data, see Section 2.4 for the experimental setup.For uniform mesh-refinement, the singularities of φ and N 1 f lead to reduced order of convergence O #E −4/7 , whereas the adaptive strategy recovers the optimal order of convergence O #E −3/2 .Even though the volume data f has a weak singularity, a reduction of order of convergence due to bad resolution of f in the uniform case is not observed demos/poisson/ folder, demonstrates how the adaptive scheme can be implemented with HILBERT.The adaptive loop still only consists of 33 command lines.

Numerical experiment
We compute the numerical solution of Example 2 from Section 2.4.3.The results of the numerical experiment are visualized in Fig. 6.For the refinement of the volume mesh we use newest-vertex bisection where marked triangles are bisected, see e.g.[37,Chapter 5].Further refinements are performed to ensure regularity of the volume mesh and the constraint T | = E .The last condition has turned out to increase stability of the integral operators as implemented in HILBERT.
In case of uniform mesh-refinements, we simply perform three bisections per triangle, i.e. all edges of T are halved.
Recall that the volume forces have a weak singularity at (0.14, 0.14).However, the data oscillations of f seem to decay fast enough in the sense that the order of convergence with respect to the number of boundary elements is not limited by resolution of volume data in practice.However, strong generic singularities of φ and N 1 f at the reentrant corner limit the convergence behavior to approximately O #E −4/7 in the case of uniform mesh-refinement.
Our adaptive algorithm now additionally measures the error due to the approximation of N 1 f and data oscillations of f over T .It recovers the mathematically optimal order of convergence O #E −3/2 .Moreover, as in the other case studies, the accuracy of the computation is increased with respect to computational time.Since Example 2 includes the computation of the Newtonian potential N 0 f , the figure showing the error and estimators with respect to memory consumption now also takes the storage requirements of N into account.It is remarkable that the adaptive approach is significantly superior to a uniform approach even then.

Features
The MATLAB program package HILBERT provides stable implementations of the discrete integral operators corresponding to the 2D Laplacian as well as many other functions necessary for an easily accessible implementation of adaptive BEM.In this paper, only a short presentation of the library is possible.However, HILBERT is distributed with a full documentation [2] and many example benchmarks.Some of these have gently been proposed at conferences by colleagues in order to stresstest our proposed adaptive approach under difficult conditions.Besides the features presented here, HILBERT comes along with many visualization tools, further error estimators, and different marking strategies.
The HILBERT program package includes several example files like e.g.adaptiveSymmMuTilde, adaptiveHypMuTilde, and adaptiveMixed VolMuTilde, to illustrate that adaptive BEM may be easily implemented by use of HILBERT.This makes the tool not only interesting for scientists, but also for lecturers planning classes on BEM or scientific computing.The software is under constant development, and updates are released continuously.In the near future, functionality for the numerical solution of transmission problems by use of adaptive FEM-BEM coupling as well as linear elasticity will be included (visit the project webpage http://www.asc.tuwien.ac.at/abem/hilbert for detailed information on the ongoing development).
The recent progress in the analytical understanding of adaptivity in general and in the context of BEM specifically, allows to implement mathematically justified adaptive algorithms, which automatically resolve the singularities of both, analytical solution and given data.We observe that our proposed algorithm -based on easyto-implement (h − h/2)-type error estimators -empirically succeeds to recover the optimal order of convergence in all benchmarks performed so far.One of the aims of the authors is to emphasize the use of adaptivity and to make the concept more accessible to practitioners.

Restrictions
HILBERT is academic code in the sense that a MATLAB implementation generically might be too slow for use in industrial applications.Moreover, HILBERT currently only provides implementations of the integral operators associated with the 2D Laplacian.It is restricted to lowest-order elements and the canonical basis functions.
The most important restriction for practitioners might be that resolution of geometry is, so far, not included in the error estimation.In particular, all analytical results as well as the implementation demand to be a polygonal domain, and boundary elements are chosen to be affine line segments.

Fig. 1
Fig. 1 Rotated L-shaped domain .The initial boundary mesh E 0 consists of 8 boundary elements.In case of non-trivial volume forces, the initial triangulation T 0 (dashed) of consists of 6 triangles.In Example 1, the solution is prescribed such that the Neumann data φ = ∂ n u has a generic singularity at the reentrant corner, the Dirichlet data g = u| has a singularity at the uppermost corner.For the mixed boundary value problem covered by Example 2, the Dirichlet boundary consists of the two boundary elements that share the reentrant corner as common point (indicated in blue).The Neumann boundary is the remaining part of (indicated in red).The corresponding data structure is presented to the right of the picture Nf of the Newtonian potential as well as the Neumann data φ = ∂ n u have a singularity at the reentrant corner.We aim at solving some mixed boundary value problem.To that end, we split the boundary into Dirichlet and Neumann part as described and shown in Fig. 1. 3 Implementation of mesh-refinement and integral operators 3.1 Discretization of the boundary Throughout, = ∂ is the piecewise affine boundary of a polygonal Lipschitz domain ⊂ R 2 .If necessary, is partitioned into finitely many relatively open and disjoint boundary pieces, e.g. in a Dirichlet boundary D and a Neumann boundary N , i.e., = D ∪ N and D ∩ N = ∅.

Fig. 2
Fig.2Speed up of operators using parallelization, see Section 2.4 for the experimental setup.The computational time is plotted in a double logarithmic scale over the number of cores.Computational were performed for a uniform mesh with 4.096 boundary elements and a volume triangulation with 24.576 triangles.For reference, the function 100/k with k = 1, 2, . . ., 6 the number of cores, is plotted

Fig. 4 5 problem 5 . 1
Fig.4 Example 1 for uniform and adaptive mesh-refinement, where the given Dirichlet data are discretized by either nodal interpolation or the L 2 -projection onto S 1 (E ).Besides the error estimator η D, , we plot the point error err , := |U (x) − u(x)| between the exact solution of (18) and its approximation (35) at the example point x = (0.2, 0.13) ∈ .For uniform mesh-refinement, the point error behaves like O #E −4/3 , i.e. twice the order of the energy error.Even for the point error, adaptivity recovers the theoretically optimal orders O #E −2 and O #E −3 for nodal interpolation resp.L 2 -projection .∈ P 0 (E ) ∩ H −1/2 *

Fig. 5
Fig.5 Example 1 computed with given Neumann and unknown Dirichlet data, see Section 2.4 for the experimental setup.For uniform mesh-refinement, the singularities of the Cauchy data lead to a reduced order of convergence O #E −2/3 , whereas the adaptive strategy recovers the optimal order of conver- +1 which is only refined locally in the sense that (at least) all elements of M are refined.A marked element E i ∈ E is bisected to certain sons e j , e k ∈ E +1 .The (N × 2)-matrix father2son provides a link between the element indices in the sense that 1/2 ( ) and φ N ∈ H −1/2 ( ) of the given data from D and N , resp., to the entire boundary .The missing boundary data, which have to be computed, are g N := u| − g D and φ D := ∂ n u − φ N .