A Multigrid Algorithm for the Mortar Finite Element Method.
The objective of this paper is to develop and analyze a multigrid algorithm for the system of equations arising from the mortar finite element discretization of second order elliptic boundary value problems. In order to establish the inf-sup condition for the saddle point formulation and to motivate the subsequent treatment of the discretizations, we first revisit briefly the theoretical concept of the mortar finite element method. Employing suitable mesh-dependent norms we verify the validity of the Ladyzhenskaya--Babuska--Brezzi (LBB) condition for the resulting mixed method and prove an L2 error estimate. This is the key for establishing a suitable approximation property for our multigrid convergence proof via a duality argument. In fact, we are able to verify optimal multigrid efficiency based on a smoother which is applied to the whole coupled system of equations. We conclude with several numerical tests of the proposed scheme which confirm the theoretical results and show the efficiency and the robustness of the method even in situations not covered by the theory.
Introduction
The mortar method as a special domain decomposition methodology appears to be
particularly attractive because different types of discretizations can be employed in
different parts of the domain. It has been analyzed in a series of papers [5, 6, 20]
mainly in connection with second order elliptic boundary value problems of the form
@n
where a(x) is a (sufficiently smooth) uniformly positive definite matrix in the bounded
domain\Omega ae R d , \Gamma D is a subset of the boundary \Gamma
Suppose
that\Omega is decomposed into non-overlapping
(\Omega\Gamma denote the usual Sobolev spaces endowed with the Sobolev norms k \Delta k
s;\Omega ,
0;D be the closure in H 1 of all C 1 -functions vanishing on \Gamma D . Although
this is not our motivation, a common approach to facilitate parallel computations is
to seek for a variational formulation of (1.1) with respect to the product space
1(\Omega
endowed with the norm
The space H 1
0;D is determined as a subspace of X ffi by appropriate linear con-
straints. Corresponding discretizations lead to saddle point problems. The central
objective of this paper is to develop a multigrid method for the efficient solution of
such indefinite systems of equations. According to standard multigrid convergence
theory the main tasks are to establish appropriate approximation properties in terms
of direct estimates as well as to design suitable smoothing procedures which give rise
to corresponding inverse estimates.
The derivation of these multigrid ingredients is, of course, based on the stability
of the discretizations which in turn hinges on a proper formulation of the continuous
problem (1.1). In this regard the subspace of those functions for which the jumps
across the interfaces of neighboring subdomains belong to the trace space H 1=2
a more suitable framework than the full space X ffi . Since H 1=2
00 is endowed with a
strictly stronger norm than H 1=2 , the concealed extensions of trace functions to the
domain require special care in the study of mortar elements. In fact, a complete
verification of the inf-sup condition is usually circumvented.
In order to establish the facts required for the above objectives in a coherent
manner we briefly reformulate in Section 2 first an analytical basis of the mortar
method. A concept with an explicit relation to the trace spaces H 1=2
00 has been given
in the recent investigation [4] which also contains a proof of the inf-sup condition.
Our present considerations, in particular a verification of the ellipticity in the discrete
case leads us to adjust the norms and to use mesh-dependent norms as in [20]. Once
a decision on the norms has been done, the analysis proceeds almost independently
of the fact whether it is done in the framework of saddle point problems or of the
theory of nonconforming elements.
The results of Sections 2 and 3 will be used in Section 4 to estimate the convergence
of the solutions of the discrete problems in the L 2 norm. This will serve
as a crucial ingredient of the multigrid convergence analysis in Section 5. Finally,
in Section 6 we present some numerical experiments, demonstrating h-independent
convergence rates for domains with one and more cross-points, distorted grids, different
mesh sizes at the interfaces and strongly varying diffusion constants.
The extension to more robust smoothers, a comparison with other transforming
iterations (cf. [19, 2]) and the application to noncontinuous cross-points and other
mortar situations (cf. [18]) will be presented in a forthcoming paper.
We will always denote by c a generic constant which does not depend on any of
the parameters involved in respective estimates but may assume different values at
each occurrence.
2 The Continuous Problem - A Characterization
of H 1
For simplicity we will assume throughout the rest of the paper that the domain
\Omega ae R d as well as the
subdomains\Omega k in (1.2) are polygonal.
and\Omega l share a
common interface, we set -
-\Omega l . The interior faces form the skeleton
k;l
will always be assumed to be the union of polygonal subsets of the
boundaries of
the\Omega k . Often such a decomposition is called geometrically conforming.
0;D is characterized as a subspace of X ffi , the space
is encountered. Here H(div; \Omega\Gamma is the space of all vector fields in L 2
d whose (weak)
divergence is in L
denotes the outward normal and q \Delta p is the standard scalar
product on R d . We recall from Proposition III.1.1 in [12] that
and the domain decomposition concept discussed in [12] is based upon (2.2). Note
that this characterization involves constraints with q 2 H 0;N (div; \Omega\Gamma which are global
in the sense that their restrictions to any interface \Gamma kl does not necessarily define
a bounded functional on the corresponding restriction of traces. By contrast the
mortar method attempts to reexpress (2.2) in terms of the jumps [v] of
across the interfaces \Gamma kl as
where
This, however, suggests that the individual terms (q
are well-defined
which requires restricting the jumps [v]. While this has been done on a discrete
level in [5, 6], a rigorous analysis has to be based on a corresponding continuous
formulation which requires suitable Sobolev spaces for the functions that live on the
skeleton S.
To this end, recall that for any (sufficiently regular) manifold \Gamma the Sobolev
spaces H s (\Gamma) can be defined by their intrinsic norms (see e.g. [14], Section 1.1.3 or
[16], Section 7.3 and [7] for a short summary tailored to the needs of finite element
discretizations) or alternatively, when \Gamma is part of a boundary, as a trace space. In
is not an integer,
kvk
@\Omega =v
kwk
is an equivalent norm for H s\Gamma1=2
(@\Omega\Gamma2 see e.g. [14], Theorem 1.5.1.2. Here and in
the sequel we will not distinguish between equivalent norms for the same space.
is a smooth subset of \Gamma, H s
consists of those elements
trivial extension ~
v by zero to all of \Gamma belongs to H s (\Gamma), cf. [16], p. 66. In the present
context the spaces H 1=2
are relevant, i.e.,
strictly contained continuously embedded subspace of H 1=2 (\Gamma 0 ), see e.g.
[14], Corollary 1.4.4.5 or [16], Theorem 11.4. By definition,
kzk
For later use we will also record the characterization of H 1=2
as an interpolation
space between L 2 (\Gamma kl ) and H 1
while
This can be realized, for instance, by the K-method [16], pp. 64-66, pp. 98-99.
Now we return to the characterization of H 1
0;D(\Omega\Gamma3 Whereas the constraints in
are non local, restricting the jumps [v] of elements in X ffi in a way that their
restriction to \Gamma kl belongs to H 1=2
indeed they become local. For this concept
the constraints are to be described through the product space
where (H 1=2
denotes the dual of H 1=2
suggests considering the
space [4]
endowed with the norm
The norms that will be introduced later for the treatment of the finite element
discretization are better understood from this norm than from (1.4).
Remark 2.1 X 00 is not a Hilbert space with respect to k \Delta k 1;ffi because H 1=2
not a closed subspace of H 1=2 (\Gamma kl ) with respect to the norm k \Delta k 1=2;\Gamma kl
Proof: By the remarks in the previous section, there exists a sequence fw n g
of C 1 functions with compact support in \Gamma kl such that kw n k 1=2;\Gamma kl
1. Since the w n admit uniformly bounded extensions
to all of
@\Omega k one can, in view of (2.5), construct a sequence of functions v n in
X 00 for which kv n k 1;ffi is bounded while kv n kX tends to infinity. Thus, the injection
Assuming that X 00 is
closed in X ffi when equipped with the weaker norm, the closed graph theorem would
lead to a contradiction. In fact, the boundedness of ' would imply the closedness of
the graph and hence the boundedness of ' \Gamma1 .
Remark 2.2 X 00 endowed with the norm k \Delta kX is a Hilbert space, and one has the
continuous embeddings
0;D is the closure with respect to k \Delta k
1;\Omega of all patch-wise smooth
globally continuous functions
on\Omega which vanish on D, the first inclusion is clear.
By the definition of X 00 we have for z 2
If fv n g is a Cauchy sequence in X 00 , it is also a Cauchy sequence in X ffi and converges
to some v 2 X ffi . From (2.11) it follows that f[v n is a Cauchy sequence in
converges to some g 2 H 1=2
00 (\Gamma kl ). By standard arguments we conclude
that
As a direct consequence we have
We now turn the problem (1.1) into a weak form based on the above characterization
of H 1
Z\Omega
Setting
Y
we consider the variational problem: find (u; -) 2 X 00 \Theta M such that
From (2.11) it follows that the operator
defined
by (Bv; -)
for any - 2 M , is bounded.
Moreover, the saddle point problem (2.14) satisfies the inf-sup condition. Indeed,
since the dual of M is M
sup
Now pick w kl 2 H 1=2
its extension to S by zero as w 0
kl . By definition of H 1=2
we have w 0
Hence assigning each \Gamma kl ae S to exactly one
(whose boundary contains \Gamma kl ), the local problems
kl
have unique solutions in H 1
defined by
z kl is easily
seen to belong to X 00 while, by construction
Since kz kl k
(see e.g. [12], p. 90) this confirms
the claim on the inf-sup condition.
Furthermore, we know from (2.12) that H 1
kvk
0;D , the bilinear form a(\Delta; \Delta) is V -elliptic.
3 The Discrete Problem
We will now turn to a conforming finite element discretization of (2.14). Throughout
the remainder of this paper we will restrict the discussion to the bivariate case
2. For each
subdomain\Omega k we choose a family of (conforming) triangulations
independently of the neighboring subdomains, i.e., the nodes in T k;h that belong
to \Gamma kl need not match with the nodes of T l;h . The corresponding spaces of piecewise
linear finite elements on T k;h are denoted by S h (T k;h ). We set
Y
i.e., the functions in X h are continuous at the cross-points of the polygonal subdo-
mains\Omega k . We associate with each interface \Gamma kl the non-mortar-side which by the
usual convention
while\Omega l is the mortar-side. Let T kl;h be the space of all continuous
piecewise linear functions on \Gamma kl on the partition induced by the triangulation
T k;h on the non-mortar-side, under the additional constraint that the elements in
T kl;h are constant on the two intervals containing the end points of \Gamma kl . Thus the dimension
of T kl;h agrees with the dimension of ~
The space of discrete multipliers is defined as
Y
Furthermore we denote the kernel of the restriction operator as
For simplicity of our notation, we assume that each \Gamma kl corresponds to one edge
of the polygonal
domain\Omega k . By including additional cross-points \Gamma kl can be divided
into parts fl i , where the mortar side of each fl i can be k or l, cf. [6]. Since this
extension does not change the analysis we will not further burden our notation with
such distinctions.
For convenience, we have labelled the finite element spaces by a global mesh
size parameter h. On the other hand, since the mortar method aims at combining
different possibly independent discretizations, we will admit different individual
discretization parameters on the subdomains. Whenever this is to be stressed, we
will denote by h k the mesh size
on\Omega k . This is actually related to the choice of the
mortar-side. In principle, one may choose the Lagrange multipliers on each part of
the skeleton from either adjacent subdomain. It will turn out, however, that the
side with the larger mesh size is the right choice if the corresponding adjacent mesh
sizes differ very much. Nevertheless, in order to avoid a severe cluttering of indices
and to keep the essence of the reasoning as transparent as possible, we will generally
suppress an explicit distinction of local mesh sizes. Instead we will always tacitly
assume the following convention.
Hypothesis (M). If each \Gamma kl ae S is labelled such
l is the mortar side, then
this choice of the mortar side has been made such that
ch k
holds with c being a constant of moderate size.
Thus the general convention will be that whenever h appears in a summand
related to \Gamma kl it is to be understood as h k , the larger of the two adjacent mesh sizes,
while otherwise h stands for the globally maximal mesh size. It will be seen that
with this choice the stability of the discretization holds uniformly for arbitrarily
varying mesh sizes.
The main handicap of all attempts to analyze the mortar element method is the
fact that we have
1;\Omega l
recall Remark 2.1, c.f. also Lemma 3.5 below. Therefore, following [20], we will use
mesh-dependent norms. Setting
let
\Gamma1=2;h :=
Whenever a distinction of local mesh sizes matters, the global h in (3.5)-(3.6) has
to be replaced by h k in the summands for \Gamma kl , i.e. in agreement with Hypothesis
(M) by the larger value of the mesh size of the neighboring subdomains.
Obviously, we have in analogy to (2.11), also by definition
In this framework we will prove that
is a stable discretization of (2.14).
In order to simplify the treatment of the trace spaces, let
be a partition of the interval [- which represents an interface \Gamma kl . We always assume
that such partitions are quasi-uniform since inverse estimates will be frequently
used. Motivated by the setting (3.2) of ~
T kl;h and T kl;h we consider two subspaces of
the space of continuous piecewise linear functions on [- h be the subspace
of those functions that vanish at the endpoints - 0 and - p , and let T h be the subspace
of those functions that are constant on the first and on the last interval. So S h and
have the same dimension p \Gamma 1.
For convenience, we suppress the explicit reference to the interval [-
there is no risk of confusion. In particular, the standard inner product on
will be denoted by (\Delta; \Delta) 0 , and the associated L 2 -norm by k \Delta k 0 .
Lemma 3.1 The projectors defined by
are uniformly bounded in L 2 , specifically
kfk
Proof: Since we are dealing with the 2-dimensional case, the proof is easy. For
be defined by v h (- i
and v h agree on [-
Z
Z
On the other hand, one obtains for the first (and last) interval
Z
Z
Z
where D := (-
Z
Z
Summing over all intervals and using Young's inequality yields
which proves (3.11).
Since the subsequent discussion involves also the special properties of the spaces
00 we briefly recall the following interpolation argument which, in principle, is
standard, see e.g. Proposition 2.5 in [16]. Suppose that Y; X are Hilbert spaces that
are continuously embedded in the Hilbert spaces Y; X , respectively, and that L is
a bounded linear operator from X to Y and from X to Y with norms C and C,
respectively. Then L is also a bounded linear operator from the interpolation space
bounded by
As a first application we state the following inverse estimate
which is perhaps only worth mentioning since a somewhat stronger norm than the
usual 1=2-norm appears on the left hand side. In fact, the standard inverse estimate
in H 1 ensures that kv h
and, by [14],
(1.4.2.1), p. 24, one has kv h k H 1= kv h k 1 one can take combined with (3.13) confirms (3.14).
We are now in a position to formulate the main result of this section.
Theorem 3.2 Assume that the triangulation in each
subdomain\Omega k is uniform and
that Hypothesis (M) is satisfied. The discretizations (3.8) based on the spaces X
defined by (3.1) and (3.2), respectively, satisfy the LBB-condition, i.e., there exists
some fi ? 0 such that
sup
holds uniformly in h.
Proof: The first part of the proof follows the lines of [20], where the LBB-condition
is proved for X h with a different norm. Given - 2 M h and \Gamma kl ae S, we define d kl
on the boundary of the subdomain on the non-mortar-side by
d kl :=
Here Q kl refers to the projector from the previous lemma when applied to L 2 (\Gamma kl ).
Since the triangulation
on\Omega k is assumed to be uniform, by Lemma 5.1 in [7], there
is an extension G k;h (d kl
Moreover, kd kl k
. Now we set
G k;h (d kl ) for x 2
From (3.12) in the proof of Lemma 3.1 we know that
Z
Z
-d kl ds - 3k-k 0;\Gamma kl kd kl k
On the other hand, the inverse inequality (3.14) yields
ch \Gamma1=2 kd kl k 0;\Gamma kl ; (3.20)
and it follows from (3.18) and (3.17) that kv kl k
ch \Gamma1=2 kd kl k 0;\Gamma kl .
we conclude, on account of (3.19) and (3.6), on one
hand, that
R
-[v kl ]ds - c k-k \Gamma1=2;h;\Gamma kl
and, on the other hand, also that
R
-[v kl ]ds - c k-k \Gamma1=2;h;\Gamma kl
which in summary yields
Z
-[v kl ]ds - ck-k \Gamma1=2;h;\Gamma kl
Now the assertion is obtained from (3.21) by summing over all \Gamma kl ae S.
Note that v kl in the preceding proof is constructed with respect to the non-
mortar-side. Therefore the larger mesh size enters into (3.20), and varying mesh
sizes are no problem in the proof.
We need two more approximation properties.
Remark 3.3 Let I h denote the Lagrange interpolation operator onto X h . Then we
have
ch 3=2 kvk
Proof: Let T be an element
in\Omega m and fl be one of its edges. By the trace theorem,
the mapping H continuous. By using the Bramble-
Hilbert lemma and the standard scaling argument we obtain
ch 3 jvj 2
By summing over all edges fl which lie on \Gamma kl , we obtain the assertion.
Remark 3.4 Given - 2 H 1=2 , there is a - h 2 T h such that
ch 1=2 k-k
Proof: To verify (3.23) consider the Lagrange interpolant L h at the interior points
to the intervals [- constants.
Since constants are reproduced, by the same arguments as those used in the proof
of the preceding remark one confirms that k- \Gamma L h -k 0 - chj-j 1 . Moreover, the L 2
projector
while also kP h . Thus, since H
applying (3.13) with
It remains to prove the ellipticity of a on the kernel V h of the operator B. In our
first approach we had a stabilizing term in the bilinear form to compensate (3.4).
This is not necessary. Recently, C. Bernardi informed us about some techniques
used now in [1]. When adapting it for obtaining a good alternative of (3.4), we
observed the analogy to Lemma 2.2 in [20]. For the reader's convenience we present
it with a short proof.
Lemma 3.5 Assume that the triangulations in each
subdomain\Omega k are shape regu-
lar. Then
1;\Omega l
kl;h be the L 2 -projector. If v h 2 V h , then by the
matching condition
R
the
jumps of H 1 -functions vanish in the sense of (2.12), we have P kl
From this, Remark 3.4, and the trace theorem we conclude that
)vj\Omega l
ch 1=2
kvj\Omega l
ch 1=2 (kvk
1;\Omega l
Finally, we divide by h 1=2 , and the proof is complete.
A direct consequence is obtained by the Cauchy-Schwarz inequality:
We note that we need the ellipticity of a h merely on V h and not on the larger set
for details cf. Lemma III.1.2 in [8]. First, the inequality
has been often used in the analysis of mortar elements. It has been proven in [6]
by a compactness argument as often found in proofs of non-standard inequalities of
Poincar'e-Friedrichs type. Moreover, from (3.25) and (3.26) it follows that
for
so that, in view of the definition (3.5), the proof of the ellipticity is complete.
Note again that the above reasoning remains valid under Hypothesis (M) concerning
the choice of the mortar side if the local mesh sizes vary significantly.
To keep the development as transparent as possible we have dispensed with
striving for utmost generality. An extension of the majority of the results to finite
element spaces with polynomials of higher degree is rather straightforward. Only a
comment on the corresponding version of Lemma 3.1 is worth mentioning. The rest
of this section is devoted to the adaptation of the lemma and may be skipped by
the reader.
Remark 3.6 Suppose that the family S h (T k;h ) in (3.8) consists of Lagrange finite
elements of degree n. One has to consider analogous spaces S h and T h consisting
of globally continuous piecewise polynomials of degree at most n with respect to the
partition (3.9). While the elements of S h vanish at the endpoints, it is consistent
with the usual setting, cf. [4], to require that the elements of T h have degree at most
on the first and last interval in order to ensure matching dimensions. The
essence of the proof for an analog to Lemma 3.1 is to show that
sup
(v
for some constant C independent of h.
So given any u h 2 S h we have to construct v h 2 T h such that (3.27) holds. As in
the proof of Lemma 3.1 we can choose v h to coincide with u h on the subset [-
so that it remains to analyse the first and last interval. By symmetry it suffices
to consider the first interval which for convenience is taken to be [0; 1]. Setting
h dx, or equivalently thatZv h
where z n is the largest zero of the n-th degree Legendre polynomial. Since v h has
degree the Gaussian quadrature formulae on [0; 1] are exact for the polynomials
Here the x k := (1 are the zeros of the n-th degree Legendre polynomial with
respect to the interval [0; 1] and the ae k 's are the weights. Since
and the latter estimate is sharp
because v h can be chosen to vanish at x Together with the inequality x 2 - x
we
This proves (3.27) with C := The remainder of the proof of Lemma
3.1 is now analogous when 3=8 is replaced by C.
A crucial ingredient of the convergence analysis of the multigrid algorithm is an
error estimate of the finite element solution. Although the existence of such an
estimate was mentioned in [4], no proof was provided there. Therefore, we will give
a proof here.
Up to now we have considered the mortar element method as a mixed method,
but one may interpret it as a non-conforming method with finite elements in
We will establish the error estimates in the framework of non-conforming elements.
On the other hand, we will do this by making use of the results for the mixed method
from the preceding section. We note that by Lemma 3.5 the norms k \Delta k 1;h and k \Delta k 1;ffi
are equivalent on V h .
In this context it is natural to assume H 2 -regularity, i.e.
Theorem 4.1 Suppose that the triangulations T h;k are shape regular and that the
variational problem (1.1) is H 2 -regular. Then the finite element solution u
satisfies
where h in the maximum of the mesh sizes of the triangulations.
Proof: By Strang's second lemma, see e.g. [8], p. 102 we have
R
a @u
@n [v h ]ds
Indeed, following [8], p. 104, integration by parts provides the following representation
of the consistency error in (4.2)
0;@\Omega
Z
Z
Z
@\Omega g(x)v h (x)ds
Z
Z
@n
Z
Z
@\Omega g(x)v h (x)ds
(a @u
@n
@n
Note that it is clear from (4.3) that the Lagrange multiplier - in (2.14) coincides
with a @u
@n . Since v h 2 V h , orthogonality (3.3) allows us to subtract an arbitrary
element h from the first factor so that
(a @u
@n
ka
@n
From the trace theorem we know that ka @u
@n
. Now
Remark 3.4 and Lemma 3.5 followed by the application of the Cauchy-Schwarz
inequality to the sum provide the estimate
jL u (v h )j - ch kuk
and the quotient in (4.2) is bounded by ch kuk
2;\Omega .
In order to establish a bound of the approximation error in (4.2), we start with
the approximation of u on each
2(\Omega\Gamma0 we can take
the interpolant in S h (T k;h ). In particular the cross-points may be chosen as nodal
points and we obtain interpolants which are continuous at the endpoints of each
. The well-known results for conforming P 1 -elements yield an estimate for the
approximation in the (larger) space
ch kuk
By Remark 3.3 and the triangle inequality, we obtain for \Gamma kl ae S:
ch 3=2 (kuk
2;\Omega l
Hence,
2;\Omega with I h u 2 X h . We do not say that I h u is contained
in V h . On the other hand, from the inf-sup condition in Theorem 3.2 above and
Remark III.4.6 in [8] we conclude that the estimate in X h yields an upper bound for
the approximation in the kernel V h ,
In particular, the general arguments above cover also the case with mesh-dependent
norms. Combining (4.5) and (4.6) yields
ch kuk
and the bound of the second term in the desired estimate (4.1) has been established.
To obtain the L 2 error estimate we move completely to the theory of nonconforming
elements. Here (4.5) and (4.7) provide the ingredients for the duality argument
in Lemma III.1.4 from [8] which yields (4.1) and completes the proof. The details
are obvious from the treatment of the Crouzeix-Raviart element in [8], p.106.
For completeness we mention that the error estimate in the energy norm can be
improved if the individual mesh sizes of the subdomains are incorporated. On the
other hand, the duality argument yields only a factor of if based on
the regularity estimate kuk
For an error estimate of the Lagrange multiplier the reader is referred to [20].
5 Multigrid Convergence Analysis
The saddle point problem (3.8) gives rise to a linear system of the form
!/
where the dimension of the vectors coincides with the dimension of the finite element
spaces X h and M h , respectively. For convenience, the same symbol is taken for
the finite element functions and their vector representations, and the index h is
suppressed whenever no confusion is possible.
We will always assume that the finite element basis functions are normalized
such that the Euclidean norm of the vectors k \Delta k ' 2
is equivalent to the L 2 -norm of
the functions, i.e.
When the equations (5.1) are to be solved by a multigrid algorithm, the design of
the smoothing procedure is the crucial point. Motivated by [10] (see also [9]) our
smoothing procedure will be based on the following concept. Suppose that C is a
preconditioner for A which, in particular, is normalized so that
and for which the linear system
e
is more easily solvable. In actual computations the vectors v; - are obtained by
implementing
where
is the Schur complement of (5.4).
In particular, we can divide the vector v into two blocks with the first block
associated to the nodal values in the interior of the subdomains and the second
block associated to the values on the skeleton. A non-diagonal preconditioner is
chosen only for the first blockB @
-C A =B @
Here, the system splits into two parts, and the Schur complement matrix
ff
Thas a simple structure. The dimension corresponds to the number of nodal points
from the non-mortar sides on the skeleton. Moreover, it consists of bands which are
coupled only through the points next to the cross-points. As a consequence, all the
other points can be eliminated with a very small fill-in.
The block matrix C 1 may be chosen as the ILU-decomposition of the corresponding
block of the given matrix A.
In order to facilitate the analysis we assume that also C 1 is a multiple of the
identity. Moreover we return to the original structure as in (5.1). Then the iteration
that will serve as a smoother in our multigrid scheme has the form
!/
where superscripts will always denote iteration indices. It is important to note that
h always satisfies the constraint, i.e.,
see [10]. Moreover (5.9) shows that the next iterate is independent of the old Lagrange
multiplier
h .
The coarse grid correction of the multigrid scheme can be performed in the
standard manner since the finite element spaces X h ae X 00 and M h ae M are nested,
see e.g. [10] or [15], p. 235. The smoothing property (5.12) below even shows that
we can abandon the transfer step proposed in [10].
As usually, the analysis of the multigrid method will be based on two different
norms. The fine topology will be measured by the norm
and the coarse one by
We recall that -
Smoothing property: Assume that - smoothing
steps of the relaxation (5.9) with C := ffI are performed, then
ch \Gamma2
Approximation property: For the coarse grid correction u 2h one has
ch
Proof of the smoothing property: We note that (5.12) is stronger than the
original version of the property in [10]. There we could estimate the jjj \Delta jjj 2 -norm
only if the Lagrange multiplier - m
h was replaced by a suitable one which could be
determined by solving once more an equation with the matrix (5.7). R. Stevenson
[17] observed that the extra step yields the same Lagrange multiplier as another
smoothing step. So the following proof is based on Stevenson's idea. As usual, only
purely algebraic properties are used.
From (5.9) we have
We may extract a recursion for the u component
with the projector
The matrix M := P
ff A)P is symmetric. Since u 1 satisfies the constraint
Since the recursion (5.14) is linear, we may assume without loss of generality that
which implies multiply (5.14) by the matrix in order to
eliminate the inverse matrix and consider the first component. Then we obtain with
or
ff
From (5.16) it follows that for m - 3
ff
ff
Note the symmetry of the matrix on the right-hand side of (5.18). Since ff -
(A), the matrix I \Gamma 1
ff A is a contraction for the ' 2 -norm. Moreover the spectrum
of M is contained in [0,1], and we conclude by the usual spectral decomposition
argument
We have assumed ff - c h \Gamma2 and obtain now from (5.18)
ch \Gamma2
The in the denominator may be replaced by m if we adapt the factor
c. This proves (5.12) for m - 3. The cases can be treated by
recalling (5.17) and
Proof of the approximation property: The proof of the approximation property
depends on the individual elliptic problem. Here it is convenient to consider
the mortar elements as nonconforming elements, and we follow the treatment of
the multigrid algorithm for the nonconforming P 1 -element in [11]. For a current
approximation the residual
is given by
We recall from (5.2) that the Riesz-Fischer representation r 2 X h ae L
2(\Omega\Gamma of d
defined by
satisfies
Let now z be the solution of the auxiliary variational problem
Since the variational problem on the
domain\Omega is assumed to be H 2 -regular, we have
are the finite element approximations for
the corresponding discretizations
respectively. The latter is true since X 2h ae X h and thus
(5.21) holds for . At this point the L 2 -estimate established in Theorem 4.1
comes into play. It is applied to the auxiliary problem (5.23) providing
0;\Omega
ch 2 kzk
ch 2 krk
ch 2 kdk ' 2
ch
which completes the proof of the approximation property.
After having established a smoothing property and the associated approximation
property it is clear from the standard multigrid theory that the W-cycles yield an
h-independent convergence rate if sufficiently many smoothing steps are chosen, see
e.g. [15], Chapter 7 or [8], pp. 222-228.
In most practical realizations a better performance is observed for multigrid
iterations with V-cycles and only a few smoothing steps. It is a common experience
that multigrid methods converge in more general situations than those assumed in
theoretical proofs. Furthermore, the bounds which are obtained for the convergence
rate strongly depend on the constants appearing in (5.12) and (5.13). Therefore a
realistic quantitative appraisal of the proposed scheme has to be based in addition
on complementary numerical experiments. Below we will report on our numerical
investigations for various typical situations of practical interest including also cases
that are not covered by our theoretical considerations.
6 Numerical Experiments
We present numerical examples implemented in the software toolbox UG [3] and its
finite element library. The efficiency of the method will first be studied in detail
for a model situation which is consistent with the assumptions of the convergence
proof. In addition, we consider the robustness of the multigrid solver for mortar
finite element problems which are not covered by our analysis.
In all examples we use a multigrid method with a V-cycle and the simple
smoother (5.9), where is the Jacobi smoother for the local problems in
the
subdomains\Omega k . We want to point out that all presented convergence rates are
asymptotic rates.
Although the Schur complement of the smoothing matrix (5.7) is small and
could be assembled without loosing the optimal complexity of the algorithm, in our
implementation the equation (5.5) is solved iteratively. Specifically, a few cg-steps,
preconditioned with a symmetric Gau-Seidel relaxation for S in (5.7) are performed.
The implementation of the Gau-Seidel procedure makes use of the graph structure
of the matrices A and B. Thus, it has complexity O(dimM h ).
Since the bilinear form b contains no differential operator, the condition number
of S is bounded independently of the mesh size h, and a bounded number of steps
for the inner iterations independent of the refinement level is sufficient. For the
numerical experiments, we prescribe an error reduction factor 0:1 - ae - 0:5 for the
approximate solution of (5.5).
In the first test series, we investigate the multigrid convergence for the model
problem
\Gammadiv a grad
in a quadrilateral
domain\Omega := (0; 1) 2 which is split into four
;\Omega 11 .
The diffusion coefficients are assumed to be constant in every subdomain and are
denoted as a 00 Figure 1).
We start with a test for the Laplace operator, i.e. a
on a regular grid (cf. Table 1). We compare the asymptotic convergence rates
of the linear multigrid method with V(1,1)-cycle for the mortar problem and the
regular grid distorted grid different step sizes

Figure

1: Grids with one cross-point
simplified problem Au additional constraints. The latter
corresponds to the decoupled boundary value problems where in every subdomain
the equation is solved with Neumann boundary conditions for the inner boundary
without requiring continuity. For the mortar problem, we use the smoother (5.9)
with an inner reduction factor ae = 0:5 and ae = 0:1, for the simplified problem we
use the Jacobi smoother diag(A). Our results on a regular grid show that the same
rates are obtained with the mortar coupling, even for only one Schur complement
iteration. The fifth column in the table refers to the configuration with a
shows that the convergence is independent of the variation
of the diffusion constants in the different subdomains for a regular grid.
level elements different diff. const. without mortar elem.
nb. of inner iter. 1 1-2 1-2

Table

1: Regular grid: asymptotic convergence rates for a linear multigrid iteration
with a V(1,1)-cycle and resulting number of inner iterations
Note that due to (5.10) the multigrid iterates are contained in V h and that the
problem (3.8) is elliptic in V h . Thus, the multigrid iteration can be accelerated by
embedding it into a cg-iteration, i.e., a cg-iteration preconditioned by a multigrid
V-cycle is performed. There is another advantage. It is often difficult to find the
optimal damping factor ff. In particular, when the discretization with the slightly
distorted grid in Figure 1 is used, the Jacobi smoother diag(A) is not convergent
without appropriate damping. Thus, the cg-method is used here for computing
the correct damping factors automatically. However, since Equation (5.5) is solved
only approximately in actual computations, the multigrid iterates are not exactly
contained in V h . In other contexts the cg-method may be very sensitive with respect
to this point, but in our tests it has turned out that 3 inner iteration steps are
sufficient in all cases (cf. the entries in the last row of Tables 1-3).
The case in which the step sizes depend strongly on the
subdomain\Omega k was next
investigated. We observe better convergence if the Lagrange parameter on each
\Gamma kl is associated to the side of \Gamma kl with the coarser mesh. This is consistent with
Hypothesis (M). Although convergence is observed for a V(1,1)-cycle in all cases,
the asymptotic rate deteriorates for more than 100000 elements in the case of the
distorted grids and large jumps of the mesh sizes. On the other hand, the V(2,2)-
cycle turns out to be a robust preconditioner for the cg-iteration also in extreme
cases.
level elements regular grid distorted grid elements various step sizes
nb. of inner iter. 1-3 1-3 1-2

Table

2: Irregular grids: asymptotic convergence rates for cg with V(2,2)-cycle
In the next example we consider a typical mortar situation with several cross-
points. In Figure 2 large bricks are separated by thin channels. Fixing the diffusion
constant for the bricks to a we test the cases where the channels have higher
or lower permeability (a We perform the cg-method
with V(1,1)-cycle and two inner iterations. We obtain stable convergence rates if
the mortar side is on the side with the smaller diffusion constant and large step size,
resp.; otherwise the method may fail. The results in Figure 2 for the case a
show clearly that the diffusion is faster in the small channels.
level elements a
number of inner iterations 1-2 1-2 1-2

Table

3: Convergence for the example with several cross-points for cg with V(1,1)-
cycle
Finally, we apply the multigrid method to an example for a rotating geometry
with two circles which occurs for time dependent problems (cf. [18]). We use a cg-
iteration with a damped Jacobi smoothing diag A. Subdomains with curved

Figure

2: Example with several cross-points
boundaries are not covered by our theory since the approximation of the curved
boundaries induces an additional consistency error. Note that the exact Lagrange
parameter is piecewise constant and discontinuous for a linear solution. Thus, linear
functions cannot be represented by the mortar ansatz space. This results in worse
convergence rates. Nevertheless, the method is stable when a cg-iteration is applied,
preconditioned by a V(3,3)-cycle and 3 inner iterations for the Schur complement
equation (cf. Table 4). On the other hand, without cg-acceleration a V(4,4)-cycle
and a strong damping for the Jacobi smoother is required.
level elements convergence rate
9 2097152 0.37

Table

4: Rotating geometry (parallel computation on 128 processors)
In summary, we have demonstrated the robustness of the method with respect to
the number of subdomains, different step sizes in the subdomains, and varying diffusion
constants. Thus, this is a very efficient solver for mortar finite elements. The
convergence rates are independent of the mesh size and the number of refinement
levels. The results show clearly that the presented smoother with inexact solution of
the corresponding Schur complement is very efficient and that no further improvement
is expected from a more accurate solution of the Schur complement equation.
In practice, for a nested multigrid cycle the accuracy of the approximation error is
obtained within one or two V(1,1)-cycles.
Of course, our tests concern the robustness with respect to different mortar
situations, whereas the equations on the subdomains are simple. For more involved
problems the smoother C has to be replaced, e.g. by a more robust ILU smoother.
Nevertheless, smoothers which are decoupled from the mortar interfaces as in (5.8)
are recommended in order to retain the low complexity.
Our numerical experiments confirm that the quality of the solver depends strongly
on the right choice of the mortar side. In extreme cases the method diverges if the
Lagrange parameter is associated with the wrong side. Apparently a good rule of
thumb is to choose the Lagrange parameter for that side for which the quotient a=h 2
attains the smaller value.



--R


A class of iterative methods for solving saddle point problems

The mortar finite element method with Lagrange multipliers
The mortar element method for three dimensional finite elements.
"Nonlinear Partial Differential Equations and Their Applications"
Iterative methods for the solution of elliptic problems in regions partitionned in substructures

A Cascade algorithm for the Stokes equation
An efficient smoother for the Stokes problem.

Mixed and Hybrid Finite Element Methods
Approximation by finite element functions using local regulariza
Elliptic Problems in Nonsmooth Domains



The coupling of mixed and conforming finite element discretizations
On the convergence of multigrid methods with transforming smoothers.
Hierarchical a posteriori error estimators for mortar finite element methods with Lagrange multipliers.
--TR

--CTR
Micol Pennacchio, The Mortar Finite Element Method for the Cardiac Bidomain Model of Extracellular Potential, Journal of Scientific Computing, v.20 n.2, p.191-210, April 2004
V. John , P. Knobloch , G. Matthies , L. Tobiska, Non-nested multi-level solvers for finite element discretisations of mixed problems, Computing, v.68 n.4, p.313-341, September 2002
Analysis of mortar-type Q1rot/Q0 element and multigrid methods for the incompressible Stokes problem, Applied Numerical Mathematics, v.57 n.5-7, p.562-576, May, 2007
O. Steinbach, A natural domain decomposition method with non-matching grids, Applied Numerical Mathematics, v.54 n.3-4, p.362-377, August 2005
D. Braess , P. Deuflhard , K. Lipnikov, A subspace cascadic multigrid method for mortar elements, Computing, v.69 n.3, p.205-225, Nov. 2002
saddle point problems;domain decomposition;mortar method;trace spaces