--- /dev/null
+\documentclass{article}
+
+\begin{document}
+
+In this example program, we will consider two aspects:
+\begin{itemize}
+\item Verification of correctness of the program;
+\item Nonhomogeneous Neumann boundary conditions for the Helmholtz equation.
+\end{itemize}
+
+\paragraph{Verification of correctness.} There has probably never been a
+non-trivial finite element program that worked right from the start. It is
+therefore necessary to find ways to verify whether a computed solution is
+correct or not. Usually, this is done by choosing the set-up of a simulation
+such that we know the exact continuous solution and evaluate the difference
+between continuous and computed discrete solution. If this difference
+converges to zero with the right order of convergence, this is already a good
+indication of correctness, although there may be other sources of error
+persisting which have only a small contribution to the total error or are of
+higher order.
+
+In this example, we will not go into the theories of systematic software
+verification which is a very complicated problem. Rather we will demonstrate
+the tools which deal.II can offer in this respect. This is basically centered
+around the functionality of a single function, \texttt{integrate\_difference}.
+This function computes the difference between a given continuous function and
+a finite element field in various norms on each cell. At present, the
+supported norms are the following, where $u$ denotes the continuous function
+and $u_h$ the finite element field, and $K$ is an element of the
+triangulation:
+\begin{eqnarray*}
+ {\| u-u_h \|}_{L_1(K)} &=& \int_K |u-u_h| \; dx,
+ \\
+ {\| u-u_h \|}_{L_2(K)} &=& \left( \int_K |u-u_h|^2 \; dx \right)^{1/2},
+ \\
+ {\| u-u_h \|}_{L_\infty(K)} &=& \max_{x \in K} |u(x) - u_h(x)|,
+ \\
+ {| u-u_h |}_{H^1(K)} &=& \left( \int_K |\nabla(u-u_h)|^2 \; dx \right)^{1/2},
+ \\
+ {\| u-u_h \|}_{H^1(K)} &=& \left( {\| u-u_h \|}^2_{L_2(K)}
+ +{| u-u_h |}^2_{H^1(K)} \right)^{1/2}.
+\end{eqnarray*}
+All these norms and seminorms can also be evaluated with weighting functions,
+for example in order to exclude singularities from the determination of the
+global error. The function also works for vector-valued functions. It should
+be noted that all these quantities are evaluated using quadrature formulas;
+the choice of the right quadrature formula is therefore crucial to the
+accurate evaluation of the error. This holds in particual for the $L_\infty$
+norm, where we evaluate the maximal deviation of numerical and exact solution
+only at the quadrature points; one should then not try to use a quadrature
+rule with points only at points where superconvergence might occur.
+
+The function \texttt{integrate\_difference} evaluates the desired norm on each
+cell $K$ of the triangulation and returns a vector which holds these
+values. From the local values, we can then obtain the global error. For
+example, if the vector $(e_i)$ contains the local $L_2$ norms, then
+$$
+ E = \| {\mathbf e} \| = \left( \sum_i e_i^2 \right)^{1/2}
+$$
+is the global $L_2$ error.
+
+In the program, we will show how to evaluate and use these quantities, and we
+will monitor their values under mesh refinement. Of course, we have to choose
+the problem at hand such that we can explicitely state the solution and its
+derivatives, but since we want to evaluate the correctness of the program,
+this is only reasonable. If we know that the program produces the correct
+solution for one (or, if one wants to be really sure: many) specifically
+chosen right hand sides, we can be rather confident that it will also compute
+the correct solution for problems where we don't know the exact values.
+
+
+\paragraph{Non-homogeneous Neumann boundary conditions.} The second, totally
+unrelated, subject of this example program is the use of non-homogeneous
+boundary conditions. These are included into the variational form using
+boundary integrals which we have to evaluate numerically when assembling the
+right hand side vector.
+
+Before we go into programming, let's have a brief look at the mathematical
+formulation. The equation which we want to solve is Helmholtz's equation
+``with the nice sign'':
+$$
+ -\Delta u + u = f,
+$$
+on the square $[-1,1]^2$, augmented by boundary conditions
+$$
+ u = g_1
+$$
+on some part $\Gamma_1$ of the boundary $\Gamma$, and
+$$
+ {\mathbf n}\cdot \nabla u = g_2
+$$
+on the rest $\Gamma_2 = \Gamma \backslash \Gamma_1$.
+
+We choose the right hand side function $f$ such that the exact solution is
+$$
+ u(x) = \sum_{i=1}^3 \exp\left(-\frac{|x-x_i|^2}{\sigma^2}\right)
+$$
+where the centers $x_i$ of the exponentials are
+ $x_1=(-\frac 12,\frac 12)$,
+ $x_2=(-\frac 12,-\frac 12)$, and
+ $x_3=(\frac 12,-\frac 12)$.
+The half width is set to $\sigma=\frac 13$.
+
+We further choose $\Gamma_1=\Gamma \cap\{\{x=1\} \cup \{y=1\}\}$, and there
+set $g_1$ such that it resembles the exact values of $u$. Likewise, we choose
+$g_2$ on the remaining portion of the boundary to be the exact normal
+derivatives of the continuous solution.
+
+Using the above definitions, we can state the weak formulation of the
+equation, which reads: find $u\in H^1_g=\{v\in H^1: v|_{\Gamma_1}=g_1\}$ such
+that
+$$
+ {(\nabla u, \nabla v)}_\Omega + {(u,v)}_\Omega
+ =
+ {(f,v)}_\Omega + {(g_2,v)}_{\Gamma_2}
+$$
+for all test functions $v\in H^1_0=\{v\in H^1: v|_{\Gamma_1}=0\}$. The
+boundary term ${(g_2,v)}_{\Gamma_2}$ has appeared by partial integration and
+using $\partial_n u=g$ on $\Gamma_2$ and $v=0$ on $\Gamma_1$. The cell
+matrices and vectors which we use to build the global matrices and right hand
+side vectors in the discrete formulation therefore look like this:
+\begin{eqnarray*}
+ A_{ij}^K &=& \left(\nabla \varphi_i, \nabla \varphi_j\right)_K
+ +\left(\varphi_i, \varphi_j\right)_K,
+ \\
+ f_i^K &=& \left(f,\varphi_i\right)_K
+ +\left(g_2, \varphi_i\right)_{\partial K\cap \Gamma_2}.
+\end{eqnarray*}
+Since the generation of the domain integrals has been shown in previous
+examples several times, only the generation of the contour integral is of
+interest here. It basically works along the following lines: for domain
+integrals we have the \texttt{FEValues} class that provides values and
+gradients of the shape values, as well as Jacobian determinants and other
+information and specified quadrature points in the cell; likewise, there is a
+class \texttt{FEFaceValues} that performs these tasks for integrations on
+faces of cells. One provides it with a quadrature formula for a mannifold with
+dimension one less than the dimension of the domain is, and the cell and the
+number of its face on which we want to perform the integration. The class will
+then compute the values, gradients, normal vectors, weights, etc at the
+quadrature points on this face, which we can then use in the same way as for
+the domain integrals. The details of how this is done are shown in the
+following program.
+
+
+\end{document}
+
+%%% Local Variables:
+%%% mode: latex
+%%% TeX-master: t
+%%% End:
projects / dealii-svn.git / commitdiff