--- /dev/null
+\documentclass{article}
+\usepackage{amsmath}
+\begin{document}
+
+\section{The maths}
+
+The Heidelberg group of Professor Rolf Rannacher, to which the three main
+authors of the deal.II library belonged for the PhD time and partly also
+afterwards, has been involved with adaptivity and error estimation for finite
+element discretizations since the mid-90ies. The main achievement is the
+development of error estimates for arbitrary functionals of the solution, and
+of optimal mesh refinement for its computation.
+
+We will not discuss the derivation of these concepts in too great detail, but
+will implement the main ideas in the present example program. For a thorough
+introduction into the general idea, we refer to the seminal work of Becker and
+Rannacher \cite{BR95,BR96r}, and the overview article of the same authors in
+Acta Numerica \cite{BR01}; the first introduces the concept of error
+estimation and adaptivity for general functional output for the Laplace
+equation, while the second gives many examples of applications of these
+concepts to a large number of other, more complicated equations. For
+applications to individual types of equations, see also the publications by
+Becker \cite{Bec95,Bec98}, Kanschat \cite{Kan96,FK97}, Suttmeier
+\cite{Sut96,RS97,RS98c,RS99}, Bangerth \cite{BR99b,Ban00w,BR01a,Ban02}, and
+Hartmann \cite{HH01,HH01a,HH01b}.
+
+The basic idea is the following: in applications, one is not usually
+interested in the solution per se, but rather in certain aspects of it. For
+example, in simulations of flow problems, one may want to know the lift or
+drag of a body emersed in the fluid; it is this quantity that we want to know
+to best accuracy, and whether the rest of the solution of the describing
+equations is well resolved is not of primary interest. Likewise, in elasticity
+one might want to know about values of the stress at certain points to guess
+whether maximal load values of joints are safe, for example. Or, in radiative
+transfer problems, mean flux intensities are of interest.
+
+In all the cases just listed, it is the evaluation of a functional $J(u)$ of
+the solution which we are interested in, rather than the values of $u$
+everywhere. Since the exact solution $u$ is not available, but only its
+numerical approximation $u_h$, it is sensible to ask whether the computed
+value $J(u_h)$ is within certain limits of the exact value $J(u)$, i.e. we
+want to bound the error with respect to this functional, $J(u)-J(u_h)$.
+
+For simplicity of exposition, we henceforth assume that both the quantity of
+interest $J$, as well as the equation are linear, and we will in particular
+show the derivation for the Laplace equation with homogeneous Dirichlet
+boundary conditions, although the concept is much more general. For this
+general case, we refer to the references listed above. The goal is to obtain
+bounds on the error, $J(e)=J(u)-J(u_h)$. For this, let us denote by $z$ the
+solution of a dual problem, defined as follows:
+\begin{gather}
+ a(\varphi,z) = J(\varphi) \qquad \forall \varphi,
+\end{gather}
+where $a(\cdot,\cdot)$ is the bilinear form associated with the differential
+equation, and the test functions are chosen from the corresponding solution
+space. Then, taking as special test function $\varphi=e$ the error, we have
+that
+\begin{gather}
+ J(e) = a(e,z)
+\end{gather}
+and we can, by Galerkin orthogonality, rewrite this as
+\begin{gather}
+ J(e) = a(e,z-\varphi_h)
+\end{gather}
+for all possible functions $\varphi_h$ from the discrete test space.
+
+Concretely, for Laplace's equation, the error identity reads
+\begin{gather}
+ J(e) = (\nabla e, \nabla(z-\varphi_h)).
+\end{gather}
+For reasons that we will not explain, we do not want to use this formula as
+is, but rather split the scalar products into terms on all cells, and
+integrate by parts on each of them:
+\begin{align*}
+ J(e)
+ &=
+ \sum_K (\nabla (u-u_h), \nabla (z-\varphi_h))_K
+ \\
+ &=
+ \sum_K (-\Delta (u-u_h), z-\varphi_h)_K
+ + (\partial_n (u-u_h), z-z_h)_{\partial K}.
+\end{align*}
+Next we use that $-\Delta u=f$, and that $\partial_n u$ is a quantity that is
+continuous almost everywhere, so the terms involving $\partial_n u$ on one
+cell cancels with that on its neighbor, where the normal vector has the
+opposite sign. At the boundary of the domain, where there is no neighbor cell
+with wich this term could cancel, the weight $z-\varphi_h$ can be chosen as
+zero, since $z$ has zero boundary values, and $\varphi_h$ can be chosen to
+have the same.
+
+Thus, we have
+\begin{align*}
+ J(e)
+ &=
+ \sum_K (f+u_h), z-\varphi_h)_K
+ - (\partial_n u_h, z-\varphi_h)_{\partial K\backslash \partial\Omega}.
+\end{align*}
+In a final step, note that when taking the normal derivative of $u_h$, we mean
+the value of this quantity as taken from this side of the cell (for the usual
+Lagrange elements, derivatives are not continuous across edges). We then
+rewrite the above formula by exchanging half of the edge integral of cell $K$
+with the neighbor cell $K'$, to obtain
+\begin{align*}
+ J(e)
+ &=
+ \sum_K (f+u_h), z-\varphi_h)_K
+ - \frac 12 (\partial_n u_h|_K + \partial_{n'} u_h|_{K'},
+ z-\varphi_h)_{\partial K\backslash \partial\Omega}.
+\end{align*}
+Using that for the normal vectors $n'=-n$ holds, we define the jump of the
+normal derivative by
+\begin{gather*}
+ [\partial_n u_h] := \partial_n u_h|_K + \partial_{n'} u_h|_{K'}
+ =
+ \partial_n u_h|_K - \partial_n u_h|_{K'},
+\end{gather*}
+and get the final form after setting the discrete function $\varphi_h$, which
+is by now still arbitrary, to the point interpolation of the dual solution,
+$\varphi_h=I_h z$:
+\begin{align*}
+ J(e)
+ &=
+ \sum_K (f+u_h), z-I_h z)_K
+ - \frac 12 ([\partial_n u_h],
+ z-I_h z)_{\partial K\backslash \partial\Omega}.
+\end{align*}
+
+With this, we have obtained an exact representation of the error of the finite
+element discretization with respect to arbitrary (linear) functionals
+$J(\cdot)$. Its structure is a weighted form of a residual estimator, as both
+$f+\Delta u_h$ and $[\partial_n u_h]$ are cell and edge residuals that vanish
+on the exact solution, and $z-I_h z$ are weights indicating how important the
+residuals on a certain cell is for the evaluation of the given functional.
+Furthermore, it is a cell-wise quantity, so we can use it as a mesh refinement
+criterion. The question, is: how to evaluate it? After all, the evaluation
+requires knowledge of the dual solution $z$, which carries the information
+about the quantity we want to know to best accuracy.
+
+In some, very special cases, this dual solution is known. For example, if the
+functional $J(\cdot)$ is the point evaluation, $J(\varphi)=\varphi(x_0)$, then
+the dual solution has to satisfy
+\begin{gather*}
+ -\Delta z = \delta(x-x_0),
+\end{gather*}
+with the Dirac delta function on the right hand side, and the dual solution is
+the Green's function with respect to the point $x_0$. For simple geometries,
+this function is analytically known, and we could insert it into the error
+representation formula.
+
+However, we do not want to restrict ourselves to such special cases. Rather,
+we will compute the dual solution numerically, and approximate $z$ by some
+numerically obtained $\tilde z$. We note that it is not sufficient to compute
+this approximation $\tilde z$ using the same method as used for the primal
+solution $u_h$, since then $\tilde z-I_h \tilde z=0$, and the overall error
+estimate would be zero. Rather, the approximation $\tilde z$ has to be from a
+larger space than the primal finite element space. There are various ways to
+obtain such an approximation (see the cited literature), and we will choose to
+compute it with a higher order finite element space. While this is certainly
+not the most efficient way, it is simple since we already have all we need to
+do that in place, and it also allows for simple experimenting. For more
+efficient methods, again refer to the given literature, in particular
+\cite{BR95,BR96r,BR01}.
+
+With this, we end the discussion of the mathematical side of this program and
+turn to the actual implementation.
+
+
+\section{The software}
+
+
+\bibliographystyle{plain}
+\bibliography{own}
+
+\end{document}
+
+
+
+
+
important for computing the point value.
</p>
+
<p>
The next point is to compare the new (duality based) mesh refinement
criterion with the old ones. These are the results:
</p>
-<h2>Outlook</h2>
+
+<h2>Conclusions and outlook</h2>
<p>
-As stated, the program is quite modular, and implementing another test
-case, or another evaluation and dual functional is simple. You are
-encouraged to take the program as a basis for your own experiments,
-and to play a little.
-</p>
\ No newline at end of file
+The results here are not too clearly indicating the superiority of the
+dual weighted error estimation approach for mesh refinement over other
+mesh refinement criteria, such as the Kelly indicator. This is due to
+the relative simplicity of the shown application. If you are not
+convinced yet that this approach is indeed superior, you are invited
+to browse through the literature indicated in the introduction, where
+plenty of examples are provided where the dual weighted approach can
+reduce the necessary numerical work by orders of magnitude, making
+this the only way to compute certain quantities to reasonable
+accuracies at all.
+</p>
+
+<p>
+Besides the objections you may raise against its use as a mesh
+refinement criterion, consider that accurate knowledge of the error in
+the quantity one might want to compute is of great use, since we can
+stop computations when we are satisfied with the accuracy. Using more
+traditional approaches, it is very difficult to get accurate estimates
+for arbitrary quantities, except for, maybe, the error in the energy
+norm, and we will then have no guarantee that the result we computed
+satisfies any requirements on its accuracy.
+</p>
+
+<p>
+Leaving these mathematical considerations, we tried to write the
+program in a modular way, such that implementing another test case, or
+another evaluation and dual functional is simple. You are encouraged
+to take the program as a basis for your own experiments, and to play a
+little.
+</p>
+
+
+
+
projects / dealii-svn.git / commitdiff