Multilevel Boundary Functionals for Least-Squares Mixed Finite Element Methods.
For least-squares mixed finite element methods for the first-order system formulation of second-order elliptic problems, a technique for the weak enforcement of boundary conditions is presented. This approach is based on least-squares boundary functionals, which are equivalent to the H-1/2 and H1/2  norms on the trace spaces of lowest-order Raviart--Thomas elements for the flux and standard continuous piecewise linear elements for the pressure, respectively. Continuity and coercivity of the resulting bilinear form is proved implying optimal order convergence of the resulting Galerkin approximation. The boundary least-squares functional is implemented using multilevel principles and the technique is tested numerically for a model problem.
Introduction
. In the context of least-squares finite element methods for first-order
systems, boundary conditions can be enforced as essential boundary conditions
in the finite element spaces. This yields optimal order convergence of the Galerkin
approximations under suitable assumptions on the regularity of the problem (see, for
example, [18, 9, 10]). However, this approach constructs an approximation which
is much more accurate on the boundary than in the interior of the domain. For
least-squares finite element methods, a natural way of treating boundary conditions
is to enforce them weakly by adding boundary functionals. The boundary functional
approach is also a natural and simple way to handle the common situation that the
boundary conditions cannot be satisfied exactly by trial functions in the finite element
spaces (see, e.g., [16, 5, 13, 1], see also [4] for a different approach to handle boundary
conditions within a least-squares method). Nonlinear boundary conditions in the
modelling of flow in porous media (see, e.g., Bear [2, Chap. 7]) may also be handled
effectively by boundary functionals.
In this paper, we derive appropriate boundary functionals for the mixed formulation
of second-order elliptic problems. These are equivalent to the H 1=2 (\Gamma D )-norm
for the Dirichlet conditions and to the H \Gamma1=2 (\Gamma N )-norm on the Neumann boundary.
From the view of trace and extension theorems of Sobolev norms this gives the proper
balance between the interior and boundary least-squares functionals. These boundary
functionals are in some sense the weakest possible without losing the optimal order of
the Galerkin approximation. Related but different boundary functional methods in
connection to least-squares finite elements were suggested before in the literature. For
pure Dirichlet problems, a weighted L 2 (\Gamma) norm was used in [13]. This approach using
weighted L 2 norms was generalized to elliptic boundary value problems of A-D-N
type in [1].
Our aim in this paper is to justify our multilevel boundary functional approach
theoretically and by computational experiments for the mixed formulation of the Poisson
equation using lowest-order Raviart-Thomas elements. Clearly, the true potential
of this methodology lies in its applicability to more complicated problems, e.g., in
flow computations in porous media. The use of the least-squares mixed finite ele-
Fachbereich 6 (Mathematik und Informatik), Universit?t-GH Essen, 45117 Essen, Germany
ment approach for nonlinear boundary value problems arising in variably saturated
subsurface flow is studied in [20]. The techniques presented in this paper can also be
extended to least-squares formulations like, e.g., Maxwell, Stokes and Navier-Stokes
equations.
In the following section, we present the least-squares formulation including the
boundary functionals and prove coercivity and continuity of the corresponding bilinear
form. Section 3 reviews some results on the approximation properties of finite
element spaces including the lowest-order Raviart-Thomas elements. In Section 4, we
are concerned with an equivalent and computable boundary functional based on multilevel
principles. Section 5 gives a study of the effect of the boundary least-squares
functionals on the accuracy of the finite element approximation by computations for
a two-dimensional model problem.
2. Least-Squares Formulation with Boundary Functionals. We consider
the first-order system formulation of Poisson's equation,
div
oe
(\Omega\Gamma in a bounded polygonal
domain\Omega ae IR 2 . The boundary
of\Omega is
divided into
are prescribed. Actually, for \Gamma ' @
may be defined as the space of traces from H
1(\Omega\Gamma and H \Gamma1=2 (\Gamma) as the corresponding
dual norm with respect to L 2 (\Gamma) (see, for example, [11, Section I.1]).
Despite the fact that we restrict ourselves to two-dimensional problems for the
purpose of exposition, all the techniques presented below can be extended to higher
dimensions. Note that these techniques can also be extended to more general diffusion
problems like those arising in porous media flow.
Clearly, each solution (u;
1(\Omega\Gamma of this boundary value problem
also minimizes the least-squares functional
0;\Omega
kn
for any ff ? 0. With the corresponding bilinear form
0;\Omega
this is equivalent to finding (u; p) 2 H(div ; \Omega\Gamma \Theta H
1(\Omega\Gamma such that
for all (v;
We have the trace inequalities
kpk
1;\Omega for all
(see [11, Theorem 1.5]) with a constant c T and
kn
\Gamma1=2;@\Omega - kuk
div;\Omega for all u
MULTILEVEL BOUNDARY FUNCTIONALS FOR LEAST-SQUARES METHODS 3
(see [11, Theorem 2.5]).
Moreover, if \Gamma D is a curve of positive measure, the generalized Poincar'e-Friedrichs
inequality in [15, Theorem 1.9] gives
kpk
for all
(\Omega\Gamma with a constant c F . If \Gamma restrict ourselves to H
normalizing p to, for example, (p; 1)
This construction is necessary in order
to ensure uniqueness of p and we have a Poincar'e-Friedrichs inequality of the form
kpk
0;\Omega for all
which satisfy (p; 1)
(cf. [3, Section II.3]).
These tools allow us to prove coercivity and continuity of the bilinear form for
any ff ? 0.
Theorem 2.1. Under the assumptions above, for any ff ? 0, the bilinear form
B(\Delta; \Delta; \Delta; \Delta) is coercive and continuous with respect to H(div; \Omega\Gamma \Theta H
and
with positive constants c S and c E (which depend on ff, c T and c F ).
Proof. The above trace inequalities and repeated use of the Cauchy-Schwarz
inequality lead to
0;\Omega kdiv vk
div;\Omega kvk
1;\Omega kqk
which proves (2.9).
For the coercivity proof we consider two separate cases: (i) \Gamma
@\Omega a curve of positive measure.
Case
Cauchy-Schwarz inequality and
ff
kn
ff
1=2;@\Omega
lead to
0;\Omega
0;\Omega
ff
1=2;@\Omega
ff
1;\Omega
for any ffi 2 (0; 1). Combined with (2.8), this implies
\Gamma1=2;@\Omega
F
1;\Omega
F
ff
0;\Omega
F
ff
Choosing
F )g.
Case positive measure. The first step is as in case (i) and using (2.6)
and (2.5) we obtain
kpk
kpk
kn
ff
ff
1;\Omega
which holds for any ffi 2 (0; ff). This leads to
0;\Omega
Combining this with (2.7) gives
0;\Omega
kuk
ff
F
Choosing
proves (2.10).
We remark that in [18], coercivity and continuity of the bilinear form
0;\Omega
is shown for
MULTILEVEL BOUNDARY FUNCTIONALS FOR LEAST-SQUARES METHODS 5
Since more emphasis is put on enforcing the boundary conditions as ff increases, this
result can be regarded as the limiting case for ff !1. Standard finite element theory
implies that the Galerkin approximation in these spaces is of optimal order. There are
many ways to achieve this optimal order convergence by boundary functionals which
are stronger than the k \Delta k \Gamma1=2;\Gamma N and k \Delta k 1=2;\Gamma D norms, respectively. The importance of
using the norms k \Delta k \Gamma1=2;\Gamma N and k \Delta k 1=2;\Gamma D in (2.2) is that the least-squares functionals
on the boundary and in the interior are properly balanced.
3. Galerkin Approximation. For the numerical approximation we consider
finite element spaces
based on a quasiuniform sequence of triangulations fT l g l=0;1;::: of \Omega\Gamma Let h l denote
the mesh-size of fT l g given, for example, by the maximal diameter of the triangles.
We compute approximate solutions u l 2 V l and p l 2 W l for u and p, respectively, such
that F(u l ; p l ; f) is minimized among all u l 2 V l and p l 2 W l . This is equivalent to
the variational problem of finding l \Theta W l such that
for all (v l ; q l l \Theta W l . The simplest choice is to use the lowest-order Raviart-
Thomas space for V l and standard continuous piecewise linear functions for W l on the
triangulation T l .
Continuity and coercivity proved in Theorem 2.1 give us the usual quasi-optimality
of the Galerkin approximation
l 2V l
q l 2W l
Let us assume that f 2 H ff
for some ff 2 (0; 1] and that g; h and the boundary
are such that p 2 H 1+ff
and, consequently, u 2 (H
details on the conditions for such regularity results). For example,
holds
if
domain\Omega with the additional property
that the interior angles at boundary points separating \Gamma N from \Gamma D are at most -=2.
In order to simplify our notation, we write - l . j l to indicate that - l - cj l holds
with a constant c which is independent of l. We also write - l h j l to indicate that
both - l . j l and j l . - l are satisfied. Thus, [8, Proposition III.3.9] implies
l 2V l
div;\Omega . h ff
l [kfk
for the approximation by lowest-order Raviart-Thomas elements and standard finite
element interpolation results (cf., e.g., [7, Chapter 4]) give
q l 2W l
1;\Omega . h ff
l
In order to compute the solution of (3.1), we use the bases f\Phi (-)
l g M l
l and
l g N l
-=1 for W l . The variational problem (3.1) may then be formulated as a linear
system of equations
- A uu A up
A pu A pp
- u
6 GERHARD STARKE
l )] -=1;:::;N l ;-=1;::: ;M l
l )] -=1;::: ;N l
l
l ) \Gamma1=2;\Gamma N
l
and l are the basis representations of u l 2
Setting up the matrices A uu and A pp involves the computation of
(n \Delta \Phi (-)
l )
l
Since we cannot compute these inner products directly, we need to replace k \Delta k \Gamma1=2;\Gamma N
by equivalent norms such that the corresponding inner products are
computable. Alternatively, one could replace these norms by stronger norms, e.g.,
norms do not give the optimal balance between accuracy on the boundary and in the
interior.
4. Multilevel Implementation of the Boundary Functionals. In this sec-
tion, we derive computable norms which are equivalent to the norms k \Delta k \Gamma1=2;\Gamma N and
using additive multilevel decompositions in the spirit of [6]. To this end,
we define the L 2 (\Gamma)-orthogonal projection Q
@\Omega onto (the trace space of)
l . We will abuse notation and denote
this trace space W l wherever it is clear from the context. With the operator
l
we have (cf. [17, Theorem 15] or [14, Corollary 3.2.4])
l
for all q l 2 W l .
For a computable expression that replaces k \Delta k \Gamma1=2;\Gamma N
we use the following result.
Theorem 4.1. Let V l and W l be the lowest-order Raviart-Thomas spaces and
standard piecewise linear continuous finite element spaces, respectively, based on a
quasiuniform sequence of triangulations. Then, for u l 2 V l ,
kn \Delta u l k \Gamma1=2;\Gamma N h sup
z l 2W l ;z l 6=0
Proof. By definition (see [11, Section I.1]),
kzk
MULTILEVEL BOUNDARY FUNCTIONALS FOR LEAST-SQUARES METHODS 7
Clearly, since W l ae H 1=2 (\Gamma N ), this implies
kn \Delta u l k \Gamma1=2;\Gamma N
z l 2W l ; z l 6=0
For the upper bound, we use the L 2 (\Gamma N )-orthogonal projection Q l;\Gamma N to (the trace
space of) W l . It is well-known (cf. [21, Section 4]) that
for all z together with the trivial inequalities
implies
l kzk 1=2;\Gamma N and kQ l;\Gamma N zk 1=2;\Gamma N . kzk
We may write
sup
kzk
kzk
kzk
Using the fact that n \Delta u l is piecewise constant along the boundary, we may estimate
the second term on the right hand side using
zk
. h l
kzk
kzk 1=2;\Gamma N . h 1=2
l kn \Delta u l k \Gamma1=2;\Gamma N
kzk
where the last step follows from the inverse inequality k- l k 0;\Gamma N . h \Gamma1=2
l k- l k \Gamma1=2;\Gamma N
for piecewise constant functions - l on the subdivision of \Gamma N corresponding to T l . This
implies
kn \Delta u l k \Gamma1=2;\Gamma N
kzk
l kn \Delta u l k
Choosing l 0 such that 2Ch 1=2
l 0
- 1, we are led to
kzk
. sup
z l 2W l ; z l 6=0
From (4.2) and (4.3) we are led to a computable expression by
kn \Delta u l k \Gamma1=2;\Gamma N h sup
z l 2W l ; z l 6=0
z l ) 1=2
z l 2W l ; z l 6=0
(n \Delta u l ; C \Gamma1=4
where C
is the adjoint of C l;\Gamma N such that
(n \Delta u l ; C \Gamma1=4
z l )
holds for all u l 2 V l ; z l 2 W l .
This gives rise to the modified bilinear form
0;\Omega
The variational problem of finding l \Theta W l such that
0;\Omega
for all (v l ; q l l \Theta W l replaces (3.1).
The computation of the contributions of the boundary functionals to the Galerkin
matrix is based on the mass matrices M j;\Gamma in W j ,
on the mass matrices ~
M j;\Gamma for coupling V j and W j ,
~
and on the restriction matrices I j
l which, for j - l, compute moments
-=1 from s l = [(s; \Psi (-)
l
(note that I l
l is just the identity matrix in IR N l ). It is easy to see that
holds. Of course, the inverse of M j;\Gamma can only be formed on the subspace of nodal
basis functions which do not vanish on \Gamma (we keep using this notation for simplicity).
For the Dirichlet part, we have
l
l ;h 0
l
l
=h l
l
l
MULTILEVEL BOUNDARY FUNCTIONALS FOR LEAST-SQUARES METHODS 9
which implies
l
h l
l
I j
l M
In order to set up the matrix entries associated with the Neumann boundary
conditions, we start from
l )
l
I j
l ]M l;\Gamma N
l I l
l
I j
l M l;\Gamma N
This implies
l )
~
with B l;\Gamma N
l I l
l
I j
l
which leads to
l ); (C \Gamma1=4
l
~
Clearly, from the quasiuniformity of the sequence of triangulations we obtain
l ); (C \Gamma1=4
l
~
~
which is the form we actually use in our implementation.
The expressions for the computation of the boundary contributions in (4.5) and
(4.6) involve the inverses of mass matrices M j;\Gamma D
and M j;\Gamma N
These inverses are dense matrices on the subspace of unknowns associated with \Gamma D
respectively. The actual computation of M \Gamma1
should be avoided,
especially for three-dimensional problems, and replaced by an inner conjugate gradient
iteration for solving linear systems with M j;\Gamma D and M j;\Gamma N . It is well-known that
the number of iterations required to achieve a certain accuracy is bounded independently
of l for these systems. Alternatively, one can use biorthogonal wavelets for
the construction of an H \Gamma1=2 (\Gamma N )-equivalent boundary functional for the Neumann
conditions. In the two-dimensional case, i.e., \Gamma N consists of a finite number of line
segments, this is described in [19, Section 5.6].
Fig. 5.1. Pressure (left) and flux (right) approximation for Example 10.5 1 ComputationalExperiments. For our computational experiments, we consider
on the rectangular domain shown in Figure 5.1. For the first
example, the boundary conditions are on the right two thirds of the upper and
on the right boundary segment, n \Delta on the left third of the upper boundary
segment and n \Delta segments. This can be viewed as a
simple model problem for flow in porous media where water infiltrates at a prescribed
rate on the left part of the upper boundary segement. The left and lower boundary
are impermeable and on the remaining part of the boundary pressure is set to zero
(hydraulic potential equals gravitational potential). The solution for the flux u is indicated
by the arrows in Figure 5.1. For this example, we have
3=2\Gamma"(\Omega\Gamma for all
due to the jump from Neumann to Dirichlet boundary condition on the upper
boundary segment.
The coarsest triangulation T 0 consists of 74 triangles which are then uniformly
refined. Figure 5.1 shows the first refinement T 1 . Table 5.1 shows the error with
respect to the "exact solution", measured in the H(div ; \Omega\Gamma \Theta H
As an "exact solution" we use the Galerkin approximation (with exact enforcement
of the boundary conditions on an adaptively refined triangulation T ?
4 which
consists of 24042 degrees of freedom for the edges and 8133 for the nodes. The triangulation
4 was constructed from T 3 using four steps of the refinement algorithm
presented in [20] based on the least-squares functional as an a posteriori error esti-
mator. The least-squares functional is reduced by a factor of more than 10 on T  compared to T 3 . This leads us to the conclusion that it is legitimate to use it as
reference solution in our experiments.

Table
Approximation error for the multilevel boundary functional for Example 1
MULTILEVEL BOUNDARY FUNCTIONALS FOR LEAST-SQUARES METHODS 11
Obviously, Table 5.1 shows that the approximation is slightly improved if the
boundary functional is weighted by instead of enforcing the boundary
conditions exactly 1). We need to compare the multilevel boundary functional
with the weighted L 2 (\Gamma) approach of [13] and [1]. In our situation, this means
replacing the least-squares functional (2.2) on V l \Theta W l by
kn
h l
This leads to the numbers listed in Table 5.2 which are slightly larger in most cases.
However, the difference between the multilevel and the weighted L 2 (\Gamma) functional is
rather marginal in this case.

Table
Approximation error for the weighted L 2 (\Gamma) boundary functional for Example 1
In our second example, the boundary conditions are chosen such that it is not
possible to satisfy them exactly with our finite element spaces. On the upper bound-
ary, we set linearly
interpolated for x 1 chosen to be smaller than the edge length of
our finest triangulation. The boundary conditions on the other boundary segments
are the same as in the first example. The numbers obtained with the multilevel and
the weighted L 2 boundary functionals are shown in Tables 5.3 and 5.4, respectively.

Table
Approximation error for the multilevel boundary functional for Example 2

Table
Approximation error for the weighted L 2 (\Gamma) boundary functional for Example 2
In conclusion, our computational results show the feasibility of using the multi-level
boundary functional approach. In particular, this method seems promising in
cases where it is not possible to enforce the boundary conditions exactly by the finite
element spaces. The construction and use of similar multilevel boundary functionals
for more complicated situations including nonlinear boundary conditions in porous
media applications is the focus of current research.

Acknowledgement

. I would like to thank Tom Manteuffel for many helpful discussions
related to the subject of this paper. The detailed comments of the anonymous
referees are also very much appreciated. In particular, I am thankful to one of them
for pointing out references [13] and [1] to me.



--R

Least squares methods for elliptic systems
Dynamics of Fluids in Porous Media
Cambridge University Press

A generalized Ritz-least-squares method for Dirichlet problems
Parallel multilevel preconditioners
The Mathematical Theory of Finite Element Methods
Mixed and Hybrid Finite Element Methods


Finite Element Methods for Navier-Stokes Equations
Elliptic Problems in Nonsmooth Domains
A least squares decomposition method for solving elliptic equations
Multilevel preconditioning - appending boundary conditions by Lagrange multi- pliers


Multilevel Finite Element Approximation



Two preconditioners based on the multi-level splitting of finite element spaces
--TR

--CTR
Huo-Yuan Duan , Guo-Ping Liang, Nonconforming elements in least-squares mixed finite element methods, Mathematics of Computation, v.73 n.245, p.1-18, January 2004
multilevel boundary functionals;least-squares mixed finite element method;raviart-thomas spaces