--- /dev/null
+# $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
+
+
--- /dev/null
+/* $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;
+}