From aefd8a82b720431ccc36d62c71032e2a2a74b378 Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Thu, 31 Jul 2014 13:23:16 +0200 Subject: [PATCH] Fixing problems with MGConstrainedDoFs --- .../deal.II/multigrid/mg_constrained_dofs.h | 105 ++- include/deal.II/multigrid/mg_tools.h | 29 +- source/multigrid/mg_tools.cc | 98 +-- source/multigrid/mg_tools.inst.in | 5 - tests/multigrid/step-16-02.cc | 6 +- tests/multigrid/step-16-bdry1.cc | 637 ++++++++++++++++++ tests/multigrid/step-16-bdry1.output | 57 ++ tests/multigrid/step-16.cc | 6 +- 8 files changed, 745 insertions(+), 198 deletions(-) create mode 100644 tests/multigrid/step-16-bdry1.cc create mode 100644 tests/multigrid/step-16-bdry1.output diff --git a/include/deal.II/multigrid/mg_constrained_dofs.h b/include/deal.II/multigrid/mg_constrained_dofs.h index 0547380d32..d7cac5aa11 100644 --- a/include/deal.II/multigrid/mg_constrained_dofs.h +++ b/include/deal.II/multigrid/mg_constrained_dofs.h @@ -42,18 +42,6 @@ class MGConstrainedDoFs : public Subscriptor public: typedef std::vector >::size_type size_dof; - /** - * Fill the internal data structures with values extracted from the dof - * handler object. - * - * This function ensures that on every level, degrees of freedom at interior - * edges of a refinement level are treated corrected but leaves degrees of - * freedom at the boundary of the domain untouched assuming that no - * Dirichlet boundary conditions for them exist. - */ - template - void initialize (const DoFHandler &dof); - /** * Fill the internal data structures with values extracted from the dof * handler object and apply the boundary values provided. @@ -74,18 +62,17 @@ public: /** * Determine whether a dof index is subject to a boundary - * constraint. (In other words, whether it is on boundary of domain.) + * constraint. */ bool is_boundary_index (const unsigned int level, const types::global_dof_index index) const; /** - * Determine whether a dof index is at an edge that is not a + * @deprecated Determine whether a dof index is at an edge that is not a * refinement edge. */ -// TODO: remove bool non_refinement_edge_index (const unsigned int level, - const types::global_dof_index index) const; + const types::global_dof_index index) const DEAL_II_DEPRECATED; /** * Determine whether a dof index is at the refinement edge. @@ -94,16 +81,19 @@ public: const types::global_dof_index index) const; /** - * Determine whether a dof index is at the refinement edge and - * subject to a boundary constraint . - * = is_boundary_index() && at_refinement_edge() + * @deprecated Use is_boundary_index() instead. The logic behind + * this function here is unclear and for practical purposes, the + * other is needed. + * + * Determine whether a dof index is subject to a boundary + * constraint. */ bool at_refinement_edge_boundary (const unsigned int level, - const types::global_dof_index index) const; + const types::global_dof_index index) const DEAL_II_DEPRECATED; /** - * Return the indices of dofs for each level that lie on the - * boundary of the domain. + * Return the indices of dofs for each level that are subject to + * boundary constraints. */ const std::vector > & get_boundary_indices () const; @@ -137,13 +127,6 @@ public: const std::vector > & get_refinement_edge_boundary_indices () const DEAL_II_DEPRECATED; - /** - * Return indices of all dofs that are on boundary faces on the given level - * if the cell has refinement edge indices (i.e. has a coarser neighbor). - */ - const IndexSet & - get_refinement_edge_boundary_indices (unsigned int level) const; - /** * @deprecated The function is_boundary_index() now returns false if * no boundary values are set. @@ -154,6 +137,21 @@ public: private: + /** + * @warning This function generates an inconsistent object if not + * called from the other initialize() in this class. + * + * Fill the internal data structures with values extracted from the dof + * handler object. + * + * This function ensures that on every level, degrees of freedom at interior + * edges of a refinement level are treated corrected but leaves degrees of + * freedom at the boundary of the domain untouched assuming that no + * Dirichlet boundary conditions for them exist. + */ + template + void initialize (const DoFHandler &dof); + /** * The indices of boundary dofs for each level. */ @@ -165,14 +163,6 @@ private: */ std::vector refinement_edge_indices; - /** - * The degrees of freedom on the refinement edge between a level and - * coarser cells, which are also on the boundary. - * - * This is a subset of #refinement_edge_indices. - */ - std::vector refinement_edge_boundary_indices; - /** * old data structure only filled on demand */ @@ -195,7 +185,6 @@ MGConstrainedDoFs::initialize(const DoFHandler &dof) boundary_indices.resize(nlevels); refinement_edge_indices.resize(nlevels); - refinement_edge_boundary_indices.resize(nlevels); refinement_edge_indices_old.clear(); refinement_edge_boundary_indices_old.clear(); for (unsigned int l=0; l &dof) boundary_indices[l].clear(); refinement_edge_indices[l] = IndexSet(dof.n_dofs(l)); - refinement_edge_boundary_indices[l] = IndexSet(dof.n_dofs(l)); } - MGTools::extract_inner_interface_dofs (dof, refinement_edge_indices, - refinement_edge_boundary_indices); + MGTools::extract_inner_interface_dofs (dof, refinement_edge_indices); } @@ -237,9 +224,6 @@ MGConstrainedDoFs::clear() for (unsigned int l=0; l > & MGConstrainedDoFs::get_refinement_edge_boundary_indices () const { - if (refinement_edge_boundary_indices_old.size()!=refinement_edge_boundary_indices.size()) - { - unsigned int n_levels = refinement_edge_boundary_indices.size(); - refinement_edge_boundary_indices_old.resize(n_levels); - for (unsigned int l=0;l - void - extract_inner_interface_dofs (const DoFHandler &mg_dof_handler, - std::vector &interface_dofs, - std::vector &boundary_interface_dofs); - - /** - * Does the same as the function above, - * but fills only the interface_dofs. + * For each level in a multigrid hierarchy, produce an IndexSet + * that indicates which of the degrees of freedom are along + * interfaces of this level to cells that only exist on coarser + * levels. */ template void diff --git a/source/multigrid/mg_tools.cc b/source/multigrid/mg_tools.cc index 9117d49916..4f942a4663 100644 --- a/source/multigrid/mg_tools.cc +++ b/source/multigrid/mg_tools.cc @@ -1418,72 +1418,7 @@ namespace MGTools } } - - - - template - void - extract_inner_interface_dofs (const DoFHandler &mg_dof_handler, - std::vector &interface_dofs) - { - Assert (interface_dofs.size() == mg_dof_handler.get_tria().n_global_levels(), - ExcDimensionMismatch (interface_dofs.size(), - mg_dof_handler.get_tria().n_global_levels())); - - for (unsigned int l=0; l &fe = mg_dof_handler.get_fe(); - - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int dofs_per_face = fe.dofs_per_face; - - std::vector local_dof_indices (dofs_per_cell); - std::vector cell_dofs(dofs_per_cell, false); - - typename DoFHandler::cell_iterator cell = mg_dof_handler.begin(), - endc = mg_dof_handler.end(); - - for (; cell!=endc; ++cell) - { - if (mg_dof_handler.get_tria().locally_owned_subdomain()!=numbers::invalid_subdomain_id - && cell->level_subdomain_id()!=mg_dof_handler.get_tria().locally_owned_subdomain()) - continue; - - std::fill (cell_dofs.begin(), cell_dofs.end(), false); - - for (unsigned int face_nr=0; face_nr::faces_per_cell; ++face_nr) - { - const typename DoFHandler::face_iterator face = cell->face(face_nr); - if (!face->at_boundary()) - { - //interior face - const typename DoFHandler::cell_iterator - neighbor = cell->neighbor(face_nr); - - if (neighbor->level() < cell->level()) - { - for (unsigned int j=0; jlevel(); - cell->get_mg_dof_indices (local_dof_indices); - - for (unsigned int i=0; i void extract_inner_interface_dofs (const DoFHandler &mg_dof_handler, @@ -1579,20 +1514,14 @@ namespace MGTools template void extract_inner_interface_dofs (const DoFHandler &mg_dof_handler, - std::vector &interface_dofs, - std::vector &boundary_interface_dofs) + std::vector &interface_dofs) { Assert (interface_dofs.size() == mg_dof_handler.get_tria().n_global_levels(), ExcDimensionMismatch (interface_dofs.size(), mg_dof_handler.get_tria().n_global_levels())); - Assert (boundary_interface_dofs.size() == mg_dof_handler.get_tria().n_global_levels(), - ExcDimensionMismatch (boundary_interface_dofs.size(), - mg_dof_handler.get_tria().n_global_levels())); std::vector > tmp_interface_dofs(interface_dofs.size()); - std::vector > - tmp_boundary_interface_dofs(interface_dofs.size()); const FiniteElement &fe = mg_dof_handler.get_fe(); @@ -1600,10 +1529,8 @@ namespace MGTools const unsigned int dofs_per_face = fe.dofs_per_face; std::vector local_dof_indices (dofs_per_cell); - std::vector face_dof_indices (dofs_per_face); std::vector cell_dofs(dofs_per_cell, false); - std::vector boundary_cell_dofs(dofs_per_cell, false); typename DoFHandler::cell_iterator cell = mg_dof_handler.begin(), endc = mg_dof_handler.end(); @@ -1617,7 +1544,6 @@ namespace MGTools bool has_coarser_neighbor = false; std::fill (cell_dofs.begin(), cell_dofs.end(), false); - std::fill (boundary_cell_dofs.begin(), boundary_cell_dofs.end(), false); for (unsigned int face_nr=0; face_nr::faces_per_cell; ++face_nr) { @@ -1642,14 +1568,7 @@ namespace MGTools if (has_coarser_neighbor == false) continue; - - for (unsigned int face_nr=0; face_nr::faces_per_cell; ++face_nr) - if (cell->at_boundary(face_nr)) - for (unsigned int j=0; jlevel(); cell->get_mg_dof_indices (local_dof_indices); @@ -1657,13 +1576,9 @@ namespace MGTools { if (cell_dofs[i]) tmp_interface_dofs[level].push_back(local_dof_indices[i]); - - if (boundary_cell_dofs[i]) - tmp_boundary_interface_dofs[level].push_back(local_dof_indices[i]); } } - for (unsigned int l=0; l&, const ComponentMask &); - template - void - extract_inner_interface_dofs (const DoFHandler &mg_dof_handler, - std::vector &interface_dofs, - std::vector &boundary_interface_dofs); template void extract_inner_interface_dofs (const DoFHandler &mg_dof_handler, diff --git a/tests/multigrid/step-16-02.cc b/tests/multigrid/step-16-02.cc index 60e181bcf9..cd4fc73163 100644 --- a/tests/multigrid/step-16-02.cc +++ b/tests/multigrid/step-16-02.cc @@ -14,11 +14,7 @@ // // --------------------------------------------------------------------- -/* - * Authors: Guido Kanschat, University of Heidelberg, 2003 - * Baerbel Janssen, University of Heidelberg, 2010 - * Wolfgang Bangerth, Texas A&M University, 2010 - */ +// Multigrid for continuous finite elements using MeshWorker #include "../tests.h" #include diff --git a/tests/multigrid/step-16-bdry1.cc b/tests/multigrid/step-16-bdry1.cc new file mode 100644 index 0000000000..95359d1765 --- /dev/null +++ b/tests/multigrid/step-16-bdry1.cc @@ -0,0 +1,637 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2003 - 2013 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// Multigrid for continuous finite elements using MeshWorker +// Investigate a seeming inconsistency in MGConstrainedDoFs at the boundary + +#include "../tests.h" +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +using namespace dealii; +using namespace LocalIntegrators; + +template +class LaplaceMatrix : public MeshWorker::LocalIntegrator +{ +public: + LaplaceMatrix(); + virtual void cell(MeshWorker::DoFInfo& dinfo, MeshWorker::IntegrationInfo& info) const; + virtual void boundary(MeshWorker::DoFInfo& dinfo, MeshWorker::IntegrationInfo& info) const; + virtual void face(MeshWorker::DoFInfo& dinfo1, MeshWorker::DoFInfo& dinfo2, + MeshWorker::IntegrationInfo& info1, MeshWorker::IntegrationInfo& info2) const; +}; + + +template +LaplaceMatrix::LaplaceMatrix() + : + MeshWorker::LocalIntegrator(true, false, false) +{} + + +template +void LaplaceMatrix::cell(MeshWorker::DoFInfo& dinfo, MeshWorker::IntegrationInfo& info) const +{ + AssertDimension (dinfo.n_matrices(), 1); + Laplace::cell_matrix(dinfo.matrix(0,false).matrix, info.fe_values(0)); +} + + +template +void LaplaceMatrix::boundary(MeshWorker::DoFInfo& /*dinfo*/, + typename MeshWorker::IntegrationInfo& /*info*/) const +{ +// const unsigned int deg = info.fe_values(0).get_fe().tensor_degree(); +// Laplace::nitsche_matrix(dinfo.matrix(0,false).matrix, info.fe_values(0), +// Laplace::compute_penalty(dinfo, dinfo, deg, deg)); +} + + +template +void LaplaceMatrix::face( + MeshWorker::DoFInfo& /*dinfo1*/, MeshWorker::DoFInfo& /*dinfo2*/, + MeshWorker::IntegrationInfo& /*info1*/, MeshWorker::IntegrationInfo& /*info2*/) const +{ +// const unsigned int deg = info1.fe_values(0).get_fe().tensor_degree(); +// Laplace::ip_matrix(dinfo1.matrix(0,false).matrix, dinfo1.matrix(0,true).matrix, +// dinfo2.matrix(0,true).matrix, dinfo2.matrix(0,false).matrix, +// info1.fe_values(0), info2.fe_values(0), +// Laplace::compute_penalty(dinfo1, dinfo2, deg, deg)); +} + +template +class LaplaceProblem +{ +public: + LaplaceProblem (const unsigned int deg); + void run (); + +private: + void setup_system (); + void assemble_system (); + void assemble_multigrid (const bool &use_mw); + void solve (); + void refine_grid (const std::string& reftype); + void output_results (const unsigned int cycle) const; + + Triangulation triangulation; + FE_Q fe; + MGDoFHandler mg_dof_handler; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + ConstraintMatrix hanging_node_constraints; + ConstraintMatrix constraints; + + Vector solution; + Vector system_rhs; + + const unsigned int degree; + LaplaceMatrix matrix_integrator; + + MGLevelObject mg_sparsity_patterns; + MGLevelObject > mg_matrices; + MGLevelObject > mg_interface_in; + MGLevelObject > mg_interface_out; + MGConstrainedDoFs mg_constrained_dofs; +}; + + +template +class Coefficient : public Function +{ +public: + Coefficient () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual void value_list (const std::vector > &points, + std::vector &values, + const unsigned int component = 0) const; +}; + + + +template +double Coefficient::value (const Point &p, + const unsigned int) const +{ +// if (p.square() < 0.5*0.5) +// return 20; +// else + return 1; +} + + + +template +void Coefficient::value_list (const std::vector > &points, + std::vector &values, + const unsigned int component) const +{ + const unsigned int n_points = points.size(); + + Assert (values.size() == n_points, + ExcDimensionMismatch (values.size(), n_points)); + + Assert (component == 0, + ExcIndexRange (component, 0, 1)); + + for (unsigned int i=0; i::value (points[i]); +} + + +template +LaplaceProblem::LaplaceProblem (const unsigned int degree) + : + triangulation (Triangulation:: + limit_level_difference_at_vertices), + fe (degree), + mg_dof_handler (triangulation), + degree(degree), + matrix_integrator() +{} + + +template +void LaplaceProblem::setup_system () +{ + mg_dof_handler.distribute_dofs (fe); + deallog << "Number of degrees of freedom: " + << mg_dof_handler.n_dofs(); + + for (unsigned int l=0; l&>(mg_dof_handler), + sparsity_pattern); + + solution.reinit (mg_dof_handler.n_dofs()); + system_rhs.reinit (mg_dof_handler.n_dofs()); + + constraints.clear (); + hanging_node_constraints.clear (); + DoFTools::make_hanging_node_constraints (mg_dof_handler, constraints); + DoFTools::make_hanging_node_constraints (mg_dof_handler, hanging_node_constraints); + typename FunctionMap::type dirichlet_boundary; + ZeroFunction homogeneous_dirichlet_bc (1); + dirichlet_boundary[0] = &homogeneous_dirichlet_bc; + MappingQ1 mapping; + VectorTools::interpolate_boundary_values (mapping, + mg_dof_handler, + dirichlet_boundary, + constraints); + constraints.close (); + hanging_node_constraints.close (); + constraints.condense (sparsity_pattern); + sparsity_pattern.compress(); + system_matrix.reinit (sparsity_pattern); + + mg_constrained_dofs.clear(); + mg_constrained_dofs.initialize(mg_dof_handler, dirichlet_boundary); + const unsigned int n_levels = triangulation.n_levels(); + + mg_interface_in.resize(0, n_levels-1); + mg_interface_in.clear (); + mg_interface_out.resize(0, n_levels-1); + mg_interface_out.clear (); + mg_matrices.resize(0, n_levels-1); + mg_matrices.clear (); + mg_sparsity_patterns.resize(0, n_levels-1); + + for (unsigned int level=0; level +void LaplaceProblem::assemble_system () +{ + const QGauss quadrature_formula(degree+1); + + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + 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(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + const Coefficient coefficient; + std::vector coefficient_values (n_q_points); + + typename MGDoFHandler::active_cell_iterator + cell = mg_dof_handler.begin_active(), + endc = mg_dof_handler.end(); + for (; cell!=endc; ++cell) + { + cell_matrix = 0; + cell_rhs = 0; + + fe_values.reinit (cell); + + coefficient.value_list (fe_values.get_quadrature_points(), + coefficient_values); + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global (cell_matrix, cell_rhs, + local_dof_indices, + system_matrix, system_rhs); + } +} + + +template +void LaplaceProblem::assemble_multigrid (const bool& use_mw) +{ + if(use_mw == true) + { + mg_matrices = 0.; + + MappingQ1 mapping; + MeshWorker::IntegrationInfoBox info_box; + UpdateFlags update_flags = update_values | update_gradients | update_hessians; + info_box.add_update_flags_all(update_flags); + info_box.initialize(fe, mapping); + + MeshWorker::DoFInfo dof_info(mg_dof_handler); + + MeshWorker::Assembler::MGMatrixSimple > assembler; + assembler.initialize(mg_constrained_dofs); + assembler.initialize(mg_matrices); + assembler.initialize_interfaces(mg_interface_in, mg_interface_out); + + MeshWorker::integration_loop ( + mg_dof_handler.begin(), mg_dof_handler.end(), + dof_info, info_box, matrix_integrator, assembler); + + const unsigned int nlevels = triangulation.n_levels(); + for (unsigned int level=0;level quadrature_formula(1+degree); + + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + 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(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + const Coefficient coefficient; + std::vector coefficient_values (n_q_points); + + std::vector > interface_dofs + = mg_constrained_dofs.get_refinement_edge_indices (); + std::vector > boundary_interface_dofs + = mg_constrained_dofs.get_refinement_edge_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) + { + cell_matrix = 0; + fe_values.reinit (cell); + + coefficient.value_list (fe_values.get_quadrature_points(), + coefficient_values); + + for (unsigned int q_point=0; q_pointget_mg_dof_indices (local_dof_indices); + + boundary_constraints[cell->level()] + .distribute_local_to_global (cell_matrix, + local_dof_indices, + mg_matrices[cell->level()]); + + for (unsigned int i=0; ilevel()][local_dof_indices[i]]==true && + interface_dofs[cell->level()][local_dof_indices[j]]==false)) + cell_matrix(i,j) = 0; + + boundary_interface_constraints[cell->level()] + .distribute_local_to_global (cell_matrix, + local_dof_indices, + mg_interface_in[cell->level()]); + } + } +} + + +template +void LaplaceProblem::solve () +{ + MGTransferPrebuilt > mg_transfer(hanging_node_constraints, mg_constrained_dofs); + mg_transfer.build_matrices(mg_dof_handler); + + FullMatrix coarse_matrix; + coarse_matrix.copy_from (mg_matrices[0]); + MGCoarseGridHouseholder<> coarse_grid_solver; + coarse_grid_solver.initialize (coarse_matrix); + + typedef PreconditionSOR > Smoother; + MGSmootherRelaxation, Smoother, Vector > + mg_smoother; + mg_smoother.initialize(mg_matrices); + mg_smoother.set_steps(2); + mg_smoother.set_symmetric(true); + + MGMatrix<> mg_matrix(&mg_matrices); + MGMatrix<> mg_interface_up(&mg_interface_in); + MGMatrix<> mg_interface_down(&mg_interface_in); + + Multigrid > mg(mg_dof_handler, + mg_matrix, + coarse_grid_solver, + 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); + + SolverControl solver_control (1000, 1e-12); + SolverCG<> cg (solver_control); + + solution = 0; + + cg.solve (system_matrix, solution, system_rhs, + preconditioner); + constraints.distribute (solution); + + deallog << " " << solver_control.last_step() + << " CG iterations needed to obtain convergence." + << std::endl; +} + + +template +void LaplaceProblem::refine_grid (const std::string& reftype) +{ + bool cell_refined = false; + if (reftype == "center" || !cell_refined) + { + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + 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/numbers::PI) + { + cell->set_refine_flag (); + cell_refined = true; + break; + } + } + } + } + if (reftype=="global" || !cell_refined) + triangulation.refine_global(1); + else + triangulation.execute_coarsening_and_refinement (); +} + + + +template +void LaplaceProblem::output_results (const unsigned int cycle) const +{ + DataOut data_out; + + data_out.attach_dof_handler (mg_dof_handler); + data_out.add_data_vector (solution, "solution"); + data_out.build_patches (); + + std::ostringstream filename; + filename << "solution-" + << cycle + << ".vtk"; + + std::ofstream output (filename.str().c_str()); + data_out.write_vtk (output); +} + + +template +void LaplaceProblem::run () +{ + for (unsigned int cycle=0; cycle<8; ++cycle) + { + deallog << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + { + GridGenerator::hyper_cube (triangulation, -1, 1, true); + +// static const HyperBallBoundary boundary; +// triangulation.set_boundary (0, boundary); + + triangulation.refine_global (1); + } + else + refine_grid ("center"); + + + deallog << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl; + + setup_system (); + + deallog << " Number of degrees of freedom: " + << mg_dof_handler.n_dofs() + << " (by level: "; + for (unsigned int level=0; level 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; +} diff --git a/tests/multigrid/step-16-bdry1.output b/tests/multigrid/step-16-bdry1.output new file mode 100644 index 0000000000..a0be674368 --- /dev/null +++ b/tests/multigrid/step-16-bdry1.output @@ -0,0 +1,57 @@ + +DEAL::Cycle 0: +DEAL:: Number of active cells: 4 +DEAL::Number of degrees of freedom: 9 L0: 4 L1: 9 +DEAL:: Number of degrees of freedom: 9 (by level: 4, 9) +DEAL:cg::Starting value 1.00000 +DEAL:cg::Convergence step 1 value 0.00000 +DEAL:: 1 CG iterations needed to obtain convergence. +DEAL::Cycle 1: +DEAL:: Number of active cells: 16 +DEAL::Number of degrees of freedom: 25 L0: 4 L1: 9 L2: 25 +DEAL:: Number of degrees of freedom: 25 (by level: 4, 9, 25) +DEAL:cg::Starting value 0.750000 +DEAL:cg::Convergence step 5 value 1.17215e-16 +DEAL:: 5 CG iterations needed to obtain convergence. +DEAL::Cycle 2: +DEAL:: Number of active cells: 28 +DEAL::Number of degrees of freedom: 41 L0: 4 L1: 9 L2: 25 L3: 25 +DEAL:: Number of degrees of freedom: 41 (by level: 4, 9, 25, 25) +DEAL:cg::Starting value 0.628894 +DEAL:cg::Convergence step 7 value 2.44070e-14 +DEAL:: 7 CG iterations needed to obtain convergence. +DEAL::Cycle 3: +DEAL:: Number of active cells: 40 +DEAL::Number of degrees of freedom: 57 L0: 4 L1: 9 L2: 25 L3: 25 L4: 25 +DEAL:: Number of degrees of freedom: 57 (by level: 4, 9, 25, 25, 25) +DEAL:cg::Starting value 0.620541 +DEAL:cg::Convergence step 7 value 5.62398e-13 +DEAL:: 7 CG iterations needed to obtain convergence. +DEAL::Cycle 4: +DEAL:: Number of active cells: 52 +DEAL::Number of degrees of freedom: 73 L0: 4 L1: 9 L2: 25 L3: 25 L4: 25 L5: 25 +DEAL:: Number of degrees of freedom: 73 (by level: 4, 9, 25, 25, 25, 25) +DEAL:cg::Starting value 0.620015 +DEAL:cg::Convergence step 8 value 4.43441e-14 +DEAL:: 8 CG iterations needed to obtain convergence. +DEAL::Cycle 5: +DEAL:: Number of active cells: 160 +DEAL::Number of degrees of freedom: 201 L0: 4 L1: 9 L2: 25 L3: 65 L4: 65 L5: 65 L6: 65 +DEAL:: Number of degrees of freedom: 201 (by level: 4, 9, 25, 65, 65, 65, 65) +DEAL:cg::Starting value 0.415999 +DEAL:cg::Convergence step 9 value 4.46942e-14 +DEAL:: 9 CG iterations needed to obtain convergence. +DEAL::Cycle 6: +DEAL:: Number of active cells: 304 +DEAL::Number of degrees of freedom: 357 L0: 4 L1: 9 L2: 25 L3: 81 L4: 81 L5: 81 L6: 81 L7: 153 +DEAL:: Number of degrees of freedom: 357 (by level: 4, 9, 25, 81, 81, 81, 81, 153) +DEAL:cg::Starting value 0.377840 +DEAL:cg::Convergence step 8 value 1.00219e-13 +DEAL:: 8 CG iterations needed to obtain convergence. +DEAL::Cycle 7: +DEAL:: Number of active cells: 676 +DEAL::Number of degrees of freedom: 761 L0: 4 L1: 9 L2: 25 L3: 81 L4: 81 L5: 81 L6: 121 L7: 209 L8: 465 +DEAL:: Number of degrees of freedom: 761 (by level: 4, 9, 25, 81, 81, 81, 121, 209, 465) +DEAL:cg::Starting value 0.377713 +DEAL:cg::Convergence step 8 value 5.29691e-13 +DEAL:: 8 CG iterations needed to obtain convergence. diff --git a/tests/multigrid/step-16.cc b/tests/multigrid/step-16.cc index ab5c4b41db..b975cb1807 100644 --- a/tests/multigrid/step-16.cc +++ b/tests/multigrid/step-16.cc @@ -14,11 +14,7 @@ // // --------------------------------------------------------------------- -/* - * Authors: Guido Kanschat, University of Heidelberg, 2003 - * Baerbel Janssen, University of Heidelberg, 2010 - * Wolfgang Bangerth, Texas A&M University, 2010 - */ +// Multigrid for continuous finite elements without MeshWorker #include "../tests.h" #include -- 2.39.5