]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Document half of the program.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 6 Jul 2009 22:18:58 +0000 (22:18 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Mon, 6 Jul 2009 22:18:58 +0000 (22:18 +0000)
git-svn-id: https://svn.dealii.org/trunk@19039 0785d39b-7218-0410-832d-ea1e28bc413d

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

index a6de2659975ddd3b15028fca00f230ecf5b76e96..472adcd01a74808a0c0679d6fd69c41be65cfddc 100644 (file)
 /*    to the file deal.II/doc/license.html for the  text  and     */
 /*    further information on this license.                        */
 
-                                // @sect3{Include files} We bundle the
-                                // "usual" deal.II include files as we
-                                // did in step-4:
+                                // @sect3{Include files}
+
+                                // As mentioned in the introduction, this
+                                // program is essentially only a slightly
+                                // revised version of step-4. As a
+                                // consequence, most of the following include
+                                // files are as used there, or at least as
+                                // used already in previous tutorial
+                                // programs:
 #include <base/logstream.h>
 #include <base/quadrature_lib.h>
 #include <base/function.h>
 #include <numerics/matrices.h>
 #include <numerics/data_out.h>
 
-                                // PETSc appears here because SLEPc
-                                // depends on this library:
+                                // PETSc appears here because SLEPc
+                                // depends on this library:
 #include <lac/petsc_sparse_matrix.h>
 #include <lac/petsc_vector.h>
 
-                                // and then we need to actually import
-                                // the interfaces for solvers that
-                                // SLEPc provides:
+                                // And then we need to actually import
+                                // the interfaces for solvers that
+                                // SLEPc provides:
 #include <lac/slepc_solver.h>
 
-                                // and some standard C++:
+                                // We also need some standard C++:
 #include <fstream>
 #include <iostream>
 
-                               // and the finally, as in previous
-                               // programs, we import all the deal.II
-                               // class and function names into the
-                               // global namespace:
+                                // Finally, as in previous programs, we
+                                // import all the deal.II class and function
+                                // names into the global namespace:
 using namespace dealii;
 
-                                // @sect1{The
-                                // <code>EigenvalueProblem</code>
-                                // class template}
+                                // @sect3{The <code>EigenvalueProblem</code> class template}
 
+                                // Following is the class declaration for the
+                                // main class template. It looks pretty much
+                                // exactly like what has already been shown
+                                // in step-4:
 template <int dim>
 class EigenvalueProblem 
 {
-public:
-  EigenvalueProblem (const std::string &prm_file);
-  void run ();
+  public:
+    EigenvalueProblem (const std::string &prm_file);
+    void run ();
   
-private:
-  void make_grid_and_dofs ();
-  void assemble_system ();
-  void solve ();
-  void output_results () const;
+  private:
+    void make_grid_and_dofs ();
+    void assemble_system ();
+    void solve ();
+    void output_results () const;
   
-  Triangulation<dim> triangulation;
-  FE_Q<dim>          fe;
-  DoFHandler<dim>    dof_handler;
-  
-                                // These are data types pertaining to
-                                // the generalized problem.
-  PETScWrappers::SparseMatrix        stiffness_matrix, mass_matrix;
-  std::vector<PETScWrappers::Vector> eigenfunctions;
-  std::vector<double>                eigenvalues;   
-  
-  ParameterHandler parameters;
+    Triangulation<dim> triangulation;
+    FE_Q<dim>          fe;
+    DoFHandler<dim>    dof_handler;
+
+                                    // With these exceptions: For our
+                                    // eigenvalue problem, we need both a
+                                    // stiffness matrix for the left hand
+                                    // side as well as a mass matrix for the
+                                    // right hand side. We also need not just
+                                    // one solution function, but a whole set
+                                    // of those for the eigenfunctions we
+                                    // want to compute, along with the
+                                    // corresponding eigenvectors:
+    PETScWrappers::SparseMatrix        stiffness_matrix, mass_matrix;
+    std::vector<PETScWrappers::Vector> eigenfunctions;
+    std::vector<double>                eigenvalues;   
+
+                                    // And then we need an object that will
+                                    // store several run-time parameters that
+                                    // we will specify in an input file:
+    ParameterHandler parameters;
+
+                                    // Finally, we will have an object that
+                                    // contains "constraints" on our degrees
+                                    // of freedom. This could include hanging
+                                    // node constraints if we had adaptively
+                                    // refined meshes (which we don't have in
+                                    // the current program). Here, we will
+                                    // store the constraints for boundary
+                                    // nodes $U_i=0$.
+    ConstraintMatrix constraints;
 };
 
-                                // @sect2{Implementation of the
-                                // <code>EigenvalueProblem</code>
-                                // class}
+                                // @sect3{Implementation of the <code>EigenvalueProblem</code> class}
 
-                                // @sect3{EigenvalueProblem::EigenvalueProblem}
+                                // @sect4{EigenvalueProblem::EigenvalueProblem}
 
+                                // First up, the constructor. The main, new
+                                // part is handling the run-time input
+                                // parameters. We need to declare their
+                                // existence first, and then read their
+                                // values from the input file whose name is
+                                // specified as an argument to this function:
 template <int dim>
 EigenvalueProblem<dim>::EigenvalueProblem (const std::string &prm_file)
-  :
-  fe (1),
-  dof_handler (triangulation)
+               :
+               fe (1),
+               dof_handler (triangulation)
 {
-
-                                // Declare some of the needed data
-                                // from a file; you can always change
-                                // this!
   parameters.declare_entry ("Global mesh refinement steps", "5",
                            Patterns::Integer (0, 20),
                            "The number number of times the 1-cell coarse mesh should "
@@ -112,13 +139,64 @@ EigenvalueProblem<dim>::EigenvalueProblem (const std::string &prm_file)
                            Patterns::Anything(),
                            "A functional description of the potential.");
   
-                                // Entries are declared, so now we
-                                // read them in...
   parameters.read_input (prm_file);
 }
 
 
-                                // @sect3{EigenvalueProblem::make_grid_and_dofs}
+                                // @sect4{EigenvalueProblem::make_grid_and_dofs}
+
+                                // The next function creates a mesh on the
+                                // domain $[-1,1]^d$, refines it as many
+                                // times as the input file calls for, and
+                                // then attaches a DoFHandler to it and
+                                // initializes the matrices and vectors to
+                                // their correct sizes. We also build the
+                                // constraints that correspond to the
+                                // boundary values $u|_{\partial\Omega}=0$.
+                                //
+                                // For the matrices, we use the PETSc
+                                // wrappers. These have the ability to
+                                // allocate memory as necessary as non-zero
+                                // entries are added. This seems inefficient:
+                                // we could as well first compute the
+                                // sparsity pattern, initialize the matrices
+                                // with it, and as we then insert entries we
+                                // can be sure that we do not need to
+                                // re-allocate memory and free the one used
+                                // previously. One way to do that would be to
+                                // use code like this:
+                                // @code
+                                //   CompressedSimpleSparsityPattern
+                                //      csp (dof_handler.n_dofs(),
+                                //           dof_handler.n_dofs());
+                                //   DoFTools::make_sparsity_pattern (dof_handler, csp);
+                                //   csp.compress ();
+                                //   stiffness_matrix.reinit (csp);
+                                //   mass_matrix.reinit (csp);
+                                // @code
+                                // instead of the two <code>reinit()</code>
+                                // calls for the stiffness and mass matrices.
+                                //
+                                // This doesn't quite work,
+                                // unfortunately. The code above may lead to
+                                // a few entries in the non-zero pattern to
+                                // which we only ever write zero entries;
+                                // most notably, this holds true for
+                                // off-diagonal entries for those rows and
+                                // columns that belong to boundary
+                                // nodes. This shouldn't be a problem, but
+                                // for whatever reason, PETSc's ILU
+                                // preconditioner, which we use to solve
+                                // linear systems in the eigenvalue solver,
+                                // doesn't like these extra entries and
+                                // aborts with an error message.
+                                //
+                                // Absent any obvious way to avoid this, we
+                                // simply settle for the second best option,
+                                // which is have PETSc allocate memory as
+                                // necessary. That said, since this is not a
+                                // time critical part, this whole affair is
+                                // of no further importance.
 template <int dim>
 void EigenvalueProblem<dim>::make_grid_and_dofs ()
 {
@@ -126,50 +204,50 @@ void EigenvalueProblem<dim>::make_grid_and_dofs ()
   triangulation.refine_global (parameters.get_integer ("Global mesh refinement steps"));
   dof_handler.distribute_dofs (fe);
 
-  CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(),
-                                              dof_handler.n_dofs());
-  DoFTools::make_sparsity_pattern (dof_handler, csp);
-  csp.compress ();
-
-  // What is going on here?
-
-  // This does not work!
-  //   stiffness_matrix.reinit (csp);
-  //   mass_matrix.reinit (csp);
+  DoFTools::make_zero_boundary_constraints (dof_handler, constraints);
+  constraints.close ();
   
-  // But this does... TODO: Fix it!
-  stiffness_matrix.reinit (dof_handler.n_dofs(), dof_handler.n_dofs(),
+  stiffness_matrix.reinit (dof_handler.n_dofs(),
+                          dof_handler.n_dofs(),
                           dof_handler.max_couplings_between_dofs());
-  mass_matrix.reinit (dof_handler.n_dofs(), dof_handler.n_dofs(),
+  mass_matrix.reinit (dof_handler.n_dofs(),
+                     dof_handler.n_dofs(),
                      dof_handler.max_couplings_between_dofs());
 
-                                // with this done we stream-out the
-                                // sparsity pattern
-  std::ofstream out ("constrained_sparsity_pattern.gpl");
-  csp.print_gnuplot (out);
-  
-                                // The next step is to take care of
-                                // the eigenspectrum. In this case,
-                                // the outputs are eigenfunctions and
-                                // eigenvalues. Set the collective
-                                // eigenfunction block to be as big as
-                                // we wanted!
+                                  // The next step is to take care of the
+                                  // eigenspectrum. In this case, the outputs
+                                  // are eigenfunctions and eigenvalues, so
+                                  // we set the size of the list of
+                                  // eigenfunctions and eigenvalues to be as
+                                  // large as asked for in the input file:
   eigenfunctions
     .resize (parameters.get_integer ("Number of eigenvalues/eigenfunctions"));
   for (unsigned int i=0; i<eigenfunctions.size (); ++i)
     eigenfunctions[i].reinit (dof_handler.n_dofs ());
 
-                                // and do the same for the eigenvalue
-                                // output, which had better be the
-                                // same size as the eigenfunction
-                                // block
-  eigenvalues
-    .resize (eigenfunctions.size ());
-
+  eigenvalues.resize (eigenfunctions.size ());
 }
 
 
-                                // @sect3{EigenvalueProblem::assemble_system}
+                                // @sect4{EigenvalueProblem::assemble_system}
+
+                                // Here, we assemble the global stiffness and
+                                // mass matrices from local contributions
+                                // $A^K_{ij} = \int_K \nabla\varphi_i(\mathbf
+                                // x) \cdot \nabla\varphi_j(\mathbf x) +
+                                // V(\mathbf x)\varphi_i(\mathbf
+                                // x)\varphi_j(\mathbf x)$. The function
+                                // should be immediately familiar if you've
+                                // seen previous tutorial programs. The only
+                                // thing new would be setting up an object
+                                // that described the potential $V(\mathbf
+                                // x)$ using the expression that we got from
+                                // the input file. We then need to evaluate
+                                // this object at the quadrature points on
+                                // each cell. If you've seen how to evaluate
+                                // function objects (see, for example the
+                                // coefficient in step-5), the code here will
+                                // also look rather familiar.
 template <int dim>
 void EigenvalueProblem<dim>::assemble_system () 
 {  
@@ -191,20 +269,9 @@ void EigenvalueProblem<dim>::assemble_system ()
   potential.initialize (FunctionParser<dim>::default_variable_names (),
                        parameters.get ("Potential"),
                        typename FunctionParser<dim>::ConstMap());
-
-                                // Here we call the actual
-                                // quadrature-point values of the
-                                // potential defined in the prm file
-                                // which we read in earlier
+  
   std::vector<double> potential_values (n_q_points);
   
-                                // Initialize objects denoting zero
-                                // boundary constraints for the
-                                // present grid.
-  ConstraintMatrix constraints;
-  constraints.clear ();
-  DoFTools::make_zero_boundary_constraints (dof_handler, constraints);
-  constraints.close ();
   
   typename DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active (),
@@ -218,12 +285,9 @@ void EigenvalueProblem<dim>::assemble_system ()
       potential.value_list (fe_values.get_quadrature_points(),
                            potential_values);
       
-      for (unsigned int q_point=0; q_point<n_q_points; 
-          ++q_point)
-       for (unsigned int i=0; i<dofs_per_cell; 
-            ++i)
-         for (unsigned int j=0; j<dofs_per_cell; 
-              ++j)
+      for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+       for (unsigned int i=0; i<dofs_per_cell; ++i)
+         for (unsigned int j=0; j<dofs_per_cell; ++j)
            {
              cell_stiffness_matrix (i, j)
                += (fe_values.shape_grad (i, q_point) *
@@ -232,18 +296,18 @@ void EigenvalueProblem<dim>::assemble_system ()
                    potential_values[q_point] *
                    fe_values.shape_value (i, q_point) *
                    fe_values.shape_value (j, q_point)
-                   ) * fe_values.JxW (q_point);
+               ) * fe_values.JxW (q_point);
              
              cell_mass_matrix (i, j)
                += (fe_values.shape_value (i, q_point) *
                    fe_values.shape_value (j, q_point) 
-                   ) * fe_values.JxW (q_point);
+               ) * fe_values.JxW (q_point);
            }
 
-                               // Now we have the local system we
-                               // transfer it into the global objects
-                               // and take care of zero boundary
-                               // constraints,
+                                      // Now that we have the local matrix
+                                      // contributions, we transfer them into
+                                      // the global objects and take care of
+                                      // zero boundary constraints:
       cell->get_dof_indices (local_dof_indices);
       
       constraints
@@ -256,55 +320,43 @@ void EigenvalueProblem<dim>::assemble_system ()
                                     mass_matrix);
     }
 
-                                // and finally, set the matrices and
-                                // vectors in an assembled state.
+                                  // At the end of the function, we tell
+                                  // PETSc that the matrices have now been
+                                  // fully assembled and that the sparse
+                                  // matrix representation can now be
+                                  // compressed as no more entries will be
+                                  // added:
   stiffness_matrix.compress();
   mass_matrix.compress();
-
-                               // make sure that the diagonal entries
-                               // of constrained degrees of freedom
-                               // are non-zero to ensure that the
-                               // matrix is actually invertible
-  //   for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
-  //     if (constraints.is_constrained(i))
-  //       {
-  //   stiffness_matrix.set (i, i, 1);
-  //   mass_matrix.set (i, i, 1);
-  //       }
-  
-                                // finally set the matrices in an
-                                // assembled state so SLEPc likes:
-  //   stiffness_matrix.compress ();
-  //   mass_matrix.compress ();  
 }
 
 
-                                // @sect3{EigenvalueProblem::solve}
-                                // Now that the system is set up, here
-                                // is a good time to actually solve
-                                // the problem: As with other examples
-                                // this is done using a "solve"
-                                // routine
+                                // @sect4{EigenvalueProblem::solve}
+                                // Now that the system is set up, here
+                                // is a good time to actually solve
+                                // the problem: As with other examples
+                                // this is done using a "solve"
+                                // routine
 template <int dim>
 void EigenvalueProblem<dim>::solve () 
 {
-                                // We start by assigning the accuracy
-                                // to which we would like to solve the
-                                // system,
+                                  // We start by assigning the accuracy
+                                  // to which we would like to solve the
+                                  // system,
   SolverControl solver_control (1000, 1e-6);
 
-                                // and assign our solver of
-                                // choice. Here we want to use the
-                                // Krylov-Schur solver, which is
-                                // pretty darn fast and robust:
+                                  // and assign our solver of
+                                  // choice. Here we want to use the
+                                  // Krylov-Schur solver, which is
+                                  // pretty darn fast and robust:
   SLEPcWrappers::SolverKrylovSchur eigensolver (solver_control);
 
-                                // Lets assign the solver which part
-                                // of the spectrum we want to solve
+                                  // Lets assign the solver which part
+                                  // of the spectrum we want to solve
   eigensolver.set_which_eigenpairs (EPS_SMALLEST_MAGNITUDE);
 
-                                // Finally, we actually solve the
-                                // generalized eigenproblem:
+                                  // Finally, we actually solve the
+                                  // generalized eigenproblem:
   eigensolver.solve (stiffness_matrix, mass_matrix, 
                     eigenvalues, eigenfunctions, 
                     eigenfunctions.size());
@@ -318,7 +370,7 @@ void EigenvalueProblem<dim>::solve ()
 }
 
 
-                                // @sect3{EigenvalueProblem::output_results}
+                                // @sect4{EigenvalueProblem::output_results}
 template <int dim>
 void EigenvalueProblem<dim>::output_results () const
 {
@@ -331,7 +383,7 @@ void EigenvalueProblem<dim>::output_results () const
                              std::string("eigenfunction_") +
                              Utilities::int_to_string(i));
 
-                                // How does this work?
+                                  // How does this work?
   Vector<double> projected_potential (dof_handler.n_dofs());
   FunctionParser<dim> potential;
   potential.initialize (FunctionParser<dim>::default_variable_names (),
@@ -347,7 +399,7 @@ void EigenvalueProblem<dim>::output_results () const
 }
 
 
-                                 // @sect3{EigenvalueProblem::run}
+                                 // @sect4{EigenvalueProblem::run}
 
                                  // This is the function which has the
                                 // top-level control over
@@ -360,9 +412,9 @@ void EigenvalueProblem<dim>::run ()
   
   make_grid_and_dofs ();
 
-                                 // While we are here, lets count the
-                                 // number of active cells and degrees
-                                 // of freedom like we always do.
+                                  // While we are here, lets count the
+                                  // number of active cells and degrees
+                                  // of freedom like we always do.
   std::cout << "   Number of active cells:       "
            << triangulation.n_active_cells()
            << std::endl
@@ -386,10 +438,11 @@ int main (int argc, char **argv)
   try
     {
 
-                                 // Here is another difference from
-                                 // other steps: We initialize the
-                                 // SLEPc work space which inherently
-                                 // initializes the PETSc work space
+                                      // Here is another difference from
+                                      // other steps: We initialize the SLEPc
+                                      // work space which inherently
+                                      // initializes the PETSc work space,
+                                      // run the whole program, ...
       SlepcInitialize (&argc,&argv,0,0);
 
       {
@@ -399,12 +452,14 @@ int main (int argc, char **argv)
        problem.run ();
       }
 
-                                 // and then unitialize the SLEPc
-                                 // work space when the job is done:
+                                      // ...and then unitialize the SLEPc
+                                      // work space when the job is done:
       SlepcFinalize ();
     }
 
-                                 // or panic if something goes wrong.
+                                  // All the while, we are watching out if
+                                  // any exceptions should have been
+                                  // generated. If that is so, we panic...
   catch (std::exception &exc)
     {
       std::cerr << std::endl << std::endl
@@ -430,10 +485,9 @@ int main (int argc, char **argv)
       return 1;
     }
   
-                                 // or show that we happy and didn't
-                                 // crap out on the calculation...
+                                  // ...or show that we are happy:
   std::cout << std::endl 
-           << "Completed
+           << "Done.
            << std::endl;
 
   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.