]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Change the commentaries
authorfrohne <frohne@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 9 Nov 2011 00:22:32 +0000 (00:22 +0000)
committerfrohne <frohne@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 9 Nov 2011 00:22:32 +0000 (00:22 +0000)
git-svn-id: https://svn.dealii.org/trunk@24739 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-41/step-41.cc

index cb8c34fb5e9db7623775ee3688228105a49f6ecf..c44c5f68baf25dd140459f96f10be2e07dece14c 100644 (file)
 
                                  // @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 <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.
+                                 // @sect3{The <code>Step41</code> class template}
+
+                                // This class supply all function and variables
+                                 // to an obstacle problem. The projection_active_set
+                                 // function and the ConstaintMatrix are important
+                                 // for the handling of the active set as we see
+                                 // later.
+
 template <int dim>
-class Step4 
+class Step41 
 {
   public:
-    Step4 ();
+    Step41 ();
     void run ();
     
   private:
@@ -119,55 +91,6 @@ class Step4
 
                                  // @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> 
 {
@@ -178,8 +101,6 @@ class RightHandSide : public Function<dim>
                          const unsigned int  component = 0) const;
 };
 
-
-
 template <int dim>
 class BoundaryValues : public Function<dim> 
 {
@@ -203,38 +124,8 @@ class Obstacle : public Function<dim>
 
 
                                 // 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++.
+                                // side function a constant force density
+                                 // like the gravitation attraction.
 template <int dim>
 double RightHandSide<dim>::value (const Point<dim> &p,
                                  const unsigned int /*component*/) const 
@@ -245,13 +136,7 @@ double RightHandSide<dim>::value (const Point<dim> &p,
 }
 
 
-                                // 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:
+                                // As boundary values, we choose the zero.
 template <int dim>
 double BoundaryValues<dim>::value (const Point<dim> &p,
                                   const unsigned int /*component*/) const 
@@ -261,6 +146,11 @@ double BoundaryValues<dim>::value (const Point<dim> &p,
   return return_value;
 }
 
+
+                                // The obstacle function describes a cascaded
+                                 // barrier. So if the gravitation attraction
+                                 // pulls the membrane down it blows over the
+                                 // steps.
 template <int dim>
 double Obstacle<dim>::value (const Point<dim> &p,
                             const unsigned int /*component*/) const 
@@ -281,86 +171,25 @@ double Obstacle<dim>::value (const Point<dim> &p,
 
 
 
-                                 // @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:
+                                 // @sect3{Implementation of the <code>Step41</code> class}
+            
+
+                                 // @sect4{Step41::Step41}
+
 template <int dim>
-Step4<dim>::Step4 ()
+Step41<dim>::Step41 ()
                :
                 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.
+                                 // @sect4{Step41::make_grid}
+
+                                 // We solve our obstacle problem on the square
+                                 // $[-1,1]\times [-1,1]$ in 2D.
 template <int dim>
-void Step4<dim>::make_grid ()
+void Step41<dim>::make_grid ()
 {
   GridGenerator::hyper_cube (triangulation, -1, 1);
   n_refinements = 6;
@@ -374,19 +203,10 @@ void Step4<dim>::make_grid ()
            << 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!
+                                 // @sect4{Step41::setup_system}
+
 template <int dim>
-void Step4<dim>::setup_system ()
+void Step41<dim>::setup_system ()
 {
   dof_handler.distribute_dofs (fe);
 
@@ -411,78 +231,27 @@ void Step4<dim>::setup_system ()
 }
 
 
-                                 // @sect4{Step4::assemble_system}
+                                 // @sect4{Step41::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.
+
+                                // At once with assembling the system matrix and
+                                 // right-hand-side we apply the constraints
+                                 // to our system. The constraint consists not
+                                 // only of the zero Dirichlet boundary values,
+                                 // in addition they contain the obstacle values.
+                                 // The projection_active_set function are used
+                                 // to fill the ConstraintMatrix.
 template <int dim>
-void Step4<dim>::assemble_system () 
+void Step41<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();
 
@@ -491,19 +260,6 @@ void Step4<dim>::assemble_system ()
 
   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();
@@ -514,28 +270,6 @@ void Step4<dim>::assemble_system ()
       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)
          {
@@ -548,57 +282,33 @@ void Step4<dim>::assemble_system ()
                            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);       
+
+      cell->get_dof_indices (local_dof_indices);
+
+                                       // This function apply the constraints
+                                       // to the system matrix and system rhs.
+                                       // The true parameter is set to make sure
+                                       // that the system rhs contains correct
+                                       // values in the rows with inhomogeneity
+                                       // constraints.
       constraints.distribute_local_to_global (cell_matrix, cell_rhs,
                                               local_dof_indices,
                                               system_matrix, system_rhs, true);
     }
 }
 
-                                 // @sect4{Step4::projection_active_set}
+                                 // @sect4{Step41::projection_active_set}
 
-                                // Projection and updating of the active set
-                                 // for the dofs which penetrates the obstacle.
+                                // Updating of the active set which means to
+                                 // set a inhomogeneity constraint in the
+                                 // ConstraintMatrix. At the same time we set
+                                 // the solution to the correct value - the obstacle value.
+                                 // To control the active set we use the vector
+                                 // active_set which contains a zero in a component
+                                 // that is not in the active set and elsewise a
+                                 // one. With the output file you can visualize it.
 template <int dim>
-void Step4<dim>::projection_active_set ()
+void Step41<dim>::projection_active_set ()
 {
   const Obstacle<dim>     obstacle;
   std::vector<bool>       vertex_touched (triangulation.n_vertices(),
@@ -624,10 +334,12 @@ void Step4<dim>::projection_active_set ()
        double obstacle_value = obstacle.value (point);
        double solution_index_x = solution (index_x);
 
-                                       // to decide which dof belongs to the
-                                       // active-set. for that we scale the
+                                       // To decide which dof belongs to the
+                                       // active-set. For that we scale the
                                        // residual-vector with the cell-size and
                                        // the diag-entry of the mass-matrix.
+
+                                       // TODO: I have to check the condition
        if ((resid_vector (index_x)*std::pow (2, 2*n_refinements)*diag_mass_matrix_vector (index_x) >= solution_index_x - obstacle_value))
        {
          constraints.add_line (index_x);
@@ -644,7 +356,7 @@ void Step4<dim>::projection_active_set ()
       }
   std::cout<< "Number of Contact-Constaints: " << counter_contact_constraints <<std::endl;
 
-                                       // to supply the boundary values of the
+                                       // To supply the boundary values of the
                                        // dirichlet-boundary in constraints
   VectorTools::interpolate_boundary_values (dof_handler,
                                            0,
@@ -653,17 +365,10 @@ void Step4<dim>::projection_active_set ()
   constraints.close ();
 }
 
-                                 // @sect4{Step4::solve}
+                                 // @sect4{Step41::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 () 
+void Step41<dim>::solve () 
 {
   ReductionControl        reduction_control (100, 1e-12, 1e-2);
   SolverCG<TrilinosWrappers::Vector>   solver (reduction_control); 
@@ -679,35 +384,13 @@ void Step4<dim>::solve ()
            << std::endl;
 }
 
-                                 // @sect4{Step4::output_results}
+                                 // @sect4{Step41::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).
+                                // We use the vtk-format for the output.
+                                 // The file contains the displacement,
+                                 // the residual and active set vectors.
 template <int dim>
-void Step4<dim>::output_results (const std::string& title) const
+void Step41<dim>::output_results (const std::string& title) const
 {
   DataOut<dim> data_out;
   
@@ -718,23 +401,22 @@ void Step4<dim>::output_results (const std::string& title) const
 
   data_out.build_patches ();
 
-  std::ofstream output_vtk (dim == 2 ?
-                           (title + ".vtk").c_str () :
-                           (title + ".vtk").c_str ());
+  std::ofstream output_vtk ((title + ".vtk").c_str ());
   data_out.write_vtk (output_vtk);
 }
 
 
 
-                                 // @sect4{Step4::run}
+                                 // @sect4{Step41::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.
+                                // top-level control over everything.
+                                 // Here the active set method is implemented.
+
+                                 // TODO: I have to compare it with the algorithm
+                                 // in the Wohlmuth-paper
 template <int dim>
-void Step4<dim>::run () 
+void Step41<dim>::run () 
 {
   std::cout << "Solving problem in " << dim << " space dimensions." << std::endl;
 
@@ -812,7 +494,7 @@ void Step4<dim>::run ()
                                 // 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
+                                // <code>Step41@<2@></code> (forcing
                                 // the compiler to compile the class template
                                 // with <code>dim</code> replaced by
                                 // <code>2</code>) and run a 2d simulation,
@@ -887,7 +569,7 @@ int main (int argc, char *argv[])
 
   Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv);
 
-  Step4<2> laplace_problem_2d;
+  Step41<2> laplace_problem_2d;
   laplace_problem_2d.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.