From c346d7718c892a589612ac22c814658bc1912b10 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 9 Feb 2010 23:59:12 +0000 Subject: [PATCH] Bring closer to what we do in step-6. git-svn-id: https://svn.dealii.org/trunk@20531 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-16/step-16.cc | 175 ++++++++-------------------- 1 file changed, 46 insertions(+), 129 deletions(-) diff --git a/deal.II/examples/step-16/step-16.cc b/deal.II/examples/step-16/step-16.cc index 13da73cfc8..93e2c348a9 100644 --- a/deal.II/examples/step-16/step-16.cc +++ b/deal.II/examples/step-16/step-16.cc @@ -31,6 +31,8 @@ #include #include #include +#include +#include #include #include @@ -40,6 +42,7 @@ #include #include +#include //These are the same include files //as in step-16 necessary for the @@ -78,8 +81,7 @@ class LaplaceProblem void assemble_system (); void assemble_multigrid (); void solve (); - void refine_local (); - void test_get_coarse_cell (); + void refine_grid (); void output_results (const unsigned int cycle) const; Triangulation triangulation; @@ -461,133 +463,35 @@ void LaplaceProblem::solve () MGTransferPrebuilt > > preconditioner(mg_dof_handler, mg, mg_transfer); - // 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); + // 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; + solution = 0; -cg.solve (system_matrix, solution, system_rhs, - preconditioner); -constraints.distribute (solution); + 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; + std::cout << " " << solver_control.last_step() + << " CG iterations needed to obtain convergence." + << std::endl; } template -void LaplaceProblem::refine_local () +void LaplaceProblem::refine_grid () { - 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(); - - + Vector estimated_error_per_cell (triangulation.n_active_cells()); + + KellyErrorEstimator::estimate (static_cast&>(mg_dof_handler), + QGauss(3), + typename FunctionMap::type(), + solution, + estimated_error_per_cell); + GridRefinement::refine_and_coarsen_fixed_number (triangulation, + estimated_error_per_cell, + 0.3, 0.03); triangulation.execute_coarsening_and_refinement (); } @@ -612,26 +516,39 @@ void LaplaceProblem::output_results (const unsigned int cycle) const template void LaplaceProblem::run () { - for (unsigned int cycle=0; cycle<9; ++cycle) + for (unsigned int cycle=0; cycle<8; ++cycle) { + std::cout << "Cycle " << cycle << ':' << std::endl; + if (cycle == 0) { - GridGenerator::hyper_cube(triangulation, -1, 1); + GridGenerator::hyper_ball (triangulation); + + static const HyperBallBoundary boundary; + triangulation.set_boundary (0, boundary); + triangulation.refine_global (1); } - refine_local (); + else + refine_grid (); - std::cout << "Cycle " << cycle - << " with " << triangulation.n_active_cells() - << " cells." + + std::cout << " Number of active cells: " + << triangulation.n_active_cells() << std::endl; setup_system (); + + std::cout << " Number of degrees of freedom: " + << mg_dof_handler.n_dofs() + << std::endl; + assemble_system (); assemble_multigrid (); + solve (); output_results (cycle); - }; + } } -- 2.39.5