From: guido Date: Thu, 29 Nov 2001 11:55:28 +0000 (+0000) Subject: multigrid cleanup branch merged X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=93ffb33f2c269b8349b8b36d362f195f2957b321;p=dealii-svn.git multigrid cleanup branch merged git-svn-id: https://svn.dealii.org/trunk@5307 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/multigrid/mg_base.h b/deal.II/deal.II/include/multigrid/mg_base.h index 3c99b73188..00e12ab33e 100644 --- a/deal.II/deal.II/include/multigrid/mg_base.h +++ b/deal.II/deal.II/include/multigrid/mg_base.h @@ -17,12 +17,8 @@ #include #include #include -#include - #include -class MGSmootherBase; - /** * An array with an object for each level. The purpose of this class @@ -105,39 +101,6 @@ class MGLevelObject : public Subscriptor }; -/** - * 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 &dst, - const Vector &src) const = 0; -}; - /** * Coarse grid solver using LAC iterative methods. @@ -150,8 +113,8 @@ class MGCoarseGridSolver : public Subscriptor * * @author Guido Kanschat, 1999 */ -template -class MGCoarseGridLACIteration : public MGCoarseGridSolver +template > +class MGCoarseGridLACIteration : public Subscriptor { public: /** @@ -171,9 +134,9 @@ class MGCoarseGridLACIteration : public MGCoarseGridSolver * matrix, vectors and * preconditioner. */ - virtual void operator() (const unsigned int level, - Vector &dst, - const Vector &src) const; + void operator() (const unsigned int level, + VECTOR &dst, + const VECTOR &src) const; private: /** * Reference to the solver. @@ -200,303 +163,63 @@ class MGCoarseGridLACIteration : public MGCoarseGridSolver * * @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 &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; -}; - - -/** - * 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 > defect; - - /** - * Auxiliary vector. - */ - MGLevelObject > solution; - - /** - * Auxiliary vector. - */ - MGLevelObject > t; - - /** - * Prolongation and restriction object. - */ - SmartPointer 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 &dst, - const Vector &src, - const Vector &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 &dst, - const Vector &src); - - /** - * Print a level vector using - * @ref{DoFHandler}. - */ - virtual void print_vector (const unsigned int level, - const Vector& 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 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 multigrid; - - /** - * The pre-smoothing object. - */ - SmartPointer pre; +// template +// 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 post; - - /** - * The coarse grid solver. - */ - SmartPointer coarse; -}; /* ------------------------------------------------------------------- */ @@ -574,10 +297,11 @@ MGLevelObject::get_maxlevel () const /* ------------------ Functions for MGCoarseGridLACIteration ------------ */ -template -MGCoarseGridLACIteration::MGCoarseGridLACIteration(SOLVER &s, - const MATRIX &m, - const PRECOND &p) +template +MGCoarseGridLACIteration +::MGCoarseGridLACIteration(SOLVER& s, + const MATRIX &m, + const PRECOND &p) : solver(s), matrix(&m), @@ -585,43 +309,17 @@ MGCoarseGridLACIteration::MGCoarseGridLACIteration(SOLVER {}; -template +template void -MGCoarseGridLACIteration::operator() (const unsigned int /* level */, - Vector &dst, - const Vector &src) const +MGCoarseGridLACIteration +::operator() (const unsigned int /* level */, + VECTOR &dst, + const VECTOR &src) const { solver.solve(*matrix, dst, src, precondition); } -/* ------------------------------------------------------------------- */ - - -template -PreconditionMG::PreconditionMG(MG &mg, - const MGSmootherBase &pre, - const MGSmootherBase &post, - const MGCoarseGridSolver &coarse) - : - multigrid(&mg), - pre(&pre), - post(&post), - coarse(&coarse) -{} - - -template -void -PreconditionMG::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 diff --git a/deal.II/deal.II/include/multigrid/mg_dof_tools.h b/deal.II/deal.II/include/multigrid/mg_dof_tools.h index 561874100b..bd3802e9dc 100644 --- a/deal.II/deal.II/include/multigrid/mg_dof_tools.h +++ b/deal.II/deal.II/include/multigrid/mg_dof_tools.h @@ -13,6 +13,11 @@ #ifndef __deal2__mg_dof_tools_h #define __deal2__mg_dof_tools_h +template class MGLevelObject; +template class MGDoFHandler; +template class Vector; +template class BlockVector; +template class FullMatrix; #include @@ -26,11 +31,11 @@ * 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: /** @@ -47,7 +52,7 @@ class MGDoFTools * hanging nodes here, since only * one level is considered. */ - template + template static void make_sparsity_pattern (const MGDoFHandler &dof_handler, SparsityPattern &sparsity, @@ -59,7 +64,7 @@ class MGDoFTools * @see make_sparsity_pattern * $see DoFTools */ - template + template static void make_flux_sparsity_pattern (const MGDoFHandler &dof_handler, SparsityPattern &sparsity, @@ -72,7 +77,7 @@ class MGDoFTools * @see{make_flux_sparsity_pattern} * @see{DoFTools} */ - template + template static void make_flux_sparsity_pattern_edge (const MGDoFHandler &dof_handler, SparsityPattern &sparsity, @@ -95,7 +100,7 @@ class MGDoFTools * for the couplings occuring in * fluxes. */ - template + template static void make_flux_sparsity_pattern (const MGDoFHandler &dof, SparsityPattern &sparsity, @@ -103,6 +108,63 @@ class MGDoFTools const FullMatrix& int_mask, const FullMatrix& 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 + static void count_dofs_per_component (const MGDoFHandler& mg_dof, + std::vector >& 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 + static void + reinit_vector (const MGDoFHandler& mg_dof, + MGLevelObject >& 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 + static void + reinit_vector (const MGDoFHandler& mg_dof, + MGLevelObject >& vector, + const vector& selected); }; diff --git a/deal.II/deal.II/include/multigrid/mg_smoother.h b/deal.II/deal.II/include/multigrid/mg_smoother.h index 715e0e35f1..0c9fbeae82 100644 --- a/deal.II/deal.II/include/multigrid/mg_smoother.h +++ b/deal.II/deal.II/include/multigrid/mg_smoother.h @@ -22,47 +22,6 @@ template 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 &u, - const Vector &rhs) const = 0; -}; - - /** * Smoother doing nothing. This class is not useful for many applications other * than for testing some multigrid procedures. Also some applications might @@ -71,12 +30,12 @@ class MGSmootherBase : public Subscriptor * * @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 @@ -84,33 +43,33 @@ class MGSmootherIdentity : public MGSmootherBase * operator equals the null * operator. */ - virtual void smooth (const unsigned int level, - Vector &u, - const Vector &rhs) const; + template + 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: @@ -122,15 +81,19 @@ class MGSmoother : public MGSmootherBase * 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); /** @@ -145,7 +108,7 @@ class MGSmoother : public MGSmootherBase * steps to be executed by @p{smooth}. */ template - MGSmoother (const MGDoFHandler &mg_dof, + MGSmootherContinuous (const MGDoFHandler &mg_dof, unsigned int steps); /** @@ -158,12 +121,15 @@ class MGSmoother : public MGSmootherBase * has no interior boundaries, this * function does nothing in this case. */ - void set_zero_interior_boundary (const unsigned int level, - Vector &u) const; + template + 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? */ @@ -174,6 +140,7 @@ class MGSmoother : public MGSmootherBase * Number of smoothing steps. */ unsigned int steps; + /** * For each level, we store a list of * degree of freedom indices which are @@ -191,27 +158,36 @@ class MGSmoother : public MGSmootherBase /** - * 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 -class MGSmootherRelaxation : public MGSmoother +template +class MGSmootherRelaxation : public MGSmootherContinuous { public: /** * Type of the smoothing * function of the matrix. */ - typedef void - (SparseMatrix::* function_ptr)(Vector&, const Vector&, - typename SparseMatrix::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. @@ -225,25 +201,26 @@ class MGSmootherRelaxation : public MGSmoother */ template - MGSmootherRelaxation(const MGDoFHandler &mg_dof, - const MGLevelObject > &matrix, - const function_ptr function, - const unsigned int steps, - const double omega = 1.); + MGSmootherRelaxation(const MGDoFHandler& mg_dof, + const MGLevelObject& matrix, + const function_ptr function, + const unsigned int steps, + const double omega = 1.); /** - * Implementation of the interface in @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 &u, - const Vector &rhs) const; + void smooth (const unsigned int level, + VECTOR& u, + const VECTOR& rhs) const; private: /** * Pointer to the matrices. */ - SmartPointer > > matrix; + SmartPointer > matrix; /** * Pointer to the relaxation function. @@ -269,16 +246,23 @@ class MGSmootherRelaxation : public MGSmoother /* ------------------------------- Inline functions -------------------------- */ +template +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; } diff --git a/deal.II/deal.II/include/multigrid/mg_smoother.templates.h b/deal.II/deal.II/include/multigrid/mg_smoother.templates.h index 6d33930854..846cd969ed 100644 --- a/deal.II/deal.II/include/multigrid/mg_smoother.templates.h +++ b/deal.II/deal.II/include/multigrid/mg_smoother.templates.h @@ -1,26 +1,28 @@ // $Id$ -template +template template -MGSmootherRelaxation::MGSmootherRelaxation (const MGDoFHandler &mg_dof, - const MGLevelObject > &matrix, - const function_ptr relaxation, - const unsigned int steps, - const double omega) +MGSmootherRelaxation +::MGSmootherRelaxation (const MGDoFHandler& mg_dof, + const MGLevelObject& 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 +template void -MGSmootherRelaxation::smooth (const unsigned int level, - Vector &u, - const Vector &rhs) const +MGSmootherRelaxation +::smooth (const unsigned int level, + VECTOR& u, + const VECTOR& rhs) const { h1.reinit(u); h2.reinit(u); diff --git a/deal.II/deal.II/include/multigrid/mg_tools.templates.h b/deal.II/deal.II/include/multigrid/mg_tools.templates.h new file mode 100644 index 0000000000..ca1aac52b2 --- /dev/null +++ b/deal.II/deal.II/include/multigrid/mg_tools.templates.h @@ -0,0 +1,471 @@ +//---------------------------------------------------------------------------- +// $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 +#include +#include +#include +#include + +#include + + +//TODO:[?] This function needs to be specially implemented, since in 2d mode we use faces +#if deal_II_dimension == 1 + +template +void +MGTransferPrebuilt::copy_to_mg (const MGDoFHandler&, + MGLevelObject >&, + const InVector&) const +{ + Assert(false, ExcNotImplemented()); +} + +#else + + +template +void +MGTransferPrebuilt::copy_to_mg (const MGDoFHandler& mg_dof_handler, + MGLevelObject >& 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 global_dof_indices (dofs_per_cell); + std::vector level_dof_indices (dofs_per_cell); + std::vector level_face_indices (dofs_per_face); + + // traverse the grid top-down + // (i.e. starting with the most + // refined grid). this way, we can + // always get that part of one + // level of the output vector which + // corresponds to a region which is + // more refined, by restriction of + // the respective vector on the + // next finer level, which we then + // already have built. + for (int level=maxlevel; level>=static_cast(minlevel); --level) + { + typename MGDoFHandler::active_cell_iterator + level_cell = mg_dof_handler.begin_active(level); + const typename MGDoFHandler::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& 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::faces_per_cell; ++face_n) + { + const typename MGDoFHandler::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(level) < maxlevel) + { + restrict_and_add (level+1, dst[level], dst[level+1]); + } + }; +}; + +#endif + + +template +void +MGTransferPrebuilt::copy_from_mg(const MGDoFHandler& mg_dof_handler, + OutVector &dst, + const MGLevelObject > &src) const +{ + const unsigned int dofs_per_cell = mg_dof_handler.get_fe().dofs_per_cell; + + std::vector global_dof_indices (dofs_per_cell); + std::vector level_dof_indices (dofs_per_cell); + + typename MGDoFHandler::active_cell_iterator + level_cell = mg_dof_handler.begin_active(); + const typename MGDoFHandler::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& 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; iset_zero(dst); +} + +template +void +MGTransferPrebuilt::copy_from_mg_add(const MGDoFHandler& mg_dof_handler, + OutVector &dst, + const MGLevelObject > &src) const +{ + const unsigned int dofs_per_cell = mg_dof_handler.get_fe().dofs_per_cell; + + std::vector global_dof_indices (dofs_per_cell); + std::vector level_dof_indices (dofs_per_cell); + + typename MGDoFHandler::active_cell_iterator + level_cell = mg_dof_handler.begin_active(); + const typename MGDoFHandler::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& 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; iset_zero(dst); +} + + +//TODO:[?] This function needs to be specially implemented, since in 2d mode we use faces +#if deal_II_dimension == 1 + +template +void +MGTransferSelect::copy_to_mg (const MGDoFHandler&, + MGLevelObject >&, + const InVector&) const +{ + Assert(false, ExcNotImplemented()); +} + +#else + + +template +void +MGTransferSelect::copy_to_mg (const MGDoFHandler& mg_dof_handler, + MGLevelObject >& dst, + const InVector& osrc) const +{ + // Make src a real finite element function + InVector src = osrc; +// constraints->distribute(src); + + const FiniteElement& 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 global_dof_indices (dofs_per_cell); + std::vector level_dof_indices (dofs_per_cell); + std::vector 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(minlevel); --level) + { + // Start of treated component + // for this level + const unsigned int level_start = mg_component_start[level][selected]; + + typename MGDoFHandler::active_cell_iterator + level_cell = mg_dof_handler.begin_active(level); + const typename MGDoFHandler::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& 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::faces_per_cell; ++face_n) + { + const typename MGDoFHandler::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(level) < maxlevel) + { + restrict_and_add (level+1, dst[level], dst[level+1]); + } + }; +}; + +#endif + + +template +void +MGTransferSelect::copy_from_mg(const MGDoFHandler& mg_dof_handler, + OutVector &dst, + const MGLevelObject > &src) const +{ + + const FiniteElement& fe = mg_dof_handler.get_fe(); + const unsigned int dofs_per_cell = fe.dofs_per_cell; + + std::vector global_dof_indices (dofs_per_cell); + std::vector level_dof_indices (dofs_per_cell); + + typename MGDoFHandler::active_cell_iterator + level_cell = mg_dof_handler.begin_active(); + const typename MGDoFHandler::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& 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; iset_zero(dst); +} + +template +void +MGTransferSelect::copy_from_mg_add(const MGDoFHandler& mg_dof_handler, + OutVector &dst, + const MGLevelObject > &src) const +{ + + const FiniteElement& fe = mg_dof_handler.get_fe(); + const unsigned int dofs_per_cell = fe.dofs_per_cell; + + std::vector global_dof_indices (dofs_per_cell); + std::vector level_dof_indices (dofs_per_cell); + + // Start of input vector inside a + // block vector. + const unsigned int start = component_start[selected]; + + typename MGDoFHandler::active_cell_iterator + level_cell = mg_dof_handler.begin_active(); + const typename MGDoFHandler::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& 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; iset_zero(dst); +} + + +#endif diff --git a/deal.II/deal.II/include/multigrid/mg_transfer.h b/deal.II/deal.II/include/multigrid/mg_transfer.h new file mode 100644 index 0000000000..1aeda705ff --- /dev/null +++ b/deal.II/deal.II/include/multigrid/mg_transfer.h @@ -0,0 +1,450 @@ +//---------------------------------------------------------------------- +// $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 +#include +#include +#include + +template 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 > +{ + public: + /** + * Destructor. + */ + virtual ~MGTransferPrebuilt (); + + /** + * Actually build the prolongation + * matrices for each level. + */ + template + void build_matrices (const MGDoFHandler &mg_dof); + + /** + * Prolongate a vector from level + * @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; + + /** + * 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; + + /** + * Transfer from a vector on the + * global grid to vectors defined + * on each of the levels + * separately, i.a. an @p{MGVector}. + */ + template + void + copy_to_mg (const MGDoFHandler& mg_dof, + MGLevelObject > &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}. In order to + * keep the result consistent, + * constrained degrees of freedom + * are set to zero. + */ + template + void + copy_from_mg (const MGDoFHandler& mg_dof, + OutVector &dst, + const MGLevelObject > &src) const; + + /** + * Add a multi-level vector to a + * normal vector. + * + * Works as the previous + * function, but probably not for + * continuous elements. + */ + template + void + copy_from_mg_add (const MGDoFHandler& mg_dof, + OutVector &dst, + const MGLevelObject > &src) const; + + private: + + /** + * Selected component. + */ + unsigned int selected; + + /** + * Sizes of the multi-level vectors. + */ + std::vector sizes; + + + std::vector prolongation_sparsities; + + /** + * The actual prolongation matrix. + * column indices belong to the + * dof indices of the mother cell, + * i.e. the coarse level. + * while row indices belong to the + * child cell, i.e. the fine level. + */ + std::vector > 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 + void build_matrices (const MGDoFHandler &mg_dof, + std::vector selected); + + /** + * Flag of selected components. + */ + std::vector selected; + + /** + * Sizes of the multi-level vectors. + */ + std::vector > sizes; + + /** + * Start index of each component. + */ + std::vector component_start; + + /** + * Start index of eah component on + * all levels. + */ + std::vector > mg_component_start; + + private: + + std::vector 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 > 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 >, + 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 + void build_matrices (const MGDoFHandler &mg_dof, + std::vector selected = std::vector()); + + /** + * 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 &dst, + const BlockVector &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 &dst, + const BlockVector &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 >, + 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 + void build_matrices (const MGDoFHandler &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 &dst, + const Vector &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 &dst, + const Vector &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 + void + copy_to_mg (const MGDoFHandler& mg_dof, + MGLevelObject > &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}. In order to + * keep the result consistent, + * constrained degrees of freedom + * are set to zero. + */ + template + void + copy_from_mg (const MGDoFHandler& mg_dof, + OutVector &dst, + const MGLevelObject > &src) const; + + /** + * Add a multi-level vector to a + * normal vector. + * + * Works as the previous + * function, but probably not for + * continuous elements. + */ + template + void + copy_from_mg_add (const MGDoFHandler& mg_dof, + OutVector &dst, + const MGLevelObject > &src) const; + + private: + + /** + * Selected component. + */ + unsigned int selected; +}; + + + +#endif diff --git a/deal.II/deal.II/include/multigrid/multigrid.h b/deal.II/deal.II/include/multigrid/multigrid.h index 02acbd1cf2..babcf2cc3c 100644 --- a/deal.II/deal.II/include/multigrid/multigrid.h +++ b/deal.II/deal.II/include/multigrid/multigrid.h @@ -27,33 +27,57 @@ /** - * 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 -class Multigrid : public MGBase +template +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 @@ -63,189 +87,439 @@ class Multigrid : public MGBase * therefore create objects of * this type as late as possible. */ -//TODO:[GK] meaning of minlevel, maxlevel? - Multigrid(const MGDoFHandler &mg_dof_handler, - const ConstraintMatrix &constraints, - const MGLevelObject &level_sparsities, - const MGLevelObject > &level_matrices, - const MGTransferBase &transfer, - const unsigned int minlevel = 0, - const unsigned int maxlevel = 1000000); + template + Multigrid(const MGDoFHandler& 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 - void copy_to_mg (const Vector &src); + template + void vcycle(const MGLevelObject& 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}. In order to - * keep the result consistent, - * constrained degrees of freedom - * are set to zero. + * Print a level vector using + * @ref{DoFHandler}. */ - template - void copy_from_mg (Vector &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 + void level_mgstep (const unsigned int level, + const MGLevelObject& 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 &dest, - const Vector &src, - const Vector &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& 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& get_matrix (const unsigned int level) const; + MGLevelObject defect; - private: /** - * Associated @p{MGDoFHandler}. + * Auxiliary vector. Contains the + * updated solution after the + * multigrid step. */ - SmartPointer > mg_dof_handler; - - -/** - * Sparsity patterns for each level. - */ - SmartPointer > level_sparsities; + MGLevelObject solution; /** - * Matrices for each level. The - * matrices are prepared by the - * constructor of @p{Multigrid} and - * can be accessed for - * assembling. + * Auxiliary vector. */ - SmartPointer > > level_matrices; + MGLevelObject t; + + + private: + /** - * Pointer to the object holding - * constraints. + * Exception. */ - SmartPointer constraints; + DeclException2(ExcSwitchedLevels, int, int, + << "minlevel and maxlevel switched, should be: " + << arg1 << "<=" << arg2); + template 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 +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& mg_dof, + Multigrid& mg, + const MGLevelObject& 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 - void build_matrices (const MGDoFHandler &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 &dst, - const Vector &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 &dst, - const Vector &src) const; - + * Not implemented, but the + * definition may be needed. + */ + void Tvmult_add (VECTOR &dst, + const VECTOR &src) const; + private: + /** + * Associated @p{MGDoFHandler}. + */ + SmartPointer > mg_dof_handler; + + /** + * The multigrid object. + */ + SmartPointer > multigrid; + + /** + * Matrices for each level. The + * matrices are prepared by the + * constructor of @p{Multigrid} and + * can be accessed for + * assembling. + */ + SmartPointer > level_matrices; - std::vector prolongation_sparsities; + /** + * Object for grid tranfer. + */ + SmartPointer transfer; + + /** + * The pre-smoothing object. + */ + SmartPointer 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 post; + + /** + * The coarse grid solver. */ - std::vector > prolongation_matrices; + SmartPointer coarse; }; + + /* --------------------------- inline functions --------------------- */ -template -inline -const SparseMatrix& -Multigrid::get_matrix (const unsigned int level) const +template +template +Multigrid::Multigrid (const MGDoFHandler& 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 +template +void +Multigrid::level_mgstep(const unsigned int level, + const MGLevelObject& matrix, + const TRANSFER& transfer, + const SMOOTHER &pre_smooth, + const SMOOTHER &post_smooth, + const COARSE &coarse_grid_solver) { - Assert((level>=minlevel) && (levelminlevel;--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 +template +void +Multigrid::vcycle(const MGLevelObject& 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 +PreconditionMG +::PreconditionMG(const MGDoFHandler& mg_dof_handler, + Multigrid& mg, + const MGLevelObject& 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 +void +PreconditionMG::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 +void +PreconditionMG::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 +void +PreconditionMG::Tvmult ( + VECTOR&, + const VECTOR&) const +{ + Assert(false, ExcNotImplemented()); +} + + +template +void +PreconditionMG::Tvmult_add ( + VECTOR&, + const VECTOR&) const +{ + Assert(false, ExcNotImplemented()); } diff --git a/deal.II/deal.II/source/dofs/dof_renumbering.cc b/deal.II/deal.II/source/dofs/dof_renumbering.cc index 33484dea77..ccf50ff2bb 100644 --- a/deal.II/deal.II/source/dofs/dof_renumbering.cc +++ b/deal.II/deal.II/source/dofs/dof_renumbering.cc @@ -276,7 +276,7 @@ void DoFRenumbering::Cuthill_McKee (MGDoFHandler &dof_handler // 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 @@ -531,6 +531,103 @@ void DoFRenumbering::component_wise (DoFHandler &dof_handle +template +void DoFRenumbering::component_wise (MGDoFHandler& dof_handler, + unsigned int level, + const std::vector &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 component_order (component_order_arg); + if (component_order.size() == 0) + for (unsigned int i=0; i local_dof_indices(dofs_per_cell); + // prebuilt list to which component + // a given dof on a cell belongs + std::vector component_list (dofs_per_cell); + for (unsigned int i=0; i > component_to_dof_map (dof_handler.get_fe().n_components()); + for (typename MGDoFHandler::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 new_indices (dof_handler.n_dofs(level), + DoFHandler::invalid_dof_index); + for (unsigned int component=0; component::const_iterator + begin_of_component = component_to_dof_map[component].begin(), + end_of_component = component_to_dof_map[component].end(); + + for (typename std::vector::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 void DoFRenumbering::sort_selected_dofs_back (DoFHandler &dof_handler, @@ -775,6 +872,11 @@ template void DoFRenumbering::component_wise (DoFHandler&, const std::vector&); +template +void DoFRenumbering::component_wise (MGDoFHandler&, + unsigned int, + const std::vector&); + template void DoFRenumbering::cell_wise_dg (DoFHandler&, diff --git a/deal.II/deal.II/source/multigrid/mg_base.cc b/deal.II/deal.II/source/multigrid/mg_base.cc index 8b50079c18..aa8d42fc47 100644 --- a/deal.II/deal.II/source/multigrid/mg_base.cc +++ b/deal.II/deal.II/source/multigrid/mg_base.cc @@ -14,139 +14,22 @@ #include #include +#include +#include #include #include -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&, - const Vector&) -{} - ////////////////////////////////////////////////////////////////////// -MGCoarseGridSolver::~MGCoarseGridSolver() -{}; - - -////////////////////////////////////////////////////////////////////// +// template +// MGTransferBase::~MGTransferBase() +// {}; -MGTransferBase::~MGTransferBase() -{}; +// // Explicit instantiations +// template class MGTransferBase >; +// template class MGTransferBase >; +// template class MGTransferBase >; +// template class MGTransferBase >; diff --git a/deal.II/deal.II/source/multigrid/mg_dof_tools.cc b/deal.II/deal.II/source/multigrid/mg_dof_tools.cc index 671bb3fb30..e84883e623 100644 --- a/deal.II/deal.II/source/multigrid/mg_dof_tools.cc +++ b/deal.II/deal.II/source/multigrid/mg_dof_tools.cc @@ -12,17 +12,29 @@ //---------------------------- mg_dof_tools.cc --------------------------- +#include +#include #include -#include -#include +#include +#include +#include +#include #include #include +#include +#include #include +#include +#include #include +//#include -template -void MGDoFTools::make_sparsity_pattern (const MGDoFHandler &dof, +#include +#include + +template +void MGTools::make_sparsity_pattern (const MGDoFHandler &dof, SparsityPattern &sparsity, const unsigned int level) { @@ -50,9 +62,9 @@ void MGDoFTools::make_sparsity_pattern (const MGDoFHandler &dof, -template +template void -MGDoFTools::make_flux_sparsity_pattern (const MGDoFHandler &dof, +MGTools::make_flux_sparsity_pattern (const MGDoFHandler &dof, SparsityPattern &sparsity, const unsigned int level) { @@ -106,9 +118,9 @@ MGDoFTools::make_flux_sparsity_pattern (const MGDoFHandler &dof, -template +template void -MGDoFTools::make_flux_sparsity_pattern_edge (const MGDoFHandler &dof, +MGTools::make_flux_sparsity_pattern_edge (const MGDoFHandler &dof, SparsityPattern &sparsity, const unsigned int level) { @@ -164,9 +176,10 @@ MGDoFTools::make_flux_sparsity_pattern_edge (const MGDoFHandler &dof, -template +//TODO: Replace FullMatrix by vector2d +template void -MGDoFTools::make_flux_sparsity_pattern (const MGDoFHandler &dof, +MGTools::make_flux_sparsity_pattern (const MGDoFHandler &dof, SparsityPattern &sparsity, const unsigned int level, const FullMatrix& int_mask, @@ -262,23 +275,281 @@ MGDoFTools::make_flux_sparsity_pattern (const MGDoFHandler &dof, +template +void +MGTools::count_dofs_per_component (const MGDoFHandler& dof_handler, + std::vector >& 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 &, - SparsityPattern &, - const unsigned int); -template void MGDoFTools::make_flux_sparsity_pattern (const MGDoFHandler &, - SparsityPattern &, - const unsigned int); + const unsigned int n_components = dof_handler.get_fe().n_components(); + for (unsigned int l=0;l > + dofs_in_component (n_components, + std::vector(dof_handler.n_dofs(l), + false)); + std::vector > + component_select (n_components, + std::vector(n_components, false)); + Threads::ThreadManager thread_manager; + for (unsigned int i=0; i &, + const std::vector &, + std::vector &) + = &DoFTools::template extract_level_dofs; + 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 +void +MGTools::reinit_vector (const MGDoFHandler& mg_dof, + MGLevelObject >& 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 +void +MGTools::reinit_vector (const MGDoFHandler& mg_dof, + MGLevelObject >& v, + const std::vector& 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 &, - SparsityPattern &, - const unsigned int); + unsigned int n_selected = std::accumulate (selected.begin(), + selected.end(), + 0U); + + std::vector > + 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 &, + SparsityPattern &, + const unsigned int); + +template void +MGTools::make_flux_sparsity_pattern (const MGDoFHandler &, + SparsityPattern &, + const unsigned int); + +template void +MGTools::make_flux_sparsity_pattern_edge (const MGDoFHandler &, + SparsityPattern &, + const unsigned int); + +template void +MGTools::make_sparsity_pattern (const MGDoFHandler &, + CompressedSparsityPattern &, + const unsigned int); + +template void +MGTools::make_flux_sparsity_pattern (const MGDoFHandler &, + CompressedSparsityPattern &, + const unsigned int); + +template void +MGTools::make_flux_sparsity_pattern_edge (const MGDoFHandler &, + CompressedSparsityPattern &, + const unsigned int); + +template void +MGTools::make_sparsity_pattern (const MGDoFHandler &, + BlockSparsityPattern &, + const unsigned int); +template void +MGTools::make_flux_sparsity_pattern (const MGDoFHandler &, + BlockSparsityPattern &, + const unsigned int); + +template void +MGTools::make_flux_sparsity_pattern_edge (const MGDoFHandler &, + BlockSparsityPattern &, + const unsigned int); + +template void +MGTools::make_sparsity_pattern (const MGDoFHandler &, + CompressedBlockSparsityPattern &, + const unsigned int); + +template void +MGTools::make_flux_sparsity_pattern (const MGDoFHandler &, + CompressedBlockSparsityPattern &, + const unsigned int); + +template void +MGTools::make_flux_sparsity_pattern_edge (const MGDoFHandler &, + CompressedBlockSparsityPattern &, + const unsigned int); #if deal_II_dimension > 1 -template void MGDoFTools::make_flux_sparsity_pattern (const MGDoFHandler &, - SparsityPattern &, - const unsigned int, - const FullMatrix&, - const FullMatrix&); +template void +MGTools::make_flux_sparsity_pattern (const MGDoFHandler &, + SparsityPattern &, + const unsigned int, + const FullMatrix&, + const FullMatrix&); + +template void +MGTools::make_flux_sparsity_pattern (const MGDoFHandler &, + CompressedSparsityPattern &, + const unsigned int, + const FullMatrix&, + const FullMatrix&); + +template void +MGTools::make_flux_sparsity_pattern (const MGDoFHandler &, + BlockSparsityPattern &, + const unsigned int, + const FullMatrix&, + const FullMatrix&); + +template void +MGTools::make_flux_sparsity_pattern (const MGDoFHandler &, + CompressedBlockSparsityPattern &, + const unsigned int, + const FullMatrix&, + const FullMatrix&); + #endif + +template void MGTools::reinit_vector (const MGDoFHandler&, + MGLevelObject >&); +template void MGTools::reinit_vector (const MGDoFHandler&, + MGLevelObject >&); + + +// template void MGTools::copy_to_mg (const MGDoFHandler&, +// const MGTransferBase >&, +// MGLevelObject >&, +// const Vector&); + +// template void MGTools::copy_from_mg(const MGDoFHandler&, +// Vector&, +// const MGLevelObject >&); + +// template void MGTools::copy_from_mg_add(const MGDoFHandler&, +// Vector&, +// const MGLevelObject >&); + +// template void MGTools::copy_to_mg (const MGDoFHandler&, +// const MGTransferBase >&, +// MGLevelObject >&, +// const Vector&); + +// template void MGTools::copy_from_mg(const MGDoFHandler&, +// Vector&, +// const MGLevelObject >&); + +// template void MGTools::copy_from_mg_add(const MGDoFHandler&, +// Vector&, +// const MGLevelObject >&); + +// template void MGTools::copy_to_mg (const MGDoFHandler&, +// const MGTransferBase >&, +// MGLevelObject >&, +// const BlockVector&); + +// template void MGTools::copy_from_mg(const MGDoFHandler&, +// BlockVector&, +// const MGLevelObject >&); + +// template void MGTools::copy_from_mg_add(const MGDoFHandler&, +// BlockVector&, +// const MGLevelObject >&); + +// template void MGTools::copy_to_mg (const MGDoFHandler&, +// const MGTransferBase >&, +// MGLevelObject >&, +// const BlockVector&); + +// template void MGTools::copy_from_mg(const MGDoFHandler&, +// BlockVector&, +// const MGLevelObject >&); + +// template void MGTools::copy_from_mg_add(const MGDoFHandler&, +// BlockVector&, +// const MGLevelObject >&); + +template void MGTools::count_dofs_per_component (const MGDoFHandler&, + std::vector >&); + diff --git a/deal.II/deal.II/source/multigrid/mg_smoother.cc b/deal.II/deal.II/source/multigrid/mg_smoother.cc index 9fd6e98053..fb505e35e1 100644 --- a/deal.II/deal.II/source/multigrid/mg_smoother.cc +++ b/deal.II/deal.II/source/multigrid/mg_smoother.cc @@ -18,10 +18,13 @@ #include #include #include -#include +//#include #include #include #include +#include +#include +#include #include @@ -29,23 +32,11 @@ ////////////////////////////////////////////////////////////////////// -MGSmootherBase::~MGSmootherBase() -{}; - - -////////////////////////////////////////////////////////////////////// - - -void -MGSmootherIdentity::smooth (const unsigned int, - Vector &, - const Vector &) const -{}; - #if deal_II_dimension == 1 -MGSmoother::MGSmoother (const MGDoFHandler<1> &/*mg_dof*/, unsigned int steps) +MGSmootherContinuous::MGSmootherContinuous (const MGDoFHandler<1> &/*mg_dof*/, + unsigned int steps) : steps(steps) { @@ -58,7 +49,8 @@ MGSmoother::MGSmoother (const MGDoFHandler<1> &/*mg_dof*/, unsigned int steps) #if deal_II_dimension > 1 template -MGSmoother::MGSmoother (const MGDoFHandler &mg_dof, unsigned int steps) +MGSmootherContinuous::MGSmootherContinuous (const MGDoFHandler &mg_dof, + unsigned int steps) : steps(steps) { @@ -121,9 +113,10 @@ MGSmoother::MGSmoother (const MGDoFHandler &mg_dof, unsigned int steps) ////////////////////////////////////////////////////////////////////// +template void -MGSmoother::set_zero_interior_boundary (const unsigned int level, - Vector &u) const +MGSmootherContinuous::set_zero_interior_boundary (const unsigned int level, + VECTOR& u) const { if (level==0) return; @@ -140,22 +133,89 @@ MGSmoother::set_zero_interior_boundary (const unsigned int level, // don't do the following instantiation in 1d, since there is a specialized // function there #if deal_II_dimension > 1 -template MGSmoother::MGSmoother (const MGDoFHandler&, unsigned int); +template MGSmootherContinuous::MGSmootherContinuous (const MGDoFHandler&, unsigned int); #endif template -MGSmootherRelaxation -::MGSmootherRelaxation(const MGDoFHandler &mg_dof, - const MGLevelObject > &matrix, - const function_ptr relaxation, - const unsigned int steps, - const double omega); +void MGSmootherContinuous::set_zero_interior_boundary (const unsigned int, + Vector&) const; + +template +void MGSmootherContinuous::set_zero_interior_boundary (const unsigned int, + Vector&) const; + +template +void MGSmootherContinuous::set_zero_interior_boundary (const unsigned int, + BlockVector&) const; + +template +void MGSmootherContinuous::set_zero_interior_boundary (const unsigned int, + BlockVector&) const; + + + +template +MGSmootherRelaxation, Vector > +::MGSmootherRelaxation(const MGDoFHandler&, + const MGLevelObject >&, + const function_ptr, + const unsigned int, + const double); + +template +MGSmootherRelaxation, Vector > +::MGSmootherRelaxation(const MGDoFHandler&, + const MGLevelObject >&, + const function_ptr, + const unsigned int, + const double); + +template +MGSmootherRelaxation, Vector > +::MGSmootherRelaxation(const MGDoFHandler&, + const MGLevelObject >&, + const function_ptr, + const unsigned int, + const double); + +template +MGSmootherRelaxation, Vector > +::MGSmootherRelaxation(const MGDoFHandler&, + const MGLevelObject >&, + const function_ptr, + const unsigned int, + const double); + + +template +MGSmootherRelaxation, BlockVector > +::MGSmootherRelaxation(const MGDoFHandler&, + const MGLevelObject >&, + const function_ptr, + const unsigned int, + const double); + +template +MGSmootherRelaxation, BlockVector > +::MGSmootherRelaxation(const MGDoFHandler&, + const MGLevelObject >&, + const function_ptr, + const unsigned int, + const double); + +template +MGSmootherRelaxation, BlockVector > +::MGSmootherRelaxation(const MGDoFHandler&, + const MGLevelObject >&, + const function_ptr, + const unsigned int, + const double); template -MGSmootherRelaxation -::MGSmootherRelaxation(const MGDoFHandler &mg_dof, - const MGLevelObject > &matrix, - const function_ptr relaxation, - const unsigned int steps, - const double omega); +MGSmootherRelaxation, BlockVector > +::MGSmootherRelaxation(const MGDoFHandler&, + const MGLevelObject >&, + const function_ptr, + const unsigned int, + const double); diff --git a/deal.II/deal.II/source/multigrid/multigrid.all_dimensions.cc b/deal.II/deal.II/source/multigrid/multigrid.all_dimensions.cc index be32a19fa5..879b8e7c32 100644 --- a/deal.II/deal.II/source/multigrid/multigrid.all_dimensions.cc +++ b/deal.II/deal.II/source/multigrid/multigrid.all_dimensions.cc @@ -13,15 +13,10 @@ -#include -#include -#include -#include -#include -#include -#include -#include #include +#include +#include +#include MGTransferPrebuilt::~MGTransferPrebuilt () @@ -50,3 +45,74 @@ void MGTransferPrebuilt::restrict_and_add (const unsigned int from_level, }; + + + +MGTransferBlock::~MGTransferBlock () +{}; + + +void MGTransferBlock::prolongate (const unsigned int to_level, + BlockVector &dst, + const BlockVector &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 &dst, + const BlockVector &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 &dst, + const Vector &src) const +{ + Assert ((to_level >= 1) && (to_level<=prolongation_matrices.size()), + ExcIndexRange (to_level, 1, prolongation_matrices.size()+1)); + + prolongation_matrices[to_level-1].block(selected, selected) + .vmult (dst, src); +}; + + +void MGTransferSelect::restrict_and_add (const unsigned int from_level, + Vector &dst, + const Vector &src) const +{ + Assert ((from_level >= 1) && (from_level<=prolongation_matrices.size()), + ExcIndexRange (from_level, 1, prolongation_matrices.size()+1)); + + prolongation_matrices[from_level-1].block(selected, selected) + .Tvmult_add (dst, src); +}; diff --git a/deal.II/deal.II/source/multigrid/multigrid.cc b/deal.II/deal.II/source/multigrid/multigrid.cc index 6861fcd100..8eeb927ee4 100644 --- a/deal.II/deal.II/source/multigrid/multigrid.cc +++ b/deal.II/deal.II/source/multigrid/multigrid.cc @@ -1,4 +1,4 @@ -//---------------------------- multigrid.cc --------------------------- +//----------------------------------------------------------------------- // $Id$ // Version: $Name$ // @@ -9,52 +9,20 @@ // to the file deal.II/doc/license.html for the text and // further information on this license. // -//---------------------------- multigrid.cc --------------------------- +//----------------------------------------------------------------------- +#include +#include +#include #include -#include -#include -#include #include +#include #include -#include -#include -#include -#include - - -/* --------------------------------- MG ------------------------------------ */ - -template -Multigrid::Multigrid (const MGDoFHandler &mg_dof_handler, - const ConstraintMatrix &constraints, - const MGLevelObject &level_sparsities, - const MGLevelObject > &level_matrices, - const MGTransferBase &transfer, - const unsigned int minlevel, - const unsigned int maxlevel) - : - MGBase(transfer, 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 -void -Multigrid::level_vmult (const unsigned int level, - Vector &result, - const Vector &u, - const Vector &/* rhs */) -{ - (*level_matrices)[level].vmult(result,u); - result.scale(-1.); -} +#include +#include +#include +#include /* ----------------------------- MGTransferPrebuilt ------------------------ */ @@ -66,18 +34,20 @@ void MGTransferPrebuilt::build_matrices (const MGDoFHandler &mg_dof) const unsigned int n_levels = mg_dof.get_tria().n_levels(); const unsigned int dofs_per_cell = mg_dof.get_fe().dofs_per_cell; + sizes.resize(n_levels); + for (unsigned int l=0;l dof_indices_mother (dofs_per_cell); + std::vector dof_indices_parent (dofs_per_cell); std::vector dof_indices_child (dofs_per_cell); // for each level: first build the sparsity @@ -87,24 +57,23 @@ void MGTransferPrebuilt::build_matrices (const MGDoFHandler &mg_dof) // level which have children for (unsigned int level=0; level::cell_iterator cell=mg_dof.begin(level); cell != mg_dof.end(level); ++cell) if (cell->has_children()) { - cell->get_mg_dof_indices (dof_indices_mother); + cell->get_mg_dof_indices (dof_indices_parent); for (unsigned int child=0; child::children_per_cell; ++child) @@ -119,20 +88,19 @@ void MGTransferPrebuilt::build_matrices (const MGDoFHandler &mg_dof) // 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()); prolongation_matrices[level].reinit (prolongation_sparsities[level]); // now actually build the matrices @@ -140,7 +108,7 @@ void MGTransferPrebuilt::build_matrices (const MGDoFHandler &mg_dof) 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::children_per_cell; ++child) @@ -161,7 +129,7 @@ void MGTransferPrebuilt::build_matrices (const MGDoFHandler &mg_dof) { prolongation_matrices[level].set (dof_indices_child[i], - dof_indices_mother[j], + dof_indices_parent[j], prolongation(i,j)); }; }; @@ -170,21 +138,205 @@ void MGTransferPrebuilt::build_matrices (const MGDoFHandler &mg_dof) }; -// explicit instatations -template class Multigrid; -template -void MGTransferPrebuilt::build_matrices (const MGDoFHandler &mg_dof); +/* ----------------------------- MGTransferBlock ------------------------ */ -template -void Multigrid::copy_to_mg (const Vector& src); -template -void Multigrid::copy_to_mg (const Vector& src); +template +void MGTransferBlockBase::build_matrices (const MGDoFHandler &mg_dof, + std::vector select) +{ + const FiniteElement& 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::copy_from_mg (Vector& 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 dof_indices_parent (dofs_per_cell); + std::vector dof_indices_child (dofs_per_cell); + + // for each level: first build the + // sparsity pattern of the matrices + // and then build the matrices + // themselves. note that we only + // need to take care of cells on + // the coarser level which have + // children + for (unsigned int level=0; level::cell_iterator cell=mg_dof.begin(level); + cell != mg_dof.end(level); ++cell) + if (cell->has_children()) + { + cell->get_mg_dof_indices (dof_indices_parent); + + for (unsigned int child=0; + child::children_per_cell; ++child) + { + // set an alias to the + // prolongation matrix for + // this child + const FullMatrix &prolongation + = mg_dof.get_fe().prolongate(child); + + cell->child(child)->get_mg_dof_indices (dof_indices_child); + + // now tag the entries in the + // matrix which will be used + // for this pair of parent/child + for (unsigned int i=0; i::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::children_per_cell; ++child) + { + // set an alias to the + // prolongation matrix for + // this child + const FullMatrix &prolongation + = mg_dof.get_fe().prolongate(child); + + cell->child(child)->get_mg_dof_indices (dof_indices_child); + + // now set the entries in the + // matrix + for (unsigned int i=0; i +void MGTransferBlock::build_matrices (const MGDoFHandler &mg_dof, + std::vector select) +{ + if (select.size() == 0) + select = std::vector (mg_dof.get_fe().n_components(), true); + + MGTransferBlockBase::build_matrices (mg_dof, select); +} + + +template +void MGTransferSelect::build_matrices (const MGDoFHandler &mg_dof, + unsigned int select) +{ + selected = select; + std::vector s(mg_dof.get_fe().n_components(), false); + s[select] = true; + + MGTransferBlockBase::build_matrices (mg_dof, s); +} + + + +// explicit instatations template -void Multigrid::copy_from_mg (Vector& src) const; +void MGTransferPrebuilt::build_matrices (const MGDoFHandler &mg_dof); +template +void MGTransferBlock::build_matrices (const MGDoFHandler &mg_dof, + std::vector); +template +void MGTransferSelect::build_matrices (const MGDoFHandler &mg_dof, + unsigned int); diff --git a/tests/multigrid/mg1.cc b/tests/multigrid/mg1.cc index 453a58ca1e..fa4c7cf818 100644 --- a/tests/multigrid/mg1.cc +++ b/tests/multigrid/mg1.cc @@ -135,7 +135,7 @@ int main() { 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]);