#include <base/subscriptor.h>
#include <base/smartpointer.h>
#include <lac/vector.h>
-#include <lac/sparsity_pattern.h>
-
#include <vector>
-class MGSmootherBase;
-
/**
* An array with an object for each level. The purpose of this 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. 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
- * @p{MGCoarseGridLACIteration}.
- *
- * Note that the information
- * about the matrix is removed to
- * that class.
- */
- virtual void operator() (const unsigned int level,
- Vector<double> &dst,
- const Vector<double> &src) const = 0;
-};
-
/**
* Coarse grid solver using LAC iterative methods.
*
* @author Guido Kanschat, 1999
*/
-template<class SOLVER, class MATRIX, class PRECOND>
-class MGCoarseGridLACIteration : public MGCoarseGridSolver
+template<class SOLVER, class MATRIX, class PRECOND, class VECTOR = Vector<double> >
+class MGCoarseGridLACIteration : public Subscriptor
{
public:
/**
* matrix, vectors and
* preconditioner.
*/
- virtual void operator() (const unsigned int level,
- Vector<double> &dst,
- const Vector<double> &src) const;
+ void operator() (const unsigned int level,
+ VECTOR &dst,
+ const VECTOR &src) const;
private:
/**
* Reference to the solver.
*
* @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
- * @p{to_level-1} to level
- * @p{to_level}. The previous
- * content of @p{dst} is
- * overwritten.
- *
- * @p{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 @p{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<double> &dst,
- const Vector<double> &src) const = 0;
-
- /**
- * Restrict a vector from level
- * @p{from_level} to level
- * @p{from_level-1} and add this
- * restriction to
- * @p{dst}. Obviously, if the
- * refined region on level
- * @p{from_level} is smaller than
- * that on level @p{from_level-1},
- * some degrees of freedom in
- * @p{dst} are not covered and will
- * not be altered. For the other
- * degress of freedom, the result
- * of the restriction is added.
- *
- * @p{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 @p{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<double> &dst,
- const Vector<double> &src) const = 0;
-};
-
-
-/**
- * Basic class for preconditioning by multigrid.
- *
- * The functionality of the multigrid method is restricted to defect
- * correction. It is @em{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 @p{precondition} is the actual multigrid method and
- * makes use of several operations to be implemented in derived
- * classes. It takes a defect in @p{src} and the result is the multigrid
- * preconditioned defect in @p{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 @p{d} is properly filled
- * with the residual in the outer
- * defect correction
- * scheme. After execution of
- * @p{vcycle()}, the result is in
- * the vector @p{s}. We propose to
- * write functions @p{copy_from_mg}
- * and @p{copy_to_mg} in derived
- * classes of @p{MGBase} to access
- * @p{d} and @p{s}.
- *
- * The actual work for this
- * function is done in
- * @p{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.
- */
- MGLevelObject<Vector<double> > defect;
-
- /**
- * Auxiliary vector.
- */
- MGLevelObject<Vector<double> > solution;
-
- /**
- * Auxiliary vector.
- */
- MGLevelObject<Vector<double> > t;
-
- /**
- * Prolongation and restriction object.
- */
- SmartPointer<const MGTransferBase> transfer;
-
- /**
- * The actual v-cycle multigrid method.
- * @p{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 @p{level-1},
- * unless we are on @p{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 @p{vmult} on all
- * cells of a level.
- * This is implemented in a
- * derived class.
- */
- virtual void level_vmult (const unsigned int level,
- Vector<double> &dst,
- const Vector<double> &src,
- const Vector<double> &rhs) = 0;
-
- /**
- * Additional multiplication on a
- * refinement edge. Using DGFEM,
- * the application of the global
- * matrix to a function on the
- * fine level produces results on
- * the surrounding cells of the
- * coarse level. Therefore,
- * additionally to the fine level
- * matrix, we need an operator
- * from the fine level to the
- * coarse level, taking care of
- * this multiplication.
- *
- * Like @ref{level_vmult}, it is
- * expected to add the negative
- * product of the matrix and
- * @p{src} to @p{dst}. Here,
- * @p{src} is a vector on the
- * fine level and @p{dst} is a
- * vector on the coarse level.
- *
- * This function has an empty
- * implementation here and must
- * be overloaded in a multigrid
- * class for discontinuous
- * methods.
- */
- virtual void edge_vmult (const unsigned int level,
- Vector<double> &dst,
- const Vector<double> &src);
-
- /**
- * Print a level vector using
- * @ref{DoFHandler}.
- */
- virtual void print_vector (const unsigned int level,
- const Vector<double>& v,
- const char* name) const = 0;
-
-};
-
-
-/**
- * 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 @p{MG} is required to inherit @p{MGBase}.
- * Furthermore, it needs functions @p{void copy_to_mg(const VECTOR&)}
- * to store @p{src} in the right hand side of the multi-level method and
- * @p{void copy_from_mg(VECTOR&)} to store the result of the v-cycle in @p{dst}.
- *
- * @author Guido Kanschat, 1999
- */
-template<class MG, class VECTOR = Vector<double> >
-class PreconditionMG : public Subscriptor
-{
- 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 @p{vcycle} function
- * of the @p{MG} object passed to
- * the constructor.
- *
- * This is the operator used by
- * LAC iterative solvers.
- */
- void vmult (VECTOR &dst,
- const VECTOR &src) const;
-
- private:
- /**
- * The multigrid object.
- */
- SmartPointer<MG> multigrid;
-
- /**
- * The pre-smoothing object.
- */
- SmartPointer<const MGSmootherBase> pre;
+// template <class VECTOR>
+// class MGTransferBase : public Subscriptor
+// {
+// public:
+// /**
+// * Destructor. Does nothing here, but
+// * needs to be declared virtual anyway.
+// */
+// virtual ~MGTransferBase();
+
+// /**
+// * Prolongate a vector from level
+// * @p{to_level-1} to level
+// * @p{to_level}. The previous
+// * content of @p{dst} is
+// * overwritten.
+// *
+// * @p{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 @p{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
+// * @p{from_level} to level
+// * @p{from_level-1} and add this
+// * restriction to
+// * @p{dst}. Obviously, if the
+// * refined region on level
+// * @p{from_level} is smaller than
+// * that on level @p{from_level-1},
+// * some degrees of freedom in
+// * @p{dst} are not covered and will
+// * not be altered. For the other
+// * degress of freedom, the result
+// * of the restriction is added.
+// *
+// * @p{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 @p{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;
+// };
- /**
- * The post-smoothing object.
- */
- SmartPointer<const MGSmootherBase> post;
-
- /**
- * The coarse grid solver.
- */
- SmartPointer<const MGCoarseGridSolver> coarse;
-};
/* ------------------------------------------------------------------- */
/* ------------------ Functions for MGCoarseGridLACIteration ------------ */
-template<class SOLVER, class MATRIX, class PRECOND>
-MGCoarseGridLACIteration<SOLVER,MATRIX,PRECOND>::MGCoarseGridLACIteration(SOLVER &s,
- const MATRIX &m,
- const PRECOND &p)
+template<class SOLVER, class MATRIX, class PRECOND, class VECTOR>
+MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND, VECTOR>
+::MGCoarseGridLACIteration(SOLVER& s,
+ const MATRIX &m,
+ const PRECOND &p)
:
solver(s),
matrix(&m),
{};
-template<class SOLVER, class MATRIX, class PRECOND>
+template<class SOLVER, class MATRIX, class PRECOND, class VECTOR>
void
-MGCoarseGridLACIteration<SOLVER,MATRIX,PRECOND>::operator() (const unsigned int /* level */,
- Vector<double> &dst,
- const Vector<double> &src) const
+MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND, VECTOR>
+::operator() (const unsigned int /* level */,
+ VECTOR &dst,
+ const VECTOR &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)
- :
- multigrid(&mg),
- pre(&pre),
- post(&post),
- coarse(&coarse)
-{}
-
-
-template<class MG, class VECTOR>
-void
-PreconditionMG<MG,VECTOR>::vmult (VECTOR &dst,
- const VECTOR &src) const
-{
- multigrid->copy_to_mg(src);
- multigrid->vcycle(*pre, *post, *coarse);
- multigrid->copy_from_mg(dst);
-}
-
-
/*---------------------------- mgbase.h ---------------------------*/
#endif
#ifndef __deal2__mg_dof_tools_h
#define __deal2__mg_dof_tools_h
+template <class Object> class MGLevelObject;
+template <int dim> class MGDoFHandler;
+template <typename number> class Vector;
+template <typename number> class BlockVector;
+template <class number> class FullMatrix;
#include <base/config.h>
* for more information.
*
* All member functions are static, so there is no need to create an
- * object of class @p{MGDoFTools}.
+ * object of class @p{MGTools}.
*
* @author Wolfgang Bangerth, Guido Kanschat, 1999, 2000
*/
-class MGDoFTools
+class MGTools
{
public:
/**
* hanging nodes here, since only
* one level is considered.
*/
- template <int dim>
+ template <int dim, class SparsityPattern>
static void
make_sparsity_pattern (const MGDoFHandler<dim> &dof_handler,
SparsityPattern &sparsity,
* @see make_sparsity_pattern
* $see DoFTools
*/
- template <int dim>
+ template <int dim, class SparsityPattern>
static void
make_flux_sparsity_pattern (const MGDoFHandler<dim> &dof_handler,
SparsityPattern &sparsity,
* @see{make_flux_sparsity_pattern}
* @see{DoFTools}
*/
- template <int dim>
+ template <int dim, class SparsityPattern>
static void
make_flux_sparsity_pattern_edge (const MGDoFHandler<dim> &dof_handler,
SparsityPattern &sparsity,
* for the couplings occuring in
* fluxes.
*/
- template <int dim>
+ template <int dim, class SparsityPattern>
static void
make_flux_sparsity_pattern (const MGDoFHandler<dim> &dof,
SparsityPattern &sparsity,
const FullMatrix<double>& int_mask,
const FullMatrix<double>& flux_mask);
+ /**
+ * Count the dofs component-wise
+ * on each level.
+ *
+ * Result is a vector containing
+ * for each level a vector
+ * containing the number of dofs
+ * for each component (access is
+ * @p{result[level][component]}).
+ */
+ template <int dim>
+ static void count_dofs_per_component (const MGDoFHandler<dim>& mg_dof,
+ std::vector<std::vector<unsigned int> >& result);
+
+
+ /**
+ * Ajust vectors on all levels to
+ * correct size. Here, we just
+ * count the numbers of degrees
+ * of freedom on each level and
+ * @p{reinit} each level vector
+ * to this length.
+ */
+ template <int dim, typename number>
+ static void
+ reinit_vector (const MGDoFHandler<dim>& mg_dof,
+ MGLevelObject<Vector<number> >& vector);
+
+ /**
+ * Adjust block-vectors on all
+ * levels to correct size. Count
+ * the numbers of degrees of
+ * freedom on each level
+ * component-wise. Then, assign
+ * each block of @p{vector} the
+ * corresponding size.
+ *
+ * The boolean field @p{selected}
+ * allows restricting this
+ * operation to certain
+ * components. In this case,
+ * @p{vector} will only have as
+ * many blocks as there are true
+ * values in @p{selected} (no
+ * blocks of length zero are
+ * padded in).
+ *
+ * Degrees of freedom must be
+ * sorted by component in order
+ * to obtain reasonable results
+ * from this function.
+ */
+ template <int dim, typename number>
+ static void
+ reinit_vector (const MGDoFHandler<dim>& mg_dof,
+ MGLevelObject<BlockVector<number> >& vector,
+ const vector<bool>& selected);
};
template <int dim> class MGDoFHandler;
-/**
- * Abstract base class for multigrid smoothers.
- * In fact, this class only provides the interface of the smoothing function.
- * Using @p{deal.II} grid handling, @p{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 @p{u} on
- * the given level. If $S$ is the
- * smoothing operator, then this
- * function should do the
- * following operation:
- * $u += S (rhs - Au)$, where @p{u} and
- * @p{rhs} are the input parameters.
- *
- * This function should keep the
- * interior level boundary
- * values, so you may want to
- * call @p{set_zero_interior_boundary}
- * in the derived class
- * @p{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;
-};
-
-
/**
* Smoother doing nothing. This class is not useful for many applications other
* than for testing some multigrid procedures. Also some applications might
*
* @author Guido Kanschat, 1999
*/
-class MGSmootherIdentity : public MGSmootherBase
+class MGSmootherIdentity : public Subscriptor
{
public:
/**
* Implementation of the
- * interface in @p{MGSmootherBase}.
+ * interface for @p{Multigrid}.
* This function does nothing,
* which by comparison with the
* definition of this function
* operator equals the null
* operator.
*/
- virtual void smooth (const unsigned int level,
- Vector<double> &u,
- const Vector<double> &rhs) const;
+ template <class VECTOR>
+ 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 @p{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.
+ * Base class for multigrid smoothers. It 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
+class MGSmootherContinuous : public Subscriptor
{
private:
/**
* Default constructor. Made private
* to prevent it being called, which
* is necessary since this could
- * circumpass the set-up of the list
+ * circumvent the set-up of the list
* if interior boundary dofs.
*/
- MGSmoother ();
+ MGSmootherContinuous ();
public:
* needed by the function
* @p{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.
+ * 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. Really amusing
+ * about this text is, that there
+ * is no 1d implementation.
*/
- MGSmoother (const MGDoFHandler<1> &mg_dof,
+ MGSmootherContinuous (const MGDoFHandler<1> &mg_dof,
unsigned int steps);
/**
* steps to be executed by @p{smooth}.
*/
template <int dim>
- MGSmoother (const MGDoFHandler<dim> &mg_dof,
+ MGSmootherContinuous (const MGDoFHandler<dim> &mg_dof,
unsigned int steps);
/**
* has no interior boundaries, this
* function does nothing in this case.
*/
- void set_zero_interior_boundary (const unsigned int level,
- Vector<double> &u) const;
+ template <class VECTOR>
+ 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?
*/
* Number of smoothing steps.
*/
unsigned int steps;
+
/**
* For each level, we store a list of
* degree of freedom indices which are
/**
- * Implementation of a smoother using matrix builtin relaxation methods.
+ * A smoother using matrix builtin relaxation methods.
+ *
+ * Additionally to smoothing with the matrix built-in relaxation
+ * scheme, hanging nodes of continuous finite elements are handled
+ * properly (at least in a future version).
+ *
+ * The library contains instantiation for @p{SparseMatrix<.>} and
+ * @p{Vector<.>}, where the template arguments are all combinations of
+ * @p{float} and @p{double}. Additional instantiations may be created
+ * by including the file mg_smoother.templates.h.
*
* @author Wolfgang Bangerth, Guido Kanschat, 1999
*/
-template<typename number>
-class MGSmootherRelaxation : public MGSmoother
+template<class MATRIX, class VECTOR>
+class MGSmootherRelaxation : public MGSmootherContinuous
{
public:
/**
* Type of the smoothing
* function of the matrix.
*/
- typedef void
- (SparseMatrix<number>::* function_ptr)(Vector<double>&, const Vector<double>&,
- typename SparseMatrix<number>::value_type) const;
+ typedef void (MATRIX::* function_ptr)(VECTOR&,
+ const VECTOR&,
+ typename MATRIX::value_type) const;
/**
* Constructor.
* The constructor uses an @p{MGDoFHandler} to initialize
* data structures for interior level boundary handling
- * in @p{MGSmoother}.
+ * in @p{MGSmootherContinuous}.
*
* Furthermore, it takes a pointer to the
* level matrices and their smoothing function.
*/
template<int dim>
- MGSmootherRelaxation(const MGDoFHandler<dim> &mg_dof,
- const MGLevelObject<SparseMatrix<number> > &matrix,
- const function_ptr function,
- const unsigned int steps,
- const double omega = 1.);
+ MGSmootherRelaxation(const MGDoFHandler<dim>& mg_dof,
+ const MGLevelObject<MATRIX>& matrix,
+ const function_ptr function,
+ const unsigned int steps,
+ const double omega = 1.);
/**
- * Implementation of the interface in @p{MGSmootherBase}.
- * We use the SOR method in @p{SparseMatrix} for the real work
- * and find a way to keep the boundary values.
+ * We use the relaxation method in
+ * @p{MATRIX} for the real
+ * work and find a way to keep
+ * the boundary values.
*/
- virtual void smooth (const unsigned int level,
- Vector<double> &u,
- const Vector<double> &rhs) const;
+ void smooth (const unsigned int level,
+ VECTOR& u,
+ const VECTOR& rhs) const;
private:
/**
* Pointer to the matrices.
*/
- SmartPointer<const MGLevelObject<SparseMatrix<number> > > matrix;
+ SmartPointer<const MGLevelObject<MATRIX> > matrix;
/**
* Pointer to the relaxation function.
/* ------------------------------- Inline functions -------------------------- */
+template <class VECTOR>
+inline void
+MGSmootherIdentity::smooth (const unsigned int, VECTOR&, const VECTOR&) const
+{};
+
+/*----------------------------------------------------------------------*/
+
inline
void
-MGSmoother::set_steps(unsigned int i)
+MGSmootherContinuous::set_steps(unsigned int i)
{
steps = i;
}
inline
unsigned int
-MGSmoother::get_steps() const
+MGSmootherContinuous::get_steps() const
{
return steps;
}
// $Id$
-template<typename number>
+template<class MATRIX, class VECTOR>
template<int dim>
-MGSmootherRelaxation<number>::MGSmootherRelaxation (const MGDoFHandler<dim> &mg_dof,
- const MGLevelObject<SparseMatrix<number> > &matrix,
- const function_ptr relaxation,
- const unsigned int steps,
- const double omega)
+MGSmootherRelaxation<MATRIX, VECTOR>
+::MGSmootherRelaxation (const MGDoFHandler<dim>& mg_dof,
+ const MGLevelObject<MATRIX>& matrix,
+ const function_ptr relaxation,
+ const unsigned int steps,
+ const double omega)
:
- MGSmoother(mg_dof, steps),
+ MGSmootherContinuous(mg_dof, steps),
matrix(&matrix),
relaxation(relaxation),
omega(omega)
{};
-template<typename number>
+template<class MATRIX, class VECTOR>
void
-MGSmootherRelaxation<number>::smooth (const unsigned int level,
- Vector<double> &u,
- const Vector<double> &rhs) const
+MGSmootherRelaxation<MATRIX, VECTOR>
+::smooth (const unsigned int level,
+ VECTOR& u,
+ const VECTOR& rhs) const
{
h1.reinit(u);
h2.reinit(u);
--- /dev/null
+//----------------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 1998, 1999, 2000, 2001 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//----------------------------------------------------------------------------
+
+#ifndef __deal2__mg_tools_templates_h
+#define __deal2__mg_tools_templates_h
+
+#include <dofs/dof_constraints.h>
+#include <numerics/data_out.h>
+#include <multigrid/multigrid.h>
+#include <algorithm>
+#include <fstream>
+
+#include <lac/sparse_matrix.h>
+
+
+//TODO:[?] This function needs to be specially implemented, since in 2d mode we use faces
+#if deal_II_dimension == 1
+
+template <int dim, class InVector>
+void
+MGTransferPrebuilt::copy_to_mg (const MGDoFHandler<dim>&,
+ MGLevelObject<Vector<double> >&,
+ const InVector&) const
+{
+ Assert(false, ExcNotImplemented());
+}
+
+#else
+
+
+template <int dim, class InVector>
+void
+MGTransferPrebuilt::copy_to_mg (const MGDoFHandler<dim>& mg_dof_handler,
+ MGLevelObject<Vector<double> >& dst,
+ const InVector& src) const
+{
+
+ // Make src a real finite element function
+// InVector src = osrc;
+// constraints->distribute(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
+ unsigned int minlevel = dst.get_minlevel();
+ unsigned int maxlevel = dst.get_maxlevel();
+
+ dst.clear();
+
+ for (unsigned int l=minlevel;l<=maxlevel;++l)
+ dst[l].reinit(sizes[l]);
+
+ std::vector<unsigned int> global_dof_indices (dofs_per_cell);
+ std::vector<unsigned int> level_dof_indices (dofs_per_cell);
+ std::vector<unsigned int> 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<int>(minlevel); --level)
+ {
+ typename MGDoFHandler<dim>::active_cell_iterator
+ level_cell = mg_dof_handler.begin_active(level);
+ const typename MGDoFHandler<dim>::active_cell_iterator
+ level_end = mg_dof_handler.end_active(level);
+
+//TODO:[?] Treat hanging nodes properly
+// The single-level vector is not an FE-function, because the values at
+// hanging nodes are set to zero. This should be treated before the restriction.
+
+ // Compute coarse level right hand side
+ // by restricting from fine level.
+ for (; level_cell!=level_end; ++level_cell)
+ {
+ DoFObjectAccessor<dim, dim>& global_cell = *level_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<dofs_per_cell; ++i)
+ dst[level](level_dof_indices[i]) = src(global_dof_indices[i]);
+
+ for (unsigned int face_n=0; face_n<GeometryInfo<dim>::faces_per_cell; ++face_n)
+ {
+ const typename MGDoFHandler<dim>::face_iterator
+ face = level_cell->face(face_n);
+ if (face->has_children())
+ {
+ face->get_mg_dof_indices(level_face_indices);
+
+
+ // Delete values on refinement edge,
+ // since restriction will add them again.
+ for (unsigned int i=0; i<dofs_per_face; ++i)
+ dst[level](level_face_indices[i])
+ = 0.;
+ };
+ };
+ };
+ // for that part of the level
+ // which is further refined:
+ // get the defect by
+ // restriction of the defect on
+ // one level higher
+ if (static_cast<unsigned int>(level) < maxlevel)
+ {
+ restrict_and_add (level+1, dst[level], dst[level+1]);
+ }
+ };
+};
+
+#endif
+
+
+template <int dim, class OutVector>
+void
+MGTransferPrebuilt::copy_from_mg(const MGDoFHandler<dim>& mg_dof_handler,
+ OutVector &dst,
+ const MGLevelObject<Vector<double> > &src) const
+{
+ const unsigned int dofs_per_cell = mg_dof_handler.get_fe().dofs_per_cell;
+
+ std::vector<unsigned int> global_dof_indices (dofs_per_cell);
+ std::vector<unsigned int> level_dof_indices (dofs_per_cell);
+
+ typename MGDoFHandler<dim>::active_cell_iterator
+ level_cell = mg_dof_handler.begin_active();
+ const typename MGDoFHandler<dim>::active_cell_iterator
+ endc = mg_dof_handler.end();
+
+ // traverse all cells and copy the
+ // data appropriately to the output
+ // vector
+
+ // Is the level monotonuosly increasing?
+
+ for (; level_cell != endc; ++level_cell)
+ {
+ DoFObjectAccessor<dim, dim>& global_cell = *level_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; i<dofs_per_cell; ++i)
+ dst(global_dof_indices[i]) = src[level](level_dof_indices[i]);
+ };
+
+ // clear constrained nodes
+//TODO: constraints->set_zero(dst);
+}
+
+template <int dim, class OutVector>
+void
+MGTransferPrebuilt::copy_from_mg_add(const MGDoFHandler<dim>& mg_dof_handler,
+ OutVector &dst,
+ const MGLevelObject<Vector<double> > &src) const
+{
+ const unsigned int dofs_per_cell = mg_dof_handler.get_fe().dofs_per_cell;
+
+ std::vector<unsigned int> global_dof_indices (dofs_per_cell);
+ std::vector<unsigned int> level_dof_indices (dofs_per_cell);
+
+ typename MGDoFHandler<dim>::active_cell_iterator
+ level_cell = mg_dof_handler.begin_active();
+ const typename MGDoFHandler<dim>::active_cell_iterator
+ endc = mg_dof_handler.end();
+
+ // traverse all cells and copy the
+ // data appropriately to the output
+ // vector
+
+ // Is the level monotonuosly increasing?
+
+ for (; level_cell != endc; ++level_cell)
+ {
+ DoFObjectAccessor<dim, dim>& global_cell = *level_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);
+
+//TODO: Probably wrong for continuous elements
+
+ // copy level-wise data to
+ // global vector
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ dst(global_dof_indices[i]) += src[level](level_dof_indices[i]);
+ };
+
+ // clear constrained nodes
+//TODO: constraints->set_zero(dst);
+}
+
+
+//TODO:[?] This function needs to be specially implemented, since in 2d mode we use faces
+#if deal_II_dimension == 1
+
+template <int dim, class InVector>
+void
+MGTransferSelect::copy_to_mg (const MGDoFHandler<dim>&,
+ MGLevelObject<Vector<double> >&,
+ const InVector&) const
+{
+ Assert(false, ExcNotImplemented());
+}
+
+#else
+
+
+template <int dim, class InVector>
+void
+MGTransferSelect::copy_to_mg (const MGDoFHandler<dim>& mg_dof_handler,
+ MGLevelObject<Vector<double> >& dst,
+ const InVector& osrc) const
+{
+ // Make src a real finite element function
+ InVector src = osrc;
+// constraints->distribute(src);
+
+ const FiniteElement<dim>& fe = mg_dof_handler.get_fe();
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ const unsigned int dofs_per_face = fe.dofs_per_face;
+
+ // set the elements of the vectors
+ // on all levels to zero
+ unsigned int minlevel = dst.get_minlevel();
+ unsigned int maxlevel = dst.get_maxlevel();
+
+ dst.clear();
+
+ for (unsigned int l=minlevel;l<=maxlevel;++l)
+ dst[l].reinit(sizes[l][selected]);
+
+ std::vector<unsigned int> global_dof_indices (dofs_per_cell);
+ std::vector<unsigned int> level_dof_indices (dofs_per_cell);
+ std::vector<unsigned int> level_face_indices (dofs_per_face);
+
+ // Start of input vector inside a
+ // block vector.
+ const unsigned int start = component_start[selected];
+
+ // 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<int>(minlevel); --level)
+ {
+ // Start of treated component
+ // for this level
+ const unsigned int level_start = mg_component_start[level][selected];
+
+ typename MGDoFHandler<dim>::active_cell_iterator
+ level_cell = mg_dof_handler.begin_active(level);
+ const typename MGDoFHandler<dim>::active_cell_iterator
+ level_end = mg_dof_handler.end_active(level);
+
+//TODO:[?] Treat hanging nodes properly
+// The single-level vector is not an FE-function, because the values at
+// hanging nodes are set to zero. This should be treated before the restriction.
+
+ // Compute coarse level right hand side
+ // by restricting from fine level.
+ for (; level_cell!=level_end; ++level_cell)
+ {
+ DoFObjectAccessor<dim, dim>& global_cell = *level_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<dofs_per_cell; ++i)
+ {
+ if (fe.system_to_component_index(i).first == selected)
+ dst[level](level_dof_indices[i] - level_start)
+ = src(global_dof_indices[i] - start);
+ }
+
+ for (unsigned int face_n=0; face_n<GeometryInfo<dim>::faces_per_cell; ++face_n)
+ {
+ const typename MGDoFHandler<dim>::face_iterator
+ face = level_cell->face(face_n);
+ if (face->has_children())
+ {
+ face->get_mg_dof_indices(level_face_indices);
+
+
+ // Delete values on refinement edge,
+ // since restriction will add them again.
+// for (unsigned int i=0; i<dofs_per_face; ++i)
+// dst[level](level_face_indices[i])
+// = 0.;
+ };
+ };
+ };
+ // for that part of the level
+ // which is further refined:
+ // get the defect by
+ // restriction of the defect on
+ // one level higher
+ if (static_cast<unsigned int>(level) < maxlevel)
+ {
+ restrict_and_add (level+1, dst[level], dst[level+1]);
+ }
+ };
+};
+
+#endif
+
+
+template <int dim, class OutVector>
+void
+MGTransferSelect::copy_from_mg(const MGDoFHandler<dim>& mg_dof_handler,
+ OutVector &dst,
+ const MGLevelObject<Vector<double> > &src) const
+{
+
+ const FiniteElement<dim>& fe = mg_dof_handler.get_fe();
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+
+ std::vector<unsigned int> global_dof_indices (dofs_per_cell);
+ std::vector<unsigned int> level_dof_indices (dofs_per_cell);
+
+ typename MGDoFHandler<dim>::active_cell_iterator
+ level_cell = mg_dof_handler.begin_active();
+ const typename MGDoFHandler<dim>::active_cell_iterator
+ endc = mg_dof_handler.end();
+
+ // Start of input vector inside a
+ // block vector.
+ const unsigned int start = component_start[selected];
+
+ // traverse all cells and copy the
+ // data appropriately to the output
+ // vector
+
+ // Is the level monotonuosly increasing?
+
+ for (; level_cell != endc; ++level_cell)
+ {
+ DoFObjectAccessor<dim, dim>& global_cell = *level_cell;
+ const unsigned int level = level_cell->level();
+
+ // Start of treated component
+ // for this level
+ const unsigned int level_start = mg_component_start[level][selected];
+
+ // 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; i<dofs_per_cell; ++i)
+ if (fe.system_to_component_index(i).first == selected)
+ dst(global_dof_indices[i]-start)
+ = src[level](level_dof_indices[i]-level_start);
+ };
+
+ // clear constrained nodes
+//TODO: constraints->set_zero(dst);
+}
+
+template <int dim, class OutVector>
+void
+MGTransferSelect::copy_from_mg_add(const MGDoFHandler<dim>& mg_dof_handler,
+ OutVector &dst,
+ const MGLevelObject<Vector<double> > &src) const
+{
+
+ const FiniteElement<dim>& fe = mg_dof_handler.get_fe();
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+
+ std::vector<unsigned int> global_dof_indices (dofs_per_cell);
+ std::vector<unsigned int> level_dof_indices (dofs_per_cell);
+
+ // Start of input vector inside a
+ // block vector.
+ const unsigned int start = component_start[selected];
+
+ typename MGDoFHandler<dim>::active_cell_iterator
+ level_cell = mg_dof_handler.begin_active();
+ const typename MGDoFHandler<dim>::active_cell_iterator
+ endc = mg_dof_handler.end();
+
+ // traverse all cells and copy the
+ // data appropriately to the output
+ // vector
+
+ // Is the level monotonuosly increasing?
+
+ for (; level_cell != endc; ++level_cell)
+ {
+ DoFObjectAccessor<dim, dim>& global_cell = *level_cell;
+ const unsigned int level = level_cell->level();
+
+ // Start of treated component
+ // for this level
+ const unsigned int level_start = mg_component_start[level][selected];
+
+ // 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);
+
+//TODO: Probably wrong for continuous elements
+
+ // copy level-wise data to
+ // global vector
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ if (fe.system_to_component_index(i).first == selected)
+ dst(global_dof_indices[i]-start)
+ += src[level](level_dof_indices[i]-level_start);
+ };
+
+ // clear constrained nodes
+//TODO: constraints->set_zero(dst);
+}
+
+
+#endif
--- /dev/null
+//----------------------------------------------------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 1998, 1999, 2000, 2001 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//----------------------------------------------------------------------
+#ifndef __deal2__mg_transfer_h
+#define __deal2__mg_transfer_h
+
+
+#include <lac/sparsity_pattern.h>
+#include <lac/block_vector.h>
+#include <lac/block_sparsity_pattern.h>
+#include <multigrid/mg_base.h>
+
+template <int dim> class MGDoFHandler;
+
+
+/**
+ * Implementation of the @p{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 Subscriptor // MGTransferBase<Vector<double> >
+{
+ public:
+ /**
+ * Destructor.
+ */
+ virtual ~MGTransferPrebuilt ();
+
+ /**
+ * Actually build the prolongation
+ * matrices for each level.
+ */
+ template <int dim>
+ void build_matrices (const MGDoFHandler<dim> &mg_dof);
+
+ /**
+ * Prolongate a vector from level
+ * @p{to_level-1} to level
+ * @p{to_level}. The previous
+ * content of @p{dst} is
+ * overwritten.
+ *
+ * @p{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 @p{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<double> &dst,
+ const Vector<double> &src) const;
+
+ /**
+ * Restrict a vector from level
+ * @p{from_level} to level
+ * @p{from_level-1} and add this
+ * restriction to
+ * @p{dst}. Obviously, if the
+ * refined region on level
+ * @p{from_level} is smaller than
+ * that on level @p{from_level-1},
+ * some degrees of freedom in
+ * @p{dst} are not covered and will
+ * not be altered. For the other
+ * degress of freedom, the result
+ * of the restriction is added.
+ *
+ * @p{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 @p{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<double> &dst,
+ const Vector<double> &src) const;
+
+ /**
+ * Transfer from a vector on the
+ * global grid to vectors defined
+ * on each of the levels
+ * separately, i.a. an @p{MGVector}.
+ */
+ template<int dim, class InVector>
+ void
+ copy_to_mg (const MGDoFHandler<dim>& mg_dof,
+ MGLevelObject<Vector<double> > &dst,
+ const InVector &src) const;
+
+ /**
+ * Transfer from multi-level vector to
+ * normal vector.
+ *
+ * Copies data from active
+ * portions of an MGVector into
+ * the respective positions of a
+ * @p{Vector<double>}. In order to
+ * keep the result consistent,
+ * constrained degrees of freedom
+ * are set to zero.
+ */
+ template<int dim, class OutVector>
+ void
+ copy_from_mg (const MGDoFHandler<dim>& mg_dof,
+ OutVector &dst,
+ const MGLevelObject<Vector<double> > &src) const;
+
+ /**
+ * Add a multi-level vector to a
+ * normal vector.
+ *
+ * Works as the previous
+ * function, but probably not for
+ * continuous elements.
+ */
+ template<int dim, class OutVector>
+ void
+ copy_from_mg_add (const MGDoFHandler<dim>& mg_dof,
+ OutVector &dst,
+ const MGLevelObject<Vector<double> > &src) const;
+
+ private:
+
+ /**
+ * Selected component.
+ */
+ unsigned int selected;
+
+ /**
+ * Sizes of the multi-level vectors.
+ */
+ std::vector<unsigned int> sizes;
+
+
+ std::vector<SparsityPattern> 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.
+ */
+ std::vector<SparseMatrix<double> > prolongation_matrices;
+};
+
+
+/**
+ * Implementation of matrix generation for @ref{MGTransferBlock} and @p{MGTransferSelected}.
+ *
+ * @author Guido Kanschat, 2001
+ */
+class MGTransferBlockBase
+{
+ protected:
+ /**
+ * Actually build the prolongation
+ * matrices for each level.
+ *
+ * If a field selected is given,
+ * only matrices for these
+ * components are to be built. By
+ * default, all matrices are
+ * built.
+ */
+ template <int dim>
+ void build_matrices (const MGDoFHandler<dim> &mg_dof,
+ std::vector<bool> selected);
+
+ /**
+ * Flag of selected components.
+ */
+ std::vector<bool> selected;
+
+ /**
+ * Sizes of the multi-level vectors.
+ */
+ std::vector<std::vector<unsigned int> > sizes;
+
+ /**
+ * Start index of each component.
+ */
+ std::vector<unsigned int> component_start;
+
+ /**
+ * Start index of eah component on
+ * all levels.
+ */
+ std::vector<std::vector<unsigned int> > mg_component_start;
+
+ private:
+
+ std::vector<BlockSparsityPattern> prolongation_sparsities;
+
+ protected:
+
+ /**
+ * 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.
+ */
+ std::vector<BlockSparseMatrix<double> > prolongation_matrices;
+};
+
+
+
+/**
+ * Implementation of the @p{MGTransferBase} interface for block
+ * matrices and block vectors. In addition to the functionality of
+ * @ref{MGTransferPrebuilt}, the operation may be restricted to
+ * certain blocks of the vector.
+ *
+ * If the restricted mode is chosen, block vectors used in the
+ * transfer routines may only have as many components as there are
+ * @p{true}s in the selected-field.
+ *
+ * @author Guido Kanschat, 2001
+ */
+class MGTransferBlock : public Subscriptor, // MGTransferBase<BlockVector<double> >,
+ private MGTransferBlockBase
+{
+ public:
+ /**
+ * Destructor.
+ */
+ virtual ~MGTransferBlock ();
+
+ /**
+ * Build the prolongation
+ * matrices for each level.
+ *
+ * If a field selected is given,
+ * only matrices for these
+ * components are to be built. By
+ * default, all matrices are
+ * built.
+ *
+ * This function is a front-end
+ * for the same function in
+ * @ref{MGTransferBlockBase}.
+ */
+ template <int dim>
+ void build_matrices (const MGDoFHandler<dim> &mg_dof,
+ std::vector<bool> selected = std::vector<bool>());
+
+ /**
+ * Prolongate a vector from level
+ * @p{to_level-1} to level
+ * @p{to_level}. The previous
+ * content of @p{dst} is
+ * overwritten.
+ *
+ * @p{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 @p{src}
+ * shall have as many elements as there
+ * are degrees of freedom on the finer
+ * level.
+ */
+ virtual void prolongate (const unsigned int to_level,
+ BlockVector<double> &dst,
+ const BlockVector<double> &src) const;
+
+ /**
+ * Restrict a vector from level
+ * @p{from_level} to level
+ * @p{from_level-1} and add this
+ * restriction to
+ * @p{dst}. Obviously, if the
+ * refined region on level
+ * @p{from_level} is smaller than
+ * that on level @p{from_level-1},
+ * some degrees of freedom in
+ * @p{dst} are not covered and will
+ * not be altered. For the other
+ * degress of freedom, the result
+ * of the restriction is added.
+ *
+ * @p{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 @p{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,
+ BlockVector<double> &dst,
+ const BlockVector<double> &src) const;
+
+ protected:
+};
+
+
+
+/**
+ * Implementation of the @p{MGTransferBase} interface for block
+ * matrices and simple vectors. This class uses @ref{MGTransferBlock}
+ * selecting a single component. The transfer operators themselves are
+ * implemented for simple vectors again.
+ *
+ * @author Guido Kanschat, 2001
+ */
+class MGTransferSelect : public Subscriptor, // MGTransferBase<Vector<double> >,
+ private MGTransferBlockBase
+{
+ public:
+ /**
+ * Destructor.
+ */
+ virtual ~MGTransferSelect ();
+
+ /**
+ * Actually build the prolongation
+ * matrices for each level.
+ *
+ * Select the component you want
+ * to apply multigrid to.
+ *
+ * This function is a front-end
+ * for the same function in
+ * @ref{MGTransferBlockBase}.
+ */
+ template <int dim>
+ void build_matrices (const MGDoFHandler<dim> &mg_dof,
+ unsigned int selected);
+
+ /**
+ * Prolongate a vector from level
+ * @p{to_level-1} to level
+ * @p{to_level}. The previous
+ * content of @p{dst} is
+ * overwritten.
+ *
+ * @p{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 @p{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<double> &dst,
+ const Vector<double> &src) const;
+
+ /**
+ * Restrict a vector from level
+ * @p{from_level} to level
+ * @p{from_level-1} and add this
+ * restriction to
+ * @p{dst}. Obviously, if the
+ * refined region on level
+ * @p{from_level} is smaller than
+ * that on level @p{from_level-1},
+ * some degrees of freedom in
+ * @p{dst} are not covered and will
+ * not be altered. For the other
+ * degress of freedom, the result
+ * of the restriction is added.
+ *
+ * @p{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 @p{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<double> &dst,
+ const Vector<double> &src) const;
+
+ /**
+ * Transfer from a vector on the
+ * global grid to vectors defined
+ * on each of the levels
+ * separately, i.a. an @p{MGVector}.
+ */
+ template<int dim, class InVector>
+ void
+ copy_to_mg (const MGDoFHandler<dim>& mg_dof,
+ MGLevelObject<Vector<double> > &dst,
+ const InVector &src) const;
+
+ /**
+ * Transfer from multi-level vector to
+ * normal vector.
+ *
+ * Copies data from active
+ * portions of an MGVector into
+ * the respective positions of a
+ * @p{Vector<double>}. In order to
+ * keep the result consistent,
+ * constrained degrees of freedom
+ * are set to zero.
+ */
+ template<int dim, class OutVector>
+ void
+ copy_from_mg (const MGDoFHandler<dim>& mg_dof,
+ OutVector &dst,
+ const MGLevelObject<Vector<double> > &src) const;
+
+ /**
+ * Add a multi-level vector to a
+ * normal vector.
+ *
+ * Works as the previous
+ * function, but probably not for
+ * continuous elements.
+ */
+ template<int dim, class OutVector>
+ void
+ copy_from_mg_add (const MGDoFHandler<dim>& mg_dof,
+ OutVector &dst,
+ const MGLevelObject<Vector<double> > &src) const;
+
+ private:
+
+ /**
+ * Selected component.
+ */
+ unsigned int selected;
+};
+
+
+
+#endif
/**
- * Implementation of multigrid with matrices.
- * While @p{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.
+ * Implementation of the multigrid method.
+ *
+ * The function actually performing a multi-level cycle,
+ * @p{level_mgstep}, as well as the function @p{vcycle}, calling it,
+ * require several helper classes handed over as template parameters.
+ * These classes have to meet the following requirements:
+ *
+ * @p{MATRIX} is a matrix as specified for LAC-solvers, that is, it
+ * has a function
+ * \begin{verbatim}
+ * void MATRIX::vmult(VECTOR& x, const VECTOR& b)
+ * \end{verbatim}
+ * performing the matrix vector product @p{x=Ab}.
+ *
+ * @p{SMOOTH} is a class with a function
+ * \begin{verbatim}
+ * void SMOOTH::smooth (unsigned int level, VECTOR& x, const VECTOR& b)
+ * \end{verbatim}
+ * modifying the vector @p{x} such that the residual @{b-Ax}
+ * on level @p{l} is smoothened. Refer to @ref{MGSmootherRelaxation}
+ * for an example.
+ *
+ * @p{COARSE} has an
+ * \begin{verbatim}
+ * void COARSE::operator() (VECTOR& x, const VECTOR& b)
+ * \end{verbatim}
+ * returning the solution to the system @p{Ax=b} on the coarsest level
+ * in @p{x}.
*
- * @author Wolfgang Bangerth, Guido Kanschat, 1999
+ * @author Guido Kanschat, 1999, 2001
*/
-template <int dim>
-class Multigrid : public MGBase
+template <class VECTOR>
+class Multigrid : public Subscriptor
{
public:
+ typedef VECTOR vector_type;
+ typedef const VECTOR const_vector_type;
+
/**
- * 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,
- * their respective sparsity
- * patterns, and an object
- * mediating the transfer of
- * solutions between different
- * levels.
+ * Constructor. The
+ * @p{MGDoFHandler} is used to
+ * determine the highest possible
+ * level. @p{transfer} is an
+ * object performing prolongation
+ * and restriction.
+ *
+ * The V-cycle will start on
+ * level @p{maxlevel} and goes
+ * down to level @p{minlevel},
+ * where the coarse grid solver
+ * will be used.
*
* This function already
* initializes the vectors which
* therefore create objects of
* this type as late as possible.
*/
-//TODO:[GK] meaning of minlevel, maxlevel?
- Multigrid(const MGDoFHandler<dim> &mg_dof_handler,
- const ConstraintMatrix &constraints,
- const MGLevelObject<SparsityPattern> &level_sparsities,
- const MGLevelObject<SparseMatrix<double> > &level_matrices,
- const MGTransferBase &transfer,
- const unsigned int minlevel = 0,
- const unsigned int maxlevel = 1000000);
+ template <int dim>
+ Multigrid(const MGDoFHandler<dim>& mg_dof_handler,
+ 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 @p{MGVector}.
+ * Execute one step of the
+ * v-cycle algorithm. This
+ * function assumes, that the
+ * vector @p{defect} is properly
+ * filled with the residual in
+ * the outer defect correction
+ * scheme (usually performed by
+ * @ref{PreconditionMG}). After
+ * execution of @p{vcycle()}, the
+ * result is in the vector
+ * @p{solution}. See
+ * @p{copy_*_mg} in class
+ * @p{MGTools} if you want to use
+ * these vectors yourself.
+ *
+ * The actual work for this
+ * function is done in
+ * @p{level_mgstep}.
*/
- template<typename number>
- void copy_to_mg (const Vector<number> &src);
+ template <class MATRIX, class TRANSFER, class SMOOTHER, class COARSE>
+ void vcycle(const MGLevelObject<MATRIX>& matrix,
+ const TRANSFER& transfer,
+ const SMOOTHER& pre_smooth,
+ const SMOOTHER& post_smooth,
+ const COARSE& cgs);
/**
- * Transfer from multi-level vector to
- * normal vector.
- *
- * Copies data from active
- * portions of an MGVector into
- * the respective positions of a
- * @p{Vector<double>}. In order to
- * keep the result consistent,
- * constrained degrees of freedom
- * are set to zero.
+ * Print a level vector using
+ * @ref{DoFHandler}.
*/
- template<typename number>
- void copy_from_mg (Vector<number> &dst) const;
+/* virtual void print_vector (const unsigned int level, */
+/* const VECTOR& v, */
+/* const char* name) const; */
+
+ private:
+
+ /**
+ * The actual v-cycle multigrid method.
+ * @p{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 @p{level-1},
+ * unless we are on @p{minlevel}
+ * where instead of recursing
+ * deeper, the coarse grid solver
+ * is used to solve the matrix of
+ * this level.
+ */
+ template <class MATRIX, class TRANSFER, class SMOOTHER, class COARSE>
+ void level_mgstep (const unsigned int level,
+ const MGLevelObject<MATRIX>& matrix,
+ const TRANSFER& transfer,
+ const SMOOTHER& pre_smooth,
+ const SMOOTHER& post_smooth,
+ const COARSE& cgs);
/**
- * Negative @p{vmult} on a level.
- * @see MGBase.
+ * Level for coarse grid solution.
*/
-//TODO:[GK] what does this function do?
- virtual void level_vmult (const unsigned int level,
- Vector<double> &dest,
- const Vector<double> &src,
- const Vector<double> &rhs);
+ unsigned int minlevel;
/**
- * Print a level vector using
- * @ref{DoFHandler}.
+ * Highest level of cells.
*/
- virtual void print_vector (const unsigned int level,
- const Vector<double>& v,
- const char* name) const;
+ unsigned int maxlevel;
/**
- * Read-only access to level matrices.
+ * Auxiliary vector. Contains the
+ * defect to be corrected before
+ * the multigrid step.
*/
- const SparseMatrix<double>& get_matrix (const unsigned int level) const;
+ MGLevelObject<VECTOR> defect;
- private:
/**
- * Associated @p{MGDoFHandler}.
+ * Auxiliary vector. Contains the
+ * updated solution after the
+ * multigrid step.
*/
- SmartPointer<const MGDoFHandler<dim> > mg_dof_handler;
-
-
-/**
- * Sparsity patterns for each level.
- */
- SmartPointer<const MGLevelObject<SparsityPattern> > level_sparsities;
+ MGLevelObject<VECTOR> solution;
/**
- * Matrices for each level. The
- * matrices are prepared by the
- * constructor of @p{Multigrid} and
- * can be accessed for
- * assembling.
+ * Auxiliary vector.
*/
- SmartPointer<const MGLevelObject<SparseMatrix<double> > > level_matrices;
+ MGLevelObject<VECTOR> t;
+
+
+ private:
+
/**
- * Pointer to the object holding
- * constraints.
+ * Exception.
*/
- SmartPointer<const ConstraintMatrix> constraints;
+ DeclException2(ExcSwitchedLevels, int, int,
+ << "minlevel and maxlevel switched, should be: "
+ << arg1 << "<=" << arg2);
+ template<int dim, class MATRIX, class VECTOR2, class TRANSFER, class SMOOTHER, class COARSE> friend class PreconditionMG;
};
/**
- * Implementation of the @p{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.
+ * 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 @p{MG} is required to inherit @p{MGBase}.
+ * Furthermore, it needs functions @p{void copy_to_mg(const VECTOR&)}
+ * to store @p{src} in the right hand side of the multi-level method and
+ * @p{void copy_from_mg(VECTOR&)} to store the result of the v-cycle in @p{dst}.
*
- * @author Wolfgang Bangerth, Guido Kanschat, 1999
+ * @author Guido Kanschat, 1999
*/
-class MGTransferPrebuilt : public MGTransferBase
+template<int dim, class MATRIX, class VECTOR, class TRANSFER, class SMOOTHER, class COARSE>
+class PreconditionMG : public Subscriptor
{
public:
/**
- * Destructor.
+ * Constructor.
+ * Arguments are the multigrid object,
+ * pre-smoother, post-smoother and
+ * coarse grid solver.
*/
- virtual ~MGTransferPrebuilt ();
+ PreconditionMG(const MGDoFHandler<dim>& mg_dof,
+ Multigrid<VECTOR>& mg,
+ const MGLevelObject<MATRIX>& matrix,
+ const TRANSFER& transfer,
+ const SMOOTHER& pre,
+ const SMOOTHER& post,
+ const COARSE& coarse);
/**
- * Actually build the prolongation
- * matrices for each level.
+ * Preconditioning operator.
+ * Calls the @p{vcycle} function
+ * of the @p{MG} object passed to
+ * the constructor.
+ *
+ * This is the operator used by
+ * LAC iterative solvers.
*/
- template <int dim>
- void build_matrices (const MGDoFHandler<dim> &mg_dof);
-
+ void vmult (VECTOR &dst,
+ const VECTOR &src) const;
+
+ /**
+ * Preconditioning operator.
+ * Calls the @p{vcycle} function
+ * of the @p{MG} object passed to
+ * the constructor.
+ */
+ void vmult_add (VECTOR &dst,
+ const VECTOR &src) const;
+
/**
- * Prolongate a vector from level
- * @p{to_level-1} to level
- * @p{to_level}. The previous
- * content of @p{dst} is
- * overwritten.
+ * Tranposed preconditioning operator.
*
- * @p{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 @p{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<double> &dst,
- const Vector<double> &src) const;
-
- /**
- * Restrict a vector from level
- * @p{from_level} to level
- * @p{from_level-1} and add this
- * restriction to
- * @p{dst}. Obviously, if the
- * refined region on level
- * @p{from_level} is smaller than
- * that on level @p{from_level-1},
- * some degrees of freedom in
- * @p{dst} are not covered and will
- * not be altered. For the other
- * degress of freedom, the result
- * of the restriction is added.
+ * Not implemented, but the
+ * definition may be needed.
+ */
+ void Tvmult (VECTOR &dst,
+ const VECTOR &src) const;
+
+ /**
+ * Tranposed preconditioning operator.
*
- * @p{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 @p{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<double> &dst,
- const Vector<double> &src) const;
-
+ * Not implemented, but the
+ * definition may be needed.
+ */
+ void Tvmult_add (VECTOR &dst,
+ const VECTOR &src) const;
+
private:
+ /**
+ * Associated @p{MGDoFHandler}.
+ */
+ SmartPointer<const MGDoFHandler<dim> > mg_dof_handler;
+
+ /**
+ * The multigrid object.
+ */
+ SmartPointer<Multigrid<VECTOR> > multigrid;
+
+ /**
+ * Matrices for each level. The
+ * matrices are prepared by the
+ * constructor of @p{Multigrid} and
+ * can be accessed for
+ * assembling.
+ */
+ SmartPointer<const MGLevelObject<MATRIX> > level_matrices;
- std::vector<SparsityPattern> prolongation_sparsities;
+ /**
+ * Object for grid tranfer.
+ */
+ SmartPointer<const TRANSFER> transfer;
+
+ /**
+ * The pre-smoothing object.
+ */
+ SmartPointer<const SMOOTHER> pre;
/**
- * 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 post-smoothing object.
+ */
+ SmartPointer<const SMOOTHER> post;
+
+ /**
+ * The coarse grid solver.
*/
- std::vector<SparseMatrix<double> > prolongation_matrices;
+ SmartPointer<const COARSE> coarse;
};
+
+
/* --------------------------- inline functions --------------------- */
-template<int dim>
-inline
-const SparseMatrix<double>&
-Multigrid<dim>::get_matrix (const unsigned int level) const
+template <class VECTOR>
+template <int dim>
+Multigrid<VECTOR>::Multigrid (const MGDoFHandler<dim>& mg_dof_handler,
+ const unsigned int min_level,
+ const unsigned int max_level)
+ :
+ minlevel(min_level),
+ maxlevel(std::min(mg_dof_handler.get_tria().n_levels()-1,
+ max_level)),
+ defect(minlevel,maxlevel),
+ solution(minlevel,maxlevel),
+ t(minlevel,maxlevel)
+{
+};
+
+
+template <class VECTOR>
+template <class MATRIX, class TRANSFER, class SMOOTHER, class COARSE>
+void
+Multigrid<VECTOR>::level_mgstep(const unsigned int level,
+ const MGLevelObject<MATRIX>& matrix,
+ const TRANSFER& transfer,
+ const SMOOTHER &pre_smooth,
+ const SMOOTHER &post_smooth,
+ const COARSE &coarse_grid_solver)
{
- Assert((level>=minlevel) && (level<maxlevel),
- ExcIndexRange(level, minlevel, maxlevel));
+#ifdef MG_DEBUG
+ char *name = new char[100];
+ sprintf(name, "MG%d-defect",level);
+ print_vector(level, defect[level], name);
+#endif
+
+ solution[level] = 0.;
+
+ if (level == minlevel)
+ {
+ coarse_grid_solver(level, solution[level], defect[level]);
+#ifdef MG_DEBUG
+ sprintf(name, "MG%d-solution",level);
+ print_vector(level, solution[level], name);
+#endif
+ return;
+ }
+
+ // smoothing of the residual by modifying s
+ pre_smooth.smooth(level, solution[level], defect[level]);
+
+#ifdef MG_DEBUG
+ sprintf(name, "MG%d-pre",level);
+ print_vector(level, solution[level], name);
+#endif
+
+ // t = -A*solution[level]
+ matrix[level].vmult(t[level], solution[level]);
+ t[level].scale(-1);
+
+ // make t rhs of lower level
+ // The non-refined parts of the
+ // coarse-level defect already contain
+ // the global defect, the refined parts
+ // its restriction.
+ for (unsigned int l = level;l>minlevel;--l)
+ {
+ t[l-1] = 0.;
+ transfer.restrict_and_add (l, t[l-1], t[l]);
+ defect[l-1] += t[l-1];
+ }
+
+ // add additional DG contribution
+// edge_vmult(level, defect[level-1], defect[level]);
+
+ // do recursion
+ level_mgstep(level-1, matrix, transfer, 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[level] = 0.;
+
+ // do coarse grid correction
+
+ transfer.prolongate(level, t[level], solution[level-1]);
+
+#ifdef MG_DEBUG
+ sprintf(name, "MG%d-cgc",level);
+ print_vector(level, t, name);
+#endif
+
+ solution[level] += t[level];
- return (*level_matrices)[level];
+ // post-smoothing
+ post_smooth.smooth(level, solution[level], defect[level]);
+
+#ifdef MG_DEBUG
+ sprintf(name, "MG%d-post",level);
+ print_vector(level, solution[level], name);
+
+ delete[] name;
+#endif
+}
+
+
+template <class VECTOR>
+template <class MATRIX, class TRANSFER, class SMOOTHER, class COARSE>
+void
+Multigrid<VECTOR>::vcycle(const MGLevelObject<MATRIX>& matrix,
+ const TRANSFER& transfer,
+ const SMOOTHER &pre_smooth,
+ const SMOOTHER &post_smooth,
+ const COARSE &coarse_grid_solver)
+{
+ // The defect vector has been
+ // initialized by copy_to_mg. Now
+ // adjust the other vectors.
+ solution.resize(minlevel, maxlevel);
+ t.resize(minlevel, maxlevel);
+
+ for (unsigned int level=minlevel; level<=maxlevel;++level)
+ {
+ solution[level].reinit(defect[level]);
+ t[level].reinit(defect[level]);
+ }
+
+ level_mgstep (maxlevel, matrix, transfer, pre_smooth, post_smooth, coarse_grid_solver);
+// abort ();
+}
+
+
+/* ------------------------------------------------------------------- */
+
+
+template<int dim, class MATRIX, class VECTOR, class TRANSFER, class SMOOTHER, class COARSE>
+PreconditionMG<dim, MATRIX, VECTOR, TRANSFER, SMOOTHER, COARSE>
+::PreconditionMG(const MGDoFHandler<dim>& mg_dof_handler,
+ Multigrid<VECTOR>& mg,
+ const MGLevelObject<MATRIX>& matrix,
+ const TRANSFER& transfer,
+ const SMOOTHER& pre,
+ const SMOOTHER& post,
+ const COARSE& coarse)
+ :
+ mg_dof_handler(&mg_dof_handler),
+ multigrid(&mg),
+ level_matrices(&matrix),
+ transfer(&transfer),
+ pre(&pre),
+ post(&post),
+ coarse(&coarse)
+{}
+
+
+template<int dim, class MATRIX, class VECTOR, class TRANSFER, class SMOOTHER, class COARSE>
+void
+PreconditionMG<dim, MATRIX, VECTOR, TRANSFER, SMOOTHER, COARSE>::vmult (
+ VECTOR& dst,
+ const VECTOR& src) const
+{
+ transfer->copy_to_mg(*mg_dof_handler,
+ multigrid->defect,
+ src);
+
+ multigrid->vcycle(*level_matrices,
+ *transfer,
+ *pre, *post, *coarse);
+
+ transfer->copy_from_mg(*mg_dof_handler,
+ dst,
+ multigrid->solution);
+}
+
+
+template<int dim, class MATRIX, class VECTOR, class TRANSFER, class SMOOTHER, class COARSE>
+void
+PreconditionMG<dim, MATRIX, VECTOR, TRANSFER, SMOOTHER, COARSE>::vmult_add (
+ VECTOR& dst,
+ const VECTOR& src) const
+{
+ transfer->copy_to_mg(*mg_dof_handler,
+ multigrid->defect,
+ src);
+
+ multigrid->vcycle(*level_matrices,
+ *transfer,
+ *pre, *post, *coarse);
+
+ transfer->copy_from_mg_add(*mg_dof_handler,
+ dst,
+ multigrid->solution);
+}
+
+
+template<int dim, class MATRIX, class VECTOR, class TRANSFER, class SMOOTHER, class COARSE>
+void
+PreconditionMG<dim, MATRIX, VECTOR, TRANSFER, SMOOTHER, COARSE>::Tvmult (
+ VECTOR&,
+ const VECTOR&) const
+{
+ Assert(false, ExcNotImplemented());
+}
+
+
+template<int dim, class MATRIX, class VECTOR, class TRANSFER, class SMOOTHER, class COARSE>
+void
+PreconditionMG<dim, MATRIX, VECTOR, TRANSFER, SMOOTHER, COARSE>::Tvmult_add (
+ VECTOR&,
+ const VECTOR&) const
+{
+ Assert(false, ExcNotImplemented());
}
// make the connection graph
SparsityPattern sparsity (dof_handler.n_dofs(level),
dof_handler.max_couplings_between_dofs());
- MGDoFTools::make_sparsity_pattern (dof_handler, sparsity, level);
+ MGTools::make_sparsity_pattern (dof_handler, sparsity, level);
const unsigned int n_dofs = sparsity.n_rows();
// store the new dof numbers; invalid_dof_index means
+template <int dim>
+void DoFRenumbering::component_wise (MGDoFHandler<dim>& dof_handler,
+ unsigned int level,
+ const std::vector<unsigned int> &component_order_arg)
+{
+ const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
+
+ // do nothing if the FE has only
+ // one component
+ if (dof_handler.get_fe().n_components() == 1)
+ return;
+
+ std::vector<unsigned int> component_order (component_order_arg);
+ if (component_order.size() == 0)
+ for (unsigned int i=0; i<dof_handler.get_fe().n_components(); ++i)
+ component_order.push_back (i);
+
+ // check whether the component list has
+ // the right length and contains all
+ // component numbers
+ Assert (component_order.size() == dof_handler.get_fe().n_components(),
+ ExcInvalidComponentOrder());
+ for (unsigned int i=0; i<dof_handler.get_fe().n_components(); ++i)
+ Assert (std::find (component_order.begin(), component_order.end(), i)
+ != component_order.end(),
+ ExcInvalidComponentOrder ());
+
+ // vector to hold the dof indices on
+ // the cell we visit at a time
+ std::vector<unsigned int> local_dof_indices(dofs_per_cell);
+ // prebuilt list to which component
+ // a given dof on a cell belongs
+ std::vector<unsigned int> component_list (dofs_per_cell);
+ for (unsigned int i=0; i<component_list.size(); ++i)
+ component_list[i] = dof_handler.get_fe().system_to_component_index(i).first;
+ // set up a map where for each component
+ // the respective degrees of freedom are
+ // collected.
+ //
+ // note that this map is sorted by component
+ // but that within each component it is NOT
+ // sorted by dof index. note also that some
+ // dof indices are entered multiply, so we
+ // will have to take care of that
+ std::vector<std::vector<unsigned int> > component_to_dof_map (dof_handler.get_fe().n_components());
+ for (typename MGDoFHandler<dim>::cell_iterator cell=dof_handler.begin(level);
+ cell!=dof_handler.end(level); ++cell)
+ {
+ // on each cell: get dof indices
+ // and insert them into the global
+ // list using their component
+ cell->get_mg_dof_indices (local_dof_indices);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ component_to_dof_map[component_list[i]].push_back (local_dof_indices[i]);
+ };
+
+ // now we've got all indices sorted into
+ // buckets labelled with their component
+ // number. we've only got to traverse this
+ // list and assign the new indices
+ //
+ // however, we first want to sort the
+ // indices entered into the buckets to
+ // preserve the order within each component
+ // and during this also remove duplicate
+ // entries
+ for (unsigned int component=0; component<dof_handler.get_fe().n_components(); ++component)
+ {
+ std::sort (component_to_dof_map[component].begin(),
+ component_to_dof_map[component].end());
+ component_to_dof_map[component].erase (std::unique (component_to_dof_map[component].begin(),
+ component_to_dof_map[component].end()),
+ component_to_dof_map[component].end());
+ };
+
+ unsigned int next_free_index = 0;
+ std::vector<unsigned int> new_indices (dof_handler.n_dofs(level),
+ DoFHandler<dim>::invalid_dof_index);
+ for (unsigned int component=0; component<dof_handler.get_fe().n_components(); ++component)
+ {
+ const typename std::vector<unsigned int>::const_iterator
+ begin_of_component = component_to_dof_map[component].begin(),
+ end_of_component = component_to_dof_map[component].end();
+
+ for (typename std::vector<unsigned int>::const_iterator dof_index = begin_of_component;
+ dof_index != end_of_component; ++dof_index)
+ new_indices[*dof_index] = next_free_index++;
+ };
+
+// Assert (next_free_index == dof_handler.n_dofs(level),
+// ExcInternalError());
+
+ dof_handler.renumber_dofs (level, new_indices);
+};
+
+
+
template <int dim>
void
DoFRenumbering::sort_selected_dofs_back (DoFHandler<dim> &dof_handler,
void DoFRenumbering::component_wise (DoFHandler<deal_II_dimension>&,
const std::vector<unsigned int>&);
+template
+void DoFRenumbering::component_wise (MGDoFHandler<deal_II_dimension>&,
+ unsigned int,
+ const std::vector<unsigned int>&);
+
template
void DoFRenumbering::cell_wise_dg (DoFHandler<deal_II_dimension>&,
#include <multigrid/mg_base.h>
#include <multigrid/mg_smoother.h>
+#include <lac/sparse_matrix.h>
+#include <lac/block_sparse_matrix.h>
#include <iostream>
#include <cmath>
-MGBase::~MGBase ()
-{}
-
-
-MGBase::MGBase(const MGTransferBase &transfer,
- const unsigned minlevel,
- const unsigned maxlevel)
- :
- maxlevel(maxlevel),
- minlevel(minlevel),
- defect(minlevel,maxlevel),
- solution(minlevel,maxlevel),
- t(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);
-// abort ();
-}
-
-
-void
-MGBase::level_mgstep(const unsigned int level,
- const MGSmootherBase &pre_smooth,
- const MGSmootherBase &post_smooth,
- const MGCoarseGridSolver &coarse_grid_solver)
-{
-#ifdef MG_DEBUG
- char *name = new char[100];
- sprintf(name, "MG%d-defect",level);
- print_vector(level, defect[level], name);
-#endif
-
- solution[level] = 0.;
-
- if (level == minlevel)
- {
- coarse_grid_solver(level, solution[level], defect[level]);
-#ifdef MG_DEBUG
- sprintf(name, "MG%d-solution",level);
- print_vector(level, solution[level], name);
-#endif
- return;
- }
-
- // smoothing of the residual by modifying s
- pre_smooth.smooth(level, solution[level], defect[level]);
-
-#ifdef MG_DEBUG
- sprintf(name, "MG%d-pre",level);
- print_vector(level, solution[level], name);
-#endif
-
- t[level] = 0.;
-
- // t = -A*solution[level]
- level_vmult(level, t[level], solution[level], defect[level]);
-
- // make t rhs of lower level
- // The non-refined parts of the
- // coarse-level defect already contain
- // the global defect, the refined parts
- // its restriction.
- for (unsigned int l = level;l>minlevel;--l)
- {
- t[l-1] = 0.;
- transfer->restrict_and_add (l, t[l-1], t[l]);
- defect[l-1] += t[l-1];
- }
-
- // add additional DG contribution
-// edge_vmult(level, defect[level-1], defect[level]);
-
- // 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[level] = 0.;
-
- // do coarse grid correction
-
- transfer->prolongate(level, t[level], solution[level-1]);
-
-#ifdef MG_DEBUG
- sprintf(name, "MG%d-cgc",level);
- print_vector(level, t, name);
-#endif
-
- solution[level] += t[level];
-
- // post-smoothing
- post_smooth.smooth(level, solution[level], defect[level]);
-
-#ifdef MG_DEBUG
- sprintf(name, "MG%d-post",level);
- print_vector(level, solution[level], name);
-
- delete[] name;
-#endif
-}
-
-
-void
-MGBase::edge_vmult (const unsigned int,
- Vector<double>&,
- const Vector<double>&)
-{}
-
//////////////////////////////////////////////////////////////////////
-MGCoarseGridSolver::~MGCoarseGridSolver()
-{};
-
-
-//////////////////////////////////////////////////////////////////////
+// template <class VECTOR>
+// MGTransferBase<VECTOR>::~MGTransferBase()
+// {};
-MGTransferBase::~MGTransferBase()
-{};
+// // Explicit instantiations
+// template class MGTransferBase<Vector<double> >;
+// template class MGTransferBase<Vector<float> >;
+// template class MGTransferBase<BlockVector<double> >;
+// template class MGTransferBase<BlockVector<float> >;
//---------------------------- mg_dof_tools.cc ---------------------------
+#include <base/multithread_info.h>
+#include <base/thread_management.h>
#include <lac/sparsity_pattern.h>
-#include <multigrid/mg_dof_handler.h>
-#include <multigrid/mg_dof_accessor.h>
+#include <lac/block_sparsity_pattern.h>
+#include <lac/compressed_sparsity_pattern.h>
+#include <lac/sparsity_pattern.h>
+#include <lac/block_vector.h>
#include <grid/tria.h>
#include <grid/tria_iterator.h>
+#include <multigrid/mg_dof_handler.h>
+#include <multigrid/mg_dof_accessor.h>
#include <multigrid/mg_dof_tools.h>
+#include <multigrid/mg_base.h>
+#include <dofs/dof_tools.h>
#include <fe/fe.h>
+//#include <multigrid/mg_tools.templates.h>
-template <int dim>
-void MGDoFTools::make_sparsity_pattern (const MGDoFHandler<dim> &dof,
+#include <algorithm>
+#include <numeric>
+
+template <int dim, class SparsityPattern>
+void MGTools::make_sparsity_pattern (const MGDoFHandler<dim> &dof,
SparsityPattern &sparsity,
const unsigned int level)
{
-template<int dim>
+template<int dim, class SparsityPattern>
void
-MGDoFTools::make_flux_sparsity_pattern (const MGDoFHandler<dim> &dof,
+MGTools::make_flux_sparsity_pattern (const MGDoFHandler<dim> &dof,
SparsityPattern &sparsity,
const unsigned int level)
{
-template<int dim>
+template<int dim, class SparsityPattern>
void
-MGDoFTools::make_flux_sparsity_pattern_edge (const MGDoFHandler<dim> &dof,
+MGTools::make_flux_sparsity_pattern_edge (const MGDoFHandler<dim> &dof,
SparsityPattern &sparsity,
const unsigned int level)
{
-template<int dim>
+//TODO: Replace FullMatrix by vector2d<bool>
+template<int dim, class SparsityPattern>
void
-MGDoFTools::make_flux_sparsity_pattern (const MGDoFHandler<dim> &dof,
+MGTools::make_flux_sparsity_pattern (const MGDoFHandler<dim> &dof,
SparsityPattern &sparsity,
const unsigned int level,
const FullMatrix<double>& int_mask,
+template <int dim>
+void
+MGTools::count_dofs_per_component (const MGDoFHandler<dim>& dof_handler,
+ std::vector<std::vector<unsigned int> >& result)
+{
+ const unsigned int nlevels = dof_handler.get_tria().n_levels();
+
+ Assert (result.size() == nlevels,
+ ExcDimensionMismatch(result.size(), nlevels));
-// explicit instantiations
-template void MGDoFTools::make_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
- SparsityPattern &,
- const unsigned int);
-template void MGDoFTools::make_flux_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
- SparsityPattern &,
- const unsigned int);
+ const unsigned int n_components = dof_handler.get_fe().n_components();
+ for (unsigned int l=0;l<nlevels;++l)
+ {
+ result[l].resize (n_components);
+
+ // special case for only one
+ // component. treat this first
+ // since it does not require any
+ // computations
+ if (n_components == 1)
+ {
+ result[l][0] = dof_handler.n_dofs(l);
+ } else {
+ // otherwise determine the number
+ // of dofs in each component
+ // separately. do so in parallel
+ std::vector<std::vector<bool> >
+ dofs_in_component (n_components,
+ std::vector<bool>(dof_handler.n_dofs(l),
+ false));
+ std::vector<std::vector<bool> >
+ component_select (n_components,
+ std::vector<bool>(n_components, false));
+ Threads::ThreadManager thread_manager;
+ for (unsigned int i=0; i<n_components; ++i)
+ {
+ void (*fun_ptr) (const unsigned int level,
+ const MGDoFHandler<dim> &,
+ const std::vector<bool> &,
+ std::vector<bool> &)
+ = &DoFTools::template extract_level_dofs<dim>;
+ component_select[i][i] = true;
+ Threads::spawn (thread_manager,
+ Threads::encapsulate (fun_ptr)
+ .collect_args (l, dof_handler, component_select[i],
+ dofs_in_component[i]));
+ };
+ thread_manager.wait();
+
+ // next count what we got
+ for (unsigned int i=0; i<n_components; ++i)
+ result[l][i] = std::count(dofs_in_component[i].begin(),
+ dofs_in_component[i].end(),
+ true);
+
+ // finally sanity check
+ Assert (std::accumulate (result[l].begin(),
+ result[l].end(), 0U)
+ ==
+ dof_handler.n_dofs(l),
+ ExcInternalError());
+ }
+ }
+}
+
+
+
+template<int dim, typename number>
+void
+MGTools::reinit_vector (const MGDoFHandler<dim>& mg_dof,
+ MGLevelObject<Vector<number> >& v)
+{
+ for (unsigned int level=v.get_minlevel();
+ level<=v.get_maxlevel();++level)
+ {
+ unsigned int n = mg_dof.n_dofs (level);
+ v[level].reinit(n);
+ }
+
+}
+
+
+
+
+
+template<int dim, typename number>
+void
+MGTools::reinit_vector (const MGDoFHandler<dim>& mg_dof,
+ MGLevelObject<BlockVector<number> >& v,
+ const std::vector<bool>& selected)
+{
+ const unsigned int ncomp = mg_dof.get_fe().n_components();
+
+ if (selected.size() == 0)
+ {
+ selected.resize(ncomp);
+ fill_n (selected.begin(), ncomp, true);
+ }
+
+ Assert (selected.size() == ncomp,
+ ExcDimensionMismatch(selected.size(), ncomp));
-template void MGDoFTools::make_flux_sparsity_pattern_edge (const MGDoFHandler<deal_II_dimension> &,
- SparsityPattern &,
- const unsigned int);
+ unsigned int n_selected = std::accumulate (selected.begin(),
+ selected.end(),
+ 0U);
+
+ std::vector<std::vector<unsigned int> >
+ ndofs(mg_dof.n_levels(), ncomp);
+
+ count_dofs_per_component (mg_dof, ndofs);
+
+ for (unsigned int level=v.get_minlevel();
+ level<=v.get_maxlevel();++level)
+ {
+ v[level].reinit(n_selected);
+ unsigned int k=0;
+ for (unsigned int i=0;i<ncomp;++i)
+ if (selected[i])
+ v[level].block(k++).reinit(ndofs[level][i]);
+ }
+}
+
+
+
+
+
+// explicit instantiations
+template void
+MGTools::make_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
+ SparsityPattern &,
+ const unsigned int);
+
+template void
+MGTools::make_flux_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
+ SparsityPattern &,
+ const unsigned int);
+
+template void
+MGTools::make_flux_sparsity_pattern_edge (const MGDoFHandler<deal_II_dimension> &,
+ SparsityPattern &,
+ const unsigned int);
+
+template void
+MGTools::make_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
+ CompressedSparsityPattern &,
+ const unsigned int);
+
+template void
+MGTools::make_flux_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
+ CompressedSparsityPattern &,
+ const unsigned int);
+
+template void
+MGTools::make_flux_sparsity_pattern_edge (const MGDoFHandler<deal_II_dimension> &,
+ CompressedSparsityPattern &,
+ const unsigned int);
+
+template void
+MGTools::make_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
+ BlockSparsityPattern &,
+ const unsigned int);
+template void
+MGTools::make_flux_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
+ BlockSparsityPattern &,
+ const unsigned int);
+
+template void
+MGTools::make_flux_sparsity_pattern_edge (const MGDoFHandler<deal_II_dimension> &,
+ BlockSparsityPattern &,
+ const unsigned int);
+
+template void
+MGTools::make_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
+ CompressedBlockSparsityPattern &,
+ const unsigned int);
+
+template void
+MGTools::make_flux_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
+ CompressedBlockSparsityPattern &,
+ const unsigned int);
+
+template void
+MGTools::make_flux_sparsity_pattern_edge (const MGDoFHandler<deal_II_dimension> &,
+ CompressedBlockSparsityPattern &,
+ const unsigned int);
#if deal_II_dimension > 1
-template void MGDoFTools::make_flux_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
- SparsityPattern &,
- const unsigned int,
- const FullMatrix<double>&,
- const FullMatrix<double>&);
+template void
+MGTools::make_flux_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
+ SparsityPattern &,
+ const unsigned int,
+ const FullMatrix<double>&,
+ const FullMatrix<double>&);
+
+template void
+MGTools::make_flux_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
+ CompressedSparsityPattern &,
+ const unsigned int,
+ const FullMatrix<double>&,
+ const FullMatrix<double>&);
+
+template void
+MGTools::make_flux_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
+ BlockSparsityPattern &,
+ const unsigned int,
+ const FullMatrix<double>&,
+ const FullMatrix<double>&);
+
+template void
+MGTools::make_flux_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
+ CompressedBlockSparsityPattern &,
+ const unsigned int,
+ const FullMatrix<double>&,
+ const FullMatrix<double>&);
+
#endif
+
+template void MGTools::reinit_vector (const MGDoFHandler<deal_II_dimension>&,
+ MGLevelObject<Vector<double> >&);
+template void MGTools::reinit_vector (const MGDoFHandler<deal_II_dimension>&,
+ MGLevelObject<Vector<float> >&);
+
+
+// template void MGTools::copy_to_mg (const MGDoFHandler<deal_II_dimension>&,
+// const MGTransferBase<Vector<double> >&,
+// MGLevelObject<Vector<double> >&,
+// const Vector<double>&);
+
+// template void MGTools::copy_from_mg(const MGDoFHandler<deal_II_dimension>&,
+// Vector<double>&,
+// const MGLevelObject<Vector<double> >&);
+
+// template void MGTools::copy_from_mg_add(const MGDoFHandler<deal_II_dimension>&,
+// Vector<double>&,
+// const MGLevelObject<Vector<double> >&);
+
+// template void MGTools::copy_to_mg (const MGDoFHandler<deal_II_dimension>&,
+// const MGTransferBase<Vector<double> >&,
+// MGLevelObject<Vector<double> >&,
+// const Vector<float>&);
+
+// template void MGTools::copy_from_mg(const MGDoFHandler<deal_II_dimension>&,
+// Vector<double>&,
+// const MGLevelObject<Vector<float> >&);
+
+// template void MGTools::copy_from_mg_add(const MGDoFHandler<deal_II_dimension>&,
+// Vector<double>&,
+// const MGLevelObject<Vector<float> >&);
+
+// template void MGTools::copy_to_mg (const MGDoFHandler<deal_II_dimension>&,
+// const MGTransferBase<Vector<double> >&,
+// MGLevelObject<Vector<double> >&,
+// const BlockVector<double>&);
+
+// template void MGTools::copy_from_mg(const MGDoFHandler<deal_II_dimension>&,
+// BlockVector<double>&,
+// const MGLevelObject<Vector<double> >&);
+
+// template void MGTools::copy_from_mg_add(const MGDoFHandler<deal_II_dimension>&,
+// BlockVector<double>&,
+// const MGLevelObject<Vector<double> >&);
+
+// template void MGTools::copy_to_mg (const MGDoFHandler<deal_II_dimension>&,
+// const MGTransferBase<Vector<double> >&,
+// MGLevelObject<Vector<double> >&,
+// const BlockVector<float>&);
+
+// template void MGTools::copy_from_mg(const MGDoFHandler<deal_II_dimension>&,
+// BlockVector<double>&,
+// const MGLevelObject<Vector<float> >&);
+
+// template void MGTools::copy_from_mg_add(const MGDoFHandler<deal_II_dimension>&,
+// BlockVector<double>&,
+// const MGLevelObject<Vector<float> >&);
+
+template void MGTools::count_dofs_per_component (const MGDoFHandler<deal_II_dimension>&,
+ std::vector<std::vector<unsigned int> >&);
+
#include <multigrid/mg_dof_accessor.h>
#include <grid/tria_iterator.h>
#include <fe/fe.h>
-#include <multigrid/multigrid.h>
+//#include <multigrid/multigrid.h>
#include <multigrid/mg_smoother.h>
#include <multigrid/mg_smoother.templates.h>
#include <lac/vector.h>
+#include <lac/block_vector.h>
+#include <lac/sparse_matrix.h>
+#include <lac/block_sparse_matrix.h>
#include <algorithm>
//////////////////////////////////////////////////////////////////////
-MGSmootherBase::~MGSmootherBase()
-{};
-
-
-//////////////////////////////////////////////////////////////////////
-
-
-void
-MGSmootherIdentity::smooth (const unsigned int,
- Vector<double> &,
- const Vector<double> &) const
-{};
-
#if deal_II_dimension == 1
-MGSmoother::MGSmoother (const MGDoFHandler<1> &/*mg_dof*/, unsigned int steps)
+MGSmootherContinuous::MGSmootherContinuous (const MGDoFHandler<1> &/*mg_dof*/,
+ unsigned int steps)
:
steps(steps)
{
#if deal_II_dimension > 1
template <int dim>
-MGSmoother::MGSmoother (const MGDoFHandler<dim> &mg_dof, unsigned int steps)
+MGSmootherContinuous::MGSmootherContinuous (const MGDoFHandler<dim> &mg_dof,
+ unsigned int steps)
:
steps(steps)
{
//////////////////////////////////////////////////////////////////////
+template <class VECTOR>
void
-MGSmoother::set_zero_interior_boundary (const unsigned int level,
- Vector<double> &u) const
+MGSmootherContinuous::set_zero_interior_boundary (const unsigned int level,
+ VECTOR& u) const
{
if (level==0)
return;
// 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<deal_II_dimension>&, unsigned int);
+template MGSmootherContinuous::MGSmootherContinuous (const MGDoFHandler<deal_II_dimension>&, unsigned int);
#endif
template
-MGSmootherRelaxation<float>
-::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension> &mg_dof,
- const MGLevelObject<SparseMatrix<float> > &matrix,
- const function_ptr relaxation,
- const unsigned int steps,
- const double omega);
+void MGSmootherContinuous::set_zero_interior_boundary (const unsigned int,
+ Vector<double>&) const;
+
+template
+void MGSmootherContinuous::set_zero_interior_boundary (const unsigned int,
+ Vector<float>&) const;
+
+template
+void MGSmootherContinuous::set_zero_interior_boundary (const unsigned int,
+ BlockVector<double>&) const;
+
+template
+void MGSmootherContinuous::set_zero_interior_boundary (const unsigned int,
+ BlockVector<float>&) const;
+
+
+
+template
+MGSmootherRelaxation<SparseMatrix<float>, Vector<float> >
+::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension>&,
+ const MGLevelObject<SparseMatrix<float> >&,
+ const function_ptr,
+ const unsigned int,
+ const double);
+
+template
+MGSmootherRelaxation<SparseMatrix<float>, Vector<double> >
+::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension>&,
+ const MGLevelObject<SparseMatrix<float> >&,
+ const function_ptr,
+ const unsigned int,
+ const double);
+
+template
+MGSmootherRelaxation<SparseMatrix<double>, Vector<float> >
+::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension>&,
+ const MGLevelObject<SparseMatrix<double> >&,
+ const function_ptr,
+ const unsigned int,
+ const double);
+
+template
+MGSmootherRelaxation<SparseMatrix<double>, Vector<double> >
+::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension>&,
+ const MGLevelObject<SparseMatrix<double> >&,
+ const function_ptr,
+ const unsigned int,
+ const double);
+
+
+template
+MGSmootherRelaxation<BlockSparseMatrix<float>, BlockVector<float> >
+::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension>&,
+ const MGLevelObject<BlockSparseMatrix<float> >&,
+ const function_ptr,
+ const unsigned int,
+ const double);
+
+template
+MGSmootherRelaxation<BlockSparseMatrix<float>, BlockVector<double> >
+::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension>&,
+ const MGLevelObject<BlockSparseMatrix<float> >&,
+ const function_ptr,
+ const unsigned int,
+ const double);
+
+template
+MGSmootherRelaxation<BlockSparseMatrix<double>, BlockVector<float> >
+::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension>&,
+ const MGLevelObject<BlockSparseMatrix<double> >&,
+ const function_ptr,
+ const unsigned int,
+ const double);
template
-MGSmootherRelaxation<double>
-::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension> &mg_dof,
- const MGLevelObject<SparseMatrix<double> > &matrix,
- const function_ptr relaxation,
- const unsigned int steps,
- const double omega);
+MGSmootherRelaxation<BlockSparseMatrix<double>, BlockVector<double> >
+::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension>&,
+ const MGLevelObject<BlockSparseMatrix<double> >&,
+ const function_ptr,
+ const unsigned int,
+ const double);
-#include <grid/tria.h>
-#include <multigrid/mg_dof_handler.h>
-#include <multigrid/mg_dof_accessor.h>
-#include <grid/tria_iterator.h>
-#include <fe/fe.h>
-#include <multigrid/multigrid.h>
-#include <multigrid/multigrid.templates.h>
-#include <multigrid/mg_smoother.h>
#include <lac/vector.h>
+#include <lac/sparse_matrix.h>
+#include <lac/block_sparse_matrix.h>
+#include <multigrid/mg_transfer.h>
MGTransferPrebuilt::~MGTransferPrebuilt ()
};
+
+
+
+MGTransferBlock::~MGTransferBlock ()
+{};
+
+
+void MGTransferBlock::prolongate (const unsigned int to_level,
+ BlockVector<double> &dst,
+ const BlockVector<double> &src) const
+{
+ Assert ((to_level >= 1) && (to_level<=prolongation_matrices.size()),
+ ExcIndexRange (to_level, 1, prolongation_matrices.size()+1));
+
+ unsigned int k=0;
+ for (unsigned int b=0; b<src.n_blocks();++b)
+ {
+ if (!selected[k])
+ ++k;
+ prolongation_matrices[to_level-1].block(k,k).vmult (dst.block(b), src.block(b));
+ ++k;
+ }
+};
+
+
+void MGTransferBlock::restrict_and_add (const unsigned int from_level,
+ BlockVector<double> &dst,
+ const BlockVector<double> &src) const
+{
+ Assert ((from_level >= 1) && (from_level<=prolongation_matrices.size()),
+ ExcIndexRange (from_level, 1, prolongation_matrices.size()+1));
+
+ unsigned int k=0;
+ for (unsigned int b=0; b<src.n_blocks();++b)
+ {
+ if (!selected[k])
+ ++k;
+ prolongation_matrices[from_level-1].block(k,k).Tvmult_add (dst.block(b), src.block(b));
+ ++k;
+ }
+};
+
+
+
+
+MGTransferSelect::~MGTransferSelect ()
+{};
+
+
+void MGTransferSelect::prolongate (const unsigned int to_level,
+ Vector<double> &dst,
+ const Vector<double> &src) const
+{
+ Assert ((to_level >= 1) && (to_level<=prolongation_matrices.size()),
+ ExcIndexRange (to_level, 1, prolongation_matrices.size()+1));
+
+ prolongation_matrices[to_level-1].block(selected, selected)
+ .vmult (dst, src);
+};
+
+
+void MGTransferSelect::restrict_and_add (const unsigned int from_level,
+ Vector<double> &dst,
+ const Vector<double> &src) const
+{
+ Assert ((from_level >= 1) && (from_level<=prolongation_matrices.size()),
+ ExcIndexRange (from_level, 1, prolongation_matrices.size()+1));
+
+ prolongation_matrices[from_level-1].block(selected, selected)
+ .Tvmult_add (dst, src);
+};
-//---------------------------- multigrid.cc ---------------------------
+//-----------------------------------------------------------------------
// $Id$
// Version: $Name$
//
// to the file deal.II/doc/license.html for the text and
// further information on this license.
//
-//---------------------------- multigrid.cc ---------------------------
+//-----------------------------------------------------------------------
+#include <lac/vector.h>
+#include <lac/sparse_matrix.h>
+#include <lac/block_sparse_matrix.h>
#include <grid/tria.h>
-#include <dofs/dof_constraints.h>
-#include <multigrid/mg_dof_handler.h>
-#include <multigrid/mg_dof_accessor.h>
#include <grid/tria_iterator.h>
+#include <dofs/dof_tools.h>
#include <fe/fe.h>
-#include <multigrid/multigrid.h>
-#include <multigrid/multigrid.templates.h>
-#include <multigrid/mg_smoother.h>
-#include <lac/vector.h>
-
-
-/* --------------------------------- MG ------------------------------------ */
-
-template <int dim>
-Multigrid<dim>::Multigrid (const MGDoFHandler<dim> &mg_dof_handler,
- const ConstraintMatrix &constraints,
- const MGLevelObject<SparsityPattern> &level_sparsities,
- const MGLevelObject<SparseMatrix<double> > &level_matrices,
- const MGTransferBase &transfer,
- const unsigned int minlevel,
- const unsigned int maxlevel)
- :
- MGBase(transfer, minlevel,
- std::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)
-{};
-
-
-template <int dim>
-void
-Multigrid<dim>::level_vmult (const unsigned int level,
- Vector<double> &result,
- const Vector<double> &u,
- const Vector<double> &/* rhs */)
-{
- (*level_matrices)[level].vmult(result,u);
- result.scale(-1.);
-}
+#include <multigrid/mg_dof_handler.h>
+#include <multigrid/mg_dof_accessor.h>
+#include <multigrid/mg_transfer.h>
+#include <multigrid/mg_dof_tools.h>
/* ----------------------------- MGTransferPrebuilt ------------------------ */
const unsigned int n_levels = mg_dof.get_tria().n_levels();
const unsigned int dofs_per_cell = mg_dof.get_fe().dofs_per_cell;
+ sizes.resize(n_levels);
+ for (unsigned int l=0;l<n_levels;++l)
+ sizes[l] = mg_dof.n_dofs(l);
+
// 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);
+ prolongation_sparsities.resize (n_levels-1);
+ prolongation_matrices.resize (n_levels-1);
-// two fields which will store the
+ // two fields which will store the
// indices of the multigrid dofs
// for a cell and one of its children
- std::vector<unsigned int> dof_indices_mother (dofs_per_cell);
+ std::vector<unsigned int> dof_indices_parent (dofs_per_cell);
std::vector<unsigned int> dof_indices_child (dofs_per_cell);
// for each level: first build the sparsity
// level which have children
for (unsigned int level=0; level<n_levels-1; ++level)
{
- prolongation_sparsities.push_back (SparsityPattern());
// reset the dimension of the structure.
// note that for the number of entries
- // per row, the number of mother dofs
+ // per row, the number of parent dofs
// coupling to a child dof is
// necessary. this, of course, is the
// number of degrees of freedom per
// cell
- prolongation_sparsities.back().reinit (mg_dof.n_dofs(level+1),
- mg_dof.n_dofs(level),
+ prolongation_sparsities[level].reinit (sizes[level+1],
+ sizes[level],
//TODO:[WB,GK] evil hack, must be corrected!
- dofs_per_cell+1);
+ dofs_per_cell+1);
for (typename MGDoFHandler<dim>::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);
+ cell->get_mg_dof_indices (dof_indices_parent);
for (unsigned int child=0;
child<GeometryInfo<dim>::children_per_cell; ++child)
// now tag the entries in the
// matrix which will be used
- // for this pair of mother/child
+ // for this pair of parent/child
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
if (prolongation(i,j) != 0)
{
prolongation_sparsities[level].add
(dof_indices_child[i],
- dof_indices_mother[j]);
+ dof_indices_parent[j]);
};
};
};
prolongation_sparsities[level].compress ();
- prolongation_matrices.push_back (SparseMatrix<double>());
prolongation_matrices[level].reinit (prolongation_sparsities[level]);
// now actually build the matrices
cell != mg_dof.end(level); ++cell)
if (cell->has_children())
{
- cell->get_mg_dof_indices (dof_indices_mother);
+ cell->get_mg_dof_indices (dof_indices_parent);
for (unsigned int child=0;
child<GeometryInfo<dim>::children_per_cell; ++child)
{
prolongation_matrices[level].set
(dof_indices_child[i],
- dof_indices_mother[j],
+ dof_indices_parent[j],
prolongation(i,j));
};
};
};
-// explicit instatations
-template class Multigrid<deal_II_dimension>;
-template
-void MGTransferPrebuilt::build_matrices (const MGDoFHandler<deal_II_dimension> &mg_dof);
+/* ----------------------------- MGTransferBlock ------------------------ */
-template
-void Multigrid<deal_II_dimension>::copy_to_mg (const Vector<float>& src);
-template
-void Multigrid<deal_II_dimension>::copy_to_mg (const Vector<double>& src);
+template <int dim>
+void MGTransferBlockBase::build_matrices (const MGDoFHandler<dim> &mg_dof,
+ std::vector<bool> select)
+{
+ const FiniteElement<dim>& fe = mg_dof.get_fe();
+ const unsigned int n_components = fe.n_components();
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+ const unsigned int n_levels = mg_dof.get_tria().n_levels();
+
+ selected = select;
+
+ Assert (selected.size() == n_components,
+ ExcDimensionMismatch(selected.size(), n_components))
-template
-void Multigrid<deal_II_dimension>::copy_from_mg (Vector<float>& src) const;
+ sizes.resize(n_levels);
+ MGTools::count_dofs_per_component(mg_dof, sizes);
+
+ // reset the size of the array of
+ // matrices
+ prolongation_sparsities.resize (n_levels-1);
+ prolongation_matrices.resize (n_levels-1);
+
+ // two fields which will store the
+ // indices of the multigrid dofs
+ // for a cell and one of its children
+ std::vector<unsigned int> dof_indices_parent (dofs_per_cell);
+ std::vector<unsigned int> 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<n_levels-1; ++level)
+ {
+ // reset the dimension of the
+ // structure. note that for
+ // the number of entries per
+ // row, the number of parent
+ // dofs coupling to a child dof
+ // is necessary. this, is the
+ // number of degrees of freedom
+ // per cell
+ prolongation_sparsities[level].reinit (n_components, n_components);
+ for (unsigned int i=0;i<n_components;++i)
+ for (unsigned int j=0;j<n_components;++j)
+ if (i==j)
+ prolongation_sparsities[level].block(i,j).reinit(
+ sizes[level+1][i],
+ sizes[level][j],
+//TODO:[GK] Split this by component to save memory
+ dofs_per_cell);
+ else
+ prolongation_sparsities[level].block(i,j).reinit(
+ sizes[level+1][i],
+ sizes[level][j],
+ 0);
+
+ prolongation_sparsities[level].collect_sizes();
+
+ for (typename MGDoFHandler<dim>::cell_iterator cell=mg_dof.begin(level);
+ cell != mg_dof.end(level); ++cell)
+ if (cell->has_children())
+ {
+ cell->get_mg_dof_indices (dof_indices_parent);
+
+ for (unsigned int child=0;
+ child<GeometryInfo<dim>::children_per_cell; ++child)
+ {
+ // set an alias to the
+ // prolongation matrix for
+ // this child
+ const FullMatrix<double> &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 parent/child
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ if (prolongation(i,j) != 0)
+ {
+ const unsigned int icomp = fe.system_to_component_index(i).first;
+ const unsigned int jcomp = fe.system_to_component_index(j).first;
+ if ((icomp==jcomp) && selected[icomp])
+ prolongation_sparsities[level].add
+ (dof_indices_child[i],
+ dof_indices_parent[j]);
+ };
+ };
+ };
+ prolongation_sparsities[level].compress ();
+
+ prolongation_matrices[level].reinit (prolongation_sparsities[level]);
+
+ // now actually build the matrices
+ for (typename MGDoFHandler<dim>::cell_iterator cell=mg_dof.begin(level);
+ cell != mg_dof.end(level); ++cell)
+ if (cell->has_children())
+ {
+ cell->get_mg_dof_indices (dof_indices_parent);
+
+ for (unsigned int child=0;
+ child<GeometryInfo<dim>::children_per_cell; ++child)
+ {
+ // set an alias to the
+ // prolongation matrix for
+ // this child
+ const FullMatrix<double> &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<dofs_per_cell; ++i)
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ if (prolongation(i,j) != 0)
+ {
+ const unsigned int icomp = fe.system_to_component_index(i).first;
+ const unsigned int jcomp = fe.system_to_component_index(j).first;
+ if ((icomp==jcomp) && selected[icomp])
+ prolongation_matrices[level].set
+ (dof_indices_child[i],
+ dof_indices_parent[j],
+ prolongation(i,j));
+ };
+ };
+ };
+ };
+ // Finally, fill some index vectors
+ // for later use.
+ mg_component_start = sizes;
+ // Compute start indices from sizes
+ for (unsigned int l=0;l<mg_component_start.size();++l)
+ {
+ unsigned int k=0;
+ for (unsigned int i=0;i<mg_component_start[l].size();++i)
+ {
+ const unsigned int t=mg_component_start[l][i];
+ mg_component_start[l][i] = k;
+ k += t;
+ }
+ }
+
+ component_start.resize(n_components);
+ DoFTools::count_dofs_per_component(mg_dof, component_start);
+
+ unsigned int k=0;
+ for (unsigned int i=0;i<component_start.size();++i)
+ {
+ const unsigned int t=component_start[i];
+ component_start[i] = k;
+ k += t;
+ }
+
+};
+
+
+template <int dim>
+void MGTransferBlock::build_matrices (const MGDoFHandler<dim> &mg_dof,
+ std::vector<bool> select)
+{
+ if (select.size() == 0)
+ select = std::vector<bool> (mg_dof.get_fe().n_components(), true);
+
+ MGTransferBlockBase::build_matrices (mg_dof, select);
+}
+
+
+template <int dim>
+void MGTransferSelect::build_matrices (const MGDoFHandler<dim> &mg_dof,
+ unsigned int select)
+{
+ selected = select;
+ std::vector<bool> s(mg_dof.get_fe().n_components(), false);
+ s[select] = true;
+
+ MGTransferBlockBase::build_matrices (mg_dof, s);
+}
+
+
+
+// explicit instatations
template
-void Multigrid<deal_II_dimension>::copy_from_mg (Vector<double>& src) const;
+void MGTransferPrebuilt::build_matrices (const MGDoFHandler<deal_II_dimension> &mg_dof);
+template
+void MGTransferBlock::build_matrices (const MGDoFHandler<deal_II_dimension> &mg_dof,
+ std::vector<bool>);
+template
+void MGTransferSelect::build_matrices (const MGDoFHandler<deal_II_dimension> &mg_dof,
+ unsigned int);
{
mgstruct[i].reinit(mgdof.n_dofs(i), mgdof.n_dofs(i),
mgdof.max_couplings_between_dofs());
- MGDoFTools::make_sparsity_pattern(mgdof, mgstruct[i], i);
+ MGTools::make_sparsity_pattern(mgdof, mgstruct[i], i);
mgstruct[i].compress();
mgA[i].reinit(mgstruct[i]);
projects / dealii.git / commitdiff