From: Wolfgang Bangerth Date: Fri, 5 Feb 2010 04:07:25 +0000 (+0000) Subject: Bring up to standard with regard to multigrid on adaptive meshes. X-Git-Tag: v8.0.0~6546 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=173676936a4a7f69e11d80710507ccff63fb0bc6;p=dealii.git Bring up to standard with regard to multigrid on adaptive meshes. git-svn-id: https://svn.dealii.org/trunk@20501 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-16/doc/intro.dox b/deal.II/examples/step-16/doc/intro.dox index 83f0c68a65..e86be0b0b6 100644 --- a/deal.II/examples/step-16/doc/intro.dox +++ b/deal.II/examples/step-16/doc/intro.dox @@ -1,3 +1,12 @@ +
+ +This program has evolved from a version originally written by Guido +Kanschat in 2003. It has undergone significant revisions by Bärbel +Janssen, Guido Kanschat and Wolfgang Bangerth in 2009 and 2010 to demonstrate +multigrid algorithms on adaptively refined meshes. + + +

Introduction

diff --git a/deal.II/examples/step-16/step-16.cc b/deal.II/examples/step-16/step-16.cc index 61f7d1f18d..2b562d310e 100644 --- a/deal.II/examples/step-16/step-16.cc +++ b/deal.II/examples/step-16/step-16.cc @@ -1,9 +1,11 @@ /* $Id$ */ -/* Author: Guido Kanschat, University of Heidelberg, 2003 */ +/* Author: Guido Kanschat, University of Heidelberg, 2003 */ +/* Baerbel Janssen, University of Heidelberg, 2010 */ +/* Wolfgang Bangerth, Texas A&M University, 2010 */ /* $Id$ */ /* */ -/* Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors */ +/* Copyright (C) 2003, 2004, 2006, 2007, 2008, 2009, 2010 by the deal.II authors */ /* */ /* This file is subject to QPL and may not be distributed */ /* without copyright and license information. Please refer */ @@ -16,42 +18,37 @@ #include #include #include +#include + +#include #include #include #include #include #include + #include #include #include -#include #include -#include + #include #include + #include #include + #include -#include #include - // These are the new include files - // required for multi-level methods. - // First, the file defining the - // multigrid method itself. + +//These are the same include files +//as in step-16 necessary for the +//multi-level methods #include - // The DoFHandler is replaced by an - // MGDoFHandler which is defined - // here. #include #include - - // Then, we need some pre-made - // transfer routines between grids. #include - // And a file in which equivalents to the - // DoFTools class are declared: #include - #include #include #include @@ -64,32 +61,25 @@ // previous programs: using namespace dealii; - // This class is based on the same - // class in step-5. Remark that we - // replaced the DoFHandler by - // MGDoFHandler. since this inherits - // from DoFHandler, the new object - // incorporates the old functionality - // plus the new functions for degrees - // of freedom on different - // levels. Furthermore, we added - // MGLevelObjects for sparsity - // patterns and matrices. + +//This class is basically the same +//class as in step-16. The only +//difference is that here we solve Laplace's +//problem on an adaptively refined grid. template class LaplaceProblem { public: - LaplaceProblem (); + LaplaceProblem (const unsigned int deg); void run (); - + private: void setup_system (); void assemble_system (); - // We add this function for - // assembling the multilevel - // matrices. void assemble_multigrid (); void solve (); + void refine_local (); + void test_get_coarse_cell (); void output_results (const unsigned int cycle) const; Triangulation triangulation; @@ -99,32 +89,34 @@ class LaplaceProblem SparsityPattern sparsity_pattern; SparseMatrix system_matrix; - // Here are the new objects for - // handling level matrices: sparsity - // patterns and matrices. We use number - // type float to save memory. It's only - // a preconditioner! - MGLevelObject mg_sparsity; - MGLevelObject > mg_matrices; - + //This object holds the information f + //or the hanging nodes. + ConstraintMatrix constraints; + + MGLevelObject mg_sparsity; + MGLevelObject > mg_matrices; + + /* The matrices at the interface + * between two refinement levels, + * coupling coarse to fine.*/ + MGLevelObject > mg_interface_matrices_up; + Vector solution; Vector system_rhs; -}; + const unsigned int degree; +}; - // This function is as before. template -LaplaceProblem::LaplaceProblem () : - fe (1), - mg_dof_handler (triangulation) +LaplaceProblem::LaplaceProblem (const unsigned int deg) : + triangulation (Triangulation::limit_level_difference_at_vertices), + fe (deg), + mg_dof_handler (triangulation), + degree(deg) {} - - // This is the function of step-5 - // augmented by the setup of the - // multi-level matrices in the end. template void LaplaceProblem::setup_system () { @@ -136,23 +128,33 @@ void LaplaceProblem::setup_system () // multilevel structure deallog << "Number of degrees of freedom: " << mg_dof_handler.n_dofs(); + for (unsigned int l=0;l&>(mg_dof_handler), - sparsity_pattern); - sparsity_pattern.compress(); - - system_matrix.reinit (sparsity_pattern); + DoFTools::make_sparsity_pattern ( + static_cast&>(mg_dof_handler), + sparsity_pattern); solution.reinit (mg_dof_handler.n_dofs()); system_rhs.reinit (mg_dof_handler.n_dofs()); + constraints.clear (); + DoFTools::make_hanging_node_constraints (mg_dof_handler, constraints); + VectorTools::interpolate_boundary_values (mg_dof_handler, + 0, + ZeroFunction(), + constraints); + constraints.close (); + constraints.condense (sparsity_pattern); + sparsity_pattern.compress(); + system_matrix.reinit (sparsity_pattern); + // The multi-level objects are // resized to hold matrices for // every level. The coarse level is @@ -167,7 +169,11 @@ void LaplaceProblem::setup_system () // SparsityPattern before it can be // destroyed. const unsigned int nlevels = triangulation.n_levels(); + + mg_interface_matrices_up.resize(0, nlevels-1); + mg_interface_matrices_up.clear (); mg_matrices.resize(0, nlevels-1); + mg_matrices.clear (); mg_sparsity.resize(0, nlevels-1); // Now, we have to build a matrix @@ -183,8 +189,23 @@ void LaplaceProblem::setup_system () mg_dof_handler.n_dofs(level), mg_dof_handler.max_couplings_between_dofs()); MGTools::make_sparsity_pattern (mg_dof_handler, mg_sparsity[level], level); + CompressedSparsityPattern ci_sparsity; + if(level>0) + { + ci_sparsity.reinit(mg_dof_handler.n_dofs(level), + mg_dof_handler.n_dofs(level)); + MGTools::make_sparsity_pattern(mg_dof_handler, ci_sparsity, level); + } + } + +//And the same for the mg matrices +//for the interface. Note that there +//is no such interface on the coarsest level + for(unsigned int level=0; level::setup_system () template void LaplaceProblem::assemble_system () { - QGauss quadrature_formula(2); + QGauss quadrature_formula(1+degree); FEValues fe_values (fe, quadrature_formula, update_values | update_gradients | - update_quadrature_points | update_JxW_values); + 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(); @@ -208,27 +229,15 @@ void LaplaceProblem::assemble_system () std::vector local_dof_indices (dofs_per_cell); - typename DoFHandler::active_cell_iterator cell = mg_dof_handler.begin_active(), - endc = mg_dof_handler.end(); + typename DoFHandler::active_cell_iterator + cell = mg_dof_handler.begin_active(), + endc = mg_dof_handler.end(); for (; cell!=endc; ++cell) { cell_matrix = 0; cell_rhs = 0; - // As before, we want the FEValues - // object to compute the quantities - // which we told him to compute in - // the constructor using the update - // flags. Then, we loop over all - // quadrature points and the local - // matrix rows and columns for - // computing the element - // contribution. This is the same as - // in step-4. For the right hand - // side, we use a constant value of - // 1. fe_values.reinit (cell); - for (unsigned int q_point=0; q_point::assemble_system () cell_rhs(i) += (fe_values.shape_value(i,q_point) * 1.0 * fe_values.JxW(q_point)); - }; + } cell->get_dof_indices (local_dof_indices); - for (unsigned int i=0; i boundary_values; - - VectorTools::interpolate_boundary_values (mg_dof_handler, - 0, - ZeroFunction(), - boundary_values); - - MatrixTools::apply_boundary_values (boundary_values, - system_matrix, - solution, - system_rhs); + constraints.distribute_local_to_global (cell_matrix, cell_rhs, + local_dof_indices, + system_matrix, system_rhs); + } } @@ -275,7 +263,10 @@ 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 correct matrix. + // into the correct matrix. Here comes + // the difference to global refinement + // into play. We have to fill the interface + // matrices correctly. // Since we only do multi-level // preconditioning, no right-hand @@ -283,26 +274,57 @@ void LaplaceProblem::assemble_system () template void LaplaceProblem::assemble_multigrid () { - QGauss quadrature_formula(2); + QGauss quadrature_formula(1+degree); FEValues fe_values (fe, quadrature_formula, update_values | update_gradients | - update_quadrature_points | update_JxW_values); + 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(); FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); std::vector local_dof_indices (dofs_per_cell); + std::vector > interface_dofs; + std::vector > boundary_interface_dofs; + for (unsigned int level = 0; level tmp (mg_dof_handler.n_dofs(level)); + interface_dofs.push_back (tmp); + boundary_interface_dofs.push_back (tmp); + } + MGTools::extract_inner_interface_dofs (mg_dof_handler, + interface_dofs, + boundary_interface_dofs); + + typename FunctionMap::type dirichlet_boundary; + ZeroFunction homogeneous_dirichlet_bc (1); + dirichlet_boundary[0] = &homogeneous_dirichlet_bc; + + std::vector > boundary_indices(triangulation.n_levels()); + MGTools::make_boundary_list (mg_dof_handler, dirichlet_boundary, + boundary_indices); + + std::vector boundary_constraints (triangulation.n_levels()); + std::vector boundary_interface_constraints (triangulation.n_levels()); + for (unsigned int level=0; level::cell_iterator cell = mg_dof_handler.begin(), endc = mg_dof_handler.end(); + for (; cell!=endc; ++cell) { - // Remember the level of the - // current cell. - const unsigned int level = cell->level(); cell_matrix = 0; // Compute the values specified @@ -321,7 +343,6 @@ void LaplaceProblem::assemble_multigrid () * fe_values.JxW(q_point); } - // Oops! This is a tiny // difference easily // forgotten. The indices we @@ -332,66 +353,26 @@ void LaplaceProblem::assemble_multigrid () // 'mg' entered into the // function call. cell->get_mg_dof_indices (local_dof_indices); + + const unsigned int level = cell->level(); + boundary_constraints[level] + .distribute_local_to_global (cell_matrix, + local_dof_indices, + mg_matrices[level]); + for (unsigned int i=0; i::type dirichlet_boundary; - ZeroFunction homogeneous_dirichlet_bc (1); - dirichlet_boundary[0] = &homogeneous_dirichlet_bc; - - // Next we generate the set of - // dof indices to be eliminated on - // each level. - std::vector > boundary_indices(triangulation.n_levels()); - MGTools::make_boundary_list (mg_dof_handler, - dirichlet_boundary, - boundary_indices); - - // And finally we eliminate these - // degrees of freedom from every - // matrix in the multilevel hierarchy. - const unsigned int nlevels = triangulation.n_levels(); - for (unsigned int level=0;level void LaplaceProblem::solve () { @@ -406,7 +387,7 @@ void LaplaceProblem::solve () // the transfer of functions // between different grid // levels. - MGTransferPrebuilt > mg_transfer; + MGTransferPrebuilt > mg_transfer(constraints); mg_transfer.build_matrices(mg_dof_handler); // Next, we need a coarse grid @@ -415,9 +396,9 @@ void LaplaceProblem::solve () // direct solver, even if its // implementation is not very // clever. - FullMatrix coarse_matrix; + FullMatrix coarse_matrix; coarse_matrix.copy_from (mg_matrices[0]); - MGCoarseGridHouseholder > mg_coarse; + MGCoarseGridHouseholder > mg_coarse; mg_coarse.initialize(coarse_matrix); // The final ingredient for the @@ -427,20 +408,22 @@ void LaplaceProblem::solve () // here. Names are getting quite // long here, so we help with // typedefs. - typedef PreconditionSOR > RELAXATION; - MGSmootherRelaxation, RELAXATION, Vector > + typedef PreconditionSOR > RELAXATION; +// typedef PreconditionJacobi > RELAXATION; +// typedef SparseILU RELAXATION; + MGSmootherRelaxation, RELAXATION, Vector > mg_smoother(vector_memory); // Initialize the smoother with our - // level matrices and the (required) + // level matrices and the required, // additional data for the - // relaxation method with default + // relaxaton method with default // values. - RELAXATION::AdditionalData smoother_data; + RELAXATION::AdditionalData smoother_data;//(0, 9,false); mg_smoother.initialize(mg_matrices, smoother_data); - + // Do two smoothing steps per level - mg_smoother.set_steps(2); + mg_smoother.set_steps(1); // Since the SOR method is not // symmetric, but we use conjugate // gradient iteration below, here @@ -449,12 +432,19 @@ void LaplaceProblem::solve () // symmetric operator even for // nonsymmetric smoothers. mg_smoother.set_symmetric(true); + mg_smoother.set_variable(true); + // We must wrap our matrices in an // object having the required // multiplication functions. - MGMatrix, Vector > + MGMatrix, Vector > mg_matrix(&mg_matrices); + //do the same for the interface matrices + MGMatrix, Vector > + mg_interface_up(&mg_interface_matrices_up); + MGMatrix, Vector > + mg_interface_down(&mg_interface_matrices_up); // Now, we are ready to set up the // V-cycle operator and the // multilevel preconditioner. @@ -464,66 +454,177 @@ void LaplaceProblem::solve () mg_transfer, mg_smoother, mg_smoother); + mg.set_edge_matrices(mg_interface_down, mg_interface_up); + PreconditionMG, MGTransferPrebuilt > > - preconditioner(mg_dof_handler, mg, mg_transfer); - - // Finally, create the solver - // object and solve the system - SolverControl solver_control (1000, 1e-12); - SolverCG<> cg (solver_control); + preconditioner(mg_dof_handler, mg, mg_transfer); - - cg.solve (system_matrix, solution, system_rhs, - preconditioner); + // Finally, create the solver + // object and solve the system +ReductionControl solver_control (100, 1.e-20, 1.e-10, true, true); +SolverCG<> cg (solver_control); + +solution = 0; + +cg.solve (system_matrix, solution, system_rhs, + preconditioner); +constraints.distribute (solution); + +std::cout << " " << solver_control.last_step() +<< " CG iterations needed to obtain convergence." +<< std::endl; } +template +void LaplaceProblem::refine_local () +{ + bool cell_refined = false; + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + { + if(false) + { + for (unsigned int vertex=0; + vertex < GeometryInfo::vertices_per_cell; + ++vertex) + { + if(cell->vertex(vertex)[dim-1]==0 && cell->vertex(vertex)[0]==0 && cell->vertex(vertex)[dim-2]==0) + { + cell->set_refine_flag (); + cell_refined = true; + break; + } + } + } + else if(true) //Kreis + { + for (unsigned int vertex=0; + vertex < GeometryInfo::vertices_per_cell; + ++vertex) + { + const Point p = cell->vertex(vertex); + const Point origin = (dim == 2 ? + Point(0,0) : + Point(0,0,0)); + const double dist = p.distance(origin); + if(dist<0.25/M_PI) + { + cell->set_refine_flag (); + cell_refined = true; + break; + } + } + } + else if(false) //linke Diagonale + { + for (unsigned int vertex=0; + vertex < GeometryInfo::vertices_per_cell; + ++vertex) + { + const Point p = cell->vertex(vertex); + if(p[0]==p[1]) + { + cell->set_refine_flag (); + cell_refined = true; + break; + } + } + } + else if(false) //inneres Quadrat + { + for (unsigned int vertex=0; + vertex < GeometryInfo::vertices_per_cell; + ++vertex) + { + const Point p = cell->vertex(vertex); + const double dist = std::max(std::fabs(p[0]),std::fabs(p[1])); + if(dist<0.5) + { + cell->set_refine_flag (); + cell_refined = true; + break; + } + } + } + else if(false) //Raute + { + for (unsigned int vertex=0; + vertex < GeometryInfo::vertices_per_cell; + ++vertex) + { + const Point p = cell->vertex(vertex); + const double dist = std::fabs(p[0])+std::fabs(p[1]); + if(dist<0.5) + { + cell->set_refine_flag (); + cell_refined = true; + break; + } + } + } + else //erster Quadrant + { + const Point p = cell->center(); + bool positive = p(0) > 0; + if (dim>1 && p(1) <= 0) + positive = false; + if (dim>2 && p(2) <= 0) + positive = false; + if (positive) + { + cell->set_refine_flag(); + cell_refined = true; + } + } + } + //Wenn nichts verfeinert wurde bisher, global verfeinern! + if(!cell_refined) + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + cell->set_refine_flag(); + + triangulation.execute_coarsening_and_refinement (); +} - // Here is the data output, which is - // a simplified version of step-5. We - // do a standard gnuplot output for - // each grid produced in the - // refinement process. template void LaplaceProblem::output_results (const unsigned int cycle) const { - // Construct and initialize a DataOut object DataOut data_out; data_out.attach_dof_handler (mg_dof_handler); data_out.add_data_vector (solution, "solution"); data_out.build_patches (); - // The following block generates - // the file name and opens the - // file: std::ostringstream filename; filename << "solution-" << cycle - << ".gnuplot"; + << ".vtk"; std::ofstream output (filename.str().c_str()); - data_out.write_gnuplot (output); + data_out.write_vtk (output); } - - template void LaplaceProblem::run () { - for (unsigned int cycle=0; cycle<6; ++cycle) + for (unsigned int cycle=0; cycle<9; ++cycle) { - deallog << "Cycle " << cycle << std::endl; - if (cycle == 0) { - // Generate a simple hyperball grid. - GridGenerator::hyper_ball(triangulation); - static const HyperBallBoundary boundary; - triangulation.set_boundary (0, boundary); + GridGenerator::hyper_cube(triangulation, -1, 1); + triangulation.refine_global (1); } - triangulation.refine_global (1); + refine_local (); + + std::cout << "Cycle " << cycle + << " with " << triangulation.n_active_cells() + << " cells." + << std::endl; + setup_system (); assemble_system (); assemble_multigrid (); @@ -532,16 +633,41 @@ void LaplaceProblem::run () }; } - - // The main function looks mostly - // like the one in the previous - // example, so we won't comment on it - // further. + int main () { - LaplaceProblem<2> laplace_problem_2d; - laplace_problem_2d.run (); - + try + { + deallog.depth_console (0); + + LaplaceProblem<2> laplace_problem(1); + laplace_problem.run (); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + return 0; }