#include <vector>
#include <base/smartpointer.h>
+
+
/**
- * 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();
{
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<double> &u,
const Vector<double> &rhs) const = 0;
-
};
{
public:
/**
- * Implementation of the interface in #MGSmootherBase#.
+ * Implementation of the
+ * interface in #MGSmootherBase#.
* This function does nothing.
*/
virtual void smooth (const unsigned int level,
+
/**
* Base class used to declare the operations needed by a concrete class
* implementing prolongation and restriction of vectors in the multigrid
* 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<double> &dst,
const Vector<double> &src) const = 0;
* 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<double> &dst,
const Vector<double> &src) const = 0;
};
/**
- * 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 VECTOR>
-class MGVector : public Subscriptor,
- public vector<VECTOR>
+template<class VECTOR = Vector<double> >
+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<T>#,
+ * #clear# will set all entries
+ * to zero, while if
+ * #VECTOR==std::vector<T>#,
+ * #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<VECTOR> vectors;
};
* 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 SOLVER, class MATRIX, class PRECOND>
* preconditioning method for later
* use.
*/
- MGCoarseGridLACIteration(SOLVER&, const MATRIX&, const PRECOND&);
+ MGCoarseGridLACIteration (SOLVER &,
+ const MATRIX &,
+ const PRECOND &);
/**
* Implementation of the abstract
/**
* Reference to the matrix.
*/
- const MATRIX& matrix;
+ const SmartPointer<const MATRIX> matrix;
/**
* Reference to the preconditioner.
+
/**
- * Basic class for multigrid preconditioning.
+ * Basic class for preconditioning by multigrid.
*
* The functionality of the multigrid method is restricted to defect
- * correction. It is <B>not</B> iterative and the start solution is
- * always zero. Since by this <I>u<SUP>E</SUP><SUB>l</SUB></I> and
- * <I>u<SUP>A</SUP><SUB>l</SUB></I> (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
* 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<Vector<double> > d;
+ unsigned int maxlevel;
/**
- * Auxiliary vector, solution.
+ * Level for coarse grid solution.
*/
- MGVector<Vector<double> > s;
+ unsigned int minlevel;
+
/**
- * Prolongation and restriction object.
+ * Auxiliary vector.
*/
- SmartPointer<const MGTransferBase> transfer;
+ MGVector<Vector<double> > defect;
- private:
/**
* Auxiliary vector.
*/
- Vector<double> t;
+ MGVector<Vector<double> > solution;
/**
- * Exception.
+ * Prolongation and restriction object.
*/
- DeclException2(ExcSwitchedLevels, int, int,
- << "minlevel and maxlevel switched, should be: "
- << arg1 << "<=" << arg2);
-
- protected:
+ SmartPointer<const MGTransferBase> 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
* This is implemented in a
* derived class.
*/
- virtual void level_vmult(unsigned int level,
- Vector<double>& dst,
- const Vector<double>& src,
- const Vector<double>& rhs) = 0;
-
+ virtual void level_vmult (const unsigned int level,
+ Vector<double> &dst,
+ const Vector<double> &src,
+ const Vector<double> &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<double> t;
};
* 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 MG, class VECTOR>
+template<class MG, class VECTOR = Vector<double> >
class PreconditionMG
{
public:
* 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<MG> multigrid;
+
/**
* The pre-smoothing object.
*/
SmartPointer<const MGSmootherBase> pre;
+
/**
* The post-smoothing object.
*/
SmartPointer<const MGSmootherBase> post;
+
/**
* The coarse grid solver.
*/
-
-template<class SOLVER, class MATRIX, class PRECOND>
-MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>::
-MGCoarseGridLACIteration(SOLVER& s, const MATRIX& m, const PRECOND& p)
- :
- solver(s), matrix(m), precondition(p)
-{}
-
-template<class SOLVER, class MATRIX, class PRECOND>
-void
-MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>::operator()
-(unsigned int, Vector<double>& dst, const Vector<double>& src) const
-{
- solver.solve(matrix, dst, src, precondition);
-}
-
-
-
-
/* ------------------------------------------------------------------- */
template<class VECTOR>
-MGVector<VECTOR>::MGVector(unsigned int min, unsigned int max)
+MGVector<VECTOR>::MGVector(const unsigned int min,
+ const unsigned int max)
:
- vector<VECTOR>(max-min+1),
- minlevel(min)
-{}
+ minlevel(min),
+ vectors(max-min+1)
+{};
+
template<class VECTOR>
-VECTOR&
-MGVector<VECTOR>::operator[](unsigned int i)
+VECTOR &
+MGVector<VECTOR>::operator[] (const unsigned int i)
{
- Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
- return vector<VECTOR>::operator[](i-minlevel);
+ Assert((i>=minlevel) && (i<minlevel+vectors.size()),
+ ExcIndexRange (i, minlevel, minlevel+vectors.size()));
+ return vectors[i-minlevel];
}
template<class VECTOR>
-const VECTOR&
-MGVector<VECTOR>::operator[](unsigned int i) const
+const VECTOR &
+MGVector<VECTOR>::operator[] (const unsigned int i) const
{
- Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
- return vector<VECTOR>::operator[](i-minlevel);
+ Assert((i>=minlevel) && (i<minlevel+vectors.size()),
+ ExcIndexRange (i, minlevel, minlevel+vectors.size()));
+ return vectors[i-minlevel];
}
+
template<class VECTOR>
void
MGVector<VECTOR>::clear()
{
typename vector<VECTOR>::iterator v;
- for (v = begin(); v != end(); ++v)
+ for (v = vectors.begin(); v != vectors.end(); ++v)
v->clear();
}
{
Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
return vector<MATRIX>::operator[](i-minlevel);
-}
+};
+
template<class MATRIX>
{
Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
return vector<MATRIX>::operator[](i-minlevel);
+};
+
+
+
+/* ------------------ Functions for MGCoarseGridLACIteration ------------ */
+
+
+template<class SOLVER, class MATRIX, class PRECOND>
+MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>::
+MGCoarseGridLACIteration(SOLVER &s,
+ const MATRIX &m,
+ const PRECOND &p)
+ :
+ solver(s),
+ matrix(&m),
+ precondition(p)
+{};
+
+
+
+template<class SOLVER, class MATRIX, class PRECOND>
+void
+MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>::
+operator() (const unsigned int /* level */,
+ Vector<double> &dst,
+ const Vector<double> &src) const
+{
+ solver.solve(*matrix, dst, src, precondition);
}
+
+
+
/* ------------------------------------------------------------------- */
template<class MG, class VECTOR>
-PreconditionMG<MG, VECTOR>::PreconditionMG(MG& mg,
- const MGSmootherBase& pre,
- const MGSmootherBase& post,
- const MGCoarseGridSolver& coarse)
+PreconditionMG<MG, VECTOR>::PreconditionMG(MG &mg,
+ const MGSmootherBase &pre,
+ const MGSmootherBase &post,
+ const MGCoarseGridSolver &coarse)
:
multigrid(&mg),
pre(&pre),
template<class MG, class VECTOR>
void
-PreconditionMG<MG, VECTOR>::operator() (VECTOR& dst,
- const VECTOR& src) const
+PreconditionMG<MG,VECTOR>::operator() (VECTOR &dst,
+ const VECTOR &src) const
{
multigrid->copy_to_mg(src);
multigrid->vcycle(*pre, *post, *coarse);
// $Id$
#include <numerics/multigrid.h>
+#include <algorithm>
#include <fstream>
-inline int minimum(int a, int b)
-{
- return ((a<b) ? a : b);
-}
+
+
template <int dim>
MG<dim>::MG(const MGDoFHandler<dim>& 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 <int dim>
template <typename number>
// set the elements of the vectors
// on all levels to zero
- d.clear();
+ defect.clear();
// hanging_nodes.condense(src);
vector<int> index(fe_dofs);
dc->get_dof_indices(index);
c->get_mg_dof_indices(mgindex);
for (unsigned int i=0;i<fe_dofs;++i)
- d[c->level()](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<dim>::faces_per_cell; ++face_n)
{
face->get_mg_dof_indices(mgfaceindex);
for (unsigned int i=0;i<face_dofs;++i)
- d[c->level()](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)
// }
}
-template <int dim> template <typename number>
+template <int dim>
+template <typename number>
void
MG<dim>::copy_from_mg(Vector<number>& dst) const
{
// }
DoFHandler<dim>::active_cell_iterator dc = dofs->DoFHandler<dim>::begin_active();
- for (MGDoFHandler<dim>::active_cell_iterator c = dofs->begin_active()
- ; c != dofs->end() ; ++c, ++dc)
+ for (MGDoFHandler<dim>::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;i<fe_dofs;++i)
- dst(index[i]) = s[c->level()](mgindex[i]);
+ dst(index[i]) = solution[c->level()](mgindex[i]);
}
hanging_nodes.set_zero(dst);
// {
// }
}
+
+
template <int dim>
void
MG<dim>::level_vmult(unsigned int l,
// $Id$
#include <numerics/multigrid.h>
+#include <algorithm>
#include <fstream>
-inline int minimum(int a, int b)
-{
- return ((a<b) ? a : b);
-}
+
+
template <int dim>
MG<dim>::MG(const MGDoFHandler<dim>& 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 <int dim>
template <typename number>
// set the elements of the vectors
// on all levels to zero
- d.clear();
+ defect.clear();
// hanging_nodes.condense(src);
vector<int> index(fe_dofs);
dc->get_dof_indices(index);
c->get_mg_dof_indices(mgindex);
for (unsigned int i=0;i<fe_dofs;++i)
- d[c->level()](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<dim>::faces_per_cell; ++face_n)
{
face->get_mg_dof_indices(mgfaceindex);
for (unsigned int i=0;i<face_dofs;++i)
- d[c->level()](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)
// }
}
-template <int dim> template <typename number>
+template <int dim>
+template <typename number>
void
MG<dim>::copy_from_mg(Vector<number>& dst) const
{
// }
DoFHandler<dim>::active_cell_iterator dc = dofs->DoFHandler<dim>::begin_active();
- for (MGDoFHandler<dim>::active_cell_iterator c = dofs->begin_active()
- ; c != dofs->end() ; ++c, ++dc)
+ for (MGDoFHandler<dim>::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;i<fe_dofs;++i)
- dst(index[i]) = s[c->level()](mgindex[i]);
+ dst(index[i]) = solution[c->level()](mgindex[i]);
}
hanging_nodes.set_zero(dst);
// {
// }
}
+
+
template <int dim>
void
MG<dim>::level_vmult(unsigned int l,
#include <iostream>
#include <cmath>
+
+
+//TODO: this function is only for debugging purposes and should be removed sometimes
template<class VECTOR>
+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<n;++i)
{
for (unsigned int j=0;j<n;++j)
-// s << ((float)i)/n << '\t' << ((float)j)/n << '\t' << v(k++) << endl;
s << '\t' << v(k++);
s << endl;
}
s << endl;
}
+
+
MGBase::~MGBase ()
{};
+//////////////////////////////////////////////////////////////////////
+
+
+MGSmootherBase::~MGSmootherBase()
+{};
+
+
+//////////////////////////////////////////////////////////////////////
+
+
+void
+MGSmootherIdentity::smooth (const unsigned int,
+ Vector<double> &,
+ const Vector<double> &) 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<double> &,
- const Vector<double> &) const
-{}
#include <vector>
#include <base/smartpointer.h>
+
+
/**
- * 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();
{
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<double> &u,
const Vector<double> &rhs) const = 0;
-
};
{
public:
/**
- * Implementation of the interface in #MGSmootherBase#.
+ * Implementation of the
+ * interface in #MGSmootherBase#.
* This function does nothing.
*/
virtual void smooth (const unsigned int level,
+
/**
* Base class used to declare the operations needed by a concrete class
* implementing prolongation and restriction of vectors in the multigrid
* 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<double> &dst,
const Vector<double> &src) const = 0;
* 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<double> &dst,
const Vector<double> &src) const = 0;
};
/**
- * 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 VECTOR>
-class MGVector : public Subscriptor,
- public vector<VECTOR>
+template<class VECTOR = Vector<double> >
+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<T>#,
+ * #clear# will set all entries
+ * to zero, while if
+ * #VECTOR==std::vector<T>#,
+ * #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<VECTOR> vectors;
};
* 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 SOLVER, class MATRIX, class PRECOND>
* preconditioning method for later
* use.
*/
- MGCoarseGridLACIteration(SOLVER&, const MATRIX&, const PRECOND&);
+ MGCoarseGridLACIteration (SOLVER &,
+ const MATRIX &,
+ const PRECOND &);
/**
* Implementation of the abstract
/**
* Reference to the matrix.
*/
- const MATRIX& matrix;
+ const SmartPointer<const MATRIX> matrix;
/**
* Reference to the preconditioner.
+
/**
- * Basic class for multigrid preconditioning.
+ * Basic class for preconditioning by multigrid.
*
* The functionality of the multigrid method is restricted to defect
- * correction. It is <B>not</B> iterative and the start solution is
- * always zero. Since by this <I>u<SUP>E</SUP><SUB>l</SUB></I> and
- * <I>u<SUP>A</SUP><SUB>l</SUB></I> (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
* 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<Vector<double> > d;
+ unsigned int maxlevel;
/**
- * Auxiliary vector, solution.
+ * Level for coarse grid solution.
*/
- MGVector<Vector<double> > s;
+ unsigned int minlevel;
+
/**
- * Prolongation and restriction object.
+ * Auxiliary vector.
*/
- SmartPointer<const MGTransferBase> transfer;
+ MGVector<Vector<double> > defect;
- private:
/**
* Auxiliary vector.
*/
- Vector<double> t;
+ MGVector<Vector<double> > solution;
/**
- * Exception.
+ * Prolongation and restriction object.
*/
- DeclException2(ExcSwitchedLevels, int, int,
- << "minlevel and maxlevel switched, should be: "
- << arg1 << "<=" << arg2);
-
- protected:
+ SmartPointer<const MGTransferBase> 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
* This is implemented in a
* derived class.
*/
- virtual void level_vmult(unsigned int level,
- Vector<double>& dst,
- const Vector<double>& src,
- const Vector<double>& rhs) = 0;
-
+ virtual void level_vmult (const unsigned int level,
+ Vector<double> &dst,
+ const Vector<double> &src,
+ const Vector<double> &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<double> t;
};
* 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 MG, class VECTOR>
+template<class MG, class VECTOR = Vector<double> >
class PreconditionMG
{
public:
* 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<MG> multigrid;
+
/**
* The pre-smoothing object.
*/
SmartPointer<const MGSmootherBase> pre;
+
/**
* The post-smoothing object.
*/
SmartPointer<const MGSmootherBase> post;
+
/**
* The coarse grid solver.
*/
-
-template<class SOLVER, class MATRIX, class PRECOND>
-MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>::
-MGCoarseGridLACIteration(SOLVER& s, const MATRIX& m, const PRECOND& p)
- :
- solver(s), matrix(m), precondition(p)
-{}
-
-template<class SOLVER, class MATRIX, class PRECOND>
-void
-MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>::operator()
-(unsigned int, Vector<double>& dst, const Vector<double>& src) const
-{
- solver.solve(matrix, dst, src, precondition);
-}
-
-
-
-
/* ------------------------------------------------------------------- */
template<class VECTOR>
-MGVector<VECTOR>::MGVector(unsigned int min, unsigned int max)
+MGVector<VECTOR>::MGVector(const unsigned int min,
+ const unsigned int max)
:
- vector<VECTOR>(max-min+1),
- minlevel(min)
-{}
+ minlevel(min),
+ vectors(max-min+1)
+{};
+
template<class VECTOR>
-VECTOR&
-MGVector<VECTOR>::operator[](unsigned int i)
+VECTOR &
+MGVector<VECTOR>::operator[] (const unsigned int i)
{
- Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
- return vector<VECTOR>::operator[](i-minlevel);
+ Assert((i>=minlevel) && (i<minlevel+vectors.size()),
+ ExcIndexRange (i, minlevel, minlevel+vectors.size()));
+ return vectors[i-minlevel];
}
template<class VECTOR>
-const VECTOR&
-MGVector<VECTOR>::operator[](unsigned int i) const
+const VECTOR &
+MGVector<VECTOR>::operator[] (const unsigned int i) const
{
- Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
- return vector<VECTOR>::operator[](i-minlevel);
+ Assert((i>=minlevel) && (i<minlevel+vectors.size()),
+ ExcIndexRange (i, minlevel, minlevel+vectors.size()));
+ return vectors[i-minlevel];
}
+
template<class VECTOR>
void
MGVector<VECTOR>::clear()
{
typename vector<VECTOR>::iterator v;
- for (v = begin(); v != end(); ++v)
+ for (v = vectors.begin(); v != vectors.end(); ++v)
v->clear();
}
{
Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
return vector<MATRIX>::operator[](i-minlevel);
-}
+};
+
template<class MATRIX>
{
Assert((i>=minlevel)&&(i<minlevel+size()),ExcIndexRange(i,minlevel,minlevel+size()));
return vector<MATRIX>::operator[](i-minlevel);
+};
+
+
+
+/* ------------------ Functions for MGCoarseGridLACIteration ------------ */
+
+
+template<class SOLVER, class MATRIX, class PRECOND>
+MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>::
+MGCoarseGridLACIteration(SOLVER &s,
+ const MATRIX &m,
+ const PRECOND &p)
+ :
+ solver(s),
+ matrix(&m),
+ precondition(p)
+{};
+
+
+
+template<class SOLVER, class MATRIX, class PRECOND>
+void
+MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>::
+operator() (const unsigned int /* level */,
+ Vector<double> &dst,
+ const Vector<double> &src) const
+{
+ solver.solve(*matrix, dst, src, precondition);
}
+
+
+
/* ------------------------------------------------------------------- */
template<class MG, class VECTOR>
-PreconditionMG<MG, VECTOR>::PreconditionMG(MG& mg,
- const MGSmootherBase& pre,
- const MGSmootherBase& post,
- const MGCoarseGridSolver& coarse)
+PreconditionMG<MG, VECTOR>::PreconditionMG(MG &mg,
+ const MGSmootherBase &pre,
+ const MGSmootherBase &post,
+ const MGCoarseGridSolver &coarse)
:
multigrid(&mg),
pre(&pre),
template<class MG, class VECTOR>
void
-PreconditionMG<MG, VECTOR>::operator() (VECTOR& dst,
- const VECTOR& src) const
+PreconditionMG<MG,VECTOR>::operator() (VECTOR &dst,
+ const VECTOR &src) const
{
multigrid->copy_to_mg(src);
multigrid->vcycle(*pre, *post, *coarse);
#include <iostream>
#include <cmath>
+
+
+//TODO: this function is only for debugging purposes and should be removed sometimes
template<class VECTOR>
+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<n;++i)
{
for (unsigned int j=0;j<n;++j)
-// s << ((float)i)/n << '\t' << ((float)j)/n << '\t' << v(k++) << endl;
s << '\t' << v(k++);
s << endl;
}
s << endl;
}
+
+
MGBase::~MGBase ()
{};
+//////////////////////////////////////////////////////////////////////
+
+
+MGSmootherBase::~MGSmootherBase()
+{};
+
+
+//////////////////////////////////////////////////////////////////////
+
+
+void
+MGSmootherIdentity::smooth (const unsigned int,
+ Vector<double> &,
+ const Vector<double> &) 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<double> &,
- const Vector<double> &) const
-{}