Error Estimator Using Higher Order FEM for an Interface Problem

A higher order finite element method is considered to treat an interface problem. The polynomial degree is allowed to be arbitrary but it is fixed for the FEM-variational formulation. We propose an error estimator which turns out to be efficient and reliable. We demonstrate upper and lower bounds of the error estimator with respect to the exact accuracy. For the transmission problem, the coefficients for the internal and external domains can be highly dissimilar. One major difficulty is the characteristic of the estimator at the interface. The a-posteriori error estimates can be computed very efficiently element by element. To corroborate the theoretical analysis, we report on a few numerical results.


Introduction
In most aspects of numerical simulation, it is desirable to provide an approximation to an unknown.But it is even more reliable to supply an additional information quantifying how accurate the approximation is.In particular for FEM (Finite Element Method), that additional value is exactly the purpose of the a-posteriori error estimator which is described subsequently.We begin by motivating the transmission problem that is based on the PBE (Poisson-Boltzmann Equation) for the interaction of solute and solvent media which are respectively denoted by int Ω and ext Ω .The surface Γ represents the solute-solvent [1] [2] [3] interface which is the molecular surface in the realistic Applied Mathematics solute is located inside the cavity int Ω .In the sequel, the whole solute-solvent domain is denoted by :=  int ext Ω Ω Ω .The nonlinear PBE [4] admits the next general expression x Ω (1)   The charge positions i x are located in the strict interior of the molecular surface Γ as illustrated in Figure 1(a) and Figure 1(b).The quantities , , , , , , unknown function is the dimensionless electrostatic potential u.The coefficients ( ) ε x and ( ) κ x are space-dependent functions related to the dielectric value and the modified Debye-Hückle parameter.Those coefficients might be discontinuous between int Ω and ext Ω but the solution u is required to be continuous everywhere.In the case 2 M = and under certain assumption on the parameters, the exponential term becomes a hyperbolic sine.By the Taylor expansion, one has ( ) ( ) ( ) x Ω which will be the purpose of this paper.We will focus on the FEM treatment of the linearized equation for which we investigate a-posteriori error estimates.Before presenting our approach, related works and our previous results are in order.Verfürth has compiled a comprehensive study [5] about APEE (a-posteriori error estimator) for which he mainly treats piecewise linear FEM.Many different a-posteriori error estimators have been proposed for the Stokes problem [6] [7] [8] [9] for isotropic grids.In the context of anisotropic meshes [10] [11] [12], there are a variety of APEEs for Poisson and reaction-diffusion problems [13] [14].An article [15] by Creusé, Kunert and Nicaise presents a survey on the residual based error estimator on anisotropic grids for the Stokes equation.An interesting APEE for two and three dimensions as well as an anisotropic adaptive mesh refinement are also detailed in [16].Basically, a-posteriori error estimators permit to evaluate the finite element errors without knowing the exact solution.That feature makes it possible to dynamically identify regions where one should have further refinements if the error there is too large.Therefore, adaptive refinements are mainly based on the quality of a-posteriori error estimators.Our approach in this paper follows the same spirit as the works in [17] [18].For the Spectral Element Method, we find in [18] an APEE for the hp-case.That is, the mesh size h is allowed to be refined in some regions while the polynomial degrees p are also variable on different elements of the mesh.The hp-case does not require that the polynomial degree or the mesh size are fixed.That has been generalized in [17] to treat hp-FEM [19] for the Poisson problem where corner singularities are allowed.We have presented in [20] an APEE based on hierarchical space enrichment on anisotropic FEM which is combined with adaptive refinements.Boundary Element Method (BEM) is very efficient [21] [22] [23] [24] [25] as far as the linear PBE is concerned because of the existence of a fundamental solution providing an explicit kernel for the integral equation formulation.In addition, BEM requires only a discretization of the surface Γ instead of the entire volume Ω (see Figure 1(c)).When treating nonlinear PBE, a solver on the volume Ω appears to be unavoidable.This paper can be viewed as a preliminary work toward nonlinear PBE.We are still reluctant to completely focus on the nonlinear PDE because the equation in (1) presents a real challenge related to the exponential nonlinear term on the left hand side and the nonsmooth Dirac functions on the right hand side.The only treatment of nonlinear PBE using BEM which we are aware of is [26] that is admittedly a very good approach.By inspecting that paper in detail, we realized that a solver in the volume Ω is also needed to construct an artificial fundamental solution.An integral equation solved by BEM is then used by applying that artificial fundamental solution in order to form a kernel.That is, a treatment exclusively on the surface Γ without recourse to a solver in Ω is so far not sufficient.Holst [27] [28] [29] [30] is one of the most prominent specialists of PBE using FEM.His work seems to be extensively based on piecewise linear variational formulation.The Finite Difference Method (FDM) is also widely used in PBE.The main reason does not seem to be attributed to its numerical efficiency but rather to code availability and to reference or comparison purpose (see Section 1.1.2 of [31]).An important component of PBE simulation is the geometric information because exact solutions of PBE are only known for very few simple geometries.Implementing a program for generating an SES (Solvent Excluded Surface) from nuclei coordinates and the Van-der-Waals radii of the atoms is not straightforward because a lot of geometric tasks [32] [33] [34] come into play.It is a long process to start from the nuclei coordinates till obtaining geometric data for computations.We have achieved a geometric task to process nuclei information in order to generate data for BEM as well as a mesh generation [35] from FEM as illustrated in Figure 1(c).Furthermore, a real chemical simulation by using wavelet BEM is described in [25] for the quantum computation.A wavelet BEM simulation using domain decomposition techniques was described in [36] which treats the case of ASM (Additive Schwarz Method).It was utilized as an efficient preconditioner for the wavelet single layer potential which is badly conditioned.A wavelet BEM formulation for computing apparent surface charge is documented in [37] for an interface problem.A simulation for chemical quantum computation using FEM is documented in [38] where we used nanotubes immersed in polymer matrices.
We consider in this work a higher order FEM formulation to treat the interface problem.That is, the polynomial degree, which is arbitrary but fixed, is uniformly constant in the entire FEM mesh.We examine the a-posteriori error estimator locally within each element and each edge.There are several types of edges: the interface edges, the interior edges, the exterior edges and the boundary edges.The difficulty for an estimator with respect to an interface edge is the discontinuous coefficients on the incident elements as well as the flux jump at the interface.In this work, we are more interested in elaborating mathematical theory than in chemical simulation.The error estimator can be efficiently computed element by element.We consider smooth load functions as right hand side of the equation.In addition to the theoretical investigation, we contribute about the numerical influence of the parameters appearing in the a-posteriori error estimator.We need numerical tests because the dependence of all various constants with respect to the problem parameters is not established theoretically.The next discussion is structured as follows.In the following section, we start by formulating the problem accurately and we introduce some important definitions as well as the expression of the estimator.That is followed by the analysis of the a-posteriori error estimator in Section 1.We report on some interesting practical results in the last section.

Problem Setting and A-Posteriori Estimators
This section describes the problem formulation, the introduction of the higher order FEM as well as the explicit expressions of the a-posteriori error estimators.
We recall also some important results related to polynomial inverse estimates.
We consider the transmission PDE: ( ) ( ) x Ω x n x x n x x x Γ (5) where ( ) T T ∈  is either empty or a common node or a complete edge,  We have the coverings:  Every node of the interface Γ is also a node of h  ,  All edges of the interface Γ are edges of the mesh h  .

For a triangle
( ) : supremum of the diameters of all balls contained in : set of elements of sharing a vertex with We use ( ) for sufficiently large i.We assume quasi-uniformity in the sense that there exists a constant 0 0 We define the following mutually disjoint subsets of edges : set of edges of on the interface : set of edges of on the boundary , : set of edges of which are not included in , : set of edges of which are included neither in nor in .
Note that an edge of h  int may have an endpoint in Γ .Likewise, an edge of h  ext is allowed to have an endpoint in Γ or ∂Ω .We introduce in addition the set of all edges For an edge For any triangle T, the affine invertible mapping from the unit reference onto T is denoted by : : where det , det That allows one to derive results on the unit reference triangle T and to carry them over to any element T in the original mesh h  .We use the standard definitions of the Sobolev spaces for ( )  Ω .From here onward, we use the usual shorthand X Y  if there is a constant c such that

X cY ≤
in which c is independent on h and p.In addition, X Y  amounts to X Y X   .
We want to consider now the Galerkin variational formulation.Denote the restriction of the solution u to the interface problem in int Ω and ext Ω by u int and u ext respectively.Due to the Green identity we have We will denote the piecewise constant function defined on Ω : ( ) ( ) The sum of ( 24) and (25) for every ( ) Taking into account the flux jump (5) in the transmission equation, we have The Galerkin weak form is therefore For a fixed polynomial degree 1 p ≥ , the finite dimensional space is { } : span : 0 The discrete Galerkin approximation is to search for ( ) Introduce the bilinear form ( ) In order to express the a-posteriori estimators, we assume that the approximated solution h u on the current discretization For a general edge e and triangle T, transformations from the reference elements are used to define e ω and T ω .For an interior element int h T  ∈ , the estimator is defined as : where T f designates the ( )

T 
-projection of the load function f onto the element T. The expression for an exterior element : For an interface edge where e Q is the ( ) where   t stands for the jump of t when evaluated from the two elements incident upon e.The orientation of the jump is irrelevant because one takes its square in the 2  -norm.The estimator corresponding to an exterior edge Since one needs computable local estimators for an element-by-element computation, an interior element T T e e e T e e T e Likewise, for an exterior element h T ∈  ext , the local estimator reads: T T e e e T e e T e The local estimators add up to the global estimator: : One has the following polynomial inverse estimates and extension properties.
The descriptions of the next lemmas are found in [17] On the 2D reference element (22), one has for every bivariate polynomial p The constants are ( ) , C C δ = which do not depend on p. Lemma 2. Consider a univariate polynomial π of degree 1 p ≥ defined on the unit reference interval [ ] ˆ0,1 e = and a parameter 0 < 1 γ ≤ .There exists a bivariate extension ( ) defined on the reference triangle T from (22) such that it has the next properties w.r.t. the weight ê α ω : ( )

Investigation of the Estimators
Theorem 1.Let e be an interface edge in h  Γ .Define for the weight exponent α :

( ) ( )
: We have for any γ as in ( 50) and (51) the next estimate using the patch ( ) Under the assumption that 2 h Cp ≤ , we have the following bound: ( ) ( ) Hence, we deduce From the equality of A and B, one obtains for the interface edge ( ) } ( ) ( ) Consequently, from (52) and the definition (37) of the interface estimator ,e α η Γ , one obtains for each fixed polynomial degree By inserting T f in the last right hand side and by multiplying with ( ) ( ) Combine with the estimation of the local residual to obtain: ( By regrouping the terms, one obtains Then, for the local load oscillation By using the assumption 2 h Cp ≤ , one concludes the last estimate in the theorem. Theorem 2. Let e be an internal edge in h  int (resp.an external edge in h  ext ).Under the assumption that 2 h Cp ≤ , we have the following bound: and respectively and proceed as in the former proof.
Theorem 3. Let the domain oscillation be: and the interface oscillation be: The next estimate holds for the weighted error estimator Proof.Due to (29), one has for ( )

( ) ( )
Use the interpolation property of the Clément interpolant [17] [19] to obtain Deduce therefore the next estimate ( ) Deduce the results for the vanishing weight 0 α ≡ : Regroup the local terms with respect to the the local edges to obtain For non-vanishing weights 0 α ≠ , one applies the polynomial inverse esti- mate (47) to the polynomial in order to obtain By using the same thing for the external domain, obtain ( ) ( ) One has for a fixed polynomial degree 1 p ≥ the following estimates Proof.The following equalities hold: ( ) ) Concerning the estimation of ( ) For the first term, apply the polynomial inverse estimate (48) to polynomial : and the affine transform : As for the second term, split A combination of (111), ( 113) and (114) yields As a consequence, one has from which (101) follows.Thus, one deduces By using the definition (35) of the interior estimator, one has The other estimates are obtained in a similar manner.

Practical Results
The computer implementation is performed by a combination of C functions and C++ classes.Some LAPACK and BLAS routines are used sometimes to perform various linear algebraic operations.

Exact Precision
In this subsection, we concentrate fully on the exact precision for the purpose of obtaining insight and confidence about the accuracy of the results of the computer implementation.That is, we do not consider yet any description of the error estimator α η .We examine several parameters comprising the polynomial degree and the problem coefficients , respectively.Since the FEM-level is used extensively, we introduce it very rapidly.For a given mesh h  on the current level  , another mesh on the next level ( ) is constructed by subdividing every edge in the middle of which one inserts a new node.Therefore, that corresponds to a global uniform refinement where every element is locally subdivided.Only the mesh on the coarsest level (level 1 =  in our entire study) is provided.Thus, one level incrementation amounts to reducing the mesh size from h to h/2.For the numerical compu-tations, we consider the unit square [ ] 2 0,1 as the entire domain Ω while the internal 1 3, 2 3 .The used mesh h  is a tensor product uniform triangular mesh.The exact solution is chosen to be ( ) ( ) ( ) for the internal domain while it is ( ) ( ) ( ) x over all samples and over all elements will then be the practical ∞  -norm.
The first configuration for the investigation of the exact precision corresponds which associates to both the ε-ratio and the μ-ratio the value of 1:1000 .We use these parameters because they highlight a situation where the internal and the external coefficients are very dissimilar.In particular, the ratio of the normal derivatives at the interface is proportional to the ε-ratio.Parameters tending to unity are very easy because that turns out to be the same as treating a problem without an interface in the whole domain Ω .The results of the computation are collected in Table 1 where we consider three fixed polynomial degrees 1, 2,3 p = .For each degree, the FEM-levels are allowed to vary from one to five.The corresponding 2  - norms, 1  -seminorms and ∞  -norms are depicted there as well as the contraction which is the ratio between two consecutive errors as the FEM-level is incremented.The 1  -seminorm errors are in general two digit worse than the 2  -norm.That is observed for different configurations related to the problem coefficients , , , and different simulation parameters.The ∞  - error is just a little worse than the 2  -error but not as imprecise as the error in the 1  -seminorm.Since the 1  -norm is quadratically the sum of the 2  - norm and the 1  -seminorm as Ω , the 1  -norm and the 1  -seminorm are practically of the same order because the 2  -norm is dominated by the 1  -seminorm.In the next description, we will be mainly interested in the contrast to the previous case, the value of ε int is currently smaller that ε ext , although the ε-ratio remains the same.Nonetheless, the general characteristic of the errors remains unchanged which demonstrates the robustness of the method when put in practice.
Our next simulation about a-priori error estimates focuses on the influence of the polynomial degree p.For that, we consider four configurations corresponding to [ ] , , , 0.5, 2, 4, 0.1 ,500, 0.5, 7 Both the 2  -error and the 1  -error decrease in exponential rate with respect to the polynomial degrees.In fact, the 2  -errors drop from an order of 0.1 to an order of 1.0 07 e − in the range of 1, , 6 p =  while the 1  -errors drop from an order of five to an order of 1.0 5 e − as the polynomial degree grows.
In addition, the decrease of the exact precision agrees well with the decrease of the error estimator α η as the FEM-levels grow.That can be very well observed for various values of the fixed polynomial degree p.In fact, the estimations by α η are somewhat influenced by the values of the chosen α .Nonetheless, the values of α η are comparatively of the same order up to some constant scaling factors.The same tests but for other configurations yield the results in Table 4 where we have   At this point, we do no know exactly which value of the parameter α to be chosen.There is no mathematical study to adjust it, let alone to find its optimal value.It is conceivable to choose the value of α adaptively but we do not have  which yield the ε-ratio and the μ-ratio of 1:10 and 1: 2.5 respectively.An observation is that for both configurations, the curves for the exact accuracies and those for the estimated ones almost agree in shape and slope.That means that the exact accuracy is comparable to the estimated precision up to a constant factor.This highlights that using the average estimator avr differ from the first test, the estimated precisions are in general somewhat above the exact precisions.

Conclusion
We considered an interface problem where the internal and external PDEs admit different coefficients.The solution to the PDE is globally continuous but it may admit highly discontinuous derivatives across the interface separating the
the two elements which are incident upon e. Designate by  e σ the extension as in (49) of the polynomial the Green identity, which describes the partial integration of a function with respect to a domain and its boundary, on e T int and e T ext yields the other hand, one has using the definition of the extension together with the properties (50) and (51), one has  of (68)-(71) with the expression of B yields ( )

ω
and apply the inverse estimate (47) to the polynomial (112) to obtain one.It is globally continuous and it admits highly discontinuous derivatives at the interface = ∂ int Γ Ω depending on the values of ε int and ε ext .The right hand side function is obtained by applying the transmission operators (2) and (3) to the exact solution.The ∞  -norm is practically obtained by considering a very dense point sample { }

1 
following ε-ratios 1: 4 , 1: 333.3 , 1: 4 , 1: 333.3 while the μ-ratios are 1: 40 , 1:14 , 1: 4 , 1:14 .The corresponding 2  -errors are depicted in Figure2(b) while the errors using the 1  -seminorm are displayed in Figure2(a).This situation confirms again the fact that the 2  -errors are about two digit worse than the -errors.In those figures, the errors are computed in function of increasing polynomial degrees which range from one to six.We observe satisfactory error decrease for all considered four configurations.

.
Comparable characteristics as in the preceding investigation are observed.That reveals again that the estimator is robust with respect to the problem configurations.

5 Nη
that utility at the time being.As shown formerly, all values of α in [ ] 0,1 provide an efficient a-posteriori error estimator.Thus, the purpose of the next study is to examine the comparison of the exact accuracy in term of the1   - norm on the the one hand with the average of the estimators α η on the other hand.In our case, the average is expressed as = while the values of α are { } 0, 0.25, 0.5, 0.75,1 .We examine in Figure3(a) and Figure3(b) the agreements of the exact precision when the FEM-level increases.In fact, we consider five levels for each test.Each curve corresponds to a fixed polynomial degree 1, 2,3, 4 p =.The problem parameters for the coefficients in the transmission PDE are the first test which is displayed in Figure3(a).Those coefficients correspond to 1:1000 and 1:13

Figure 3 .
Figure 3.Comparison of the exact accuracy and the average avr η in function of the FEM-levels:

η
makes sense to evaluate the desired precision.The constant factors are problem dependent in our case as they are affected by the problem parameters , it can be observed in Figure 3(a), the estimated accuracies are in general somewhat smaller than the exact ones.In contrast, for the second test in Figure 3(b) where only the coefficients 2

ext int ext int ext int int int int ext ext ext ext int int int ext ext ext int int int ext ext ext
)

Table 1 .
1000 and 1:13 .The outcome of the errors is collected in Table2.In Exact precision and contraction ratio for

Table 2 .
Exact precision and contraction ratio for

Table 3 .
Exact and estimated accuracies for Estimated accuracy α

Table 4 .
Exact and estimated accuracies for