In most aspects of numerical simulation, it is desirable to provide an approxi- mation 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-Boltz- mann Equation) for the interaction of solute and solvent media which are respectively denoted by and . The surface represents the solute-solvent    interface which is the molecular surface in the realistic case. The solvent is represented by a continuous dielectric medium while the solute is located inside the cavity . In the sequel, the whole solute-solvent domain is denoted by . The nonlinear PBE  admits the next general expression
Figure 1. Motivation: (a) Ionic solution; (b) Connolly surface: composition of trimmed toroidal and spherical surfaces; (c) Tetrahedral mesh of the solvent.
unknown function is the dimensionless electrostatic potential u. The coefficients and are space-dependent functions related to the dielectric value and the modified Debye-Hückle parameter. Those coefficients might be discontinuous between and but the solution u is required to be continuous everywhere. In the case and under certain assumption on the parameters, the exponential term becomes a hyperbolic sine. By the Taylor expansion, one has . The linear PBE considers only the first term of the above expansion to obtain for 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  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     for isotropic grids. In the context of anisotropic meshes    , there are a variety of APEEs for Poisson and reaction-diffusion problems   . An article  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  . 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   . For the Spectral Element Method, we find in  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  to treat hp-FEM  for the Poisson problem where corner singularities are allowed. We have presented in  an APEE based on hierarchical space enrichment on anisotropic FEM which is combined with adaptive refinements. Boundary Element Method (BEM) is very efficient      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  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     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  ). 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    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  from FEM as illustrated in Figure 1(c). Furthermore, a real chemical simulation by using wavelet BEM is described in  for the quantum computation. A wavelet BEM simulation using domain decomposition techniques was described in  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  for an interface problem. A simulation for chemical quantum computation using FEM is documented in  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.
2. 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:
where designates the normal vector at pointing outward of . The load function and the flux jump are given while the interface and the boundary are polygons in . We consider a mesh of the entire domain such that the restrictions of the mesh in and are respectively denoted by and . The mesh is composed of triangles admitting the next properties:
• The intersection of two different elements 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 ,
• All edges of the interface are edges of the mesh .
For a triangle , we denote
We use to denote for sufficiently large i. We assume quasi-uniformity in the sense that there exists a constant such that
We define the following mutually disjoint subsets of edges
Note that an edge of may have an endpoint in . Likewise, an edge of is allowed to have an endpoint in or . We introduce in addition the set of all edges
For an edge , we denote
The direction of the normal vector is outward if the edge while it is pointed toward the exterior of if the edge . For all other edges in , the normal vectors are pointed in an arbitrary but fixed orientation.
For any triangle T, the affine invertible mapping from the unit reference
onto T is denoted by in which
That allows one to derive results on the unit reference triangle and to carry them over to any element T in the original mesh . We use the standard definitions of the Sobolev spaces for , and . From here onward, we use the usual shorthand if there is a constant c such that in which c is independent on h and p. In addition, amounts to .
We want to consider now the Galerkin variational formulation. Denote the restriction of the solution u to the interface problem in and by and 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 yields
Taking into account the flux jump (5) in the transmission equation, we have
The Galerkin weak form is therefore
For a fixed polynomial degree , the finite dimensional space is
The discrete Galerkin approximation is to search for such that
Introduce the bilinear form
In order to express the a-posteriori estimators, we assume that the appro- ximated solution on the current discretization is available and we consider a parameter . The 1D-weight on and the 2D-weight on are respectively
For a general edge e and triangle T, transformations from the reference elements are used to define and . For an interior element , the estimator is defined as
where designates the -projection of the load function onto the element T. The expression for an exterior element is similar:
For an interface edge , one introduces
where is the -projection of the flux jump Q onto the edge e. The estimator for an interior edge having a normal vector is defined as
where 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 -norm. The estimator corresponding to an exterior edge is
Since one needs computable local estimators for an element-by-element computation, an interior element is introduced
Likewise, for an exterior element , the local estimator reads:
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       .
Lemma 1. Given such that and some . For every univariate polynomial of degree on the 1D reference element , one has
On the 2D reference element (22), one has for every bivariate polynomial of degree
The constants are , , which do not depend on p.
Lemma 2. Consider a univariate polynomial of degree defined on the unit reference interval and a parameter . There exists a bivariate extension defined on the reference triangle from (22) such that it has the next properties w.r.t. the weight :
3. Investigation of the Estimators
Theorem 1. Let be an interface edge in . Define for the weight exponent :
We have for any as in (50) and (51) the next estimate using the patch :
Under the assumption that , we have the following bound:
Proof. Let and be the two elements which are incident upon e. Designate by the extension as in (49) of the polynomial
over the whole patch such that we have the restriction and such that we have the boundary value . An application of the Green identity, which describes the partial integration of a function with respect to a domain and its boundary, on and yields
Introduce the expressions:
Apply now (56) and (57) to the last identity in order to obtain where
By adding on both sides of , one has on the one hand
On the other hand, one has
Hence, we deduce
By using the definition of the extension together with the properties (50) and (51), one has
A combination of (68)-(71) with the expression of yields
From the equality of A and B, one obtains for the interface edge :
Consequently, from (52) and the definition (37) of the interface estimator , one obtains for each fixed polynomial degree :
By inserting in the last right hand side and by multiplying with , one deduces
Combine with the estimation of the local residual to obtain:
By regrouping the terms, one obtains
Then, for the local load oscillation and the local weighted flux oscillation one has
By using the assumption , one concludes the last estimate in the theorem.
Theorem 2. Let be an internal edge in (resp. an external edge in ). Under the assumption that , we have the following bound:
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 from the bilinear form (33)
Every side of is represented as an edge of which is categorized in , , , . As a consequence,
where because . Therefore,
In particular, for where is the Clément interpolant in   . Since , one has .
Use the interpolation property of the Clément interpolant   to obtain
Deduce therefore the next estimate
Deduce the results for the vanishing weight :
Regroup the local terms with respect to the the local edges to obtain
For non-vanishing weights , 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 . Hence,
Theorem 4. For an element , respectively , introduce:
One has for a fixed polynomial degree the following estimates
Thus, the local expressions and verify:
Proof. The following equalities hold:
Consequently, one obtains
Concerning the estimation of , one has
For the first term, apply the polynomial inverse estimate (48) to the polynomial
and the affine transform in order to obtain
As for the second term, split into and , use the boundedness of and apply the inverse estimate (47) to the polynomial (112) to obtain
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.
4. 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.
4.1. 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 . The ε-ratio and the μ-ratio are the positive constants and respectively. Since the FEM-level is used extensively, we introduce it very rapidly. For a given mesh 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 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 as the entire domain while the internal domain is . The used mesh is a tensor product uniform triangular mesh. The exact solution is chosen to be for the internal domain while it is for the external one. It is globally continuous and it admits highly discontinuous derivatives at the interface depending on the values of and . 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 inside every element . The maximum values of over all samples and over all elements will then be the practical -norm.
The first configuration for the investigation of the exact precision corresponds to which associates to both the ε-ratio and the μ-ratio the value of . 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 . For each degree, the FEM-levels are allowed to vary from one to five. The corresponding - norms, -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 -seminorm errors are in general two digit worse than the -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 -error but not as imprecise as the error in the -seminorm. Since the -norm is quadratically the sum of the - norm and the -seminorm as , the -norm and the -seminorm are practically of the same order because the -norm is dominated by the -seminorm. In the next description, we will be mainly interested in the -norm or practically the -seminorm. We observe practically constant contraction numbers depending on the used polynomial degree. The same test has been conducted for another configuration of whose ε-ratio and μ-ratio are respec- tively and . The outcome of the errors is collected in Table 2. In
Table 1. Exact precision and contraction ratio for .
Table 2. Exact precision and contraction ratio for .
contrast to the previous case, the value of is currently smaller that , 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 , , , . They have respectively the following e-ratios , , , while the μ-ratios are , , , . The corresponding -errors are depicted in Figure 2(b) while the errors using the -seminorm are displayed in Figure 2(a). This situation confirms again the fact that the -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. Both the -error and the -error decrease in exponential rate with respect to the polynomial degrees. In fact, the -errors drop from an order of 0.1 to an order of in the range of while the -errors drop from an order of five to an order of as the polynomial degree grows.
Figure 2. A-priori error in function of the increasing polynomial degree for various values of : (a) -error, (b) -error.
4.2. Practical A-Posteriori Error Estimation
Now that we have gained insight about the a-priori accuracy, we want to turn our attention to the a-posteriori estimator . The principal role of an a-posteriori error estimator is twofold. For one, it serves as gaining some idea of whether to continue or to abort a simulation. The computation is aborted when the desired precision is provided by the estimator. Since the exact solution u is not known for real applications, an estimator is needed. A further purpose of the error estimator is to identify the regions within the domain where the precision is unsatisfactory. Mesh refinements are therefore applied at those regions to improve the local precision. The proposed estimator can be used for both purposes. But we have currently only the first utility in our computer implementation. We use different values of which comprise , , , , . The results of the computer simulation is collected in Table 3 where we examine four polynomial degrees . For every fixed polynomial degree p, the value of the FEM-level is allowed to
Table 3. Exact and estimated accuracies for .
range from one to five. We consider a configuration where the problem coefficients are that corresponds to an e-ratio and μ-ratio of and . According to the proposed theory, it is sufficient to conduct a comparison between the -error and the error estimator . For all values of , Table 3 highlights that the estimator provides an efficient estimation of the exact precision as predicted theoretically. 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 . Comparable charac- teristics as in the preceding investigation are observed. That reveals again that the estimator is robust with respect to the problem configurations.
Table 4. Exact and estimated accuracies for .
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 that utility at the time being. As shown formerly, all values of in 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 the - norm on the the one hand with the average of the estimators on the other
hand. In our case, the average is expressed as where we set
while the values of are . We examine in Figure 3(a) and Figure 3(b) the agreements of the exact precision and estimated precision when the FEM-level increases. In fact, we consider five levels for each test. Each curve corresponds to a fixed polynomial degree . The problem parameters for the coefficients in the transmission PDE are for the first test which is displayed in Figure 3(a). Those coefficients correspond to and
Figure 3. Comparison of the exact accuracy and the average in function of the FEM-levels: (a) ; (b) .
as ε-ratio and μ-ratio. As for the second test, we consider the problem coefficients which yield the ε-ratio and the μ-ratio of and 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 makes sense to evaluate the desired precision. The constant factors are problem dependent in our case as they are affected by the problem parameters . As 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 differ from the first test, the estimated precisions are in general somewhat above the exact precisions.
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 internal and the external domains. We proposed an a-posteriori error estimator which has been theoretically demonstrated to be equivalent to the exact precision. The performance of the estimator has also been investigated practically. For every fixed polynomial degree, it provides satisfactory precisions which are comparable to the order of the exact accuracies on different FEM-levels. Several tests have been performed where one varies the problem configurations and the simulation parameters.
 Cancès, E., Mennucci, B. and Tomasi, J. (1998) A New Integral Equation Formalism for the Polarizable Continuum Model, Theoretical Background and Applications to Isotropic and Anisotropic Dielectrics. The Journal of Chemical Physics, 107, 3032-3041.
 Teso, A., Filho, E. and Neto, A. (1997) Solution of the Poisson-Boltzmann Equation for a System with Four Ionic Species. Journal of Mathematical Biology, 35, 814-824.
 Xie, D. and Zhou, S. (2007) A New Minimization Protocol for Solving Nonlinear Poisson-Boltzmann Mortar Finite Element Equation. BIT Numerical Mathematics, 47, 853-871.
 Verfürth, R. (1994) A-Posteriori Error Estimation and Adaptive Mesh Refinement Techniques. Journal of Computational and Applied Mathematics, 50, 67-83. https://doi.org/10.1016/0377-0427(94)90290-9
 Jou, J. and Liu, J. (2000) An A-Posteriori Finite Element Analysis for the Stokes Equations. Journal of Computational and Applied Mathematics, 114, 333-343.
 Kunert, G. and Verfürth, R. (2000) Edge Residuals Dominate A-Posteriori Error Estimates for Linear Finite Element Methods on Anisotropic Triangular and Tetrahedral Meshes. Numerische Mathematik, 86, 283-303. https://doi.org/10.1007/PL00005407
 Kunert, G. (2001) Robust a Posteriori Error Estimation for a Singularly Perturbed Reaction-Diffusion Equation on Anisotropic Tetrahedral Meshes. Advances in Computational Mathematics, 15, 237-259.
 Creusé, E., Kunert, G. and Nicaise, S. (2004) A-Posteriori Error Estimation for the Stokes Problem, Anisotropic and Isotropic Discretizations. Mathematical Models and Methods in Applied Sciences, 14, 1297-1341. https://doi.org/10.1142/S0218202504003635
 Randrianarivony, M. (2004) Anisotropic Finite Elements for the Stokes Problem, A-Posteriori Error Estimator and Adaptive Mesh. Journal of Computational and Applied Mathematics, 169, 255-275. https://doi.org/10.1016/j.cam.2003.12.025
 Costabel, M. and Stephan, E. (1985) A Direct Boundary Integral Equation Method for Transmission Problems. Journal of Mathematical Analysis and Applications, 106, 367-413. https://doi.org/10.1016/0022-247X(85)90118-0
 Holst, M., McCammon, J., Yu, Z., Zhou, Y. and Zhu, Y. (2012) Adaptive Finite Element Modeling Techniques for the Poisson-Boltzmann Equation. Communications in Computational Physics, 11, 179-214. https://doi.org/10.4208/cicp.081009.130611a
 Lu, B., Zhou, Y., Holst, M. and McCammon, J. (2008) Recent Progress in Numerical Methods for the Poisson-Boltzmann Equation in Biophysical Applications. Computer Physics Communications, 3, 973-1009.
 Baker, N., Sept, D., Holst, M. and McCammon, J. (2001) The Adaptive Multilevel Finite Element Solution of the Poisson-Boltzmann Equation on Massively Parallel Computers. IBM Journal of Research and Development, 45, 427-438. https://doi.org/10.1147/rd.453.0427
 Chen, L., Holst, M. and Xu, J. (2007) The Finite Element Approximation of the Nonlinear Poisson-Boltzmann Equation. SIAM Journal on Numerical Analysis, 45, 2298-2320.
 Randrianarivony, M. and Brunnett, G. (2008) Preparation of CAD and Molecular Surfaces for Meshfree Solvers. Lecture Notes in Computational Science and Engineering, 65, 231-245. https://doi.org/10.1007/978-3-540-79994-8_14
 Diedrich, C., Dijkstra, D., Hamaekers, J., Henninger, B. and Randrianarivony, M. (2015) A Finite Element Study on the Effect of Curvature on the Reinforcement of Matrices by Randomly Distributed and Curved Nanotubes. Journal of Computational and Theoretical Nanoscience, 12, 2108-2116. https://doi.org/10.1166/jctn.2015.3995