Local Defect Correction for Boundary Integral Equation Methods

G.Kakuba and M.J.H.Anthonissen

1 Introduction

Often boundary value problems have small localised regions of high activity where the solution varies very rapidly compared to the rest of the domain.This behaviour maybe either due to boundary conditions or due to an irregular boundary.One therefore has to use relatively fine meshes to capture the high activity.Since the activity is localised,one may also choose to solve on a uniform structured grid.That is,instead of a uniform global grid,the solution is approximated using several uniform grids with different grid sizes that cover different parts of the domain.The size of each grid is chosen in agreement with the activity of the solution in that part of the domain.This refinement strategy is calledlocal uniform grid refinement[Ferket and Reusken(1996)].The solution is approximated on a composite grid which is the union of the various uniform local grids.One way of approximating this composite grid solution that is simple and less complex is byLocal Defect Correction(LDC).

In LDC,at least one grid,theglobal coarse grid,covers the entire domain.Then a uniformlocal fine gridis used in a small part of the domain containing the high activity.In[Ferket and Reusken(1996);Hackbusch(1984)],LDC has been shown to be a useful way of approximating the composite grid solution in which a global coarse grid solution is improved by a local fine grid solution through a process whereby the right hand side of the global coarse grid problem iscorrectedby thedefectof a local fine grid approximation.This method has been well explored for other numerical methods such as finite differences and finite volumes,see[Ferket and Reusken(1996);Hackbusch(1984);Anthonissen(2001);Minero,Anthonissen,and Mattheij(2006)].In this paper we explore potential analogues and develop an LDC strategy for boundary integral equation(BIE)methods,in particular,the boundary element method(BEM).The BEM is now a well established technique for potential problems as it leads to a reduction in the dimensionality of the problem,due to the need to discretise only the domain boundary.Moreover it provides accurate solutions due to the use of fundamental solutions and,in electrostatics,has many numerical advantages over other methods as finite differences methods as discussed in Liang and Subramaniam(1997).However,although reduction in dimensionality is guaranteed,the accuracy depends a lot on the mesh used[Guiggiani(1990);Liapis(1994)].Although a lot has been done for the finite element methods,BEM’s main competitor,comparatively little has been done on adaptive mesh refinement in BEM[Carsten and Stephan(1995,1996)].This paper therefore seeks to add to a gradually growing literature on mesh refinement in BEM.

The initial attempts on LDC for BEM were in[Kakuba,Mattheij,and Anthonissen(2006)]and[Kakuba and Mattheij(2007)]where an algebraic approach was suggested and studied.The algorithm in[Kakuba,Mattheij,and Anthonissen(2006)]was based on block decomposition and manipulation of the coefficient matrix and right hand side of the BEM equations in three dimension space.That formulation ignores the inherent nature of the Boundary Integral Equation(BIE)that it is glob-al,that is,each node of the grid contributes to all the others in the grid through integration.In this paper we consider a new and better approach to the LDC algorithm that is based on the global integration.This approach takes into account the global integral properties of BEM in the computation of the defects and in the formulation of the local problem.Since in BEM we discretise the boundary,we are concerned with problems in which the high activity occurs at the boundary.

One ofthe most important steps in adaptive refinementis error estimation.One way of estimating the error in collocation BEM is by using higher order interpolation at the elements to estimate the exact solutions and then use an appropriate norm of the difference between the BEM solution and the estimated exact solution to measure the error[Kita and Kamiya(2001)].Since the focus of this paper is to develop the essential steps of the LDC algorithm in BEM,we use simple examples whose exact solutions are known.We also consider high activity due to boundary conditions.Then we use the infinity norm||·||∞of the difference between the exact and the BEM solutions as the measure for the error.

The paper is organised as follows:First,we give a brief introduction to the BIE and its discretisation using BEMin Section 2.In Section 3,we develop an LDCstrategy for BEM alongside an example.In Section 4 we test the strategy on a typical example and discuss results.In Section 5 we give the advantages of the algorithm and finish in Section 6 with a summary of the concepts and results presented in the paper.

2 The Boundary Integral Equation and the Boundary Element Method discretisation

The BEM results from a numerical discretisation of a BIE.Consider a domain ?with boundary?? on which we have the following Laplace problem:

where??1∪??2=?? and??1∩??2=/0 andgandhare given functions,see Fig.1.If??1≡?? we have a Dirichlet problem and if??2≡?? we have a Neumann problem,otherwise we have a mixed problem.The boundary integral equation gives relations for the potentialu(s)at different locations of the domain ?.These relations have been abundantly derived in literature and are readily available in various books on boundary element methods such as[Paris and Canas(1997),Katsikadelis(2002),Pozrikidis(1992)].The BIE is

Figure 1:Domain illustration for a Dirichlet or Neumann problem

α(s)is the internal angle ats,?c=R2? is the complement of ? andv(s;r)is the fundamental solution of the Laplace equation.

At the boundary,the discretised BIE leads to the linear system of equations

We also have introduced the vectors

The solution of the system(8)gives a BEM approximation of the unknowns in x in the grid nodesri.We denote by xLa BEM approximation on a grid of sizeL.ThusuLj(orqLj)is a BEM approximation ofuj(orqj)using a grid of sizeL.Solving(8)gives the unknown boundary quantities ofuandq.Therefore we now have all the boundary quantities.The solutionuiat any pointri∈? can then be computed using

3 Local Defect Correction

Consider the Neumann problem(10)whose domain is a unit square in two dimension.That is,

The Neumann problem(10)results in a Fredholm integral equation of the second kind[Atkinson(1997);Hackbusch(1995);Kress(1989)].It is singular and hence has no unique solution[Atkinson(1967)].To ensure a unique solution in the discretised problem,a value ofuis prescribed in one of the nodes[Chen and Zhou(1992);Strese(1984)].In the implementations in this paper,the value ofuin thefirst element node is prescribed.This will also help us compare the numerical solutions directly with the exact solutions for error measurement.The solution in ?,shown in Fig.2,has a small area close to the boundary where it changes rapidly.As a result,the solutionu(r)in the boundary has a region of high activity in a small part of the boundary,see Fig.2.Therefore we can identify a small region inside ? which contains the high activity.This region we call thelocal domainand denote it by ?local,see Fig.3.Its boundary Γlocal,which we will call thelocal boundary,consists of two parts:a part Γactivethat is also part of the global boundary and a part Γinsidethat is contained in the global domain ?,Fig.3(b).We will call the part Γactivethelocal active boundary.For instance,in the problem corresponding to the solutions shown in Fig.2,the boundary Γactive,may be identified as Γactive={(x,y):y=0,x∈[0.2,0.8]}.The part of the global boundary Γ that is outside the active region Γactivewill be denoted Γc,that is,Γc:= ΓΓactive.

The interest in BEM is to compute a numerical approximation ofu(r)at the boundary as accurately as possible.For such kind of multiscaled variations one is faced with the option of using a global uniform grid with a mesh of relatively small sizelin order to capture the high activity.This would result in very large systems which are computationally expensive since BEM matrices are full matrices.Besides,outside the local active boundary Γactive,the variation of the solution is smooth and a relatively coarse grid would suffice.The other option is to use a uniform structured grid designed to capture the different activities.This would be a composite grid with a relatively fine mesh of sizelin the local active region and a coarse grid of sizeLelsewhere.

Figure 2:Solution in the domain and solution at the boundary 0≤x≤1,y=0,for problem(10).

With LDC we approximate the solution on a composite grid in an iterative way that involves solving a so calledlocal problemwhich is a boundary value problem defined on the local domain.The local problem is solved on a fine mesh whose size is chosen in agreement with the local activity.The solution on the local fine grid is combined with the solution on the global coarse grid throughdefect correctionto obtain a composite grid solution on Γ.The advantage of this approach is that instead of solving a large composite grid system,two smaller systems;a global coarse grid system and a local fine grid system,are solved independently.For problemswith various local activities the local problems can be solved separately in parallel giving a tremendously cheaper way of obtaining a composite grid solution other than solving directly on the composite grid.

Let Γ be the numerical representation of?? in BEM.Theglobal coarse gridΓLis a uniform mesh ofNelements each of sizeLcovering the whole of Γ,that is,

Figure 3:An example of a multiscaled solution with localised high activity in 3(a)and,in 3(b),an illustration of a local problem domain.The boundary of ?local is

where|ΓLj|=Lfor allj.Thelocal fine gridΓllocalis a uniform mesh ofNlocalelements each of sizelcovering Γlocal,that is,

where|Γllocal,i|=lfor alli=1,2,...,Nlocal.The size of the local fine gridlis chosen in agreement with the activity of the solution in Γactive.Since the solution varies much more rapidly in Γactivethan elsewhere,we expectlto be much smaller thanL.Part of the grid Γllocalbelongs to Γactiveand part belongs to Γinside.The part that belongs to Γactiveis denoted Γlactiveand that that belongs to Γinsideis denoted

Figure 4:Global coarse and local fine grids.The small dots are the nodes rl localof the local fine grid Γl local and the big circles are the nodes rL of the global coarse grid ΓL.Node 2 belongs to rL∩rl active.

First we discretise the BIE on ΓLto yield

which gives the initial global coarse grid system of equations

Once we have solved(18),the next step is to use the initial global coarse grid solution uL0to formulate a local problem on ?local.This local problem on ?localsatisfies the same operator as in the global problem.The boundary conditions on Γactiveare the same as those in the global problem,that is,q(r)=h(r),since Γactive?Γ.On Γinsidewe prescribe an artificial boundary condition?g(r)defined below.So we have

The functiong?(r)is piecewise constant on Γinsideand takes on values ofuinside(ri)whereriis a node of Γlinside,i,that is,

Equation(21)means that we use the solution of the initial global coarse grid problem to obtain artificial Dirichlet boundary conditions at Γinside.Since at Γactiveqis known,the local problem is mixed and the BIE for(19)is,forr,r(χ)∈Γlocal,

Discretising(22)on the local fine grid defined in(13)and(14)we have

The solution ul0activeis expected to be more accurate than the coarse grid solution uL0in Γactive.The next step of LDC is to use the local fine grid solution to update the global coarse grid problem.In updating,the right hand side of the global coarse grid problem is corrected by the defect of the local fine grid approximation,we will call this step thedefect correctionstep.The two approximations are then used to define a composite grid approximation ofu(r).The question now is:how do we compute the defect?

Consider the coarse grid discretisation(17).If we knew the exact continuous functionu(r)and hence the exact solutionuj:=u(rj)in the nodes we would use it in(17)to obtain

Subtracting(27)from(26)gives

Therefore if we know the exact continuous functionu(r)we can compute the local defectdiL,add it to the right hand side of(17)and solve for the exact solutionujon each element.Howeveru(r)is not known and therefore we cannot compute the defect using(29).All we can do is estimatediLjas accurately as possible using the best solution available,which is

So for elements in the high activity region we have the fine grid solution which we can use to estimate the local defect as follows.

Figure 5:A coarse element that is refined into three elements in the local fine grid

Therefore,using the initial fine grid solution,we have the following best approximation of the initial defect per element

By default,integration in the BIE is global.Each node of the global coarse grid communicateswith the active region through integration.So although the activity is local,it affetcs the entire system of equations.The defectdL0iis therefore computed for all nodes of the global coarse grid to generate the local defect vector

The next step now is the updating step.The global coarse grid discretisation is updated with the defect of the local fine grid solution.So we have

Solving(36)gives the updated coarse grid solution uL1.At this stage we use the fine grid solution on Γlactiveand the global coarse grid solution to form a composite grid solutionul,Las

The composite grid solution(37)can now be used to compute better boundary conditions on Γinsideand then form and solve the updated fine grid problem

Then we obtain the updated composite grid solution given by

This step marks the end of one complete cycle of the LDC algorithm.The iteration process is summarised in the following algorithm:

BEM-LDC Algorithm

(i)Initialisation

–Solve the global coarse grid system

–Solve the initial fine grid system

–Compute the initial defect

(ii)Fori=1,2,...

4 Examples and results

The LDC procedure above has been used to solve the problem(10).The results are shown in Fig.6 and Fig.7.The figures show the solution on the sidey=0 of the unit square.The coarse grid used is of sizeL=0.2 and the fine grid is of

Figure 6:Results of a typical LDC process for a Neumann problem in one iteration.The solid line is the exact solution.

sizel=0.2/9.Fig.6 shows how the initial results compare with the exact solution(the solid line),the initial coarse grid solution in Fig.6(a)and the initial fine grid solution in Fig.6(b).Fig.7 shows the results after the first update,the updated coarse grid solution in Fig.7(a)is better than the initial one in Fig.6(a)and the updated local fine grid solution is better than the initial one in Fig.6(b).Fig.8 shows how fast the global error decreases.Basically the algorithm has converged already in the first iteration since the error reduction between successive iterations after the first one is much smaller compared to that in the first iteration.

5 Advantages of LDC of the BEM-LDC algorithm

In brief,the LDC iterative process involves the following:solve the global coarse grid problem on Γ,computeuon Γinside,solve a fine grid problem on Γlocaland then update the global coarse grid problem.

Suppose we haveplocally active small regions and thusplocal problems.Let,for each local problem,Mlbe the number of elements on ΓlocalandMinthe number of elements on Γinside.For instance in the illustration in Fig.4,Ml=3 andMin=9.We can increaseMlwithout necessarily increasingMinsince the activity is on Γlocal.LetMinbe so small compared toMlthat the size of the local problem isM≈Ml.LetNbe the size of the global problem andNlLocalthe number of global elements in Γlocal.We assume Γlocalis such a small part of the global boundary thatN?NlLocal≈N.Then the equivalent composite grid on Γ would be of sizepM+N.The operational count for LU-decomposition isN3/3 for a sizeNmatrix.So the complexity of the equivalent composite grid problem would be

Figure 7:Results of a typical LDC process for a Neumann problem in one iteration.The solid line is the exact solution.

Figure 8:Graph of the global coarse grid error||u??uLi||∞,i=0,1,2,...,12 where u?is the exact solution.A logarithmic scale is used on the error axis.

The BEM-LDC algorithm converges in one step which involves solving two coarse grid problems andplocal problems and so has total complexity

6 Conclusions

We have presented the technique of Local Defect Correction for BEM,a technique for solving problems with high local activity in the boundary using BEM.The focus of the paper has been to develop the LDC algorith for BEM that takes into account the global integral nature of the BIE,a property that was ignored in Kakuba,Mattheij,and Anthonissen(2006).This technique offers an alternative to solving directly on a composite grid or a uniform fine grid both of which would result in large matrices that are more expensive than using LDC.What is also interesting to note is that one iteration of the algorithm suffices.

Anthonissen,M.(2001):Local Defect Correction Techniques:Analysis and Ap-plication to Combustion.PhD thesis,Eindhoven University of Technology,Eindhoven,The Netherlands,2001.

Atkinson,K.E.(1967): The solution on non-unique linear integral equations.Numerische Mathematik,vol.10,pp.117–124.

Atkinson,K.E.(1997):The Numerical Solution of Integral Equations of the Second Kind.Cambridge University Press.

Carsten,C.;Stephan,E.P.(1995): A posteriori error estimates for boundary element methods.Mathematics of Computation.,vol.64,pp.483–500.

Carsten,C.;Stephan,E.P.(1996): Adaptive boundary element methods for some first kind integral equations.SIAM Journal of Numerical Analysis.,vol.33,pp.2166–2183.

Chen,G.;Zhou,J.(1992):Boundary Element Methods.Academic press,London,San Diego.

Ferket,P.J.J.;Reusken,A.A.(1996): Further Analysis of the Local Defect Correction Method.Computing,vol.56,pp.117–139.

Guiggiani,M.(1990):Errorindicators foradaptive mesh refinementin the boundary element method-a new approach.International Journal for Numerical Methods in Engineering.,vol.29,pp.1247–1269.

Hackbusch,W.(1984):Local defect correction method and domain decomposition technique.Computing[Suppl.],vol.5,pp.89–113.

Hackbusch,W.(1995):Integral Equations:Theory and Numerical Treatment.Birkhauser Verlag,Basel.

Kakuba,G.;Mattheij,R.M.M.(2007):Convergence analysis of Local Defect Correction for the BEM.Advances in Boundary Integral Methods.Proceedings of the sixth UK conference on Boundary Integral Methods,vol.15,pp.251–261.

Kakuba,G.;Mattheij,R.M.M.;Anthonissen,M.J.(2006): Local Defect Correction for the Boundary Element Method.Computer Modeling in Engineering Sciences,vol.15,pp.127–135.

Katsikadelis,J.T.(2002):Boundary Elements theory and applications.Elsevier.

Kita,E.;Kamiya,N.(2001): Error estimation and adaptive mesh refinement in boundary element method,an overview.Engineering Analysis with Boundary Elements,vol.25,pp.479–495.

Kress,R.(1989):Linear Integral Equations.Springer-Verlag.

Liang,J.;Subramaniam,S.(1997): Computation of molecular electrostatics with boundary element methods.Biophysical Journal,vol.13,pp.1830–1841.

Liapis,S.(1994): A review of error estimation and adaptivity in the boundary element method.Engineering analysis with boundary elements.,vol.14,pp.315–323.

Minero,R.;Anthonissen,M.J.H.;Mattheij,R.M.M.(2006): A Local Defect Correction Technique for Time-Dependent Problems.Numerical Methods for Partial Differential Equations,vol.22,pp.128–144.

Paris,F.;Canas,J.(1997):Boundary Element Method:Fundamentals and Applications.Oxford University Press,Oxford.

Pozrikidis,C.(1992):Boundary integral and singularity methods for linearised viscous flow.Cambridge University Press.

Strese,H.(1984):Remarks concerning the boundary element method in potential theory.Appl.Math.Modelling,vol.8,pp.40–44.