<TABLE WIDTH="60%" ALIGN="center">
<tr>
<td ALIGN="center">
- @image html step-1.grid-1.jpg
+ @image html step-1.grid-1.png
</td>
<td ALIGN="center">
- @image html step-1.grid-2.jpg
+ @image html step-1.grid-2.png
</td>
</tr>
</table>
<a name="Intro"></a>
<h1>Introduction</h1>
-<p>
-<acronym>deal.II</acronym> has a unique feature which we call
+
+deal.II has a unique feature which we call
``dimension independent programming''. You may have noticed in the
previous examples that many classes had a number in angle brackets
suffixed to them. This is to indicate that for example the
triangulation in two and three space dimensions are different, but
related data types. We could as well have called them
<code>Triangulation2d</code> and <code>Triangulation3d</code> instead
-of <code>Triangulation<2></code> and
-<code>Triangulation<3></code> to name the two classes, but this
+of <code>Triangulation@<2@></code> and
+<code>Triangulation@<3@></code> to name the two classes, but this
has an important drawback: assume you have a function which does
exactly the same functionality, but on 2d or 3d triangulations,
depending on which dimension we would like to solve the equation in
<code>Triangulation3d</code>. This is an unnecessary obstacle in
programming and leads to a nuisance to keep the two function in synch
(at best) or difficult to find errors if the two versions get out of
-synch (at worst; this would probably the more common case).
-</p>
+sync (at worst; this would probably the more common case).
+
+
-<p>
Such obstacles can be circumvented by using some template magic as
provided by the C++ language: templatized classes and functions are
not really classes or functions but only a pattern depending on an
the template can depend on the template parameters, and they will be
resolved at the time of compilation for a specific template
parameter. For example, consider the following piece of code:
-<pre><code>
- template <int dim>
- void make_grid (Triangulation<dim> &triangulation)
+@code
+ template <int dim>
+ void make_grid (Triangulation<dim> &triangulation)
{
GridGenerator::hyper_cube (triangulation, -1, 1);
};
-</code></pre>
-</p>
+@endcode
+
+
-<p>
At the point where the compiler sees this function, it does not know
-anything about the actual value of ``dim''. The only thing the compiler has is
+anything about the actual value of <code>dim</code>. The only thing the compiler has is
a template, i.e. a blueprint, to generate
functions <code>make_grid</code> if given a particular value of
-``dim''. Since ``dim'' has an unknown value, there is no
+<code>dim</code>. Since <code>dim</code> has an unknown value, there is no
code the compiler can generate for the moment.
-</p>
-<p>
+
+
However, if later down the compiler would encounter code that looks, for
example, like this,
-<pre><code>
- Triangulation<2> triangulation;
+@code
+ Triangulation<2> triangulation;
make_grid (triangulation);
-</code></pre>
+@endcode
then the compiler will deduce that the function <code>make_grid</code> for
<code>dim==2</code> was
requested and will compile the template above into a function with dim replaced
by 2 everywhere, i.e. it will compile the function as if it were defined
as
-<pre><code>
- void make_grid (Triangulation<2> &triangulation)
+@code
+ void make_grid (Triangulation<2> &triangulation)
{
GridGenerator::hyper_cube (triangulation, -1, 1);
};
-</code></pre>
-</p>
+@endcode
+
+
-<p>
However, it is worth to note that the function
<code>GridGenerator::hyper_cube</code> depends on the dimension as
well, so in this case, the compiler will call the function
-<code>GridGenerator::hyper_cube<2></code> while if dim were 3,
-it would call <code>GridGenerator::hyper_cube<3></code> which
+<code>GridGenerator::hyper_cube@<2@></code> while if dim were 3,
+it would call <code>GridGenerator::hyper_cube@<3@></code> which
might be (and actually is) a totally unrelated function.
-</p>
-<p>
+
+
The same can be done with member variables. Consider the following
function, which might in turn call the above one:
-<pre><code>
- template <int dim>
- void make_grid_and_dofs (Triangulation<dim> &triangulation)
+@code
+ template <int dim>
+ void make_grid_and_dofs (Triangulation<dim> &triangulation)
{
make_grid (triangulation);
- DoFHandler<dim> dof_handler(triangulation);
+ DoFHandler<dim> dof_handler(triangulation);
...
};
-</code></pre>
+@endcode
This function has a member variable of type
-<code>DoFHandler<dim></code>. Again, the compiler can't
+<code>DoFHandler@<dim@></code>. Again, the compiler can't
compile this function until it knows for which dimension. If you call
this function for a specific dimension as above, the compiler will
take the template, replace all occurences of dim by the dimension for
which it was called, and compile it. If you call the function several
times for different dimensions, it will compile it several times, each
-time calling the right ``make_grid'' function and reserving the right
+time calling the right <code>make_grid</code> function and reserving the right
amount of memory for the member variable; note that the size of a
<code>DoFHandler</code> might, and indeed does, depend on the space dimension.
-</p>
-<p>
-The <acronym>deal.II</acronym> library is build around this concept
+
+
+The deal.II library is build around this concept
of dimension-independent programming, and therefore allows you to program in
a way that will not need to
distinguish between the space dimensions. It should be noted that in
only a very few places is it necessary to actually compare the
-dimension using ``if''s or ``switch''es. However, since the compiler
+dimension using <code>if</code>s or <code>switch</code>es. However, since the compiler
has to compile each function for each dimension separately, even there
-it knows the value of ``dim'' at the time of compilation and will
-therefore be able to optimize away the ``if'' statement along with the
+it knows the value of <code>dim</code> at the time of compilation and will
+therefore be able to optimize away the <code>if</code> statement along with the
unused branch.
-</p>
-<p>
+
+
In this example program, we will show how to program dimension
independently (which in fact is even simpler than if you had to take
care about the dimension) and we will extend the Laplace problem of
the last example to a program that runs in two and three space
dimensions at the same time. Other extensions are the use of a
non-constant right hand side function and of non-zero boundary values.
-</p>
+
<a name="Intro"></a>
<h1>Introduction</h1>
-<p>
+
This example does not show revolutionary new things, but it shows many
small improvements over the previous examples, and also many small
things that can usually be found in finite element programs. Among
preconditioned iterative solvers for the linear systems of
equations.
</ul>
-</p>
Number of degrees of freedom: 20609
182 CG iterations needed to obtain convergence.
--------------------------------------------------------
-An error occurred in line <273> of file <step-5.cc> in function
- void Coefficient<dim>::value_list(const std::vector<Point<dim>, std::allocator<Point<dim> > >&, std::vector<double, std::allocator<double> >&, unsigned int)
+An error occurred in line <273> of file <step-5.cc> in function
+ void Coefficient<dim>::value_list(const std::vector<Point<dim>, std::allocator<Point<dim> > >&, std::vector<double, std::allocator<double> >&, unsigned int)
const [with int dim = 2]
The violated condition was:
values.size() == points.size()
Stacktrace:
-----------
-#0 ./step-5: Coefficient<2>::value_list(std::vector<Point<2>, std::allocator<Point<2> > > const&, std::vector<double, std::allocator<double> >&, unsigned) const
+#0 ./step-5: Coefficient<2>::value_list(std::vector<Point<2>, std::allocator<Point<2> > > const&, std::vector<double, std::allocator<double> >&, unsigned) const
#1 ./step-5: main
--------------------------------------------------------
make: *** [run] Aborted
As for the error — let's look at it again:
@code
--------------------------------------------------------
-An error occurred in line <273> of file <step-5.cc> in function
- void Coefficient<dim>::value_list(const std::vector<Point<dim>, std::allocator<Point<dim> > >&, std::vector<double, std::allocator<double> >&, unsigned int)
+An error occurred in line <273> of file <step-5.cc> in function
+ void Coefficient<dim>::value_list(const std::vector<Point<dim>, std::allocator<Point<dim> > >&, std::vector<double, std::allocator<double> >&, unsigned int)
const [with int dim = 2]
The violated condition was:
values.size() == points.size()
Stacktrace:
-----------
-#0 ./step-5: Coefficient<2>::value_list(std::vector<Point<2>, std::allocator<Point<2> > > const&, std::vector<double, std::allocator<double> >&, unsigned) const
+#0 ./step-5: Coefficient<2>::value_list(std::vector<Point<2>, std::allocator<Point<2> > > const&, std::vector<double, std::allocator<double> >&, unsigned) const
#1 ./step-5: main
--------------------------------------------------------
make: *** [run] Aborted
<a name="Intro"></a>
<h1>Introduction</h1>
-<p>
+
The main emphasis in this example is the handling of locally refined
grids. The approach to adaptivity chosen in deal.II us to use grids in which
neighboring cells may be refined a different number of times. This then
results in nodes on the interfaces of cells which belong to one
side, but are unbalanced on the other. The common term for these is
-“hanging nodes”.
-</p>
+“hanging nodes”.
+
+
-<p>
To guarantee that the global solution is continuous at these nodes as
well, we have to state some additional constraints on the values of
the solution at these nodes. In the program below, we will show how we
can get these constraints from the library, and how to use them in the
solution of the linear system of equations.
-</p>
-<p>
+
+
The locally refined grids are produced using an error estimator class
which estimates the energy error with respect to the Laplace
operator. This error estimator, although developed for Laplace's
equations, it is often a good way to quickly produce meshes that are
well adapted to the features of solutions, such as regions of great
variation or discontinuities. Since it was developed by Kelly and
-co-workers, we often refer to it as the “Kelly refinement
-indicator” in the library, documentation, and mailing list. The
+co-workers, we often refer to it as the “Kelly refinement
+indicator” in the library, documentation, and mailing list. The
class that implements it is called
-“KellyErrorEstimator”. Although the error estimator (and
+<code>KellyErrorEstimator</code>. Although the error estimator (and
its
implementation in the deal.II library) is capable of handling variable
coefficients in the equation, we will not use this feature since we
are only interested in a quick and simple way to generate locally
refined grids.
-</p>
-<p>
+
+
Since the concepts used for locally refined grids are so important,
we do not show much additional new stuff in this example. The most
important exception is that we show how to use biquadratic elements
instead of the bilinear ones which we have used in all previous
examples. In fact, The use of higher order elements is accomplished by
only replacing three lines of the program, namely the declaration of
-the ``fe'' variable, and the use of an appropriate quadrature formula
+the <code>fe</code> variable, and the use of an appropriate quadrature formula
in two places. The rest of the program is unchanged.
-</p>
-<p>
+
+
The only other new thing is a method to catch exceptions in the
-“main” function in order to output some information in case the
+<code>main</code> function in order to output some information in case the
program crashes for some reason.
-</p>
+
@code
--------------------------------------------------------
-An error occurred in line <79> of file <source/subscriptor.cc> in function
+An error occurred in line <79> of file <source/subscriptor.cc> in function
virtual Subscriptor::~Subscriptor()
The violated condition was:
counter == 0
The name and call sequence of the exception was:
- ExcInUse(counter, object_info->name(), infostring)
+ ExcInUse(counter, object_info->name(), infostring)
Additional Information:
Object of class 4FE_QILi2EE is still used by 1 other objects.
from Subscriber 10DoFHandlerILi2EE
Stacktrace:
-----------
#0 /u/bangerth/p/deal.II/1/deal.II/lib/libbase.g.so: Subscriptor::~Subscriptor()
-#1 /u/bangerth/p/deal.II/1/deal.II/lib/libdeal_II_2d.g.so: FiniteElement<2>::~FiniteElement()
-#2 ./step-6: FE_Poly<TensorProductPolynomials<2>, 2>::~FE_Poly()
-#3 ./step-6: FE_Q<2>::~FE_Q()
-#4 ./step-6: LaplaceProblem<2>::~LaplaceProblem()
+#1 /u/bangerth/p/deal.II/1/deal.II/lib/libdeal_II_2d.g.so: FiniteElement<2>::~FiniteElement()
+#2 ./step-6: FE_Poly<TensorProductPolynomials<2>, 2>::~FE_Poly()
+#3 ./step-6: FE_Q<2>::~FE_Q()
+#4 ./step-6: LaplaceProblem<2>::~LaplaceProblem()
#5 ./step-6: main
--------------------------------------------------------
make: *** [run] Aborted
<a name="Results"></a>
<h1>Results</h1>
-<p>
+
The program generates two kinds of output. The first are the output
files <code>solution-adaptive-q1.gmv</code>,
<code>solution-global-q1.gmv</code>, and
<code>solution-global-q2.gmv</code>. We show the latter in a 3d view
here:
-</p>
-<p align="center">
- <img src="step-7.data/solution-3d.jpg" alt="solution" width="75%">
-</p>
+
+@image html step-7.solution.png
+
+
-<p>
Secondly, the program writes tables not only to disk, but also to the
screen while running:
-</p>
-<pre>
-<code>
+
+@code
examples/step-7> make run
============================ Running step-7
Solving with Q1 elements, adaptive refinement
4 1024 1.571e-02 1.99 1.265e-04 7.88 2.98
5 4096 3.937e-03 2.00 1.587e-05 7.97 2.99
6 16384 9.847e-04 2.00 1.986e-06 7.99 3.00
-</code>
-</pre>
+@endcode
-<p>
One can see the error reduction upon grid refinement, and for the
cases where global refinement was performed, also the convergence
rates can be seen. The linear and quadratic convergence rates of Q1
-and Q2 elements in the H<sup>1</sup> norm can clearly be seen, as
-are the quadratic and cubic rates in the L<sub>2</sub> norm.
-</p>
+and Q2 elements in the $H^1$ norm can clearly be seen, as
+are the quadratic and cubic rates in the $L_2$ norm.
+
+
-<p>
Finally, the program generated various LaTeX tables. We show here
the convergence table of the Q2 element with global refinement, after
converting the format to HTML:
-</p>
-<DIV ALIGN="CENTER">
+
<TABLE CELLPADDING=3 BORDER="1" ALIGN="CENTER">
-<TR><TD ALIGN="CENTER" COLSPAN=2><P>
+<TR><TD ALIGN="CENTER" COLSPAN=2>
n cells</TD>
<TD ALIGN="CENTER" COLSPAN=2><I>H</I><SUP>1</SUP>-error</TD>
<TD ALIGN="CENTER" COLSPAN=3><I>L</I><SUP>2</SUP>-error</TD>
<TD ALIGN="CENTER">7.99</TD>
<TD ALIGN="CENTER">3.00</TD>
</TR>
-</TABLE></DIV>
-<BR>
+</TABLE>
+