]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Clean up a tiny bit.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 16 May 2012 12:47:00 +0000 (12:47 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 16 May 2012 12:47:00 +0000 (12:47 +0000)
git-svn-id: https://svn.dealii.org/trunk@25511 0785d39b-7218-0410-832d-ea1e28bc413d

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

index c5bfa569b8f2b3a945a69ec7e9c4c9316de8b80c..e0d0d63225fd9003baf8f35d16712bce307e8f8d 100644 (file)
@@ -1,9 +1,9 @@
 /* $Id$ */
-/* Author: Wolfgang Bangerth, University of Heidelberg, 2000 */
+/* Author: Sven Wetterauer, University of Heidelberg, 2012 */
 
 /*    $Id$       */
 /*                                                                */
-/*    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2010, 2011 by the deal.II authors */
+/*    Copyright (C) 2012 by the deal.II authors */
 /*                                                                */
 /*    This file is subject to QPL and may not be  distributed     */
 /*    without copyright and license information. Please refer     */
 
                                  // @sect3{Include files}
 
-                                // The first few files have already
-                                // been covered in previous examples
-                                // and will thus not be further
-                                // commented on.
+                                 // The first few files have already
+                                 // been covered in previous examples
+                                 // and will thus not be further
+                                 // commented on.
 #include <deal.II/base/quadrature_lib.h>
 #include <deal.II/base/function.h>
 #include <deal.II/base/logstream.h>
 #include <deal.II/base/utilities.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/precondition.h>
 #include <deal.II/lac/constraint_matrix.h>
+
 #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/grid/tria_boundary_lib.h>
 #include <deal.II/grid/grid_refinement.h>
+
+#include <deal.II/dofs/dof_handler.h>
 #include <deal.II/dofs/dof_accessor.h>
 #include <deal.II/dofs/dof_tools.h>
+
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/fe/fe_q.h>
+
 #include <deal.II/numerics/vectors.h>
 #include <deal.II/numerics/matrices.h>
 #include <deal.II/numerics/data_out.h>
 #include <fstream>
 #include <iostream>
 
-                               // We will use adaptive mesh refinement between Newton
-                               // interations. To do so, we need to be able to work
-                               // with a solution on the new mesh, although it was
-                               // computed on the old one. The SolutionTransfer
-                               // class transfers the solution to the new mesh.
+                                 // We will use adaptive mesh refinement between Newton
+                                 // interations. To do so, we need to be able to work
+                                 // with a solution on the new mesh, although it was
+                                 // computed on the old one. The SolutionTransfer
+                                 // class transfers the solution to the new mesh.
 
 #include <deal.II/numerics/solution_transfer.h>
 
-                               // In this tutorial, we can't use the CG-method as a solver, as
-                               // described above, but we use the minimal residual method, which
-                               // is included with this file.
-
-#include <deal.II/lac/solver_minres.h>
+                                 // In this tutorial, we can't use the CG-method as a solver, as
+                                 // described above, but we use the minimal residual method, which
+                                 // is included with this file.
+#include <deal.II/lac/solver_cg.h>
+//#include <deal.II/lac/solver_minres.h>
 
-                               // As in previous programs:
-
-using namespace dealii;
+                                 // As in previous programs:
+namespace Step15
+{
+  using namespace dealii;
 
 
-                                                               // @sect3{The <code>Step15</code> class template}
+                                  // @sect3{The <code>MinimalSurfaceProblem</code> class template}
 
-                               // The class template is basically the same as in step 6.
-                               // Four additions are made: There are two solution vectors,
-                               // one for the Newton update, and one for the solution of
-                               // the original pde. Also we need a double for the residual
-                               // of the Newton method, an integer, which counts the mesh
-                               // refinements and a bool for the boundary condition in the first
-                               // Newton step.
+                                  // The class template is basically the same as in step 6.
+                                  // Four additions are made: There are two solution vectors,
+                                  // one for the Newton update, and one for the solution of
+                                  // the original pde. Also we need a double for the residual
+                                  // of the Newton method, an integer, which counts the mesh
+                                  // refinements and a bool for the boundary condition in the first
+                                  // Newton step.
 
-template <int dim>
-class Step15
-{
-  public:
-    Step15 ();
-    ~Step15 ();
+  template <int dim>
+  class MinimalSurfaceProblem
+  {
+    public:
+      MinimalSurfaceProblem ();
+      ~MinimalSurfaceProblem ();
 
-    void run ();
+      void run ();
 
-  private:
-    void setup_system ();
-    void assemble_system ();
-    void solve ();
-    void refine_grid ();
-    void set_boundary_values ();
-    double compute_residual (const double alpha) const;
-    double determine_step_length() const;
+    private:
+      void setup_system ();
+      void assemble_system ();
+      void solve ();
+      void refine_grid ();
+      void set_boundary_values ();
+      double compute_residual (const double alpha) const;
+      double determine_step_length() const;
 
-    Triangulation<dim>   triangulation;
+      Triangulation<dim>   triangulation;
 
-    DoFHandler<dim>      dof_handler;
-    FE_Q<dim>            fe;
+      DoFHandler<dim>      dof_handler;
+      FE_Q<dim>            fe;
 
-    ConstraintMatrix     hanging_node_constraints;
+      ConstraintMatrix     hanging_node_constraints;
 
-    SparsityPattern      sparsity_pattern;
-    SparseMatrix<double> system_matrix;
+      SparsityPattern      sparsity_pattern;
+      SparseMatrix<double> system_matrix;
 
-    Vector<double>       present_solution;
-    Vector<double>       newton_update;
-    Vector<double>       system_rhs;
+      Vector<double>       present_solution;
+      Vector<double>       newton_update;
+      Vector<double>       system_rhs;
 
 
 
-    unsigned int         refinement;
+      unsigned int         refinement;
 
-                       // As described in the Introduction, the first Newton iteration
-                       // is special, because of the boundary condition. To implement
-                       // these correctly, there is a bool, which is true in the first
-                       // step and false ever after.
-    bool                                first_step;
-};
+                                      // As described in the Introduction, the first Newton iteration
+                                      // is special, because of the boundary condition. To implement
+                                      // these correctly, there is a bool, which is true in the first
+                                      // step and false ever after.
+      bool                                 first_step;
+  };
 
-                                                               // @sect3{Boundary condition}
+                                  // @sect3{Boundary condition}
 
-                               // The boundary condition is implemented just like in step 4.
-                               // It was chosen as $g(x,y)=sin(2 \pi (x+y))$ in this example.
+                                  // The boundary condition is implemented just like in step 4.
+                                  // It was chosen as $g(x,y)=sin(2 \pi (x+y))$ in this example.
 
-template <int dim>
-class BoundaryValues : public Function<dim>
-{
-  public:
-    BoundaryValues () : Function<dim>() {}
+  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;
-};
+      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
-{
-  double return_value=sin(2*M_PI*(p[0]+p[1]));
-  return return_value;
-}
+  template <int dim>
+  double BoundaryValues<dim>::value (const Point<dim> &p,
+                                    const unsigned int /*component*/) const
+  {
+    double return_value=sin(2*M_PI*(p[0]+p[1]));
+    return return_value;
+  }
 
-                                                               // @sect3{The <code>Step15</code> class implementation}
+                                  // @sect3{The <code>MinimalSurfaceProblem</code> class implementation}
 
-                                                               // @sect4{Step15::Step15}
+                                  // @sect4{MinimalSurfaceProblem::MinimalSurfaceProblem}
 
-                               // The constructor and destructor of the class are the same
-                               // as in the first few tutorials.
+                                  // The constructor and destructor of the class are the same
+                                  // as in the first few tutorials.
 
-template <int dim>
-Step15<dim>::Step15 ()
-               :
-               dof_handler (triangulation),
-                fe (2)
-{}
+  template <int dim>
+  MinimalSurfaceProblem<dim>::MinimalSurfaceProblem ()
+                 :
+                 dof_handler (triangulation),
+                 fe (2)
+  {}
 
 
 
-                               //
-template <int dim>
-Step15<dim>::~Step15 ()
-{
-  dof_handler.clear ();
-}
-
-                                                               // @sect4{Step15::setup_system}
-
-                               // As always in the setup-system function, we setup the variables
-                               // of the finite element method. There are same differences to
-                               // step 6, because we don't have to solve one pde over all,
-                               // but one in every Newton step. Also the starting function
-                               // has to be setup in the first step.
+                                  //
+  template <int dim>
+  MinimalSurfaceProblem<dim>::~MinimalSurfaceProblem ()
+  {
+    dof_handler.clear ();
+  }
 
-template <int dim>
-void Step15<dim>::setup_system ()
-{
+                                  // @sect4{MinimalSurfaceProblem::setup_system}
 
-                               // This function will be called, every time we refine the mesh
-                               // to resize the system matrix,  Newton update - and right hand
-                               // side vector and to set the right values of hanging nodes to
-                               // get a continuous solution.
-                               // But only the first time, the starting solution has to be
-                               // initialized. Also the vector of the solution will be
-                               // resized in the <code>refine_grid</code> function, while the
-                               // vector is transfered to the new mesh.
+                                  // As always in the setup-system function, we setup the variables
+                                  // of the finite element method. There are same differences to
+                                  // step 6, because we don't have to solve one pde over all,
+                                  // but one in every Newton step. Also the starting function
+                                  // has to be setup in the first step.
 
-  if(first_step)
+  template <int dim>
+  void MinimalSurfaceProblem<dim>::setup_system ()
   {
+
+                                    // This function will be called, every time we refine the mesh
+                                    // to resize the system matrix,  Newton update - and right hand
+                                    // side vector and to set the right values of hanging nodes to
+                                    // get a continuous solution.
+                                    // But only the first time, the starting solution has to be
+                                    // initialized. Also the vector of the solution will be
+                                    // resized in the <code>refine_grid</code> function, while the
+                                    // vector is transfered to the new mesh.
+
+    if (first_step)
+      {
        dof_handler.distribute_dofs (fe);
        present_solution.reinit (dof_handler.n_dofs());
        for(unsigned int i=0; i<dof_handler.n_dofs();++i)
-       {
-         present_solution(i)=0;
-       }
-                               // The constraint matrix, holding a list of the hanging nodes,
-                               // will be setup in the <code>refine_grid</code> function
-                               // after refining the mesh.
+         {
+           present_solution(i)=0;
+         }
+                                        // The constraint matrix, holding a list of the hanging nodes,
+                                        // will be setup in the <code>refine_grid</code> function
+                                        // after refining the mesh.
 
        hanging_node_constraints.clear ();
        DoFTools::make_hanging_node_constraints (dof_handler,
-                                  hanging_node_constraints);
+                                                hanging_node_constraints);
        hanging_node_constraints.close ();
+      }
+
+
+                                    // The remaining parts of the function are the same as in step 6.
+
+    newton_update.reinit (dof_handler.n_dofs());
+    system_rhs.reinit (dof_handler.n_dofs());
+
+    CompressedSparsityPattern c_sparsity(dof_handler.n_dofs());
+    DoFTools::make_sparsity_pattern (dof_handler, c_sparsity);
+
+    hanging_node_constraints.condense (c_sparsity);
+
+    sparsity_pattern.copy_from(c_sparsity);
+    system_matrix.reinit (sparsity_pattern);
   }
 
+                                  // @sect4{MinimalSurfaceProblem::assemble_system}
 
-                           // The remaining parts of the function are the same as in step 6.
+                                  // This function does the same as in the previous tutorials.
+                                  // The only additional step is the correct implementation of
+                                  // the boundary condition and the usage of the gradients of
+                                  // the old solution.
 
-  newton_update.reinit (dof_handler.n_dofs());
-  system_rhs.reinit (dof_handler.n_dofs());
+  template <int dim>
+  void MinimalSurfaceProblem<dim>::assemble_system ()
+  {
+    const QGauss<dim>  quadrature_formula(3);
 
-  CompressedSparsityPattern c_sparsity(dof_handler.n_dofs());
-  DoFTools::make_sparsity_pattern (dof_handler, c_sparsity);
+    system_matrix = 0;
+    system_rhs = 0;
 
-  hanging_node_constraints.condense (c_sparsity);
+    FEValues<dim> fe_values (fe, quadrature_formula,
+                            update_gradients |
+                            update_quadrature_points  |  update_JxW_values);
 
-  sparsity_pattern.copy_from(c_sparsity);
-  system_matrix.reinit (sparsity_pattern);
-}
+    const unsigned int   dofs_per_cell = fe.dofs_per_cell;
+    const unsigned int   n_q_points    = quadrature_formula.size();
 
-                                                               // @sect4{Step15::assemble_system}
+    FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+    Vector<double>       cell_rhs (dofs_per_cell);
 
-                               // This function does the same as in the previous tutorials.
-                               // The only additional step is the correct implementation of
-                               // the boundary condition and the usage of the gradients of
-                               // the old solution.
+    std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 
-template <int dim>
-void Step15<dim>::assemble_system ()
-{
-  const QGauss<dim>  quadrature_formula(3);
+    typename DoFHandler<dim>::active_cell_iterator
+      cell = dof_handler.begin_active(),
+      endc = dof_handler.end();
+    for (; cell!=endc; ++cell)
+      {
+       cell_matrix = 0;
+       cell_rhs = 0;
 
-  system_matrix = 0;
-  system_rhs = 0;
+       fe_values.reinit (cell);
 
-  FEValues<dim> fe_values (fe, quadrature_formula,
-                          update_gradients |
-                          update_quadrature_points  |  update_JxW_values);
 
-  const unsigned int   dofs_per_cell = fe.dofs_per_cell;
-  const unsigned int   n_q_points    = quadrature_formula.size();
+       for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
 
-  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
-  Vector<double>       cell_rhs (dofs_per_cell);
+                                          // To setup up the linear system, the gradient of the old solution
+                                          // in the quadrature points is needed. For this purpose there is
+                                          // is a function, which will write these gradients in a vector,
+                                          // where every component of the vector is a vector itself:
 
-  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+         std::vector<Tensor<1, dim> > gradients(n_q_points);
+         fe_values.get_function_gradients(present_solution, gradients);
 
-  typename DoFHandler<dim>::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
-  for (; cell!=endc; ++cell)
-    {
-      cell_matrix = 0;
-      cell_rhs = 0;
-
-      fe_values.reinit (cell);
-
-
-      for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
-
-                 // To setup up the linear system, the gradient of the old solution
-                 // in the quadrature points is needed. For this purpose there is
-                 // is a function, which will write these gradients in a vector,
-                 // where every component of the vector is a vector itself:
-
-                       std::vector<Tensor<1, dim> > gradients(n_q_points);
-                       fe_values.get_function_gradients(present_solution, gradients);
-
-                         // Having the gradients of the old solution in the quadrature
-                         // points, we are able to compute the coefficients $a_{n}$
-                         // in these points.
-
-                       const double coeff = 1/sqrt(1 + gradients[q_point] * gradients[q_point]);
-
-                         // The assembly of the system then is the same as always, except
-                         // of the damping parameter of the Newton method, which we set on
-                         // 0.1 in this case.
-
-                       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)
-                                                       * coeff
-                                                       * (fe_values.shape_grad(j, q_point)
-                                                               - coeff * coeff
-                                                               * (fe_values.shape_grad(j, q_point)
-                                                               * gradients[q_point])
-                                                               * gradients[q_point])
-                                                       * fe_values.JxW(q_point));
-                               }
-
-                               cell_rhs(i) -= (fe_values.shape_grad(i, q_point) * coeff
-                                                               * gradients[q_point] * fe_values.JxW(q_point));
-                       }
-               }
-
-      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);
+                                          // Having the gradients of the old solution in the quadrature
+                                          // points, we are able to compute the coefficients $a_{n}$
+                                          // in these points.
+
+         const double coeff = 1/sqrt(1 + gradients[q_point] * gradients[q_point]);
+
+                                          // The assembly of the system then is the same as always, except
+                                          // of the damping parameter of the Newton method, which we set on
+                                          // 0.1 in this case.
+
+         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)
+                                   * coeff
+                                   * (fe_values.shape_grad(j, q_point)
+                                      - coeff * coeff
+                                      * (fe_values.shape_grad(j, q_point)
+                                         * gradients[q_point])
+                                      * gradients[q_point])
+                                   * fe_values.JxW(q_point));
+           }
+
+           cell_rhs(i) -= (fe_values.shape_grad(i, q_point) * coeff
+                           * gradients[q_point] * fe_values.JxW(q_point));
+         }
        }
-    }
-  hanging_node_constraints.condense (system_matrix);
-  hanging_node_constraints.condense (system_rhs);
-  std::map<unsigned int,double> boundary_values;
+
+       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);
+         }
+      }
+    hanging_node_constraints.condense (system_matrix);
+    hanging_node_constraints.condense (system_rhs);
+    std::map<unsigned int,double> boundary_values;
 
 
     VectorTools::interpolate_boundary_values (dof_handler,
-                                             0,
-                                             ZeroFunction<dim>(),
-                                             boundary_values);
-
-  MatrixTools::apply_boundary_values (boundary_values,
-                                     system_matrix,
-                                     newton_update,
-                                     system_rhs);
-}
+                                             0,
+                                             ZeroFunction<dim>(),
+                                             boundary_values);
+
+    MatrixTools::apply_boundary_values (boundary_values,
+                                       system_matrix,
+                                       newton_update,
+                                       system_rhs);
+  }
 
 
-template <int dim>
-double Step15<dim>::compute_residual (const double alpha) const
-{
-  const QGauss<dim>  quadrature_formula(3);
+  template <int dim>
+  double MinimalSurfaceProblem<dim>::compute_residual (const double alpha) const
+  {
+    const QGauss<dim>  quadrature_formula(3);
 
-  Vector<double> residual (dof_handler.n_dofs());
+    Vector<double> residual (dof_handler.n_dofs());
 
-  Vector<double> linearization_point (dof_handler.n_dofs());
-  linearization_point = present_solution;
-  linearization_point.add (alpha, newton_update);
+    Vector<double> linearization_point (dof_handler.n_dofs());
+    linearization_point = present_solution;
+    linearization_point.add (alpha, newton_update);
 
-  FEValues<dim> fe_values (fe, quadrature_formula,
-                           update_gradients |
-                           update_quadrature_points  |  update_JxW_values);
+    FEValues<dim> fe_values (fe, quadrature_formula,
+                            update_gradients |
+                            update_quadrature_points  |  update_JxW_values);
 
-  const unsigned int   dofs_per_cell = fe.dofs_per_cell;
-  const unsigned int   n_q_points    = quadrature_formula.size();
+    const unsigned int   dofs_per_cell = fe.dofs_per_cell;
+    const unsigned int   n_q_points    = quadrature_formula.size();
 
-  Vector<double>       cell_rhs (dofs_per_cell);
+    Vector<double>       cell_rhs (dofs_per_cell);
 
-  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+    std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 
-  typename DoFHandler<dim>::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
-  for (; cell!=endc; ++cell)
-    {
-      cell_rhs = 0;
+    typename DoFHandler<dim>::active_cell_iterator
+      cell = dof_handler.begin_active(),
+      endc = dof_handler.end();
+    for (; cell!=endc; ++cell)
+      {
+       cell_rhs = 0;
 
-      fe_values.reinit (cell);
+       fe_values.reinit (cell);
 
 
-      for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
+       for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) {
 
-                  // To setup up the linear system, the gradient of the old solution
-                  // in the quadrature points is needed. For this purpose there is
-                  // is a function, which will write these gradients in a vector,
-                  // where every component of the vector is a vector itself:
+                                          // To setup up the linear system, the gradient of the old solution
+                                          // in the quadrature points is needed. For this purpose there is
+                                          // is a function, which will write these gradients in a vector,
+                                          // where every component of the vector is a vector itself:
 
-                        std::vector<Tensor<1, dim> > gradients(n_q_points);
-                        fe_values.get_function_gradients(linearization_point, gradients);
+         std::vector<Tensor<1, dim> > gradients(n_q_points);
+         fe_values.get_function_gradients(linearization_point, gradients);
 
-                          // Having the gradients of the old solution in the quadrature
-                          // points, we are able to compute the coefficients $a_{n}$
-                          // in these points.
+                                          // Having the gradients of the old solution in the quadrature
+                                          // points, we are able to compute the coefficients $a_{n}$
+                                          // in these points.
 
-                        const double coeff = 1/sqrt(1 + gradients[q_point] * gradients[q_point]);
+         const double coeff = 1/sqrt(1 + gradients[q_point] * gradients[q_point]);
 
-                          // The assembly of the system then is the same as always, except
-                          // of the damping parameter of the Newton method, which we set on
-                          // 0.1 in this case.
+                                          // The assembly of the system then is the same as always, except
+                                          // of the damping parameter of the Newton method, which we set on
+                                          // 0.1 in this case.
 
-                        for (unsigned int i = 0; i < dofs_per_cell; ++i) {
-                                cell_rhs(i) -= (fe_values.shape_grad(i, q_point) * coeff
-                                                                * gradients[q_point] * fe_values.JxW(q_point));
-                        }
-                }
+         for (unsigned int i = 0; i < dofs_per_cell; ++i) {
+           cell_rhs(i) -= (fe_values.shape_grad(i, q_point) * coeff
+                           * gradients[q_point] * fe_values.JxW(q_point));
+         }
+       }
 
-      cell->get_dof_indices (local_dof_indices);
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-          residual(local_dof_indices[i]) += cell_rhs(i);
-    }
-  hanging_node_constraints.condense (residual);
-
-  std::map<unsigned int,double> boundary_values;
-  VectorTools::interpolate_boundary_values (dof_handler,
-                                              0,
-                                              ZeroFunction<dim>(),
-                                              boundary_values);
-  for (std::map<unsigned int,double>::const_iterator p = boundary_values.begin();
-      p != boundary_values.end(); ++p)
-    residual(p->first) = 0;
-
-  return residual.l2_norm();
-}
+       cell->get_dof_indices (local_dof_indices);
+       for (unsigned int i=0; i<dofs_per_cell; ++i)
+         residual(local_dof_indices[i]) += cell_rhs(i);
+      }
+    hanging_node_constraints.condense (residual);
 
-                                                               // @sect4{Step15::solve}
+    std::map<unsigned int,double> boundary_values;
+    VectorTools::interpolate_boundary_values (dof_handler,
+                                             0,
+                                             ZeroFunction<dim>(),
+                                             boundary_values);
+    for (std::map<unsigned int,double>::const_iterator p = boundary_values.begin();
+        p != boundary_values.end(); ++p)
+      residual(p->first) = 0;
+
+    return residual.l2_norm();
+  }
 
-                               // The solve function is the same as always, we just have to
-                               // implement the minimal residual method as a solver and
-                               // apply the Newton update to the solution.
+                                  // @sect4{MinimalSurfaceProblem::solve}
 
-template <int dim>
-void Step15<dim>::solve ()
-{
-  SolverControl           solver_control (1000, system_rhs.l2_norm()*1e-6);
-  SolverMinRes<>              solver (solver_control);
+                                  // The solve function is the same as always, we just have to
+                                  // implement the minimal residual method as a solver and
+                                  // apply the Newton update to the solution.
+
+  template <int dim>
+  void MinimalSurfaceProblem<dim>::solve ()
+  {
+    SolverControl solver_control (system_rhs.size(),
+                                 system_rhs.l2_norm()*1e-6);
+    SolverCG<>    solver (solver_control);
 
-  PreconditionSSOR<> preconditioner;
-  preconditioner.initialize(system_matrix, 1.2);
+    PreconditionSSOR<> preconditioner;
+    preconditioner.initialize(system_matrix, 1.2);
 
-  solver.solve (system_matrix, newton_update, system_rhs,
-               preconditioner);
+    solver.solve (system_matrix, newton_update, system_rhs,
+                 preconditioner);
 
-  hanging_node_constraints.distribute (newton_update);
+    hanging_node_constraints.distribute (newton_update);
 
-                                 // In this step, the old solution is updated to the new one:
-  const double alpha = determine_step_length();
-  std::cout << "  step length alpha=" << alpha << std::endl;
-  present_solution.add (alpha, newton_update);
-}
+                                    // In this step, the old solution is updated to the new one:
+    const double alpha = determine_step_length();
+    std::cout << "  step length alpha=" << alpha << std::endl;
+    present_solution.add (alpha, newton_update);
+  }
 
 
-template <int dim>
-double Step15<dim>::determine_step_length() const
-{
-  return 0.1;
-}
-                                                               // @sect4{Step15::refine_grid}
+  template <int dim>
+  double MinimalSurfaceProblem<dim>::determine_step_length() const
+  {
+    return 0.1;
+  }
+                                  // @sect4{MinimalSurfaceProblem::refine_grid}
 
-                               // The first part of this function is the same as in step 6.
-                               // But after refining the mesh we have to transfer the old
-                               // solution to the new one, which is done with the help of
-                               // the SolutionTransfer class.
+                                  // The first part of this function is the same as in step 6.
+                                  // But after refining the mesh we have to transfer the old
+                                  // solution to the new one, which is done with the help of
+                                  // the SolutionTransfer class.
 
 
-template <int dim>
-void Step15<dim>::refine_grid ()
-{
-  Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
-
-  KellyErrorEstimator<dim>::estimate (dof_handler,
-                                     QGauss<dim-1>(3),
-                                     typename FunctionMap<dim>::type(),
-                                     present_solution,
-                                     estimated_error_per_cell);
-
-  GridRefinement::refine_and_coarsen_fixed_number (triangulation,
-                                                  estimated_error_per_cell,
-                                                  0.3, 0.03);
-
-                       // Then we need an additional step: if, for example,
-                               // you flag a cell that is once more refined than its neighbor,
-                               // and that neighbor is not flagged for refinement, we would end
-                               // up with a jump of two refinement levels across a cell interface.
-                               // To avoid these situations, the library will
-                               // silently also have to refine the neighbor cell once. It does so
-                               // by calling the Triangulation::prepare_coarsening_and_refinement
-                               // function before actually doing the refinement and coarsening.
-                               // This function flags a set of additional cells for refinement or
-                               // coarsening, to enforce rules like the one-hanging-node rule.
-                               // The cells that are flagged for refinement and coarsening after
-                               // calling this function are exactly the ones that will actually
-                               // be refined or coarsened. Since the SolutionTransfer class needs
-                               // this information in order to store the data from the old mesh
-                               // and transfer to the new one.
-
-  triangulation.prepare_coarsening_and_refinement ();
-
-                               // With this out of the way, we initialize a SolutionTransfer
-                               // object with the present DoFHandler and attach the solution
-                               // vector to it:
-
-  SolutionTransfer<dim> solution_transfer(dof_handler);
-  solution_transfer.prepare_for_coarsening_and_refinement(present_solution);
-
-                               // Then we do the actual refinement, and distribute degrees
-                               // of freedom on the new mesh:
-
-  triangulation.execute_coarsening_and_refinement();
-  dof_handler.distribute_dofs(fe);
-
-                               // Finally, we retrieve the old solution interpolated to the new
-                               // mesh. Since the SolutionTransfer function does not actually
-                               // store the values of the old solution, but rather indices, we
-                               // need to preserve the old solution vector until we have gotten
-                               // the new interpolated values. Thus, we have the new values
-                               // written into a temporary vector, and only afterwards write
-                               // them into the solution vector object:
-
-  Vector<double> tmp(dof_handler.n_dofs());
-  solution_transfer.interpolate(present_solution,tmp);
-  present_solution=tmp;
-
-  set_boundary_values ();
-
-                               // On the new mesh, there are different hanging nodes, which shall
-                               // be enlisted in a matrix like before. To ensure there are no
-                               // hanging nodes of the old mesh in the matrix, it's first cleared:
-  hanging_node_constraints.clear();
-
-                               // After doing so, the hanging nodes of the new mesh can be
-                               // enlisted in the matrix, like before. Calling the
-                               // <code>setup_system</code> function in the <code>run</code>
-                               // function again after this, the hanging nodes don't have to
-                               // be enlisted there once more.
-
-  DoFTools::make_hanging_node_constraints(dof_handler, hanging_node_constraints);
-  hanging_node_constraints.close();
-  hanging_node_constraints.distribute(present_solution);
-}
+  template <int dim>
+  void MinimalSurfaceProblem<dim>::refine_grid ()
+  {
+    Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+
+    KellyErrorEstimator<dim>::estimate (dof_handler,
+                                       QGauss<dim-1>(3),
+                                       typename FunctionMap<dim>::type(),
+                                       present_solution,
+                                       estimated_error_per_cell);
+
+    GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+                                                    estimated_error_per_cell,
+                                                    0.3, 0.03);
+
+                                    // Then we need an additional step: if, for example,
+                                    // you flag a cell that is once more refined than its neighbor,
+                                    // and that neighbor is not flagged for refinement, we would end
+                                    // up with a jump of two refinement levels across a cell interface.
+                                    // To avoid these situations, the library will
+                                    // silently also have to refine the neighbor cell once. It does so
+                                    // by calling the Triangulation::prepare_coarsening_and_refinement
+                                    // function before actually doing the refinement and coarsening.
+                                    // This function flags a set of additional cells for refinement or
+                                    // coarsening, to enforce rules like the one-hanging-node rule.
+                                    // The cells that are flagged for refinement and coarsening after
+                                    // calling this function are exactly the ones that will actually
+                                    // be refined or coarsened. Since the SolutionTransfer class needs
+                                    // this information in order to store the data from the old mesh
+                                    // and transfer to the new one.
+
+    triangulation.prepare_coarsening_and_refinement ();
+
+                                    // With this out of the way, we initialize a SolutionTransfer
+                                    // object with the present DoFHandler and attach the solution
+                                    // vector to it:
+
+    SolutionTransfer<dim> solution_transfer(dof_handler);
+    solution_transfer.prepare_for_coarsening_and_refinement(present_solution);
+
+                                    // Then we do the actual refinement, and distribute degrees
+                                    // of freedom on the new mesh:
+
+    triangulation.execute_coarsening_and_refinement();
+    dof_handler.distribute_dofs(fe);
+
+                                    // Finally, we retrieve the old solution interpolated to the new
+                                    // mesh. Since the SolutionTransfer function does not actually
+                                    // store the values of the old solution, but rather indices, we
+                                    // need to preserve the old solution vector until we have gotten
+                                    // the new interpolated values. Thus, we have the new values
+                                    // written into a temporary vector, and only afterwards write
+                                    // them into the solution vector object:
+
+    Vector<double> tmp(dof_handler.n_dofs());
+    solution_transfer.interpolate(present_solution,tmp);
+    present_solution=tmp;
+
+    set_boundary_values ();
+
+                                    // On the new mesh, there are different hanging nodes, which shall
+                                    // be enlisted in a matrix like before. To ensure there are no
+                                    // hanging nodes of the old mesh in the matrix, it's first cleared:
+    hanging_node_constraints.clear();
+
+                                    // After doing so, the hanging nodes of the new mesh can be
+                                    // enlisted in the matrix, like before. Calling the
+                                    // <code>setup_system</code> function in the <code>run</code>
+                                    // function again after this, the hanging nodes don't have to
+                                    // be enlisted there once more.
+
+    DoFTools::make_hanging_node_constraints(dof_handler, hanging_node_constraints);
+    hanging_node_constraints.close();
+    hanging_node_constraints.distribute(present_solution);
+  }
 
 
-template <int dim>
-void Step15<dim>::set_boundary_values ()
-{
-    // Having refined the mesh, there might be new nodal points on
-    // the boundary. These have just interpolated values, but
-    // not the right boundary values. This is fixed up, by
-    // setting all boundary nodals explicit to the right value:
+  template <int dim>
+  void MinimalSurfaceProblem<dim>::set_boundary_values ()
+  {
+                                    // Having refined the mesh, there might be new nodal points on
+                                    // the boundary. These have just interpolated values, but
+                                    // not the right boundary values. This is fixed up, by
+                                    // setting all boundary nodals explicit to the right value:
 
     std::map<unsigned int, double> boundary_values2;
     VectorTools::interpolate_boundary_values(dof_handler, 0,
-        BoundaryValues<dim>(), boundary_values2);
+                                            BoundaryValues<dim>(), boundary_values2);
     for (std::map<unsigned int, double>::const_iterator p =
-        boundary_values2.begin(); p != boundary_values2.end(); ++p)
+          boundary_values2.begin(); p != boundary_values2.end(); ++p)
       present_solution(p->first) = p->second;
-}
-                                                               // @sect4{Step15::run}
+  }
+                                  // @sect4{MinimalSurfaceProblem::run}
 
-                               // In the run function, the first grid is build. Also in this
-                               // function, the Newton iteration is implemented.
+                                  // In the run function, the first grid is build. Also in this
+                                  // function, the Newton iteration is implemented.
 
-template <int dim>
-void Step15<dim>::run ()
-{
+  template <int dim>
+  void MinimalSurfaceProblem<dim>::run ()
+  {
 
-                               // The integer refinement counts the mesh refinements. Obviously
-                               // starting the program, it should be zero.
-       refinement=0;
-       first_step=true;
+                                    // The integer refinement counts the mesh refinements. Obviously
+                                    // starting the program, it should be zero.
+    refinement=0;
+    first_step=true;
 
-                               // As described in the introduction, the domain is a unitball around
-                               // the origin. The Mesh is globally refined two times, not to start
-                               // on the coarse mesh, which consists only of five cells.
+                                    // As described in the introduction, the domain is a unitball around
+                                    // the origin. The Mesh is globally refined two times, not to start
+                                    // on the coarse mesh, which consists only of five cells.
 
-       GridGenerator::hyper_ball (triangulation);
-       static const HyperBallBoundary<dim> boundary;
-       triangulation.set_boundary (0, boundary);
-       triangulation.refine_global(2);
+    GridGenerator::hyper_ball (triangulation);
+    static const HyperBallBoundary<dim> boundary;
+    triangulation.set_boundary (0, boundary);
+    triangulation.refine_global(2);
 
-                               // The Newton iteration starts here. During the first step, there is
-                               // no residual computed, so the bool is needed here to enter the
-                               // iteration scheme. Later the Newton method will continue until the
-                               // residual is less than $10^{-3}$.
+                                    // The Newton iteration starts here. During the first step, there is
+                                    // no residual computed, so the bool is needed here to enter the
+                                    // iteration scheme. Later the Newton method will continue until the
+                                    // residual is less than $10^{-3}$.
 
-       double previous_res = 0;
-       while(first_step || (previous_res>1e-3))
-       {
+    double previous_res = 0;
+    while(first_step || (previous_res>1e-3))
+      {
 
-                               // In the first step, we compute the solution on the two times globally 
-                               // refined mesh. After that the mesh will be refined 
-                               // adaptively, in order to not get too many cells. The refinement 
-                               // is the first thing done every time we restart the process in the while-loop.
-               if(!first_step)
-               {
-                       refine_grid();
+                                        // In the first step, we compute the solution on the two times globally
+                                        // refined mesh. After that the mesh will be refined
+                                        // adaptively, in order to not get too many cells. The refinement
+                                        // is the first thing done every time we restart the process in the while-loop.
+       if(!first_step)
+         {
+           refine_grid();
 
-                       std::cout<<"********mesh-refinement:"<<refinement+1<<" ********"<<std::endl;
-                       refinement++;
-               }
+           std::cout<<"********mesh-refinement:"<<refinement+1<<" ********"<<std::endl;
+           refinement++;
+         }
 
 
-                               // First thing to do after refining the mesh, is to setup the vectors,
-                               // matrices, etc., which is done in the <code>setup_system</code>
-                               // function.
+                                        // First thing to do after refining the mesh, is to setup the vectors,
+                                        // matrices, etc., which is done in the <code>setup_system</code>
+                                        // function.
 
-               setup_system();
+       setup_system();
 
-               if (first_step)
-                 set_boundary_values ();
+       if (first_step)
+         set_boundary_values ();
 
-                               // On every mesh there are done five Newton steps, in order to get a
-                               // better solution, before the mesh gets too fine and the computations
-                               // take more time.
-                std::cout<<"initial residual:"<<compute_residual(0)<<std::endl;
+                                        // On every mesh there are done five Newton steps, in order to get a
+                                        // better solution, before the mesh gets too fine and the computations
+                                        // take more time.
+       std::cout<<"initial residual:"<<compute_residual(0)<<std::endl;
 
-               for(unsigned int i=0; i<5;++i)
-               {
+       for(unsigned int i=0; i<5;++i)
+         {
 
-                               // In every Newton step the system matrix and the right hand side
-                               // have to be computed.
+                                            // In every Newton step the system matrix and the right hand side
+                                            // have to be computed.
 
-                       assemble_system ();
-                        previous_res = system_rhs.l2_norm();
+           assemble_system ();
+           previous_res = system_rhs.l2_norm();
 
-                       solve ();
-                       first_step=false;
-                       std::cout<<"residual:"<<compute_residual(0)<<std::endl;
-               }
+           solve ();
+           first_step=false;
+           std::cout<<"residual:"<<compute_residual(0)<<std::endl;
+         }
 
-                           // The fifth solution, as well as the Newton update,
-                               // on every mesh will be written in a vtk-file,
-                               // in order to show the convergence of the solution.
+                                        // The fifth solution, as well as the Newton update,
+                                        // on every mesh will be written in a vtk-file,
+                                        // in order to show the convergence of the solution.
 
-               Assert (refinement < 100, ExcNotImplemented());
+       Assert (refinement < 100, ExcNotImplemented());
 
-               DataOut<dim> data_out;
+       DataOut<dim> data_out;
 
-               data_out.attach_dof_handler (dof_handler);
-               data_out.add_data_vector (newton_update, "update");
-               data_out.add_data_vector (present_solution, "solution");
-               data_out.build_patches ();
-               const std::string filename = "solution-" + Utilities::int_to_string (refinement, 2) + ".vtk";
-               std::ofstream output (filename.c_str());
-               data_out.write_vtk (output);
+       data_out.attach_dof_handler (dof_handler);
+       data_out.add_data_vector (newton_update, "update");
+       data_out.add_data_vector (present_solution, "solution");
+       data_out.build_patches ();
+       const std::string filename = "solution-" + Utilities::int_to_string (refinement, 2) + ".vtk";
+       std::ofstream output (filename.c_str());
+       data_out.write_vtk (output);
 
-       }
+      }
+  }
 }
 
-                                                               // @ sect4{The main function}
+                                // @ sect4{The main function}
 
-                               // Finally the main function, this follows the scheme of all other main
-                               // functions:
+                                // Finally the main function, this follows the scheme of all other main
+                                // functions:
 
 int main ()
 {
-
   try
     {
+      using namespace dealii;
+      using namespace Step15;
+
       deallog.depth_console (0);
 
-      Step15<2> laplace_problem_2d;
+      MinimalSurfaceProblem<2> laplace_problem_2d;
       laplace_problem_2d.run ();
     }
   catch (std::exception &exc)

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.