-\documentclass{article}
-\usepackage{amsmath}
+<a name="Intro"></a>
+<h1>Introduction</h1>
-\begin{document}
The problem we will be considering is the solution of Laplace's problem with
Neumann boundary conditions only:
-\begin{align*}
- -\Delta u &= f &&\text{in $\Omega$},
+@f{eqnarray*}
+ -\Delta u &=& f \mathrm{in}\ \Omega,
\\
- \partial_n u &= g && \text{on $\partial\Omega$}.
-\end{align*}
+ \partial_n u &=& g \mathrm{on}\ \partial\Omega.
+@f}
It is well known that if this problem is to have a solution, then the forces
need to satisfy the compatibility condition
-\begin{gather*}
+@f[
\int_\Omega f\; dx + \int_{\partial\Omega} g\; ds = 0.
-\end{gather*}
+@f]
We will consider the special case that $\Omega$ is the circle of radius 1
around the origin, and $f=-2$, $g=1$. This choice satisfies the compatibility
condition.
which fixes this constant.
For this, there are various possibilities:
-\begin{itemize}
-\item Fix one node of the discretization to zero or any other fixed value.
+<ol>
+<li> Fix one node of the discretization to zero or any other fixed value.
This amounts to an additional condition $u_h(x_0)=0$. Although this is
common practice, it is not necessarily a good idea, since we know that the
solutions of Laplace's equation are only in $H^1$, which does not allow for
often see this in a resulting error spike at this point in the numerical
solution.
-\item Fixing the mean value over the domain to zero or any other value. This
+<li> Fixing the mean value over the domain to zero or any other value. This
is allowed on the continuous level, since $H^1(\Omega)\subset L^1(\Omega)$
by Sobolev's inequality, and thus also on the discrete level since we
there only consider subsets of $H^1$.
-\item Fixing the mean value over the boundary of the domain to zero or any
+<li> Fixing the mean value over the boundary of the domain to zero or any
other value. This is also allowed on the continuous level, since
$H^{1/2}(\partial\Omega)\subset L^1(\partial\Omega)$, again by Sobolev's
inequality.
-\end{itemize}
+</ol>
We will choose the last possibility, since we want to demonstrate another
technique with it.
consider the objects describing mappings as a black box which you need not
worry about, because their only uses seem to be to be passed to places deep
inside the library where functions know how to handle them (i.e. in the
-\texttt{FEValues} classes and their descendents).
+<code>FEValues</code> classes and their descendents).
The tricky point in this program is the use of the mean value
constraint. Fortunately, there is a class in the library which knows how to
handle such constraints, and we have used it quite often already, without
mentioning its generality. Note that if we assume that the boundary nodes are
spaced equally along the boundary, then the mean value constraint
-\begin{gather*}
+@f[
\int_{\partial \Omega} u(x) \; ds = 0
-\end{gather*}
+@f]
can be written as
-\begin{gather*}
+@f[
\sum_{i\in\partial\Omega_h} u_i = 0,
-\end{gather*}
+@f]
where the sum shall run over all degree of freedom indices which are located
on the boundary of the computational domain. Let us denote by $i_0$ that index
on the boundary with the lowest number (or any other conveniently chosen
index), then the constraint can also be represented by
-\begin{gather*}
+@f[
u_{i_0} = \sum_{i\in\partial\Omega_h\backslash i_0} -u_i.
-\end{gather*}
+@f]
This, luckily, is exactly the form of constraints for which the
-\texttt{ConstraintMatrix} class was designed. Note that we have used this
+<code>ConstraintMatrix</code> class was designed. Note that we have used this
class in several previous examples for the representation of hanging nodes
constraints, which also have this form: there, the middle vertex shall have
the mean of the values of the adjacent vertices. In general, the
-\texttt{ConstraintMatrix} class is designed to handle homogeneous constraints
+<code>ConstraintMatrix</code> class is designed to handle homogeneous constraints
of the form
-\begin{gather*}
+@f[
CU = 0
-\end{gather*}
+@f]
where $C$ denotes a matrix, and $U$ the vector of nodal values.
In this example, the mean value along the boundary allows just such a
representation, with $C$ being a matrix with just one row (i.e. there is only
one constraint). In the implementation, we will create a
-\texttt{ConstraintMatrix} object, add one constraint (i.e. add another row to
+<code>ConstraintMatrix</code> object, add one constraint (i.e. add another row to
the matrix) referring to the first boundary node $i_0$, and insert the weights
with which all the other nodes contribute, which in this example happens to be
just $-1$.
of the rows of the matrix. We will show how we can use an intermediate object
to work around this problem.
-But now on for the implementation of the program solving this problem...
-
-\end{document}
+But now on to the implementation of the program solving this problem...
+++ /dev/null
-<a name="Intro"></a>
-<h1>Introduction</h1>
-
-<p>
-In this example, we demonstrate solving a simple Laplace problem using
-higher order mappings, as well as using some stranger sort of
-constraints and a variant of the <code>SparsityPattern</code>
-class. As was noted in the introduction to the <a href="step-10.html"
-target="body">step-10</a> example program, it is difficult to find
-problems where it actually makes a difference in the convergence order
-whether we use simple linear mappings of the unit cells to the cells
-in real space, or higher order mappings. Such problems are not exotic
-(in fact, only the work with such problems generated the interest in
-implementing higher order mappings in deal.II), and include for
-example the numerical solution of the Euler equations of inviscid gas
-flow, but they are too complicated to be made the subject of an
-example program. We therefore restrict our attention to a problem
-which exhibits the same order of convergence, but profits
-significantly from higher order mappings nevertheless, by a reduced
-size of the error.
-</p>
-
-<p>
-The problem we will be considering is the solution of Laplace's problem with
-Neumann boundary conditions only:
-<BR>
-<IMG
- WIDTH="263" HEIGHT="40" ALIGN="BOTTOM" BORDER="0"
- SRC="step-11.data/img1.gif"
- ALT="\begin{align*}-\Delta u &= f &&\text{in $\Omega$ },
-\\
-\partial_n u &= g && \text{on $\partial\Omega$ }.
-\end{align*}">
-<BR>
-It is well known that if this problem is to have a solution, then the forces
-need to satisfy the compatibility condition
-<BR>
-<IMG
- WIDTH="162" HEIGHT="37" ALIGN="BOTTOM" BORDER="0"
- SRC="step-11.data/img2.gif"
- ALT="\begin{gather*}\int_\Omega f\; dx + \int_{\partial\Omega} g\; ds = 0.
-\end{gather*}">
-<BR>
-We will consider the special case that <IMG
- WIDTH="15" HEIGHT="14" ALIGN="BOTTOM" BORDER="0"
- SRC="step-11.data/img3.gif"
- ALT="$\Omega$">
-is the circle of radius 1
-around the origin, and <I>f</I>=-2, <I>g</I>=1. This choice satisfies the compatibility
-condition.
-
-<P>
-The compatibility condition allows a solution of the above equation, but it
-nevertheless retains an ambiguity: since only derivatives of the solution
-appear in the equations, the solution is only determined up to a constant. For
-this reason, we have to pose another condition for the numerical solution,
-which fixes this constant.
-
-<P>
-For this, there are various possibilities:
-<UL>
-<LI>Fix one node of the discretization to zero or any other fixed value.
- This amounts to an additional condition
-<!-- MATH: $u_h(x_0)=0$ -->
-<I>u</I><SUB><I>h</I></SUB>(<I>x</I><SUB>0</SUB>)=0. Although this is
- common practice, it is not necessarily a good idea, since we know that the
- solutions of Laplace's equation are only in <I>H</I><SUP>1</SUP>, which does not allow for
- the definition of point values because it is not a subset of the continuous
- functions. Therefore, even though fixing one node is allowed for
- discretitized functions, it is not for continuous functions, and one can
- often see this in a resulting error spike at this point in the numerical
- solution.
-<LI>Fixing the mean value over the domain to zero or any other value. This
- is allowed on the continuous level, since
-<!-- MATH: $H^1(\Omega)\subset L^1(\Omega)$ -->
-<IMG
- WIDTH="112" HEIGHT="33" ALIGN="MIDDLE" BORDER="0"
- SRC="step-11.data/img4.gif"
- ALT="$H^1(\Omega)\subset L^1(\Omega)$">
-by Sobolev's inequality, and thus also on the discrete level since we
- there only consider subsets of <I>H</I><SUP>1</SUP>.
-<LI>Fixing the mean value over the boundary of the domain to zero or any
- other value. This is also allowed on the continuous level, since
-
-<!-- MATH: $H^{1/2}(\partial\Omega)\subset L^1(\partial\Omega)$ -->
-<IMG
- WIDTH="143" HEIGHT="36" ALIGN="MIDDLE" BORDER="0"
- SRC="step-11.data/img5.gif"
- ALT="$H^{1/2}(\partial\Omega)\subset L^1(\partial\Omega)$">,
-again by Sobolev's
- inequality.
-</UL>We will choose the last possibility, since we want to demonstrate another
-technique with it.
-
-<P>
-While this describes the problem to be solved, we still have to figure out how
-to implement it. Basically, except for the additional mean value constraint,
-we have solved this problem several times, using Dirichlet boundary values,
-and we only need to drop the treatment of Dirichlet boundary nodes. The use of
-higher order mappings is also rather trivial and will be explained at the
-various places where we use it; in almost all conceivable cases, you will only
-consider the objects describing mappings as a black box which you need not
-worry about, because their only uses seem to be to be passed to places deep
-inside the library where functions know how to handle them (i.e. in the
-<TT>FEValues</TT> classes and their descendents).
-
-<P>
-The tricky point in this program is the use of the mean value
-constraint. Fortunately, there is a class in the library which knows how to
-handle such constraints, and we have used it quite often already, without
-mentioning its generality. Note that if we assume that the boundary nodes are
-spaced equally along the boundary, then the mean value constraint
-<BR>
-<IMG
- WIDTH="109" HEIGHT="37" ALIGN="BOTTOM" BORDER="0"
- SRC="step-11.data/img7.gif"
- ALT="\begin{gather*}\int_{\partial \Omega} u(x) \; ds = 0
-\end{gather*}">
-<BR>
-can be written as
-<BR>
-<IMG
- WIDTH="87" HEIGHT="37" ALIGN="BOTTOM" BORDER="0"
- SRC="step-11.data/img8.gif"
- ALT="\begin{gather*}\sum_{i\in\partial\Omega_h} u_i = 0,
-\end{gather*}">
-<BR>
-where the sum shall run over all degree of freedom indices which are located
-on the boundary of the computational domain. Let us denote by <I>i</I><SUB>0</SUB> that index
-on the boundary with the lowest number (or any other conveniently chosen
-index), then the constraint can also be represented by
-<BR>
-<IMG
- WIDTH="130" HEIGHT="38" ALIGN="BOTTOM" BORDER="0"
- SRC="step-11.data/img9.gif"
- ALT="\begin{gather*}u_{i_0} = \sum_{i\in\partial\Omega_h\backslash i_0} -u_i.
-\end{gather*}">
-<BR>
-This, luckily, is exactly the form of constraints for which the
-<TT>ConstraintMatrix</TT> class was designed. Note that we have used this
-class in several previous examples for the representation of hanging nodes
-constraints, which also have this form: there, the middle vertex shall have
-the mean of the values of the adjacent vertices. In general, the
-<TT>ConstraintMatrix</TT> class is designed to handle homogeneous constraints
-of the form
-<BR>
-<IMG
- WIDTH="55" HEIGHT="13" ALIGN="BOTTOM" BORDER="0"
- SRC="step-11.data/img10.gif"
- ALT="\begin{gather*}CU = 0
-\end{gather*}">
-<BR>
-where <I>C</I> denotes a matrix, and <I>U</I> the vector of nodal values.
-
-<P>
-In this example, the mean value along the boundary allows just such a
-representation, with <I>C</I> being a matrix with just one row (i.e. there is only
-one constraint). In the implementation, we will create a
-<TT>ConstraintMatrix</TT> object, add one constraint (i.e. add another row to
-the matrix) referring to the first boundary node <I>i</I><SUB>0</SUB>, and insert the weights
-with which all the other nodes contribute, which in this example happens to be
-just -1.
-
-<P>
-Later, we will use this object to eliminate the first boundary node from the
-linear system of equations, reducing it to one which has a solution without
-the ambiguity of the constant shift value. One of the problems of the
-implementation will be that the explicit elimination of this node results in a
-number of additional elements in the matrix, of which we do not know in
-advance where they are located and how many additional entries will be in each
-of the rows of the matrix. We will show how we can use an intermediate object
-to work around this problem.
-
-<P>
-But now on for the implementation of the program solving this problem...
-</p>