]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Almost finish everything.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 16 Apr 2002 08:34:13 +0000 (08:34 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 16 Apr 2002 08:34:13 +0000 (08:34 +0000)
git-svn-id: https://svn.dealii.org/trunk@5655 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/tutorial/chapter-2.step-by-step/step-13.data/intro.html
deal.II/doc/tutorial/chapter-2.step-by-step/step-13.data/results.html

index c59cfdf6092aac8d5c180aab3ea01f0da3d3390c..0aba96ec67b6b8ed1cb4338daedef857385eb93b 100644 (file)
@@ -1,2 +1,177 @@
 <a name="Intro"></a>
 <h1>Introduction</h1>
+
+<h3>Background and purpose</h3>
+
+<p>
+In this example program, we will not so much be concerned with
+describing new ways how to use deal.II and its facilities, but rather
+with presenting methods of writing modular and extensible finite
+element programs. The main reason for this is the size and complexity
+of modern research software: applications implementing modern error
+estimation concepts and adaptive solution methods tend to become
+rather large. For example, the three largest applications by the main
+authors of deal.II, are at the time of writing of this example
+program:
+<ol>
+<li> a program for solving conservation hyperbolic equations by the
+     Discontinuous Galerkin Finite Element method: 33,775 lines of
+     code; 
+<li> a parameter estimation program: 28,980 lines of code;
+<li> a wave equation solver: 21,020 lines of code.
+</ol>
+(The library proper - without example programs and
+test suite - has slightly more than 150,000 lines of code.) In
+the opinion of the author of this example program, the sizes of these
+applications are at the edge of what one person, even an experienced
+programmer, can manage. 
+</p>
+
+<p>
+The numbers above make one thing rather clear: monolithic programs that
+are not broken up into smaller, mostly independent pieces have no way
+of surviving, since even the author will quickly lose the overview of
+the various dependencies between different parts of a program. Only
+data encapsulation, for example using object oriented programming
+methods, and modularization by defining small but fixed interfaces can
+help structure data flow and mutual interdependencies. It is also an
+absolute prerequisite if more than one person is developing a program,
+since otherwise confusion will quickly prevail as one developer
+would need to know if another changed something about the internals of
+a different module if they were not cleanly separated.
+</p>
+
+<p>
+In previous examples, you have seen how the library itself is broken
+up into several complexes each building atop the underying ones, but
+relatively independent of the other ones:
+<ol>
+<li>the triangulation class complex, with associated iterator classes;
+<li>the finite element classes;
+<li>the DoFHandler class complex, with associated iterators, built on
+    the triangulation and finite element classes;
+<li>the classes implementing mappings between unit and real cells;
+<li>the FEValues class complex, built atop the finite elements and
+    mappings.
+</ol>
+Besides these, and a large number of smaller classes, there are of
+course the following ``tool'' modules:
+<ol>
+<li>output in various graphical formats;
+<li>linear algebra classes.
+</ol>
+</p>
+
+
+<p>
+The goal of this program is now to give an example of how a relatively
+simple finite element program could be structured such that we end up
+with a set of modules that are as independent of each other as
+possible. This allows to change the program at one end, without having to
+worry that it might break at the other, as long as we do not touch the
+interface through which the two ends communicate. The interface in
+C++, of course, is the declaration of abstract base classes.
+</p>
+
+<p>
+Here, we will implement (again) a Laplace solver, although with a
+number of differences compared to previous example programs:
+<ol>
+<li>The classes that implement the process of numerically solving the
+    equation are no more responsible for driving the process of
+    ``solving-estimating error-refining-solving again'', but we delegate
+    this to external functions. This allows first to use it as a
+    building block in a larger context, where the solution of a
+    Laplace equation might only be one part (for example, in a
+    nonlinear problem, where Laplace equations might have to be solved
+    in each nonlinear step). It would also allow to build a framework
+    around this class that would allow using solvers for other
+    equations (but with the same external interface) instead, in case
+    some techniques shall be evaluated for different types of partial
+    differential equations.
+<li>It splits the process of evaluating the computed solution to a
+    separate set of classes. The reason is that one is usually not
+    interested in the solution of a PDE per se, but rather in certain
+    aspects of it. For example, one might wish to compute the traction
+    at a certain boundary in elastic computations, or in the signal of
+    a seismic wave at a receiver position at a given
+    location. Sometimes, one might have an interest in several of
+    these aspects. Since the evaluation of a solution is something
+    that does not usually affect the process of solution, we split it
+    off into a separate module, to allow for the development of such
+    evaluation filters independently of the development of the solver
+    classes.
+<li>Separate the classes that implement mesh refinement from the
+    classes that compute the solution.
+<li>Separate the description of the test case with which we will
+    present the program, from the rest of the program.
+</ol>
+</p>
+
+<p>
+The things the program does are not new. In fact, this is more like a
+melange of previous programs, cannibalizing various parts and
+functions from earlier examples. It is the way they are arranged in
+this program that should be the focus of the reader.
+</p>
+
+<p>
+Once you have worked through the program, you will remark that it is
+already somewhat complex in its structure. Nevertheless, it
+only has about 850 lines of code, without comments. In real
+applications, there would of course be comments and class
+documentation, which would bring that to maybe 1200 lines. Yet, compared to
+the applications listed above, this is still small, as they are 20 to
+25 times as large. For programs as large, a proper design right from
+the start is thus indispensible. Otherwise, it will have to be
+redesigned at one point in its life, once it becomes too large to be
+manageable.
+</p>
+
+<p>
+Despite of this, all three programs listed above have undergone major
+revisions, or even rewrites. The wave program, for example, was once
+entirely teared to parts when it was still significantly smaller, just
+to assemble it again in a more modular form. By that time, it had
+become impossible to add functionality without affecting older parts
+of the code (the main problem with the code was the data flow: in time
+dependent application, the major concern is when to store data to disk
+and when to reload it again; if this is not done in an organized
+fashion, then you end up with data released too early, loaded too
+late, or not released at all). Although the present example program
+thus draws from sevelar years of experience, it is certainly not
+without flaws in its design, and in particular might not be suited for
+an application where the objective is different. It should serve as an
+inspiration for writing your own application in a modular way, to
+avoid the pitfalls of too closely coupled codes.
+</p>
+
+
+<h3>What the program does</h3>
+
+<p>
+What the program actually does is not even the main point of this
+program, the structure of the program is more important. However, in a
+few words, a description would be: solve the Laplace equation for a
+given right hand side such that the solution is the function 
+  <em> u(x,t)=exp(x+sin(10y+5x<sup>2</sup>)) </em>. The goal of the
+computation is to get the value of the solution at the point
+<em>x<sub>0</sub>=(0.5,0.5)</em>, and to compare the accuracy with
+which we resolve this value for two refinement criteria, namely global
+refinement and refinement by the error indicator by Kelly et al. which
+we have already used in previous examples.
+</p>
+
+<p>
+The results will, as usual, be discussed in the respective section of
+this document. In doing so, we will find a slightly irritating
+observation about the relative performance of the two refinement
+criteria. In a later example program, building atop this one, we will
+devise a different method that should hopefully perform better than
+the techniques discussed here.
+</p>
+
+<p>
+So much now for all the theoretical and anecdotal background. The best
+way of learning about a program is to look at it, so here it is:
+</p>
index 51041370d36a6f0346d2a5ed3349d39b9c6af197..81ccd6d393c6e284b3ef298674ffd68fa1a7d1c9 100644 (file)
@@ -2,38 +2,62 @@
 <h1>Results</h1>
 
 
+<p>
+The results of this program are not that interesting - after all
+its purpose was not to demonstrate some new mathematical idea, and
+also not how to program with deal.II, but rather to use the material
+which we have developed in the previous examples to form something
+which demonstrates a way to build modern finite element software in a
+modular and extensible way.
+</p>
+
+<p>
+Nevertheless, we of course show the results of the program. Of
+foremost interest is the point value computation, for which we had
+implemented the corresponding evaluation class. The results (i.e. the
+output) of the program looks as follows:
 <code>
 <pre>
-Running tests with "global" refinement criterion:
--------------------------------------------------
-Refinement cycle: 0 1 2 3 4 5 6
-DoFs  u(x_0)
-   25 1.2868
-   81 1.6945
-  289 1.4658
- 1089 1.5679
- 4225 1.5882
-16641 1.5932
-66049 1.5945
-
-Running tests with "kelly" refinement criterion:
-------------------------------------------------
-Refinement cycle: 0 1 2 3 4 5 6 7 8 9 10 11
-DoFs  u(x_0)
-   25 1.2868
-   47 0.8775
-   89 1.5365
-  165 1.2974
-  316 1.6442
-  589 1.5221
- 1090 1.5724
- 2035 1.5622
- 3754 1.5916
- 7100 1.5876
-13059 1.5942
-24749 1.5933
+    Running tests with "global" refinement criterion:
+    -------------------------------------------------
+    Refinement cycle: 0 1 2 3 4 5 6
+    DoFs  u(x_0)
+       25 1.2868
+       81 1.6945
+      289 1.4658
    1089 1.5679
    4225 1.5882
+    16641 1.5932
+    66049 1.5945
+
+    Running tests with "kelly" refinement criterion:
+    ------------------------------------------------
+    Refinement cycle: 0 1 2 3 4 5 6 7 8 9 10 11
+    DoFs  u(x_0)
+       25 1.2868
+       47 0.8775
+       89 1.5365
+      165 1.2974
+      316 1.6442
+      589 1.5221
    1090 1.5724
    2035 1.5622
    3754 1.5916
    7100 1.5876
+    13059 1.5942
+    24749 1.5933
 </pre>
 </code>
+</p>
+
+<p>
+What surprises here is that the the exact value is 1.59492, and that
+it is therefore suprisingly complicated to compute the solution to
+only one per cent accuracy, although the solution is smooth (in fact
+infinite often differentiable). This smoothness is shown in the
+graphical output generated by the program, here coarse grid and the
+first 9 refinement steps of the Kelly refinement indicator:
+</p>
 
 <table width="80%" align="center">
   <tr>
@@ -91,3 +115,94 @@ DoFs  u(x_0)
     </td>
   </tr>
 </table>
+
+<p>
+While we're already at watching pictures, this is the eighth grid, as
+viewed from top:
+</p>
+
+<p align="center">
+  <img src="step-13.data/pix-kelly/grid-kelly-8.gif" alt="Kelly, grid 8">
+</p>
+
+
+
+<p>
+However, we are not yet finished with evaluation the point value
+computation. In fact, plotting the error
+<em>e=|u(x<sub>h</sub>)-u<sub>h</sub>(x<sub>0</sub>)|</em> for the two
+refinement criteria yields the following picture:
+</p>
+
+<p align="center">
+  <img src="step-13.data/error.gif" alt="error" width="80%">
+</p>
+
+
+<p>
+What <em>is</em> disturbing about this picture is that not only is the
+adaptive mesh refinement not better than global refinement as one
+would usually expect, it is even significantly worse since its
+convergence is irregular, preventing all extrapolation techniques when
+using the values of subsequent meshes! On the other hand, global
+refinement provides a perfect <em>1/N</em> or <em>h<sup>-2</sup></em>
+convergence history and provides every opportunity to even improve on
+the point values by extrapolation. Global mesh refinement must
+therefore be considered superior in this example! This is even more
+surprising as the evaluation point is not somewhere in the left part
+where the mesh is coarse, but rather to the right and the adaptive
+refinement should refine the mesh around the evaluation point as well.
+</p>
+
+<p>
+We thus close the discussion of this example program with a question:
+</p>
+<p align="center">
+  <strong><em>What is wrong with adaptivity if it is not better than
+  global refinement?</em></strong>
+</p>
+
+
+<p>
+<em>Exercise at the end of this example:</em> There is a simple reason
+for the bad and irregular behavior of the adapted mesh solutions. It
+is simple to find out by looking at the mesh around the evaluation
+point in each of the steps - the data for this is in the output files
+of the program. An exercise would therefore be to modify the mesh
+refinement routine such that the problem (once you remark it) is
+avoided. The second exercise is to check whether the results are then
+better than global refinement, and if so if even a better order of
+convergence (in terms of the number of degrees of freedom) is
+achieved, or only by a better constant.
+</p>
+
+<p>
+(<em>Very brief answers for the impatient:</em> at steps with larger
+errors, the mesh is not regular at the point of evaluation, i.e. some
+of the adjacent cells have hanging nodes; this destroys some
+superapproximation effects of which the globally refined mesh can
+profit. Answer 2: this quick hack
+<code><pre>
+    bool refinement_indicated = false;
+    typename Triangulation<dim>::active_cell_iterator cell;
+    for (cell=triangulation->begin_active();
+        cell!=triangulation->end(); ++cell)
+      for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+       if (cell->vertex(v) == Point<dim>(.5,.5))
+         {
+           cell->clear_coarsen_flag();
+           refinement_indicated |= cell->refine_flag_set();
+         };
+    if (refinement_indicated)
+      for (cell=triangulation->begin_active();
+          cell!=triangulation->end(); ++cell)
+       for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+         if (cell->vertex(v) == Point<dim>(.5,.5))
+           cell->set_refine_flag ();
+</pre></code>
+in the refinement function of the Kelly refinement class right before
+executing refinement would improve the results (exercise: what does
+the code do?), making them consistently better than global
+refinement. Behavior is still irregular, though, so no results about
+an order of convergence are possible.)
+</p>

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.