From: guido Date: Sun, 1 Feb 2004 22:27:46 +0000 (+0000) Subject: step 16 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a9b6118c803e296a396e7f1cb182690ffe45092e;p=dealii-svn.git step 16 git-svn-id: https://svn.dealii.org/trunk@8379 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-16.data/intro.html b/deal.II/doc/tutorial/chapter-2.step-by-step/step-16.data/intro.html index c59cfdf609..39d27081c8 100644 --- a/deal.II/doc/tutorial/chapter-2.step-by-step/step-16.data/intro.html +++ b/deal.II/doc/tutorial/chapter-2.step-by-step/step-16.data/intro.html @@ -1,2 +1,34 @@

Introduction

+ +

+This example shows the basic usage of the multilevel functions in +deal.II. 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. +

+ +

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 +

+

+ +

+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. +

+ +

+The multilevel method in deal.II 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. +

\ No newline at end of file diff --git a/deal.II/examples/step-16/Makefile b/deal.II/examples/step-16/Makefile index 433c82b233..dae158698b 100644 --- a/deal.II/examples/step-16/Makefile +++ b/deal.II/examples/step-16/Makefile @@ -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. # diff --git a/deal.II/examples/step-16/step-16.cc b/deal.II/examples/step-16/step-16.cc index 1c6a1d1fd1..73f3c5f51b 100644 --- a/deal.II/examples/step-16/step-16.cc +++ b/deal.II/examples/step-16/step-16.cc @@ -53,6 +53,10 @@ // DoFTools class are declared: #include +#include +#include +#include + // This is C++ ... see step 5 for // further comments. #include @@ -151,8 +155,14 @@ void LaplaceProblem::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::assemble_system () for (unsigned int i=0; i::assemble_system () }; // Again use zero boundary values: - std::map boundary_values; - VectorTools::interpolate_boundary_values (mg_dof_handler, - 0, - ZeroFunction(), - boundary_values); - MatrixTools::apply_boundary_values (boundary_values, - system_matrix, - solution, - system_rhs); +// std::map boundary_values; +// VectorTools::interpolate_boundary_values (mg_dof_handler, +// 0, +// ZeroFunction(), +// boundary_values); +// MatrixTools::apply_boundary_values (boundary_values, +// system_matrix, +// solution, +// system_rhs); } @@ -280,7 +292,7 @@ void LaplaceProblem::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::assemble_multigrid () QGauss2 quadrature_formula; FEValues 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::assemble_multigrid () std::vector local_dof_indices (dofs_per_cell); // - typename DoFHandler::cell_iterator cell = mg_dof_handler.begin(), - endc = mg_dof_handler.end(); + typename MGDoFHandler::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::assemble_multigrid () for (unsigned int i=0; i::assemble_multigrid () }; // Again use zero boundary values: - std::map boundary_values; - VectorTools::interpolate_boundary_values (mg_dof_handler, - 0, - ZeroFunction(), - boundary_values); - MatrixTools::apply_boundary_values (boundary_values, - system_matrix, - solution, - system_rhs); +// std::map boundary_values; +// VectorTools::interpolate_boundary_values (mg_dof_handler, +// 0, +// ZeroFunction(), +// boundary_values); +// MatrixTools::apply_boundary_values (boundary_values, +// system_matrix, +// solution, +// system_rhs); } @@ -374,48 +389,84 @@ void LaplaceProblem::assemble_multigrid () template void LaplaceProblem::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 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 coarse_matrix; + coarse_matrix.copy_from (mg_matrices[0]); + MGCoarseGridHouseholder > 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 > RELAXATION; + MGSmootherRelaxation, RELAXATION, Vector > + 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, Vector > + mg_matrix(&mg_matrices); + // Now, we are ready to set up the + // V-cycle operator and the + // multilevel preconditioner. + Multigrid > mg(mg_dof_handler, + mg_matrix, + mg_coarse, + mg_transfer, + mg_smoother, + mg_smoother); + PreconditionMG, MGTransferPrebuilt > + 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::run () setup_system (); assemble_system (); + assemble_multigrid (); solve (); output_results (cycle); };