+<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 continuous 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 to
+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>
--- /dev/null
+\documentclass{article}
+\usepackage{amsmath}
+
+\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$},
+ \\
+ \partial_n u &= g && \text{on $\partial\Omega$}.
+\end{align*}
+It is well known that if this problem is to have a solution, then the forces
+need to satisfy the compatibility condition
+\begin{gather*}
+ \int_\Omega f\; dx + \int_{\partial\Omega} g\; ds = 0.
+\end{gather*}
+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.
+
+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.
+
+For this, there are various possibilities:
+\begin{itemize}
+\item 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
+ 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.
+
+\item 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 continuous 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
+ 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}
+We will choose the last possibility, since we want to demonstrate another
+technique with it.
+
+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 to
+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).
+
+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*}
+ \int_{\partial \Omega} u(x) \; ds = 0
+\end{gather*}
+can be written as
+\begin{gather*}
+ \sum_{i\in\partial\Omega_h} u_i = 0,
+\end{gather*}
+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*}
+ u_{i_0} = \sum_{i\in\partial\Omega_h\backslash i_0} -u_i.
+\end{gather*}
+This, luckily, is exactly the form of constraints for which the
+\texttt{ConstraintMatrix} 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
+of the form
+\begin{gather*}
+ CU = 0
+\end{gather*}
+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
+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$.
+
+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.
+
+But now on for the implementation of the program solving this problem...
+
+\end{document}
projects / dealii-svn.git / commitdiff
