Serge Nicaise
email: serge.nicaise@univ-valenciennes.fr

Hengguang Li

Anna Mazzucato

Regularity and a priori error analysis of a Ventcel problem in polyhedral domains

Keywords: AMS (MOS) subject classification (2010): 35J25, 65N30, 46E35, 52B70, 58J05 Elliptic boundary value problems, Ventcel boundary conditions, polyhedral domains, weighted Sobolev spaces, Finite Element, anistropic meshes

We consider the regularity of a mixed boundary value problem for the Laplace operator on a polyhedral domain, where Ventcel boundary conditions are imposed on one face of the polyhedron and Dirichlet boundary conditions are imposed on the complement of that face in the boundary. We establish improved regularity estimates for the trace of the variational solution on the Ventcel face, and use them to derive a decomposition of the solution into a regular and a singular part that belongs to suitable weighted Sobolev spaces. This decomposition, in turn, via interpolation estimates both in the interior as well as on the Ventcel face, allows us to perform an a priori error analysis for the Finite Element approximation of the solution on anisotropic graded meshes. Numerical tests support the theoretical analysis.

Introduction

This article concerns the regularity of solutions to an elliptic boundary-value problem for the Laplace operator on a polyhedral domain in R 3 under so-called Ventcel or Wentzell boundary conditions. The regularity result we establish in weighted Sobolev spaces gives rise, in turn, to a priori error estimates for the Finite Element Method (FEM) on a suitable anisotropic mesh.

We first introduce the Ventcel boundary-value problem. Let Ω be a bounded domain of R 3 with Lipschitz boundary Γ and let Γ V be an open subset of Γ with positive measure. We denote by Γ D = Γ \ Γ V the complement of Γ V , which we assume has also positive measure.

We consider the following mixed boundary-value problem:

-∆u = f, in Ω, (1a) u = 0 on Γ D , (1b)
-∆ LB u + ∂ ν u = g, on Γ V , (1c)
which can be thought of as a generalized Robin-type boundary condition. The more general problem is well posed only under conditions on the sign of α and β. Ventcel boundary conditions arise naturally in many contexts. In the context of multidimentional diffusion processes, Ventcel boundary conditions were introduced in the pioneering work of Ventcel [START_REF] Ventcel | Semigroups of operators that correspond to a generalized differential operator of second order[END_REF][START_REF] Ventcel | On boundary conditions for multi-dimensional diffusion processes[END_REF] (see also the work of Feller for one-dimensional processes [START_REF] Feller | The parabolic differential equations and the associated semi-groups of transformations[END_REF][START_REF] Feller | Generalized second order differential operators and their lateral conditions[END_REF]). They can model heat conduction in materials for which the boundary can store, but not absorb or transmit heat. They can also be derived as approximate boundary conditions in asymptotic problems or artificial boundary conditions in exterior problems (see e.g. [START_REF] Bonnaillie-Noël | On generalized Ventcel's type boundary conditions for Laplace operator in a bounded domain[END_REF][START_REF] Bonnaillie-Noël | Artificial conditions for the linear elasticity equations[END_REF][START_REF] Nicaise | Polynomial stabilization of the wave equation with Ventcel's boundary conditions[END_REF] and references therein), in particular in fluid-structure interaction problems.

Problem ( 1) is known to have a unique variational solution if f and g are in the appropriate Sobolev space as recalled in Section 2. We are concerned here with the higher regularity for solutions to this problem when the data is also regular, in the case that the domain Ω is a polyhedral domain in R 3 . It is well known that, due to the presence of edges and corners at the boundary of Ω, even when Γ V is empty, elliptic regularity does not hold, and the solution is not smooth even if the data is smooth. This loss of regularity affects the rate of convergence of the Finite Element approximation to the solution if uniform meshes are used.

By using weighted Sobolev spaces, where the weights are the distance to the edges and vertices, respectively, one can characterize precisely the behavior of the variational solution near the singular set in terms of singular function and singular exponents (Theorem 2.5). In turn, the decomposition of the solution into a regular and a singular part, together with interpolation estimates (Theorem 3.1), leads to establishing a priori error estimates for the Finite Element approximation (Corollary 3.4), where the elements are given on an anisotropic mesh that exploits the improved regularity of the solution along the edges versus the corners of the polyhedron. There is a well established literature on this approach for mixed Dirichlet, Neumann, and even standard Robin boundary condition (see for example [START_REF] Apel | The finite element method with anisotropic mesh grading for elliptic problems in domains with corners and edges[END_REF][START_REF] Bacuta | Improving the rate of convergence of high-order finite elements on polyhedra. II. Mesh refinements and interpolation[END_REF][START_REF] Schötzau | hp-dGFEM for second-order mixed elliptic problems in polyhedra[END_REF]). There are also several works in the literature concerning the Ventcel boundary-value problems on singular domains (see in particular [START_REF] Lemrabet | Problème aux limites de Ventcel dans un domaine non régulier[END_REF][START_REF] Popivanov | On Ventcel's type boundary condition for Laplace operator in a sector[END_REF]), and their implementation of the FEM (see [START_REF] Kashiwabara | Well-posedness, regularity, and convergence analysis of the finite element approximation of a generalized Robin boundary value problem[END_REF] and references therein). The novelty of this work consists in extending the approach using weighted spaces and anisotropic meshes to the Ventcel boundary conditions, which include tangential differential operators at the boundary of the same order as the main operator in Ω. As a matter of fact, the main difficulty in considering such boundary conditions lies in establishing the needed regularity of the traces on the faces of the polyhedron. For simplicity, we restrict here to the case where the Ventcel condition is imposed on only one face of the polyhedron. If the Ventcel condition is imposed on adjacent faces, one would expect higher regularity to hold for the solution on these faces, under suitable transmission conditions at the common edges. however, capturing this behavior entails studying weighted Sobolev spaces for which the weight is the distance to the boundary and not the distance to the singular set (as those arising from the analysis of equations with degenerate coefficients). We reserve to address this problem in future works.

The paper is organized as follows. In Section 2, we recall the variational formulation for Problem [START_REF] Mazzucato | Nonlinear waves in networks[END_REF], and prove our main regularity result for the solution in weighted spaces. In Section 3, we introduce the anisotropic mesh and the associated Finite Element discretization of the problem, and derive a priori error estimates. Section 4 contains some refined 2D interpolation estimates valid on the polyhedral faces, needed for the error analysis. We close in Section 5 by presenting some numerical examples to validate the theoretical analysis.

We end this Introduction with some needed notation.

If Ω is a domain of R n , n ≥ 1, we employ the standard notation H m (Ω) to denote the Sobolev space that consists of functions whose ith derivatives, for 0 ≤ i ≤ m, are square-integrable. The L 2 (Ω)-inner product (resp. norm) will be denoted by (•, •) Ω (resp. • Ω ). The usual norm and semi-norm in H s (Ω), for s ≥ 0, are denoted by • s,Ω and | • | s,Ω , respectively. The trace operator from H 1 (Ω) into H 1 2 (∂Ω) will be denoted by γ. We also introduce the space:

H 1 Γ D (Ω) = {u ∈ H 1 (Ω) : γu = 0 on Γ D }, which is clearly a closed subspace of H 1 (Ω). If v is a d-dimensional

Some regularity results

In this section we recall needed facts about the well-posedness of the Ventcel Problem (1), and establish regularity estimates for its variational solution in weighted spaces. The variation formulation of (1) is well known (see [START_REF] Mazzucato | Nonlinear waves in networks[END_REF][START_REF] Lemrabet | Problème aux limites de Ventcel dans un domaine non régulier[END_REF][START_REF] Kashiwabara | Well-posedness, regularity, and convergence analysis of the finite element approximation of a generalized Robin boundary value problem[END_REF]). We let

V := {u ∈ H 1 Γ D (Ω) : γu ∈ H 1 0 (Γ V )},
which is a Hilbert spaces equipped with the natural norm

u 2 V := u 2 1,Ω + |γu| 1,Γ V , ∀u ∈ V.
We further introduce the bilinear form

a(u, v) = Ω ∇u • ∇v dx + Γ V ∇ T (γu) • ∇ T (γv) dσ(x), ∀u, v ∈ V.
As this bilinear form is continuous and coercive in V , by the Lax-Milgram lemma, for any f ∈ L 2 (Ω) and g ∈ L 2 (Γ V ), there exists a unique solution u ∈ V of (2)

a(u, v) = Ω f v dx + Γ V gv dσ(x), ∀v ∈ V.
It was shown in [START_REF] Kashiwabara | Well-posedness, regularity, and convergence analysis of the finite element approximation of a generalized Robin boundary value problem[END_REF]Thm 3.3] that if Γ D is empty and if Γ is C 1,1 , then u belongs to H 2 (Ω) and γu belongs to H 2 (Γ V ). This regularity is no longer valid if Ω is a non-convex polyhedral domain, and the main purpose of this section is to describe the behavior of the solution near the singular set, which consists of the edges and vertices of the boundary faces of the polyhedron, and characterize the regularity of boundary traces of u and its derivatives.

To this end, we will employ anisotropic weighted Sobolev spaces, for which the weights are (variants of) the distance to the edges and vertices, respectively. There is a vast literature concerning the use of weighted Sobolev spaces in the analysis of singular domains (we refer for instance to [START_REF] Dauge | Elliptic boundary value problems on corner domains[END_REF][START_REF] Kozlov | Elliptic boundary value problems in domains with point singularities[END_REF][START_REF] Maz | Elliptic boundary value problems on manifolds with singularities[END_REF][START_REF] Maz | Weighted L p estimates of solutions to boundary value problems for second order elliptic systems in polyhedral domains[END_REF][START_REF] Nazarov | Elliptic problems in domains with piecewise smooth boundaries[END_REF] and references therein). In the context of the analysis of Dirichlet/Neumann boundary conditions, anisotropic Sobolev spaces were used in [START_REF] Apel | Anisotropic interpolation with applications to the finite element method[END_REF][START_REF] Apel | The finite element method with anisotropic mesh grading for elliptic problems in domains with corners and edges[END_REF][START_REF] Apel | Anisotropic mesh refinement in polyhedral domains: error estimates with data in L 2 (Ω)[END_REF][START_REF] Bacuta | Improving the rate of convergence of high-order finite elements on polyhedra. II. Mesh refinements and interpolation[END_REF].

From now on we assume that Ω is a polyhedral domain of the space and that Γ V is reduced to one face F of the boundary.

By a face, we mean an open face on the boundary. Let S and E be the set of vertices and the set of open edges of Ω, respectively.

On the polygonal face F , we denote its set of vertices by S F . Given a vertex S ∈ S F , we denote by (r S , θ S ) the radial distance and angular component of the local polar coordinate system centered at S on the plane containing F . In addition, we let ω F,S be the interior angle on the face F associated with the vertex S.

Following [START_REF] Apel | The finite element method with anisotropic mesh grading for elliptic problems in domains with corners and edges[END_REF], we consider a triangulation {Λ } L =1 of the domain Ω that consists of disjoint tetrahedra Λ . We will refer to each tetrahedron Λ as a macro element, to distinguish it from the elements of the mesh utilized in the analysis of the FEM in Section 3. The purpose of the macro elements is to localize the construction and the regularity estimates near edges and vertices of Ω. We will also refer to any edge or vertex of an element Λ as a singular edge or singular vertex, if that edge or vertex lies along a true edge or is a true vertex of Ω and the solution is not in H 2 near that true edge or vertex.

We will assume that each Λ contains at most one singular edge and at most one singular vertex. If Λ contains both a singular edge and a singular vertex, that vertex belongs to that edge. We will also assume that all Λ are shape regular with diameter of order O(1). In each macro element Λ , we introduce a local Cartesian coordinate system x ( ) = (x

( ) 1 , x ( ) 2 , x ( )
3 ), such that the singular vertex, if it exists, is at the origin, and the singular edge, if it exists, lies along the x ( ) 3 -axis. We then define the distance functions to the set of singular edges and singular vertices, respectively, as follows:

r ( ) (x ( ) ) = (x ( ) 1 ) 2 + (x ( ) 2 ) 2 , (3a) R ( ) (x ( ) ) = (x ( ) 1 ) 2 + (x ( ) 2 ) 2 + (x ( ) 3 ) 2 , (3b)
and introduce the auxiliary function

θ ( ) (x ( ) ) = r ( ) (x ( ) )/R ( ) (x ( ) ). (3c)
We observe that r ( ) , and R ( ) extend as continuous functions to the closure of the macro element Λ , while θ ( ) extends as a bounded function.

In what follows, we will omit the sup-index ( ) in these distance functions and in x, when there is no confusion about the underlying macro element. Given a subdomain Λ ⊂ Ω, we define the following weighted Sobolev space for k ∈ N and β, δ ∈ R:

V k β,δ (Λ) := {v ∈ D (Λ), v V k β,δ (Λ) < ∞}, where v 2 V k β,δ (Λ) = |α|≤k R β-k+|α| θ δ-k+|α| ∂ |α| v 2 L 2 (Λ) ,
and R(x) = R ( ) (x ( ) ) and θ(x) = θ ( ) (x ( ) ), given in (3b) and in (3c), if x ∈ Λ is represented by x ( ) in local coordinates. We will also need to define spaces on the faces of Ω. To this end, given G a bounded polygonal domain in R 2 , we also define

V k γ (G) := {v ∈ D (G), v V k γ (G) < ∞}, where v 2 V k γ (G) = |α|≤k ρ γ-k+|α| ∂ |α| v 2 L 2 (G) .
Above, ρ is the distance function to the set of vertices of G, defined in a manner similar to R above.

We further classify the initial macro elements Λ into four types as follows:

Type 1. Λ ∩ (S ∪ E) = ∅; Type 2. Λ ∩ S = ∅ but Λ ∩ E = ∅; Type 3. Λ ∩ E = ∅ but Λ ∩ S = ∅; Type 4. Λ ∩ E = ∅ and Λ ∩ S = ∅.
We first start with an improved regularity of ∂ ν u on Γ V . In what follows, for ease of notation we will let u F be the trace γu of u on the face F . Furthermore, for a two-dimensional domain D, we define the space Hs (D), 0 < s < 1, as the closure of

C ∞ c (D) in H s (D). Lemma 2.1 If D ⊂ R 2 is a two-dimensional domain with Lipschitz boundary, then for any h ∈ ( H 1 2 (D)) , the unique solution w ∈ H 1 0 (D) of -∆w = h in D, belongs to H 1+ε (D) ∩ H 1 0 (D) for any ε ∈ (0, 1 2 ). Proof. We fix ε ∈ (0, 1 2 ). Since H 1-ε 0 (D) = H1-ε (D) is continuously and densely embedded into H 1 2 (D), by duality we obtain that ( H 1 2 (D)) is continuously embedded into (H 1-ε 0 (D)) = H -1+ε (D).
Hence, w can be seen as a solution of the Laplace equation with datum in H -1+ε (D). Owing to Theorem 18.13 and Remark 18.17/2 in [START_REF] Dauge | Elliptic boundary value problems on corner domains[END_REF], w belongs to H 1+ε (D). Lemma 2.2 Let u ∈ V be the solution of (2), then we have

∂ ν u ∈ L 2 (Γ V ).
Proof. We first observe that, by Theorem 2.8 of [START_REF] Nicaise | Exact controllability of a pluridimensional coupled problem[END_REF],

∂ ν u ∈ ( H 1 2 (Γ V ))
. Then, we may interpret u F ∈ H 1 0 (F ) as the unique variational solution of (4)

∆ LB u F = -g + ∂ ν u ∈ ( H 1 2 (F )) .
By Lemma 2.1, we deduce that u F belongs to H 1+ε (F ) ∩ H 1 0 (F ) for any ε ∈ (0, 1 2 ). We now fix ε ∈ (0, 1 2 ) small enough that the mapping

H 3 2 +ε (Ω) ∩ H 1 0 (Ω) → H -1 2 +ε (Ω) : v → ∆v,
is an isomorphism (see [START_REF] Dauge | Elliptic boundary value problems on corner domains[END_REF]Thm 18.13]). Now, by applying the trace theorem from [START_REF] Grisvard | Théorèmes de traces relatifs à un polyèdre[END_REF], there exists w ∈ H

3 2 +ε (Ω) such that γw = 0 on Γ D , (5)
γw = u F on F. (6)
This implies, again by uniqueness, that v := u -w ∈ H 1 0 (Ω) is the solution of

-∆v = f + ∆w ∈ H -1 2 +ε (Ω).
We therefore deduce that v belongs to H 3 2 +ε (Ω) and, hence, u belongs to this space as well. By a standard trace theorem, we finally conclude that ∇u ∈ H ε (Γ) 3 .

Thus, we have the following decomposition of the singular solution u F on the polygonal face.

Corollary 2.3 Let again u ∈ V be the solution of (2). Then, it holds

(7) u F = u F,R + S∈S F :ω F,S >π c S r π ω F,S S sin πθ S ω F,S ,
where u F,R ∈ H 2 (F ) and c S ∈ C.

Proof. As ∂ ν u belongs to L 2 (F ) by Lemma 2.2, the right-hand side in ( 4) is now in L 2 (F ) and therefore Theorem 4.4.3.7 of [START_REF] Grisvard | Elliptic problems in nonsmooth domains[END_REF] yields [START_REF] Bonnaillie-Noël | On generalized Ventcel's type boundary conditions for Laplace operator in a bounded domain[END_REF].

We will refer to u F,R as the regular part of u F , hence the subscript, as it has the expected regularity from elliptic theory. We will consequently call u F -u F,R the singular part of u F .

For the regularity of the solution in the interior of the domain Ω, we first have the following lifting estimate based on the trace theorem.

Lemma 2.4 Given u F ∈ H 3 2 +ε (F )∩H 1 0 (F ) for some ε ∈ (0, 1 2
), there exists a lifting w ∈ H 2+ε (Ω) satisfying (5)- [START_REF] Bacuta | Improving the rate of convergence of high-order finite elements on polyhedra. II. Mesh refinements and interpolation[END_REF].

Proof. The idea is to use again the trace theorem from [START_REF] Grisvard | Théorèmes de traces relatifs à un polyèdre[END_REF] with s = 2 + ε and the operator [START_REF] Bonnaillie-Noël | Artificial conditions for the linear elasticity equations[END_REF] u

→ (f j,0 , f j,1 ) N j=1 := u |Γj , (∂ νj u) |Γj ) N j=1
, where Γ j are the faces of Ω and ν j the outward normal vector along Γ j . As 1+ε is not an integer, this trace operator (8) is surjective from H 2+ (Ω) onto the subspace of

N j=1 (H 3 2 +ε (Γ j ) × H 1 2 +ε (Γ j
)) that satisfies the compatibility conditions (C 1 ) of [START_REF] Grisvard | Théorèmes de traces relatifs à un polyèdre[END_REF]. If we assume that Γ 1 = F , it is therefore sufficient to show that there exist f j,1 ∈ H

1 2 +ε (Γ j ), j = 1, • • • , N such that (u F , f 1,1 ) × (0, f j,1 ) N j=2
satisfies these conditions (C 1 ). Since such conditions are quite technical to check, as in [START_REF] Grisvard | Théorèmes de traces relatifs à un polyèdre[END_REF] we can reduce to check such conditions in the case where Ω is the trihedral x i > 0, i = 1, 2, 3 and F is the face x 1 = 0 (and hence N = 3 with Γ 2 ≡ x 2 = 0 and Γ 3 ≡ x 3 = 0), by means of a localization argument and a linear change of variables. In such a case, the conditions (C 1 ) of [START_REF] Grisvard | Théorèmes de traces relatifs à un polyèdre[END_REF] for (u F , f 1,1 ) × (0, f j,1 ) 3 j=2 take the form:

u F = 0 on A 1,2 ∪ A 2,3 , (9a)
f 1,1 = 0 on A 1,3 , (9b)
∂ 2 u F = 0 on A 1,3 , (9c)
∂ 3 u F = f 3,1 on A 1,3 , (9d)
f 1,1 = 0 on A 1,2 , (9e)
∂ 2 u F = f 2,1 on A 1,2 , (9f)
∂ 3 u F = 0 on A 1,2 , (9g)
f 2,1 = 0 on A 2,3 , (9h)
f 3,1 = 0 on A 2,3 , (9i)
where A i,j = Γi ∩ Γj . The first condition trivially holds as u F belongs to H 1 0 (F ), and similarly (9c) (resp. (9g)) because ∂ 2 u F (resp. ∂ 3 u F ) is the tangential derivatives of u F on A 1,3 (resp. A 1,2 ). To satisfy the second and fourth conditions we simply take f 1,1 = 0. Hence it remains to verify the conditions (9f) and (9h) (resp. (9d) and (9i)) that can be interpreted as constraints on f 2,1 and f 3,1 , respectively. In other words, we look for

f 2,1 ∈ H 1 2 +ε (Γ 2 ) (resp. f 3,1 ∈ H 1 2 +ε (Γ 3
)) satisfying the boundary conditions (9f) and (9h) (resp. (9d) and (9i)). Such a solution f 2,1 (and similarly f 3,1 ) exists by applying Theorem 1.5.1.2 of [START_REF] Grisvard | Elliptic problems in nonsmooth domains[END_REF] (valid for a quarter plane), because the function defined by ∂ 2 u F on A 1,2 and 0 on A 2,3 belongs to H ε (Γ 2 ).

For a vertex v ∈ S, let C v be the infinite polyhedral cone that coincides with Ω in the neighborhood of v. Let G v = C v ∩ S 2 (v) be the intersection of C v and the unit sphere centered at v. For an edge e ∈ E, let ω e be the interior angle between the two faces of Ω that contain e. Then, for v ∈ S and for e ∈ E, respectively, we define the following parameters associated to the singularities in the solution near v and e: [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF] λ

v := - 1 2 + λ v,1 + 1 4 , λ e := π/ω e ,
where λ v,1 is the smallest positive eigenvalue of the Laplace-Beltrami operator on G s with Dirichlet boundary conditions. We observe that a vertex v is singular if λ v < 1/2 and an edge e is singular if λ e < 1. For a given macro element Λ , we set λ Theorem 2.5 Let u ∈ V be again the solution of (2). We have:

( ) v = λ v if Λ contains one singular vertex v of Ω and λ ( ) v = ∞
(11) u = u R + u S ,
where u R ∈ H 2 (Ω) and u S ∈ H 1 (Ω) satisfies, for all ∈ {1, . . . , L}, ∂u S ∂x

( ) j ∈ V 1 β,δ (Λ ), j = 1, 2, ( 12
)
∂u S ∂x ( ) 3 ∈ V 1 β,0 (Λ ), (13)
for any β, δ ≥ 0 such that

β > 1 2 -λ ( ) v , δ > 1 -λ ( ) e ,
Again, the subscripts refer to the fact that u R has the expected regularity, and hence it will be called the regular part of the solution, while u S = u -u R represents the singular part.

Proof. The decomposition [START_REF] Bonnaillie-Noël | On generalized Ventcel's type boundary conditions for Laplace operator in a bounded domain[END_REF] implies that there exists ε ∈ (0, 1 2 ) small enough such that

u F ∈ H 3 2 +ε (F ) ∩ H 1 0 (F ).
Hence by Lemma 2.4, there exists a lifting w ∈ H 2+ε (Ω) satisfying ( 5)- [START_REF] Bacuta | Improving the rate of convergence of high-order finite elements on polyhedra. II. Mesh refinements and interpolation[END_REF]. With this lifting at hands, we consider u -w, which belongs to H 1 0 (Ω) and is the weak solution of -∆(u -w) = f + ∆w.

As f + ∆w belongs to L 2 (Ω), we can apply Theorem 2.10 of [START_REF] Apel | The finite element method with anisotropic mesh grading for elliptic problems in domains with corners and edges[END_REF] to u -w, which gives the decomposition:

u -w = u R + u S ,
with u R ∈ H 2 (Ω) and u S satisfying ( 12)- [START_REF] Feller | Generalized second order differential operators and their lateral conditions[END_REF]. Finally, the result follows by setting u R = w + u r .

Remark 2.6 Theorem 2.5 shows that for the solution to (2) with the Ventcel boundary condition, its regularity in Ω, determined by the geometry of the domain, is similar to the regularity of the Poisson equation with the Dirichlet boundary condition. Meanwhile, the trace of the solution u on the face F is the solution of a two-dimensional elliptic problem with the Dirichlet boundary condition. Corollary 2.3 implies that the regularity of the trace depends on the interior angles of the polygon F .

Finite element approximation

We consider an (anisotropic) triangulation T h = {T i } N i=1 of Ω as in Section 3 of [START_REF] Apel | The finite element method with anisotropic mesh grading for elliptic problems in domains with corners and edges[END_REF] or in Section 2 of [START_REF] Apel | Anisotropic mesh refinement in polyhedral domains: error estimates with data in L 2 (Ω)[END_REF], consisting of tetrahedra with refinement parameters µ and ν . We assume the general conditions for a triangulation of the domain (see e.g. [START_REF] Brenner | The mathematical theory of finite element methods[END_REF][START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]) and that the number of tetrahedra m satisfies N ∼ h -3 , where h is the global mesh size. In addition, we assume that the initial subdomains Λ are resolved exactly, namely, Λ = ∪ i∈L Ti , where = 1, • • • , L and L ⊂ {1, • • • , m} is the index set of the tetrahedra included in Λ .

In each Λ , the parameters µ , ν ∈ (0, 1] determine the anisotropic mesh refinement close to edges and vertices, respectively as indicated in ( 14) below. When µ = 1 or ν = 1, there will be no graded refinement in Λ . We recall the local Cartesian coordinate system (x

( ) 1 , x ( ) 2 , x ( )
3 ) in each of the subdomain Λ , which is such that the singular vertex is at the origin and the singular edge is along the x 3 -axis, if they exist. Then, for each element T i ⊂ Λ of the triangularization, we let

r i := inf x∈Ti [(x ( ) 1 ) 2 + (x ( ) 2 ) 2 ] 1/2 , R i := inf x∈Ti [(x ( ) 1 ) 2 + (x ( ) 2 ) 2 + (x ( ) 3 ) 2 ] 1/2 ,
be the distance of T i to the origin and the x 3 -axis, respectively. We then introduce local, anisotropic mesh parameters in T i as follows:

h i := h 1/µ if r i = 0, hr 1-µ i if r i > 0, H i := h 1/ν if 0 ≤ R i h 1/ν , hR 1-ν i if R i h 1/ν , ( 14
)
We also introduce the actual mesh sizes hj,i , which are the lengths of the projections of T i ⊂ Λ on the x ( ) j -axis, 1 ≤ j ≤ 3. Then, there exists a triangulation T h satisfying the following conditions:

1. If µ < 1, then hj,i ∼ h i , j = 1, 2, h3,i H i , and h3,i ∼ H i if r i = 0.
2. The number of tetrahedra in Λ with r i = 0 is of order h -1 .

3. The number of tetrahedra in Λ such that 0 ≤ R i h 1/ν is bounded by h 2µ /ν -2 , and there is only one tetrahedral element with R i = 0.

4. If µ < 1, then µ ≤ ν for 1 ≤ ≤ L.
We refer to [START_REF] Apel | The finite element method with anisotropic mesh grading for elliptic problems in domains with corners and edges[END_REF] or a detailed description of these conditions. It is clear that this triangulation T h induces an exact triangulation F h of the face F , the elements of which are simply given by T ∩ F for T ∈ T h .

Based on these triangulations, we introduce the approximation space V h of V as follows:

V h := {u h ∈ V : u |T ∈ P 1 (T ), ∀T ∈ T h },
where P m , m ∈ Z + , denotes the space of all polynomials of degree ≤ m. This is clearly a closed subspace of V . Then, the Finite Element approximation of Problem (2) consists of looking for a solution

u h ∈ V h of (15) a(u h , v h ) = Ω f v h dx + Γ V gv h dσ(x), ∀v h ∈ V h . By Céa's lemma, we have u -u h V inf v h ∈V h u -v h V ,
Hence an error estimate will be available if we can built an appropriate approximation v h of u. This is the purpose of the next theorems in this section.

Theorem 3.1 Recall the parameters in [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF] and in [START_REF] Grisvard | Théorèmes de traces relatifs à un polyèdre[END_REF]. Assume that for all = 1, • • • , L, we have:

µ < λ ( ) e , (16a) ν < λ ( ) v + 1 2 (16b) 1 ν + 1 µ (λ ( ) v - 1 2 ) > 1. (16c)
Then, there exists v h ∈ V h such that

(17) u -v h 1,Ω h.
Proof. The proof of Theorem 2.5 furnishes the splitting of u as

u = ũ + w,
where ũ = u -w and w ∈ H 2+ε (Ω) with ε ∈ (0, 1 2 ). Hence we define an interpolant I h u of u as follows:

(18)

I h u := ũI + D h (ũ -ũI ) + L h w,
where D h is the interpolant introduced in [START_REF] Apel | Anisotropic mesh refinement in polyhedral domains: error estimates with data in L 2 (Ω)[END_REF], ũI is the Lagrange interpolant of ũ with respect to the partition {Λ }, while L h w is the standard Lagrange interpolant of w, which consists of piecewise polynomials of degree 1. Then, using the regularity estimate in Theorem 2.5 and applying Theorem 3.11 of [START_REF] Apel | Anisotropic mesh refinement in polyhedral domains: error estimates with data in L 2 (Ω)[END_REF], we have

(19) ũ -ũI + D h (ũ -ũI ) 1,Ω h.
On the other hand, as w belongs to H 2+ε (Ω) and

H 2+ε (Ω) is continuously embedded into W 2,p (Ω) with p ∈ (2, 6 3-2-ε ), by the estimate (5.6) of [4] for a fixed p ∈ (2, 6 3-2-ε ), we deduce that w -L h w 1,Ω h w W 2,p (Ω) h w 2+ε,Ω .
This estimate and [START_REF] Maz | Elliptic boundary value problems on manifolds with singularities[END_REF] prove the estimate [START_REF] Kozlov | Elliptic boundary value problems in domains with point singularities[END_REF].

We observe that I h u = L h w = L h u F on the face F , since ũI and D h (ũ -ũI ) vanish on F . We next state and prove an error estimate for the Finite Element approximation on the face F . Theorem 3.2 For a macro element Λ such that Λ ∩ F = ∅, let ω F,v, be the interior angle of F associated with the vertex v ∈ V. Assume that the conditions

ν < π ω F,v, , (20a)
1 ν + 1 µ ( π ω F,v, -1) > 1, (20b)
are satisfied. Then, it holds:

(21) u F -L h u F 1,F h.
Proof. We will prove that for all = 1, • • • , L, we have

(22) |u F -L h u F | 1,F ∩ Λ h.
Hence, summing on , we find that

|u F -L h u F | 1,F h,
and the conclusion of the theorem follows from Poincaré's inequality.

To prove [START_REF] Nicaise | Exact controllability of a pluridimensional coupled problem[END_REF], we distinguish different cases:

1. F ∩ Λ contains no singular vertex or singular edge: In this case, u F belongs to H 2 (F ∩ Λ ) and the mesh on F ∩ Λ is quasi-uniform. Thus, the estimate ( 22) is standard.

2. F ∩ Λ contains a singular vertex v but no singular edge: Thus, u F belongs to

V 2 γ (F ∩ Λ ) = {v ∈ L 2 loc (F ∩ Λ ) : R γ+|β|-2 D β v ∈ L 2 (F ∩ Λ ), ∀|β| ≤ 2}, with γ > 1 -π ω F,v,
, and the estimate ( 22) is also standard, since the triangulation in F ∩ Λ is isotropic (see for instance [START_REF] Raugel | Résolution numérique par une méthode d'éléments finis du problème de Dirichlet pour le laplacien dans un polygone[END_REF], [15, §8.4]).

3. F ∩ Λ contains a singular edge: Then, the mesh on F ∩ Λ is anisotropic. There are two possible situations: (S1) F ∩ Λ contains no singular vertex; and (S2) F ∩ Λ also contains a singular vertex v. Due to Corollary 2.3, for (S1), u F belongs to H 2 (F ∩ Λ ), while for (S2), u F belongs to V 2 γ (F ∩ Λ ). Now for any triangle T i in F ∩ Λ , we will prove that

(23) |u F -L h u F | 1,Ti h|u F | 2,γ,Ti ,
with γ = 0 for (S1) and γ > 1 -π ω F,v, for (S2), where

|u| 2 2,γ,T = |α|=2 T R 2γ |D α u| 2 dx.
If this estimate is valid, then summing on T i , we get [START_REF] Nicaise | Exact controllability of a pluridimensional coupled problem[END_REF].

To prove [START_REF] Nicaise | Polynomial stabilization of the wave equation with Ventcel's boundary conditions[END_REF], we distinguish two cases.

i. If T i is far from the singular corner, then we know that u F belongs to H 2 (T i ), and, by using Estimate (25) below, we have:

|u F -L h u F | 1,Ti h i |∂ 1 u F | 1,Ti + H i |∂ 3 u F | 1,Ti (24)
H i |u F | 2,Ti .
If Λ is of Type 3, then u F belongs to H 2 (F ∩ Λ ), but as ν = 1, by the assumptions on the mesh we have H i h, hence the estimate (24) directly yields [START_REF] Nicaise | Polynomial stabilization of the wave equation with Ventcel's boundary conditions[END_REF]. If Λ is of Type 4, we again distinguish two cases:

a) If R i h 1 ν , then H i ∼ hR 1-ν i , and therefore, |u F -L h u F | 1,Ti hR 1-ν i |u F | 2,Ti .
This yields [START_REF] Nicaise | Polynomial stabilization of the wave equation with Ventcel's boundary conditions[END_REF] by our assumption (20a). b) If 0 < R i h 1 ν , then as H i ∼ h 1 ν , the estimate (24) becomes:

|u F -L h u F | 1,Ti h 1 ν |u F | 2,Ti .
But from Lemma 4.5 below, we know that R i h 1 µ and, therefore,

|u F -L h u F | 1,Ti h 1 ν -γ µ R γ i |u F | 2,Ti .
This yields [START_REF] Nicaise | Polynomial stabilization of the wave equation with Ventcel's boundary conditions[END_REF] by our assumption (20b).

ii. If T i is near a singular corner (i.e., R i = 0), then applying Lemma 4.3 we have:

|u F -L h u F | 1,Ti h 1 ν -γ µ |u F | 2,γ,Ti .
Again we get [START_REF] Nicaise | Polynomial stabilization of the wave equation with Ventcel's boundary conditions[END_REF] owing to our assumption (20b).

The proof is now complete. Theorems 3.1 and 3.2 directly lead to the following a priori global interpolation estimate on u and error estimate on the Finite Element solution u h . Corollary 3.3 Assume that for all = 1, • • • , L, (16a), (16b), (16c), (20a) and (20b) hold. Then there exists

v h ∈ V h such that u -v h V h.
Corollary 3.4 Under the assumption of Corollary 3.3, if u ∈ V is the solution of (2) and u h ∈ V h the solution of (15), then u -u h V h.

Anisotropic error estimates in two dimension

To complete the proof of Theorem 3.2 we need some interpolation estimates in two space dimensions. In this section, T i will be a triangle in the triangulation F h of the face F , which is induced by the triangulation T h of Ω. We will need the two-dimensional version of Theorem 4.10 of [START_REF] Apel | The finite element method with anisotropic mesh grading for elliptic problems in domains with corners and edges[END_REF], given below.

Theorem 4.1 Assume that Λ is of Type 3 or 4. Suppose that F ∩ Λ contains the singular edge.

Recall the local Cartesian coordinate system (x 1 , x 2 , x 3 ) for Λ , for which the singular edge is on the x 3 -axis. Let F ∩ Λ be in the plane given by x 2 = 0. Let T i ⊂ F ∩ Λ be a triangle in the triangulation F h . Then, for v ∈ H 2 (T i ), we have

(25) |v -L h v| 1,Ti h i |∂ 1 v| 1,Ti + H i |∂ 3 v| 1,Ti ,
where h i and H i are defined in [START_REF] Grisvard | Théorèmes de traces relatifs à un polyèdre[END_REF].

Proof. Let h1,i and h3,i be the lengths of the projections of T i on the x 1 -and x 3 -axis, respectively. We distinguish between the case h3,i h i or not. This estimate implies [START_REF] Raugel | Résolution numérique par une méthode d'éléments finis du problème de Dirichlet pour le laplacien dans un polygone[END_REF], because hi,1 h i and h3,i H i (see assumption (B1) in [START_REF] Apel | Anisotropic interpolation with applications to the finite element method[END_REF], recalling that µ < 1 if a macro element is of Type 3 or 4).

With this choice, we see that T (∂ 1 (u -q))(x) dx = 0, and therefore by [START_REF] Ventcel | On boundary conditions for multi-dimensional diffusion processes[END_REF] we obtain:

(30) ∂ 1 (u -q) 1,γ |∂ 1 (u -q)| 1,γ = |∂ 1 u| 1,γ .
As q -Lu is linear, ∂ 1 (q -Lu) is constant, and we can write

∂ 1 (q -Lu) 1,γ |F (∂ 1 (q -Lu))| = |F (∂ 1 (q -u))|.
By (29), we deduce that

∂ 1 (q -Lu) 1,γ ∂ 1 (q -u) 1,γ .
This estimate and the triangle inequality imply that

∂ 1 (u -Lu) 0, T ≤ ∂ 1 (u -Lu) 1,γ ≤ ∂ 1 (u -q) 1,γ + ∂ 1 (q -Lu) 1,γ ∂ 1 (q -u) 1,γ
and the conclusion for ∂ 1 follows from (30). The estimate for ∂ 2 follows in an analogous manner.

Then, we are ready to derive the interpolation error estimate near a singular corner of F .

Lemma 4.3 Assume that Λ is of Type 3 or 4. Let 0 ≤ γ < 1. If T i is near a singular corner (i.e., R i = 0), then for any u F ∈ V 2 γ (T i ), we have |u F -L h u F | 1,Ti h 1 ν -γ µ |u F | 2,γ,Ti .
Proof. The result follows by mapping T i to T as in Lemma 4.8 of [START_REF] Apel | The finite element method with anisotropic mesh grading for elliptic problems in domains with corners and edges[END_REF], by using the estimate (26), and then mapping back to T i by using the properties (3.2) and (3.3) in [START_REF] Apel | The finite element method with anisotropic mesh grading for elliptic problems in domains with corners and edges[END_REF] and the fact that R h [4, p. 538]). Remark 4.4 If T i is isotropic, the previous Lemma is well known and can be found in [START_REF] Raugel | Résolution numérique par une méthode d'éléments finis du problème de Dirichlet pour le laplacien dans un polygone[END_REF] (see also [15, §8.4]). Lemma 4.5 Assume that Λ is of Type 4. Let T i be a triangle belonging to Λ ∩ F such that

-1 µ R i (see
R i > 0, then R i h 1 µ .
Proof. Without loss of generality, by a relabeling, we can always assume that T 0 is the triangle that contains the singular vertex v . Then, it has two edges that contain v , the first one is the edge in the x 1 -axis and is of length ∼ h 1 µ , while the other one is of length ∼ h 1 µ . Moreover, as the angle between these two edges is independent of the mesh, the ball of center v and radius ch

Numerical examples

In this section, we present some numerical examples to illustrate the theory presented in the previous sections.

We will solve the boundary-value problem (1) using the FEM with linear elements on a polyhedral domain. The domain is given as follows. We let T be the triangle with vertices (0, 0), (1, 0), and (0.5, 0.5), and let the domain be the prism Ω := (0, 1) 2 \ T × (0, 1). We refer to the labeling in Figure 1 in what follows. We will solve (1) in variational form (2) with data f = 1 and g = 0.

The interior angle between the two faces that contain the edge e := v 2 v 7 is 3π/2. Based on the estimates in (10) and Theorem 2.5, e is the singular edge; and the solution u admits a decomposition into the singular and regular parts with regularity determined by λ e = 2/3. By Theorem 2.5, the location of the face F , where the Ventcel boundary condition is imposed does not drastically affect the regularity of the solution.

To verify our theory, we implement two sets of numerical tests regarding different locations of the special boundary face F = Γ V : (I) F is the bottom face of prism Ω, with vertices v 1 , v 2 , v 3 , v 4 , and v 5 ; (II) F is a face that contains the singular edge with vertices v 2 , v 3 , v 7 , and v 8 .

For both cases, the singular parts u S of the solution have anisotropic exponents and belong to the same weighted space. Moreover, by Corollary 3.4, it is sufficient to choose the parameters in [START_REF] Grisvard | Théorèmes de traces relatifs à un polyèdre[END_REF] corresponding to the singular edge such that µ < 2/3 and ν = 1, in order to achieve the optimal (first-order) convergence rate.

In Table 1, we list the convergence rates of the numerical solution for the aforementioned model problems with ν = 1, but with different values of the mesh grading parameter µ . We let N be the number of degrees of freedom in the discrete system. Then, the mesh size satisfies h ∼ N -1/3 . Since the exact solution is not known, the convergence rate is computed using the numerical solutions for successive mesh refinements, u 2h u h , and u h/2 , as the convergence rate = log 2 (

u h -u 2h V u h/2 -u h V ), (31)
where u 2h and u h/2 are the finite element solutions with mesh parameters 2h and h/2, respectively. Therefore, as h decreases, the asymptotic convergence rate in (31) is a reasonable indicator of the actual convergence rate for the Finite Element solution.

It is clear from the table that for both cases, the first-order convergence rate is obtained for µ = 0.58 < 2/3, while we lose the optimal convergence rate if µ = 0.76, 1.00, both larger than the critical value 2/3. When µ = 0.76, that is, 2/3 < µ < mesh graded toward the singular edge, but the grading is insufficient to resolve the singularity in the solution, and hence does not give rise to the predicted first-order convergence rate. These results are in strong agreement with the theoretical results in Sections 3 and 4.

  otherwise. Similarly, we set λ ( ) e = λ e if Λ contains on singular edge e of Ω and λ ( ) e = ∞ otherwise. Then, the following decomposition for the variational solution of (1) holds.

1 µ

 1 intersects only T 0 by choosing c small enough. The estimate follows from the definition of the distance.

Figure 1 :

 1 Figure 1: The computational domain Ω (left) and the macro elements (right).

  vector, we will write v ∈ H s (Ω) d , although for ease of notation, we may write H s (Ω) simply for H s (Ω) d . Lastly, we employ the standard notation D (Ω)) to denote the space of distributions on Ω.Throughout, the notation A B is used for the estimate A ≤ C B, where C is a generic constant that does not depend on A and B. The notation A ∼ B means that both A B and B A hold. We will also employ standard multi-index notation for partial derivatives in R d , i.e., ∂ α = ∂ α1 x1 . . . ∂ α d x d where α = (α 1 , . . . , α d ) ∈ Z d + and |α| = α 1 + . . . + α d .

	Acknowledgements: The second author was partially supported by the US National Science
	Foundation (NSF) grant DMS-1418853. The third author was partially supported by NSF grant
	DMS-1312727. The visit of the first author to Penn State University and Wayne State University,
	where part of this work was conducted, was partially supported through NSF grant DMS-1312727
	and the Wayne State University Grants Plus Program.

1 .

 1 If h3,i h i , then diam T i h i (see the proof of Theorem 4.10 of [4]) and owing to Theorem 2 of [2], we have |v -L h v| 1,Ti h i |v| 2,Ti , and (25) holds since h i H i . 2. If h3,i h i , then Theorem 1 of [2] on the reference element and Lemma 4.8 of [4] yield |v -L h v| 1,Ti hi,1 |∂ 1 v| 1,Ti + h3,i |∂ 3 v| 1,Ti .

Table 1 :

 1 1, this choice still leads to an anisotropic Convergence rates for different values of µ : (I) F is the bottom face (left); (II) F is a side face containing the singular edge (right).

	h\µ	0.58	0.76	1.00	0.58	0.76	1.00
	2 -3 0.834 0.843 0.825	0.821 0.833 0.825
	2 -4 0.938 0.930 0.890	0.936 0.896 0.889
	2 -5 0.977 0.960 0.894	0.977 0.899 0.890
	2 -6 0.991 0.968 0.871	0.990 0.876 0.866
	2 -7 0.995 0.968 0.837	1.000 0.842 0.831

We continue with an anisotropic error estimate in weighted Sobolev spaces (compare with Theorem 1 of [START_REF] Apel | Anisotropic interpolation with applications to the finite element method[END_REF] for two-dimensional triangles in standard Sobolev spaces and Theorem 4.5 of [START_REF] Apel | The finite element method with anisotropic mesh grading for elliptic problems in domains with corners and edges[END_REF] for three-dimensional tetrahedra in weighted Sobolev spaces). Theorem 4.2 Let T be the standard reference element of vertices (0, 0), (1, 0) and (0, 1). Denote by R the distance to (0, 0). Let 0 ≤ γ < 1. Then for all u ∈ V 2 γ ( T ), and i = 1 or 2, we have:

where Lu is the Lagrange interpolant of u.

Proof. We first remark that Lemma 8.4.1.2 of [START_REF] Grisvard | Elliptic problems in nonsmooth domains[END_REF] shows that V 2 γ ( T ) is continuously embedded into C( T ), hence the Lagrange interpolant Lu of u is well-defined. We define the space:

which is an Hilbert space equipped with its natural norm • 1,γ . We will also use the semi-norm:

Then by the proof of Lemma 8.4.1.2 of [START_REF] Grisvard | Elliptic problems in nonsmooth domains[END_REF], we know that H 1,γ ( T ) is embedded into W 1,p ( T ) for all 1 < p < 2 1+γ , and hence compactly embedded into L 2 ( T ). The first embedding and a trace theorem also guarantee that any v ∈ H 1,γ ( T ) satisfies [START_REF] Ventcel | Semigroups of operators that correspond to a generalized differential operator of second order[END_REF] v

for any edge ê of T . The second embedding implies that

for all v ∈ H 1,γ ( T ) such that T v dx = 0. Now we follow the arguments of Lemma 3 and Theorem 1 of [START_REF] Apel | Anisotropic interpolation with applications to the finite element method[END_REF]. We will first prove the estimate for ∂ 1 . We observe that [START_REF] Ventcel | Semigroups of operators that correspond to a generalized differential operator of second order[END_REF] implies that the functional

where ê1 is the edge of T parallel to the x1 axis, is well defined and continuous on H 1,γ ( T ):

Next, we note note that F (∂ 1 (u -Lu)) = 0, if u ∈ V 2 γ ( T ). We then define the polynomial q of degree 1 by

where c = 2 T (∂ 1 u)(x) dx.