]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
undo r24637
authorheister <heister@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 19 Oct 2011 16:04:43 +0000 (16:04 +0000)
committerheister <heister@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 19 Oct 2011 16:04:43 +0000 (16:04 +0000)
git-svn-id: https://svn.dealii.org/trunk@24643 0785d39b-7218-0410-832d-ea1e28bc413d

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

diff --git a/deal.II/examples/step-41/step-41.cc b/deal.II/examples/step-41/step-41.cc
new file mode 100644 (file)
index 0000000..6ae88cb
--- /dev/null
@@ -0,0 +1,940 @@
+/* $Id: step-4.cc 24093 2011-08-16 13:58:12Z bangerth $ */
+/* Author: Wolfgang Bangerth, University of Heidelberg, 1999 */
+
+/*    $Id: step-4.cc 24093 2011-08-16 13:58:12Z bangerth $       */
+/*                                                                */
+/*    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 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 <deal.II/grid/tria.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+#include <deal.II/numerics/vectors.h>
+#include <deal.II/numerics/matrices.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_sparsity_pattern.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_bicgstab.h>
+#include <deal.II/lac/precondition.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <fstream>
+#include <iostream>
+#include <list>
+
+                                // 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 <deal.II/base/logstream.h>
+
+                                // The final step, as in previous
+                                // programs, is to import all the
+                                // deal.II class and function names
+                                // into the global namespace:
+using namespace dealii;
+
+                                 // @sect3{The <code>Step4</code> class template}
+
+                                // This is again the same
+                                // <code>Step4</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 Step4 
+{
+  public:
+    Step4 ();
+    void run ();
+    
+  private:
+    void make_grid ();
+    void setup_system();
+    void assemble_system ();
+    void projection_active_set ();
+    void solve ();
+    void output_results (Vector<double> vector_to_plot, const std::string& title) const;
+
+    Triangulation<dim>   triangulation;
+    FE_Q<dim>            fe;
+    DoFHandler<dim>      dof_handler;
+
+    ConstraintMatrix     constraints;
+    
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+    SparseMatrix<double> system_matrix_complete;
+
+    Vector<double>       solution;
+    Vector<double>       system_rhs;
+    Vector<double>       system_rhs_complete;
+    Vector<double>       resid_vector;
+    Vector<double>       active_set;
+
+    std::map<unsigned int,double> boundary_values;
+};
+
+
+                                 // @sect3{Right hand side and boundary values}
+
+                                // In the following, we declare two more
+                                // classes denoting the right hand side and
+                                // the non-homogeneous Dirichlet boundary
+                                // values. Both are functions of a
+                                // dim-dimensional space variable, so we
+                                // declare them as templates as well.
+                                //
+                                // Each of these classes is derived from a
+                                // common, abstract base class Function,
+                                // which declares the common interface which
+                                // all functions have to follow. In
+                                // particular, concrete classes have to
+                                // overload the <code>value</code> function,
+                                // which takes a point in dim-dimensional
+                                // space as parameters and shall return the
+                                // value at that point as a
+                                // <code>double</code> variable.
+                                //
+                                // The <code>value</code> function takes a
+                                // second argument, which we have here named
+                                // <code>component</code>: This is only meant
+                                // for vector valued functions, where you may
+                                // want to access a certain component of the
+                                // vector at the point
+                                // <code>p</code>. However, our functions are
+                                // scalar, so we need not worry about this
+                                // parameter and we will not use it in the
+                                // implementation of the functions. Inside
+                                // the library's header files, the Function
+                                // base class's declaration of the
+                                // <code>value</code> function has a default
+                                // value of zero for the component, so we
+                                // will access the <code>value</code>
+                                // function of the right hand side with only
+                                // one parameter, namely the point where we
+                                // want to evaluate the function. A value for
+                                // the component can then simply be omitted
+                                // for scalar functions.
+                                //
+                                // Note that the C++ language forces
+                                // us to declare and define a
+                                // constructor to the following
+                                // classes even though they are
+                                // empty. This is due to the fact
+                                // that the base class has no default
+                                // constructor (i.e. one without
+                                // arguments), even though it has a
+                                // constructor which has default
+                                // values for all arguments.
+template <int dim>
+class RightHandSide : public Function<dim> 
+{
+  public:
+    RightHandSide () : Function<dim>() {}
+    
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component = 0) const;
+};
+
+
+
+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>
+class Obstacle : public Function<dim> 
+{
+  public:
+    Obstacle () : Function<dim>() {}
+    
+    virtual double value (const Point<dim>   &p,
+                         const unsigned int  component = 0) const;
+};
+
+
+
+                                // For this example, we choose as right hand
+                                // side function to function $4(x^4+y^4)$ in
+                                // 2D, or $4(x^4+y^4+z^4)$ in 3D. We could
+                                // write this distinction using an
+                                // if-statement on the space dimension, but
+                                // here is a simple way that also allows us
+                                // to use the same function in 1D (or in 4D,
+                                // if you should desire to do so), by using a
+                                // short loop.  Fortunately, the compiler
+                                // knows the size of the loop at compile time
+                                // (remember that at the time when you define
+                                // the template, the compiler doesn't know
+                                // the value of <code>dim</code>, but when it later
+                                // encounters a statement or declaration
+                                // <code>RightHandSide@<2@></code>, it will take the
+                                // template, replace all occurrences of dim
+                                // by 2 and compile the resulting function);
+                                // in other words, at the time of compiling
+                                // this function, the number of times the
+                                // body will be executed is known, and the
+                                // compiler can optimize away the overhead
+                                // needed for the loop and the result will be
+                                // as fast as if we had used the formulas
+                                // above right away.
+                                //
+                                // The last thing to note is that a
+                                // <code>Point@<dim@></code> denotes a point in
+                                // dim-dimensionsal space, and its individual
+                                // components (i.e. $x$, $y$,
+                                // ... coordinates) can be accessed using the
+                                // () operator (in fact, the [] operator will
+                                // work just as well) with indices starting
+                                // at zero as usual in C and C++.
+template <int dim>
+double RightHandSide<dim>::value (const Point<dim> &p,
+                                 const unsigned int /*component*/) const 
+{
+  double return_value = 0;
+  // for (unsigned int i=0; i<dim; ++i)
+  //   return_value += 4*std::pow(p(i), 4);
+
+  return return_value;
+}
+
+
+                                // As boundary values, we choose x*x+y*y in
+                                // 2D, and x*x+y*y+z*z in 3D. This happens to
+                                // be equal to the square of the vector from
+                                // the origin to the point at which we would
+                                // like to evaluate the function,
+                                // irrespective of the dimension. So that is
+                                // what we return:
+template <int dim>
+double BoundaryValues<dim>::value (const Point<dim> &p,
+                                  const unsigned int /*component*/) const 
+{
+  double return_value = 0;
+
+  return return_value;
+}
+
+template <int dim>
+double Obstacle<dim>::value (const Point<dim> &p,
+                            const unsigned int /*component*/) const 
+{
+  return 2.0*p.square() - 0.5;
+}
+
+
+
+                                 // @sect3{Implementation of the <code>Step4</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>Step4@<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>Step4@<3@></code> in
+                                 // <code>main()</code>. This will cause it to
+                                 // again go back to the general
+                                 // <code>Step4@<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>Step4@<2@></code> and
+                                 // <code>Step4@<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{Step4::Step4}
+
+                                // After this introduction, here is the
+                                // constructor of the <code>Step4</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>
+Step4<dim>::Step4 ()
+               :
+                fe (1),
+               dof_handler (triangulation)
+{}
+
+
+                                 // @sect4{Step4::make_grid}
+
+                                // 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]\times [-1,1]$
+                                // in 2D, or on the cube $[-1,1] \times
+                                // [-1,1] \times [-1,1]$ in 3D; both can be
+                                // termed GridGenerator::hyper_cube(), 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.
+template <int dim>
+void Step4<dim>::make_grid ()
+{
+  GridGenerator::hyper_cube (triangulation, -1, 1);
+  triangulation.refine_global (6);
+  
+  std::cout << "   Number of active cells: "
+           << triangulation.n_active_cells()
+           << std::endl
+           << "   Total number of cells: "
+           << triangulation.n_cells()
+           << std::endl;
+}
+
+                                 // @sect4{Step4::setup_system}
+
+                                // This function 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 Step4<dim>::setup_system ()
+{
+  dof_handler.distribute_dofs (fe);
+
+  std::cout << "   Number of degrees of freedom: "
+           << dof_handler.n_dofs()
+           << std::endl;
+
+  CompressedSparsityPattern c_sparsity(dof_handler.n_dofs());
+  DoFTools::make_sparsity_pattern (dof_handler, c_sparsity, constraints, false);
+//   c_sparsity.compress ();
+  sparsity_pattern.copy_from(c_sparsity);
+  
+  system_matrix.reinit (sparsity_pattern);
+  system_matrix_complete.reinit (sparsity_pattern);
+  
+  solution.reinit (dof_handler.n_dofs());
+  system_rhs.reinit (dof_handler.n_dofs());
+  system_rhs_complete.reinit (dof_handler.n_dofs());
+  resid_vector.reinit (dof_handler.n_dofs());
+  active_set.reinit (dof_handler.n_dofs());
+}
+
+
+                                 // @sect4{Step4::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 Step4<dim>::assemble_system () 
+{  
+  QGauss<dim>  quadrature_formula(2);
+
+                                  // We wanted to have a non-constant right
+                                  // hand side, so we use an object of the
+                                  // class declared above to generate the
+                                  // necessary data. Since this right hand
+                                  // side object is only used locally in the
+                                  // present function, we declare it here as
+                                  // a local variable:
+  const RightHandSide<dim> right_hand_side;
+
+                                  // Compared to the previous example, in
+                                  // order to evaluate the non-constant right
+                                  // hand side function we now also need the
+                                  // quadrature points on the cell we are
+                                  // presently on (previously, we only
+                                  // required values and gradients of the
+                                  // shape function from the
+                                  // FEValues object, as well as
+                                  // the quadrature weights,
+                                  // FEValues::JxW() ). We can tell the
+                                  // FEValues object to do for
+                                  // us by also giving it the
+                                  // #update_quadrature_points
+                                  // flag:
+  FEValues<dim> fe_values (fe, quadrature_formula, 
+                          update_values   | update_gradients |
+                           update_quadrature_points | update_JxW_values);
+
+                                  // We then again define a few
+                                  // abbreviations. The values of these
+                                  // variables of course depend on the
+                                  // dimension which we are presently
+                                  // using. However, the FE and Quadrature
+                                  // classes do all the necessary work for
+                                  // you and you don't have to care about the
+                                  // dimension dependent parts:
+  const unsigned int   dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int   n_q_points    = quadrature_formula.size();
+
+  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+  Vector<double>       cell_rhs (dofs_per_cell);
+
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+                                   // Next, we again have to loop over all
+                                  // cells and assemble local contributions.
+                                  // Note, that a cell is a quadrilateral in
+                                  // two space dimensions, but a hexahedron
+                                  // in 3D. In fact, the
+                                  // <code>active_cell_iterator</code> data
+                                  // type is something different, depending
+                                  // on the dimension we are in, but to the
+                                  // outside world they look alike and you
+                                  // will probably never see a difference
+                                  // although the classes that this typedef
+                                  // stands for are in fact completely
+                                  // unrelated:
+  typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+  
+  for (; cell!=endc; ++cell)
+    {
+      fe_values.reinit (cell);
+      cell_matrix = 0;
+      cell_rhs = 0;
+
+                                      // Now we have to assemble the
+                                      // local matrix and right hand
+                                      // side. This is done exactly
+                                      // like in the previous
+                                      // example, but now we revert
+                                      // the order of the loops
+                                      // (which we can safely do
+                                      // since they are independent
+                                      // of each other) and merge the
+                                      // loops for the local matrix
+                                      // and the local vector as far
+                                      // as possible to make
+                                      // things a bit faster.
+                                       //
+                                       // Assembling the right hand side
+                                       // presents the only significant
+                                       // difference to how we did things in
+                                       // step-3: Instead of using a constant
+                                       // right hand side with value 1, we use
+                                       // the object representing the right
+                                       // hand side and evaluate it at the
+                                       // quadrature points:
+      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+       for (unsigned int i=0; i<dofs_per_cell; ++i)
+         {
+           for (unsigned int j=0; j<dofs_per_cell; ++j)
+             cell_matrix(i,j) += (fe_values.shape_grad (i, q_point) *
+                                  fe_values.shape_grad (j, q_point) *
+                                  fe_values.JxW (q_point));
+
+           cell_rhs(i) += (fe_values.shape_value (i, q_point) *
+                           right_hand_side.value (fe_values.quadrature_point (q_point)) *
+                           fe_values.JxW (q_point));
+         }
+                                       // As a final remark to these loops:
+                                       // when we assemble the local
+                                       // contributions into
+                                       // <code>cell_matrix(i,j)</code>, we
+                                       // have to multiply the gradients of
+                                       // shape functions $i$ and $j$ at point
+                                       // q_point and multiply it with the
+                                       // scalar weights JxW. This is what
+                                       // actually happens:
+                                       // <code>fe_values.shape_grad(i,q_point)</code>
+                                       // returns a <code>dim</code>
+                                       // dimensional vector, represented by a
+                                       // <code>Tensor@<1,dim@></code> object,
+                                       // and the operator* that multiplies it
+                                       // with the result of
+                                       // <code>fe_values.shape_grad(j,q_point)</code>
+                                       // makes sure that the <code>dim</code>
+                                       // components of the two vectors are
+                                       // properly contracted, and the result
+                                       // is a scalar floating point number
+                                       // that then is multiplied with the
+                                       // weights. Internally, this operator*
+                                       // makes sure that this happens
+                                       // correctly for all <code>dim</code>
+                                       // components of the vectors, whether
+                                       // <code>dim</code> be 2, 3, or any
+                                       // other space dimension; from a user's
+                                       // perspective, this is not something
+                                       // worth bothering with, however,
+                                       // making things a lot simpler if one
+                                       // wants to write code dimension
+                                       // independently.
+      
+                                      // With the local systems assembled,
+                                      // the transfer into the global matrix
+                                      // and right hand side is done exactly
+                                      // as before, but here we have again
+                                      // merged some loops for efficiency:
+      cell->get_dof_indices (local_dof_indices);
+//       for (unsigned int i=0; i<dofs_per_cell; ++i)
+//     {
+//       for (unsigned int j=0; j<dofs_per_cell; ++j)
+//         system_matrix.add (local_dof_indices[i],
+//                            local_dof_indices[j],
+//                            cell_matrix(i,j));
+//       
+//       system_rhs(local_dof_indices[i]) += cell_rhs(i);
+//     }
+       
+      constraints.distribute_local_to_global (cell_matrix, cell_rhs,
+                                              local_dof_indices,
+                                              system_matrix, system_rhs);
+    }
+  
+//                                // As the final step in this function, we
+//                                // wanted to have non-homogeneous boundary
+//                                // values in this example, unlike the one
+//                                // before. This is a simple task, we only
+//                                // have to replace the
+//                                // ZeroFunction used there by
+//                                // an object of the class which describes
+//                                // the boundary values we would like to use
+//                                // (i.e. the <code>BoundaryValues</code>
+//                                // class declared above):
+// 
+//   MatrixTools::apply_boundary_values (boundary_values,
+//                                   system_matrix,
+//                                   solution,
+//                                   system_rhs);
+}
+
+                                 // @sect4{Step4::projection_active_set}
+
+                                // Projection and updating of the active set
+                                 // for the dofs which penetrates the obstacle.
+template <int dim>
+void Step4<dim>::projection_active_set ()
+{
+//   const Obstacle<dim>     obstacle;
+//   std::vector<bool>       vertex_touched (triangulation.n_vertices(),
+//                                 false);
+// 
+//   boundary_values.clear ();
+//   VectorTools::interpolate_boundary_values (dof_handler,
+//                                         0,
+//                                         BoundaryValues<dim>(),
+//                                         boundary_values);
+// 
+//   typename DoFHandler<dim>::active_cell_iterator
+//     cell = dof_handler.begin_active(),
+//     endc = dof_handler.end();
+// 
+//   active_set = 0;
+//   unsigned int n = 0;
+//   for (; cell!=endc; ++cell)
+//     for (unsigned int v=0; v<GeometryInfo<2>::vertices_per_cell; ++v)
+//       {
+//     if (vertex_touched[cell->vertex_index(v)] == false)
+//       {
+//         vertex_touched[cell->vertex_index(v)] = true;
+//         unsigned int index_x = cell->vertex_dof_index (v,0);
+//         // unsigned int index_y = cell->vertex_dof_index (v,1);
+// 
+//         Point<dim> point (cell->vertex (v)[0], cell->vertex (v)[1]);
+//         double obstacle_value = obstacle.value (point);
+//         if (solution (index_x) >= obstacle_value && resid_vector (index_x) <= 0)
+//           {
+//             solution (index_x) = obstacle_value;
+//             boundary_values.insert (std::pair<unsigned int, double>(index_x, obstacle_value));
+//             active_set (index_x) = 1;
+//             n += 1;
+//           }
+//       }
+//       }
+//   std::cout<< "Number of active contraints: " << n <<std::endl;
+  
+  const Obstacle<dim>     obstacle;
+  std::vector<bool>       vertex_touched (triangulation.n_vertices(),
+                                   false);
+  typename DoFHandler<dim>::active_cell_iterator
+  cell = dof_handler.begin_active(),
+  endc = dof_handler.end();
+
+  constraints.clear();
+  active_set = 0;
+  for (; cell!=endc; ++cell)
+    for (unsigned int v=0; v<GeometryInfo<2>::vertices_per_cell; ++v)
+      {
+       unsigned int index_x = cell->vertex_dof_index (v,0);
+
+       Point<dim> point (cell->vertex (v)[0], cell->vertex (v)[1]);
+       double obstacle_value = obstacle.value (point);
+       if (solution (index_x) >= obstacle_value && resid_vector (index_x) <= 0)
+       {
+         constraints.add_line (index_x);
+         constraints.set_inhomogeneity (index_x, obstacle_value);
+         solution (index_x) = 0;
+         active_set (index_x) = 1;
+       }
+      }
+      
+  VectorTools::interpolate_boundary_values (dof_handler,
+                                           0,
+                                           BoundaryValues<dim>(),
+                                           constraints);
+  constraints.close ();
+}
+
+                                 // @sect4{Step4::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 Step4<dim>::solve () 
+{
+  ReductionControl        reduction_control (100, 1e-12, 1e-2);
+  SolverCG<>              solver (reduction_control);
+  SolverBicgstab<>        solver_bicgstab (reduction_control);
+  PreconditionSSOR<SparseMatrix<double> > precondition;
+  precondition.initialize (system_matrix, 1.2);
+
+  solver.solve (system_matrix, solution, system_rhs, precondition);
+
+  std::cout << "Initial error: " << reduction_control.initial_value() <<std::endl;
+  std::cout << "   " << reduction_control.last_step()
+           << " CG iterations needed to obtain convergence with an error: "
+           <<  reduction_control.last_value()
+           << std::endl;
+
+  constraints.distribute (solution);
+}
+
+                                 // @sect4{Step4::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
+                                 // VTK format, rather than for gnuplot. VTK
+                                 // format is currently the most widely used
+                                 // one and is supported by a number of
+                                 // visualization programs such as Visit and
+                                 // Paraview (for ways to obtain these
+                                 // programs 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_vtk</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
+                                 // <code>dim</code> can have other values
+                                 // than 2 or 3, but we neglect this here for
+                                 // the sake of brevity).
+template <int dim>
+void Step4<dim>::output_results (Vector<double> vector_to_plot, const std::string& title) const
+{
+  DataOut<dim> data_out;
+
+  data_out.attach_dof_handler (dof_handler);
+  data_out.add_data_vector (vector_to_plot, "vector_to_plot");
+
+  data_out.build_patches ();
+
+  std::ofstream output_vtk (dim == 2 ?
+                           (title + ".vtk").c_str () :
+                           (title + ".vtk").c_str ());
+  data_out.write_vtk (output_vtk);
+
+  std::ofstream output_gnuplot (dim == 2 ?
+                               (title + ".gp").c_str () :
+                               (title + ".gp").c_str ());
+  data_out.write_gnuplot (output_gnuplot);
+}
+
+
+
+                                 // @sect4{Step4::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 Step4<dim>::run () 
+{
+  std::cout << "Solving problem in " << dim << " space dimensions." << std::endl;
+
+  make_grid();
+  setup_system ();
+
+  constraints.clear ();
+  VectorTools::interpolate_boundary_values (dof_handler,
+                                           0,
+                                           BoundaryValues<dim>(),
+                                           constraints);
+  constraints.close ();
+  ConstraintMatrix constraints_complete (constraints);
+  assemble_system ();
+
+  system_matrix_complete.copy_from (system_matrix);
+  system_rhs_complete = system_rhs;
+
+  std::cout<< "Update Active Set:" <<std::endl;
+  solution = 0;
+  resid_vector = 0;
+  projection_active_set ();
+
+  for (unsigned int i=0; i<solution.size (); i++)
+    {
+//       std::ostringstream filename_matrix;
+//       filename_matrix << "system_matrix_";
+//       filename_matrix << i;
+//       filename_matrix << ".dat";
+//       std::ofstream matrix (filename_matrix.str ().c_str());
+
+      std::cout<< "Assemble System:" <<std::endl;
+      system_matrix = 0;
+      system_rhs = 0;
+      assemble_system ();
+//       constraints.print (matrix);
+//       system_matrix.print (matrix);
+//       for (unsigned int k=0; k<solution.size (); k++)
+//     std::cout<< system_rhs (k) << ", "
+//              << solution (k) << ", "
+//              << system_rhs.l2_norm ()
+//              <<std::endl;
+      std::cout<< "Solve System" <<std::endl;
+      solve ();
+
+      std::ostringstream filename_solution;
+      filename_solution << "solution_";
+      filename_solution << i;
+      output_results (solution, filename_solution.str ());
+
+      resid_vector = 0;
+      resid_vector -= system_rhs_complete;
+      system_matrix_complete.vmult_add  (resid_vector, solution);
+
+      for (unsigned int k = 0; k<solution.size (); k++)
+       if (resid_vector (k) < 0)
+         resid_vector (k) = 0;
+
+      std::ostringstream filename_residuum;
+      filename_residuum << "residuum_";
+      filename_residuum << i;
+      output_results (resid_vector, filename_residuum.str ());
+
+      std::ostringstream filename_active_set;
+      filename_active_set << "active_set_";
+      filename_active_set << i;
+      output_results (active_set, filename_active_set.str ());
+
+      double resid = resid_vector.l2_norm ();
+      std::cout<< i << ". Residuum = " << resid <<std::endl;
+      if (resid < 1e-10)
+       {
+         break;
+       }
+
+      std::cout<< "Update Active Set:"<<std::endl;
+      projection_active_set ();
+    }
+}
+
+
+                                 // @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>Step4@<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);
+  {
+    Step4<2> laplace_problem_2d;
+    laplace_problem_2d.run ();
+  }
+  
+  // {
+  //   Step4<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.