]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
step 16
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 1 Feb 2004 22:27:46 +0000 (22:27 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 1 Feb 2004 22:27:46 +0000 (22:27 +0000)
git-svn-id: https://svn.dealii.org/trunk@8379 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/tutorial/chapter-2.step-by-step/step-16.data/intro.html
deal.II/examples/step-16/Makefile
deal.II/examples/step-16/step-16.cc

index c59cfdf6092aac8d5c180aab3ea01f0da3d3390c..39d27081c8010009c1c0d523d1cb26c0bad56cc9 100644 (file)
@@ -1,2 +1,34 @@
 <a name="Intro"></a>
 <h1>Introduction</h1>
+
+<p>
+This example shows the basic usage of the multilevel functions in
+<tt>deal.II</tt>. It solves Helmholtz equation with Neumann boundary conditions
+to avoid additional complications due to Dirichlet boundary conditions (there,
+some library functionsare missing). In all other respects, it is similar to
+step 5.
+</p>
+
+<p> In order to allow sufficient flexibility in conjunction with systems of
+differential equations and block preconditioners, quite a few different objects
+have to be created before starting the multilevel method. These are
+<ul>
+<li>MGTransfer, the object handling transfer between grids
+<li>MGCoarse, the solver on the coarsest level
+<li>MGSmoother, the smoother on all other levels
+<li>MGMatrix, the matrix object having a special level multipplication.
+</ul>
+</p>
+
+<p>
+These objects are combined in an object of type Multigrid, containing the
+implementation of the V-cycle, which is in turn used by the preconditioner
+PreconditionMG, ready for plug-in into a linear slver of the LAC library.
+</p>
+
+<p>
+The multilevel method in <tt>deal.II</tt> follows in many respects the outlines
+of the various publications by James Bramble, Joseph Pasciak and ??? Xu. In
+order to understand many of the options, a rough familiarity with their work is
+quite helpful.
+</p>
\ No newline at end of file
index 433c82b233dad1b64781e624af8cb4791f946717..dae158698b9557610c74b83bd5727f4731733202 100644 (file)
@@ -1,9 +1,5 @@
 # $Id$
 
-# Protect the file from being made automatically as long as it is not tested
-default:
-       @echo Example not ready for compiling
-
 # For the small projects Makefile, you basically need to fill in only
 # four fields.
 #
index 1c6a1d1fd175f2dabb401baefa3b2cfe12e15108..73f3c5f51b58f442dd6ecd812b8c586c39600c1b 100644 (file)
                                  // DoFTools class are declared:
 #include <multigrid/mg_dof_tools.h>
 
+#include <multigrid/mg_coarse.h>
+#include <multigrid/mg_smoother.h>
+#include <multigrid/mg_matrix.h>
+
                                 // This is C++ ... see step 5 for
                                 // further comments.
 #include <fstream>
@@ -151,8 +155,14 @@ void LaplaceProblem<dim>::setup_system ()
                                   // revision). Remark, that the
                                   // finest level is nlevels-1.
   const unsigned int nlevels = triangulation.n_levels();
-  mg_sparsity.resize(0, nlevels-1);
+                                  // We first have to resize the
+                                  // container holding the
+                                  // SparseMatrix classes, since they
+                                  // have to release their
+                                  // SparsityPattern before it can be
+                                  // destroyed.
   mg_matrices.resize(0, nlevels-1);
+  mg_sparsity.resize(0, nlevels-1);
 
                                   // Now, we have to build a matrix
                                   // on each level. Technically, we
@@ -238,9 +248,11 @@ void LaplaceProblem<dim>::assemble_system ()
        for (unsigned int i=0; i<dofs_per_cell; ++i)
          {
            for (unsigned int j=0; j<dofs_per_cell; ++j)
-             cell_matrix(i,j) += (fe_values.shape_grad(i,q_point)
-                                  * fe_values.shape_grad(j,q_point)
-                                  * fe_values.JxW(q_point));
+             cell_matrix(i,j) += ((fe_values.shape_grad(i,q_point)
+                                   * fe_values.shape_grad(j,q_point))
+                                   + fe_values.shape_value(i,q_point)
+                                   * fe_values.shape_value(j,q_point))
+                                  * fe_values.JxW(q_point);
 
                                             // For the right hand
                                             // side, a constant value
@@ -263,15 +275,15 @@ void LaplaceProblem<dim>::assemble_system ()
     };
 
                                   // Again use zero boundary values:
-  std::map<unsigned int,double> boundary_values;
-  VectorTools::interpolate_boundary_values (mg_dof_handler,
-                                           0,
-                                           ZeroFunction<dim>(),
-                                           boundary_values);
-  MatrixTools::apply_boundary_values (boundary_values,
-                                     system_matrix,
-                                     solution,
-                                     system_rhs);
+//   std::map<unsigned int,double> boundary_values;
+//   VectorTools::interpolate_boundary_values (mg_dof_handler,
+//                                         0,
+//                                         ZeroFunction<dim>(),
+//                                         boundary_values);
+//   MatrixTools::apply_boundary_values (boundary_values,
+//                                   system_matrix,
+//                                   solution,
+//                                   system_rhs);
 }
 
 
@@ -280,7 +292,7 @@ void LaplaceProblem<dim>::assemble_system ()
                                 // the same as above. Only the loop
                                 // goes over all existing cells now
                                 // and the results must be entered
-                                // into the right matrix.
+                                // into the correct matrix.
 
                                 // Since we only do multi-level
                                 // preconditioning, no right-hand
@@ -291,7 +303,8 @@ void LaplaceProblem<dim>::assemble_multigrid ()
   QGauss2<dim>  quadrature_formula;
 
   FEValues<dim> fe_values (fe, quadrature_formula, 
-                          UpdateFlags(update_gradients |
+                          UpdateFlags(update_values |
+                                      update_gradients |
                                       update_q_points  |
                                       update_JxW_values));
 
@@ -303,8 +316,8 @@ void LaplaceProblem<dim>::assemble_multigrid ()
   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 
                                   // 
-  typename DoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
-                                         endc = mg_dof_handler.end();
+  typename MGDoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
+                                           endc = mg_dof_handler.end();
   for (; cell!=endc; ++cell)
     {
                                       // Remember the level of the
@@ -323,9 +336,11 @@ void LaplaceProblem<dim>::assemble_multigrid ()
        for (unsigned int i=0; i<dofs_per_cell; ++i)
          {
            for (unsigned int j=0; j<dofs_per_cell; ++j)
-             cell_matrix(i,j) += (fe_values.shape_grad(i,q_point)
-                                  * fe_values.shape_grad(j,q_point)
-                                  * fe_values.JxW(q_point));
+             cell_matrix(i,j) += ((fe_values.shape_grad(i,q_point)
+                                   * fe_values.shape_grad(j,q_point))
+                                  + fe_values.shape_value(i,q_point)
+                                  * fe_values.shape_value(j,q_point))
+                                 * fe_values.JxW(q_point);
          };
 
 
@@ -352,15 +367,15 @@ void LaplaceProblem<dim>::assemble_multigrid ()
     };
 
                                   // Again use zero boundary values:
-  std::map<unsigned int,double> boundary_values;
-  VectorTools::interpolate_boundary_values (mg_dof_handler,
-                                           0,
-                                           ZeroFunction<dim>(),
-                                           boundary_values);
-  MatrixTools::apply_boundary_values (boundary_values,
-                                     system_matrix,
-                                     solution,
-                                     system_rhs);
+//   std::map<unsigned int,double> boundary_values;
+//   VectorTools::interpolate_boundary_values (mg_dof_handler,
+//                                         0,
+//                                         ZeroFunction<dim>(),
+//                                         boundary_values);
+//   MatrixTools::apply_boundary_values (boundary_values,
+//                                   system_matrix,
+//                                   solution,
+//                                   system_rhs);
 }
 
 
@@ -374,48 +389,84 @@ void LaplaceProblem<dim>::assemble_multigrid ()
 template <int dim>
 void LaplaceProblem<dim>::solve () 
 {
+                                  // Create a memory handler for
+                                  // regular vectors. Note, that
+                                  // GrowingVectorMemory is more time
+                                  // efficient than
+                                  // PrimitiveVectorMemory.
+  GrowingVectorMemory<>   vector_memory;
+
+                                  // Now, create an object handling
+                                  // the transfer of functions
+                                  // between different grid
+                                  // levels.
+  MGTransferPrebuilt<double> mg_transfer;
+  mg_transfer.build_matrices(mg_dof_handler);
+
+                                  // Next, we need a coarse grid
+                                  // solver. Since our coarse grid is
+                                  // VERY coarse, we decide for a
+                                  // direct solver, even if its
+                                  // implementation is not very
+                                  // clever.
+  FullMatrix<float> coarse_matrix;
+  coarse_matrix.copy_from (mg_matrices[0]);
+  MGCoarseGridHouseholder<float, Vector<double> > mg_coarse;
+  mg_coarse.initialize(coarse_matrix);
+
+                                  // The final ingredient for the
+                                  // multilevel preconditioner is the
+                                  // smoother. It is very customary
+                                  // to use a relaxation method
+                                  // here. Names are getting quite
+                                  // long here, so we help with
+                                  // typedefs.
+  typedef PreconditionSOR<SparseMatrix<float> > RELAXATION;
+  MGSmootherRelaxation<SparseMatrix<float>, RELAXATION, Vector<double> >
+    mg_smoother(vector_memory);
+
+                                  // Initialize the smoother with our
+                                  // level matrices and the required,
+                                  // additional data for the
+                                  // relaxaton method with default
+                                  // values.
+  RELAXATION::AdditionalData smoother_data;
+  mg_smoother.initialize(mg_matrices, smoother_data);
+  
+                                  // Do two smoothing steps per level
+  mg_smoother.set_steps(2);
+                                  // Since the SOR method is not
+                                  // symmetric, but we use conjugate
+                                  // gradient iteration below, here
+                                  // is a trick to make the
+                                  // multilevel preconditioner a
+                                  // symmetric operator even for
+                                  // nonsymmetric smoothers.
+  mg_smoother.set_symmetric(true);
+
+                                  // We must wrap our matrices in an
+                                  // object having the required
+                                  // multiplication functions.
+  MGMatrix<SparseMatrix<float>, Vector<double> >
+    mg_matrix(&mg_matrices);
+                                  // Now, we are ready to set up the
+                                  // V-cycle operator and the
+                                  // multilevel preconditioner.
+  Multigrid<Vector<double> > mg(mg_dof_handler,
+                               mg_matrix,
+                               mg_coarse,
+                               mg_transfer,
+                               mg_smoother,
+                               mg_smoother);
+  PreconditionMG<dim, Vector<double>, MGTransferPrebuilt<double> >
+    preconditioner(mg_dof_handler, mg, mg_transfer);
+  
+                                  // Finally, create the solver
+                                  // object and solve the system
   SolverControl           solver_control (1000, 1e-12);
-  PrimitiveVectorMemory<> vector_memory;
   SolverCG<>              cg (solver_control, vector_memory);
 
-                                  // The only thing we have to alter
-                                  // is that we need an object which
-                                  // will act as a preconditioner. We
-                                  // will use SSOR (symmetric
-                                  // successive overrelaxation), with
-                                  // a relaxation factor of 1.2. For
-                                  // this purpose, the SparseMatrix
-                                  // class has a function which does
-                                  // one SSOR step, and we need to
-                                  // package the address of this
-                                  // function together with the
-                                  // matrix on which it should act
-                                  // (which is the matrix to be
-                                  // inverted) and the relaxation
-                                  // factor into one object. This can
-                                  // be done like this:
-  PreconditionSSOR<> preconditioner;
-  preconditioner.initialize(system_matrix, 1.2);
-                                  // (Note that we did not have to
-                                  // explicitely pass the address of
-                                  // the SSOR function of the matrix
-                                  // to this objects, rather it is
-                                  // hardcoded into the object, thus
-                                  // the name.)
-                                  //
-                                  // The default template parameters
-                                  // of the ``PreconditionRelaxation''
-                                  // class is the matrix type, which
-                                  // defaults to the types used in
-                                  // this program.
-
-                                  // Calling the solver now looks
-                                  // mostly like in the example
-                                  // before, but where there was an
-                                  // object of type
-                                  // PreconditionIdentity before,
-                                  // there now is the newly generated
-                                  // preconditioner object.
+  
   cg.solve (system_matrix, solution, system_rhs,
            preconditioner);
 
@@ -503,6 +554,7 @@ void LaplaceProblem<dim>::run ()
 
       setup_system ();
       assemble_system ();
+      assemble_multigrid ();
       solve ();
       output_results (cycle);
     };

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.