From: wolf Date: Wed, 9 Feb 2000 20:01:36 +0000 (+0000) Subject: Collect all multigrid files in one directory, rename MG{Matrix,Vector} to MGLevelObje... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e88634ff8e36f247a05f334fc1028cff9bfae4d0;p=dealii-svn.git Collect all multigrid files in one directory, rename MG{Matrix,Vector} to MGLevelObject, and small other changes. git-svn-id: https://svn.dealii.org/trunk@2356 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/Attic/examples/multigrid/multigrid.cc b/deal.II/deal.II/Attic/examples/multigrid/multigrid.cc index e48fd7e987..7769168471 100644 --- a/deal.II/deal.II/Attic/examples/multigrid/multigrid.cc +++ b/deal.II/deal.II/Attic/examples/multigrid/multigrid.cc @@ -2,712 +2,514 @@ /* Copyright W. Bangerth, University of Heidelberg, 1998 */ -#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 +#include +#include +#include +#include +#include -#include -#include -#include -#include -#include - - - - - -template -class PoissonEquation : public Equation { - public: - PoissonEquation (const Function &rhs) : - Equation(1), - right_hand_side (rhs) {}; - - virtual void assemble (FullMatrix &cell_matrix, - Vector &rhs, - const FEValues &fe_values, - const DoFHandler::cell_iterator &cell) const; - virtual void assemble (FullMatrix &cell_matrix, - const FEValues &fe_values, - const DoFHandler::cell_iterator &cell) const; - virtual void assemble (Vector &rhs, - const FEValues &fe_values, - const DoFHandler::cell_iterator &cell) const; - protected: - const Function &right_hand_side; -}; - +#include +#include +#include +#include +#include +#include +#include +#include template -class PoissonProblem { +class LaplaceProblem +{ public: - /** - * Declare a data type which denotes a - * mapping between a boundary indicator - * and the function denoting the boundary - * values on this part of the boundary. - * Only one boundary function may be given - * for each boundary indicator, which is - * guaranteed by the #map# data type. - * - * See the general documentation of this - * class for more detail. - */ - typedef map*> FunctionMap; - /** - * Typdedef an iterator which assembles - * matrices and vectors. - */ - typedef TriaActiveIterator > active_assemble_iterator; - - PoissonProblem (unsigned int order); - virtual ~PoissonProblem (); + LaplaceProblem (); + ~LaplaceProblem (); + void run (); - void clear (); - void create_new (); + private: + void setup_system (); + void assemble_system (); void solve (); - - /** - * Initiate the process of assemblage of - * vectors and system matrix. Use the - * given equation object and the given - * quadrature formula. Also use the list - * of dirichlet boundary value functions - * (by default, no dirichlet bc are assumed - * which means that all bc are included - * into the weak formulation). - * - * For what exactly happens here, refer to - * the general doc of this class. - */ - virtual void assemble (const Equation &equation, - const Quadrature &q, - const UpdateFlags update_flags, - const FunctionMap &dirichlet_bc = FunctionMap()); - - int run (unsigned int level); - void print_history (string filename) const; - - /** - * Exception - */ - DeclException0 (ExcNoTriaSelected); - - protected: - Triangulation *tria; - MGDoFHandler *dof; - - - /** - * Sparsity pattern of the system matrix. - */ - SparsityPattern system_sparsity; - - /** - * System matrix. - */ - SparseMatrix system_matrix; + void refine_grid (); + void output_results (const unsigned int cycle) const; - /** - * Vector storing the right hand side. - */ - Vector right_hand_side; + Triangulation triangulation; + MGDoFHandler mg_dof_handler; - /** - * Solution vector. - */ - Vector solution; + FEQ1 fe; - /** - * List of constraints introduced by - * hanging nodes. - */ - ConstraintMatrix constraints; + ConstraintMatrix hanging_node_constraints; - Function *rhs; - Function *boundary_values; + SparsityPattern global_sparsity_pattern; + SparseMatrix global_system_matrix; - vector l1_error, l2_error, linfty_error, h1_seminorm_error, h1_error; - vector n_dofs; - - unsigned int order; + MGLevelObject level_sparsity_patterns; + MGLevelObject > level_system_matrices; + + Vector solution; + Vector system_rhs; }; - - -/** - Right hand side constructed such that the exact solution is - $sin(2 pi x) + sin(2 pi y)$ - */ template -class RHSPoly : public Function { +class Coefficient : public Function +{ public: - /** - * Return the value of the function - * at the given point. - */ - virtual double value (const Point &p, - const unsigned int component) const; + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual void value_list (const vector > &points, + vector &values, + const unsigned int component = 0) const; }; template -class Solution : public Function { - public: - /** - * Return the value of the function - * at the given point. - */ - virtual double value (const Point &p, - const unsigned int component) const; - /** - * Return the gradient of the function - * at the given point. - */ - virtual Tensor<1,dim> gradient (const Point &p, - const unsigned int component) const; +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 vector > &points, + vector &values, + const unsigned int component) const +{ + const unsigned int n_points = points.size(); -template <> -double RHSPoly<2>::value (const Point<2> &p, - const unsigned int component) const { - Assert (component==0, ExcIndexRange (component, 0, 1)); + Assert (values.size() == n_points, + ExcVectorHasWrongSize (values.size(), n_points)); + + Assert (component == 0, + ExcWrongComponent (component, 1)); - const double x = p(0), - y = p(1); - const double pi= 3.1415926536; - return 4*pi*pi*(sin(2*pi*x)+sin(2*pi*y)); + for (unsigned int i=0; i -double Solution<2>::value (const Point<2> &p, - const unsigned int component) const { - Assert (component==0, ExcIndexRange (component, 0, 1)); - const double x = p(0), - y = p(1); - const double pi= 3.1415926536; - return sin(2*pi*x)+sin(2*pi*y); +class MGSmootherLAC : public MGSmootherBase +{ + private: + SmartPointer > >matrices; + public: + MGSmootherLAC(MGLevelObject >&); + + virtual void smooth (const unsigned int level, + Vector &u, + const Vector &rhs) const; }; -template <> -Tensor<1,2> Solution<2>::gradient (const Point<2> &p, - const unsigned int component) const { - Assert (component==0, ExcIndexRange (component, 0, 1)); +MGSmootherLAC::MGSmootherLAC(MGLevelObject >& matrix) + : + matrices(&matrix) +{} - const double x = p(0), - y = p(1); - const double pi= 3.1415926536; - return Point<2> (2*pi*cos(2*pi*x), - 2*pi*cos(2*pi*y)); -}; - +void +MGSmootherLAC::smooth (const unsigned int level, + Vector &u, + const Vector &rhs) const +{ + SolverControl control(2,1.e-300,false,false); + PrimitiveVectorMemory<> mem; + SolverRichardson<> rich(control, mem); + PreconditionRelaxation<> + prec((*matrices)[level], &SparseMatrix ::template precondition_SSOR, 1.); + rich.solve((*matrices)[level], u, rhs, prec); +} -template <> -void PoissonEquation<2>::assemble (FullMatrix &cell_matrix, - Vector &rhs, - const FEValues<2> &fe_values, - const DoFHandler<2>::cell_iterator &) const { - for (unsigned int point=0; point +LaplaceProblem::LaplaceProblem () : + mg_dof_handler (triangulation) +{}; template -void PoissonEquation::assemble (FullMatrix &, - const FEValues &, - const DoFHandler::cell_iterator &) const { - Assert (false, ExcPureVirtualFunctionCalled()); +LaplaceProblem::~LaplaceProblem () +{ + mg_dof_handler.clear (); }; template -void PoissonEquation::assemble (Vector &, - const FEValues &, - const DoFHandler::cell_iterator &) const { - Assert (false, ExcPureVirtualFunctionCalled()); -}; - +void LaplaceProblem::setup_system () +{ + mg_dof_handler.distribute_dofs (fe); + hanging_node_constraints.clear (); + DoFTools::make_hanging_node_constraints (mg_dof_handler, + hanging_node_constraints); + hanging_node_constraints.close (); + global_sparsity_pattern.reinit (mg_dof_handler.DoFHandler::n_dofs(), + mg_dof_handler.DoFHandler::n_dofs(), + mg_dof_handler.max_couplings_between_dofs()); + DoFTools::make_sparsity_pattern (mg_dof_handler, global_sparsity_pattern); + hanging_node_constraints.condense (global_sparsity_pattern); + global_sparsity_pattern.compress(); + global_system_matrix.reinit (global_sparsity_pattern); + solution.reinit (mg_dof_handler.DoFHandler::n_dofs()); + system_rhs.reinit (mg_dof_handler.DoFHandler::n_dofs()); + const unsigned int n_levels = triangulation.n_levels(); + level_system_matrices.resize (0, n_levels); + level_sparsity_patterns.resize (0, n_levels); + + for (unsigned int level=0; level -PoissonProblem::PoissonProblem (unsigned int order) : - tria(0), dof(0), rhs(0), - boundary_values(0), order(order) {}; +void LaplaceProblem::assemble_system () +{ + const Coefficient coefficient; + QGauss2 quadrature_formula; + FEValues fe_values (fe, quadrature_formula, + UpdateFlags(update_values | + update_gradients | + update_q_points | + update_JxW_values)); -template -PoissonProblem::~PoissonProblem () -{ - clear (); -}; + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.n_quadrature_points; + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + vector local_dof_indices (dofs_per_cell); -template -void PoissonProblem::clear () { - if (dof != 0) { - delete dof; - dof = 0; - }; - - if (tria != 0) { - delete tria; - tria = 0; - }; - - if (rhs != 0) - { - delete rhs; - rhs = 0; - }; + // FIX + vector coefficient_values (n_q_points, 1.0); - if (boundary_values != 0) + // not only active cells + MGDoFHandler::cell_iterator cell = mg_dof_handler.begin(), + endc = mg_dof_handler.end(); + for (; cell!=endc; ++cell) { - delete boundary_values; - boundary_values = 0; + cell_matrix.clear (); + cell_rhs.clear (); + + fe_values.reinit (cell); + const FullMatrix + & shape_values = fe_values.get_shape_values(); + const vector > > + & shape_grads = fe_values.get_shape_grads(); + const vector + & JxW_values = fe_values.get_JxW_values(); + const vector > + & q_points = fe_values.get_quadrature_points(); + + // FIX +// coefficient.value_list (q_points, coefficient_values); + + for (unsigned int q_point=0; q_pointget_mg_dof_indices (local_dof_indices); + const unsigned int level = cell->level(); + for (unsigned int i=0; iactive()) + { + 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, +// global_system_matrix, +// solution, +// system_rhs); }; - template -void PoissonProblem::create_new () { - clear (); - - tria = new Triangulation(); - dof = new MGDoFHandler (*tria); -}; +void LaplaceProblem::solve () +{ + { + SolverControl solver_control (1000, 1e-12); + PrimitiveVectorMemory<> vector_memory; + SolverCG<> cg (solver_control, vector_memory); + SolverControl coarse_grid_solver_control (1000, 1e-12); + PrimitiveVectorMemory<> coarse_grid_vector_memory; + + SolverCG<> coarse_grid_cg (coarse_grid_solver_control, + coarse_grid_vector_memory); + +// PreconditionRelaxation<> +// coarse_grid_solver_preconditioner(level_system_matrices[level_system_matrices.get_minlevel()], +// &SparseMatrix::template precondition_SSOR, +// 1.2); + PreconditionIdentity coarse_grid_solver_preconditioner; + + MGCoarseGridLACIteration, SparseMatrix, PreconditionIdentity> + coarse_grid_solver (coarse_grid_cg, + level_system_matrices[level_system_matrices.get_minlevel()], + coarse_grid_solver_preconditioner); + + MGSmootherLAC smoother (level_system_matrices); + MGTransferPrebuilt grid_transfer; + grid_transfer.build_matrices (mg_dof_handler); + + Multigrid<2> multigrid (mg_dof_handler, + hanging_node_constraints, + level_sparsity_patterns, + level_system_matrices, + grid_transfer); + + PreconditionMG > + mg_precondition (multigrid, smoother, smoother, coarse_grid_solver); -template -void PoissonProblem::assemble (const Equation &equation, - const Quadrature &quadrature, - const UpdateFlags update_flags, - const FunctionMap &dirichlet_bc) { - Assert ((tria!=0) && (dof!=0), ExcNoTriaSelected()); - - system_sparsity.reinit (dof->DoFHandler::n_dofs(), - dof->DoFHandler::n_dofs(), - dof->max_couplings_between_dofs()); - right_hand_side.reinit (dof->DoFHandler::n_dofs()); - - // make up sparsity pattern and - // compress with constraints - constraints.clear (); - DoFTools::make_hanging_node_constraints (*dof, constraints); - constraints.close (); - DoFTools::make_sparsity_pattern (*dof, system_sparsity); - constraints.condense (system_sparsity); - - MGTransferPrebuilt p; - p.build_matrices (*dof); -// MGSmoother smoother(*dof); - + solution.clear (); + cg.solve (global_system_matrix, solution, system_rhs, + mg_precondition); - // reinite system matrix - system_matrix.reinit (system_sparsity); - // reinit solution vector, preset - // with zeroes. - solution.reinit (dof->DoFHandler::n_dofs()); + cout << " MG Outer iterations: " << solver_control.last_step() + << endl; + + cout << " MG Total inner iterations: " << coarse_grid_solver_control.last_step() + << endl; + }; - // create assembler - Assembler::AssemblerData data (*dof, - true, true, //assemble matrix and rhs - system_matrix, - right_hand_side, - quadrature, - update_flags); - active_assemble_iterator assembler (tria, - tria->begin_active()->level(), - tria->begin_active()->index(), - &data); - // loop over all cells, fill matrix and rhs - do { - assembler->assemble (equation); - } - while ((++assembler).state() == valid); - - // condense system matrix in-place - constraints.condense (system_matrix); - - // condense right hand side in-place - constraints.condense (right_hand_side); - - // apply Dirichlet bc as described - // in the docs - map boundary_value_list; - - for (FunctionMap::const_iterator bc=dirichlet_bc.begin(); - bc != dirichlet_bc.end(); ++bc) - VectorTools::interpolate_boundary_values (*dof, - bc->first, - *bc->second, - boundary_value_list); - MatrixTools::apply_boundary_values (boundary_value_list, - system_matrix, solution, - right_hand_side); + SolverControl solver_control (1000, 1e-12); + PrimitiveVectorMemory<> vector_memory; + SolverCG<> cg (solver_control, vector_memory); + + PreconditionRelaxation<> + preconditioner(global_system_matrix, + &SparseMatrix::template precondition_SSOR, + 1.2); + + solution.clear (); + cg.solve (global_system_matrix, solution, system_rhs, + preconditioner); + + cout << " CG Outer iterations: " << solver_control.last_step() + << endl; + }; + + hanging_node_constraints.distribute (solution); +}; + + +template +void LaplaceProblem::refine_grid () +{ + Vector estimated_error_per_cell (triangulation.n_active_cells()); + + KellyErrorEstimator::FunctionMap neumann_boundary; + KellyErrorEstimator::estimate (mg_dof_handler, + QGauss3(), + neumann_boundary, + solution, + estimated_error_per_cell); + + triangulation.refine_and_coarsen_fixed_number (estimated_error_per_cell, + 0.3, 0.03); + triangulation.execute_coarsening_and_refinement (); }; template -void PoissonProblem::solve () { - Assert ((tria!=0) && (dof!=0), ExcNoTriaSelected()); +void LaplaceProblem::output_results (const unsigned int cycle) const +{ + string filename = "grid-"; + filename += ('0' + cycle); + Assert (cycle < 10, ExcInternalError()); - SolverControl control(4000, 1e-16); - PrimitiveVectorMemory > memory; - SolverCG,Vector > cg(control,memory); - - // solve - cg.solve (system_matrix, solution, right_hand_side, - PreconditionIdentity()); - // distribute solution - constraints.distribute (solution); + filename += ".eps"; + ofstream output (filename.c_str()); + + GridOut grid_out; + grid_out.write_eps (triangulation, output); }; template -int PoissonProblem::run (const unsigned int level) { - create_new (); - - cout << "Refinement level = " << level - << ", using elements of type <"; - switch (order) +void LaplaceProblem::run () +{ + for (unsigned int cycle=0; cycle<8; ++cycle) { - case 0: - cout << "criss-cross"; - break; - default: - cout << "Lagrange-" << order; - break; - }; - cout << ">" << endl; - - cout << " Making grid... "; - GridGenerator::hyper_cube (*tria); - tria->refine_global (level+1); - tria->begin_active()->set_refine_flag(); - (++(++(tria->begin_active())))->set_refine_flag(); - tria->execute_coarsening_and_refinement (); - cout << tria->n_active_cells() << " active cells." << endl; - - rhs = new RHSPoly(); - boundary_values = new Solution (); - + cout << "Cycle " << cycle << ':' << endl; - FiniteElement *fe; - PoissonEquation equation (*rhs); - Quadrature *quadrature; - Quadrature *boundary_quadrature; - switch (order) { - case 0: - fe = new FECrissCross(); - quadrature = new QCrissCross1(); - boundary_quadrature = new QGauss2(); - break; - case 1: - fe = new FEQ1(); - quadrature = new QGauss3(); - boundary_quadrature = new QGauss2(); - break; - case 2: - fe = new FEQ2(); - quadrature = new QGauss4(); - boundary_quadrature = new QGauss3(); - break; - case 3: - fe = new FEQ3(); - quadrature = new QGauss5(); - boundary_quadrature = new QGauss4(); - break; - case 4: - fe = new FEQ4(); - quadrature = new QGauss6(); - boundary_quadrature = new QGauss5(); - break; - default: - return 100000; - }; - - cout << " Distributing dofs... "; - dof->distribute_dofs (*fe); - cout << dof->DoFHandler::n_dofs() << " degrees of freedom." << endl; - n_dofs.push_back (dof->DoFHandler::n_dofs()); - - cout << " Assembling matrices..." << endl; - UpdateFlags update_flags = UpdateFlags(update_values | update_q_points | - update_gradients | update_JxW_values); - - map*> dirichlet_bc; - dirichlet_bc[0] = boundary_values; - assemble (equation, *quadrature, update_flags, dirichlet_bc); + if (cycle == 0) + { + GridGenerator::hyper_cube (triangulation); + triangulation.refine_global (1); + } + else + { + refine_grid (); + }; + - cout << " Solving..." << endl; - solve (); + cout << " Number of active cells: " + << triangulation.n_active_cells() + << endl; - Solution sol; - Vector l1_error_per_cell, l2_error_per_cell, linfty_error_per_cell; - Vector h1_seminorm_error_per_cell, h1_error_per_cell; - - cout << " Calculating L1 error... "; - VectorTools::integrate_difference (*dof, - solution, sol, - l1_error_per_cell, - *quadrature, L1_norm); - cout << l1_error_per_cell.l1_norm() << endl; - l1_error.push_back (l1_error_per_cell.l1_norm()); - - cout << " Calculating L2 error... "; - VectorTools::integrate_difference (*dof, - solution, sol, - l2_error_per_cell, - *quadrature, L2_norm); - cout << l2_error_per_cell.l2_norm() << endl; - l2_error.push_back (l2_error_per_cell.l2_norm()); - - cout << " Calculating L-infinity error... "; - VectorTools::integrate_difference (*dof, - solution, sol, - linfty_error_per_cell, - *quadrature, Linfty_norm); - cout << linfty_error_per_cell.linfty_norm() << endl; - linfty_error.push_back (linfty_error_per_cell.linfty_norm()); - - cout << " Calculating H1-seminorm error... "; - VectorTools::integrate_difference (*dof, - solution, sol, - h1_seminorm_error_per_cell, - *quadrature, H1_seminorm); - cout << h1_seminorm_error_per_cell.l2_norm() << endl; - h1_seminorm_error.push_back (h1_seminorm_error_per_cell.l2_norm()); - - cout << " Calculating H1 error... "; - VectorTools::integrate_difference (*dof, - solution, sol, - h1_error_per_cell, - *quadrature, H1_norm); - cout << h1_error_per_cell.l2_norm() << endl; - h1_error.push_back (h1_error_per_cell.l2_norm()); - - if (dof->DoFHandler::n_dofs()<=5000) - { - Vector l1_error_per_dof (dof->DoFHandler::n_dofs()); - Vector l2_error_per_dof (dof->DoFHandler::n_dofs()); - Vector linfty_error_per_dof (dof->DoFHandler::n_dofs()); - Vector h1_seminorm_error_per_dof (dof->DoFHandler::n_dofs()); - Vector h1_error_per_dof (dof->DoFHandler::n_dofs()); - DoFTools::distribute_cell_to_dof_vector (*dof, l1_error_per_cell, l1_error_per_dof); - DoFTools::distribute_cell_to_dof_vector (*dof, l2_error_per_cell, l2_error_per_dof); - DoFTools::distribute_cell_to_dof_vector (*dof, linfty_error_per_cell, - linfty_error_per_dof); - DoFTools::distribute_cell_to_dof_vector (*dof, h1_seminorm_error_per_cell, - h1_seminorm_error_per_dof); - DoFTools::distribute_cell_to_dof_vector (*dof, h1_error_per_cell, h1_error_per_dof); - -// Vector projected_solution; -// ConstraintMatrix constraints; -// constraints.close (); -// VectorTools::project (*dof, constraints, *fe, -// StraightBoundary(), *quadrature, -// sol, projected_solution, false, -// *boundary_quadrature); -// cout << " Calculating L2 error of projected solution... "; -// VectorTools::integrate_difference (*dof, -// projected_solution, sol, -// l2_error_per_cell, -// *quadrature, *fe, L2_norm); -// cout << l2_error_per_cell.l2_norm() << endl; - - - string filename; - filename = ('0'+order); - filename += "."; - filename += ('0'+level); - filename += ".ucd"; - cout << " Writing error plots to <" << filename << ">..." << endl; + setup_system (); + + cout << " Number of degrees of freedom: " + << mg_dof_handler.DoFHandler::n_dofs() + << endl; - DataOut out; - ofstream o(filename.c_str()); - out.attach_dof_handler (*dof); - - out.add_data_vector (solution, "u"); - out.add_data_vector (l1_error_per_dof, "L1_Error"); - out.add_data_vector (l2_error_per_dof, "L2_Error"); - out.add_data_vector (linfty_error_per_dof, "Linfty_Error"); - out.add_data_vector (h1_seminorm_error_per_dof, "H1_seminorm_Error"); - out.add_data_vector (h1_error_per_dof, "H1_Error"); - out.write_ucd (o); - o.close (); - } - else - cout << " Not writing error as grid." << endl; - - cout << endl; - - const unsigned int n_dofs = dof->DoFHandler::n_dofs(); - dof->clear (); - tria->set_boundary (0); - delete fe; - delete quadrature; - delete boundary_quadrature; - - return n_dofs; -}; + assemble_system (); + solve (); + output_results (cycle); + DataOut::EpsFlags eps_flags; + eps_flags.z_scaling = 4; + + DataOut data_out; + data_out.set_flags (eps_flags); -template -void PoissonProblem::print_history (string filename) const { - ofstream out(filename.c_str()); - out << "# n_dofs l1_error l2_error linfty_error h1_seminorm_error h1_error" - << endl; - for (unsigned int i=0; ih/2:" << endl; - cout << " L1 error : " << 1./average_l1 << endl - << " L2 error : " << 1./average_l2 << endl - << " Linfty error : " << 1./average_linfty << endl - << " H1 seminorm error: " << 1./average_h1_semi << endl - << " H1 error : " << 1./average_h1 << endl; - cout << "==========================================================\n"; - cout << "==========================================================\n"; }; + +int main () +{ + try + { + deallog.depth_console (0); - -int main () { - for (unsigned int order=0; order<5; ++order) + LaplaceProblem<2> laplace_problem_2d; + laplace_problem_2d.run (); + } + catch (exception &exc) { - PoissonProblem<2> problem (order); - - unsigned int level=0; - unsigned int n_dofs; - do - n_dofs = problem.run (level++); - while (n_dofs<25000); - - string filename; - switch (order) - { - case 0: - filename = "criss_cross"; - break; - case 1: - filename = "linear"; - break; - case 2: - filename = "quadratic"; - break; - case 3: - filename = "cubic"; - break; - case 4: - filename = "quartic"; - break; - }; - filename += ".history"; - - cout << endl << "Printing convergence history to <" - << filename << ">..." << endl; - problem.print_history (filename); - cout << endl << endl << endl; + cerr << endl << endl + << "----------------------------------------------------" + << endl; + cerr << "Exception on processing: " << endl + << exc.what() << endl + << "Aborting!" << endl + << "----------------------------------------------------" + << endl; + return 1; + } + catch (...) + { + cerr << endl << endl + << "----------------------------------------------------" + << endl; + cerr << "Unknown exception!" << endl + << "Aborting!" << endl + << "----------------------------------------------------" + << endl; + return 1; }; - + return 0; }; diff --git a/deal.II/deal.II/include/dofs/mg_dof_accessor.h b/deal.II/deal.II/include/dofs/mg_dof_accessor.h deleted file mode 100644 index 4b30be0c8a..0000000000 --- a/deal.II/deal.II/include/dofs/mg_dof_accessor.h +++ /dev/null @@ -1,670 +0,0 @@ -/*---------------------------- mg_dof_accessor.h ---------------------------*/ -/* $Id$ */ -/* Copyright W. Bangerth, University of Heidelberg, 1998 */ -#ifndef __mg_dof_accessor_H -#define __mg_dof_accessor_H -/*---------------------------- mg_dof_accessor.h ---------------------------*/ - - -#include - - - - -/** - * Define the basis for accessors to the degrees of freedom for - * a multigrid DoF handler object. - * - * Note that it is allowed to construct an object of which the - * #mg_dof_handler# pointer is a Null pointer. Such an object would - * result in a strange kind of behaviour, though every reasonable - * operating system should disallow access through that pointer. - * The reason we do not check for the null pointer in the - * constructor which gets passed the #DoFHandler# pointer is that - * if we did we could not make dof iterators member of other classes - * (like in the #FEValues# class) if we did not know about the - * #DoFHandler# object to be used upon construction of that object. - * Through the way this class is implemented here, we allow the - * creation of a kind of virgin object which only gets useful if - * assigned to from another object before first usage. - * - * Opposite to construction, it is not possible to copy an object - * which has an invalid dof handler pointer. This is to guarantee - * that every iterator which is once assigned to is a valid - * object. However, this assertion only holds in debug mode, when - * the #Assert# macro is switched on. - * - * @author Wolfgang Bangerth, 1998 - */ -template -class MGDoFAccessor { - public: - /** - * Constructor - */ - MGDoFAccessor () : mg_dof_handler(0) { - Assert (false, ExcInvalidObject()); - }; - - /** - * This should be the default constructor. - * We cast away the #const#ness of the - * pointer which clearly is EVIL but - * we can't help without making all - * functions which could somehow use - * iterators (directly or indirectly) make - * non-const, even if they preserve - * constness. - */ - MGDoFAccessor (const MGDoFHandler *mg_dof_handler) : - mg_dof_handler(const_cast*>(mg_dof_handler)) {}; - - /** - * Reset the DoF handler pointer. - */ - void set_mg_dof_handler (MGDoFHandler *dh) { - Assert (dh != 0, DoFAccessor::ExcInvalidObject()); - mg_dof_handler = dh; - }; - - /** - * Exception for child classes - */ - DeclException0 (ExcInvalidObject); - - /** - * Copy operator. - */ - MGDoFAccessor & operator = (const MGDoFAccessor &da) { - set_dof_handler (da.mg_dof_handler); - return *this; - }; - - protected: - /** - * Store the address of the #MGDoFHandler# object - * to be accessed. - */ - MGDoFHandler *mg_dof_handler; -}; - - -/* -------------------------------------------------------------------------- */ - -/** - * This is a switch class which only declares a #typdef#. It is meant to - * determine which class a #MGDoFObjectAccessor# class is to be derived - * from. By default, #MGDoFObjectAccessor# derives from - * the #typedef# in the general #MGDoFObjectAccessor_Inheritance# - * class, which is #DoFObjectAccessor#, - * but if #celldim==dim#, then the specialization #MGDoFObjectAccessor_Inheritance# - * is used which declares its local type to be #DoFCellAccessor#. Therefore, - * the inheritance is automatically chosen to be from #DoFCellAccessor# if the - * object under consideration has full dimension, i.e. constitutes a cell. - * - * @author Wolfgang Bangerth, 1999 - */ -template -class MGDoFObjectAccessor_Inheritance -{ - /** - * Declaration of the #typedef#. - * See the full documentation for - * more information. - */ - typedef DoFObjectAccessor BaseClass; -}; - - - -/** - * This is a switch class which only declares a #typdef#. It is meant to - * determine which class a #DoFObjectAccessor# class is to be derived - * from. By default, #DoFObjectAccessor# derives from - * the #typedef# in the general #DoFObjectAccessor_Inheritance# - * class, which is #TriaObjectAccessor#, - * but if #celldim==dim#, then the specialization #TriaObjectAccessor# - * is used which declares its local type to be #CellAccessor#. Therefore, - * the inheritance is automatically chosen to be from #CellAccessor# if the - * object under consideration has full dimension, i.e. constitutes a cell. - * - * @author Wolfgang Bangerth, 1999 - */ -template -class MGDoFObjectAccessor_Inheritance -{ - /** - * Declaration of the #typedef#. - * See the full documentation for - * more information. - */ - typedef DoFCellAccessor BaseClass; -}; - - - - -/* -------------------------------------------------------------------------- */ - - - -/** - * Common template for line, quad, hex. - * - * Internal: inheritance is necessary for the general template due to - * a compiler error. - * @author Guido Kanschat, 1999 - */ -template -class MGDoFObjectAccessor : public MGDoFAccessor, - public MGDoFObjectAccessor_Inheritance::BaseClass -{}; - - - -/** - * Closure class. - */ -template -class MGDoFObjectAccessor<0, dim> -{ - public: - typedef void AccessorData; - - /** - * Constructor. Should never be called - * and thus throws an error. - */ - MGDoFObjectAccessor (Triangulation *, - const int, - const int, - const AccessorData *) - { - Assert (false, ExcInternalError()); - } -}; - - - -/** - * Grant access to the degrees of freedom located on lines. - * This class follows mainly the route laid out by the accessor library - * declared in the triangulation library (\Ref{TriaAccessor}). It enables - * the user to access the degrees of freedom on the lines (there are similar - * versions for the DoFs on quads, etc), where the dimension of the underlying - * triangulation does not really matter (i.e. this accessor works with the - * lines in 1D-, 2D-, etc dimensions). - * - * - * \subsection{Usage} - * - * The \Ref{DoFDimensionInfo} classes inherited by the \Ref{DoFHandler} classes - * declare typedefs to iterators using the accessors declared in this class - * hierarchy tree. Usage is best to happens through these typedefs, since they - * are more secure to changes in the class naming and template interface as well - * as they provide easier typing (much less complicated names!). - * - * - * \subsection{Notes about the class hierarchy structure} - * - * Inheritance from #MGDoFObjectAccessor_Inheritance<1,dim>::BaseClass# yields - * inheritance from #DoFCellAccessor<1># if #dim==1# and from - * #TriaObjectAccessor<1,dim># for all other #dim# values. Thus, an object - * of this class shares all features of cells in one dimension, but behaves - * like an ordinary line in all other cases. - * - * @author Wolfgang Bangerth, 1998 - */ -template -class MGDoFObjectAccessor<1, dim> : public MGDoFAccessor, - public MGDoFObjectAccessor_Inheritance<1, dim>::BaseClass -{ - public: - /** - * Declare the data type that this accessor - * class expects to get passed from the - * iterator classes. - */ - typedef MGDoFHandler AccessorData; - - /** - * Default constructor, unused thus - * not implemented. - */ - MGDoFObjectAccessor (); - - /** - * Constructor. The #local_data# - * argument is assumed to be a pointer - * to a #MGDoFHandler# object. - */ - MGDoFObjectAccessor (Triangulation *tria, - const int level, - const int index, - const AccessorData *local_data); - - /** - * Return the index of the #i#th degree - * of freedom of this line on the level - * this line lives on. - */ - int mg_dof_index (const unsigned int i) const; - - /** - * Set the index of the #i#th degree - * of freedom of this line on the - * level this line lives on to #index#. - */ - void set_mg_dof_index (const unsigned int i, const int index) const; - - /** - * Return the index of the #i#th degree - * on the #vertex#th vertex for the - * level this line lives on. - */ - int mg_vertex_dof_index (const unsigned int vertex, - const unsigned int i) const; - - /** - * Set the index of the #i#th degree - * on the #vertex#th vertex to #index# - * for the level this line lives on. - */ - void set_mg_vertex_dof_index (const unsigned int vertex, - const unsigned int i, - const int index) const; - - /** - * Return the indices of the dofs of this - * line in the standard ordering: dofs - * on vertex 0, dofs on vertex 1, - * dofs on line. - * - * It is assumed that the vector already - * has the right size beforehand. The - * indices refer to the local numbering - * for the level this line lives on. - */ - void get_mg_dof_indices (vector &dof_indices) const; - - /** - * Return the value of the given vector - * restricted to the dofs of this - * cell in the standard ordering: dofs - * on vertex 0, dofs on vertex 1, etc, - * dofs on line 0, dofs on line 1, etc, - * dofs on quad 0, etc. - * - * It is assumed that the vector already - * has the right size beforehand. The - * indices refer to the multilevel - * numbering local to the present - * level of this cell. The vector shall - * therefore have the same number of - * entries as there are degrees of - * freedom on this level. - */ - template - void get_mg_dof_values (const Vector &values, - Vector &dof_values) const; - - /** - * Return the #i#th child as a MGDoF line - * iterator. This function is needed since - * the child function of the base - * class returns a line accessor without - * access to the MG data. - */ - TriaIterator > child (const unsigned int) const; - - /** - * Implement the copy operator needed - * for the iterator classes. - */ - void copy_from (const MGDoFObjectAccessor<1, dim> &a); -}; - - - - - - -/** - * Grant access to the multilevel degrees of freedom located on quads. - * - * @see DoFLineAccessor - */ -template -class MGDoFObjectAccessor<2, dim> : public MGDoFAccessor, - public MGDoFObjectAccessor_Inheritance<2, dim>::BaseClass -{ - public: - /** - * Declare the data type that this accessor - * class expects to get passed from the - * iterator classes. - */ - typedef MGDoFHandler AccessorData; - - /** - * Default constructor, unused thus - * not implemented. - */ - MGDoFObjectAccessor (); - - /** - * Constructor. The #local_data# - * argument is assumed to be a pointer - * to a #DoFHandler# object. - */ - MGDoFObjectAccessor (Triangulation *tria, - const int level, - const int index, - const AccessorData *local_data); - - /** - * Return the index of the #i#th degree - * of freedom of this quad on the level - * this quad lives on. - */ - int mg_dof_index (const unsigned int i) const; - - /** - * Set the index of the #i#th degree - * of freedom of this quad on the - * level this quad lives on to #index#. - */ - void set_mg_dof_index (const unsigned int i, const int index) const; - - /** - * Return the index of the #i#th degree - * on the #vertex#th vertex for the level - * this quad lives on. - */ - int mg_vertex_dof_index (const unsigned int vertex, - const unsigned int i) const; - - /** - * Set the index of the #i#th degree - * on the #vertex#th vertex for the - * level this quad lives on to #index#. - */ - void set_mg_vertex_dof_index (const unsigned int vertex, - const unsigned int i, - const int index) const; - - /** - * Return the indices of the dofs of this - * quad in the standard ordering: dofs - * on vertex 0, dofs on vertex 1, etc, - * dofs on line 0, dofs on line 1, etc, - * dofs on quad 0, etc. - * - * It is assumed that the vector already - * has the right size beforehand. The - * indices refer to the local numbering - * for the level this quad lives on. - */ - void get_mg_dof_indices (vector &dof_indices) const; - - /** - * Return the value of the given vector - * restricted to the dofs of this - * cell in the standard ordering: dofs - * on vertex 0, dofs on vertex 1, etc, - * dofs on line 0, dofs on line 1, etc, - * dofs on quad 0, etc. - * - * It is assumed that the vector already - * has the right size beforehand. The - * indices refer to the multilevel - * numbering local to the present - * level of this cell. The vector shall - * therefore have the same number of - * entries as there are degrees of - * freedom on this level. - */ - template - void get_mg_dof_values (const Vector &values, - Vector &dof_values) const; - - /** - * Return a pointer to the #i#th line - * bounding this #Quad#. - */ - TriaIterator > - line (const unsigned int i) const; - - /** - * Return the #i#th child as a DoF quad - * iterator. This function is needed since - * the child function of the base - * class returns a quad accessor without - * access to the DoF data. - */ - TriaIterator > - child (const unsigned int) const; - - /** - * Implement the copy operator needed - * for the iterator classes. - */ - void copy_from (const MGDoFObjectAccessor<2, dim> &a); -}; - - - - -/** - * Grant access to the multilevel degrees of freedom located on hexhedra. - * - * @see DoFLineAccessor - */ -template -class MGDoFObjectAccessor<3, dim> : public MGDoFAccessor, - public MGDoFObjectAccessor_Inheritance<3, dim>::BaseClass -{ - public: - /** - * Declare the data type that this accessor - * class expects to get passed from the - * iterator classes. - */ - typedef MGDoFHandler AccessorData; - - /** - * Default constructor, unused thus - * not implemented. - */ - MGDoFObjectAccessor (); - - /** - * Constructor. The #local_data# - * argument is assumed to be a pointer - * to a #DoFHandler# object. - */ - MGDoFObjectAccessor (Triangulation *tria, - const int level, - const int index, - const AccessorData *local_data); - - /** - * Return the index of the #i#th degree - * of freedom of this hex on the level - * this quad lives on. - */ - int mg_dof_index (const unsigned int i) const; - - /** - * Set the index of the #i#th degree - * of freedom of this hex on the - * level this hex lives on to #index#. - */ - void set_mg_dof_index (const unsigned int i, const int index) const; - - /** - * Return the index of the #i#th degree - * on the #vertex#th vertex for the level - * this hex lives on. - */ - int mg_vertex_dof_index (const unsigned int vertex, - const unsigned int i) const; - - /** - * Set the index of the #i#th degree - * on the #vertex#th vertex for the - * level this hex lives on to #index#. - */ - void set_mg_vertex_dof_index (const unsigned int vertex, - const unsigned int i, - const int index) const; - - /** - * Return the indices of the dofs of this - * hex in the standard ordering: dofs - * on vertex 0, dofs on vertex 1, etc, - * dofs on line 0, dofs on line 1, etc, - * dofs on quad 0, etc. - * - * It is assumed that the vector already - * has the right size beforehand. The - * indices refer to the local numbering - * for the level this hex lives on. - */ - void get_mg_dof_indices (vector &dof_indices) const; - - /** - * Return the value of the given vector - * restricted to the dofs of this - * cell in the standard ordering: dofs - * on vertex 0, dofs on vertex 1, etc, - * dofs on line 0, dofs on line 1, etc, - * dofs on quad 0, etc. - * - * It is assumed that the vector already - * has the right size beforehand. The - * indices refer to the multilevel - * numbering local to the present - * level of this cell. The vector shall - * therefore have the same number of - * entries as there are degrees of - * freedom on this level. - */ - template - void get_mg_dof_values (const Vector &values, - Vector &dof_values) const; - - /** - * Return a pointer to the #i#th line - * bounding this #Hex#. - */ - TriaIterator > - line (const unsigned int i) const; - - /** - * Return a pointer to the #i#th quad - * bounding this #Hex#. - */ - TriaIterator > - quad (const unsigned int i) const; - - /** - * Return the #i#th child as a DoF quad - * iterator. This function is needed since - * the child function of the base - * class returns a hex accessor without - * access to the DoF data. - */ - TriaIterator > child (const unsigned int) const; - - /** - * Implement the copy operator needed - * for the iterator classes. - */ - void copy_from (const MGDoFObjectAccessor<3, dim> &a); -}; - - - -/** - * Grant access to the degrees of freedom on a cell. In fact, since all - * access to the degrees of freedom has been enabled by the classes - * #DoFObjectAccessor<1, 1># and #DoFObjectAccessor<2, 2># for the space dimension - * one and two, respectively, this class only collects the pieces - * together by deriving from the appropriate #DoF*Accessor# and the - * right #CellAccessor# and finally adding two functions which give - * access to the neighbors and children as #DoFCellAccessor# objects - * rather than #CellAccessor# objects (the latter function was inherited - * from the #CellAccessor# class). - * - * @author Wolfgang Bangerth, 1998 - */ -template -class MGDoFCellAccessor : public MGDoFObjectAccessor { - public: - /** - * Type of faces. - */ - typedef - TriaIterator > - face_iterator; - /** - * Declare the data type that this accessor - * class expects to get passed from the - * iterator classes. - */ - typedef typename MGDoFObjectAccessor::AccessorData AccessorData; - - /** - * Constructor - */ - MGDoFCellAccessor (Triangulation *tria, - const int level, - const int index, - const AccessorData *local_data) : - MGDoFObjectAccessor (tria,level,index,local_data) {}; - - /** - * Return the #i#th neighbor as a MGDoF cell - * iterator. This function is needed since - * the neighbor function of the base - * class returns a cell accessor without - * access to the MGDoF data. - */ - TriaIterator > neighbor (const unsigned int) const; - - /** - * Return the #i#th child as a MGDoF cell - * iterator. This function is needed since - * the child function of the base - * class returns a cell accessor without - * access to the DoF data. - */ - TriaIterator > child (const unsigned int) const; - - /** - * Return an iterator to the #i#th face - * of this cell. - * - * This function is not implemented in 1D, - * and maps to MGDoFObjectAccessor<2, dim>::line in 2D. - */ - face_iterator - face (const unsigned int i) const; - - /** - * Exception - */ - DeclException0 (ExcNotUsefulForThisDimension); -}; - - - - - - - -/*---------------------------- mg_dof_accessor.h ---------------------------*/ -/* end of #ifndef __mg_dof_accessor_H */ -#endif -/*---------------------------- mg_dof_accessor.h ---------------------------*/ diff --git a/deal.II/deal.II/include/dofs/mg_dof_handler.h b/deal.II/deal.II/include/dofs/mg_dof_handler.h deleted file mode 100644 index 5dda5971b8..0000000000 --- a/deal.II/deal.II/include/dofs/mg_dof_handler.h +++ /dev/null @@ -1,1009 +0,0 @@ -/*---------------------------- mg_dof.h ---------------------------*/ -/* $Id$ */ -/* Copyright W. Bangerth, University of Heidelberg, 1998 */ -#ifndef __mg_dof_H -#define __mg_dof_H -/*---------------------------- mg_dof.h ---------------------------*/ - - -#include - - - -/** - * Define some types which differ between the dimensions. This class - * is analogous to the \Ref{TriaDimensionInfo} class hierarchy. - * - * @see MGDoFDimensionInfo<1> - * @see MGDoFDimensionInfo<2> - */ -template -class MGDoFDimensionInfo; - - - - - -/** - * Define some types for the DoF handling in one dimension. - * - * The types have the same meaning as those declared in \Ref{TriaDimensionInfo<2>}. - */ -class MGDoFDimensionInfo<1> { - public: - typedef TriaRawIterator<1,MGDoFCellAccessor<1> > raw_line_iterator; - typedef TriaIterator<1,MGDoFCellAccessor<1> > line_iterator; - typedef TriaActiveIterator<1,MGDoFCellAccessor<1> > active_line_iterator; - - typedef void * raw_quad_iterator; - typedef void * quad_iterator; - typedef void * active_quad_iterator; - - typedef void * raw_hex_iterator; - typedef void * hex_iterator; - typedef void * active_hex_iterator; - - typedef raw_line_iterator raw_cell_iterator; - typedef line_iterator cell_iterator; - typedef active_line_iterator active_cell_iterator; - - typedef void * raw_face_iterator; - typedef void * face_iterator; - typedef void * active_face_iterator; -}; - - - - - - -/** - * Define some types for the DoF handling in two dimensions. - * - * The types have the same meaning as those declared in \Ref{TriaDimensionInfo<2>}. - */ -class MGDoFDimensionInfo<2> { - public: - typedef TriaRawIterator<2,MGDoFObjectAccessor<1, 2> > raw_line_iterator; - typedef TriaIterator<2,MGDoFObjectAccessor<1, 2> > line_iterator; - typedef TriaActiveIterator<2,MGDoFObjectAccessor<1, 2> > active_line_iterator; - - typedef TriaRawIterator<2,MGDoFCellAccessor<2> > raw_quad_iterator; - typedef TriaIterator<2,MGDoFCellAccessor<2> > quad_iterator; - typedef TriaActiveIterator<2,MGDoFCellAccessor<2> > active_quad_iterator; - - typedef void * raw_hex_iterator; - typedef void * hex_iterator; - typedef void * active_hex_iterator; - - typedef raw_quad_iterator raw_cell_iterator; - typedef quad_iterator cell_iterator; - typedef active_quad_iterator active_cell_iterator; - - typedef raw_line_iterator raw_face_iterator; - typedef line_iterator face_iterator; - typedef active_line_iterator active_face_iterator; -}; - - - - -/** - * Define some types for the DoF handling in two dimensions. - * - * The types have the same meaning as those declared in \Ref{TriaDimensionInfo<2>}. - */ -class MGDoFDimensionInfo<3> { - public: - typedef TriaRawIterator<3,MGDoFObjectAccessor<1, 3> > raw_line_iterator; - typedef TriaIterator<3,MGDoFObjectAccessor<1, 3> > line_iterator; - typedef TriaActiveIterator<3,MGDoFObjectAccessor<1, 3> > active_line_iterator; - - typedef TriaRawIterator<3,MGDoFObjectAccessor<2, 3> > raw_quad_iterator; - typedef TriaIterator<3,MGDoFObjectAccessor<2, 3> > quad_iterator; - typedef TriaActiveIterator<3,MGDoFObjectAccessor<2, 3> > active_quad_iterator; - - typedef TriaRawIterator<3,MGDoFCellAccessor<3> > raw_hex_iterator; - typedef TriaIterator<3,MGDoFCellAccessor<3> > hex_iterator; - typedef TriaActiveIterator<3,MGDoFCellAccessor<3> > active_hex_iterator; - - typedef raw_hex_iterator raw_cell_iterator; - typedef hex_iterator cell_iterator; - typedef active_hex_iterator active_cell_iterator; - - typedef raw_quad_iterator raw_face_iterator; - typedef quad_iterator face_iterator; - typedef active_quad_iterator active_face_iterator; -}; - - - - -/** - * This class manages degrees of freedom for a multilevel hierarchy of - * grids. It does mostly the same as does the #DoDHandler# class, but - * it uses a separate enumeration of the degrees of freedom on each - * level. For example, a vertex has several DoF numbers, one for each - * level of the triangulation on which it exists. - * - * At present, multilevel algorithms are not fully functional, so this - * documentation is still very brief. -//TODO: Extend MGDoFHandler doc - * - * @author Wolfgang Bangerth, 1998, 1999 - */ -template -class MGDoFHandler : public DoFHandler -{ - public: - typedef typename MGDoFDimensionInfo::raw_line_iterator raw_line_iterator; - typedef typename MGDoFDimensionInfo::line_iterator line_iterator; - typedef typename MGDoFDimensionInfo::active_line_iterator active_line_iterator; - - typedef typename MGDoFDimensionInfo::raw_quad_iterator raw_quad_iterator; - typedef typename MGDoFDimensionInfo::quad_iterator quad_iterator; - typedef typename MGDoFDimensionInfo::active_quad_iterator active_quad_iterator; - - typedef typename MGDoFDimensionInfo::raw_hex_iterator raw_hex_iterator; - typedef typename MGDoFDimensionInfo::hex_iterator hex_iterator; - typedef typename MGDoFDimensionInfo::active_hex_iterator active_hex_iterator; - - typedef typename MGDoFDimensionInfo::raw_cell_iterator raw_cell_iterator; - typedef typename MGDoFDimensionInfo::cell_iterator cell_iterator; - typedef typename MGDoFDimensionInfo::active_cell_iterator active_cell_iterator; - - typedef typename MGDoFDimensionInfo::raw_face_iterator raw_face_iterator; - typedef typename MGDoFDimensionInfo::face_iterator face_iterator; - typedef typename MGDoFDimensionInfo::active_face_iterator active_face_iterator; - - /** - * Constructor. Take #tria# as the - * triangulation to work on. - */ - MGDoFHandler (Triangulation &tria); - - /** - * Destructor - */ - virtual ~MGDoFHandler (); - - /** - * Go through the triangulation and - * distribute the degrees of freedoms - * needed for the given finite element - * according to the given distribution - * method. We first call the #DoFHandler#'s - * function and then distribute the - * levelwise numbers. - * - * A copy of the transferred finite - * element is stored. - * - * This function uses the user flags of the - * triangulation object, so make sure you - * don't need them after calling this - * function, or if so store them. - */ - virtual void distribute_dofs (const FiniteElement &, unsigned int offset=0); - - /** - * Actually do the renumbering based on - * a list of new dof numbers for all the - * dofs. - * - * #new_numbers# is an array of integers - * with size equal to the number of dofs - * on the present level. It stores the new - * indices after renumbering in the - * order of the old indices. - */ - void renumber_dofs (const unsigned int level, - const vector &new_numbers); - - /*--------------------------------------*/ - - /** - * @name Cell iterator functions - */ - /*@{*/ - /** - * Iterator to the first cell, used - * or not, on level #level#. If a level - * has no cells, a past-the-end iterator - * is returned. - * - * This function calls #begin_raw_line# - * in 1D and #begin_raw_quad# in 2D. - */ - raw_cell_iterator begin_raw (const unsigned int level = 0) const; - - /** - * Iterator to the first used cell - * on level #level#. - * - * This function calls #begin_line# - * in 1D and #begin_quad# in 2D. - */ - cell_iterator begin (const unsigned int level = 0) const; - - /** - * Iterator to the first active - * cell on level #level#. - * - * This function calls #begin_active_line# - * in 1D and #begin_active_quad# in 2D. - */ - active_cell_iterator begin_active(const unsigned int level = 0) const; - - /** - * Iterator past the end; this - * iterator serves for comparisons of - * iterators with past-the-end or - * before-the-beginning states. - * - * This function calls #end_line# - * in 1D and #end_quad# in 2D. - */ - raw_cell_iterator end () const; - - /** - * Return an iterator which is the first - * iterator not on level. If #level# is - * the last level, then this returns - * #end()#. - */ - cell_iterator end (const unsigned int level) const; - - /** - * Return a raw iterator which is the first - * iterator not on level. If #level# is - * the last level, then this returns - * #end()#. - */ - raw_cell_iterator end_raw (const unsigned int level) const; - - /** - * Return an active iterator which is the - * first iterator not on level. If #level# - * is the last level, then this returns - * #end()#. - */ - active_cell_iterator end_active (const unsigned int level) const; - - /** - * Return an iterator pointing to the - * last cell, used or not. - * - * This function calls #last_raw_line# - * in 1D and #last_raw_quad# in 2D. - */ - raw_cell_iterator last_raw () const; - - /** - * Return an iterator pointing to the last - * cell of the level #level#, used or not. - * - * This function calls #last_raw_line# - * in 1D and #last_raw_quad# in 2D. - */ - raw_cell_iterator last_raw (const unsigned int level) const; - - /** - * Return an iterator pointing to the last - * used cell. - * - * This function calls #last_line# - * in 1D and #last_quad# in 2D. - */ - cell_iterator last () const; - - /** - * Return an iterator pointing to the last - * used cell on level #level#. - * - * This function calls #last_line# - * in 1D and #last_quad# in 2D. - */ - cell_iterator last (const unsigned int level) const; - - /** - * Return an iterator pointing to the last - * active cell. - * - * This function calls #last_active_line# - * in 1D and #last_active_quad# in 2D. - */ - active_cell_iterator last_active () const; - - /** - * Return an iterator pointing to the last - * active cell on level #level#. - * - * This function calls #last_active_line# - * in 1D and #last_active_quad# in 2D. - */ - active_cell_iterator last_active (const unsigned int level) const; - //@} - - /*---------------------------------------*/ - - /** - * @name Face iterator functions - */ - /*@{*/ - /** - * Iterator to the first face, used - * or not, on level #level#. If a level - * has no faces, a past-the-end iterator - * is returned. - * - * This function calls #begin_raw_line# - * in 2D and #begin_raw_quad# in 3D. - */ - raw_face_iterator begin_raw_face (const unsigned int level = 0) const; - - /** - * Iterator to the first used face - * on level #level#. - * - * This function calls #begin_line# - * in 2D and #begin_quad# in 3D. - */ - face_iterator begin_face (const unsigned int level = 0) const; - - /** - * Iterator to the first active - * face on level #level#. - * - * This function calls #begin_active_line# - * in 2D and #begin_active_quad# in 3D. - */ - active_face_iterator begin_active_face(const unsigned int level = 0) const; - - /** - * Iterator past the end; this - * iterator serves for comparisons of - * iterators with past-the-end or - * before-the-beginning states. - * - * This function calls #end_line# - * in 2D and #end_quad# in 3D. - */ - raw_face_iterator end_face () const; - - /** - * Return an iterator which is the first - * iterator not on level. If #level# is - * the last level, then this returns - * #end()#. - */ - face_iterator end_face (const unsigned int level) const; - - /** - * Return a raw iterator which is the first - * iterator not on level. If #level# is - * the last level, then this returns - * #end()#. - */ - raw_face_iterator end_raw_face (const unsigned int level) const; - - /** - * Return an active iterator which is the - * first iterator not on level. If #level# - * is the last level, then this returns - * #end()#. - */ - active_face_iterator end_active_face (const unsigned int level) const; - - /** - * Return an iterator pointing to the - * last face, used or not. - * - * This function calls #last_raw_line# - * in 2D and #last_raw_quad# in 3D. - */ - raw_face_iterator last_raw_face () const; - - /** - * Return an iterator pointing to the last - * face of the level #level#, used or not. - * - * This function calls #last_raw_line# - * in 2D and #last_raw_quad# in 3D. - */ - raw_face_iterator last_raw_face (const unsigned int level) const; - - /** - * Return an iterator pointing to the last - * used face. - * - * This function calls #last_line# - * in 2D and #last_quad# in 3D. - */ - face_iterator last_face () const; - - /** - * Return an iterator pointing to the last - * used face on level #level#. - * - * This function calls #last_line# - * in 2D and #last_quad# in 3D. - */ - face_iterator last_face (const unsigned int level) const; - - /** - * Return an iterator pointing to the last - * active face. - * - * This function calls #last_active_line# - * in 2D and #last_active_quad# in 3D. - */ - active_face_iterator last_active_face () const; - - /** - * Return an iterator pointing to the last - * active face on level #level#. - * - * This function calls #last_active_line# - * in 2D and #last_active_quad# in 3D. - */ - active_face_iterator last_active_face (const unsigned int level) const; - //@} - - - /*---------------------------------------*/ - - /** - * @name Line iterator functions - */ - /*@{*/ - /** - * Iterator to the first line, used - * or not, on level #level#. If a level - * has no lines, a past-the-end iterator - * is returned. - */ - raw_line_iterator begin_raw_line (const unsigned int level = 0) const; - - /** - * Iterator to the first used line - * on level #level#. - */ - line_iterator begin_line (const unsigned int level = 0) const; - - /** - * Iterator to the first active - * line on level #level#. - */ - active_line_iterator begin_active_line(const unsigned int level = 0) const; - - /** - * Iterator past the end; this - * iterator serves for comparisons of - * iterators with past-the-end or - * before-the-beginning states. - */ - raw_line_iterator end_line () const; - - /** - * Return an iterator which is the first - * iterator not on level. If #level# is - * the last level, then this returns - * #end()#. - */ - line_iterator end_line (const unsigned int level) const; - - /** - * Return a raw iterator which is the first - * iterator not on level. If #level# is - * the last level, then this returns - * #end()#. - */ - raw_line_iterator end_raw_line (const unsigned int level) const; - - /** - * Return an active iterator which is the - * first iterator not on level. If #level# - * is the last level, then this returns - * #end()#. - */ - active_line_iterator end_active_line (const unsigned int level) const; - - - /** - * Return an iterator pointing to the - * last line, used or not. - */ - raw_line_iterator last_raw_line () const; - - /** - * Return an iterator pointing to the last - * line of the level #level#, used or not. - - */ - raw_line_iterator last_raw_line (const unsigned int level) const; - - /** - * Return an iterator pointing to the last - * used line. - */ - line_iterator last_line () const; - - /** - * Return an iterator pointing to the last - * used line on level #level#. - */ - line_iterator last_line (const unsigned int level) const; - - /** - * Return an iterator pointing to the last - * active line. - */ - active_line_iterator last_active_line () const; - - /** - * Return an iterator pointing to the last - * active line on level #level#. - */ - active_line_iterator last_active_line (const unsigned int level) const; - /*@}*/ - - /*---------------------------------------*/ - - /** - * @name Quad iterator functions*/ - /*@{ - */ - /** - * Iterator to the first quad, used - * or not, on level #level#. If a level - * has no quads, a past-the-end iterator - * is returned. - */ - raw_quad_iterator begin_raw_quad (const unsigned int level = 0) const; - - /** - * Iterator to the first used quad - * on level #level#. - */ - quad_iterator begin_quad (const unsigned int level = 0) const; - - /** - * Iterator to the first active - * quad on level #level#. - */ - active_quad_iterator begin_active_quad(const unsigned int level = 0) const; - - /** - * Iterator past the end; this - * iterator serves for comparisons of - * iterators with past-the-end or - * before-the-beginning states. - */ - raw_quad_iterator end_quad () const; - - /** - * Return an iterator which is the first - * iterator not on level. If #level# is - * the last level, then this returns - * #end()#. - */ - quad_iterator end_quad (const unsigned int level) const; - - /** - * Return a raw iterator which is the first - * iterator not on level. If #level# is - * the last level, then this returns - * #end()#. - */ - raw_quad_iterator end_raw_quad (const unsigned int level) const; - - /** - * Return an active iterator which is the - * first iterator not on level. If #level# - * is the last level, then this returns - * #end()#. - */ - active_quad_iterator end_active_quad (const unsigned int level) const; - - - /** - * Return an iterator pointing to the - * last quad, used or not. - */ - raw_quad_iterator last_raw_quad () const; - - /** - * Return an iterator pointing to the last - * quad of the level #level#, used or not. - - */ - raw_quad_iterator last_raw_quad (const unsigned int level) const; - - /** - * Return an iterator pointing to the last - * used quad. - */ - quad_iterator last_quad () const; - - /** - * Return an iterator pointing to the last - * used quad on level #level#. - */ - quad_iterator last_quad (const unsigned int level) const; - - /** - * Return an iterator pointing to the last - * active quad. - */ - active_quad_iterator last_active_quad () const; - - /** - * Return an iterator pointing to the last - * active quad on level #level#. - */ - active_quad_iterator last_active_quad (const unsigned int level) const; - /*@}*/ - - /*---------------------------------------*/ - - /** - * @name Hex iterator functions*/ - /*@{ - */ - /** - * Iterator to the first hex, used - * or not, on level #level#. If a level - * has no hexs, a past-the-end iterator - * is returned. - */ - raw_hex_iterator begin_raw_hex (const unsigned int level = 0) const; - - /** - * Iterator to the first used hex - * on level #level#. - */ - hex_iterator begin_hex (const unsigned int level = 0) const; - - /** - * Iterator to the first active - * hex on level #level#. - */ - active_hex_iterator begin_active_hex(const unsigned int level = 0) const; - - /** - * Iterator past the end; this - * iterator serves for comparisons of - * iterators with past-the-end or - * before-the-beginning states. - */ - raw_hex_iterator end_hex () const; - - /** - * Return an iterator which is the first - * iterator not on level. If #level# is - * the last level, then this returns - * #end()#. - */ - hex_iterator end_hex (const unsigned int level) const; - - /** - * Return a raw iterator which is the first - * iterator not on level. If #level# is - * the last level, then this returns - * #end()#. - */ - raw_hex_iterator end_raw_hex (const unsigned int level) const; - - /** - * Return an active iterator which is the - * first iterator not on level. If #level# - * is the last level, then this returns - * #end()#. - */ - active_hex_iterator end_active_hex (const unsigned int level) const; - - - /** - * Return an iterator pointing to the - * last hex, used or not. - */ - raw_hex_iterator last_raw_hex () const; - - /** - * Return an iterator pointing to the last - * hex of the level #level#, used or not. - - */ - raw_hex_iterator last_raw_hex (const unsigned int level) const; - - /** - * Return an iterator pointing to the last - * used hex. - */ - hex_iterator last_hex () const; - - /** - * Return an iterator pointing to the last - * used hex on level #level#. - */ - hex_iterator last_hex (const unsigned int level) const; - - /** - * Return an iterator pointing to the last - * active hex. - */ - active_hex_iterator last_active_hex () const; - - /** - * Return an iterator pointing to the last - * active hex on level #level#. - */ - active_hex_iterator last_active_hex (const unsigned int level) const; - /*@}*/ - - /*---------------------------------------*/ - - /** - * Return the number of degrees of freedom - * on the specified level. - * Included in this number are those - * DoFs which are constrained by - * hanging nodes. - */ - unsigned int n_dofs (const unsigned int level) const; - - /** - * Exception. - */ - DeclException1 (ExcInvalidLevel, - int, - << "The level index " << arg1 << "was not in the valid range."); - - private: - - /** We need for each vertex a list of the - * degree of freedom indices on each of the - * levels this vertex lives on. Since most - * vertices live only on a few levels, it - * is not economical to reserve indices for - * all the levels there are; rather, we - * create an object which holds the indices - * on those levels only where the vertex - * lives. To construct such an array, it - * is necessary to know beforehand which - * is the coarsest level the vertex lives - * on, how many levels it lives on and - * how many dofs there are on each vertex. - * If we have this information, we can - * allocate exactly the amount of memory - * which is needed and need not handle - * growing arrays and the like. - */ - class MGVertexDoFs { - public: - - /** - * Constructor. This one is empty - * because it is difficult to make it - * efficient to use vector<>'s and - * still construct the object using - * the constructor. Use the #init# - * function to really allocate - * memory. - */ - MGVertexDoFs (); - - /** - * Allocate memory and - * set all indices to #-1#. - */ - void init (const unsigned int coarsest_level, - const unsigned int finest_level, - const unsigned int dofs_per_vertex); - - /** - * Destructor - */ - ~MGVertexDoFs (); - - /** - * Set the index with number - * #dof_number# of this vertex on - * #level# to the given index. To - * compute the position in the array, - * one has to specify how many - * dofs per vertex there are. It is - * not checked that the level number - * is below the number of the finest - * level this vertex lives on. - * - * The function is inline, so should - * be reasonably fast. - */ - void set_index (const unsigned int level, - const unsigned int dof_number, - const unsigned int dofs_per_vertex, - const unsigned int index); - - /** - * Return the index with number - * #dof_number# of this vertex on - * #level#. To - * compute the position in the array, - * one has to specify how many - * dofs per vertex there are. It is - * not checked that the level number - * is below the number of the finest - * level this vertex lives on. - * - * The function is inline, so should - * be reasonably fast. - */ - int get_index (const unsigned int level, - const unsigned int dof_number, - const unsigned int dofs_per_vertex) const; - - /** - * Return the index of the coarsest - * level this vertex lives on. - */ - unsigned int get_coarsest_level () const; - - /** - * Return the index of the finest - * level this vertex lives on. - */ - unsigned int get_finest_level () const; - - /** - * Exception. - */ - DeclException0 (ExcNoMemory); - /** - * Exception. - */ - DeclException1 (ExcInvalidLevel, - int, - << "The given level number " << arg1 << " is outside " - << "the range of levels this vertex lives on."); - /** - * Exception. - */ - DeclException0 (ExcInternalError); - - private: - /** - * Store the coarsest level this - * vertex lives on. This is used - * as an offset when accessing the - * dofs of a specific level. - */ - unsigned int coarsest_level; - - /** - * Finest level this level lives on. - * This is mostly used for error - * checking but can also be accessed - * by the function #get_finest_level#. - */ - unsigned int finest_level; - - /** - * Array holding the indices. - */ - int *indices; - }; - - - /** - * Distribute dofs on the given cell, - * with new dofs starting with index - * #next_free_dof#. Return the next - * unused index number. The finite - * element used is the one given to - * #distribute_dofs#, which is copied - * to #selected_fe#. - * - * This function is excluded from the - * #distribute_dofs# function since - * it can not be implemented dimension - * independent. - * - * Note that unlike for the usual dofs, - * here all cells and not only active - * ones are allowed. - */ - unsigned int distribute_dofs_on_cell (cell_iterator &cell, - unsigned int next_free_dof); - - /** - * Reserve enough space for the MG dof - * indices for a given triangulation. - */ - void reserve_space (); - - /** - * Space to store the DoF numbers for the - * different levels. Unlike the #levels# - * object in the #DoFHandler#, these are - * not global numbers but rather are - * numbers which start from zero on each - * level. - */ - vector*> mg_levels; - - /** - * For each vertex there is a list of - * indices of the degrees of freedom indices - * on the different levels it lives on and - * which are these levels. - */ - vector mg_vertex_dofs; - - /** - * Vectors storing the number of degrees of - * freedom on each level. - */ - vector mg_used_dofs; - -#if (__GNUC__==2) && (__GNUC_MINOR__ < 95) - // this seems to be disallowed - // by the standard, so gcc2.95 - // does not accept it - friend class MGDoFObjectAccessor<1, dim>; - friend class MGDoFObjectAccessor<2, dim>; - friend class MGDoFObjectAccessor<3, dim>; -#else - // this, however, may grant - // access to too many classes, - // but ... - template friend class MGDoFObjectAccessor; -#endif -}; - - - - - - - -/* ----------------------- Inline functions of MGVertexDoFs -------------------*/ - - -template -inline -void MGDoFHandler::MGVertexDoFs::set_index (const unsigned int level, - const unsigned int dof_number, - const unsigned int dofs_per_vertex, - const unsigned int index) { - Assert ((level >= coarsest_level) && (level <= finest_level), - ExcInvalidLevel(level)); - Assert (dof_number < dofs_per_vertex, - ExcIndexRange(dof_number, 0, dofs_per_vertex)); - - indices[(level-coarsest_level)*dofs_per_vertex + dof_number] = index; -}; - - - - -template -inline -int MGDoFHandler::MGVertexDoFs::get_index (const unsigned int level, - const unsigned int dof_number, - const unsigned int dofs_per_vertex) const { - Assert ((level >= coarsest_level) && (level <= finest_level), - ExcInvalidLevel(level)); - Assert (dof_number < dofs_per_vertex, - ExcIndexRange (dof_number, 0, dofs_per_vertex)); - - return indices[(level-coarsest_level)*dofs_per_vertex + dof_number]; -}; - - - -/*---------------------------- mg_dof.h ---------------------------*/ -/* end of #ifndef __mg_dof_H */ -#endif -/*---------------------------- mg_dof.h ---------------------------*/ diff --git a/deal.II/deal.II/include/dofs/mg_dof_tools.h b/deal.II/deal.II/include/dofs/mg_dof_tools.h deleted file mode 100644 index c97e55e53c..0000000000 --- a/deal.II/deal.II/include/dofs/mg_dof_tools.h +++ /dev/null @@ -1,61 +0,0 @@ -/*---------------------------- mg_dof_tools.h ---------------------------*/ -/* $Id$ */ -#ifndef __mg_dof_tools_H -#define __mg_dof_tools_H -/*---------------------------- mg_dof_tools.h ---------------------------*/ - - -#include -#include - - -/** - * This is a collection of functions operating on, and manipulating - * the numbers of degrees of freedom in a multilevel triangulation. It - * is similar in purpose and function to the #DoFTools# class, but - * operates on #MGDoFHandler# objects instead of #DoFHandler# - * objects. See there and the documentation of the member functions - * for more information. - * - * All member functions are static, so there is no need to create an - * object of class #DoFTools#. - * - * @author Wolfgang Bangerth and others, 1999 - */ -class MGDoFTools -{ - public: - /** - * Write the sparsity structure - * of the matrix belonging to the - * specified #level# including - * constrained degrees of freedom - * into the matrix structure. The - * sparsity pattern does not - * include entries introduced by - * the elimination of constrained - * nodes. The sparsity pattern - * is not compressed, since if - * you want to call - * #ConstraintMatrix::condense(1)# - * afterwards, new entries have - * to be added. However, if you - * don't want to call - * #ConstraintMatrix::condense(1)#, - * you have to compress the - * matrix yourself, using - * #SparsityPattern::compress()#. - */ - template - static void - make_sparsity_pattern (const MGDoFHandler &dof_handler, - const unsigned int level, - SparsityPattern &sparsity); -}; - - - -/*---------------------------- mg_dof_tools.h ---------------------------*/ -/* end of #ifndef __mg_dof_tools_H */ -#endif -/*---------------------------- mg_dof_tools.h ---------------------------*/ diff --git a/deal.II/deal.II/include/multigrid/mg_base.h b/deal.II/deal.II/include/multigrid/mg_base.h index 3e92c8a903..496e9bd76d 100644 --- a/deal.II/deal.II/include/multigrid/mg_base.h +++ b/deal.II/deal.II/include/multigrid/mg_base.h @@ -1,20 +1,107 @@ -/*---------------------------- mg_base.h ---------------------------*/ +/*---------------------------- mgbase.h ---------------------------*/ /* $Id$ */ -#ifndef __mg_base_H -#define __mg_base_H -/*---------------------------- mg_base.h ---------------------------*/ +#ifndef __mgbase_H +#define __mgbase_H +/*---------------------------- mgbase.h ---------------------------*/ #include #include -#include #include +#include +#include + #include + +/** + * An array with an object for each level. The purpose of this class + * is mostly to allow access by level number, even if the lower levels + * are not used and therefore have no object at all; this is done by + * simply shifting the given index by the minimum level we have + * stored. + * + * In most cases, the objects which are stored on each levels, are + * either matrices or vectors. + * + * @author Wolfgang Bangerth, Guido Kanschat, 1999 + */ +template +class MGLevelObject : public Subscriptor +{ + public: + /** + * Constructor allowing to + * initialize the number of + * levels. By default, the object + * is created empty. + */ + MGLevelObject (const unsigned int minlevel = 0, + const unsigned int maxlevel = 0); + + /** + * Access object on level #level#. + */ + Object & operator[] (const unsigned int level); + + /** + * Access object on level + * #level#. Constant version. + */ + const Object & operator[] (const unsigned int level) const; + + /** + * Delete all previous contents + * of this object and reset its + * size according to the values + * of #new_minlevel# and + * #new_maxlevel#. + */ + void resize (const unsigned int new_minlevel, + const unsigned int new_maxlevel); + + /** + * Call #clear# on all objects + * stored by this object. This + * function is only implemented + * for some #Object# classes, + * most notably for vectors and + * matrices. Note that if + * #Object==Vector#, #clear# + * will set all entries to zero, + * while if + * #Object==std::vector#, + * #clear# deletes the elements + * of the vectors. This class + * might therefore not be useful + * for STL vectors. + */ + void clear(); + + ////// + unsigned int get_minlevel () const; + unsigned int get_maxlevel () const; + + private: + /** + * Level of first component. + */ + unsigned int minlevel; + + /** + * Array of the objects to be held. + */ + vector objects; +}; + + + + + /** * Abstract base class for coarse grid solvers. The interface of a * function call operator is defined to execute coarse grid solution @@ -50,78 +137,62 @@ class MGCoarseGridSolver : public Subscriptor - /** - * Abstract base class for multigrid smoothers. - * In fact, this class only provides the interface of the smoothing function. - * Using #deal.II# grid handling, #MGSmoother# is a good point to start. + * Coarse grid solver using LAC iterative methods. + * This is a little wrapper, transforming a triplet of iterative + * solver, matrix and preconditioner into a coarse grid solver. * - * @author Wolfgang Bangerth, Guido Kanschat, 1999 + * The type of the matrix (i.e. the template parameter #MATRIX#) + * should be derived from #Subscriptor# to allow for the use of a + * smart pointer to it. + * + * @author Guido Kanschat, 1999 */ -class MGSmootherBase : public Subscriptor -{ +template +class MGCoarseGridLACIteration : public MGCoarseGridSolver +{ public: /** - * Virtual destructor. Does - * nothing in particular, but - * needs to be declared anyhow. + * Constructor. + * Store solver, matrix and + * preconditioning method for later + * use. */ - virtual ~MGSmootherBase(); + MGCoarseGridLACIteration (SOLVER &, + const MATRIX &, + const PRECOND &); /** - * Smooth the residual of #u# on - * the given level. If $S$ is the - * smoothing operator, then this - * function should do the - * following operation: - * $u += S (rhs - Au)$, where #u# and - * #rhs# are the input parameters. - * - * This function should keep the - * interior level boundary - * values, so you may want to - * call #set_zero_interior_boundary# - * in the derived class - * #MGSmoother# somewhere in your - * derived function, or another - * function doing similar things. - */ - virtual void smooth (const unsigned int level, - Vector &u, - const Vector &rhs) const = 0; + * Implementation of the abstract + * function. + * Calls the solver method with + * matrix, vectors and + * preconditioner. + */ + virtual void operator() (const unsigned int level, + Vector &dst, + const Vector &src) const; + private: + /** + * Reference to the solver. + */ + SOLVER& solver; + + /** + * Reference to the matrix. + */ + const SmartPointer matrix; + + /** + * Reference to the preconditioner. + */ + const PRECOND& precondition; }; -/** - * Smoother doing nothing. This class is not useful for many applications other - * than for testing some multigrid procedures. Also some applications might - * get convergence without smoothing and then this class brings you the - * cheapest possible multigrid. - * - * @author Guido Kanschat, 1999 - */ -class MGSmootherIdentity : public MGSmootherBase -{ - public: - /** - * Implementation of the - * interface in #MGSmootherBase#. - * This function does nothing, - * which by comparison with the - * definition of this function - * means that the the smoothing - * operator equals the null - * operator. - */ - virtual void smooth (const unsigned int level, - Vector &u, - const Vector &rhs) const; -}; - - /** @@ -190,154 +261,7 @@ class MGTransferBase : public Subscriptor -/** - * An array with a vector for each level. The purpose of this class - * is mostly to allow access by level number, even if the lower levels - * are not used and therefore have no vector at all; this is done by - * simply shifting the given index by the minimum level we have - * stored. - * - * @author Guido Kanschat, 1999 - */ -template > -class MGVector : public Subscriptor -{ - public: - /** - * Constructor allowing to - * initialize the number of - * levels. - */ - MGVector (const unsigned int minlevel, - const unsigned int maxlevel); - - /** - * Access vector on level #level#. - */ - VECTOR & operator[] (const unsigned int level); - - /** - * Access vector on level - * #level#. Constant version. - */ - const VECTOR & operator[] (const unsigned int level) const; - - /** - * Call #clear# on all vectors - * stored by this object. Note - * that if #VECTOR==Vector#, - * #clear# will set all entries - * to zero, while if - * #VECTOR==std::vector#, - * #clear# deletes the elements - * of the vectors. This class - * might therefore not be useful - * for STL vectors. - */ - void clear(); - - private: - /** - * Level of first component. - */ - const unsigned int minlevel; - - /** - * Array of the vectors to be held. - */ - vector vectors; -}; - - - -/** - * An array of matrices for each level. - * This class provides the functionality of #vector# combined - * with a subscriptor for smart pointers. - * @author Guido Kanschat, 1999 - */ -template > -class MGMatrix : public Subscriptor, - public vector -{ - public: - /** - * Constructor allowing to initialize the number of levels. - */ - MGMatrix(unsigned int minlevel, unsigned int maxlevel); - - /** - * Safe access operator. - */ - MATRIX& operator[](unsigned int); - - /** - * Safe access operator. - */ - const MATRIX& operator[](unsigned int) const; - - private: - /** - * Level of first component. - */ - unsigned int minlevel; -}; - - - - -/** - * Coarse grid solver using LAC iterative methods. - * This is a little wrapper, transforming a triplet of iterative - * solver, matrix and preconditioner into a coarse grid solver. - * - * The type of the matrix (i.e. the template parameter #MATRIX#) - * should be derived from #Subscriptor# to allow for the use of a - * smart pointer to it. - * - * @author Guido Kanschat, 1999 - */ -template -class MGCoarseGridLACIteration : public MGCoarseGridSolver -{ - public: - /** - * Constructor. - * Store solver, matrix and - * preconditioning method for later - * use. - */ - MGCoarseGridLACIteration (SOLVER &, - const MATRIX &, - const PRECOND &); - - /** - * Implementation of the abstract - * function. - * Calls the solver method with - * matrix, vectors and - * preconditioner. - */ - virtual void operator() (const unsigned int level, - Vector &dst, - const Vector &src) const; - private: - /** - * Reference to the solver. - */ - SOLVER& solver; - - /** - * Reference to the matrix. - */ - const SmartPointer matrix; - - /** - * Reference to the preconditioner. - */ - const PRECOND& precondition; -}; @@ -432,12 +356,12 @@ class MGBase : public Subscriptor /** * Auxiliary vector. */ - MGVector > defect; + MGLevelObject > defect; /** * Auxiliary vector. */ - MGVector > solution; + MGLevelObject > solution; /** * Prolongation and restriction object. @@ -549,74 +473,77 @@ class PreconditionMG /* ------------------------------------------------------------------- */ -template -MGVector::MGVector(const unsigned int min, - const unsigned int max) + +template +MGLevelObject::MGLevelObject(const unsigned int min, + const unsigned int max) : - minlevel(min), - vectors(max-min+1) -{}; + minlevel(0) +{ + resize (min, max); +}; -template -VECTOR & -MGVector::operator[] (const unsigned int i) +template +Object & +MGLevelObject::operator[] (const unsigned int i) { - Assert((i>=minlevel) && (i=minlevel) && (i -const VECTOR & -MGVector::operator[] (const unsigned int i) const +template +const Object & +MGLevelObject::operator[] (const unsigned int i) const { - Assert((i>=minlevel) && (i=minlevel) && (i + +template void -MGVector::clear() +MGLevelObject::resize (const unsigned int new_minlevel, + const unsigned int new_maxlevel) { - typename vector::iterator v; - for (v = vectors.begin(); v != vectors.end(); ++v) - v->clear(); - -} + Assert (new_minlevel <= new_maxlevel, ExcInternalError()); + objects.clear (); + minlevel = new_minlevel; + objects.resize (new_maxlevel - new_minlevel + 1); +}; -/* ------------------------------------------------------------------- */ -template -MGMatrix::MGMatrix(unsigned int min, unsigned int max) - : - vector(max-min+1), - minlevel(min) -{} +template +void +MGLevelObject::clear () +{ + typename vector::iterator v; + for (v = objects.begin(); v != objects.end(); ++v) + v->clear(); +}; -template -MATRIX& -MGMatrix::operator[](unsigned int i) +template +unsigned int +MGLevelObject::get_minlevel () const { - Assert((i>=minlevel)&&(i::operator[](i-minlevel); + return minlevel; }; -template -const MATRIX& -MGMatrix::operator[](unsigned int i) const +template +unsigned int +MGLevelObject::get_maxlevel () const { - Assert((i>=minlevel)&&(i::operator[](i-minlevel); + return minlevel + objects.size() - 1; }; @@ -681,7 +608,7 @@ PreconditionMG::operator() (VECTOR &dst, -/*---------------------------- mg_base.h ---------------------------*/ -/* end of #ifndef __mg_base_H */ +/*---------------------------- mgbase.h ---------------------------*/ +/* end of #ifndef __mgbase_H */ #endif -/*---------------------------- mg_base.h ---------------------------*/ +/*---------------------------- mgbase.h ---------------------------*/ diff --git a/deal.II/deal.II/include/multigrid/mg_dof_tools.h b/deal.II/deal.II/include/multigrid/mg_dof_tools.h index c97e55e53c..ee586c0aa5 100644 --- a/deal.II/deal.II/include/multigrid/mg_dof_tools.h +++ b/deal.II/deal.II/include/multigrid/mg_dof_tools.h @@ -44,7 +44,7 @@ class MGDoFTools * #ConstraintMatrix::condense(1)#, * you have to compress the * matrix yourself, using - * #SparsityPattern::compress()#. + * #SparseMatrixStruct::compress()#. */ template static void diff --git a/deal.II/deal.II/include/multigrid/mg_smoother.h b/deal.II/deal.II/include/multigrid/mg_smoother.h index 0c73624f05..dce540bdca 100644 --- a/deal.II/deal.II/include/multigrid/mg_smoother.h +++ b/deal.II/deal.II/include/multigrid/mg_smoother.h @@ -6,12 +6,85 @@ #include -#include +#include #include #include #include + +/** + * Abstract base class for multigrid smoothers. + * In fact, this class only provides the interface of the smoothing function. + * Using #deal.II# grid handling, #MGSmoother# is a good point to start. + * + * @author Wolfgang Bangerth, Guido Kanschat, 1999 + */ +class MGSmootherBase : public Subscriptor +{ + public: + /** + * Virtual destructor. Does + * nothing in particular, but + * needs to be declared anyhow. + */ + virtual ~MGSmootherBase(); + + /** + * Smooth the residual of #u# on + * the given level. If $S$ is the + * smoothing operator, then this + * function should do the + * following operation: + * $u += S (rhs - Au)$, where #u# and + * #rhs# are the input parameters. + * + * This function should keep the + * interior level boundary + * values, so you may want to + * call #set_zero_interior_boundary# + * in the derived class + * #MGSmoother# somewhere in your + * derived function, or another + * function doing similar things. + */ + virtual void smooth (const unsigned int level, + Vector &u, + const Vector &rhs) const = 0; +}; + + + + + +/** + * Smoother doing nothing. This class is not useful for many applications other + * than for testing some multigrid procedures. Also some applications might + * get convergence without smoothing and then this class brings you the + * cheapest possible multigrid. + * + * @author Guido Kanschat, 1999 + */ +class MGSmootherIdentity : public MGSmootherBase +{ + public: + /** + * Implementation of the + * interface in #MGSmootherBase#. + * This function does nothing, + * which by comparison with the + * definition of this function + * means that the the smoothing + * operator equals the null + * operator. + */ + virtual void smooth (const unsigned int level, + Vector &u, + const Vector &rhs) const; +}; + + + /** * Base class for multigrid smoothers. It declares the interface * to smoothers by inheriting #MGSmootherBase# and implements some functionality for setting @@ -111,6 +184,9 @@ class MGSmoother : public MGSmootherBase vector > interior_boundary_dofs; }; + + + /** * Implementation of a smoother using matrix builtin relaxation methods. * @@ -146,11 +222,11 @@ class MGSmootherRelaxation : public MGSmoother */ template - MGSmootherRelaxation(const MGDoFHandler &mg_dof, - const MGMatrix >& matrix, - function_ptr function, - unsigned int steps, - double omega = 1.); + MGSmootherRelaxation(const MGDoFHandler &mg_dof, + const MGLevelObject > &matrix, + const function_ptr function, + const unsigned int steps, + const double omega = 1.); /** * Implementation of the interface in #MGSmootherBase#. @@ -164,25 +240,33 @@ class MGSmootherRelaxation : public MGSmoother /** * Pointer to the matrices. */ - SmartPointer< const MGMatrix< SparseMatrix > > matrix; + SmartPointer > > matrix; + /** * Pointer to the relaxation function. */ function_ptr relaxation; + /** * Relaxation parameter. */ double omega; + /** * Auxiliary vector. */ mutable Vector h1; + /** * Auxiliary vector. */ mutable Vector h2; }; + + +/* ------------------------------- Inline functions -------------------------- */ + inline void MGSmoother::set_steps(unsigned int i) diff --git a/deal.II/deal.II/include/multigrid/mg_smoother.templates.h b/deal.II/deal.II/include/multigrid/mg_smoother.templates.h index f72af10b02..4c7637a43a 100644 --- a/deal.II/deal.II/include/multigrid/mg_smoother.templates.h +++ b/deal.II/deal.II/include/multigrid/mg_smoother.templates.h @@ -1,25 +1,28 @@ // $Id$ + template template MGSmootherRelaxation:: -MGSmootherRelaxation(const MGDoFHandler &mg_dof, - const MGMatrix >& matrix, - function_ptr relaxation, - unsigned int steps, - double omega) +MGSmootherRelaxation (const MGDoFHandler &mg_dof, + const MGLevelObject > &matrix, + const function_ptr relaxation, + const unsigned int steps, + const double omega) : MGSmoother(mg_dof, steps), matrix(&matrix), relaxation(relaxation), omega(omega) -{} +{}; + + template void MGSmootherRelaxation::smooth (const unsigned int level, - Vector &u, - const Vector &rhs) const + Vector &u, + const Vector &rhs) const { h1.reinit(u); h2.reinit(u); diff --git a/deal.II/deal.II/include/multigrid/multigrid.h b/deal.II/deal.II/include/multigrid/multigrid.h index 16f83a9083..602e11bafa 100644 --- a/deal.II/deal.II/include/multigrid/multigrid.h +++ b/deal.II/deal.II/include/multigrid/multigrid.h @@ -8,8 +8,9 @@ #include #include #include +#include #include -#include +#include #include @@ -17,14 +18,14 @@ /** * Implementation of multigrid with matrices. - * While #MGBase# was only an abstract framework for the v-cycle, + * While #MGBase# was only an abstract framework for the V-cycle, * here we have the implementation of the pure virtual functions defined there. * Furthermore, level information is obtained from a triangulation object. * * @author Wolfgang Bangerth, Guido Kanschat, 1999 */ -template -class MG : public MGBase +template +class Multigrid : public MGBase { public: /** @@ -37,10 +38,12 @@ class MG : public MGBase * (for example for hanging * nodes), a set of matrices * which describe the equation - * discretized on each level, and - * an object mediating the - * transfer of solutions between - * different levels. + * discretized on each level, + * their respective sparsity + * patterns, and an object + * mediating the transfer of + * solutions between different + * levels. * * This function already * initializes the vectors which @@ -51,12 +54,13 @@ class MG : public MGBase * this type as late as possible. */ //TODO: minlevel, maxlevel? - MG(const MGDoFHandler &mg_dof_handler, - ConstraintMatrix &constraints, - const MGMatrix > &level_matrices, - const MGTransferBase &transfer, - const unsigned int minlevel = 0, - const unsigned int maxlevel = 1000000); + Multigrid(const MGDoFHandler &mg_dof_handler, + const ConstraintMatrix &constraints, + const MGLevelObject &level_sparsities, + const MGLevelObject > &level_matrices, + const MGTransferBase &transfer, + const unsigned int minlevel = 0, + const unsigned int maxlevel = 1000000); /** * Transfer from a vector on the @@ -76,7 +80,8 @@ class MG : public MGBase * the respective positions of a * #Vector#. In order to * keep the result consistent, - * constrained degrees of freedom are set to zero. + * constrained degrees of freedom + * are set to zero. */ template void copy_from_mg (Vector &dst) const; @@ -102,13 +107,20 @@ class MG : public MGBase */ SmartPointer > mg_dof_handler; + + /** + * Sparsity patterns for each level. + */ + SmartPointer > level_sparsities; + /** - * Matrices for each level. - * The matrices are prepared by - * the constructor of #MG# and can - * be accessed for assembling. + * Matrices for each level. The + * matrices are prepared by the + * constructor of #Multigrid# and + * can be accessed for + * assembling. */ - SmartPointer > > level_matrices; + SmartPointer > > level_matrices; /** * Pointer to the object holding @@ -193,21 +205,16 @@ class MGTransferPrebuilt : public MGTransferBase private: - /** - * Sparsity patterns for the - * prolongation matrices on each - * level. - */ - vector prolongation_sparsities; + vector prolongation_sparsities; - /** - * The actual prolongation matrix. - * column indices belong to the - * dof indices of the mother cell, - * i.e. the coarse level. - * while row indices belong to the - * child cell, i.e. the fine level. - */ + /** + * The actual prolongation matrix. + * column indices belong to the + * dof indices of the mother cell, + * i.e. the coarse level. + * while row indices belong to the + * child cell, i.e. the fine level. + */ vector > prolongation_matrices; }; @@ -220,7 +227,7 @@ class MGTransferPrebuilt : public MGTransferBase template inline const SparseMatrix& -MG::get_matrix(unsigned int level) const +Multigrid::get_matrix (const unsigned int level) const { Assert((level>=minlevel) && (level -#include +#include #include #include +#include +//TODO: This function needs to be specially implemented, since in 2d mode we use faces +#if deal_II_dimension != 1 template template void -MG::copy_to_mg(const Vector& src) +Multigrid::copy_to_mg(const Vector& src) { const unsigned int dofs_per_cell = mg_dof_handler->get_fe().dofs_per_cell; const unsigned int dofs_per_face = mg_dof_handler->get_fe().dofs_per_face; @@ -27,9 +30,19 @@ MG::copy_to_mg(const Vector& src) // constraints->condense(src); vector global_dof_indices (dofs_per_cell); - vector level_dof_indices (dofs_per_cell); + vector level_dof_indices (dofs_per_cell); vector level_face_indices (dofs_per_face); + // initialize the objects with + // their correct size + for (unsigned int level=minlevel; level<=maxlevel ; ++level) + { + const unsigned int system_size = (*level_matrices)[level].m(); + + solution[level].reinit(system_size); + defect[level].reinit(system_size); + }; + // traverse the grid top-down // (i.e. starting with the most // refined grid). this way, we can @@ -88,12 +101,13 @@ MG::copy_to_mg(const Vector& src) }; }; +#endif template template void -MG::copy_from_mg(Vector &dst) const +Multigrid::copy_from_mg(Vector &dst) const { const unsigned int dofs_per_cell = mg_dof_handler->get_fe().dofs_per_cell; diff --git a/deal.II/deal.II/include/numerics/mg_smoother.h b/deal.II/deal.II/include/numerics/mg_smoother.h deleted file mode 100644 index 0c73624f05..0000000000 --- a/deal.II/deal.II/include/numerics/mg_smoother.h +++ /dev/null @@ -1,204 +0,0 @@ -/*---------------------------- mg_smoother.h ---------------------------*/ -/* $Id$ */ -#ifndef __mg_smoother_H -#define __mg_smoother_H -/*---------------------------- mg_smoother.h ---------------------------*/ - - -#include -#include -#include -#include -#include - - -/** - * Base class for multigrid smoothers. It declares the interface - * to smoothers by inheriting #MGSmootherBase# and implements some functionality for setting - * the values - * of vectors at interior boundaries (i.e. boundaries between differing - * levels of the triangulation) to zero, by building a list of these degrees - * of freedom's indices at construction time. - * - * @author Wolfgang Bangerth, Guido Kanschat, 1999 - */ -class MGSmoother : public MGSmootherBase -{ - private: - /** - * Default constructor. Made private - * to prevent it being called, which - * is necessary since this could - * circumpass the set-up of the list - * if interior boundary dofs. - */ - MGSmoother (); - - public: - - /** - * Constructor. This one collects - * the indices of the degrees of freedom - * on the interior boundaries between - * the different levels, which are - * needed by the function - * #set_zero_interior_boundaries#. - * - * Since this function is implemented - * a bit different in 1d (there are no - * faces of cells, just vertices), - * there are actually two sets of - * constructors, namely this one for 1d - * and the following one for all other - * dimensions. - */ - MGSmoother (const MGDoFHandler<1> &mg_dof, - unsigned int steps); - - /** - * Constructor. This one collects - * the indices of the degrees of freedom - * on the interior boundaries between - * the different levels, which are - * needed by the function - * #set_zero_interior_boundaries#. - * - * The parameter steps indicates the number of smoothing - * steps to be executed by #smooth#. - */ - template - MGSmoother (const MGDoFHandler &mg_dof, - unsigned int steps); - - /** - * Reset the values of the degrees of - * freedom on interior boundaries between - * different levels to zero in the given - * data vector #u#. - * - * Since the coarsest level (#level==0#) - * has no interior boundaries, this - * function does nothing in this case. - */ - void set_zero_interior_boundary (const unsigned int level, - Vector &u) const; - /** - * Modify the number of smoothing steps. - */ - void set_steps(unsigned int steps); - /** - * How many steps should be used? - */ - unsigned int get_steps() const; - - private: - /** - * Number of smoothing steps. - */ - unsigned int steps; - /** - * For each level, we store a list of - * degree of freedom indices which are - * located on interior boundaries between - * differing levels of the triangulation. - * Since the coarsest level has no - * interior boundary dofs, the first - * entry refers to the second level. - * - * These arrays are set by the constructor. - * The entries for each level are sorted ascendingly. - */ - vector > interior_boundary_dofs; -}; - -/** - * Implementation of a smoother using matrix builtin relaxation methods. - * - * @author Wolfgang Bangerth, Guido Kanschat, 1999 - */ -template -class MGSmootherRelaxation : public MGSmoother -{ - public: - /** - * Type of the smoothing - * function of the matrix. - */ - typedef void - (SparseMatrix::* function_ptr)(Vector&, const Vector&, - typename SparseMatrix::value_type) const; - - /** - * Constructor. - * The constructor uses an #MGDoFHandler# to initialize - * data structures for interior level boundary handling - * in #MGSmoother#. - * - * Furthermore, it takes a pointer to the - * level matrices and their smoothing function. - * This function must perform one relaxation step - * like #SparseMatrix::SOR_step# does. Do not - * use the preconditioning methods because they apply - * a preconditioning matrix to the residual. - * - * The final two parameters are the number of relaxation - * steps and the relaxation parameter. - */ - - template - MGSmootherRelaxation(const MGDoFHandler &mg_dof, - const MGMatrix >& matrix, - function_ptr function, - unsigned int steps, - double omega = 1.); - - /** - * Implementation of the interface in #MGSmootherBase#. - * We use the SOR method in #SparseMatrix# for the real work - * and find a way to keep the boundary values. - */ - virtual void smooth (const unsigned int level, - Vector &u, - const Vector &rhs) const; - private: - /** - * Pointer to the matrices. - */ - SmartPointer< const MGMatrix< SparseMatrix > > matrix; - /** - * Pointer to the relaxation function. - */ - function_ptr relaxation; - /** - * Relaxation parameter. - */ - double omega; - /** - * Auxiliary vector. - */ - mutable Vector h1; - /** - * Auxiliary vector. - */ - mutable Vector h2; -}; - -inline -void -MGSmoother::set_steps(unsigned int i) -{ - steps = i; -} - -inline -unsigned int -MGSmoother::get_steps() const -{ - return steps; -} - -/*---------------------------- mg_smoother.h ---------------------------*/ -/* end of #ifndef __mg_smoother_H */ -#endif -/*---------------------------- mg_smoother.h ---------------------------*/ - diff --git a/deal.II/deal.II/include/numerics/mg_smoother.templates.h b/deal.II/deal.II/include/numerics/mg_smoother.templates.h deleted file mode 100644 index f72af10b02..0000000000 --- a/deal.II/deal.II/include/numerics/mg_smoother.templates.h +++ /dev/null @@ -1,34 +0,0 @@ -// $Id$ - -template -template -MGSmootherRelaxation:: -MGSmootherRelaxation(const MGDoFHandler &mg_dof, - const MGMatrix >& matrix, - function_ptr relaxation, - unsigned int steps, - double omega) - : - MGSmoother(mg_dof, steps), - matrix(&matrix), - relaxation(relaxation), - omega(omega) -{} - -template -void -MGSmootherRelaxation::smooth (const unsigned int level, - Vector &u, - const Vector &rhs) const -{ - h1.reinit(u); - h2.reinit(u); - for(unsigned i=0;i -#include -#include -#include -#include -#include - -#include - - - -/** - * Implementation of multigrid with matrices. - * While #MGBase# was only an abstract framework for the v-cycle, - * here we have the implementation of the pure virtual functions defined there. - * Furthermore, level information is obtained from a triangulation object. - * - * @author Wolfgang Bangerth, Guido Kanschat, 1999 - */ -template -class MG : public MGBase -{ - public: - /** - * Constructor. Take the object - * describing the multilevel - * hierarchy of degrees of - * freedom, an object describing - * constraints to degrees of - * freedom on the global grid - * (for example for hanging - * nodes), a set of matrices - * which describe the equation - * discretized on each level, and - * an object mediating the - * transfer of solutions between - * different levels. - * - * This function already - * initializes the vectors which - * will be used later on in the - * course of the - * computations. You should - * therefore create objects of - * this type as late as possible. - */ -//TODO: minlevel, maxlevel? - MG(const MGDoFHandler &mg_dof_handler, - ConstraintMatrix &constraints, - const MGMatrix > &level_matrices, - const MGTransferBase &transfer, - const unsigned int minlevel = 0, - const unsigned int maxlevel = 1000000); - - /** - * Transfer from a vector on the - * global grid to vectors defined - * on each of the levels - * separately, i.a. an #MGVector#. - */ - template - void copy_to_mg (const Vector &src); - - /** - * Transfer from multi-level vector to - * normal vector. - * - * Copies data from active - * portions of an MGVector into - * the respective positions of a - * #Vector#. In order to - * keep the result consistent, - * constrained degrees of freedom are set to zero. - */ - template - void copy_from_mg (Vector &dst) const; - - /** - * Negative #vmult# on a level. -//TODO: what does this function do? - * @see MGBase. - */ - virtual void level_vmult (const unsigned int level, - Vector &dest, - const Vector &src, - const Vector &rhs); - - /** - * Read-only access to level matrices. - */ - const SparseMatrix& get_matrix (const unsigned int level) const; - - private: - /** - * Associated #MGDoFHandler#. - */ - SmartPointer > mg_dof_handler; - - /** - * Matrices for each level. - * The matrices are prepared by - * the constructor of #MG# and can - * be accessed for assembling. - */ - SmartPointer > > level_matrices; - - /** - * Pointer to the object holding - * constraints. - */ - SmartPointer constraints; -}; - - - - - -/** - * Implementation of the #MGTransferBase# interface for which the transfer - * operations are prebuilt upon construction of the object of this class as - * matrices. This is the fast way, since it only needs to build the operation - * once by looping over all cells and storing the result in a matrix for - * each level, but requires additional memory. - * - * @author Wolfgang Bangerth, Guido Kanschat, 1999 - */ -class MGTransferPrebuilt : public MGTransferBase -{ - public: - /** - * Destructor. - */ - virtual ~MGTransferPrebuilt (); - - /** - * Actually build the prolongation - * matrices for each level. - */ - template - void build_matrices (const MGDoFHandler &mg_dof); - - /** - * Prolongate a vector from level - * #to_level-1# to level - * #to_level#. The previous - * content of #dst# is - * overwritten. - * - * #src# is assumed to be a vector with - * as many elements as there are degrees - * of freedom on the coarser level of - * the two involved levels, while #src# - * shall have as many elements as there - * are degrees of freedom on the finer - * level. - */ - virtual void prolongate (const unsigned int to_level, - Vector &dst, - const Vector &src) const; - - /** - * Restrict a vector from level - * #from_level# to level - * #from_level-1# and add this - * restriction to - * #dst#. Obviously, if the - * refined region on level - * #from_level# is smaller than - * that on level #from_level-1#, - * some degrees of freedom in - * #dst# are not covered and will - * not be altered. For the other - * degress of freedom, the result - * of the restriction is added. - * - * #src# is assumed to be a vector with - * as many elements as there are degrees - * of freedom on the finer level of - * the two involved levels, while #src# - * shall have as many elements as there - * are degrees of freedom on the coarser - * level. - */ - virtual void restrict_and_add (const unsigned int from_level, - Vector &dst, - const Vector &src) const; - - private: - - /** - * Sparsity patterns for the - * prolongation matrices on each - * level. - */ - vector prolongation_sparsities; - - /** - * The actual prolongation matrix. - * column indices belong to the - * dof indices of the mother cell, - * i.e. the coarse level. - * while row indices belong to the - * child cell, i.e. the fine level. - */ - vector > prolongation_matrices; -}; - - - - -/* --------------------------- inline functions --------------------- */ - - -template -inline -const SparseMatrix& -MG::get_matrix(unsigned int level) const -{ - Assert((level>=minlevel) && (level -#include -#include -#include - - - - -template -template -void -MG::copy_to_mg(const Vector& src) -{ - const unsigned int dofs_per_cell = mg_dof_handler->get_fe().dofs_per_cell; - const unsigned int dofs_per_face = mg_dof_handler->get_fe().dofs_per_face; - - // set the elements of the vectors - // on all levels to zero - defect.clear(); -// constraints->condense(src); - - vector global_dof_indices (dofs_per_cell); - vector level_dof_indices (dofs_per_cell); - vector level_face_indices (dofs_per_face); - - // traverse the grid top-down - // (i.e. starting with the most - // refined grid). this way, we can - // always get that part of one - // level of the output vector which - // corresponds to a region which is - // more refined, by restriction of - // the respective vector on the - // next finer level, which we then - // already have built. - for (int level=maxlevel; level>=static_cast(minlevel); --level) - { - DoFHandler::active_cell_iterator - global_cell = mg_dof_handler->DoFHandler::begin_active(level); - MGDoFHandler::active_cell_iterator - level_cell = mg_dof_handler->begin_active(level); - const MGDoFHandler::active_cell_iterator - level_end = mg_dof_handler->end_active(level); - - for (; level_cell!=level_end; ++level_cell, ++global_cell) - { - // get the dof numbers of - // this cell for the global - // and the level-wise - // numbering - global_cell->get_dof_indices(global_dof_indices); - level_cell->get_mg_dof_indices (level_dof_indices); - - // transfer the global - // defect in the vector - // into the level-wise one - for (unsigned int i=0; i::faces_per_cell; ++face_n) - { - const MGDoFHandler::face_iterator face = level_cell->face(face_n); - if (face->has_children()) - { - face->get_mg_dof_indices(level_face_indices); - for (unsigned int i=0; i(level) < maxlevel) - transfer->restrict_and_add (level+1, defect[level], defect[level+1]); - }; -}; - - - -template -template -void -MG::copy_from_mg(Vector &dst) const -{ - const unsigned int dofs_per_cell = mg_dof_handler->get_fe().dofs_per_cell; - - vector global_dof_indices (dofs_per_cell); - vector level_dof_indices (dofs_per_cell); - - DoFHandler::active_cell_iterator - global_cell = mg_dof_handler->DoFHandler::begin_active(); - MGDoFHandler::active_cell_iterator - level_cell = mg_dof_handler->begin_active(); - const MGDoFHandler::active_cell_iterator - endc = mg_dof_handler->end(); - - // traverse all cells and copy the - // data appropriately to the output - // vector - for (; level_cell != endc; ++level_cell, ++global_cell) - { - const unsigned int level = level_cell->level(); - - // get the dof numbers of - // this cell for the global - // and the level-wise - // numbering - global_cell->get_dof_indices (global_dof_indices); - level_cell->get_mg_dof_indices(level_dof_indices); - - // copy level-wise data to - // global vector - for (unsigned int i=0; iset_zero(dst); -} - - - - -/*---------------------------- multigrid.templates.h ---------------------------*/ -/* end of #ifndef __multigrid.templates_H */ -#endif -/*---------------------------- multigrid.templates.h ---------------------------*/ diff --git a/deal.II/deal.II/source/dofs/dof_renumbering.cc b/deal.II/deal.II/source/dofs/dof_renumbering.cc index 576a66fa1f..2e95645a98 100644 --- a/deal.II/deal.II/source/dofs/dof_renumbering.cc +++ b/deal.II/deal.II/source/dofs/dof_renumbering.cc @@ -4,10 +4,10 @@ #include #include #include -#include #include #include -#include +#include +#include #include #include diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc index 95faa5439e..05268f83fe 100644 --- a/deal.II/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/deal.II/source/dofs/dof_tools.cc @@ -4,10 +4,10 @@ #include #include #include -#include #include -#include #include +#include +#include #include #include #include diff --git a/deal.II/deal.II/source/dofs/mg_dof_accessor.cc b/deal.II/deal.II/source/dofs/mg_dof_accessor.cc deleted file mode 100644 index 9b0b78cf2e..0000000000 --- a/deal.II/deal.II/source/dofs/mg_dof_accessor.cc +++ /dev/null @@ -1,776 +0,0 @@ -/* $Id$ */ -/* Copyright W. Bangerth, University of Heidelberg, 1998 */ - - - -#include -#include -#include -#include -#include -#include - -#include - - - - - - - - -/* ------------------------ MGDoFLineAccessor --------------------------- */ - -template -MGDoFObjectAccessor<1, dim>::MGDoFObjectAccessor (Triangulation *tria, - const int level, - const int index, - const AccessorData *local_data) : - MGDoFAccessor (local_data), - MGDoFObjectAccessor_Inheritance<1,dim>::BaseClass(tria,level,index,local_data) -{}; - - - -template -int MGDoFObjectAccessor<1, dim>::mg_dof_index (const unsigned int i) const -{ - Assert (dof_handler != 0, DoFAccessor::ExcInvalidObject()); - Assert (mg_dof_handler != 0, DoFAccessor::ExcInvalidObject()); - // make sure a FE has been selected - // and enough room was reserved - Assert (&dof_handler->get_fe() != 0, DoFAccessor::ExcInvalidObject()); - Assert (iget_fe().dofs_per_line, - ExcIndexRange (i, 0, dof_handler->get_fe().dofs_per_line)); - - return mg_dof_handler->mg_levels[present_level] - ->line_dofs[present_index*dof_handler->get_fe().dofs_per_line+i]; -}; - - - - - -template -void MGDoFObjectAccessor<1, dim>::set_mg_dof_index (const unsigned int i, - const int index) const -{ - Assert (dof_handler != 0, DoFAccessor::ExcInvalidObject()); - Assert (mg_dof_handler != 0, DoFAccessor::ExcInvalidObject()); - // make sure a FE has been selected - // and enough room was reserved - Assert (&dof_handler->get_fe() != 0, DoFAccessor::ExcInvalidObject()); - Assert (iget_fe().dofs_per_line, - ExcIndexRange (i, 0, dof_handler->get_fe().dofs_per_line)); - - mg_dof_handler->mg_levels[present_level] - ->line_dofs[present_index*dof_handler->get_fe().dofs_per_line+i] = index; -}; - - - - -template -int MGDoFObjectAccessor<1, dim>::mg_vertex_dof_index (const unsigned int vertex, - const unsigned int i) const -{ - Assert (dof_handler != 0, DoFAccessor::ExcInvalidObject()); - Assert (mg_dof_handler != 0, DoFAccessor::ExcInvalidObject()); - Assert (&dof_handler->get_fe() != 0, DoFAccessor::ExcInvalidObject()); - Assert (vertex<2, ExcIndexRange (i,0,2)); - Assert (iget_fe().dofs_per_vertex, - ExcIndexRange (i, 0, dof_handler->get_fe().dofs_per_vertex)); - - return (mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)] - .get_index (present_level, i, dof_handler->get_fe().dofs_per_vertex)); -}; - - - -template -void MGDoFObjectAccessor<1, dim>::set_mg_vertex_dof_index (const unsigned int vertex, - const unsigned int i, - const int index) const -{ - Assert (dof_handler != 0, DoFAccessor::ExcInvalidObject()); - Assert (mg_dof_handler != 0, DoFAccessor::ExcInvalidObject()); - Assert (&dof_handler->get_fe() != 0, DoFAccessor::ExcInvalidObject()); - Assert (vertex<2, ExcIndexRange (i,0,2)); - Assert (iget_fe().dofs_per_vertex, - ExcIndexRange (i, 0, dof_handler->get_fe().dofs_per_vertex)); - - mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)] - .set_index (present_level, i, dof_handler->get_fe().dofs_per_vertex, index); -}; - - - -template -void -MGDoFObjectAccessor<1, dim>::get_mg_dof_indices (vector &dof_indices) const -{ - Assert (dof_handler != 0, DoFAccessor::ExcInvalidObject()); - Assert (mg_dof_handler != 0, DoFAccessor::ExcInvalidObject()); - Assert (&dof_handler->get_fe() != 0, DoFAccessor::ExcInvalidObject()); - Assert (dof_indices.size() == (2*dof_handler->get_fe().dofs_per_vertex + - dof_handler->get_fe().dofs_per_line), - DoFAccessor::ExcVectorDoesNotMatch()); - - const unsigned int dofs_per_vertex = dof_handler->get_fe().dofs_per_vertex, - dofs_per_line = dof_handler->get_fe().dofs_per_line; - vector::iterator next = dof_indices.begin(); - for (unsigned int vertex=0; vertex<2; ++vertex) - for (unsigned int d=0; d -template -void -MGDoFObjectAccessor<1,dim>::get_mg_dof_values (const Vector &values, - Vector &dof_values) const -{ - Assert (dof_handler != 0, DoFAccessor<1>::DoFAccessor<1>::ExcInvalidObject()); - Assert (mg_dof_handler != 0, DoFAccessor<1>::DoFAccessor<1>::ExcInvalidObject()); - Assert (&dof_handler->get_fe() != 0, DoFAccessor<1>::DoFAccessor<1>::ExcInvalidObject()); - Assert (dof_values.size() == dof_handler->get_fe().dofs_per_cell, - DoFAccessor<1>::ExcVectorDoesNotMatch()); - Assert (values.size() == dof_handler->n_dofs(), - DoFAccessor<1>::ExcVectorDoesNotMatch()); - - const unsigned int dofs_per_vertex = dof_handler->get_fe().dofs_per_vertex, - dofs_per_line = dof_handler->get_fe().dofs_per_line; - vector::iterator next_dof_value=dof_values.begin(); - for (unsigned int vertex=0; vertex<2; ++vertex) - for (unsigned int d=0; d -TriaIterator > -MGDoFObjectAccessor<1, dim>::child (const unsigned int i) const -{ - TriaIterator > q (tria, - present_level+1, - child_index (i), - mg_dof_handler); - -#ifdef DEBUG - if (q.state() != past_the_end) - Assert (q->used(), typename TriaAccessor::ExcUnusedCellAsChild()); -#endif - return q; -}; - - - -template -void -MGDoFObjectAccessor<1, dim>::copy_from (const MGDoFObjectAccessor<1, dim> &a) -{ - DoFObjectAccessor<1, dim>::copy_from (a); - set_mg_dof_handler (a.mg_dof_handler); -}; - - - - -/* ------------------------ MGDoFQuadAccessor --------------------------- */ - -template -MGDoFObjectAccessor<2, dim>::MGDoFObjectAccessor (Triangulation *tria, - const int level, - const int index, - const AccessorData *local_data) : - MGDoFAccessor (local_data), - MGDoFObjectAccessor_Inheritance<2,dim>::BaseClass(tria,level,index,local_data) -{}; - - - -template -inline -int MGDoFObjectAccessor<2, dim>::mg_dof_index (const unsigned int i) const -{ - Assert (dof_handler != 0, DoFAccessor::ExcInvalidObject()); - Assert (mg_dof_handler != 0, DoFAccessor::ExcInvalidObject()); - // make sure a FE has been selected - // and enough room was reserved - Assert (&dof_handler->get_fe() != 0, DoFAccessor::ExcInvalidObject()); - Assert (iget_fe().dofs_per_quad, - ExcIndexRange (i, 0, dof_handler->get_fe().dofs_per_quad)); - - return mg_dof_handler->mg_levels[present_level] - ->quad_dofs[present_index*dof_handler->get_fe().dofs_per_quad+i]; -}; - - - -template -void MGDoFObjectAccessor<2, dim>::set_mg_dof_index (const unsigned int i, - const int index) const -{ - Assert (dof_handler != 0, DoFAccessor::ExcInvalidObject()); - Assert (mg_dof_handler != 0, DoFAccessor::ExcInvalidObject()); - // make sure a FE has been selected - // and enough room was reserved - Assert (&dof_handler->get_fe() != 0, DoFAccessor::ExcInvalidObject()); - Assert (iget_fe().dofs_per_quad, - ExcIndexRange (i, 0, dof_handler->get_fe().dofs_per_quad)); - - mg_dof_handler->mg_levels[present_level] - ->quad_dofs[present_index*dof_handler->get_fe().dofs_per_quad+i] = index; -}; - - - -template -inline -int MGDoFObjectAccessor<2, dim>::mg_vertex_dof_index (const unsigned int vertex, - const unsigned int i) const -{ - Assert (dof_handler != 0, DoFAccessor::ExcInvalidObject()); - Assert (mg_dof_handler != 0, DoFAccessor::ExcInvalidObject()); - Assert (&dof_handler->get_fe() != 0, DoFAccessor::ExcInvalidObject()); - Assert (vertex<4, ExcIndexRange (i,0,4)); - Assert (iget_fe().dofs_per_vertex, - ExcIndexRange (i, 0, dof_handler->get_fe().dofs_per_vertex)); - - return (mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)] - .get_index (present_level, i, dof_handler->get_fe().dofs_per_vertex)); -}; - - - -template -void MGDoFObjectAccessor<2, dim>::set_mg_vertex_dof_index (const unsigned int vertex, - const unsigned int i, - const int index) const -{ - Assert (dof_handler != 0, DoFAccessor::ExcInvalidObject()); - Assert (mg_dof_handler != 0, DoFAccessor::ExcInvalidObject()); - Assert (&dof_handler->get_fe() != 0, DoFAccessor::ExcInvalidObject()); - Assert (vertex<4, ExcIndexRange (i,0,4)); - Assert (iget_fe().dofs_per_vertex, - ExcIndexRange (i, 0, dof_handler->get_fe().dofs_per_vertex)); - - mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)] - .set_index (present_level, i, dof_handler->get_fe().dofs_per_vertex, index); -}; - - - -template -void -MGDoFObjectAccessor<2, dim>::get_mg_dof_indices (vector &dof_indices) const -{ - Assert (dof_handler != 0, DoFAccessor::ExcInvalidObject()); - Assert (mg_dof_handler != 0, DoFAccessor::ExcInvalidObject()); - Assert (&dof_handler->get_fe() != 0, DoFAccessor::ExcInvalidObject()); - Assert (dof_indices.size() == (4*dof_handler->get_fe().dofs_per_vertex + - 4*dof_handler->get_fe().dofs_per_line + - dof_handler->get_fe().dofs_per_quad), - DoFAccessor<2>::ExcVectorDoesNotMatch()); - - const unsigned int dofs_per_vertex = dof_handler->get_fe().dofs_per_vertex, - dofs_per_line = dof_handler->get_fe().dofs_per_line, - dofs_per_quad = dof_handler->get_fe().dofs_per_quad; - vector::iterator next = dof_indices.begin(); - for (unsigned int vertex=0; vertex<4; ++vertex) - for (unsigned int d=0; dline(line)->mg_dof_index(d); - for (unsigned int d=0; d -template -void -MGDoFObjectAccessor<2,dim>::get_mg_dof_values (const Vector &values, - Vector &dof_values) const -{ - Assert (dof_handler != 0, DoFAccessor<2>::ExcInvalidObject()); - Assert (mg_dof_handler != 0, DoFAccessor<2>::ExcInvalidObject()); - Assert (&dof_handler->get_fe() != 0, DoFAccessor<2>::ExcInvalidObject()); - Assert (dof_values.size() == dof_handler->get_fe().dofs_per_cell, - DoFAccessor<2>::ExcVectorDoesNotMatch()); - Assert (values.size() == mg_dof_handler->n_dofs(present_level), - DoFAccessor<2>::ExcVectorDoesNotMatch()); - - const unsigned int dofs_per_vertex = dof_handler->get_fe().dofs_per_vertex, - dofs_per_line = dof_handler->get_fe().dofs_per_line, - dofs_per_quad = dof_handler->get_fe().dofs_per_quad; - vector::iterator next_dof_value=dof_values.begin(); - for (unsigned int vertex=0; vertex<4; ++vertex) - for (unsigned int d=0; dline(line)->mg_dof_index(d)); - for (unsigned int d=0; d -TriaIterator > -MGDoFObjectAccessor<2, dim>::line (const unsigned int i) const -{ - Assert (i<4, ExcIndexRange (i, 0, 4)); - - return TriaIterator > - ( - tria, - present_level, - line_index (i), - mg_dof_handler - ); -}; - - - -template -TriaIterator > -MGDoFObjectAccessor<2, dim>::child (const unsigned int i) const -{ - TriaIterator > q (tria, - present_level+1, - child_index (i), - mg_dof_handler); - -#ifdef DEBUG - if (q.state() != past_the_end) - Assert (q->used(), typename TriaAccessor::ExcUnusedCellAsChild()); -#endif - return q; -}; - - - -template -void -MGDoFObjectAccessor<2, dim>::copy_from (const MGDoFObjectAccessor<2, dim> &a) -{ - DoFObjectAccessor<2, dim>::copy_from (a); - set_mg_dof_handler (a.mg_dof_handler); -}; - - - - -/* ------------------------ MGDoFHexAccessor --------------------------- */ - -template -MGDoFObjectAccessor<3, dim>::MGDoFObjectAccessor (Triangulation *tria, - const int level, - const int index, - const AccessorData *local_data) : - MGDoFAccessor (local_data), - MGDoFObjectAccessor_Inheritance<3,dim>::BaseClass(tria,level,index,local_data) -{}; - - - -template -inline -int MGDoFObjectAccessor<3, dim>::mg_dof_index (const unsigned int i) const -{ - Assert (dof_handler != 0, DoFAccessor::ExcInvalidObject()); - Assert (mg_dof_handler != 0, DoFAccessor::ExcInvalidObject()); - // make sure a FE has been selected - // and enough room was reserved - Assert (&dof_handler->get_fe() != 0, DoFAccessor::ExcInvalidObject()); - Assert (iget_fe().dofs_per_hex, - ExcIndexRange (i, 0, dof_handler->get_fe().dofs_per_hex)); - - return mg_dof_handler->mg_levels[present_level] - ->hex_dofs[present_index*dof_handler->get_fe().dofs_per_hex+i]; -}; - - - -template -void MGDoFObjectAccessor<3, dim>::set_mg_dof_index (const unsigned int i, - const int index) const -{ - Assert (dof_handler != 0, DoFAccessor::ExcInvalidObject()); - Assert (mg_dof_handler != 0, DoFAccessor::ExcInvalidObject()); - // make sure a FE has been selected - // and enough room was reserved - Assert (&dof_handler->get_fe() != 0, DoFAccessor::ExcInvalidObject()); - Assert (iget_fe().dofs_per_hex, - ExcIndexRange (i, 0, dof_handler->get_fe().dofs_per_hex)); - - mg_dof_handler->mg_levels[present_level] - ->hex_dofs[present_index*dof_handler->get_fe().dofs_per_hex+i] = index; -}; - - - -template -inline -int MGDoFObjectAccessor<3, dim>::mg_vertex_dof_index (const unsigned int vertex, - const unsigned int i) const -{ - Assert (dof_handler != 0, DoFAccessor::ExcInvalidObject()); - Assert (mg_dof_handler != 0, DoFAccessor::ExcInvalidObject()); - Assert (&dof_handler->get_fe() != 0, DoFAccessor::ExcInvalidObject()); - Assert (vertex<8, ExcIndexRange (i,0,8)); - Assert (iget_fe().dofs_per_vertex, - ExcIndexRange (i, 0, dof_handler->get_fe().dofs_per_vertex)); - - return (mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)] - .get_index (present_level, i, dof_handler->get_fe().dofs_per_vertex)); -}; - - - -template -void MGDoFObjectAccessor<3, dim>::set_mg_vertex_dof_index (const unsigned int vertex, - const unsigned int i, - const int index) const -{ - Assert (dof_handler != 0, DoFAccessor::ExcInvalidObject()); - Assert (mg_dof_handler != 0, DoFAccessor::ExcInvalidObject()); - Assert (&dof_handler->get_fe() != 0, DoFAccessor::ExcInvalidObject()); - Assert (vertex<4, ExcIndexRange (i,0,4)); - Assert (iget_fe().dofs_per_vertex, - ExcIndexRange (i, 0, dof_handler->get_fe().dofs_per_vertex)); - - mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)] - .set_index (present_level, i, dof_handler->get_fe().dofs_per_vertex, index); -}; - - - -template -void -MGDoFObjectAccessor<3, dim>::get_mg_dof_indices (vector &dof_indices) const -{ - Assert (dof_handler != 0, DoFAccessor::ExcInvalidObject()); - Assert (mg_dof_handler != 0, DoFAccessor::ExcInvalidObject()); - Assert (&dof_handler->get_fe() != 0, DoFAccessor::ExcInvalidObject()); - Assert (dof_indices.size() == (8*dof_handler->get_fe().dofs_per_vertex + - 12*dof_handler->get_fe().dofs_per_line + - 6*dof_handler->get_fe().dofs_per_quad + - dof_handler->get_fe().dofs_per_hex), - DoFAccessor<3>::ExcVectorDoesNotMatch()); - - const unsigned int dofs_per_vertex = dof_handler->get_fe().dofs_per_vertex, - dofs_per_line = dof_handler->get_fe().dofs_per_line, - dofs_per_quad = dof_handler->get_fe().dofs_per_quad, - dofs_per_hex = dof_handler->get_fe().dofs_per_hex; - vector::iterator next = dof_indices.begin(); - for (unsigned int vertex=0; vertex<8; ++vertex) - for (unsigned int d=0; dline(line)->mg_dof_index(d); - for (unsigned int quad=0; quad<12; ++quad) - for (unsigned int d=0; dquad(quad)->mg_dof_index(d); - for (unsigned int d=0; d -template -void -MGDoFObjectAccessor<3,dim>::get_mg_dof_values (const Vector &values, - Vector &dof_values) const -{ - Assert (dof_handler != 0, DoFAccessor<2>::ExcInvalidObject()); - Assert (mg_dof_handler != 0, DoFAccessor<2>::ExcInvalidObject()); - Assert (&dof_handler->get_fe() != 0, DoFAccessor<2>::ExcInvalidObject()); - Assert (dof_values.size() == dof_handler->get_fe().dofs_per_cell, - DoFAccessor<3>::ExcVectorDoesNotMatch()); - Assert (values.size() == mg_dof_handler->n_dofs(present_level), - DoFAccessor<3>::ExcVectorDoesNotMatch()); - - const unsigned int dofs_per_vertex = dof_handler->get_fe().dofs_per_vertex, - dofs_per_line = dof_handler->get_fe().dofs_per_line, - dofs_per_quad = dof_handler->get_fe().dofs_per_quad, - dofs_per_hex = dof_handler->get_fe().dofs_per_hex; - vector::iterator next_dof_value=dof_values.begin(); - for (unsigned int vertex=0; vertex<8; ++vertex) - for (unsigned int d=0; dline(line)->mg_dof_index(d)); - for (unsigned int quad=0; quad<6; ++quad) - for (unsigned int d=0; dquad(quad)->mg_dof_index(d)); - for (unsigned int d=0; d -TriaIterator > -MGDoFObjectAccessor<3, dim>::line (const unsigned int i) const { - Assert (i<12, ExcIndexRange (i, 0, 12)); - - return TriaIterator > - ( - tria, - present_level, - line_index (i), - mg_dof_handler - ); -}; - - - -template -TriaIterator > -MGDoFObjectAccessor<3, dim>::quad (const unsigned int i) const { - Assert (i<12, ExcIndexRange (i, 0, 6)); - - return TriaIterator > - ( - tria, - present_level, - quad_index (i), - mg_dof_handler - ); -}; - - - -template -TriaIterator > -MGDoFObjectAccessor<3, dim>::child (const unsigned int i) const { - TriaIterator > q (tria, - present_level+1, - child_index (i), - mg_dof_handler); - -#ifdef DEBUG - if (q.state() != past_the_end) - Assert (q->used(), typename TriaAccessor::ExcUnusedCellAsChild()); -#endif - return q; -}; - - - -template -void -MGDoFObjectAccessor<3, dim>::copy_from (const MGDoFObjectAccessor<3, dim> &a) { - DoFObjectAccessor<3, dim>::copy_from (a); - set_mg_dof_handler (a.mg_dof_handler); -}; - - - - - -/*------------------------- Functions: MGDoFCellAccessor -----------------------*/ - -template -TriaIterator > -MGDoFCellAccessor::neighbor (const unsigned int i) const { - TriaIterator > q (tria, - neighbor_level (i), - neighbor_index (i), - mg_dof_handler); - -#ifdef DEBUG - if (q.state() != past_the_end) - Assert (q->used(), typename TriaAccessor::ExcUnusedCellAsNeighbor()); -#endif - return q; -}; - - - -template -TriaIterator > -MGDoFCellAccessor::child (const unsigned int i) const { - TriaIterator > q (tria, - present_level+1, - child_index (i), - mg_dof_handler); - -#ifdef DEBUG - if (q.state() != past_the_end) - Assert (q->used(), typename TriaAccessor::ExcUnusedCellAsChild()); -#endif - return q; -}; - - - -#if deal_II_dimension == 1 - -template <> -MGDoFCellAccessor<1>::face_iterator -MGDoFCellAccessor<1>::face (const unsigned int) const -{ - Assert (false, ExcNotUsefulForThisDimension()); - return face_iterator(); -}; - -#endif - - - -#if deal_II_dimension == 2 - -template <> -MGDoFCellAccessor<2>::face_iterator -MGDoFCellAccessor<2>::face (const unsigned int i) const -{ - return line(i); -}; - -#endif - - - -#if deal_II_dimension == 3 - -template <> -MGDoFCellAccessor<3>::face_iterator -MGDoFCellAccessor<3>::face (const unsigned int i) const -{ - return quad(i); -}; - -#endif - - - - - -// explicit instantiations - -template -void -MGDoFObjectAccessor<1,deal_II_dimension>:: -get_mg_dof_values (const Vector &values, - Vector &dof_values) const; - -template -void -MGDoFObjectAccessor<1,deal_II_dimension>:: -get_mg_dof_values (const Vector &values, - Vector &dof_values) const; - - - -#if deal_II_dimension >= 2 - -template -void -MGDoFObjectAccessor<2,deal_II_dimension>:: -get_mg_dof_values (const Vector &values, - Vector &dof_values) const; - -template -void -MGDoFObjectAccessor<2,deal_II_dimension>:: -get_mg_dof_values (const Vector &values, - Vector &dof_values) const; - -#endif - - -#if deal_II_dimension >= 3 - -template -void -MGDoFObjectAccessor<3,deal_II_dimension>:: -get_mg_dof_values (const Vector &values, - Vector &dof_values) const; - -template -void -MGDoFObjectAccessor<3,deal_II_dimension>:: -get_mg_dof_values (const Vector &values, - Vector &dof_values) const; - -#endif - - -#if deal_II_dimension == 1 -template class MGDoFObjectAccessor<1, 1>; -template class MGDoFCellAccessor<1>; - -template class TriaRawIterator<1,MGDoFCellAccessor<1> >; -template class TriaIterator<1,MGDoFCellAccessor<1> >; -template class TriaActiveIterator<1,MGDoFCellAccessor<1> >; -#endif - -#if deal_II_dimension == 2 -template class MGDoFObjectAccessor<1, 2>; -template class MGDoFObjectAccessor<2, 2>; -template class MGDoFCellAccessor<2>; - -template class TriaRawIterator <2,MGDoFObjectAccessor<1, 2> >; -template class TriaIterator <2,MGDoFObjectAccessor<1, 2> >; -template class TriaActiveIterator<2,MGDoFObjectAccessor<1, 2> >; - -template class TriaRawIterator <2,MGDoFCellAccessor<2> >; -template class TriaIterator <2,MGDoFCellAccessor<2> >; -template class TriaActiveIterator<2,MGDoFCellAccessor<2> >; -#endif - - -#if deal_II_dimension == 3 -template class MGDoFObjectAccessor<1, 3>; -template class MGDoFObjectAccessor<2, 3>; -template class MGDoFObjectAccessor<3, 3>; -template class MGDoFCellAccessor<3>; - -template class TriaRawIterator <3,MGDoFObjectAccessor<1, 3> >; -template class TriaIterator <3,MGDoFObjectAccessor<1, 3> >; -template class TriaActiveIterator<3,MGDoFObjectAccessor<1, 3> >; - -template class TriaRawIterator <3,MGDoFObjectAccessor<2, 3> >; -template class TriaIterator <3,MGDoFObjectAccessor<2, 3> >; -template class TriaActiveIterator<3,MGDoFObjectAccessor<2, 3> >; - -template class TriaRawIterator <3,MGDoFCellAccessor<3> >; -template class TriaIterator <3,MGDoFCellAccessor<3> >; -template class TriaActiveIterator<3,MGDoFCellAccessor<3> >; -#endif - - diff --git a/deal.II/deal.II/source/dofs/mg_dof_handler.cc b/deal.II/deal.II/source/dofs/mg_dof_handler.cc deleted file mode 100644 index d880f442b1..0000000000 --- a/deal.II/deal.II/source/dofs/mg_dof_handler.cc +++ /dev/null @@ -1,1931 +0,0 @@ -/* $Id$ */ -/* Copyright W. Bangerth, University of Heidelberg, 1998 */ - - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include - - - -/* ------------------------ MGVertexDoFs ----------------------------------- */ - -template -MGDoFHandler::MGVertexDoFs::MGVertexDoFs () : - coarsest_level (1<<30), - finest_level (0), - indices (0) -{}; - - - -template -void MGDoFHandler::MGVertexDoFs::init (const unsigned int cl, - const unsigned int fl, - const unsigned int dofs_per_vertex) { - Assert (indices == 0, ExcInternalError()); - - coarsest_level = cl; - finest_level = fl; - - const unsigned int n_levels = finest_level-coarsest_level + 1; - - indices = new int[n_levels * dofs_per_vertex]; - Assert (indices != 0, ExcNoMemory ()); - - for (unsigned int i=0; i -MGDoFHandler::MGVertexDoFs::~MGVertexDoFs () { - delete[] indices; -}; - - - -template -inline -unsigned int MGDoFHandler::MGVertexDoFs::get_coarsest_level () const { - return coarsest_level; -}; - - - -template -inline -unsigned int MGDoFHandler::MGVertexDoFs::get_finest_level () const { - return finest_level; -}; - - - - - - -/* ------------------------ MGDoFHandler ------------------------------------- */ - -template -MGDoFHandler::MGDoFHandler (Triangulation &tria) : - DoFHandler (tria) -{}; - - - -template -MGDoFHandler::~MGDoFHandler () {}; - - - - - - -#if deal_II_dimension == 1 - -template <> -MGDoFHandler<1>::raw_cell_iterator -MGDoFHandler<1>::begin_raw (const unsigned int level) const { - return begin_raw_line (level); -}; - - - -template <> -MGDoFHandler<1>::cell_iterator -MGDoFHandler<1>::begin (const unsigned int level) const { - return begin_line (level); -}; - - - -template <> -MGDoFHandler<1>::active_cell_iterator -MGDoFHandler<1>::begin_active (const unsigned int level) const { - return begin_active_line (level); -}; - - - -template <> -MGDoFHandler<1>::raw_cell_iterator -MGDoFHandler<1>::end () const { - return end_line (); -}; - - - - -template <> -MGDoFHandler<1>::raw_cell_iterator -MGDoFHandler<1>::last_raw () const { - return last_raw_line (); -}; - - - -template <> -MGDoFHandler<1>::raw_cell_iterator -MGDoFHandler<1>::last_raw (const unsigned int level) const { - return last_raw_line (level); -}; - - - -template <> -MGDoFHandler<1>::cell_iterator -MGDoFHandler<1>::last () const { - return last_line (); -}; - - - -template <> -MGDoFHandler<1>::cell_iterator -MGDoFHandler<1>::last (const unsigned int level) const { - return last_line (level); -}; - - - -template <> -MGDoFHandler<1>::active_cell_iterator -MGDoFHandler<1>::last_active () const { - return last_active_line (); -}; - - - -template <> -MGDoFHandler<1>::active_cell_iterator -MGDoFHandler<1>::last_active (const unsigned int level) const { - return last_active_line (level); -}; - - - -template <> -MGDoFDimensionInfo<1>::raw_face_iterator -MGDoFHandler<1>::begin_raw_face (const unsigned int) const { - Assert (false, ExcFunctionNotUseful()); - return 0; -}; - - - -template <> -MGDoFDimensionInfo<1>::face_iterator -MGDoFHandler<1>::begin_face (const unsigned int) const { - Assert (false, ExcFunctionNotUseful()); - return 0; -}; - - - -template <> -MGDoFDimensionInfo<1>::active_face_iterator -MGDoFHandler<1>::begin_active_face (const unsigned int) const { - Assert (false, ExcFunctionNotUseful()); - return 0; -}; - - - -template <> -MGDoFDimensionInfo<1>::raw_face_iterator -MGDoFHandler<1>::end_face () const { - Assert (false, ExcFunctionNotUseful()); - return 0; -}; - - - -template <> -MGDoFDimensionInfo<1>::raw_face_iterator -MGDoFHandler<1>::last_raw_face () const { - Assert (false, ExcFunctionNotUseful()); - return 0; -}; - - - -template <> -MGDoFDimensionInfo<1>::raw_face_iterator -MGDoFHandler<1>::last_raw_face (const unsigned int) const { - Assert (false, ExcFunctionNotUseful()); - return 0; -}; - - - -template <> -MGDoFDimensionInfo<1>::face_iterator -MGDoFHandler<1>::last_face () const { - Assert (false, ExcFunctionNotUseful()); - return 0; -}; - - - -template <> -MGDoFDimensionInfo<1>::face_iterator -MGDoFHandler<1>::last_face (const unsigned int) const { - Assert (false, ExcFunctionNotUseful()); - return 0; -}; - - - -template <> -MGDoFDimensionInfo<1>::active_face_iterator -MGDoFHandler<1>::last_active_face () const { - Assert (false, ExcFunctionNotUseful()); - return 0; -}; - - - -template <> -MGDoFDimensionInfo<1>::active_face_iterator -MGDoFHandler<1>::last_active_face (const unsigned int) const { - Assert (false, ExcFunctionNotUseful()); - return 0; -}; - - - - -template <> -MGDoFHandler<1>::raw_quad_iterator -MGDoFHandler<1>::begin_raw_quad (const unsigned int) const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - - -template <> -MGDoFHandler<1>::quad_iterator -MGDoFHandler<1>::begin_quad (const unsigned int) const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - - -template <> -MGDoFHandler<1>::active_quad_iterator -MGDoFHandler<1>::begin_active_quad (const unsigned int) const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - - -template <> -MGDoFHandler<1>::raw_quad_iterator -MGDoFHandler<1>::end_quad () const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - - -template <> -MGDoFHandler<1>::raw_quad_iterator -MGDoFHandler<1>::last_raw_quad (const unsigned int) const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - -template <> -MGDoFHandler<1>::quad_iterator -MGDoFHandler<1>::last_quad (const unsigned int) const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - - -template <> -MGDoFHandler<1>::active_quad_iterator -MGDoFHandler<1>::last_active_quad (const unsigned int) const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - - -template <> -MGDoFHandler<1>::raw_quad_iterator -MGDoFHandler<1>::last_raw_quad () const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - - -template <> -MGDoFHandler<1>::quad_iterator -MGDoFHandler<1>::last_quad () const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - - -template <> -MGDoFHandler<1>::active_quad_iterator -MGDoFHandler<1>::last_active_quad () const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - -template <> -MGDoFHandler<1>::raw_hex_iterator -MGDoFHandler<1>::begin_raw_hex (const unsigned int) const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - - -template <> -MGDoFHandler<1>::hex_iterator -MGDoFHandler<1>::begin_hex (const unsigned int) const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - - -template <> -MGDoFHandler<1>::active_hex_iterator -MGDoFHandler<1>::begin_active_hex (const unsigned int) const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - - -template <> -MGDoFHandler<1>::raw_hex_iterator -MGDoFHandler<1>::end_hex () const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - - -template <> -MGDoFHandler<1>::raw_hex_iterator -MGDoFHandler<1>::last_raw_hex (const unsigned int) const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - -template <> -MGDoFHandler<1>::hex_iterator -MGDoFHandler<1>::last_hex (const unsigned int) const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - - -template <> -MGDoFHandler<1>::active_hex_iterator -MGDoFHandler<1>::last_active_hex (const unsigned int) const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - - -template <> -MGDoFHandler<1>::raw_hex_iterator -MGDoFHandler<1>::last_raw_hex () const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - - -template <> -MGDoFHandler<1>::hex_iterator -MGDoFHandler<1>::last_hex () const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - - -template <> -MGDoFHandler<1>::active_hex_iterator -MGDoFHandler<1>::last_active_hex () const { - Assert (false, ExcNotImplemented()); - return 0; -}; - -#endif - - - -#if deal_II_dimension == 2 - -template <> -MGDoFHandler<2>::raw_cell_iterator -MGDoFHandler<2>::begin_raw (const unsigned int level) const { - return begin_raw_quad (level); -}; - - - -template <> -MGDoFHandler<2>::cell_iterator -MGDoFHandler<2>::begin (const unsigned int level) const { - return begin_quad (level); -}; - - - -template <> -MGDoFHandler<2>::active_cell_iterator -MGDoFHandler<2>::begin_active (const unsigned int level) const { - return begin_active_quad (level); -}; - - - -template <> -MGDoFHandler<2>::raw_cell_iterator -MGDoFHandler<2>::end () const { - return end_quad (); -}; - - - -template <> -MGDoFHandler<2>::raw_cell_iterator -MGDoFHandler<2>::last_raw () const { - return last_raw_quad (); -}; - - - -template <> -MGDoFHandler<2>::raw_cell_iterator -MGDoFHandler<2>::last_raw (const unsigned int level) const { - return last_raw_quad (level); -}; - - - -template <> -MGDoFHandler<2>::cell_iterator -MGDoFHandler<2>::last () const { - return last_quad (); -}; - - - -template <> -MGDoFHandler<2>::cell_iterator -MGDoFHandler<2>::last (const unsigned int level) const { - return last_quad (level); -}; - - - -template <> -MGDoFHandler<2>::active_cell_iterator -MGDoFHandler<2>::last_active () const { - return last_active_quad (); -}; - - - -template <> -MGDoFHandler<2>::active_cell_iterator -MGDoFHandler<2>::last_active (const unsigned int level) const { - return last_active_quad (level); -}; - - -template <> -MGDoFDimensionInfo<2>::raw_face_iterator -MGDoFHandler<2>::begin_raw_face (const unsigned int level) const { - return begin_raw_line (level); -}; - - - -template <> -MGDoFDimensionInfo<2>::face_iterator -MGDoFHandler<2>::begin_face (const unsigned int level) const { - return begin_line (level); -}; - - - -template <> -MGDoFDimensionInfo<2>::active_face_iterator -MGDoFHandler<2>::begin_active_face (const unsigned int level) const { - return begin_active_line (level); -}; - - - -template <> -MGDoFDimensionInfo<2>::raw_face_iterator -MGDoFHandler<2>::end_face () const { - return end_line (); -}; - - - -template <> -MGDoFDimensionInfo<2>::raw_face_iterator -MGDoFHandler<2>::last_raw_face () const { - return last_raw_line (); -}; - - - -template <> -MGDoFDimensionInfo<2>::raw_face_iterator -MGDoFHandler<2>::last_raw_face (const unsigned int level) const { - return last_raw_line (level); -}; - - - -template <> -MGDoFDimensionInfo<2>::face_iterator -MGDoFHandler<2>::last_face () const { - return last_line (); -}; - - - -template <> -MGDoFDimensionInfo<2>::face_iterator -MGDoFHandler<2>::last_face (const unsigned int level) const { - return last_line (level); -}; - - - -template <> -MGDoFDimensionInfo<2>::active_face_iterator -MGDoFHandler<2>::last_active_face () const { - return last_active_line (); -}; - - - -template <> -MGDoFDimensionInfo<2>::active_face_iterator -MGDoFHandler<2>::last_active_face (const unsigned int level) const { - return last_active_line (level); -}; - - - -template <> -MGDoFHandler<2>::raw_hex_iterator -MGDoFHandler<2>::begin_raw_hex (const unsigned int) const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - - -template <> -MGDoFHandler<2>::hex_iterator -MGDoFHandler<2>::begin_hex (const unsigned int) const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - - -template <> -MGDoFHandler<2>::active_hex_iterator -MGDoFHandler<2>::begin_active_hex (const unsigned int) const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - - -template <> -MGDoFHandler<2>::raw_hex_iterator -MGDoFHandler<2>::end_hex () const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - - -template <> -MGDoFHandler<2>::raw_hex_iterator -MGDoFHandler<2>::last_raw_hex (const unsigned int) const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - -template <> -MGDoFHandler<2>::hex_iterator -MGDoFHandler<2>::last_hex (const unsigned int) const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - - -template <> -MGDoFHandler<2>::active_hex_iterator -MGDoFHandler<2>::last_active_hex (const unsigned int) const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - - -template <> -MGDoFHandler<2>::raw_hex_iterator -MGDoFHandler<2>::last_raw_hex () const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - - -template <> -MGDoFHandler<2>::hex_iterator -MGDoFHandler<2>::last_hex () const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - - -template <> -MGDoFHandler<2>::active_hex_iterator -MGDoFHandler<2>::last_active_hex () const { - Assert (false, ExcNotImplemented()); - return 0; -}; - - - -#endif - - - - -#if deal_II_dimension == 3 - -template <> -MGDoFHandler<3>::raw_cell_iterator -MGDoFHandler<3>::begin_raw (const unsigned int level) const { - return begin_raw_hex (level); -}; - - - -template <> -MGDoFHandler<3>::cell_iterator -MGDoFHandler<3>::begin (const unsigned int level) const { - return begin_hex (level); -}; - - - -template <> -MGDoFHandler<3>::active_cell_iterator -MGDoFHandler<3>::begin_active (const unsigned int level) const { - return begin_active_hex (level); -}; - - - -template <> -MGDoFHandler<3>::raw_cell_iterator -MGDoFHandler<3>::end () const { - return end_hex (); -}; - - - -template <> -MGDoFHandler<3>::raw_cell_iterator -MGDoFHandler<3>::last_raw () const { - return last_raw_hex (); -}; - - - -template <> -MGDoFHandler<3>::raw_cell_iterator -MGDoFHandler<3>::last_raw (const unsigned int level) const { - return last_raw_hex (level); -}; - - - -template <> -MGDoFHandler<3>::cell_iterator -MGDoFHandler<3>::last () const { - return last_hex (); -}; - - - -template <> -MGDoFHandler<3>::cell_iterator -MGDoFHandler<3>::last (const unsigned int level) const { - return last_hex (level); -}; - - - -template <> -MGDoFHandler<3>::active_cell_iterator -MGDoFHandler<3>::last_active () const { - return last_active_hex (); -}; - - - -template <> -MGDoFHandler<3>::active_cell_iterator -MGDoFHandler<3>::last_active (const unsigned int level) const { - return last_active_hex (level); -}; - - -template <> -MGDoFDimensionInfo<3>::raw_face_iterator -MGDoFHandler<3>::begin_raw_face (const unsigned int level) const { - return begin_raw_quad (level); -}; - - - -template <> -MGDoFDimensionInfo<3>::face_iterator -MGDoFHandler<3>::begin_face (const unsigned int level) const { - return begin_quad (level); -}; - - - -template <> -MGDoFDimensionInfo<3>::active_face_iterator -MGDoFHandler<3>::begin_active_face (const unsigned int level) const { - return begin_active_quad (level); -}; - - - -template <> -MGDoFDimensionInfo<3>::raw_face_iterator -MGDoFHandler<3>::end_face () const { - return end_quad (); -}; - - - -template <> -MGDoFDimensionInfo<3>::raw_face_iterator -MGDoFHandler<3>::last_raw_face () const { - return last_raw_quad (); -}; - - - -template <> -MGDoFDimensionInfo<3>::raw_face_iterator -MGDoFHandler<3>::last_raw_face (const unsigned int level) const { - return last_raw_quad (level); -}; - - - -template <> -MGDoFDimensionInfo<3>::face_iterator -MGDoFHandler<3>::last_face () const { - return last_quad (); -}; - - - -template <> -MGDoFDimensionInfo<3>::face_iterator -MGDoFHandler<3>::last_face (const unsigned int level) const { - return last_quad (level); -}; - - - -template <> -MGDoFDimensionInfo<3>::active_face_iterator -MGDoFHandler<3>::last_active_face () const { - return last_active_quad (); -}; - - - -template <> -MGDoFDimensionInfo<3>::active_face_iterator -MGDoFHandler<3>::last_active_face (const unsigned int level) const { - return last_active_quad (level); -}; - -#endif - - - - - -template -typename MGDoFHandler::raw_line_iterator -MGDoFHandler::begin_raw_line (const unsigned int level) const { - return raw_line_iterator (tria, - tria->begin_raw_line(level)->level(), - tria->begin_raw_line(level)->index(), - this); -}; - - - -template -typename MGDoFHandler::line_iterator -MGDoFHandler::begin_line (const unsigned int level) const { - return line_iterator (tria, - tria->begin_line(level)->level(), - tria->begin_line(level)->index(), - this); -}; - - - -template -typename MGDoFHandler::active_line_iterator -MGDoFHandler::begin_active_line (const unsigned int level) const { - return active_line_iterator (tria, - tria->begin_active_line(level)->level(), - tria->begin_active_line(level)->index(), - this); -}; - - - -template -typename MGDoFHandler::raw_quad_iterator -MGDoFHandler::begin_raw_quad (const unsigned int level) const { - return raw_quad_iterator (tria, - tria->begin_raw_quad(level)->level(), - tria->begin_raw_quad(level)->index(), - this); -}; - - - -template -typename MGDoFHandler::quad_iterator -MGDoFHandler::begin_quad (const unsigned int level) const { - return quad_iterator (tria, - tria->begin_quad(level)->level(), - tria->begin_quad(level)->index(), - this); -}; - - - -template -typename MGDoFHandler::active_quad_iterator -MGDoFHandler::begin_active_quad (const unsigned int level) const { - return active_quad_iterator (tria, - tria->begin_active_quad(level)->level(), - tria->begin_active_quad(level)->index(), - this); -}; - - - -template -typename MGDoFHandler::raw_hex_iterator -MGDoFHandler::begin_raw_hex (const unsigned int level) const { - return raw_hex_iterator (tria, - tria->begin_raw_hex(level)->level(), - tria->begin_raw_hex(level)->index(), - this); -}; - - - -template -typename MGDoFHandler::hex_iterator -MGDoFHandler::begin_hex (const unsigned int level) const { - return hex_iterator (tria, - tria->begin_hex(level)->level(), - tria->begin_hex(level)->index(), - this); -}; - - - -template -typename MGDoFHandler::active_hex_iterator -MGDoFHandler::begin_active_hex (const unsigned int level) const { - return active_hex_iterator (tria, - tria->begin_active_hex(level)->level(), - tria->begin_active_hex(level)->index(), - this); -}; - - - -template -typename MGDoFHandler::raw_line_iterator -MGDoFHandler::end_line () const { - return raw_line_iterator (tria, -1, -1, this); -}; - - - -template -typename MGDoFHandler::raw_quad_iterator -MGDoFHandler::end_quad () const { - return raw_quad_iterator (tria, -1, -1, this); -}; - - - -template -typename MGDoFHandler::raw_hex_iterator -MGDoFHandler::end_hex () const { - return raw_hex_iterator (tria, -1, -1, this); -}; - - - -template -typename MGDoFDimensionInfo::raw_cell_iterator -MGDoFHandler::end_raw (const unsigned int level) const { - return (level == mg_levels.size()-1 ? - end() : - begin_raw (level+1)); -}; - - - -template -typename MGDoFDimensionInfo::cell_iterator -MGDoFHandler::end (const unsigned int level) const { - return (level == mg_levels.size()-1 ? - cell_iterator(end()) : - begin (level+1)); -}; - - - -template -typename MGDoFDimensionInfo::active_cell_iterator -MGDoFHandler::end_active (const unsigned int level) const { - return (level == mg_levels.size()-1 ? - active_cell_iterator(end()) : - begin_active (level+1)); -}; - - - -template -typename MGDoFDimensionInfo::raw_face_iterator -MGDoFHandler::end_raw_face (const unsigned int level) const { - return (level == mg_levels.size()-1 ? - end_face() : - begin_raw_face (level+1)); -}; - - - -template -typename MGDoFDimensionInfo::face_iterator -MGDoFHandler::end_face (const unsigned int level) const { - return (level == mg_levels.size()-1 ? - face_iterator(end_face()) : - begin_face (level+1)); -}; - - - -template -typename MGDoFDimensionInfo::active_face_iterator -MGDoFHandler::end_active_face (const unsigned int level) const { - return (level == mg_levels.size()-1 ? - active_face_iterator(end_face()) : - begin_active_face (level+1)); -}; - - - -template -typename MGDoFDimensionInfo::raw_line_iterator -MGDoFHandler::end_raw_line (const unsigned int level) const { - return (level == mg_levels.size()-1 ? - end_line() : - begin_raw_line (level+1)); -}; - - - -template -typename MGDoFDimensionInfo::line_iterator -MGDoFHandler::end_line (const unsigned int level) const { - return (level == mg_levels.size()-1 ? - line_iterator(end_line()) : - begin_line (level+1)); -}; - - - -template -typename MGDoFDimensionInfo::active_line_iterator -MGDoFHandler::end_active_line (const unsigned int level) const { - return (level == mg_levels.size()-1 ? - active_line_iterator(end_line()) : - begin_active_line (level+1)); -}; - - - - -template -typename MGDoFDimensionInfo::raw_quad_iterator -MGDoFHandler::end_raw_quad (const unsigned int level) const { - return (level == mg_levels.size()-1 ? - end_quad() : - begin_raw_quad (level+1)); -}; - - - -template -typename MGDoFDimensionInfo::quad_iterator -MGDoFHandler::end_quad (const unsigned int level) const { - return (level == mg_levels.size()-1 ? - quad_iterator(end_quad()) : - begin_quad (level+1)); -}; - - - - -template -typename MGDoFDimensionInfo::active_quad_iterator -MGDoFHandler::end_active_quad (const unsigned int level) const { - return (level == mg_levels.size()-1 ? - active_quad_iterator(end_quad()) : - begin_active_quad (level+1)); -}; - - - -template -typename MGDoFDimensionInfo::raw_hex_iterator -MGDoFHandler::end_raw_hex (const unsigned int level) const { - return (level == mg_levels.size()-1 ? - end_hex() : - begin_raw_hex (level+1)); -}; - - - -template -typename MGDoFDimensionInfo::hex_iterator -MGDoFHandler::end_hex (const unsigned int level) const { - return (level == mg_levels.size()-1 ? - hex_iterator(end_hex()) : - begin_hex (level+1)); -}; - - - - -template -typename MGDoFDimensionInfo::active_hex_iterator -MGDoFHandler::end_active_hex (const unsigned int level) const { - return (level == mg_levels.size()-1 ? - active_hex_iterator(end_hex()) : - begin_active_hex (level+1)); -}; - - - -template -typename MGDoFHandler::raw_line_iterator -MGDoFHandler::last_raw_line (const unsigned int level) const { - return raw_line_iterator (tria, - tria->last_raw_line(level)->level(), - tria->last_raw_line(level)->index(), - this); -}; - - - -template -typename MGDoFHandler::line_iterator -MGDoFHandler::last_line (const unsigned int level) const { - return line_iterator (tria, - tria->last_line(level)->level(), - tria->last_line(level)->index(), - this); -}; - - - -template -typename MGDoFHandler::active_line_iterator -MGDoFHandler::last_active_line (const unsigned int level) const { - return active_line_iterator (tria, - tria->last_active_line(level)->level(), - tria->last_active_line(level)->index(), - this); -}; - - - -template -typename MGDoFHandler::raw_quad_iterator -MGDoFHandler::last_raw_quad (const unsigned int level) const { - return raw_quad_iterator (tria, - tria->last_raw_quad(level)->level(), - tria->last_raw_quad(level)->index(), - this); -}; - - - - -template -typename MGDoFHandler::quad_iterator -MGDoFHandler::last_quad (const unsigned int level) const { - return quad_iterator (tria, - tria->last_quad(level)->level(), - tria->last_quad(level)->index(), - this); -}; - - - - -template -typename MGDoFHandler::active_quad_iterator -MGDoFHandler::last_active_quad (const unsigned int level) const { - return active_quad_iterator (tria, - tria->last_active_quad(level)->level(), - tria->last_active_quad(level)->index(), - this); -}; - - - - -template -typename MGDoFHandler::raw_hex_iterator -MGDoFHandler::last_raw_hex (const unsigned int level) const { - return raw_hex_iterator (tria, - tria->last_raw_hex(level)->level(), - tria->last_raw_hex(level)->index(), - this); -}; - - - - -template -typename MGDoFHandler::hex_iterator -MGDoFHandler::last_hex (const unsigned int level) const { - return hex_iterator (tria, - tria->last_hex(level)->level(), - tria->last_hex(level)->index(), - this); -}; - - - - -template -typename MGDoFHandler::active_hex_iterator -MGDoFHandler::last_active_hex (const unsigned int level) const { - return active_hex_iterator (tria, - tria->last_active_hex(level)->level(), - tria->last_active_hex(level)->index(), - this); -}; - - - - -template -typename MGDoFHandler::raw_line_iterator -MGDoFHandler::last_raw_line () const { - return last_raw_line (mg_levels.size()-1); -}; - - - -template -typename MGDoFHandler::raw_quad_iterator -MGDoFHandler::last_raw_quad () const { - return last_raw_quad (mg_levels.size()-1); -}; - - - -template -typename MGDoFHandler::raw_hex_iterator -MGDoFHandler::last_raw_hex () const { - return last_raw_hex (mg_levels.size()-1); -}; - - - -template -typename MGDoFHandler::line_iterator -MGDoFHandler::last_line () const { - return last_line (mg_levels.size()-1); -}; - - - -template -typename MGDoFHandler::quad_iterator -MGDoFHandler::last_quad () const { - return last_quad (mg_levels.size()-1); -}; - - - -template -typename MGDoFHandler::hex_iterator -MGDoFHandler::last_hex () const { - return last_hex (mg_levels.size()-1); -}; - - - -template -typename MGDoFHandler::active_line_iterator -MGDoFHandler::last_active_line () const { - return last_active_line (mg_levels.size()-1); -}; - - - - -template -typename MGDoFHandler::active_quad_iterator -MGDoFHandler::last_active_quad () const { - return last_active_quad (mg_levels.size()-1); -}; - - - -template -typename MGDoFHandler::active_hex_iterator -MGDoFHandler::last_active_hex () const { - return last_active_hex (mg_levels.size()-1); -}; - - - - - - -//------------------------------------------------------------------ - - - - - - - -template -void MGDoFHandler::distribute_dofs (const FiniteElement &fe, - unsigned int offset) -{ - // first distribute global dofs - DoFHandler::distribute_dofs (fe, offset); - - - // reserve space for the MG dof numbers - reserve_space (); - mg_used_dofs.resize (tria->n_levels(), 0); - - // clear user flags because we will - // need them - tria->clear_user_flags (); - - // now distribute indices on each level - // separately - for (unsigned int level=0; leveln_levels(); ++level) - { - unsigned int next_free_dof = offset; - cell_iterator cell = begin(level), - endc = end(level); - - for (; cell != endc; ++cell) - next_free_dof = distribute_dofs_on_cell (cell, next_free_dof); - - mg_used_dofs[level] = next_free_dof; - }; -}; - - - -#if deal_II_dimension == 1 - -template <> -unsigned int -MGDoFHandler<1>::distribute_dofs_on_cell (cell_iterator &cell, - unsigned int next_free_dof) { - - // distribute dofs of vertices - if (selected_fe->dofs_per_vertex > 0) - for (unsigned int v=0; v::vertices_per_cell; ++v) - { - cell_iterator neighbor = cell->neighbor(v); - - if (neighbor.state() == valid) - { - // has neighbor already been processed? - if (neighbor->user_flag_set() && - (neighbor->level() == cell->level())) - // copy dofs if the neighbor is on - // the same level (only then are - // mg dofs the same) - { - if (v==0) - for (unsigned int d=0; ddofs_per_vertex; ++d) - cell->set_mg_vertex_dof_index (0, d, - neighbor->mg_vertex_dof_index (1, d)); - else - for (unsigned int d=0; ddofs_per_vertex; ++d) - cell->set_mg_vertex_dof_index (1, d, - neighbor->mg_vertex_dof_index (0, d)); - - // next neighbor - continue; - }; - }; - - // otherwise: create dofs newly - for (unsigned int d=0; ddofs_per_vertex; ++d) - cell->set_mg_vertex_dof_index (v, d, next_free_dof++); - }; - - // dofs of line - if (selected_fe->dofs_per_line > 0) - for (unsigned int d=0; ddofs_per_line; ++d) - cell->set_mg_dof_index (d, next_free_dof++); - - // note that this cell has been processed - cell->set_user_flag (); - - return next_free_dof; -}; - -#endif - - -#if deal_II_dimension == 2 - -template <> -unsigned int -MGDoFHandler<2>::distribute_dofs_on_cell (cell_iterator &cell, - unsigned int next_free_dof) { - if (selected_fe->dofs_per_vertex > 0) - // number dofs on vertices - for (unsigned int vertex=0; vertex::vertices_per_cell; ++vertex) - // check whether dofs for this - // vertex have been distributed - // (only check the first dof) - if (cell->mg_vertex_dof_index(vertex, 0) == -1) - for (unsigned int d=0; ddofs_per_vertex; ++d) - cell->set_mg_vertex_dof_index (vertex, d, next_free_dof++); - - // for the four sides - if (selected_fe->dofs_per_line > 0) - for (unsigned int side=0; side::faces_per_cell; ++side) - { - line_iterator line = cell->line(side); - - // distribute dofs if necessary: - // check whether line dof is already - // numbered (check only first dof) - if (line->mg_dof_index(0) == -1) - // if not: distribute dofs - for (unsigned int d=0; ddofs_per_line; ++d) - line->set_mg_dof_index (d, next_free_dof++); - }; - - - // dofs of quad - if (selected_fe->dofs_per_quad > 0) - for (unsigned int d=0; ddofs_per_quad; ++d) - cell->set_mg_dof_index (d, next_free_dof++); - - - // note that this cell has been processed - cell->set_user_flag (); - - return next_free_dof; -}; - -#endif - - - -#if deal_II_dimension == 3 - -template <> -unsigned int -MGDoFHandler<3>::distribute_dofs_on_cell (cell_iterator &cell, - unsigned int next_free_dof) { - if (selected_fe->dofs_per_vertex > 0) - // number dofs on vertices - for (unsigned int vertex=0; vertex::vertices_per_cell; ++vertex) - // check whether dofs for this - // vertex have been distributed - // (only check the first dof) - if (cell->mg_vertex_dof_index(vertex, 0) == -1) - for (unsigned int d=0; ddofs_per_vertex; ++d) - cell->set_mg_vertex_dof_index (vertex, d, next_free_dof++); - - // for the lines - if (selected_fe->dofs_per_line > 0) - for (unsigned int l=0; l::lines_per_cell; ++l) - { - line_iterator line = cell->line(l); - - // distribute dofs if necessary: - // check whether line dof is already - // numbered (check only first dof) - if (line->mg_dof_index(0) == -1) - // if not: distribute dofs - for (unsigned int d=0; ddofs_per_line; ++d) - line->set_mg_dof_index (d, next_free_dof++); - }; - - // for the quads - if (selected_fe->dofs_per_quad > 0) - for (unsigned int q=0; q::lines_per_cell; ++q) - { - quad_iterator quad = cell->quad(q); - - // distribute dofs if necessary: - // check whether line dof is already - // numbered (check only first dof) - if (quad->mg_dof_index(0) == -1) - // if not: distribute dofs - for (unsigned int d=0; ddofs_per_line; ++d) - quad->set_mg_dof_index (d, next_free_dof++); - }; - - - // dofs of cell - if (selected_fe->dofs_per_hex > 0) - for (unsigned int d=0; ddofs_per_hex; ++d) - cell->set_mg_dof_index (d, next_free_dof++); - - - // note that this cell has been processed - cell->set_user_flag (); - - return next_free_dof; -}; - -#endif - - - - -template -unsigned int MGDoFHandler::n_dofs (const unsigned int level) const { - Assert (level < mg_used_dofs.size(), ExcInvalidLevel(level)); - - return mg_used_dofs[level]; -}; - - - -#if deal_II_dimension == 1 - -template <> -void MGDoFHandler<1>::renumber_dofs (const unsigned int level, - const vector &new_numbers) { - Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete()); - - // note that we can not use cell iterators - // in this function since then we would - // renumber the dofs on the interface of - // two cells more than once. Anyway, this - // ways it's not only more correct but also - // faster - for (vector::iterator i=mg_vertex_dofs.begin(); - i!=mg_vertex_dofs.end(); ++i) - // if the present vertex lives on - // the present level - if ((i->get_coarsest_level() <= level) && - (i->get_finest_level() >= level)) - for (unsigned int d=0; ddofs_per_vertex; ++d) - i->set_index (level, d, selected_fe->dofs_per_vertex, - new_numbers[i->get_index (level, d, - selected_fe->dofs_per_vertex)]); - - for (vector::iterator i=mg_levels[level]->line_dofs.begin(); - i!=mg_levels[level]->line_dofs.end(); ++i) - { - Assert (*i != -1, ExcInternalError()); - *i = new_numbers[*i]; - }; -}; - -#endif - - - - -#if deal_II_dimension == 2 - -template <> -void MGDoFHandler<2>::renumber_dofs (const unsigned int level, - const vector &new_numbers) { - Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete()); - - for (vector::iterator i=mg_vertex_dofs.begin(); - i!=mg_vertex_dofs.end(); ++i) - // if the present vertex lives on - // the present level - if ((i->get_coarsest_level() <= level) && - (i->get_finest_level() >= level)) - for (unsigned int d=0; ddofs_per_vertex; ++d) - i->set_index (level, d, selected_fe->dofs_per_vertex, - new_numbers[i->get_index (level, d, - selected_fe->dofs_per_vertex)]); - - for (vector::iterator i=mg_levels[level]->line_dofs.begin(); - i!=mg_levels[level]->line_dofs.end(); ++i) - { - Assert (*i != -1, ExcInternalError()); - *i = new_numbers[*i]; - }; - - for (vector::iterator i=mg_levels[level]->quad_dofs.begin(); - i!=mg_levels[level]->quad_dofs.end(); ++i) - { - Assert (*i != -1, ExcInternalError()); - *i = new_numbers[*i]; - }; -}; - -#endif - - - -#if deal_II_dimension == 3 - -template <> -void MGDoFHandler<3>::renumber_dofs (const unsigned int level, - const vector &new_numbers) { - Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete()); - - for (vector::iterator i=mg_vertex_dofs.begin(); - i!=mg_vertex_dofs.end(); ++i) - // if the present vertex lives on - // the present level - if ((i->get_coarsest_level() <= level) && - (i->get_finest_level() >= level)) - for (unsigned int d=0; ddofs_per_vertex; ++d) - i->set_index (level, d, selected_fe->dofs_per_vertex, - new_numbers[i->get_index (level, d, - selected_fe->dofs_per_vertex)]); - - for (vector::iterator i=mg_levels[level]->line_dofs.begin(); - i!=mg_levels[level]->line_dofs.end(); ++i) - { - Assert (*i != -1, ExcInternalError()); - *i = new_numbers[*i]; - }; - - for (vector::iterator i=mg_levels[level]->quad_dofs.begin(); - i!=mg_levels[level]->quad_dofs.end(); ++i) - { - Assert (*i != -1, ExcInternalError()); - *i = new_numbers[*i]; - }; - - for (vector::iterator i=mg_levels[level]->hex_dofs.begin(); - i!=mg_levels[level]->hex_dofs.end(); ++i) - { - Assert (*i != -1, ExcInternalError()); - *i = new_numbers[*i]; - }; -}; - -#endif - - - - -#if deal_II_dimension == 1 - -template <> -void MGDoFHandler<1>::reserve_space () { - const unsigned int dim = 1; - - Assert (selected_fe != 0, ExcNoFESelected()); - Assert (tria->n_levels() > 0, ExcInvalidTriangulation()); - - ////////////////////////// - // DESTRUCTION - - // delete all levels and set them up - // newly, since vectors are - // troublesome if you want to change - // their size - for (unsigned int i=0; in_levels(); ++i) - { - mg_levels.push_back (new DoFLevel<1>); - - mg_levels.back()->line_dofs = vector(tria->levels[i]->lines.lines.size() * - selected_fe->dofs_per_line, - -1); - }; - - // now allocate space for the - // vertices. To this end, we need - // to construct as many objects as - // there are vertices and let them - // allocate enough space for their - // vertex indices on the levels they - // live on. We need therefore to - // count to how many levels a cell - // belongs to, which we do by looping - // over all cells and storing the - // maximum and minimum level each - // vertex we pass by belongs to - mg_vertex_dofs.resize (tria->vertices.size()); - - vector min_level (tria->vertices.size(), tria->n_levels()); - vector max_level (tria->vertices.size(), 0); - - Triangulation::cell_iterator cell = tria->begin(), - endc = tria->end(); - for (; cell!=endc; ++cell) - for (unsigned int vertex=0; vertex::vertices_per_cell; - ++vertex) - { - const unsigned int vertex_index = cell->vertex_index(vertex); - if (min_level[vertex_index] > static_cast(cell->level())) - min_level[vertex_index] = cell->level(); - if (max_level[vertex_index] < static_cast(cell->level())) - max_level[vertex_index] = cell->level(); - }; - - // now allocate the needed space - for (unsigned int vertex=0; vertexvertices.size(); ++vertex) - { - Assert (min_level[vertex] < tria->n_levels(), ExcInternalError()); - Assert (max_level[vertex] >= min_level[vertex], ExcInternalError()); - - mg_vertex_dofs[vertex].init (min_level[vertex], - max_level[vertex], - selected_fe->dofs_per_vertex); - }; -}; - -#endif - - - -#if deal_II_dimension == 2 - -template <> -void MGDoFHandler<2>::reserve_space () { - const unsigned int dim = 2; - - Assert (selected_fe != 0, ExcNoFESelected()); - Assert (tria->n_levels() > 0, ExcInvalidTriangulation()); - - //////////////////////////// - // DESTRUCTION - - // delete all levels and set them up - // newly, since vectors are - // troublesome if you want to change - // their size - for (unsigned int i=0; in_levels(); ++i) - { - mg_levels.push_back (new DoFLevel<2>); - - mg_levels.back()->line_dofs = vector (tria->levels[i]->lines.lines.size() * - selected_fe->dofs_per_line, - -1); - mg_levels.back()->quad_dofs = vector (tria->levels[i]->quads.quads.size() * - selected_fe->dofs_per_quad, - -1); - }; - - - // now allocate space for the - // vertices. To this end, we need - // to construct as many objects as - // there are vertices and let them - // allocate enough space for their - // vertex indices on the levels they - // live on. We need therefore to - // count to how many levels a cell - // belongs to, which we do by looping - // over all cells and storing the - // maximum and minimum level each - // vertex we pass by belongs to - mg_vertex_dofs.resize (tria->vertices.size()); - - vector min_level (tria->vertices.size(), tria->n_levels()); - vector max_level (tria->vertices.size(), 0); - - Triangulation::cell_iterator cell = tria->begin(), - endc = tria->end(); - for (; cell!=endc; ++cell) - for (unsigned int vertex=0; vertex::vertices_per_cell; - ++vertex) - { - const unsigned int vertex_index = cell->vertex_index(vertex); - if (min_level[vertex_index] > static_cast(cell->level())) - min_level[vertex_index] = cell->level(); - if (max_level[vertex_index] < static_cast(cell->level())) - max_level[vertex_index] = cell->level(); - }; - - - // now allocate the needed space - for (unsigned int vertex=0; vertexvertices.size(); ++vertex) - { - Assert (min_level[vertex] < tria->n_levels(), ExcInternalError()); - Assert (max_level[vertex] >= min_level[vertex], ExcInternalError()); - - mg_vertex_dofs[vertex].init (min_level[vertex], - max_level[vertex], - selected_fe->dofs_per_vertex); - }; -}; - -#endif - - - - -#if deal_II_dimension == 3 - -template <> -void MGDoFHandler<3>::reserve_space () { - const unsigned int dim = 3; - - Assert (selected_fe != 0, ExcNoFESelected()); - Assert (tria->n_levels() > 0, ExcInvalidTriangulation()); - - //////////////////////////// - // DESTRUCTION - - // delete all levels and set them up - // newly, since vectors are - // troublesome if you want to change - // their size - for (unsigned int i=0; in_levels(); ++i) - { - mg_levels.push_back (new DoFLevel<3>); - - mg_levels.back()->line_dofs = vector (tria->levels[i]->lines.lines.size() * - selected_fe->dofs_per_line, - -1); - mg_levels.back()->quad_dofs = vector (tria->levels[i]->quads.quads.size() * - selected_fe->dofs_per_quad, - -1); - mg_levels.back()->hex_dofs = vector (tria->levels[i]->hexes.hexes.size() * - selected_fe->dofs_per_hex, - -1); - }; - - - // now allocate space for the - // vertices. To this end, we need - // to construct as many objects as - // there are vertices and let them - // allocate enough space for their - // vertex indices on the levels they - // live on. We need therefore to - // count to how many levels a cell - // belongs to, which we do by looping - // over all cells and storing the - // maximum and minimum level each - // vertex we pass by belongs to - mg_vertex_dofs.resize (tria->vertices.size()); - - vector min_level (tria->vertices.size(), tria->n_levels()); - vector max_level (tria->vertices.size(), 0); - - Triangulation::cell_iterator cell = tria->begin(), - endc = tria->end(); - for (; cell!=endc; ++cell) - for (unsigned int vertex=0; vertex::vertices_per_cell; - ++vertex) - { - const unsigned int vertex_index = cell->vertex_index(vertex); - if (min_level[vertex_index] > static_cast(cell->level())) - min_level[vertex_index] = cell->level(); - if (max_level[vertex_index] < static_cast(cell->level())) - max_level[vertex_index] = cell->level(); - }; - - - // now allocate the needed space - for (unsigned int vertex=0; vertexvertices.size(); ++vertex) - { - Assert (min_level[vertex] < tria->n_levels(), ExcInternalError()); - Assert (max_level[vertex] >= min_level[vertex], ExcInternalError()); - - mg_vertex_dofs[vertex].init (min_level[vertex], - max_level[vertex], - selected_fe->dofs_per_vertex); - }; -}; - -#endif - - -// explicite instantiations -template class MGDoFHandler; diff --git a/deal.II/deal.II/source/dofs/mg_dof_tools.cc b/deal.II/deal.II/source/dofs/mg_dof_tools.cc deleted file mode 100644 index 2fbb4da48c..0000000000 --- a/deal.II/deal.II/source/dofs/mg_dof_tools.cc +++ /dev/null @@ -1,42 +0,0 @@ -/* $Id$ */ - -#include -#include -#include -#include -#include -#include - - -template -void MGDoFTools::make_sparsity_pattern (const MGDoFHandler &mg_dof_handler, - const unsigned int level, - SparsityPattern &sparsity) -{ - Assert (sparsity.n_rows() == mg_dof_handler.n_dofs(level), - ExcDimensionMismatch (sparsity.n_rows(), mg_dof_handler.n_dofs(level))); - Assert (sparsity.n_cols() == mg_dof_handler.n_dofs(level), - ExcDimensionMismatch (sparsity.n_cols(), mg_dof_handler.n_dofs(level))); - - const unsigned int dofs_per_cell = mg_dof_handler.get_fe().dofs_per_cell; - vector dofs_on_this_cell(dofs_per_cell); - MGDoFHandler::cell_iterator cell = mg_dof_handler.begin(level), - endc = mg_dof_handler.end(level); - for (; cell!=endc; ++cell) - { - cell->get_mg_dof_indices (dofs_on_this_cell); - // make sparsity pattern for this cell - for (unsigned int i=0; i &, - const unsigned int, - SparsityPattern &); - diff --git a/deal.II/deal.II/source/grid/intergrid_map.cc b/deal.II/deal.II/source/grid/intergrid_map.cc index 4142e5e096..0299d848f3 100644 --- a/deal.II/deal.II/source/grid/intergrid_map.cc +++ b/deal.II/deal.II/source/grid/intergrid_map.cc @@ -2,10 +2,10 @@ #include #include -#include #include #include -#include +#include +#include #include #include diff --git a/deal.II/deal.II/source/multigrid/mg_base.cc b/deal.II/deal.II/source/multigrid/mg_base.cc index 25fa926268..06f3843ee3 100644 --- a/deal.II/deal.II/source/multigrid/mg_base.cc +++ b/deal.II/deal.II/source/multigrid/mg_base.cc @@ -1,6 +1,7 @@ // $Id$ -#include +#include +#include #include #include @@ -27,32 +28,14 @@ void print_vector(ostream& s, const VECTOR& v) -MGBase::~MGBase () -{}; - - -////////////////////////////////////////////////////////////////////// - - -MGSmootherBase::~MGSmootherBase() -{}; - - ////////////////////////////////////////////////////////////////////// -void -MGSmootherIdentity::smooth (const unsigned int, - Vector &, - const Vector &) const +MGBase::~MGBase () {}; -////////////////////////////////////////////////////////////////////// - - - MGBase::MGBase(const MGTransferBase &transfer, const unsigned minlevel, const unsigned maxlevel) diff --git a/deal.II/deal.II/source/multigrid/mg_dof_accessor.cc b/deal.II/deal.II/source/multigrid/mg_dof_accessor.cc index 9b0b78cf2e..df234c2e30 100644 --- a/deal.II/deal.II/source/multigrid/mg_dof_accessor.cc +++ b/deal.II/deal.II/source/multigrid/mg_dof_accessor.cc @@ -3,16 +3,16 @@ -#include -#include #include +#include +#include #include #include #include #include - - +#include +#include @@ -665,6 +665,7 @@ MGDoFCellAccessor<2>::face (const unsigned int i) const + #if deal_II_dimension == 3 template <> diff --git a/deal.II/deal.II/source/multigrid/mg_dof_handler.cc b/deal.II/deal.II/source/multigrid/mg_dof_handler.cc index d880f442b1..937ef0b0e1 100644 --- a/deal.II/deal.II/source/multigrid/mg_dof_handler.cc +++ b/deal.II/deal.II/source/multigrid/mg_dof_handler.cc @@ -2,16 +2,17 @@ /* Copyright W. Bangerth, University of Heidelberg, 1998 */ -#include #include -#include #include +#include +#include #include #include #include #include #include #include +#include #include diff --git a/deal.II/deal.II/source/multigrid/mg_dof_tools.cc b/deal.II/deal.II/source/multigrid/mg_dof_tools.cc index 2fbb4da48c..899113d635 100644 --- a/deal.II/deal.II/source/multigrid/mg_dof_tools.cc +++ b/deal.II/deal.II/source/multigrid/mg_dof_tools.cc @@ -1,10 +1,11 @@ /* $Id$ */ + #include -#include -#include +#include +#include #include -#include +#include #include diff --git a/deal.II/deal.II/source/multigrid/mg_smoother.cc b/deal.II/deal.II/source/multigrid/mg_smoother.cc index 9b6a5b4e97..f2f8bd4e40 100644 --- a/deal.II/deal.II/source/multigrid/mg_smoother.cc +++ b/deal.II/deal.II/source/multigrid/mg_smoother.cc @@ -1,18 +1,39 @@ /* $Id$ */ + #include -#include -#include +#include +#include #include #include -#include -#include +#include +#include +#include #include -#include #include + +////////////////////////////////////////////////////////////////////// + + +MGSmootherBase::~MGSmootherBase() +{}; + + +////////////////////////////////////////////////////////////////////// + + +void +MGSmootherIdentity::smooth (const unsigned int, + Vector &, + const Vector &) const +{}; + + + + #if deal_II_dimension == 1 MGSmoother::MGSmoother (const MGDoFHandler<1> &/*mg_dof*/, unsigned int steps) @@ -105,7 +126,6 @@ MGSmoother::set_zero_interior_boundary (const unsigned int level, ////////////////////////////////////////////////////////////////////// -#include // explicit instantiations // don't do the following instantiation in 1d, since there is a specialized @@ -116,16 +136,17 @@ template MGSmoother::MGSmoother (const MGDoFHandler&, unsigne template MGSmootherRelaxation -::MGSmootherRelaxation(const MGDoFHandler &mg_dof, - const MGMatrix >& matrix, - function_ptr relaxation, - unsigned int steps, - double omega); +::MGSmootherRelaxation(const MGDoFHandler &mg_dof, + const MGLevelObject > &matrix, + const function_ptr relaxation, + const unsigned int steps, + const double omega); + template MGSmootherRelaxation -::MGSmootherRelaxation(const MGDoFHandler &mg_dof, - const MGMatrix >& matrix, - function_ptr relaxation, - unsigned int steps, - double omega); +::MGSmootherRelaxation(const MGDoFHandler &mg_dof, + const MGLevelObject > &matrix, + const function_ptr relaxation, + const unsigned int steps, + const double omega); diff --git a/deal.II/deal.II/source/multigrid/multigrid.cc b/deal.II/deal.II/source/multigrid/multigrid.cc index 5ac6408d95..99f0ae2bbd 100644 --- a/deal.II/deal.II/source/multigrid/multigrid.cc +++ b/deal.II/deal.II/source/multigrid/multigrid.cc @@ -1,50 +1,46 @@ // $Id$ -// Copyright Guido Kanschat, Universitaet Heidelberg, 1999 + #include -#include -#include +#include +#include #include #include -#include -#include -#include +#include +#include +#include #include -#include + /* --------------------------------- MG ------------------------------------ */ template -MG::MG(const MGDoFHandler &mg_dof_handler, - ConstraintMatrix &constraints, - const MGMatrix > &level_matrices, - const MGTransferBase &transfer, - unsigned int minl, unsigned int maxl) +Multigrid::Multigrid (const MGDoFHandler &mg_dof_handler, + const ConstraintMatrix &constraints, + const MGLevelObject &level_sparsities, + const MGLevelObject > &level_matrices, + const MGTransferBase &transfer, + const unsigned int minlevel, + const unsigned int maxlevel) : - MGBase(transfer, minl, - min(mg_dof_handler.get_tria().n_levels()-1, maxl)), + MGBase(transfer, minlevel, + min(mg_dof_handler.get_tria().n_levels()-1, + maxlevel)), mg_dof_handler(&mg_dof_handler), + level_sparsities(&level_sparsities), level_matrices(&level_matrices), constraints(&constraints) -{ - // initialize the objects with - // their correct size - for (unsigned int level=minlevel; level<=maxlevel ; ++level) - { - solution[level].reinit(level_matrices[level].m()); - defect[level].reinit(level_matrices[level].m()); - }; -}; +{}; template void -MG::level_vmult (const unsigned int level, - Vector &result, - const Vector &u, - const Vector &/* rhs */) +Multigrid::level_vmult (const unsigned int level, + Vector &result, + const Vector &u, + const Vector &/* rhs */) { (*level_matrices)[level].vmult(result,u); result.scale(-1.); @@ -91,7 +87,7 @@ void MGTransferPrebuilt::build_matrices (const MGDoFHandler &mg_dof) // level which have children for (unsigned int level=0; level; +template class Multigrid; template void MGTransferPrebuilt::build_matrices (const MGDoFHandler &mg_dof); +template +void Multigrid::copy_to_mg (const Vector& src); + +template +void Multigrid::copy_to_mg (const Vector& src); + +template +void Multigrid::copy_from_mg (Vector& src) const; + +template +void Multigrid::copy_from_mg (Vector& src) const; + + diff --git a/deal.II/deal.II/source/numerics/mg_smoother.cc b/deal.II/deal.II/source/numerics/mg_smoother.cc deleted file mode 100644 index 9b6a5b4e97..0000000000 --- a/deal.II/deal.II/source/numerics/mg_smoother.cc +++ /dev/null @@ -1,131 +0,0 @@ -/* $Id$ */ - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include - - -#if deal_II_dimension == 1 - -MGSmoother::MGSmoother (const MGDoFHandler<1> &/*mg_dof*/, unsigned int steps) - : - steps(steps) -{ - Assert (false, ExcNotImplemented()); -}; - -#endif - - -#if deal_II_dimension > 1 - -template -MGSmoother::MGSmoother (const MGDoFHandler &mg_dof, unsigned int steps) - : - steps(steps) -{ - const unsigned int n_levels = mg_dof.get_tria().n_levels(); - - // allocate the right number of - // elements - interior_boundary_dofs.resize (n_levels-1); - - // use a temporary to store the - // indices. this allows to not allocate - // to much memory by later copying the - // content of this vector to its final - // destination - vector boundary_dofs; - - // temporary to hold the dof indices - // on a face between to levels - vector dofs_on_face (mg_dof.get_fe().dofs_per_face); - - for (unsigned int level=1; level::cell_iterator cell=mg_dof.begin(level); - cell != mg_dof.end(level); ++cell) - for (unsigned int face=0; face::faces_per_cell; ++face) - if ((cell->neighbor(face).state() == valid) && - (static_cast(cell->neighbor(face)->level()) - == level-1)) - { - // get indices of this face - cell->face(face)->get_mg_dof_indices (dofs_on_face); - // append them to the levelwise - // list - boundary_dofs.insert (boundary_dofs.end(), - dofs_on_face.begin(), - dofs_on_face.end()); - }; - - // now sort the list of interior boundary - // dofs and eliminate duplicates - sort (boundary_dofs.begin(), boundary_dofs.end()); - boundary_dofs.erase (unique (boundary_dofs.begin(), - boundary_dofs.end()), - boundary_dofs.end()); - - // now finally copy the result - // for this level its destination - interior_boundary_dofs[level-1] = boundary_dofs; - }; -}; - -#endif - -////////////////////////////////////////////////////////////////////// - -void -MGSmoother::set_zero_interior_boundary (const unsigned int level, - Vector &u) const -{ - if (level==0) - return; - else - for (vector::const_iterator p=interior_boundary_dofs[level-1].begin(); - p!=interior_boundary_dofs[level-1].end(); ++p) - u(*p) = 0; -}; - -////////////////////////////////////////////////////////////////////// - -#include - -// explicit instantiations -// don't do the following instantiation in 1d, since there is a specialized -// function there -#if deal_II_dimension > 1 -template MGSmoother::MGSmoother (const MGDoFHandler&, unsigned int); -#endif - -template -MGSmootherRelaxation -::MGSmootherRelaxation(const MGDoFHandler &mg_dof, - const MGMatrix >& matrix, - function_ptr relaxation, - unsigned int steps, - double omega); -template -MGSmootherRelaxation -::MGSmootherRelaxation(const MGDoFHandler &mg_dof, - const MGMatrix >& matrix, - function_ptr relaxation, - unsigned int steps, - double omega); - diff --git a/deal.II/deal.II/source/numerics/multigrid.cc b/deal.II/deal.II/source/numerics/multigrid.cc deleted file mode 100644 index 5ac6408d95..0000000000 --- a/deal.II/deal.II/source/numerics/multigrid.cc +++ /dev/null @@ -1,207 +0,0 @@ -// $Id$ -// Copyright Guido Kanschat, Universitaet Heidelberg, 1999 - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - - -/* --------------------------------- MG ------------------------------------ */ - -template -MG::MG(const MGDoFHandler &mg_dof_handler, - ConstraintMatrix &constraints, - const MGMatrix > &level_matrices, - const MGTransferBase &transfer, - unsigned int minl, unsigned int maxl) - : - MGBase(transfer, minl, - min(mg_dof_handler.get_tria().n_levels()-1, maxl)), - mg_dof_handler(&mg_dof_handler), - level_matrices(&level_matrices), - constraints(&constraints) -{ - // initialize the objects with - // their correct size - for (unsigned int level=minlevel; level<=maxlevel ; ++level) - { - solution[level].reinit(level_matrices[level].m()); - defect[level].reinit(level_matrices[level].m()); - }; -}; - - - -template -void -MG::level_vmult (const unsigned int level, - Vector &result, - const Vector &u, - const Vector &/* rhs */) -{ - (*level_matrices)[level].vmult(result,u); - result.scale(-1.); -} - - - - - - -/* ----------------------------- MGTransferPrebuilt ------------------------ */ - - - -MGTransferPrebuilt::~MGTransferPrebuilt () -{}; - - - -template -void MGTransferPrebuilt::build_matrices (const MGDoFHandler &mg_dof) -{ - const unsigned int n_levels = mg_dof.get_tria().n_levels(); - const unsigned int dofs_per_cell = mg_dof.get_fe().dofs_per_cell; - - // reset the size of the array of - // matrices - prolongation_sparsities.clear (); - prolongation_matrices.clear (); - prolongation_sparsities.reserve (n_levels-1); - prolongation_matrices.reserve (n_levels-1); - - - // two fields which will store the - // indices of the multigrid dofs - // for a cell and one of its children - vector dof_indices_mother (dofs_per_cell); - vector dof_indices_child (dofs_per_cell); - - // for each level: first build the sparsity - // pattern of the matrices and then build the - // matrices themselves. note that we only - // need to take care of cells on the coarser - // level which have children - for (unsigned int level=0; level::cell_iterator cell=mg_dof.begin(level); - cell != mg_dof.end(level); ++cell) - if (cell->has_children()) - { - cell->get_mg_dof_indices (dof_indices_mother); - - for (unsigned int child=0; - child::children_per_cell; ++child) - { - // set an alias to the - // prolongation matrix for - // this child - const FullMatrix &prolongation - = mg_dof.get_fe().prolongate(child); - - cell->child(child)->get_mg_dof_indices (dof_indices_child); - - // now tag the entries in the - // matrix which will be used - // for this pair of mother/child - for (unsigned int i=0; i()); - prolongation_matrices[level].reinit (prolongation_sparsities[level]); - - - // now actually build the matrices - for (MGDoFHandler::cell_iterator cell=mg_dof.begin(level); - cell != mg_dof.end(level); ++cell) - if (cell->has_children()) - { - cell->get_mg_dof_indices (dof_indices_mother); - - for (unsigned int child=0; - child::children_per_cell; ++child) - { - // set an alias to the - // prolongation matrix for - // this child - const FullMatrix &prolongation - = mg_dof.get_fe().prolongate(child); - - cell->child(child)->get_mg_dof_indices (dof_indices_child); - - // now set the entries in the - // matrix - for (unsigned int i=0; i &dst, - const Vector &src) const -{ - Assert ((to_level >= 1) && (to_level<=prolongation_matrices.size()), - ExcIndexRange (to_level, 1, prolongation_matrices.size()+1)); - - prolongation_matrices[to_level-1].vmult (dst, src); -}; - - - -void MGTransferPrebuilt::restrict_and_add (const unsigned int from_level, - Vector &dst, - const Vector &src) const -{ - Assert ((from_level >= 1) && (from_level<=prolongation_matrices.size()), - ExcIndexRange (from_level, 1, prolongation_matrices.size()+1)); - - prolongation_matrices[from_level-1].Tvmult_add (dst, src); -}; - - - -// explicit instatations -template class MG; -template -void MGTransferPrebuilt::build_matrices (const MGDoFHandler &mg_dof); - diff --git a/deal.II/lac/include/lac/mgbase.h b/deal.II/lac/include/lac/mgbase.h deleted file mode 100644 index 6012ba2b64..0000000000 --- a/deal.II/lac/include/lac/mgbase.h +++ /dev/null @@ -1,687 +0,0 @@ -/*---------------------------- mgbase.h ---------------------------*/ -/* $Id$ */ -#ifndef __mgbase_H -#define __mgbase_H -/*---------------------------- mgbase.h ---------------------------*/ - - - -#include -#include -#include -#include - -#include - - - -/** - * Abstract base class for coarse grid solvers. The interface of a - * function call operator is defined to execute coarse grid solution - * in a derived class. - * - * @author Guido Kanschat, 1999 - */ -class MGCoarseGridSolver : public Subscriptor -{ - public: - /** - * Virtual destructor. Does - * nothing in particular, but - * needs to be declared anyhow. - */ - virtual ~MGCoarseGridSolver(); - - /** - * Coarse grid solving method. - * This is only the interface for - * a function defined in a derived - * class like - * #MGCoarseGridLACIteration#. - * - * Note that the information - * about the matrix is removed to - * that class. - */ - virtual void operator() (const unsigned int level, - Vector &dst, - const Vector &src) const = 0; -}; - - - - -/** - * Abstract base class for multigrid smoothers. - * In fact, this class only provides the interface of the smoothing function. - * Using #deal.II# grid handling, #MGSmoother# is a good point to start. - * - * @author Wolfgang Bangerth, Guido Kanschat, 1999 - */ -class MGSmootherBase : public Subscriptor -{ - public: - /** - * Virtual destructor. Does - * nothing in particular, but - * needs to be declared anyhow. - */ - virtual ~MGSmootherBase(); - - /** - * Smooth the residual of #u# on - * the given level. If $S$ is the - * smoothing operator, then this - * function should do the - * following operation: - * $u += S (rhs - Au)$, where #u# and - * #rhs# are the input parameters. - * - * This function should keep the - * interior level boundary - * values, so you may want to - * call #set_zero_interior_boundary# - * in the derived class - * #MGSmoother# somewhere in your - * derived function, or another - * function doing similar things. - */ - virtual void smooth (const unsigned int level, - Vector &u, - const Vector &rhs) const = 0; -}; - - - - - -/** - * Smoother doing nothing. This class is not useful for many applications other - * than for testing some multigrid procedures. Also some applications might - * get convergence without smoothing and then this class brings you the - * cheapest possible multigrid. - * - * @author Guido Kanschat, 1999 - */ -class MGSmootherIdentity : public MGSmootherBase -{ - public: - /** - * Implementation of the - * interface in #MGSmootherBase#. - * This function does nothing, - * which by comparison with the - * definition of this function - * means that the the smoothing - * operator equals the null - * operator. - */ - virtual void smooth (const unsigned int level, - Vector &u, - const Vector &rhs) const; -}; - - - - -/** - * Base class used to declare the operations needed by a concrete class - * implementing prolongation and restriction of vectors in the multigrid - * context. This class is an abstract one and has no implementations of - * possible algorithms for these operations. - * - * @author Wolfgang Bangerth, Guido Kanschat, 1999 - */ -class MGTransferBase : public Subscriptor -{ - public: - /** - * Destructor. Does nothing here, but - * needs to be declared virtual anyway. - */ - virtual ~MGTransferBase(); - - /** - * Prolongate a vector from level - * #to_level-1# to level - * #to_level#. The previous - * content of #dst# is - * overwritten. - * - * #src# is assumed to be a vector with - * as many elements as there are degrees - * of freedom on the coarser level of - * the two involved levels, while #src# - * shall have as many elements as there - * are degrees of freedom on the finer - * level. - */ - virtual void prolongate (const unsigned int to_level, - Vector &dst, - const Vector &src) const = 0; - - /** - * Restrict a vector from level - * #from_level# to level - * #from_level-1# and add this - * restriction to - * #dst#. Obviously, if the - * refined region on level - * #from_level# is smaller than - * that on level #from_level-1#, - * some degrees of freedom in - * #dst# are not covered and will - * not be altered. For the other - * degress of freedom, the result - * of the restriction is added. - * - * #src# is assumed to be a vector with - * as many elements as there are degrees - * of freedom on the finer level of - * the two involved levels, while #src# - * shall have as many elements as there - * are degrees of freedom on the coarser - * level. - */ - virtual void restrict_and_add (const unsigned int from_level, - Vector &dst, - const Vector &src) const = 0; -}; - - - -/** - * An array with a vector for each level. The purpose of this class - * is mostly to allow access by level number, even if the lower levels - * are not used and therefore have no vector at all; this is done by - * simply shifting the given index by the minimum level we have - * stored. - * - * @author Guido Kanschat, 1999 - */ -template > -class MGVector : public Subscriptor -{ - public: - /** - * Constructor allowing to - * initialize the number of - * levels. - */ - MGVector (const unsigned int minlevel, - const unsigned int maxlevel); - - /** - * Access vector on level #level#. - */ - VECTOR & operator[] (const unsigned int level); - - /** - * Access vector on level - * #level#. Constant version. - */ - const VECTOR & operator[] (const unsigned int level) const; - - /** - * Call #clear# on all vectors - * stored by this object. Note - * that if #VECTOR==Vector#, - * #clear# will set all entries - * to zero, while if - * #VECTOR==std::vector#, - * #clear# deletes the elements - * of the vectors. This class - * might therefore not be useful - * for STL vectors. - */ - void clear(); - - private: - /** - * Level of first component. - */ - const unsigned int minlevel; - - /** - * Array of the vectors to be held. - */ - vector vectors; -}; - - - - -/** - * An array of matrices for each level. - * This class provides the functionality of #vector# combined - * with a subscriptor for smart pointers. - * @author Guido Kanschat, 1999 - */ -template > -class MGMatrix : public Subscriptor, - public vector -{ - public: - /** - * Constructor allowing to initialize the number of levels. - */ - MGMatrix(unsigned int minlevel, unsigned int maxlevel); - - /** - * Safe access operator. - */ - MATRIX& operator[](unsigned int); - - /** - * Safe access operator. - */ - const MATRIX& operator[](unsigned int) const; - - private: - /** - * Level of first component. - */ - unsigned int minlevel; -}; - - - - -/** - * Coarse grid solver using LAC iterative methods. - * This is a little wrapper, transforming a triplet of iterative - * solver, matrix and preconditioner into a coarse grid solver. - * - * The type of the matrix (i.e. the template parameter #MATRIX#) - * should be derived from #Subscriptor# to allow for the use of a - * smart pointer to it. - * - * @author Guido Kanschat, 1999 - */ -template -class MGCoarseGridLACIteration : public MGCoarseGridSolver -{ - public: - /** - * Constructor. - * Store solver, matrix and - * preconditioning method for later - * use. - */ - MGCoarseGridLACIteration (SOLVER &, - const MATRIX &, - const PRECOND &); - - /** - * Implementation of the abstract - * function. - * Calls the solver method with - * matrix, vectors and - * preconditioner. - */ - virtual void operator() (const unsigned int level, - Vector &dst, - const Vector &src) const; - private: - /** - * Reference to the solver. - */ - SOLVER& solver; - - /** - * Reference to the matrix. - */ - const SmartPointer matrix; - - /** - * Reference to the preconditioner. - */ - const PRECOND& precondition; -}; - - - - -/** - * Basic class for preconditioning by multigrid. - * - * The functionality of the multigrid method is restricted to defect - * correction. It is {\bf not} iterative and the start solution is - * always zero. Since by this $u^E_l$ and $u^A_l$ (see report on - * multigrid) are always zero, restriction is simplified a lot and - * maybe even the seam condition on grids is oblivious. Still, I am - * not sure that these restrictions on the class might cause numerical - * inefficiencies. - * - * The function #precondition# is the actual multigrid method and - * makes use of several operations to be implemented in derived - * classes. It takes a defect in #src# and the result is the multigrid - * preconditioned defect in #dst#. - * - * @author Guido Kanschat, 1999 - */ -class MGBase : public Subscriptor -{ - private: - /** - * Copy constructor. Made private - * to prevent use. - */ - MGBase(const MGBase&); - - /** - * Copy operator. Made private - * to prevent use. - */ - const MGBase& operator=(const MGBase&); - - public: - /** - * Constructor, subject to change. - */ - MGBase (const MGTransferBase &transfer, - const unsigned int minlevel, - const unsigned int maxlevel); - - /** - * Virtual destructor. - */ - virtual ~MGBase(); - - /** - * Execute one step of the - * v-cycle algorithm. This - * function assumes, that the - * vector #d# is properly filled - * with the residual in the outer - * defect correction - * scheme. After execution of - * #vcycle()#, the result is in - * the vector #s#. We propose to - * write functions #copy_from_mg# - * and #copy_to_mg# in derived - * classes of #MGBase# to access - * #d# and #s#. - * - * The actual work for this - * function is done in - * #level_mgstep#. - */ - void vcycle(const MGSmootherBase &pre_smooth, - const MGSmootherBase &post_smooth, - const MGCoarseGridSolver &cgs); - - /** - * Exception. - */ - DeclException2(ExcSwitchedLevels, int, int, - << "minlevel and maxlevel switched, should be: " - << arg1 << "<=" << arg2); - - protected: - /** - * Highest level of cells. - */ - unsigned int maxlevel; - - /** - * Level for coarse grid solution. - */ - unsigned int minlevel; - - /** - * Auxiliary vector. - */ - MGVector > defect; - - /** - * Auxiliary vector. - */ - MGVector > solution; - - /** - * Prolongation and restriction object. - */ - SmartPointer transfer; - - /** - * The actual v-cycle multigrid method. - * #level# is the level the - * function should work on. It - * will usually be called for the - * highest level from outside, - * but will then call itself - * recursively for #level-1#, - * unless we are on #minlevel# - * where instead of recursing - * deeper, the coarse grid solver - * is used to solve the matrix of - * this level. - */ - void level_mgstep (const unsigned int level, - const MGSmootherBase &pre_smooth, - const MGSmootherBase &post_smooth, - const MGCoarseGridSolver &cgs); - - /** - * Apply negative #vmult# on all - * cells of a level. - * This is implemented in a - * derived class. - */ - virtual void level_vmult (const unsigned int level, - Vector &dst, - const Vector &src, - const Vector &rhs) = 0; - - private: - /** - * Auxiliary vector. - */ - Vector t; -}; - - - -/** - * Multi-level preconditioner. - * Here, we collect all information needed for multi-level preconditioning - * and provide the standard interface for LAC iterative methods. - * - * The template parameter class #MG# is required to inherit #MGBase#. - * Furthermore, it needs functions #void copy_to_mg(const VECTOR&)# - * to store #src# in the right hand side of the multi-level method and - * #void copy_from_mg(VECTOR&)# to store the result of the v-cycle in #dst#. - * - * @author Guido Kanschat, 1999 - */ -template > -class PreconditionMG -{ - public: - /** - * Constructor. - * Arguments are the multigrid object, - * pre-smoother, post-smoother and - * coarse grid solver. - */ - PreconditionMG(MG &mg, - const MGSmootherBase &pre, - const MGSmootherBase &post, - const MGCoarseGridSolver &coarse); - - /** - * Preconditioning operator, - * calling the #vcycle# function - * of the #MG# object passed to - * the constructor. - * - * This is the operator used by - * LAC iterative solvers. - */ - void operator() (VECTOR &dst, - const VECTOR &src) const; - - private: - /** - * The multigrid object. - */ - SmartPointer multigrid; - - /** - * The pre-smoothing object. - */ - SmartPointer pre; - - /** - * The post-smoothing object. - */ - SmartPointer post; - - /** - * The coarse grid solver. - */ - SmartPointer coarse; -}; - - - - -/* ------------------------------------------------------------------- */ - -template -MGVector::MGVector(const unsigned int min, - const unsigned int max) - : - minlevel(min), - vectors(max-min+1) -{}; - - - -template -VECTOR & -MGVector::operator[] (const unsigned int i) -{ - Assert((i>=minlevel) && (i -const VECTOR & -MGVector::operator[] (const unsigned int i) const -{ - Assert((i>=minlevel) && (i -void -MGVector::clear() -{ - typename vector::iterator v; - for (v = vectors.begin(); v != vectors.end(); ++v) - v->clear(); - -} - - -/* ------------------------------------------------------------------- */ - - -template -MGMatrix::MGMatrix(unsigned int min, unsigned int max) - : - vector(max-min+1), - minlevel(min) -{} - - -template -MATRIX& -MGMatrix::operator[](unsigned int i) -{ - Assert((i>=minlevel)&&(i::operator[](i-minlevel); -}; - - - -template -const MATRIX& -MGMatrix::operator[](unsigned int i) const -{ - Assert((i>=minlevel)&&(i::operator[](i-minlevel); -}; - - - -/* ------------------ Functions for MGCoarseGridLACIteration ------------ */ - - -template -MGCoarseGridLACIteration:: -MGCoarseGridLACIteration(SOLVER &s, - const MATRIX &m, - const PRECOND &p) - : - solver(s), - matrix(&m), - precondition(p) -{}; - - - -template -void -MGCoarseGridLACIteration:: -operator() (const unsigned int /* level */, - Vector &dst, - const Vector &src) const -{ - solver.solve(*matrix, dst, src, precondition); -} - - - - - - -/* ------------------------------------------------------------------- */ - - -template -PreconditionMG::PreconditionMG(MG &mg, - const MGSmootherBase &pre, - const MGSmootherBase &post, - const MGCoarseGridSolver &coarse) - : - multigrid(&mg), - pre(&pre), - post(&post), - coarse(&coarse) -{} - - - -template -void -PreconditionMG::operator() (VECTOR &dst, - const VECTOR &src) const -{ - multigrid->copy_to_mg(src); - multigrid->vcycle(*pre, *post, *coarse); - multigrid->copy_from_mg(dst); -} - - - -/*---------------------------- mgbase.h ---------------------------*/ -/* end of #ifndef __mgbase_H */ -#endif -/*---------------------------- mgbase.h ---------------------------*/ diff --git a/deal.II/lac/source/mgbase.cc b/deal.II/lac/source/mgbase.cc deleted file mode 100644 index 25fa926268..0000000000 --- a/deal.II/lac/source/mgbase.cc +++ /dev/null @@ -1,142 +0,0 @@ -// $Id$ - -#include -#include -#include - - - -//TODO: this function is only for debugging purposes and should be removed sometimes -template -static -void print_vector(ostream& s, const VECTOR& v) -{ - const unsigned int n = (unsigned int)(sqrt(v.size())+.3); - unsigned int k=0; - - // write the vector entries in a - // kind of square - for (unsigned int i=0;i &, - const Vector &) const -{}; - - - -////////////////////////////////////////////////////////////////////// - - - -MGBase::MGBase(const MGTransferBase &transfer, - const unsigned minlevel, - const unsigned maxlevel) - : - maxlevel(maxlevel), - minlevel(minlevel), - defect(minlevel,maxlevel), - solution(minlevel,maxlevel), - transfer(&transfer) -{ - Assert(minlevel <= maxlevel, - ExcSwitchedLevels(minlevel, maxlevel)); -}; - - - -void -MGBase::vcycle(const MGSmootherBase &pre_smooth, - const MGSmootherBase &post_smooth, - const MGCoarseGridSolver &coarse_grid_solver) -{ - level_mgstep(maxlevel, pre_smooth, post_smooth, coarse_grid_solver); -}; - - - -void -MGBase::level_mgstep(const unsigned int level, - const MGSmootherBase &pre_smooth, - const MGSmootherBase &post_smooth, - const MGCoarseGridSolver &coarse_grid_solver) -{ - solution[level] = 0.; - - if (level == minlevel) - { - coarse_grid_solver(level, solution[level], defect[level]); - return; - } - - // smoothing of the residual by modifying s - pre_smooth.smooth(level, solution[level], defect[level]); - // t = d-As - t.reinit(solution[level].size()); - level_vmult(level, t, solution[level], defect[level]); - - // make t rhs of lower level -//TODO: this function adds the restricted t to defect[level-1]. -//TODO: why don't we have to clear it before? - transfer->restrict_and_add (level, defect[level-1], t); - - // do recursion - level_mgstep(level-1, pre_smooth, post_smooth, coarse_grid_solver); - - // reset size of the auxiliary - // vector, since it has been - // resized in the recursive call to - // level_mgstep directly above - t.reinit(solution[level].size()); - - // do coarse grid correction - transfer->prolongate(level, t, solution[level-1]); - solution[level] += t; - - // smoothing (modify solution again) -//TODO: what happens here? smooth overwrites the solution[level], -//TODO: so the previous two statements should have no effect. No? - post_smooth.smooth(level, solution[level], defect[level]); -}; - - -////////////////////////////////////////////////////////////////////// - -MGCoarseGridSolver::~MGCoarseGridSolver() -{}; - - - -////////////////////////////////////////////////////////////////////// - -MGTransferBase::~MGTransferBase() -{}; - - - - - diff --git a/tests/big-tests/multigrid/multigrid.cc b/tests/big-tests/multigrid/multigrid.cc index e48fd7e987..7769168471 100644 --- a/tests/big-tests/multigrid/multigrid.cc +++ b/tests/big-tests/multigrid/multigrid.cc @@ -2,712 +2,514 @@ /* Copyright W. Bangerth, University of Heidelberg, 1998 */ -#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 +#include +#include +#include +#include +#include -#include -#include -#include -#include -#include - - - - - -template -class PoissonEquation : public Equation { - public: - PoissonEquation (const Function &rhs) : - Equation(1), - right_hand_side (rhs) {}; - - virtual void assemble (FullMatrix &cell_matrix, - Vector &rhs, - const FEValues &fe_values, - const DoFHandler::cell_iterator &cell) const; - virtual void assemble (FullMatrix &cell_matrix, - const FEValues &fe_values, - const DoFHandler::cell_iterator &cell) const; - virtual void assemble (Vector &rhs, - const FEValues &fe_values, - const DoFHandler::cell_iterator &cell) const; - protected: - const Function &right_hand_side; -}; - +#include +#include +#include +#include +#include +#include +#include +#include template -class PoissonProblem { +class LaplaceProblem +{ public: - /** - * Declare a data type which denotes a - * mapping between a boundary indicator - * and the function denoting the boundary - * values on this part of the boundary. - * Only one boundary function may be given - * for each boundary indicator, which is - * guaranteed by the #map# data type. - * - * See the general documentation of this - * class for more detail. - */ - typedef map*> FunctionMap; - /** - * Typdedef an iterator which assembles - * matrices and vectors. - */ - typedef TriaActiveIterator > active_assemble_iterator; - - PoissonProblem (unsigned int order); - virtual ~PoissonProblem (); + LaplaceProblem (); + ~LaplaceProblem (); + void run (); - void clear (); - void create_new (); + private: + void setup_system (); + void assemble_system (); void solve (); - - /** - * Initiate the process of assemblage of - * vectors and system matrix. Use the - * given equation object and the given - * quadrature formula. Also use the list - * of dirichlet boundary value functions - * (by default, no dirichlet bc are assumed - * which means that all bc are included - * into the weak formulation). - * - * For what exactly happens here, refer to - * the general doc of this class. - */ - virtual void assemble (const Equation &equation, - const Quadrature &q, - const UpdateFlags update_flags, - const FunctionMap &dirichlet_bc = FunctionMap()); - - int run (unsigned int level); - void print_history (string filename) const; - - /** - * Exception - */ - DeclException0 (ExcNoTriaSelected); - - protected: - Triangulation *tria; - MGDoFHandler *dof; - - - /** - * Sparsity pattern of the system matrix. - */ - SparsityPattern system_sparsity; - - /** - * System matrix. - */ - SparseMatrix system_matrix; + void refine_grid (); + void output_results (const unsigned int cycle) const; - /** - * Vector storing the right hand side. - */ - Vector right_hand_side; + Triangulation triangulation; + MGDoFHandler mg_dof_handler; - /** - * Solution vector. - */ - Vector solution; + FEQ1 fe; - /** - * List of constraints introduced by - * hanging nodes. - */ - ConstraintMatrix constraints; + ConstraintMatrix hanging_node_constraints; - Function *rhs; - Function *boundary_values; + SparsityPattern global_sparsity_pattern; + SparseMatrix global_system_matrix; - vector l1_error, l2_error, linfty_error, h1_seminorm_error, h1_error; - vector n_dofs; - - unsigned int order; + MGLevelObject level_sparsity_patterns; + MGLevelObject > level_system_matrices; + + Vector solution; + Vector system_rhs; }; - - -/** - Right hand side constructed such that the exact solution is - $sin(2 pi x) + sin(2 pi y)$ - */ template -class RHSPoly : public Function { +class Coefficient : public Function +{ public: - /** - * Return the value of the function - * at the given point. - */ - virtual double value (const Point &p, - const unsigned int component) const; + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual void value_list (const vector > &points, + vector &values, + const unsigned int component = 0) const; }; template -class Solution : public Function { - public: - /** - * Return the value of the function - * at the given point. - */ - virtual double value (const Point &p, - const unsigned int component) const; - /** - * Return the gradient of the function - * at the given point. - */ - virtual Tensor<1,dim> gradient (const Point &p, - const unsigned int component) const; +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 vector > &points, + vector &values, + const unsigned int component) const +{ + const unsigned int n_points = points.size(); -template <> -double RHSPoly<2>::value (const Point<2> &p, - const unsigned int component) const { - Assert (component==0, ExcIndexRange (component, 0, 1)); + Assert (values.size() == n_points, + ExcVectorHasWrongSize (values.size(), n_points)); + + Assert (component == 0, + ExcWrongComponent (component, 1)); - const double x = p(0), - y = p(1); - const double pi= 3.1415926536; - return 4*pi*pi*(sin(2*pi*x)+sin(2*pi*y)); + for (unsigned int i=0; i -double Solution<2>::value (const Point<2> &p, - const unsigned int component) const { - Assert (component==0, ExcIndexRange (component, 0, 1)); - const double x = p(0), - y = p(1); - const double pi= 3.1415926536; - return sin(2*pi*x)+sin(2*pi*y); +class MGSmootherLAC : public MGSmootherBase +{ + private: + SmartPointer > >matrices; + public: + MGSmootherLAC(MGLevelObject >&); + + virtual void smooth (const unsigned int level, + Vector &u, + const Vector &rhs) const; }; -template <> -Tensor<1,2> Solution<2>::gradient (const Point<2> &p, - const unsigned int component) const { - Assert (component==0, ExcIndexRange (component, 0, 1)); +MGSmootherLAC::MGSmootherLAC(MGLevelObject >& matrix) + : + matrices(&matrix) +{} - const double x = p(0), - y = p(1); - const double pi= 3.1415926536; - return Point<2> (2*pi*cos(2*pi*x), - 2*pi*cos(2*pi*y)); -}; - +void +MGSmootherLAC::smooth (const unsigned int level, + Vector &u, + const Vector &rhs) const +{ + SolverControl control(2,1.e-300,false,false); + PrimitiveVectorMemory<> mem; + SolverRichardson<> rich(control, mem); + PreconditionRelaxation<> + prec((*matrices)[level], &SparseMatrix ::template precondition_SSOR, 1.); + rich.solve((*matrices)[level], u, rhs, prec); +} -template <> -void PoissonEquation<2>::assemble (FullMatrix &cell_matrix, - Vector &rhs, - const FEValues<2> &fe_values, - const DoFHandler<2>::cell_iterator &) const { - for (unsigned int point=0; point +LaplaceProblem::LaplaceProblem () : + mg_dof_handler (triangulation) +{}; template -void PoissonEquation::assemble (FullMatrix &, - const FEValues &, - const DoFHandler::cell_iterator &) const { - Assert (false, ExcPureVirtualFunctionCalled()); +LaplaceProblem::~LaplaceProblem () +{ + mg_dof_handler.clear (); }; template -void PoissonEquation::assemble (Vector &, - const FEValues &, - const DoFHandler::cell_iterator &) const { - Assert (false, ExcPureVirtualFunctionCalled()); -}; - +void LaplaceProblem::setup_system () +{ + mg_dof_handler.distribute_dofs (fe); + hanging_node_constraints.clear (); + DoFTools::make_hanging_node_constraints (mg_dof_handler, + hanging_node_constraints); + hanging_node_constraints.close (); + global_sparsity_pattern.reinit (mg_dof_handler.DoFHandler::n_dofs(), + mg_dof_handler.DoFHandler::n_dofs(), + mg_dof_handler.max_couplings_between_dofs()); + DoFTools::make_sparsity_pattern (mg_dof_handler, global_sparsity_pattern); + hanging_node_constraints.condense (global_sparsity_pattern); + global_sparsity_pattern.compress(); + global_system_matrix.reinit (global_sparsity_pattern); + solution.reinit (mg_dof_handler.DoFHandler::n_dofs()); + system_rhs.reinit (mg_dof_handler.DoFHandler::n_dofs()); + const unsigned int n_levels = triangulation.n_levels(); + level_system_matrices.resize (0, n_levels); + level_sparsity_patterns.resize (0, n_levels); + + for (unsigned int level=0; level -PoissonProblem::PoissonProblem (unsigned int order) : - tria(0), dof(0), rhs(0), - boundary_values(0), order(order) {}; +void LaplaceProblem::assemble_system () +{ + const Coefficient coefficient; + QGauss2 quadrature_formula; + FEValues fe_values (fe, quadrature_formula, + UpdateFlags(update_values | + update_gradients | + update_q_points | + update_JxW_values)); -template -PoissonProblem::~PoissonProblem () -{ - clear (); -}; + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.n_quadrature_points; + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + vector local_dof_indices (dofs_per_cell); -template -void PoissonProblem::clear () { - if (dof != 0) { - delete dof; - dof = 0; - }; - - if (tria != 0) { - delete tria; - tria = 0; - }; - - if (rhs != 0) - { - delete rhs; - rhs = 0; - }; + // FIX + vector coefficient_values (n_q_points, 1.0); - if (boundary_values != 0) + // not only active cells + MGDoFHandler::cell_iterator cell = mg_dof_handler.begin(), + endc = mg_dof_handler.end(); + for (; cell!=endc; ++cell) { - delete boundary_values; - boundary_values = 0; + cell_matrix.clear (); + cell_rhs.clear (); + + fe_values.reinit (cell); + const FullMatrix + & shape_values = fe_values.get_shape_values(); + const vector > > + & shape_grads = fe_values.get_shape_grads(); + const vector + & JxW_values = fe_values.get_JxW_values(); + const vector > + & q_points = fe_values.get_quadrature_points(); + + // FIX +// coefficient.value_list (q_points, coefficient_values); + + for (unsigned int q_point=0; q_pointget_mg_dof_indices (local_dof_indices); + const unsigned int level = cell->level(); + for (unsigned int i=0; iactive()) + { + 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, +// global_system_matrix, +// solution, +// system_rhs); }; - template -void PoissonProblem::create_new () { - clear (); - - tria = new Triangulation(); - dof = new MGDoFHandler (*tria); -}; +void LaplaceProblem::solve () +{ + { + SolverControl solver_control (1000, 1e-12); + PrimitiveVectorMemory<> vector_memory; + SolverCG<> cg (solver_control, vector_memory); + SolverControl coarse_grid_solver_control (1000, 1e-12); + PrimitiveVectorMemory<> coarse_grid_vector_memory; + + SolverCG<> coarse_grid_cg (coarse_grid_solver_control, + coarse_grid_vector_memory); + +// PreconditionRelaxation<> +// coarse_grid_solver_preconditioner(level_system_matrices[level_system_matrices.get_minlevel()], +// &SparseMatrix::template precondition_SSOR, +// 1.2); + PreconditionIdentity coarse_grid_solver_preconditioner; + + MGCoarseGridLACIteration, SparseMatrix, PreconditionIdentity> + coarse_grid_solver (coarse_grid_cg, + level_system_matrices[level_system_matrices.get_minlevel()], + coarse_grid_solver_preconditioner); + + MGSmootherLAC smoother (level_system_matrices); + MGTransferPrebuilt grid_transfer; + grid_transfer.build_matrices (mg_dof_handler); + + Multigrid<2> multigrid (mg_dof_handler, + hanging_node_constraints, + level_sparsity_patterns, + level_system_matrices, + grid_transfer); + + PreconditionMG > + mg_precondition (multigrid, smoother, smoother, coarse_grid_solver); -template -void PoissonProblem::assemble (const Equation &equation, - const Quadrature &quadrature, - const UpdateFlags update_flags, - const FunctionMap &dirichlet_bc) { - Assert ((tria!=0) && (dof!=0), ExcNoTriaSelected()); - - system_sparsity.reinit (dof->DoFHandler::n_dofs(), - dof->DoFHandler::n_dofs(), - dof->max_couplings_between_dofs()); - right_hand_side.reinit (dof->DoFHandler::n_dofs()); - - // make up sparsity pattern and - // compress with constraints - constraints.clear (); - DoFTools::make_hanging_node_constraints (*dof, constraints); - constraints.close (); - DoFTools::make_sparsity_pattern (*dof, system_sparsity); - constraints.condense (system_sparsity); - - MGTransferPrebuilt p; - p.build_matrices (*dof); -// MGSmoother smoother(*dof); - + solution.clear (); + cg.solve (global_system_matrix, solution, system_rhs, + mg_precondition); - // reinite system matrix - system_matrix.reinit (system_sparsity); - // reinit solution vector, preset - // with zeroes. - solution.reinit (dof->DoFHandler::n_dofs()); + cout << " MG Outer iterations: " << solver_control.last_step() + << endl; + + cout << " MG Total inner iterations: " << coarse_grid_solver_control.last_step() + << endl; + }; - // create assembler - Assembler::AssemblerData data (*dof, - true, true, //assemble matrix and rhs - system_matrix, - right_hand_side, - quadrature, - update_flags); - active_assemble_iterator assembler (tria, - tria->begin_active()->level(), - tria->begin_active()->index(), - &data); - // loop over all cells, fill matrix and rhs - do { - assembler->assemble (equation); - } - while ((++assembler).state() == valid); - - // condense system matrix in-place - constraints.condense (system_matrix); - - // condense right hand side in-place - constraints.condense (right_hand_side); - - // apply Dirichlet bc as described - // in the docs - map boundary_value_list; - - for (FunctionMap::const_iterator bc=dirichlet_bc.begin(); - bc != dirichlet_bc.end(); ++bc) - VectorTools::interpolate_boundary_values (*dof, - bc->first, - *bc->second, - boundary_value_list); - MatrixTools::apply_boundary_values (boundary_value_list, - system_matrix, solution, - right_hand_side); + SolverControl solver_control (1000, 1e-12); + PrimitiveVectorMemory<> vector_memory; + SolverCG<> cg (solver_control, vector_memory); + + PreconditionRelaxation<> + preconditioner(global_system_matrix, + &SparseMatrix::template precondition_SSOR, + 1.2); + + solution.clear (); + cg.solve (global_system_matrix, solution, system_rhs, + preconditioner); + + cout << " CG Outer iterations: " << solver_control.last_step() + << endl; + }; + + hanging_node_constraints.distribute (solution); +}; + + +template +void LaplaceProblem::refine_grid () +{ + Vector estimated_error_per_cell (triangulation.n_active_cells()); + + KellyErrorEstimator::FunctionMap neumann_boundary; + KellyErrorEstimator::estimate (mg_dof_handler, + QGauss3(), + neumann_boundary, + solution, + estimated_error_per_cell); + + triangulation.refine_and_coarsen_fixed_number (estimated_error_per_cell, + 0.3, 0.03); + triangulation.execute_coarsening_and_refinement (); }; template -void PoissonProblem::solve () { - Assert ((tria!=0) && (dof!=0), ExcNoTriaSelected()); +void LaplaceProblem::output_results (const unsigned int cycle) const +{ + string filename = "grid-"; + filename += ('0' + cycle); + Assert (cycle < 10, ExcInternalError()); - SolverControl control(4000, 1e-16); - PrimitiveVectorMemory > memory; - SolverCG,Vector > cg(control,memory); - - // solve - cg.solve (system_matrix, solution, right_hand_side, - PreconditionIdentity()); - // distribute solution - constraints.distribute (solution); + filename += ".eps"; + ofstream output (filename.c_str()); + + GridOut grid_out; + grid_out.write_eps (triangulation, output); }; template -int PoissonProblem::run (const unsigned int level) { - create_new (); - - cout << "Refinement level = " << level - << ", using elements of type <"; - switch (order) +void LaplaceProblem::run () +{ + for (unsigned int cycle=0; cycle<8; ++cycle) { - case 0: - cout << "criss-cross"; - break; - default: - cout << "Lagrange-" << order; - break; - }; - cout << ">" << endl; - - cout << " Making grid... "; - GridGenerator::hyper_cube (*tria); - tria->refine_global (level+1); - tria->begin_active()->set_refine_flag(); - (++(++(tria->begin_active())))->set_refine_flag(); - tria->execute_coarsening_and_refinement (); - cout << tria->n_active_cells() << " active cells." << endl; - - rhs = new RHSPoly(); - boundary_values = new Solution (); - + cout << "Cycle " << cycle << ':' << endl; - FiniteElement *fe; - PoissonEquation equation (*rhs); - Quadrature *quadrature; - Quadrature *boundary_quadrature; - switch (order) { - case 0: - fe = new FECrissCross(); - quadrature = new QCrissCross1(); - boundary_quadrature = new QGauss2(); - break; - case 1: - fe = new FEQ1(); - quadrature = new QGauss3(); - boundary_quadrature = new QGauss2(); - break; - case 2: - fe = new FEQ2(); - quadrature = new QGauss4(); - boundary_quadrature = new QGauss3(); - break; - case 3: - fe = new FEQ3(); - quadrature = new QGauss5(); - boundary_quadrature = new QGauss4(); - break; - case 4: - fe = new FEQ4(); - quadrature = new QGauss6(); - boundary_quadrature = new QGauss5(); - break; - default: - return 100000; - }; - - cout << " Distributing dofs... "; - dof->distribute_dofs (*fe); - cout << dof->DoFHandler::n_dofs() << " degrees of freedom." << endl; - n_dofs.push_back (dof->DoFHandler::n_dofs()); - - cout << " Assembling matrices..." << endl; - UpdateFlags update_flags = UpdateFlags(update_values | update_q_points | - update_gradients | update_JxW_values); - - map*> dirichlet_bc; - dirichlet_bc[0] = boundary_values; - assemble (equation, *quadrature, update_flags, dirichlet_bc); + if (cycle == 0) + { + GridGenerator::hyper_cube (triangulation); + triangulation.refine_global (1); + } + else + { + refine_grid (); + }; + - cout << " Solving..." << endl; - solve (); + cout << " Number of active cells: " + << triangulation.n_active_cells() + << endl; - Solution sol; - Vector l1_error_per_cell, l2_error_per_cell, linfty_error_per_cell; - Vector h1_seminorm_error_per_cell, h1_error_per_cell; - - cout << " Calculating L1 error... "; - VectorTools::integrate_difference (*dof, - solution, sol, - l1_error_per_cell, - *quadrature, L1_norm); - cout << l1_error_per_cell.l1_norm() << endl; - l1_error.push_back (l1_error_per_cell.l1_norm()); - - cout << " Calculating L2 error... "; - VectorTools::integrate_difference (*dof, - solution, sol, - l2_error_per_cell, - *quadrature, L2_norm); - cout << l2_error_per_cell.l2_norm() << endl; - l2_error.push_back (l2_error_per_cell.l2_norm()); - - cout << " Calculating L-infinity error... "; - VectorTools::integrate_difference (*dof, - solution, sol, - linfty_error_per_cell, - *quadrature, Linfty_norm); - cout << linfty_error_per_cell.linfty_norm() << endl; - linfty_error.push_back (linfty_error_per_cell.linfty_norm()); - - cout << " Calculating H1-seminorm error... "; - VectorTools::integrate_difference (*dof, - solution, sol, - h1_seminorm_error_per_cell, - *quadrature, H1_seminorm); - cout << h1_seminorm_error_per_cell.l2_norm() << endl; - h1_seminorm_error.push_back (h1_seminorm_error_per_cell.l2_norm()); - - cout << " Calculating H1 error... "; - VectorTools::integrate_difference (*dof, - solution, sol, - h1_error_per_cell, - *quadrature, H1_norm); - cout << h1_error_per_cell.l2_norm() << endl; - h1_error.push_back (h1_error_per_cell.l2_norm()); - - if (dof->DoFHandler::n_dofs()<=5000) - { - Vector l1_error_per_dof (dof->DoFHandler::n_dofs()); - Vector l2_error_per_dof (dof->DoFHandler::n_dofs()); - Vector linfty_error_per_dof (dof->DoFHandler::n_dofs()); - Vector h1_seminorm_error_per_dof (dof->DoFHandler::n_dofs()); - Vector h1_error_per_dof (dof->DoFHandler::n_dofs()); - DoFTools::distribute_cell_to_dof_vector (*dof, l1_error_per_cell, l1_error_per_dof); - DoFTools::distribute_cell_to_dof_vector (*dof, l2_error_per_cell, l2_error_per_dof); - DoFTools::distribute_cell_to_dof_vector (*dof, linfty_error_per_cell, - linfty_error_per_dof); - DoFTools::distribute_cell_to_dof_vector (*dof, h1_seminorm_error_per_cell, - h1_seminorm_error_per_dof); - DoFTools::distribute_cell_to_dof_vector (*dof, h1_error_per_cell, h1_error_per_dof); - -// Vector projected_solution; -// ConstraintMatrix constraints; -// constraints.close (); -// VectorTools::project (*dof, constraints, *fe, -// StraightBoundary(), *quadrature, -// sol, projected_solution, false, -// *boundary_quadrature); -// cout << " Calculating L2 error of projected solution... "; -// VectorTools::integrate_difference (*dof, -// projected_solution, sol, -// l2_error_per_cell, -// *quadrature, *fe, L2_norm); -// cout << l2_error_per_cell.l2_norm() << endl; - - - string filename; - filename = ('0'+order); - filename += "."; - filename += ('0'+level); - filename += ".ucd"; - cout << " Writing error plots to <" << filename << ">..." << endl; + setup_system (); + + cout << " Number of degrees of freedom: " + << mg_dof_handler.DoFHandler::n_dofs() + << endl; - DataOut out; - ofstream o(filename.c_str()); - out.attach_dof_handler (*dof); - - out.add_data_vector (solution, "u"); - out.add_data_vector (l1_error_per_dof, "L1_Error"); - out.add_data_vector (l2_error_per_dof, "L2_Error"); - out.add_data_vector (linfty_error_per_dof, "Linfty_Error"); - out.add_data_vector (h1_seminorm_error_per_dof, "H1_seminorm_Error"); - out.add_data_vector (h1_error_per_dof, "H1_Error"); - out.write_ucd (o); - o.close (); - } - else - cout << " Not writing error as grid." << endl; - - cout << endl; - - const unsigned int n_dofs = dof->DoFHandler::n_dofs(); - dof->clear (); - tria->set_boundary (0); - delete fe; - delete quadrature; - delete boundary_quadrature; - - return n_dofs; -}; + assemble_system (); + solve (); + output_results (cycle); + DataOut::EpsFlags eps_flags; + eps_flags.z_scaling = 4; + + DataOut data_out; + data_out.set_flags (eps_flags); -template -void PoissonProblem::print_history (string filename) const { - ofstream out(filename.c_str()); - out << "# n_dofs l1_error l2_error linfty_error h1_seminorm_error h1_error" - << endl; - for (unsigned int i=0; ih/2:" << endl; - cout << " L1 error : " << 1./average_l1 << endl - << " L2 error : " << 1./average_l2 << endl - << " Linfty error : " << 1./average_linfty << endl - << " H1 seminorm error: " << 1./average_h1_semi << endl - << " H1 error : " << 1./average_h1 << endl; - cout << "==========================================================\n"; - cout << "==========================================================\n"; }; + +int main () +{ + try + { + deallog.depth_console (0); - -int main () { - for (unsigned int order=0; order<5; ++order) + LaplaceProblem<2> laplace_problem_2d; + laplace_problem_2d.run (); + } + catch (exception &exc) { - PoissonProblem<2> problem (order); - - unsigned int level=0; - unsigned int n_dofs; - do - n_dofs = problem.run (level++); - while (n_dofs<25000); - - string filename; - switch (order) - { - case 0: - filename = "criss_cross"; - break; - case 1: - filename = "linear"; - break; - case 2: - filename = "quadratic"; - break; - case 3: - filename = "cubic"; - break; - case 4: - filename = "quartic"; - break; - }; - filename += ".history"; - - cout << endl << "Printing convergence history to <" - << filename << ">..." << endl; - problem.print_history (filename); - cout << endl << endl << endl; + cerr << endl << endl + << "----------------------------------------------------" + << endl; + cerr << "Exception on processing: " << endl + << exc.what() << endl + << "Aborting!" << endl + << "----------------------------------------------------" + << endl; + return 1; + } + catch (...) + { + cerr << endl << endl + << "----------------------------------------------------" + << endl; + cerr << "Unknown exception!" << endl + << "Aborting!" << endl + << "----------------------------------------------------" + << endl; + return 1; }; - + return 0; }; diff --git a/tests/deal.II/helmholtz.h b/tests/deal.II/helmholtz.h index f51156022a..5088d78f94 100644 --- a/tests/deal.II/helmholtz.h +++ b/tests/deal.II/helmholtz.h @@ -1,6 +1,6 @@ // $Id$ -#include +#include /** * Generator for system matrix and right hand side. @@ -16,7 +16,7 @@ class Helmholtz const Function& rhs); template - void build_mgmatrix(MGMatrix& A, + void build_mgmatrix(MGLevelObject& A, const MGDoFHandler& dof, const Quadrature& quadrature); }; diff --git a/tests/deal.II/helmholtz1mg.th b/tests/deal.II/helmholtz1mg.th index 5bf4392710..4487979877 100644 --- a/tests/deal.II/helmholtz1mg.th +++ b/tests/deal.II/helmholtz1mg.th @@ -2,7 +2,7 @@ template void -Helmholtz::build_mgmatrix(MGMatrix& A, +Helmholtz::build_mgmatrix(MGLevelObject& A, const MGDoFHandler& dof, const Quadrature& quadrature) { diff --git a/tests/deal.II/mg.cc b/tests/deal.II/mg.cc index 47698cab7b..bace5c9231 100644 --- a/tests/deal.II/mg.cc +++ b/tests/deal.II/mg.cc @@ -15,21 +15,27 @@ #include #include #include -#include #include #include -#include #include #include #include -#include -#include -#include #include -#include + +#include +#include +#include +#include +#include +#include #include "helmholtz.h" +#include + + + + template class RHSFunction : public Function { @@ -42,9 +48,9 @@ class RHSFunction : public Function class MGSmootherLAC : public MGSmootherBase { private: - SmartPointer > >matrices; + SmartPointer > >matrices; public: - MGSmootherLAC(MGMatrix >&); + MGSmootherLAC(MGLevelObject >&); virtual void smooth (const unsigned int level, Vector &u, @@ -115,8 +121,8 @@ int main() solver.solve(A,u,f,precondition); u = 0.; - vector mgstruct(tr.n_levels()); - MGMatrix > mgA(0,tr.n_levels()-1); + MGLevelObject mgstruct(0, tr.n_levels()-1); + MGLevelObject > mgA(0,tr.n_levels()-1); for (unsigned int i=0;i multigrid(mgdof, hanging_nodes, mgA, transfer); - PreconditionMG > + Multigrid<2> multigrid(mgdof, hanging_nodes, mgstruct, mgA, transfer); + PreconditionMG > mgprecondition(multigrid, smoother, smoother, coarse); solver.solve(A, u, f, mgprecondition); @@ -158,7 +164,7 @@ RHSFunction::value (const Point&, return 1.; } -MGSmootherLAC::MGSmootherLAC(MGMatrix >& matrix) +MGSmootherLAC::MGSmootherLAC(MGLevelObject >& matrix) : matrices(&matrix) {} diff --git a/tests/deal.II/mglocal.cc b/tests/deal.II/mglocal.cc index 5dec88651c..650f5a790b 100644 --- a/tests/deal.II/mglocal.cc +++ b/tests/deal.II/mglocal.cc @@ -15,19 +15,18 @@ #include #include #include -#include +#include #include #include -#include +#include #include #include #include #include -#include -#include -#include +#include +#include #include -#include +#include #include @@ -130,8 +129,8 @@ int main() SolverControl control(20, 1.e-12, true); SolverRichardson<> solver(control, mem); - vector mgstruct(tr.n_levels()); - MGMatrix > mgA(0,tr.n_levels()-1); + MGLevelObject mgstruct(0, tr.n_levels()-1); + MGLevelObject > mgA(0,tr.n_levels()-1); for (unsigned int i=0;i multigrid(mgdof, hanging_nodes, mgA, transfer, tr.n_levels()-2); - PreconditionMG > + Multigrid<2> multigrid(mgdof, hanging_nodes, mgstruct, mgA, transfer, tr.n_levels()-2); + PreconditionMG > mgprecondition(multigrid, smoother, smoother, coarse); u = 0.; diff --git a/tests/lac/Makefile b/tests/lac/Makefile index c8c0a58a0f..c4fa2492a9 100644 --- a/tests/lac/Makefile +++ b/tests/lac/Makefile @@ -24,9 +24,11 @@ debug-mode = on ############################################################ libs.g = $(lib-lac.g) \ - $(lib-base.g) + $(lib-base.g) \ + $(lib-deal2-2d.g) libs = $(lib-lac.o) \ - $(lib-base.o) + $(lib-base.o) \ + $(lib-deal2-2d.o) ############################################################ diff --git a/tests/lac/mg.cc b/tests/lac/mg.cc index 0279dfb26d..a5f39f0538 100644 --- a/tests/lac/mg.cc +++ b/tests/lac/mg.cc @@ -15,7 +15,9 @@ #include #include #include -#include +#include +#include + #define TYPE double #define ACCURACY 1.e-5 @@ -44,9 +46,9 @@ class FDMG /** * Pointer to the level matrices. */ - SmartPointer > >matrices; + SmartPointer > >matrices; public: - FDMG(unsigned int maxlevel, MGMatrix >& matrices, + FDMG(unsigned int maxlevel, MGPerLevelObject >& matrices, FDMGTransfer& transfer); @@ -64,9 +66,9 @@ class MGSmootherLAC public MGSmootherBase { private: - SmartPointer > >matrices; + SmartPointer > >matrices; public: - MGSmootherLAC(MGMatrix >&); + MGSmootherLAC(MGPerLevelObject >&); virtual void smooth (const unsigned int level, Vector &u, @@ -113,7 +115,7 @@ int main() // Make matrix vector structure(maxlevel+1); - MGMatrix > A(minlevel,maxlevel); + MGPerLevelObject > A(minlevel,maxlevel); FDMatrix testproblem(size, size); @@ -157,7 +159,7 @@ int main() } } -FDMG::FDMG(unsigned int maxlevel, MGMatrix >& matrices, +FDMG::FDMG(unsigned int maxlevel, MGPerLevelObject >& matrices, FDMGTransfer& transfer) : MGBase(transfer, 0, maxlevel), @@ -193,7 +195,7 @@ FDMG::copy_from_mg(Vector& v) } -MGSmootherLAC::MGSmootherLAC(MGMatrix >& matrix) +MGSmootherLAC::MGSmootherLAC(MGPerLevelObject >& matrix) : matrices(&matrix) {} diff --git a/tests/lac/mgbase.cc b/tests/lac/mgbase.cc index 984437ea8f..882d62039e 100644 --- a/tests/lac/mgbase.cc +++ b/tests/lac/mgbase.cc @@ -3,7 +3,8 @@ // deal_II_libraries.g=-ldeal_II_2d.g // deal_II_libraries=-ldeal_II_2d -#include +#include +#include #include #include diff --git a/tests/lac/testmatrix.h b/tests/lac/testmatrix.h index a81ca7210b..5c1b392b77 100644 --- a/tests/lac/testmatrix.h +++ b/tests/lac/testmatrix.h @@ -2,7 +2,7 @@ #include #include -#include +#include /** * Finite difference matrix on uniform grid.