FETI-DP for the three-dimensional Virtual Element Method
Silvia Bertoluzza, Micol Pennacchio, Daniele Prada
TL;DR
This paper extends the FETI-DP preconditioner analysis to 3D virtual element methods, proving polylogarithmic condition number bounds that are independent of subdomain count, mesh size, and coefficient jumps, with numerical validation.
Contribution
It generalizes FETI-DP preconditioning results to three-dimensional VEM, establishing bounds independent of key problem parameters.
Findings
Polylogarithmic condition number bounds proven for 3D VEM.
Bounds are independent of subdomain number, mesh size, and coefficient jumps.
Numerical experiments confirm theoretical results.
Abstract
We deal with the Finite Element Tearing and Interconnecting Dual Primal (FETI-DP) preconditioner for elliptic problems discretized by the virtual element method (VEM). We extend the result of [22] to the three dimensional case. We prove polylogarithmic condition number bounds, independent of the number of subdomains, the mesh size, and jumps in the diffusion coefficients. Numerical experiments validate the theory
| Mesh | Mesh | ||||||
|---|---|---|---|---|---|---|---|
| oct1 | voro1 | ||||||
| oct2 | voro2 | ||||||
| oct3 | voro3 | ||||||
| oct4 | voro4 | ||||||
| oct5 | voro5 | ||||||
| oct6 | voro6 | ||||||
| oct7 | voro7 | ||||||
| oct8 | voro8 |
| Octahedra | Voronoi | Coarse | |||||||
|---|---|---|---|---|---|---|---|---|---|
| d.o.f. | d.o.f. | V | E | F | VE | VF | |||
| V | E | F | VE | VF | ||||||
|---|---|---|---|---|---|---|---|---|---|---|
| it | it | it | it | it | ||||||
| V | E | F | VE | VF | ||||||
|---|---|---|---|---|---|---|---|---|---|---|
| it | it | it | it | it | ||||||
| V | E | F | VE | VF | ||||||
|---|---|---|---|---|---|---|---|---|---|---|
| it | it | it | it | it | ||||||
| V | E | F | VE | VF | ||||||
|---|---|---|---|---|---|---|---|---|---|---|
| it | it | it | it | it | ||||||
| Mesh | Npol | |||
|---|---|---|---|---|
| voro-b1 | ||||
| voro-b2 | ||||
| voro-b3 | ||||
| voro-b4 |
| d.o.f. | it | d.o.f. | it | d.o.f. | it | ||||
|---|---|---|---|---|---|---|---|---|---|
| d.o.f. | it | d.o.f. | it | d.o.f. | it | ||||
|---|---|---|---|---|---|---|---|---|---|
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
\headers
FETI-DP for 3D VEMS. Bertoluzza, M. Pennacchio, and D. Prada
FETI-DP for the three-dimensional Virtual Element Method††thanks:
Submitted December 17, 2018, revised February 27, 2020. \fundingThe work of the authors was realized in the framework of ERC Project CHANGE, which received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program grant 694515.
Silvia Bertoluzza
Istituto di Matematica Applicata e Tecnologie Informatiche “E. Magenes”, IMATI-CNR, Pavia, PV, 27100, Italy
(, ,
[email protected]
Micol Pennacchio22footnotemark: 2
Daniele Prada22footnotemark: 2
Abstract
We deal with the finite element tearing and interconnecting dual primal (FETI-DP) preconditioner for elliptic problems discretized by the virtual element method (VEM). We extend the result of [22] to the three dimensional case. We prove polylogarithmic condition number bounds, independent of the number of subdomains, the mesh size, and jumps in the diffusion coefficients. Numerical experiments validate the theory.
keywords:
Virtual Element Method, Domain decomposition methods, Substructuring preconditioners
{textblock*}
14cm(3.5cm,24.8cm) This is the final draft of the paper published in SIAM J. Numer. Anal. Vol. 58, No. 3, pp. 1556-1591 (2020), doi: 10.1137/18M1233303. The original publication is available at https://epubs.siam.org.
{AMS}
65N30, 65N55
1 Introduction
Methods for the solution of PDEs based on polytopal meshes have recently attracted an increasing attention, mainly due to the necessity of tackling what is nowadays a bottleneck in the overall process of simulating real life phenomena, namely the task of mesh generation. Several methods have been recently introduced which allow for quite general polygonal or polyhedral elements, such as mimetic finite differences [13, 26], discontinuous Galerkin-finite element method (DG-FEM) [2, 30], hybridizable and hybrid high-order methods [33, 35], weak galerkin method [59], BEM-based FEM [53] and polygonal FEM [55] to name a few.
Here we deal with the virtual element method (VEM) [6], a discretization technique that can be considered as an extension of the FEM to polytopal tessellations. In such a method, local approximation spaces containing polynomial functions are defined and assembled in a global conforming approximation space, but the explicit construction and integration of the associated shape functions are avoided, whence the name virtual [6]. The evaluation of the operators and matrices needed in the implementation of the method is carried out by relying only on an implicit knowledge of the local shape functions, as described in [8] (see also [5, 11, 48], where the and versions of the method are discussed and analyzed). Though introduced fairly recently, such a method has already been applied and extended to a wide variety of different model problems; we recall applications to: parabolic problems [58], Cahn-Hilliard, Stokes, Navier-Stokes and Helmholtz equations [3, 4, 16, 17, 51], linear and nonlinear elasticity problems [14, 7, 39], general elliptic problems in mixed form [9], fracture networks [18], Laplace-Beltrami equation [38].
In this paper, we focus on the linear system of equations associated with the VEM discretization. As happens in the case of finite elements, the efficient solution of such a linear system is of paramount importance to fully exploit the potential of the method. Little work has been done on this issue up to now, all limited to the spatial dimension two. First works in the literature tackled the increase of the condition number appearing already at the level of the elementary stiffness matrix, due either to a degradation of the quality of the tessellation and/or to the increase in the polynomial order of the method [11, 34, 47, 19]. If we consider rather the increase of the condition number resulting from refining the discretization, to the best of our knowledge the approaches considered up to now are domain decomposition ([28, 27, 22, 52]) and multigrid ([5], for refinement). In the present paper we extend to the three dimensional case the results obtained in [22]. More precisely, we focus on one of the most efficient preconditioning techniques: the dual-primal finite element tearing and interconnecting (FETI-DP) [37, 56], a non overlapping domain decomposition method where the problem is reformulated as a constrained optimization problem and solved by iterating on the set of Lagrange multipliers representing the fluxes across the interface between the non overlapping subdomains. The FETI-DP method has been already extensively studied in the context of many different discretization methods – spectral elements [49, 41], mortar discretizations [40], NURBS discretizations in isogeometric analysis [50].
Following the approach presented in [22] for two dimensional domains, which mainly relies on the properties of the trace of the discrete space on the interface of the domain decomposition, we prove that the properties of scalability, quasi-optimality and independence on the discontinuities of the elliptic operator coefficients across subdomain interfaces, that are known for the finite element case, still hold when dealing with VEM. More specifically, we show that the condition number of the preconditioned matrix is bounded by a constant times the factor where and are, respectively, mesh-size of the subdomain decomposition and of the tessellation, see Theorem 4.6. In order to do so, we need to prove several inequalities related to the VEM approximation space, by only relying on the implicit definition of the discrete functions, which, we recall, are not explicitly known.
We observe that, since we are in the framework of [46], the equivalence of the BDDC (balancing domain decomposition by constraint) and the FETI-DP preconditioners holds. Therefore the bound for the condition number obtained here also yields an estimate on the BDDC preconditioner for VEM.
The paper is organized as follows. The basic notation, functional setting and the description of the VEM are given in Section 2. The dual-primal preconditioner is introduced and analyzed in Section 4, whereas some relevant properties of the virtual element (VE) discretization space, mainly used for the proof, are presented in Section 3. The analysis of the preconditioner, with the proof of the estimate for the condition number (Theorem 4.6), is carried out in Section 5 where we also give some detail specific to its implementation in the VEM framework. Numerical experiments that validate the theory are presented in Section 6.
2 The VEM
We start by recalling the definition and the main properties of the virtual element method [1]. To fix the ideas, we focus on the following elliptic model problem:
[TABLE]
with , where is (for simplicity) a convex polyhedron. We assume that the coefficient is a scalar function of such that for almost all , for two constants .
The tessellation
The virtual element method looks for an approximation to the solution of (2.1) in a conforming subspace of constructed on a polyhedral tessellation of . Let us then start by introducing the assumptions on the tessellation. We consider a family of tessellations of into a finite number of polyhedra , which we assume to be shape regular according to the following definition (quite standard in the theoretical study of VEM).
Definition 2.1*.*
We say that a polyhedron is shape regular of diameter with constants and if satisfies the following assumptions ([44]):
has at most faces and edges; 2. 2.
for every face and every edge we have:
[TABLE] 3. 3.
for each face , there exists a point such that is star-shaped with respect to every point in the disk of radius centered at ; 4. 4.
there exists a point such that is star-shaped with respect to every point in the sphere of radius centered at ; 5. 5.
for every face of , there exists a pyramid contained in such that its base equals to , its height equals and the projection of its vertex onto is .
We recall ([44]) that is shape regular if and only if it admits a conformal decomposition that is made of less than shape regular tetrahedra and includes all its vertices.
We assume the tessellation to be quasi-uniform and uniformly shape regular. More precisely we make the following assumption.
Assumption 2.1**.**
We assume that there exist two constants and such that the tessellation verifies the following assumptions:
* is geometrically conforming, that is for each , if a vertex, edge or face of is contained in , it is also, respectively, a vertex, edge, or face of ;* 2. 2.
all are shape regular of diameter with constants and ; 3. 3.
the tessellation is quasi uniform, that is there exists an such that for all , .
As we are interested here in explicitly studying the dependence of the estimates that we are going to prove on the number and size of the subdomains and the number and size of the elements of the tessellations, throughout the paper we will employ the notation (resp. ) to say that the quantity is bounded from above (resp. from below) by , with a constant independent of the diffusion coefficient in the PDE (and in particular, independent of , , and of their jump across the interface of the decomposition), that can depend on the polynomial order of the VEM method, and depending on the tessellation and on the domain decomposition only via the constants in Assumptions 2.1 and 4.1. The expression will stand for .
The local face space
The order virtual element discretization space on is defined element by element starting from the edges of the tessellation, where the discrete functions are defined as degree polynomials, then subsequently defining it on the faces and in the interior of the polyhedra. As the space on the faces, and the degrees of freedom identifying the elements of such a space, will play a key role in the definition and analysis of the FETI preconditioner, it is worth recalling its definition in some detail.
On the boundary of each face we introduce the space:
[TABLE]
where, for any one, two or three dimensional domain , denotes the set of order polynomials on . In order to define the discrete face space , being a face of one of the polyhedra of the tessellation, we start by introducing an auxiliary space , defined as
[TABLE]
It is known [6] that an element in is uniquely identified once its trace on and its moments up to order (or, equivalently, the scalar products with the elements of a basis for the space ) are known.
The local face space is a suitable subspace of . In order to define it, we introduce the projection , orthogonal to the scalar product
[TABLE]
To define we choose two bases for the two embedded subspaces of polynomials of order less than or equal to, respectively, and (where, for we write conventionally ):
[TABLE]
where, for , denotes the dimension of the space . We assume that is a suitably scaled Riesz basis. More precisely, we make the following assumption.
Assumption 2.2**.**
For all polynomial it holds that
[TABLE]
We then define as
[TABLE]
We immediately see that, as is a projector, for all polynomials , , and hence, as and , it holds that .
It is also easy to see that a function in is uniquely determined by
- .1)
the values of at the vertices of (vertex degrees of freedom); 2. .2)
the values of at the interior nodes of the points Gauss-Lobatto integration rule (edge degrees of freedom); 3. .3)
the value of the scalar products (face degrees of freedom)
[TABLE]
In fact, the degrees of freedom .1)-.2) uniquely determine the trace of on , which is an element of . Moreover, the values of .1)–.3) uniquely determine , as they are sufficient to compute for all . Indeed [6], for , Green’s formula yields
[TABLE]
where is the outer unit normal to . The values of .1)–.3) give access to on and to its scalar product with any polynomial of order up to , thus allowing us to compute the right hand side of (2.6), as . In view of the definition of , this implies, among other things, that the knowledge of the values of .1)–.3) allows us to compute exactly the scalar product of a function in times any polynomial of order less than or equal to .
Observe that the space does not depend on the choice of the basis functions in , but it does depend on the choice of the basis functions in . Usually, the , , are chosen to be the scaled monomials of order greater than and less than or equal to , though other choices are possible (and even advisable, for ).
Remark 2.2*.*
Requiring that Assumption 2.2 holds is equivalent to requiring that for all vectors it holds that
[TABLE]
(see [32]). It is a quite natural assumption, satisfied by many possible choices of . Of course, (2.4) holds for any suitably normalized orthogonal basis for . In [31], the authors prove that it is indeed satisfied for the basis of rescaled monomials , , , where we set
[TABLE]
A similar argument shows that such a property is satisfied by the rescaled version of any basis for
[TABLE]
Remark 2.3*.*
For the sake of simplicity, the basis functions are assumed to be normalized in such a way that (2.4) holds rather than the maybe more natural
[TABLE]
(that would hold, for instance, in the case of orthonormal bases). Remark, however, that the normalization of the degrees of freedom is transparent to the action of the preconditioner, so that both the design of the preconditioner itself and the theoretical bounds on the resulting condition number turn out not to depend on it.
The element space and the global space
In order to build the discrete space on a polyhedron , we assemble the local face spaces to build a local space defined on the boundary of :
[TABLE]
The local space on is finally defined as
[TABLE]
As contains the subspace of continuous piecewise polynomials of order up to , it is immediate to see that contains the polynomials of order up to .
A function in is uniquely determined by the value of the vertex, edge and face degrees of freedom .1) – .3) corresponding to all the faces of , plus
- ()
the values of the moments of in , up to order (interior degrees of freedom).
Such degrees of freedom are sufficient to compute the action of the projection , orthogonal to the scalar product
[TABLE]
In fact, for , integrating by parts we obtain
[TABLE]
where is the outer unit normal to . The right hand side can be computed since, on the one hand, is in and, on the other hand, since is in for all faces , as observed in the previous section, the degrees of freedom – enable us to compute the second integral on the right hand side.
We finally define the global virtual element space on by gluing the local spaces:
[TABLE]
Remark 2.4*.*
The definition of can appear quite cumbersome, and one might think that there are simpler options for defining the space . Two possibilities come to mind. The first one is to define with
[TABLE]
as is done in the definition of the simplest form of the two dimensional virtual element method. However, if we define in this way, the knowledge of the values of the degrees of freedom .1)–.3) would not be sufficient to compute the scalar product of a function in with an order polynomial. This is needed to evaluate the projector , which, as we will see, is necessary to build the discretization of the PDE that we want to solve. The second possibility would be to let . This choice would indeed allow the computation of but it would lead to a dramatic increase in the number of degrees of freedom (an extra per face), without any relevant increase in the order of approximation.
The global degrees of freedom
It is immediate to see that a function in is uniquely determined by the values of the degrees of freedom .1)–.3) and K) for all vertices, edges, faces, and polyhedra of the tessellation. Requiring the continuity of reduces to asking that all .1), .2) and .3) are single valued, that is, that they produce the same result when evaluated on the restriction of to two polyhedra sharing, respectively, a vertex, edge or face.
It will be convenient, in the following, to explicitly introduce functionals , with denoting the dimension of the global space , mapping the elements of to the value of the corresponding -th degree of freedom. In other words, a function is uniquely determined by the vector where is either one of the following:
is the value of at one of the vertices of the tessellation (vertex degrees of freedom); 2. 2.
is the value of at one of the interior nodes of the points Gauss-Lobatto quadrature formula (edge degrees of freedom); 3. 3.
is the scalar product , , with being one of the faces of the tessellation (face degrees of freedom); 4. 4.
is one of the moments of order up to of over one of the elements (interior degrees of freedom).
Note that edges and vertex degrees of freedom are “nodal”, that is they correspond to a node, that we will denote .
We let denote the set of all the degrees of freedom of ,
[TABLE]
We recall that for we can define its support as follows [36]:
[TABLE]
where the annihilator is the largest open subset of where vanishes (where, by definition, a distribution vanishes on an open subset of , if and only if for all ). It is easy to see that the support of a vertex or edge degree of freedom is the corresponding node, whereas the support of a face degree of freedom is the closure of the corresponding face, and the support of an interior degree of freedom is the closure of the corresponding element. Observe that
[TABLE]
Remark 2.5*.*
Also for the interior degrees of freedom ) we could consider a basis for the space different from the monomial basis. For increasing this can help to increase the stability of the method (and in particular it allows better constants in the stabilization bound (2.12)). Remark, however, that the choice of such degrees of freedom does not directly enter in the definition of the preconditioner, as they are locally eliminated as a preliminary step, and it affects the theoretical analysis only through such constants.
Remark 2.6*.*
It is not difficult to realize that the value of a function at any vertex, edge or face of the tessellation is uniquely determined by the values for supported, respectively, in the vertex, in the closure of the edge and in the closure of the face. This is trivial for the values at the vertices, but is easily checked also for the value at the edges (a polynomial of degree is, of course uniquely determined by its value at nodes of the point Gauss-Lobatto integration rule) and on the faces (uniquely determined by the values of .1)–.3), as observed previously when introducing the local face space). Consequently, if is obtained as the union of some vertices, edges and/or faces, is uniquely determined by those degrees of freedom such that .
The VE Discretization
The Virtual Element Method stems from a Galerkin discretization of problem (2.1), starting from its variational formulation, which reads
[TABLE]
with
[TABLE]
As the functions in are not known in closed form, it is not possible to directly evaluate the bilinear form on two of such functions (this would imply solving Poisson equations on each face and in each element). Clearly we have
[TABLE]
with local counterpart of . The virtual element method is obtained by replacing the operator in the second term on the right hand side, which cannot be computed exactly, with an “equivalent” computable operator . We then define
[TABLE]
where is any continuous symmetric bilinear form satisfying
[TABLE]
Equations (2.11) and (2.12) immediately yield
[TABLE]
Finally, we let be defined by
[TABLE]
and we consider the following discrete problem:
Problem 2.7*.*
Find such that
[TABLE]
where we let denote the projection of onto the space of polynomials of order less than or equal to .
Different choices of the stabilization bilinear forms are proposed in the literature. The simplest, widely used choice, is directly expressed in terms of the degrees of freedom as the suitably scaled euclidean scalar product [12]
[TABLE]
An alternative which gives sharper bounds in (2.12) as increases is the so called “diagonal recipe”:
[TABLE]
where is the canonical basis function corresponding to the degree of freedom , that is the unique function in such that , denoting the Kronecker delta. It has also been proposed to replace the sum on the right hand side of (2.14) or (2.15) by a sum over those such that . Other recipes (proposed in two dimensions, but also valid in three dimensions), include suitably scaled versions of the scalar product, and of the bilinear form corresponding to the Laplace-Beltrami operator on [15].
For the study of the convergence, stability and robustness properties of the method we refer to [6, 10, 24].
Remark 2.8*.*
Different alternatives have been proposed for defining the constant components and entering the definition of the scalar products and and, consequently, of the projectors and . In particular, for we can consider
[TABLE]
while for we can set
[TABLE]
Remark that the choice of does not have any effect on the matrix resulting from discretizing the bilinear form, as, in its definition, the constant component of and are canceled by the operator. On the contrary, the choice of does have a direct influence on the definition of the projection of onto , which enters the definition of through the second term on the right hand side of (2.9).
3 Some relevant bounds
In this section we present some bounds that will play a role in the forthcoming analysis. More precisely, letting be a shape regular polyhedron of diameter , and letting be any face of , we have the following bounds.
Agmon inequality
For all functions in , it holds that [15, 44]
[TABLE]
Inverse estimates
The following two inverse inequalities hold for all [15, 29, 31, 57]:
[TABLE]
Riesz basis property
The following lemma, related to the equivalence between the norm of a function in and the euclidean norm of the vector of its degrees of freedom, holds.
Lemma 3.1*.*
Let be a face of the tessellation. For all we have that
[TABLE]
Bound (3.4) has been proven in [31] with a different choice of the edge degrees of freedom (namely, the moments of order up to ) and for being the basis of the rescaled monomials.
It is not difficult to verify that the proof therein still holds, with some minor changes, for our choice of degrees of freedom, provided Assumption 2.2 holds, and provided that for we have
[TABLE]
where we recall that, for , denotes the node corresponding to the degree of freedom, which can be either a vertex of or one of the nodes of the points Gauss-Lobatto quadrature formula on an edge of . The norm equivalence (3.5) does indeed hold with our choice of degrees of freedom, and can be proven, using, edge by edge, a scaling argument, and the equivalence of the norm of a polynomial in and the norm of the vector of its values at the nodes of the points Gauss-Lobatto quadrature formula.
For the sake of completeness let us sketch the main steps of a proof of Lemma 3.1, which is slightly different from the one proposed in [31], yielding a sharper bound for suitable choices of the face degrees of freedom (such as, for instance, when is chosen to be an orthogonal basis for ).
Letting denote the orthogonal projection, and combining (2.4) with the definition of orthogonal projection, we can easily see that for all we have
[TABLE]
Then, proving that the right hand side of (3.4) is bounded by a constant times the left hand side is not difficult. In fact we have that
[TABLE]
We can bound the first term on the right hand side by combining (3.5), (3.1), and (3.3) and the second term by using (3.6) and the boundedness of the orthogonal projection.
To prove the converse inequality, as in the proof of [31], combining (3.5) with the maximum principle, it is not difficult to prove the following proposition.
Proposition 3.2*.*
Let verify and in . Then
[TABLE]
Additionally, we have the following proposition.
Proposition 3.3*.*
For all it holds that
[TABLE]
Proof 3.4*.*
Using the Poincaré inequality, integrating by parts, applying a Cauchy–Schwartz inequality, and the inverse bounds (3.3) and (3.2) we obtain
[TABLE]
We conclude by dividing both sides by .
Letting now denote the harmonic lifting of for all we then have
[TABLE]
where we used Proposition 3.3, (3.6), and Proposition 3.2. We immediately see that if is orthogonal to , , we have that
[TABLE]
Let now denote the function in orthogonal to for , and verify
[TABLE]
We easily see that, given any , we have that . Then, for we can write
[TABLE]
Now, as under our assumptions the operator is bounded in uniformly in [31], we have
[TABLE]
Bounding the right hand side by (3.8) and combining with (3.9) we get the thesis.
4 Domain decomposition for the Virtual Element Method
4.1 The subdomain decomposition
We assume that can be split as , inducing a decomposition of as the union of disjoint polyhedral subdomains :
[TABLE]
We make the following assumptions.
Assumption 4.1**.**
There exist a constant and such that the subdomain decomposition satisfies the following properties:
it is geometrically conforming, that is, for all , if a vertex, edge, or face of is contained in , it is also, respectively, a vertex, edge, or face of ; 2. 2.
*the subdomains are shape regular (in the sense of Definition 2.1) of diameter with constants and ; * 3. 3.
*for all , there exists a scalar such that ; * 4. 4.
*the decomposition is quasi uniform: there exists an such that for all we have . *
We will refer to the edges and faces of the subdomains as macro edges and macro faces. We let denote the skeleton of the decomposition, and denote, respectively, the set of macro edges and of macro faces of the subdomain decomposition interior to , and and denote the set of, respectively, macro faces and macro edges of the subdomain .
We also let denote the wirebasket (that is, the union of all the edges) of the decomposition. Finally, we let denote the set of the vertices of the decomposition (the cross points), see Figure 1.
We would like to point out that Assumption 4.1, which is quite standard in the framework of domain decomposition methods, is actually also an assumption on the tessellation , satisfied, for instance, if is built by first introducing the subdomains and then refining them. This is, of course, quite strong an assumption which might leave out some relevant discretization. In this respect, a couple of observations are in order. On the one hand, also in the finite element framework, the literature on FETI methods in irregular domains is quite scarce and, to our knowledge, limited to two dimensions [43]. On the other hand, while the theory relies on Assumption 4.1, the design of the preconditioner itself does not, so that in practice, the FETI-DP is also available for subdomains decomposition with irregularly shaped subdomains such as those obtained by graph partitioning algorithms starting from a fine tessellation.
Notation: Global versus local degrees of freedom
In the following we will need to single out different subsets of the set of degrees of freedom. To this aim we start by letting
[TABLE]
denote the set of indexes such that the degree of freedom is supported on the interface of the decomposition. For each subdomain we let
[TABLE]
denote the set of indexes of those degrees of freedom supported in .
Analogously, for each macro edge and for each macro face we let
[TABLE]
denote the set of indices of those degrees of freedom supported on and . Moreover, we let be the set of indices of those degrees of freedom that are, respectively, supported on the wirebasket ,
[TABLE]
and we let
[TABLE]
respectively, denote the set of indices of those degrees of freedom that are supported on the edges of and on .
Finally, we let and denote, respectively, the set of indices of those degrees of freedom supported in the set of cross points (vertices of the subdomain decomposition), respectively, in the set of vertices of the subdomain :
[TABLE]
Conversely, for each degree of freedom supported on the skeleton of the decomposition we let denote the set of indices of those subdomains whose boundary the support of lies on:
[TABLE]
Observe that, for we can, in a natural way, give meaning to the expression for all . Finally, for each vertex , edge , and face we can also define the set , , and of the indices of those subdomains that share, respectively, , and as a vertex, edge or face:
[TABLE]
Remark that for all we have .
Notation: Scaled norms and seminorms
In the following we will make use of suitably scaled norms for the Sobolev spaces defined on faces and edges of the subdomains. More precisely, letting denote any dimensional domain () we set
[TABLE]
For bounded domain, we will also consider the spaces () and of those functions whose extension by zero is in () and , respectively, which we will equip with the norm
[TABLE]
Domain decomposition and FETI-DP preconditioner
The subdomain spaces and bilinear forms are defined, as usual, as
[TABLE]
In view of (2.13) we immediately obtain that for all
[TABLE]
Solving Problem 2.7 is then reduced to finding minimizing
[TABLE]
subject to a continuity constraint across the interface.
The FETI-DP preconditioner is then constructed according to the same strategy used in the finite element case, which we recall, mainly to fix some notation. We let
[TABLE]
and
[TABLE]
Moreover, for each macro face and macro edge we let
[TABLE]
denote, respectively, the trace on and of . On we define a norm and a seminorm:
[TABLE]
As usual, we define local discrete lifting operators as
[TABLE]
The following proposition holds.
Proposition 4.1*.*
is well defined, and it verifies
[TABLE]
Proof 4.2*.*
We start by recalling that there exists a linear operator such that, for , , one has ([29])
[TABLE]
(note that we are using scaled norms, see (4.2)–(4.4)). Moreover is constructed in such a way that if one has on . In view of this result it is not difficult to construct an operator satisfying and
[TABLE]
can be for instance defined as applied to the harmonic lifting of ; (4.8) follows then from the stability of the harmonic lifting and of , exactly as in the two dimensional case [22]. We can then write
[TABLE]
After dividing both sides by , in view of (4.8) it is easy to conclude by applying a triangular inequality.
Let now denote the subset of functions which are single valued across :
[TABLE]
For we let , so that is split as
[TABLE]
We next define the bilinear form as
[TABLE]
The proof of the following proposition, where the factor in bounds (4.11) stems from using norms scaled as in (4.3), is trivial.
Proposition 4.3*.*
For all we have
[TABLE]
The solution of Problem 2.7 is then split as , where and are the solutions of the following two independent problems.
Problem 4.4*.*
Find such that for all
[TABLE]
Problem 4.5*.*
Find such that for all
[TABLE]
Exactly as in the finite element case, the design of different versions of the FETI-DP method will rely on the choice of a subspace of whose elements have some degree of continuity on , ensuring that the restriction to of the bilinear form is coercive (which is equivalent to asking that the seminorm is a norm on ). We recall that, while in two dimensions the space can be defined as the subspace of functions continuous at the vertices of the subdomains, it is known that this is not sufficient to get a quasi-optimal result in the three dimensional case, so that we also need to impose continuity of either edge or face averages (or both). While later in the paper we will analyze the different possible choices, for now we will only assume that is coercive on (or, equivalently, that is coercive on ), and that .
Following the approach of [25], we introduce the operators and defined, respectively, as
[TABLE]
and we let denote the natural injection operator. We observe that
[TABLE]
Problem 4.5 becomes
[TABLE]
As in [45, 42], we define a scalar product in terms of the degrees of freedom , as
[TABLE]
where, for , the scaling coefficient is defined as
[TABLE]
Next, we introduce the projection operator , orthogonal with respect to the scalar product :
[TABLE]
Following [25], we introduce the quotient space . We let denote its dual and we let denote the quotient mapping, defined as
[TABLE]
Observe that two elements and are representative of the same equivalence class if and only if they have the same jump across the interface: for all macro facs with , . We can then identify with the set of jumps of elements of . The quotient map can then be interpreted as the operator that maps an element of to its jump on the interface. Clearly
[TABLE]
where is defined as . Problem 2.7 is then equivalent to the following saddle point problem: find , solution to
[TABLE]
Using the first equation in (4.17) to express as a function of , we eliminate the former unknown and finally reduce the solution of Problem 4.5 to the solution of a problem in the unknown of the form
[TABLE]
Letting now be any right inverse of (for , is any element in such that ), we set
[TABLE]
where denotes the identity operator in the space . Recall that, as in the two dimensional case, the definition of is independent of the actual choice of the operator . Indeed implies and hence .
We finally let the FETI-DP preconditioner be defined by
[TABLE]
Let us now come to the choice of the space , or, equivalently, to the choice of the so called primal degrees of freedom (d.o.f.s), which are taken as single valued, thus directly imposing the corresponding degree of continuity across the interface, whereas continuity for the remaining dual (d.o.f.s) is imposed via Lagrange multipliers in . Exactly as in the finite element case there are several possibilities. Letting
[TABLE]
denote the subset of of traces of functions which, respectively, are continuous at cross-points, have the same average at all edges, and have same average at all faces, different choices for can be considered, resulting in different versions of the FETI-DP algorithms. More precisely we have
- E:
(primal d.o.f.s are the edge averages);
- F:
\widetilde{W}_{h}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(F)}} (primal d.o.f.s are the face averages);
- VE:
\widetilde{W}_{h}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(V)}}\cap{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(E)}} (primal d.o.f.s are the values at cross points and the edge averages);
- VF:
\widetilde{W}_{h}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(V)}}\cap{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(F)}} (primal d.o.f.s are the values at cross points and the face averages);
- EF:
\widetilde{W}_{h}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(E)}}\cap{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(F)}} (primal d.o.f.s are edge and face averages);
- VEF:
\widetilde{W}_{h}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(V)}}\cap{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(E)}}\cap{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(F)}} (primal d.o.f.s are the values at cross points and the face and edge averages).
4.2 The main theorem
We can now state the main result of this paper.
Theorem 4.6*.*
Letting denote the condition number of the matrix corresponding to the operator , depending on the choice of the space we have the following bounds:
Algorithms E/EF:
If \widetilde{W}_{h}\subseteq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(E)}} then
[TABLE]
Algorithm F:
If \widetilde{W}_{h}\subseteq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(F)}} then
[TABLE]
Algorithms VE/VEF:
If \widetilde{W}_{h}\subseteq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(V)}}\cap{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(E)}} then
[TABLE]
Algorithm VF:
If \widetilde{W}_{h}\subseteq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(V)}}\cap{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(F)}} then
[TABLE]
where , , and are constants depending on the diffusion coefficient that satisfy
[TABLE]
As it happens for the finite element case, all the above choices lead to a quasi-optimal (in terms of dependence on and ) preconditioning, whereas robustness with respect to the jumps in the coefficient is achieved for algorithms VE and VEF.
Remark 4.7*.*
The precise definition of the constants , and (which will be detailed in Section 5) is the same as in the finite element case [56], and it is quite technical. We would like to remark that the bounds (4.22) are often quite pessimistic, as we will see in Section 6.
Remark 4.8*.*
We made the assumption that the subdomains are shape regular (point 2 of Assumption 4.1) in order to have some of the inequalities needed in the proof of Theorem 4.6. More precisely, such an assumption yields the validity of the Poincaré inequality on the macro faces, of the trace inequalities on the subdomains and on the macro faces, and of the injection bound (5.3). While Poincaré and trace inequalities are valid for a much larger class of domains, including the ones obtained by applying graph partitioning to a fine tessellation, in such a case the constants are a priori dependent on . In order to extend the theoretical result to the irregular case, we would need to provide a partitioning recipe, resulting in a subdomain decomposition for which such bounds can be proven uniformly in . A relevant case in which this is indeed possible, is the case where the “roughness” of the subdomains derives from properties of the continuous problems, such as, for instance, when the subdomains are irregular in order to fit a (fixed) irregular boundary or interior interface. More in general, we believe this is still possible if the partitioning is performed in such a way that the rough edges and faces are “as straight as possible”, however, it is outside the scope of this paper to deal with such an issue, which we plan to address in the future.
Remark 4.9*.*
The assumption that the tessellation and the domain decomposition are quasi–uniform can be dropped out, provided that we replace the ratio with the worst possible instance of such a ratio, that is, . Conversely, the assumption that the elements of the tessellation are shape regular (and more precisely, that all the faces of the tessellation lying on the interface are shape regular) plays a key role in the proof, as it is needed for the bounds in Section 3 to hold.
4.3 Implementation of the FETI method in the VEM context
Contrary to what happens for finite elements, VE discrete functions are explicitly known only on the edges of the subdomains, where they are piecewise degree polynomials. On the faces and within the subdomains the VEM basis functions are not explicitly known, and all quantities needed for the implementation of the preconditioner (but also of the stiffness matrix and load vector) have to be retrieved in terms of the d.o.f.s .1), .2), .3), and K), by exploiting the definition of the spaces and (in the terminology of VEM, they must be computable). Extended details on the computability of the stiffness matrix and the right hand side can be found in [8]. Let us concentrate here on the quantities needed for implementing the FETI-DP preconditioner in the VEM context.
As in the finite element case, the FETI-DP method is carried out by, at first, assembling local systems , where the matrix is the Schur complement obtained from the local stiffness matrix by eliminating the interior degrees of freedom. The efficient implementation of the FETI-DP method relies on expressing the elements of in a suitable basis, singling the primal degrees of freedom out, so that imposing reduces to requiring that the corresponding entries in the unknown vector (the primal unknowns) are single valued. Equation (4.18) then becomes
[TABLE]
where has entries in , and where the matrix and the vector are obtained from the local matrices and load vectors by partial assembly in the primal unknowns. The operators involved in equation (4.18) and in the definition (4.20) of the preconditioner can be expressed in a block form, easily allowing us to reduce the evaluation of their action (in particular the action of ) to the solution of local problems, plus a small coarse problem involving the primal degrees of freedom. The linear system (4.23) is solved by a preconditioned conjugate gradient method, with preconditioner given by (4.20), whose action involves the suitably combined (see (4.20) and (4.19)) actions of , any right pseudoinverse of , and the algebraic realization of the projector , whose explicit expression is given by (5.15) in Section 5.
Implementing the change of basis from the canonical VEM basis to the new basis entails the need to evaluate the primal degrees of freedom for any given discrete functions, and to explicitly express them as a function of the degrees of freedom .1), .2) and .3). The primal degrees of freedom involved in the definition of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(V)}} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(E)}} (vertex values and edge averages) are easily computed, as discrete functions are explicitly known on the edges of the subdomains. Let us then consider the face averages, involved in the definition of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(F)}}. In the case , in order to compute the face average of functions we need to resort to the definition (2.5) of the space , by which we have . We recall that, thanks to (2.6), is computable in terms of the (known) value of on . For the computation is simpler, as the degrees of freedom .3) provide direct access to the integral over any face of times any polynomial of degree less than or equal to , and, in particular, times the constants. Observe that, in view of the FETI implementation, it might be advisable to include a constant function among the elements of the basis .
5 Proof of Theorem 4.6
We start by remarking that, by construction, we are in the framework of [46]. In particular we have the identity
[TABLE]
In fact, it is not difficult to see that for all we have that and that we have
[TABLE]
Then we have [46]
[TABLE]
In order to have a bound on the condition number, we then only need to bound in terms of .
To start, let us recall some functional inequalities that will be useful in the following [21, 20]. Let be a shape regular polygon (in the following will be a face of one of the subdomains). Then, for all , , we have, uniformly in , the following trace inequalities,
[TABLE]
On the other hand, we recall that for the two spaces and coincide, and the two corresponding norms are equivalent. However, the constant in the equivalence depends on and it explodes as converges to . More precisely, for all , , and for all it holds that
[TABLE]
once again uniformly in (recall that we are using scaled norms, as defined in Section 2, so that the bounds are uniform in ).
We also observe that, thanks to the inverse inequality (3.3) and to the scaling of the norms, by using a standard space interpolation technique it is not difficult to prove that for all with , and for all , (resp. ) it holds that
[TABLE]
An analogous bound holds for the norms in the spaces and (with the usual care when either or is equal to ), provided . In particular, in such cases we have
[TABLE]
The following proposition holds.
Proposition 5.1*.*
Let be a shape regular subdomain and let be a face of . Then for we have
[TABLE]
Proof 5.2*.*
Using inequalities (5.4) and (5.2) we can write, for arbitrary,
[TABLE]
where the last bound is obtained by choosing .
We now prove the following lemma, which is the equivalent, for the Virtual Element Method, of Lemma 5.6 of [54] and Lemma 4.3 of [23].
Lemma 5.3*.*
Let and let be defined by for all , for all . Then, for all faces of it holds that and
[TABLE]
Moreover, if is constant on then
[TABLE]
Proof 5.4*.*
We let be defined as
[TABLE]
It is easy to see that for all (degrees of freedom supported on ) implies on and, hence, . Let and denote the -projection onto and onto , respectively. Recall that for all , using (4.7) we have
[TABLE]
and, since implies that , we also have for all
[TABLE]
Let now be defined by for all , so that, on , . Remark that, thanks to Lemma 3.1, we have
[TABLE]
Consider now the operator obtained by first projecting onto and then setting the values at nodes on to zero. We will prove that the restriction of to is uniformly bounded for all , that is, that for all we have
[TABLE]
with a constant independent of . Then, for arbitrary, using (5.5) and (5.3), we can write
[TABLE]
which, by choosing , yields
[TABLE]
Observing that for we have
[TABLE]
we immediately get the thesis (the result for constant on is obtained by setting in (5.3)).
Let us then prove (5.9). We easily see that is bounded: for all ,
[TABLE]
On the other hand, observing that , using (5.8) we see that, for ,
[TABLE]
By adding and subtracting in the second term on the right hand side and using (5.6) and (5.7), we obtain, for ,
[TABLE]
This allows us to prove, by a standard argument, that is -bounded. In fact, letting denote the projection, for , using (3.3), adding and subtracting and then using an approximation bound, we have
[TABLE]
The bound (5.9) follows by a standard space interpolation argument.
Let us now consider the projector . We have the following lemma.
Lemma 5.5*.*
For all it holds that
[TABLE]
with
[TABLE]
where
[TABLE]
Proof 5.6*.*
It is not difficult to check that, for , we have
[TABLE]
and that these relations completely define , as it is uniquely determined by the value of the degrees of freedom , , that is, by the degrees of freedom supported on the interface of the domain decomposition.
Let now both and be split as the sum of the contributions of the degrees of freedom supported on the wirebasket, which we will denote by and , respectively, and the contribution of remaining degrees of freedom, which we will denote by and , respectively. More precisely, we let and be defined by
[TABLE]
and we set and . Remark that
[TABLE]
To start, let us consider the contribution of the faces. We have
[TABLE]
We recall that for and we have . Let be a common face of the subdomains and . On we have , where . Let denote the function whose degrees of freedom assume the following values:
[TABLE]
where stands for the action of the functional on the constant unit function. We can write
[TABLE]
We now apply Lemma 5.3, which, thanks to the Poincaré inequality, gives us
[TABLE]
Adding up over all subdomains and over all faces (each face is counted twice) we obtain
[TABLE]
We now consider the contribution of the wirebasket. We recall that for ( is supported on the wirebasket ), the degree of freedom corresponds to a node and takes the form . Then, for all , (5.15) can be rewritten as
[TABLE]
Using (5.4) and (3.5) we can write
[TABLE]
and
[TABLE]
Plugging (5.17) into (5.16) and adding up over all we obtain
[TABLE]
Now, given , for and we can write
[TABLE]
yielding
[TABLE]
where we used once again (3.5) and the fact that, under the assumptions made on the tessellation, we have that . Then
[TABLE]
We conclude by observing that
[TABLE]
(we used that minimizes ). Applying Proposition 5.1 we then obtain
[TABLE]
Observe that , , and vanish provided that belongs to , , and , respectively, so that, depending on the choice of , some of the terms on the right hand side of (5.10) disappear. In order to get a bound for for the different choices of , we then need to bound the remaining terms in the different cases. We start by comparing, for a given function , the average over a face with the average over one of its edges.
Lemma 5.7*.*
For edge of , it holds for defined in (5.14):
[TABLE]
Proof 5.8*.*
Let denote the projection onto the space of constant functions, which is defined by
[TABLE]
Trivially, such an operator preserves the constants. Then, by using Proposition 5.1 and a Poincaré type inequality, allowing to bound the norm with the seminorm for the average free function , we can write
[TABLE]
We can now bound in terms of either or .
Lemma 5.9*.*
The following inequalities hold:
[TABLE]
with , constants depending on the diffusion coefficient , satisfying .
Proof 5.10*.*
We start by proving (5.21). Let denote one of the cross points, let (with ) denote the corresponding d.o.f., and let . Assume at first that and share an edge having as one of the vertices. Adding and subtracting , using Proposition 5.1 as well as a Poincaré–inequality for function with vanishing average in a portion of the boundary (allowing us to bound the norm with the seminorm), we can write
[TABLE]
Let now be two subdomains sharing the vertex but not an edge. In this case we bound by adding and subtracting a suitable sequence of values in such a way that we fall back into the previous case. To this aim we start by introducing the following definitions:
- •
a path of length is any sequence of subdomains such that for all , and share at least a vertex;
- •
for a given path we set , ;
- •
a path connects and via edges (resp. via faces) if , and for all the subdomains and share an edge (resp. a face).
Letting be the maximum number of subdomains sharing a vertex, we denote by (resp. ) the set of paths of length connecting and via edges (resp. via faces). For all paths we can bound
[TABLE]
and, using the bound for subdomains sharing an edge, we obtain (5.21) with
[TABLE]
where
[TABLE]
Bound (5.22) with and
[TABLE]
is obtained by a similar argument. As (if two subdomains share a face they also share a vertex) we easily get that .
Finally, we bound and in terms of each other.
Lemma 5.11*.*
The following bounds hold:
[TABLE]
with a constant that satisfies .
Proof 5.12*.*
Let us consider the first inequality. Let be a face common to the subdomains and , and let be any of the edges of . By Lemma 5.7 we have
[TABLE]
In view of the definition of and , the bound (5.23) follows up adding the contribution of all faces.
Let us now consider the second bound. Let be an edge, and let . Assume at first that and share a face . Then, as we did for the previous bound, it is not difficult to prove that
[TABLE]
Let us now assume that and do not share a face. Proceeding as in the proof of the previous lemma it is easy to see that bound (5.24) holds with
[TABLE]
As two subdomains that share an edge also share a vertex, we have that .
Now we have all the ingredients to prove Theorem 4.6. From (5.1) we get that to bound the condition number, we only need to bound in terms of . By using Lemma (5.21), (5.22), (5.23), and (5.24) we can easily obtain that
[TABLE]
As, by (4.11), we have that and , (5.25)–(5.27) yield the thesis.
Remark 5.13*.*
Observe that, since we are in the framework of [46], we have the equivalence of the BDDC preconditioner with the FETI-DP preconditioner. Therefore the analysis presented also yields an estimate on the BDDC preconditioner for the Virtual Element Method.
6 Numerical tests
We consider the model problem
[TABLE]
In all experiments is divided into cubic subdomains with side length , which are then each discretized by a rescaled version of a reference tessellation of the unit cube made of either truncated octahedra or Voronoi cells (see Figure 2). In order to guarantee conformity, the mappings from the different subdomains to the reference unit cube are defined in such a way that each macro face is a symmetry plane for the tessellations on the two sides of . Table 1 lists the values of the following geometrical parameters for the reference meshes used in the experiments:
[TABLE]
where is the diameter of element , is the minimum distance between any two vertices of and is the parameter given in Definition 2.1. From Table 1 we see that Voronoi meshes satisfy Assumption 2.1 but with worse constants than the octahedra ones. For all tests the stabilization bilinear form is the simplest one, defined in (2.14).
Concerning the space , the following choices from Section 5 are tested:
E: \widetilde{W}_{h}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(E)}}, F: \widetilde{W}_{h}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(F)}}, VE: \widetilde{W}_{h}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(V)}}\cap{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(E)}}, VF: \widetilde{W}_{h}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(V)}}\cap{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(F)}}.
For the sake of completeness, we also test V: \widetilde{W}_{h}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(V)}}, for which we expect a worse behaviour, as for the finite element case, where the condition number increases as (see [56], Remark 6.39).
In order to analyze the performance of these algorithms, we carry out two series of experiments:
Test 1:
FETI-DP scalability. We increase the number of subdomains and the overall problem size, while keeping the reference mesh (and, consequently, the ratio ), fixed. Table 2 shows the dimension of the primal spaces for the different algorithms. According to Theorem 4.6, we expect the condition number for the FETI-DP preconditioner to remain asymptotically constant.
Test 2:
FETI-DP quasi-optimality. We fix the number of subdomains , so that is kept constant, and increase the size of the local problems by choosing finer and finer reference meshes (from oct1 to oct8, or from voro1 to voro8 in Table 1), thereby incrementing the overall problem size. This results in a decreasing and, consequently, increasing ratio . According to Theorem 4.6, we now expect the condition number for the FETI-DP preconditioner to asymptotically exhibit a polylogarithmic behavior.
To test the robustness of FETI-DP, each test is run with two types of data (unless otherwise stated):
- •
in and ;
- •
constant in each subdomain, assuming values and in a (3D) checkerboard distribution, so that there is a coefficient jump of along . The right hand side is implicitly chosen by choosing a right hand side vector with values uniformly randomly distributed in .
Unless otherwise stated, we use the conjugate gradient with the FETI-DP Dirichlet preconditioner with zero initial guess and, as a stopping criterion, the relative reduction of the dual residual by either or when using the first or the second type of data, respectively. MATLAB R2016b is used as the subdomain and coarse sparse direct solver. If not otherwise specified, results are obtained on a machine equipped with a processor Intel Core i7-7820HQ, operating system Ubuntu Linux 16.04 LTS, memory 64 GB, 2400 MHz DDR4 Non-ECC SDRAM.
6.1 FETI-DP scalability
We test the scalability for both the constant coefficient case and the varying coefficients case on the lowest order method (). We fix the number of subdomains to be (). Results for the first series of experiments (Test 1) with constant coefficients are reported in Table 3 for meshes of truncated octahedra, and in Table 4 for Voronoi meshes. The results are in accordance with the theoretical bounds for both sets of meshes. Tables 3 and 4 show that, when no jumps in the coefficients are present, the results for E, F, VE, and VF are similar. In switching from the octahedra to the Voronoi mesh, which satisfies Assumption 2.1 but with worse constants, we observe an increase in the number of iterations and in the condition number which, however, still displays the expected behavior. As one expects, the choice V is, instead, not competitive.
The numerical tests with varying coefficients are displayed in Tables 5 and 6 for octahedra and Voronoi meshes, respectively. With highly oscillating coefficients, when only face averages are selected as primal degrees of freedom (Algorithm F), the preconditioner performs very poorly on both types of meshes.
This is not in contradiction with our theoretical results. Indeed, with coefficient jumps of in a checkerboard distribution, we have 1\text{\cdot}{10}^{10}, and $\tau_{EF}=$1\text{\cdot}{10}^{10}, in agreement with the better performance of E with respect to F and VF, as shown in Tables 5, 6. The independence from the jumps of the coefficients is shown by VE as predicted by the bound of Theorem 4.6.
6.2 FETI-DP quasi-optimality
We test quasi optimality for both the constant coefficient case and the varying coefficients case, again on the lowest order method (). Results for our second set of runs (Test 2) with both smooth and random data are shown in Figures 3 for meshes of truncated octahedra, and in Figures 4 for Voronoi meshes. In both cases, the results are in agreement with the estimates of Theorem 4.6. In particular, the experiments show that FETI-DP achieves quasi-optimality for our model problem if VE is used, i.e. \widetilde{W}_{h}\subset{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(V)}}\cap{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{W}_{h}^{(E)}}.
We conclude by a first test of the FETI-DP preconditioner on order VEM, for (more extensive testing will be carried out in the future). More precisely, for (we include the lowest order case for comparison), we test the most efficient choice of primal space, namely, Algorithm VE, with, once again, two types of data: constant coefficients and the 3D checkerboard coefficient distribution . For this last test, the forcing term and Dirichlet boundary condition on are chosen so that the exact solution is . The basis is chosen as the rescaled monomial basis and for the stabilization we use the diagonal recipe (2.15).
For this test the stopping criterion is the relative reduction of the dual residual by . In a C implementation, the PETSc built-in LU solver is used as subdomain solver, whereas we use GMRES with absolute convergence tolerance of as the coarse solver. These last experiments are run on the EOS cluster of the University of Pavia (http://matematica.unipv.it/it/cluster-eos). We used nodes equipped with processors Intel(R) Xeon(R) Gold 6130 CPU @ 2.10 GHz, operating system GNU/Linux 3.10.0-957.1.3.el7.x86_64.
Table 7 lists the values of the geometrical parameters for the reference meshes used in the experiments, where Npol denotes the number of polyhedra. Table 8 shows, for the case of constant, the number of iterations and the estimated condition number for the FETI-DP preconditioned system with a fixed number of subdomains () but varying both the polynomial degree and the reference mesh. The numerical tests with a 3D checkerboard distribution of are displayed in Table 9. Focusing, for fixed, on the number of iterations rather that on the condition numbers (which are only estimated), we observe that the results are in agreement with the theory. As is to be expected, the performance deteriorates as increases. Indeed, several of the inequalities on which the theoretical estimate relies on are (possibly heavily) depending on . In particular it is well known that, as increases, the choice of the stabilization bilinear form and of the face and interior degrees of freedom become crucial for the performance of the method [47, 34]. It is out of the scope of this paper to study the dependence of the performance of the preconditioner on such choices, and on the order of the method, a question that we plan to address in the future.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] B. Ahmad, A. Alsaedi, F. Brezzi, L. D. Marini, and A. Russo , Equivalent projectors for virtual element methods , Comput. Math. Appl., 66 (2013), pp. 376–391.
- 2[2] P. F. Antonietti, A. Cangiani, J. Collis, Z. Dong, E. H. Georgoulis, S. Giani, and P. Houston , Review of discontinuous Galerkin finite element methods for partial differential equations on complicated domains , in Connections and Challenges in Modern Approaches to Numerical Partial Differential Equations, Springer, 2016, pp. 279–308.
- 3[3] P. F. Antonietti, L. B. da Veiga, D. Mora, and M. Verani , A stream virtual element formulation of the Stokes problem on polygonal meshes , SIAM J. Numer. Anal., 52 (2014), pp. 386–404.
- 4[4] P. F. Antonietti, L. B. da Veiga, S. Scacchi, and M. Verani , A C 1 superscript 𝐶 1 C^{1} virtual element method for the Cahn-Hilliard equation with polygonal meshes , SIAM J. Numer. Anal., 54 (2016), pp. 34–56.
- 5[5] P. F. Antonietti, L. Mascotto, and M. Verani , A multigrid algorithm for the p 𝑝 p -version of the virtual element method , ESAIM Math. Model. Numer. Anal., 52 (2018), pp. 337–364.
- 6[6] L. Beirão da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. D. Marini, and A. Russo , Basic principles of virtual element methods , Math. Models Methods Appl. Sci., 23 (2013), pp. 199–214.
- 7[7] L. Beirão da Veiga, F. Brezzi, and L. D. Marini , Virtual elements for linear elasticity problems , SIAM J. Numer. Anal., 51 (2013), pp. 794–812.
- 8[8] L. Beirão da Veiga, F. Brezzi, L. D. Marini, and A. Russo , The hitchhiker’s guide to the virtual element method , Math. Models Methods Appl. Sci., 24 (2014), pp. 1541–1573.