]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add step-26
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 4 Sep 2006 02:39:18 +0000 (02:39 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 4 Sep 2006 02:39:18 +0000 (02:39 +0000)
git-svn-id: https://svn.dealii.org/trunk@13829 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-26/Makefile [new file with mode: 0644]
deal.II/examples/step-26/step-26.cc [new file with mode: 0644]

diff --git a/deal.II/examples/step-26/Makefile b/deal.II/examples/step-26/Makefile
new file mode 100644 (file)
index 0000000..7f74cb9
--- /dev/null
@@ -0,0 +1,154 @@
+# $Id: Makefile 11779 2005-11-23 15:54:25Z wolf $
+
+
+# For the small projects Makefile, you basically need to fill in only
+# four fields.
+#
+# The first is the name of the application. It is assumed that the
+# application name is the same as the base file name of the single C++
+# file from which the application is generated.
+target = $(basename $(shell echo step-*.cc))
+
+# The second field determines whether you want to run your program in
+# debug or optimized mode. The latter is significantly faster, but no
+# run-time checking of parameters and internal states is performed, so
+# you should set this value to `on' while you develop your program,
+# and to `off' when running production computations.
+debug-mode = off
+
+
+# As third field, we need to give the path to the top-level deal.II
+# directory. You need to adjust this to your needs. Since this path is
+# probably the most often needed one in the Makefile internals, it is
+# designated by a single-character variable, since that can be
+# reference using $D only, i.e. without the parentheses that are
+# required for most other parameters, as e.g. in $(target).
+D = ../../
+
+
+# The last field specifies the names of data and other files that
+# shall be deleted when calling `make clean'. Object and backup files,
+# executables and the like are removed anyway. Here, we give a list of
+# files in the various output formats that deal.II supports.
+clean-up-files = *gmv *gnuplot *gpl *eps *pov
+
+
+
+
+#
+#
+# Usually, you will not need to change something beyond this point.
+#
+#
+# The next statement tell the `make' program where to find the
+# deal.II top level directory and to include the file with the global
+# settings
+include $D/common/Make.global_options
+
+
+# Since the whole project consists of only one file, we need not
+# consider difficult dependencies. We only have to declare the
+# libraries which we want to link to the object file, and there need
+# to be two sets of libraries: one for the debug mode version of the
+# application and one for the optimized mode. Here we have selected
+# the versions for 2d. Note that the order in which the libraries are
+# given here is important and that your applications won't link
+# properly if they are given in another order.
+#
+# You may need to augment the lists of libraries when compiling your
+# program for other dimensions, or when using third party libraries
+libs.g   = $(lib-deal2-3d.g) \
+          $(lib-lac.g)      \
+           $(lib-base.g)
+libs.o   = $(lib-deal2-3d.o) \
+          $(lib-lac.o)      \
+           $(lib-base.o)
+
+
+# We now use the variable defined above which switch between debug and
+# optimized mode to select the set of libraries to link with. Included
+# in the list of libraries is the name of the object file which we
+# will produce from the single C++ file. Note that by default we use
+# the extension .g.o for object files compiled in debug mode and .o for
+# object files in optimized mode (or whatever the local default on your
+# system is instead of .o).
+ifeq ($(debug-mode),on)
+  libraries = $(target).g.$(OBJEXT) $(libs.g)
+else
+  libraries = $(target).$(OBJEXT) $(libs.o)
+endif
+
+
+# Now comes the first production rule: how to link the single object
+# file produced from the single C++ file into the executable. Since
+# this is the first rule in the Makefile, it is the one `make' selects
+# if you call it without arguments.
+$(target) : $(libraries)
+       @echo ============================ Linking $@
+       @$(CXX) -o $@$(EXEEXT) $^ $(LIBS) $(LDFLAGS)
+
+
+# To make running the application somewhat independent of the actual
+# program name, we usually declare a rule `run' which simply runs the
+# program. You can then run it by typing `make run'. This is also
+# useful if you want to call the executable with arguments which do
+# not change frequently. You may then want to add them to the
+# following rule:
+run: $(target)
+       @echo ============================ Running $<
+       @./$(target)$(EXEEXT)
+
+
+# As a last rule to the `make' program, we define what to do when
+# cleaning up a directory. This usually involves deleting object files
+# and other automatically created files such as the executable itself,
+# backup files, and data files. Since the latter are not usually quite
+# diverse, you needed to declare them at the top of this file.
+clean:
+       -rm -f *.$(OBJEXT) *~ Makefile.dep $(target)$(EXEEXT) $(clean-up-files)
+
+
+# Since we have not yet stated how to make an object file from a C++
+# file, we should do so now. Since the many flags passed to the
+# compiler are usually not of much interest, we suppress the actual
+# command line using the `at' sign in the first column of the rules
+# and write the string indicating what we do instead.
+./%.g.$(OBJEXT) :
+       @echo ==============debug========= $(<F)
+       @$(CXX) $(CXXFLAGS.g) -c $< -o $@
+./%.$(OBJEXT) :
+       @echo ==============optimized===== $(<F)
+       @$(CXX) $(CXXFLAGS.o) -c $< -o $@
+
+
+# The following statement tells make that the rules `run' and `clean'
+# are not expected to produce files of the same name as Makefile rules
+# usually do.
+.PHONY: run clean
+
+
+# Finally there is a rule which you normally need not care much about:
+# since the executable depends on some include files from the library,
+# besides the C++ application file of course, it is necessary to
+# re-generate the executable when one of the files it depends on has
+# changed. The following rule to created a dependency file
+# `Makefile.dep', which `make' uses to determine when to regenerate
+# the executable. This file is automagically remade whenever needed,
+# i.e. whenever one of the cc-/h-files changed. Make detects whether
+# to remake this file upon inclusion at the bottom of this file.
+#
+# If the creation of Makefile.dep fails, blow it away and fail
+Makefile.dep: $(target).cc Makefile \
+              $(shell echo $D/*/include/*/*.h)
+       @echo ============================ Remaking $@
+       @$D/common/scripts/make_dependencies  $(INCLUDE) -B. $(target).cc \
+               > $@ \
+         || (rm -f $@ ; false)
+       @if test -s $@ ; then : else rm $@ ; fi
+
+
+# To make the dependencies known to `make', we finally have to include
+# them:
+include Makefile.dep
+
+
diff --git a/deal.II/examples/step-26/step-26.cc b/deal.II/examples/step-26/step-26.cc
new file mode 100644 (file)
index 0000000..0cffe51
--- /dev/null
@@ -0,0 +1,759 @@
+/* $Id: step-4.cc 13345 2006-07-06 13:44:43Z guido $ */
+/* Author: Wolfgang Bangerth, University of Heidelberg, 1999 */
+
+/*    $Id: step-4.cc 13345 2006-07-06 13:44:43Z guido $       */
+/*    Version: $Name$                                          */
+/*                                                                */
+/*    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors */
+/*                                                                */
+/*    This file is subject to QPL and may not be  distributed     */
+/*    without copyright and license information. Please refer     */
+/*    to the file deal.II/doc/license.html for the  text  and     */
+/*    further information on this license.                        */
+
+                                 // @sect3{Include files}
+
+                                // The first few (many?) include
+                                // files have already been used in
+                                // the previous example, so we will
+                                // not explain their meaning here
+                                // again.
+#include <grid/tria.h>
+#include <dofs/dof_handler.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_accessor.h>
+#include <fe/fe_q.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_values.h>
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <lac/sparse_matrix.h>
+#include <lac/solver_cg.h>
+#include <lac/precondition.h>
+
+#include <numerics/data_out.h>
+#include <fstream>
+#include <iostream>
+
+                                // This is new, however: in the previous
+                                // example we got some unwanted output from
+                                // the linear solvers. If we want to suppress
+                                // it, we have to include this file and add a
+                                // single line somewhere to the program (see
+                                // the main() function below for that):
+#include <base/logstream.h>
+
+
+#include <algorithm>
+#include <numeric>
+#include <grid/tria_boundary.h>
+
+
+class PointDefinedSurface : public StraightBoundary<3>
+{
+  public:
+    PointDefinedSurface (const std::string &filename);
+
+    Point<3> closest_point (const Point<3> &p) const;
+    
+                                    /**
+                                     * Let the new point be the
+                                     * arithmetic mean of the two
+                                     * vertices of the line.
+                                     *
+                                     * Refer to the general
+                                     * documentation of this class
+                                     * and the documentation of the
+                                     * base class for more
+                                     * information.
+                                     */
+    virtual Point<3>
+    get_new_point_on_line (const Triangulation<3>::line_iterator &line) const;
+
+                                    /**
+                                     * Let the new point be the
+                                     * arithmetic mean of the four
+                                     * vertices of this quad and the
+                                     * four midpoints of the lines,
+                                     * which are already created at
+                                     * the time of calling this
+                                     * function.
+                                     *
+                                     * Refer to the general
+                                     * documentation of this class
+                                     * and the documentation of the
+                                     * base class for more
+                                     * information.
+                                     */
+    virtual Point<3>
+    get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const;
+
+                                    /**
+                                     * Gives <tt>n=points.size()</tt>
+                                     * points that splits the
+                                     * p{StraightBoundary} line into
+                                     * p{n+1} partitions of equal
+                                     * lengths.
+                                     *
+                                     * Refer to the general
+                                     * documentation of this class
+                                     * and the documentation of the
+                                     * base class.
+                                     */
+    virtual void
+    get_intermediate_points_on_line (const Triangulation<3>::line_iterator &line,
+                                    std::vector<Point<3> > &points) const;
+
+                                    /**
+                                     * Gives <tt>n=points.size()=m*m</tt>
+                                     * points that splits the
+                                     * p{StraightBoundary} quad into
+                                     * <tt>(m+1)(m+1)</tt> subquads of equal
+                                     * size.
+                                     *
+                                     * Refer to the general
+                                     * documentation of this class
+                                     * and the documentation of the
+                                     * base class.
+                                     */
+    virtual void
+    get_intermediate_points_on_quad (const Triangulation<3>::quad_iterator &quad,
+                                    std::vector<Point<3> > &points) const;
+  private:
+    std::vector<Point<3> > point_list;
+};
+
+
+PointDefinedSurface::PointDefinedSurface (const std::string &filename)
+{
+                                  // first read in all the points
+  {
+    std::ifstream in (filename.c_str());
+    AssertThrow (in, ExcIO());
+    
+    while (in)
+      {
+       Point<3> p;
+       in >> p;
+       point_list.push_back (p);
+      }
+
+    AssertThrow (point_list.size() > 1, ExcIO());
+  }
+  
+                                  // next fit a linear model through the data
+                                  // cloud to rectify it in a local
+                                  // coordinate system
+                                  //
+                                  // the first step is to move the center of
+                                  // mass of the points to the origin
+  {
+    const Point<3> c_o_m = std::accumulate (point_list.begin(),
+                                           point_list.end(),
+                                           Point<3>()) /
+                            point_list.size();
+    for (unsigned int i=0; i<point_list.size(); ++i)
+      point_list[i] -= c_o_m;
+  }
+
+                                  // next do a least squares fit to the
+                                  // function ax+by. this leads to the
+                                  // following equations:
+  
+                                  // min f(a,b) = sum_i (zi-a xi - b yi)^2 / 2
+                                  //
+                                  // f_a = sum_i (zi - a xi - b yi) xi = 0
+                                  // f_b = sum_i (zi - a xi - b yi) yi = 0
+                                  //
+                                  // f_a = (sum_i zi xi) - (sum xi^2) a - (sum xi yi) b = 0
+                                  // f_a = (sum_i zi yi) - (sum xi yi) a - (sum yi^2) b = 0
+  {
+    double A[2][2] = {{0,0},{0,0}};
+    double B[2] = {0,0};
+
+    for (unsigned int i=0; i<point_list.size(); ++i)
+      {
+       A[0][0] += point_list[i][0] * point_list[i][0];
+       A[0][1] += point_list[i][0] * point_list[i][1];
+       A[1][1] += point_list[i][1] * point_list[i][1];
+
+       B[0] += point_list[i][0] * point_list[i][2];
+       B[1] += point_list[i][1] * point_list[i][2];
+      }
+
+    const double det = A[0][0]*A[1][1]-2*A[0][1];
+    const double a = (A[1][1] * B[0] - A[0][1] * B[1]) / det;
+    const double b = (A[0][0] * B[1] - A[0][1] * B[0]) / det;
+
+
+                                    // with this information, we can rotate
+                                    // the points so that the corresponding
+                                    // least-squares fit would be the x-y
+                                    // plane
+    const Point<2> gradient_direction
+      = Point<2>(a,b) / std::sqrt(a*a+b*b);
+    const Point<2> orthogonal_direction
+      = Point<2>(-b,a) / std::sqrt(a*a+b*b);
+
+    const double stretch_factor = std::sqrt(1.+a*a+b*b);
+    
+    for (unsigned int i=0; i<point_list.size(); ++i)
+      {
+                                        // we can do that by, for each point,
+                                        // first subtract the points in the
+                                        // plane:
+       point_list[i][2] -= a*point_list[i][0] + b*point_list[i][1];
+
+                                      // we made a mistake here, though:
+                                      // we've shrunk the plan in the
+                                      // direction parallel to the
+                                      // gradient. we will have to correct
+                                      // for this:
+       const Point<2> xy (point_list[i][0],
+                          point_list[i][1]);
+       const double grad_distance = xy * gradient_direction;
+       const double orth_distance = xy * orthogonal_direction;
+
+                                        // we then have to stretch the points
+                                        // in the gradient direction. the
+                                        // stretch factor is defined above
+                                        // (zero if the original plane was
+                                        // already the xy plane, infinity if
+                                        // it was vertical)
+       const Point<2> new_xy
+         = (grad_distance * stretch_factor * gradient_direction +
+            orth_distance * orthogonal_direction);
+       point_list[i][0] = new_xy[0];
+       point_list[i][1] = new_xy[1];   
+      }
+  }
+}
+
+
+Point<3>
+PointDefinedSurface::closest_point (const Point<3> &p) const
+{
+  double distance = p.distance (point_list[0]);
+  Point<3> point = point_list[0];
+  
+  for (std::vector<Point<3> >::const_iterator i=point_list.begin();
+       i != point_list.end(); ++i)
+    {
+      const double d = p.distance (*i);
+      if (d < distance)
+       {
+         distance = d;
+         point = *i;
+       }
+    }
+
+  return point;
+}
+
+  
+Point<3>
+PointDefinedSurface::
+get_new_point_on_line (const Triangulation<3>::line_iterator &line) const
+{
+  return closest_point (StraightBoundary<3>::get_new_point_on_line (line));
+}
+
+
+
+Point<3>
+PointDefinedSurface::
+get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const
+{
+  return closest_point (StraightBoundary<3>::get_new_point_on_quad (quad));
+}
+
+
+
+void
+PointDefinedSurface::
+get_intermediate_points_on_line (const Triangulation<3>::line_iterator &line,
+                                std::vector<Point<3> > &points) const
+{
+  StraightBoundary<3>::get_intermediate_points_on_line (line,
+                                                       points);
+  for (unsigned int i=0; i<points.size(); ++i)
+    points[i] = closest_point(points[i]);
+}
+
+
+
+void
+PointDefinedSurface::
+get_intermediate_points_on_quad (const Triangulation<3>::quad_iterator &quad,
+                                std::vector<Point<3> > &points) const
+{
+  StraightBoundary<3>::get_intermediate_points_on_quad (quad,
+                                                       points);
+  for (unsigned int i=0; i<points.size(); ++i)
+    points[i] = closest_point(points[i]);
+}
+
+
+
+PointDefinedSurface pds("surface-points");
+
+
+
+
+
+
+
+
+                                 // @sect3{The <code>LaplaceProblem</code> class template}
+
+                                // This is again the same
+                                // <code>LaplaceProblem</code> class as in the
+                                // previous example. The only
+                                // difference is that we have now
+                                // declared it as a class with a
+                                // template parameter, and the
+                                // template parameter is of course
+                                // the spatial dimension in which we
+                                // would like to solve the Laplace
+                                // equation. Of course, several of
+                                // the member variables depend on
+                                // this dimension as well, in
+                                // particular the Triangulation
+                                // class, which has to represent
+                                // quadrilaterals or hexahedra,
+                                // respectively. Apart from this,
+                                // everything is as before.
+template <int dim>
+class LaplaceProblem 
+{
+  public:
+    LaplaceProblem ();
+    void run ();
+    
+  private:
+    void make_grid_and_dofs ();
+    void assemble_system ();
+    void solve ();
+    void output_results () const;
+
+    Triangulation<dim>   triangulation;
+    FE_Q<dim>            fe;
+    DoFHandler<dim>      dof_handler;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+
+    Vector<double>       solution;
+    Vector<double>       system_rhs;
+};
+
+
+                                 // @sect3{Right hand side and boundary values}
+
+
+
+
+template <int dim>
+class BoundaryValues : public Function<dim> 
+{
+  public:
+    BoundaryValues () : Function<dim>() {};
+    
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component = 0) const;
+};
+
+
+
+template <int dim>
+double BoundaryValues<dim>::value (const Point<dim> &p,
+                                  const unsigned int /*component*/) const 
+{
+  return std::max(p[dim-1], -5.);
+}
+
+
+
+                                 // @sect3{Implementation of the <code>LaplaceProblem</code> class}
+
+                                 // Next for the implementation of the class
+                                 // template that makes use of the functions
+                                 // above. As before, we will write everything
+                                 // as templates that have a formal parameter
+                                 // <code>dim</code> that we assume unknown at the time
+                                 // we define the template functions. Only
+                                 // later, the compiler will find a
+                                 // declaration of <code>LaplaceProblem@<2@></code> (in
+                                 // the <code>main</code> function, actually) and
+                                 // compile the entire class with <code>dim</code>
+                                 // replaced by 2, a process referred to as
+                                 // `instantiation of a template'. When doing
+                                 // so, it will also replace instances of
+                                 // <code>RightHandSide@<dim@></code> by
+                                 // <code>RightHandSide@<2@></code> and instantiate the
+                                 // latter class from the class template.
+                                 //
+                                 // In fact, the compiler will also find a
+                                 // declaration <code>LaplaceProblem@<3@></code> in
+                                 // <code>main()</code>. This will cause it to again go
+                                 // back to the general
+                                 // <code>LaplaceProblem@<dim@></code> template, replace
+                                 // all occurrences of <code>dim</code>, this time by
+                                 // 3, and compile the class a second
+                                 // time. Note that the two instantiations
+                                 // <code>LaplaceProblem@<2@></code> and
+                                 // <code>LaplaceProblem@<3@></code> are completely
+                                 // independent classes; their only common
+                                 // feature is that they are both instantiated
+                                 // from the same general template, but they
+                                 // are not convertible into each other, for
+                                 // example, and share no code (both
+                                 // instantiations are compiled completely
+                                 // independently).
+
+
+                                 // @sect4{LaplaceProblem::LaplaceProblem}
+
+                                // After this introduction, here is the
+                                // constructor of the <code>LaplaceProblem</code>
+                                // class. It specifies the desired polynomial
+                                // degree of the finite elements and
+                                // associates the DoFHandler to the
+                                // triangulation just as in the previous
+                                // example program, step-3:
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem () :
+                fe (1),
+               dof_handler (triangulation)
+{}
+
+
+                                 // @sect4{LaplaceProblem::make_grid_and_dofs}
+
+                                // Grid creation is something
+                                // inherently dimension
+                                // dependent. However, as long as the
+                                // domains are sufficiently similar
+                                // in 2D or 3D, the library can
+                                // abstract for you. In our case, we
+                                // would like to again solve on the
+                                // square [-1,1]x[-1,1] in 2D, or on
+                                // the cube [-1,1]x[-1,1]x[-1,1] in
+                                // 3D; both can be termed
+                                // <code>hyper_cube</code>, so we may use the
+                                // same function in whatever
+                                // dimension we are. Of course, the
+                                // functions that create a hypercube
+                                // in two and three dimensions are
+                                // very much different, but that is
+                                // something you need not care
+                                // about. Let the library handle the
+                                // difficult things.
+                                //
+                                // Likewise, associating a degree of freedom
+                                // with each vertex is something which
+                                // certainly looks different in 2D and 3D,
+                                // but that does not need to bother you
+                                // either. This function therefore looks
+                                // exactly like in the previous example,
+                                // although it performs actions that in their
+                                // details are quite different if <code>dim</code>
+                                // happens to be 3. The only significant
+                                // difference from a user's perspective is
+                                // the number of cells resulting, which is
+                                // much higher in three than in two space
+                                // dimensions!
+template <int dim>
+void LaplaceProblem<dim>::make_grid_and_dofs ()
+{
+  GridGenerator::hyper_cube (triangulation, -30, 30);
+
+  for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+    if (triangulation.begin()->face(f)->center()[2] > 15)
+      {
+       triangulation.begin()->face(f)->set_boundary_indicator (1);
+       for (unsigned int i=0; i<GeometryInfo<dim>::lines_per_face; ++i)
+         triangulation.begin()->face(f)->line(i)->set_boundary_indicator (1);
+       break;
+      }
+  triangulation.set_boundary (1, pds);
+  
+  
+  for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+    if (triangulation.begin()->vertex(v)[2] > 0)
+      triangulation.begin()->vertex(v)
+       = pds.closest_point (Point<3>(triangulation.begin()->vertex(v)[0],
+                                     triangulation.begin()->vertex(v)[1],
+                                     0));
+       
+  for (unsigned int i=0; i<7; ++i)
+    {
+      for (typename Triangulation<dim>::active_cell_iterator
+            cell = triangulation.begin_active();
+          cell != triangulation.end(); ++cell)
+       for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+         if (cell->face(f)->boundary_indicator() == 1)
+           cell->set_refine_flag ();
+      
+      triangulation.execute_coarsening_and_refinement ();
+
+      std::cout << "Refinement cycle " << i << std::endl
+               << "   Number of active cells: "
+               << triangulation.n_active_cells()
+               << std::endl
+               << "   Total number of cells: "
+               << triangulation.n_cells()
+               << std::endl;
+
+    }
+  
+  
+  dof_handler.distribute_dofs (fe);
+
+  std::cout << "   Number of degrees of freedom: "
+           << dof_handler.n_dofs()
+           << std::endl;
+
+  sparsity_pattern.reinit (dof_handler.n_dofs(),
+                          dof_handler.n_dofs(),
+                          dof_handler.max_couplings_between_dofs());
+  DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+  sparsity_pattern.compress();
+
+  system_matrix.reinit (sparsity_pattern);
+
+  solution.reinit (dof_handler.n_dofs());
+  system_rhs.reinit (dof_handler.n_dofs());
+}
+
+
+                                 // @sect4{LaplaceProblem::assemble_system}
+
+                                // Unlike in the previous example, we
+                                // would now like to use a
+                                // non-constant right hand side
+                                // function and non-zero boundary
+                                // values. Both are tasks that are
+                                // readily achieved with a only a few
+                                // new lines of code in the
+                                // assemblage of the matrix and right
+                                // hand side.
+                                //
+                                // More interesting, though, is the
+                                // way we assemble matrix and right
+                                // hand side vector dimension
+                                // independently: there is simply no
+                                // difference to the 
+                                // two-dimensional case. Since the
+                                // important objects used in this
+                                // function (quadrature formula,
+                                // FEValues) depend on the dimension
+                                // by way of a template parameter as
+                                // well, they can take care of
+                                // setting up properly everything for
+                                // the dimension for which this
+                                // function is compiled. By declaring
+                                // all classes which might depend on
+                                // the dimension using a template
+                                // parameter, the library can make
+                                // nearly all work for you and you
+                                // don't have to care about most
+                                // things.
+template <int dim>
+void LaplaceProblem<dim>::assemble_system () 
+{  
+  MatrixTools::create_laplace_matrix (dof_handler,
+                                     QGauss<dim>(2),
+                                     system_matrix);
+  system_rhs = 0;
+  
+  std::map<unsigned int,double> boundary_values;
+  VectorTools::interpolate_boundary_values (dof_handler,
+                                           0,
+                                           BoundaryValues<dim>(),
+                                           boundary_values);
+  MatrixTools::apply_boundary_values (boundary_values,
+                                     system_matrix,
+                                     solution,
+                                     system_rhs);
+}
+
+
+                                 // @sect4{LaplaceProblem::solve}
+
+                                // Solving the linear system of
+                                // equations is something that looks
+                                // almost identical in most
+                                // programs. In particular, it is
+                                // dimension independent, so this
+                                // function is copied verbatim from the
+                                // previous example.
+template <int dim>
+void LaplaceProblem<dim>::solve () 
+{
+  SolverControl           solver_control (1000, 1e-12);
+  SolverCG<>              cg (solver_control);
+
+  PreconditionSSOR<> preconditioner;
+  preconditioner.initialize(system_matrix, 1.2);
+
+  cg.solve (system_matrix, solution, system_rhs,
+           preconditioner);
+}
+
+
+                                 // @sect4{LaplaceProblem::output_results}
+
+                                // This function also does what the
+                                // respective one did in step-3. No changes
+                                // here for dimension independence either.
+                                 //
+                                 // The only difference to the previous
+                                 // example is that we want to write output in
+                                 // GMV format, rather than for gnuplot (GMV
+                                 // is another graphics program that, contrary
+                                 // to gnuplot, shows data in nice colors,
+                                 // allows rotation of geometries with the
+                                 // mouse, and generates reasonable
+                                 // representations of 3d data; for ways to
+                                 // obtain it see the ReadMe file of
+                                 // deal.II). To write data in this format, we
+                                 // simply replace the
+                                 // <code>data_out.write_gnuplot</code> call by
+                                 // <code>data_out.write_gmv</code>.
+                                 //
+                                 // Since the program will run both 2d and 3d
+                                 // versions of the laplace solver, we use the
+                                 // dimension in the filename to generate
+                                 // distinct filenames for each run (in a
+                                 // better program, one would check whether
+                                 // `dim' can have other values than 2 or 3,
+                                 // but we neglect this here for the sake of
+                                 // brevity).
+template <int dim>
+void LaplaceProblem<dim>::output_results () const
+{
+  DataOut<dim> data_out;
+
+  data_out.attach_dof_handler (dof_handler);
+  data_out.add_data_vector (solution, "solution");
+
+  data_out.build_patches ();
+
+  std::ofstream output (dim == 2 ?
+                       "solution-2d.gmv" :
+                       "solution-3d.gmv");
+  data_out.write_gmv (output);
+}
+
+
+
+                                 // @sect4{LaplaceProblem::run}
+
+                                 // This is the function which has the
+                                // top-level control over
+                                // everything. Apart from one line of
+                                // additional output, it is the same
+                                // as for the previous example.
+template <int dim>
+void LaplaceProblem<dim>::run () 
+{
+  std::cout << "Solving problem in " << dim << " space dimensions." << std::endl;
+  
+  make_grid_and_dofs();
+  assemble_system ();
+  solve ();
+  output_results ();
+}
+
+
+                                 // @sect3{The <code>main</code> function}
+
+                                // And this is the main function. It also
+                                // looks mostly like in step-3, but if you
+                                // look at the code below, note how we first
+                                // create a variable of type
+                                // <code>LaplaceProblem@<2@></code> (forcing the
+                                // compiler to compile the class template
+                                // with <code>dim</code> replaced by <code>2</code>) and run a
+                                // 2d simulation, and then we do the whole
+                                // thing over in 3d.
+                                //
+                                // In practice, this is probably not what you
+                                // would do very frequently (you probably
+                                // either want to solve a 2d problem, or one
+                                // in 3d, but not both at the same
+                                // time). However, it demonstrates the
+                                // mechanism by which we can simply change
+                                // which dimension we want in a single place,
+                                // and thereby force the compiler to
+                                // recompile the dimension independent class
+                                // templates for the dimension we
+                                // request. The emphasis here lies on the
+                                // fact that we only need to change a single
+                                // place. This makes it rather trivial to
+                                // debug the program in 2d where computations
+                                // are fast, and then switch a single place
+                                // to a 3 to run the much more computing
+                                // intensive program in 3d for `real'
+                                // computations.
+                                //
+                                // Each of the two blocks is enclosed in
+                                // braces to make sure that the
+                                // <code>laplace_problem_2d</code> variable goes out
+                                // of scope (and releases the memory it
+                                // holds) before we move on to allocate
+                                // memory for the 3d case. Without the
+                                // additional braces, the
+                                // <code>laplace_problem_2d</code> variable would only
+                                // be destroyed at the end of the function,
+                                // i.e. after running the 3d problem, and
+                                // would needlessly hog memory while the 3d
+                                // run could actually use it.
+                                 //
+                                 // Finally, the first line of the function is
+                                 // used to suppress some output.  Remember
+                                 // that in the previous example, we had the
+                                 // output from the linear solvers about the
+                                 // starting residual and the number of the
+                                 // iteration where convergence was
+                                 // detected. This can be suppressed through
+                                 // the <code>deallog.depth_console(0)</code> call.
+                                 //
+                                 // The rationale here is the following: the
+                                 // deallog (i.e. deal-log, not de-allog)
+                                 // variable represents a stream to which some
+                                 // parts of the library write output. It
+                                 // redirects this output to the console and
+                                 // if required to a file. The output is
+                                 // nested in a way so that each function can
+                                 // use a prefix string (separated by colons)
+                                 // for each line of output; if it calls
+                                 // another function, that may also use its
+                                 // prefix which is then printed after the one
+                                 // of the calling function. Since output from
+                                 // functions which are nested deep below is
+                                 // usually not as important as top-level
+                                 // output, you can give the deallog variable
+                                 // a maximal depth of nested output for
+                                 // output to console and file. The depth zero
+                                 // which we gave here means that no output is
+                                 // written. By changing it you can get more
+                                 // information about the innards of the
+                                 // library.
+int main () 
+{
+  deallog.depth_console (0);
+  {
+    LaplaceProblem<3> laplace_problem_3d;
+    laplace_problem_3d.run ();
+  }
+  
+  return 0;
+}

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.