From e7a7f11e886d9010a1c3e76fbefb77d7b9964980 Mon Sep 17 00:00:00 2001 From: wolf Date: Fri, 19 Nov 1999 08:12:28 +0000 Subject: [PATCH] Large clean-ups of multigrid things (part 1). git-svn-id: https://svn.dealii.org/trunk@1893 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/multigrid/mg_base.h | 397 +++++++++++------- .../include/multigrid/multigrid.templates.h | 38 +- .../include/numerics/multigrid.templates.h | 38 +- deal.II/deal.II/source/multigrid/mg_base.cc | 103 +++-- deal.II/lac/include/lac/mgbase.h | 397 +++++++++++------- deal.II/lac/source/mgbase.cc | 103 +++-- 6 files changed, 664 insertions(+), 412 deletions(-) diff --git a/deal.II/deal.II/include/multigrid/mg_base.h b/deal.II/deal.II/include/multigrid/mg_base.h index 1a5d713443..4214e4269b 100644 --- a/deal.II/deal.II/include/multigrid/mg_base.h +++ b/deal.II/deal.II/include/multigrid/mg_base.h @@ -14,16 +14,22 @@ #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. + * 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. + * Virtual destructor. Does + * nothing in particular, but + * needs to be declared anyhow. */ virtual ~MGCoarseGridSolver(); @@ -57,24 +63,33 @@ class MGSmootherBase : public Subscriptor { public: /** - * Virtual destructor. + * Virtual destructor. Does + * nothing in particular, but + * needs to be declared anyhow. */ virtual ~MGSmootherBase(); /** - * Smooth the residual of #u# on the given - * level. 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. + * 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; - }; @@ -93,7 +108,8 @@ class MGSmootherIdentity : public MGSmootherBase { public: /** - * Implementation of the interface in #MGSmootherBase#. + * Implementation of the + * interface in #MGSmootherBase#. * This function does nothing. */ virtual void smooth (const unsigned int level, @@ -103,6 +119,7 @@ class MGSmootherIdentity : public MGSmootherBase + /** * Base class used to declare the operations needed by a concrete class * implementing prolongation and restriction of vectors in the multigrid @@ -132,7 +149,7 @@ class MGTransferBase : public Subscriptor * are degrees of freedom on the finer * level. */ - virtual void prolongate (const unsigned int to_level, + virtual void prolongate (const unsigned int to_level, Vector &dst, const Vector &src) const = 0; @@ -149,7 +166,7 @@ class MGTransferBase : public Subscriptor * are degrees of freedom on the coarser * level. */ - virtual void restrict (const unsigned int from_level, + virtual void restrict (const unsigned int from_level, Vector &dst, const Vector &src) const = 0; }; @@ -157,37 +174,61 @@ class MGTransferBase : public Subscriptor /** - * An array with a vector for each level. - * The purpose of this class is mostly to provide range checking that is missing in - * those ultra-fashionable STL vectors. + * 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 vector +template > +class MGVector : public Subscriptor { public: /** - * Constructor allowing to initialize the number of levels. + * Constructor allowing to + * initialize the number of + * levels. */ - MGVector(unsigned int minlevel, unsigned int maxlevel); + MGVector (const unsigned int minlevel, + const unsigned int maxlevel); + /** - * Safe access operator. + * Access vector on level #level#. */ - VECTOR& operator[](unsigned int); + VECTOR & operator[] (const unsigned int level); + /** - * Safe access operator. + * Access vector on level + * #level#. Constant version. */ - const VECTOR& operator[](unsigned int) const; + const VECTOR & operator[] (const unsigned int level) const; + /** - * Delete all entries. + * 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. */ - unsigned int minlevel; + const unsigned int minlevel; + + /** + * Array of the vectors to be held. + */ + vector vectors; }; @@ -234,6 +275,10 @@ class MGMatrix : public Subscriptor, * 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 @@ -246,7 +291,9 @@ class MGCoarseGridLACIteration : public MGCoarseGridSolver * preconditioning method for later * use. */ - MGCoarseGridLACIteration(SOLVER&, const MATRIX&, const PRECOND&); + MGCoarseGridLACIteration (SOLVER &, + const MATRIX &, + const PRECOND &); /** * Implementation of the abstract @@ -267,7 +314,7 @@ class MGCoarseGridLACIteration : public MGCoarseGridSolver /** * Reference to the matrix. */ - const MATRIX& matrix; + const SmartPointer matrix; /** * Reference to the preconditioner. @@ -277,16 +324,16 @@ class MGCoarseGridLACIteration : public MGCoarseGridSolver + /** - * Basic class for multigrid preconditioning. + * Basic class for preconditioning by multigrid. * * The functionality of the multigrid method is restricted to defect - * correction. It is not iterative and the start solution is - * always zero. Since by this uEl and - * uAl (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 + * 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 @@ -295,65 +342,109 @@ class MGCoarseGridLACIteration : public MGCoarseGridSolver * 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&); - protected: + public: /** - * Highest level of cells. + * Constructor, subject to change. */ - unsigned int maxlevel; + 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); /** - * Level for coarse grid solution. + * Exception. */ - unsigned int minlevel; + DeclException2(ExcSwitchedLevels, int, int, + << "minlevel and maxlevel switched, should be: " + << arg1 << "<=" << arg2); + + protected: /** - * Auxiliary vector, defect. + * Highest level of cells. */ - MGVector > d; + unsigned int maxlevel; /** - * Auxiliary vector, solution. + * Level for coarse grid solution. */ - MGVector > s; + unsigned int minlevel; + /** - * Prolongation and restriction object. + * Auxiliary vector. */ - SmartPointer transfer; + MGVector > defect; - private: /** * Auxiliary vector. */ - Vector t; + MGVector > solution; /** - * Exception. + * Prolongation and restriction object. */ - DeclException2(ExcSwitchedLevels, int, int, - << "minlevel and maxlevel switched, should be: " - << arg1 << "<=" << arg2); - - protected: + SmartPointer transfer; /** * The actual v-cycle multigrid method. - * This function is called on the - * highest level and recursively - * invokes itself down to the - * coarsest. There, it calls - * #coarse_grid_solver# and - * proceeds back up. - */ - void level_mgstep(unsigned int level, - const MGSmootherBase& pre_smooth, - const MGSmootherBase& post_smooth, - const MGCoarseGridSolver& cgs); + * #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 @@ -361,39 +452,16 @@ class MGBase : public Subscriptor * This is implemented in a * derived class. */ - virtual void level_vmult(unsigned int level, - Vector& dst, - const Vector& src, - const Vector& rhs) = 0; - + virtual void level_vmult (const unsigned int level, + Vector &dst, + const Vector &src, + const Vector &rhs) = 0; - public: - /** - * Constructor, subject to change. - */ - MGBase(const MGTransferBase& transfer, - unsigned int minlevel, unsigned int maxlevel); - /** - * Virtual destructor. - */ - virtual ~MGBase(); + private: /** - * Execute 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#. + * Auxiliary vector. */ - void vcycle(const MGSmootherBase& pre_smooth, - const MGSmootherBase& post_smooth, - const MGCoarseGridSolver& cgs); + Vector t; }; @@ -407,9 +475,10 @@ class MGBase : public Subscriptor * 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 +template > class PreconditionMG { public: @@ -419,29 +488,39 @@ class PreconditionMG * pre-smoother, post-smoother and * coarse grid solver. */ - PreconditionMG(MG& mg, - const MGSmootherBase& pre, - const MGSmootherBase& post, - const MGCoarseGridSolver& coarse); + PreconditionMG(MG &mg, + const MGSmootherBase &pre, + const MGSmootherBase &post, + const MGCoarseGridSolver &coarse); + /** - * Preconditioning operator. - * This is the operator used by LAC - * iterative solvers. + * 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; + void operator() (VECTOR &dst, + const VECTOR &src) const; + private: /** - * The mulrigrid object. + * The multigrid object. */ SmartPointer multigrid; + /** * The pre-smoothing object. */ SmartPointer pre; + /** * The post-smoothing object. */ SmartPointer post; + /** * The coarse grid solver. */ @@ -451,58 +530,44 @@ class PreconditionMG - -template -MGCoarseGridLACIteration:: -MGCoarseGridLACIteration(SOLVER& s, const MATRIX& m, const PRECOND& p) - : - solver(s), matrix(m), precondition(p) -{} - -template -void -MGCoarseGridLACIteration::operator() -(unsigned int, Vector& dst, const Vector& src) const -{ - solver.solve(matrix, dst, src, precondition); -} - - - - /* ------------------------------------------------------------------- */ template -MGVector::MGVector(unsigned int min, unsigned int max) +MGVector::MGVector(const unsigned int min, + const unsigned int max) : - vector(max-min+1), - minlevel(min) -{} + minlevel(min), + vectors(max-min+1) +{}; + template -VECTOR& -MGVector::operator[](unsigned int i) +VECTOR & +MGVector::operator[] (const unsigned int i) { - Assert((i>=minlevel)&&(i::operator[](i-minlevel); + Assert((i>=minlevel) && (i -const VECTOR& -MGVector::operator[](unsigned int i) const +const VECTOR & +MGVector::operator[] (const unsigned int i) const { - Assert((i>=minlevel)&&(i::operator[](i-minlevel); + Assert((i>=minlevel) && (i void MGVector::clear() { typename vector::iterator v; - for (v = begin(); v != end(); ++v) + for (v = vectors.begin(); v != vectors.end(); ++v) v->clear(); } @@ -525,7 +590,8 @@ MGMatrix::operator[](unsigned int i) { Assert((i>=minlevel)&&(i::operator[](i-minlevel); -} +}; + template @@ -534,18 +600,49 @@ 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) +PreconditionMG::PreconditionMG(MG &mg, + const MGSmootherBase &pre, + const MGSmootherBase &post, + const MGCoarseGridSolver &coarse) : multigrid(&mg), pre(&pre), @@ -557,8 +654,8 @@ PreconditionMG::PreconditionMG(MG& mg, template void -PreconditionMG::operator() (VECTOR& dst, - const VECTOR& src) const +PreconditionMG::operator() (VECTOR &dst, + const VECTOR &src) const { multigrid->copy_to_mg(src); multigrid->vcycle(*pre, *post, *coarse); diff --git a/deal.II/deal.II/include/multigrid/multigrid.templates.h b/deal.II/deal.II/include/multigrid/multigrid.templates.h index 7caabae233..5d01b9b16a 100644 --- a/deal.II/deal.II/include/multigrid/multigrid.templates.h +++ b/deal.II/deal.II/include/multigrid/multigrid.templates.h @@ -1,12 +1,11 @@ // $Id$ #include +#include #include -inline int minimum(int a, int b) -{ - return ((a MG::MG(const MGDoFHandler& dofs, @@ -16,17 +15,19 @@ MG::MG(const MGDoFHandler& dofs, unsigned int minl, unsigned int maxl) : MGBase(transfer, minl, - minimum(dofs.get_tria().n_levels()-1, maxl)), + min(dofs.get_tria().n_levels()-1, maxl)), dofs(&dofs), matrices(&matrices), hanging_nodes(hanging_nodes) { for (unsigned int level = minlevel; level<=maxlevel ; ++level) { - s[level].reinit(matrices[level].m()); - d[level].reinit(matrices[level].m()); - } -} + solution[level].reinit(matrices[level].m()); + defect[level].reinit(matrices[level].m()); + }; +}; + + template template @@ -38,7 +39,7 @@ MG::copy_to_mg(const Vector& src) // set the elements of the vectors // on all levels to zero - d.clear(); + defect.clear(); // hanging_nodes.condense(src); vector index(fe_dofs); @@ -65,7 +66,7 @@ MG::copy_to_mg(const Vector& src) dc->get_dof_indices(index); c->get_mg_dof_indices(mgindex); for (unsigned int i=0;ilevel()](mgindex[i]) = src(index[i]); + defect[c->level()](mgindex[i]) = src(index[i]); // Delete values on refinement edge for (unsigned int face_n = 0; face_n< GeometryInfo::faces_per_cell; ++face_n) @@ -75,14 +76,14 @@ MG::copy_to_mg(const Vector& src) { face->get_mg_dof_indices(mgfaceindex); for (unsigned int i=0;ilevel()](mgfaceindex[i]) = 0.; + defect[c->level()](mgfaceindex[i]) = 0.; } } } // if (level > (int) minlevel) // transfer->restrict(level, d[level-1], d[level]); if (level < (int) maxlevel) - transfer->restrict(level+1, d[level], d[level+1]); + transfer->restrict(level+1, defect[level], defect[level+1]); } // for (unsigned int i=minlevel; i<= maxlevel; ++i) @@ -94,7 +95,8 @@ MG::copy_to_mg(const Vector& src) // } } -template template +template +template void MG::copy_from_mg(Vector& dst) const { @@ -112,13 +114,13 @@ MG::copy_from_mg(Vector& dst) const // } DoFHandler::active_cell_iterator dc = dofs->DoFHandler::begin_active(); - for (MGDoFHandler::active_cell_iterator c = dofs->begin_active() - ; c != dofs->end() ; ++c, ++dc) + for (MGDoFHandler::active_cell_iterator c = dofs->begin_active(); + c != dofs->end(); ++c, ++dc) { dc->get_dof_indices(index); c->get_mg_dof_indices(mgindex); for (unsigned int i=0;ilevel()](mgindex[i]); + dst(index[i]) = solution[c->level()](mgindex[i]); } hanging_nodes.set_zero(dst); // { @@ -130,6 +132,8 @@ MG::copy_from_mg(Vector& dst) const // } } + + template void MG::level_vmult(unsigned int l, diff --git a/deal.II/deal.II/include/numerics/multigrid.templates.h b/deal.II/deal.II/include/numerics/multigrid.templates.h index 7caabae233..5d01b9b16a 100644 --- a/deal.II/deal.II/include/numerics/multigrid.templates.h +++ b/deal.II/deal.II/include/numerics/multigrid.templates.h @@ -1,12 +1,11 @@ // $Id$ #include +#include #include -inline int minimum(int a, int b) -{ - return ((a MG::MG(const MGDoFHandler& dofs, @@ -16,17 +15,19 @@ MG::MG(const MGDoFHandler& dofs, unsigned int minl, unsigned int maxl) : MGBase(transfer, minl, - minimum(dofs.get_tria().n_levels()-1, maxl)), + min(dofs.get_tria().n_levels()-1, maxl)), dofs(&dofs), matrices(&matrices), hanging_nodes(hanging_nodes) { for (unsigned int level = minlevel; level<=maxlevel ; ++level) { - s[level].reinit(matrices[level].m()); - d[level].reinit(matrices[level].m()); - } -} + solution[level].reinit(matrices[level].m()); + defect[level].reinit(matrices[level].m()); + }; +}; + + template template @@ -38,7 +39,7 @@ MG::copy_to_mg(const Vector& src) // set the elements of the vectors // on all levels to zero - d.clear(); + defect.clear(); // hanging_nodes.condense(src); vector index(fe_dofs); @@ -65,7 +66,7 @@ MG::copy_to_mg(const Vector& src) dc->get_dof_indices(index); c->get_mg_dof_indices(mgindex); for (unsigned int i=0;ilevel()](mgindex[i]) = src(index[i]); + defect[c->level()](mgindex[i]) = src(index[i]); // Delete values on refinement edge for (unsigned int face_n = 0; face_n< GeometryInfo::faces_per_cell; ++face_n) @@ -75,14 +76,14 @@ MG::copy_to_mg(const Vector& src) { face->get_mg_dof_indices(mgfaceindex); for (unsigned int i=0;ilevel()](mgfaceindex[i]) = 0.; + defect[c->level()](mgfaceindex[i]) = 0.; } } } // if (level > (int) minlevel) // transfer->restrict(level, d[level-1], d[level]); if (level < (int) maxlevel) - transfer->restrict(level+1, d[level], d[level+1]); + transfer->restrict(level+1, defect[level], defect[level+1]); } // for (unsigned int i=minlevel; i<= maxlevel; ++i) @@ -94,7 +95,8 @@ MG::copy_to_mg(const Vector& src) // } } -template template +template +template void MG::copy_from_mg(Vector& dst) const { @@ -112,13 +114,13 @@ MG::copy_from_mg(Vector& dst) const // } DoFHandler::active_cell_iterator dc = dofs->DoFHandler::begin_active(); - for (MGDoFHandler::active_cell_iterator c = dofs->begin_active() - ; c != dofs->end() ; ++c, ++dc) + for (MGDoFHandler::active_cell_iterator c = dofs->begin_active(); + c != dofs->end(); ++c, ++dc) { dc->get_dof_indices(index); c->get_mg_dof_indices(mgindex); for (unsigned int i=0;ilevel()](mgindex[i]); + dst(index[i]) = solution[c->level()](mgindex[i]); } hanging_nodes.set_zero(dst); // { @@ -130,6 +132,8 @@ MG::copy_from_mg(Vector& dst) const // } } + + template void MG::level_vmult(unsigned int l, diff --git a/deal.II/deal.II/source/multigrid/mg_base.cc b/deal.II/deal.II/source/multigrid/mg_base.cc index 9df47ff594..3662347101 100644 --- a/deal.II/deal.II/source/multigrid/mg_base.cc +++ b/deal.II/deal.II/source/multigrid/mg_base.cc @@ -4,104 +4,129 @@ #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, - unsigned minlevel, unsigned maxlevel) + +////////////////////////////////////////////////////////////////////// + + + +MGBase::MGBase(const MGTransferBase &transfer, + const unsigned minlevel, + const unsigned maxlevel) : - maxlevel(maxlevel), minlevel(minlevel), - d(minlevel,maxlevel), - s(minlevel,maxlevel), + maxlevel(maxlevel), + minlevel(minlevel), + defect(minlevel,maxlevel), + solution(minlevel,maxlevel), transfer(&transfer) { - Assert(minlevel <= maxlevel, ExcSwitchedLevels(minlevel, maxlevel)); -} + Assert(minlevel <= maxlevel, + ExcSwitchedLevels(minlevel, maxlevel)); +}; + + void -MGBase::vcycle(const MGSmootherBase& pre_smooth, - const MGSmootherBase& post_smooth, - const MGCoarseGridSolver& coarse_grid_solver) +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(unsigned level, - const MGSmootherBase& pre_smooth, - const MGSmootherBase& post_smooth, - const MGCoarseGridSolver& coarse_grid_solver) +MGBase::level_mgstep(const unsigned int level, + const MGSmootherBase &pre_smooth, + const MGSmootherBase &post_smooth, + const MGCoarseGridSolver &coarse_grid_solver) { - s[level] = 0.; + solution[level] = 0.; if (level == minlevel) { - coarse_grid_solver(level, s[level], d[level]); + coarse_grid_solver(level, solution[level], defect[level]); return; } // smoothing of the residual by modifying s - pre_smooth.smooth(level, s[level], d[level]); + pre_smooth.smooth(level, solution[level], defect[level]); // t = d-As - t.reinit(s[level].size()); - level_vmult(level, t, s[level], d[level]); + t.reinit(solution[level].size()); + level_vmult(level, t, solution[level], defect[level]); // make t rhs of lower level - transfer->restrict(level, d[level-1], t); + transfer->restrict(level, defect[level-1], t); // do recursion level_mgstep(level-1, pre_smooth, post_smooth, coarse_grid_solver); // do coarse grid correction - t.reinit(s[level].size()); - transfer->prolongate(level, t, s[level-1]); - s[level].add(t); + t.reinit(solution[level].size()); + transfer->prolongate(level, t, solution[level-1]); + solution[level].add(t); // smoothing (modify s again) - post_smooth.smooth(level, s[level], d[level]); + post_smooth.smooth(level, solution[level], defect[level]); } ////////////////////////////////////////////////////////////////////// MGCoarseGridSolver::~MGCoarseGridSolver() -{} +{}; + ////////////////////////////////////////////////////////////////////// MGTransferBase::~MGTransferBase() -{} - -MGSmootherBase::~MGSmootherBase() -{} - +{}; -////////////////////////////////////////////////////////////////////// -void -MGSmootherIdentity::smooth (const unsigned int, - Vector &, - const Vector &) const -{} diff --git a/deal.II/lac/include/lac/mgbase.h b/deal.II/lac/include/lac/mgbase.h index d1fd1fe0c3..9456174df4 100644 --- a/deal.II/lac/include/lac/mgbase.h +++ b/deal.II/lac/include/lac/mgbase.h @@ -14,16 +14,22 @@ #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. + * 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. + * Virtual destructor. Does + * nothing in particular, but + * needs to be declared anyhow. */ virtual ~MGCoarseGridSolver(); @@ -57,24 +63,33 @@ class MGSmootherBase : public Subscriptor { public: /** - * Virtual destructor. + * Virtual destructor. Does + * nothing in particular, but + * needs to be declared anyhow. */ virtual ~MGSmootherBase(); /** - * Smooth the residual of #u# on the given - * level. 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. + * 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; - }; @@ -93,7 +108,8 @@ class MGSmootherIdentity : public MGSmootherBase { public: /** - * Implementation of the interface in #MGSmootherBase#. + * Implementation of the + * interface in #MGSmootherBase#. * This function does nothing. */ virtual void smooth (const unsigned int level, @@ -103,6 +119,7 @@ class MGSmootherIdentity : public MGSmootherBase + /** * Base class used to declare the operations needed by a concrete class * implementing prolongation and restriction of vectors in the multigrid @@ -132,7 +149,7 @@ class MGTransferBase : public Subscriptor * are degrees of freedom on the finer * level. */ - virtual void prolongate (const unsigned int to_level, + virtual void prolongate (const unsigned int to_level, Vector &dst, const Vector &src) const = 0; @@ -149,7 +166,7 @@ class MGTransferBase : public Subscriptor * are degrees of freedom on the coarser * level. */ - virtual void restrict (const unsigned int from_level, + virtual void restrict (const unsigned int from_level, Vector &dst, const Vector &src) const = 0; }; @@ -157,37 +174,61 @@ class MGTransferBase : public Subscriptor /** - * An array with a vector for each level. - * The purpose of this class is mostly to provide range checking that is missing in - * those ultra-fashionable STL vectors. + * 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 vector +template > +class MGVector : public Subscriptor { public: /** - * Constructor allowing to initialize the number of levels. + * Constructor allowing to + * initialize the number of + * levels. */ - MGVector(unsigned int minlevel, unsigned int maxlevel); + MGVector (const unsigned int minlevel, + const unsigned int maxlevel); + /** - * Safe access operator. + * Access vector on level #level#. */ - VECTOR& operator[](unsigned int); + VECTOR & operator[] (const unsigned int level); + /** - * Safe access operator. + * Access vector on level + * #level#. Constant version. */ - const VECTOR& operator[](unsigned int) const; + const VECTOR & operator[] (const unsigned int level) const; + /** - * Delete all entries. + * 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. */ - unsigned int minlevel; + const unsigned int minlevel; + + /** + * Array of the vectors to be held. + */ + vector vectors; }; @@ -234,6 +275,10 @@ class MGMatrix : public Subscriptor, * 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 @@ -246,7 +291,9 @@ class MGCoarseGridLACIteration : public MGCoarseGridSolver * preconditioning method for later * use. */ - MGCoarseGridLACIteration(SOLVER&, const MATRIX&, const PRECOND&); + MGCoarseGridLACIteration (SOLVER &, + const MATRIX &, + const PRECOND &); /** * Implementation of the abstract @@ -267,7 +314,7 @@ class MGCoarseGridLACIteration : public MGCoarseGridSolver /** * Reference to the matrix. */ - const MATRIX& matrix; + const SmartPointer matrix; /** * Reference to the preconditioner. @@ -277,16 +324,16 @@ class MGCoarseGridLACIteration : public MGCoarseGridSolver + /** - * Basic class for multigrid preconditioning. + * Basic class for preconditioning by multigrid. * * The functionality of the multigrid method is restricted to defect - * correction. It is not iterative and the start solution is - * always zero. Since by this uEl and - * uAl (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 + * 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 @@ -295,65 +342,109 @@ class MGCoarseGridLACIteration : public MGCoarseGridSolver * 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&); - protected: + public: /** - * Highest level of cells. + * Constructor, subject to change. */ - unsigned int maxlevel; + 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); /** - * Level for coarse grid solution. + * Exception. */ - unsigned int minlevel; + DeclException2(ExcSwitchedLevels, int, int, + << "minlevel and maxlevel switched, should be: " + << arg1 << "<=" << arg2); + + protected: /** - * Auxiliary vector, defect. + * Highest level of cells. */ - MGVector > d; + unsigned int maxlevel; /** - * Auxiliary vector, solution. + * Level for coarse grid solution. */ - MGVector > s; + unsigned int minlevel; + /** - * Prolongation and restriction object. + * Auxiliary vector. */ - SmartPointer transfer; + MGVector > defect; - private: /** * Auxiliary vector. */ - Vector t; + MGVector > solution; /** - * Exception. + * Prolongation and restriction object. */ - DeclException2(ExcSwitchedLevels, int, int, - << "minlevel and maxlevel switched, should be: " - << arg1 << "<=" << arg2); - - protected: + SmartPointer transfer; /** * The actual v-cycle multigrid method. - * This function is called on the - * highest level and recursively - * invokes itself down to the - * coarsest. There, it calls - * #coarse_grid_solver# and - * proceeds back up. - */ - void level_mgstep(unsigned int level, - const MGSmootherBase& pre_smooth, - const MGSmootherBase& post_smooth, - const MGCoarseGridSolver& cgs); + * #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 @@ -361,39 +452,16 @@ class MGBase : public Subscriptor * This is implemented in a * derived class. */ - virtual void level_vmult(unsigned int level, - Vector& dst, - const Vector& src, - const Vector& rhs) = 0; - + virtual void level_vmult (const unsigned int level, + Vector &dst, + const Vector &src, + const Vector &rhs) = 0; - public: - /** - * Constructor, subject to change. - */ - MGBase(const MGTransferBase& transfer, - unsigned int minlevel, unsigned int maxlevel); - /** - * Virtual destructor. - */ - virtual ~MGBase(); + private: /** - * Execute 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#. + * Auxiliary vector. */ - void vcycle(const MGSmootherBase& pre_smooth, - const MGSmootherBase& post_smooth, - const MGCoarseGridSolver& cgs); + Vector t; }; @@ -407,9 +475,10 @@ class MGBase : public Subscriptor * 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 +template > class PreconditionMG { public: @@ -419,29 +488,39 @@ class PreconditionMG * pre-smoother, post-smoother and * coarse grid solver. */ - PreconditionMG(MG& mg, - const MGSmootherBase& pre, - const MGSmootherBase& post, - const MGCoarseGridSolver& coarse); + PreconditionMG(MG &mg, + const MGSmootherBase &pre, + const MGSmootherBase &post, + const MGCoarseGridSolver &coarse); + /** - * Preconditioning operator. - * This is the operator used by LAC - * iterative solvers. + * 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; + void operator() (VECTOR &dst, + const VECTOR &src) const; + private: /** - * The mulrigrid object. + * The multigrid object. */ SmartPointer multigrid; + /** * The pre-smoothing object. */ SmartPointer pre; + /** * The post-smoothing object. */ SmartPointer post; + /** * The coarse grid solver. */ @@ -451,58 +530,44 @@ class PreconditionMG - -template -MGCoarseGridLACIteration:: -MGCoarseGridLACIteration(SOLVER& s, const MATRIX& m, const PRECOND& p) - : - solver(s), matrix(m), precondition(p) -{} - -template -void -MGCoarseGridLACIteration::operator() -(unsigned int, Vector& dst, const Vector& src) const -{ - solver.solve(matrix, dst, src, precondition); -} - - - - /* ------------------------------------------------------------------- */ template -MGVector::MGVector(unsigned int min, unsigned int max) +MGVector::MGVector(const unsigned int min, + const unsigned int max) : - vector(max-min+1), - minlevel(min) -{} + minlevel(min), + vectors(max-min+1) +{}; + template -VECTOR& -MGVector::operator[](unsigned int i) +VECTOR & +MGVector::operator[] (const unsigned int i) { - Assert((i>=minlevel)&&(i::operator[](i-minlevel); + Assert((i>=minlevel) && (i -const VECTOR& -MGVector::operator[](unsigned int i) const +const VECTOR & +MGVector::operator[] (const unsigned int i) const { - Assert((i>=minlevel)&&(i::operator[](i-minlevel); + Assert((i>=minlevel) && (i void MGVector::clear() { typename vector::iterator v; - for (v = begin(); v != end(); ++v) + for (v = vectors.begin(); v != vectors.end(); ++v) v->clear(); } @@ -525,7 +590,8 @@ MGMatrix::operator[](unsigned int i) { Assert((i>=minlevel)&&(i::operator[](i-minlevel); -} +}; + template @@ -534,18 +600,49 @@ 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) +PreconditionMG::PreconditionMG(MG &mg, + const MGSmootherBase &pre, + const MGSmootherBase &post, + const MGCoarseGridSolver &coarse) : multigrid(&mg), pre(&pre), @@ -557,8 +654,8 @@ PreconditionMG::PreconditionMG(MG& mg, template void -PreconditionMG::operator() (VECTOR& dst, - const VECTOR& src) const +PreconditionMG::operator() (VECTOR &dst, + const VECTOR &src) const { multigrid->copy_to_mg(src); multigrid->vcycle(*pre, *post, *coarse); diff --git a/deal.II/lac/source/mgbase.cc b/deal.II/lac/source/mgbase.cc index 9df47ff594..3662347101 100644 --- a/deal.II/lac/source/mgbase.cc +++ b/deal.II/lac/source/mgbase.cc @@ -4,104 +4,129 @@ #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, - unsigned minlevel, unsigned maxlevel) + +////////////////////////////////////////////////////////////////////// + + + +MGBase::MGBase(const MGTransferBase &transfer, + const unsigned minlevel, + const unsigned maxlevel) : - maxlevel(maxlevel), minlevel(minlevel), - d(minlevel,maxlevel), - s(minlevel,maxlevel), + maxlevel(maxlevel), + minlevel(minlevel), + defect(minlevel,maxlevel), + solution(minlevel,maxlevel), transfer(&transfer) { - Assert(minlevel <= maxlevel, ExcSwitchedLevels(minlevel, maxlevel)); -} + Assert(minlevel <= maxlevel, + ExcSwitchedLevels(minlevel, maxlevel)); +}; + + void -MGBase::vcycle(const MGSmootherBase& pre_smooth, - const MGSmootherBase& post_smooth, - const MGCoarseGridSolver& coarse_grid_solver) +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(unsigned level, - const MGSmootherBase& pre_smooth, - const MGSmootherBase& post_smooth, - const MGCoarseGridSolver& coarse_grid_solver) +MGBase::level_mgstep(const unsigned int level, + const MGSmootherBase &pre_smooth, + const MGSmootherBase &post_smooth, + const MGCoarseGridSolver &coarse_grid_solver) { - s[level] = 0.; + solution[level] = 0.; if (level == minlevel) { - coarse_grid_solver(level, s[level], d[level]); + coarse_grid_solver(level, solution[level], defect[level]); return; } // smoothing of the residual by modifying s - pre_smooth.smooth(level, s[level], d[level]); + pre_smooth.smooth(level, solution[level], defect[level]); // t = d-As - t.reinit(s[level].size()); - level_vmult(level, t, s[level], d[level]); + t.reinit(solution[level].size()); + level_vmult(level, t, solution[level], defect[level]); // make t rhs of lower level - transfer->restrict(level, d[level-1], t); + transfer->restrict(level, defect[level-1], t); // do recursion level_mgstep(level-1, pre_smooth, post_smooth, coarse_grid_solver); // do coarse grid correction - t.reinit(s[level].size()); - transfer->prolongate(level, t, s[level-1]); - s[level].add(t); + t.reinit(solution[level].size()); + transfer->prolongate(level, t, solution[level-1]); + solution[level].add(t); // smoothing (modify s again) - post_smooth.smooth(level, s[level], d[level]); + post_smooth.smooth(level, solution[level], defect[level]); } ////////////////////////////////////////////////////////////////////// MGCoarseGridSolver::~MGCoarseGridSolver() -{} +{}; + ////////////////////////////////////////////////////////////////////// MGTransferBase::~MGTransferBase() -{} - -MGSmootherBase::~MGSmootherBase() -{} - +{}; -////////////////////////////////////////////////////////////////////// -void -MGSmootherIdentity::smooth (const unsigned int, - Vector &, - const Vector &) const -{} -- 2.39.5