#include <vector>
+
+
/**
* Implementation of multigrid with matrices.
* While #MGBase# was only an abstract framework for the v-cycle,
* @author Wolfgang Bangerth, Guido Kanschat, 1999
*/
template<int dim>
-class MG
- :
- public MGBase
+class MG : public MGBase
{
public:
- MG(const MGDoFHandler<dim>&,
- ConstraintMatrix& hanging_nodes,
- const MGMatrix<SparseMatrix<double> >&,
- const MGTransferBase& transfer,
- unsigned int minlevel = 0, unsigned int maxlevel = 10000);
-
- /** Transfer from dVector to
- * MGVector.
+ /**
+ * Constructor. Take the object
+ * describing the multilevel
+ * hierarchy of degrees of
+ * freedom, an object describing
+ * constraints to degrees of
+ * freedom on the global grid
+ * (for example for hanging
+ * nodes), a set of matrices
+ * which describe the equation
+ * discretized on each level, and
+ * an object mediating the
+ * transfer of solutions between
+ * different levels.
*
- * This function copies data from a
- * dVector, that is, data on the
- * locally finest level, into the
- * corresponding levels of an
- * MGVector.
+ * This function already
+ * initializes the vectors which
+ * will be used later on in the
+ * course of the
+ * computations. You should
+ * therefore create objects of
+ * this type as late as possible.
+ */
+//TODO: minlevel, maxlevel?
+ MG(const MGDoFHandler<dim> &mg_dof_handler,
+ ConstraintMatrix &constraints,
+ const MGMatrix<SparseMatrix<double> > &level_matrices,
+ const MGTransferBase &transfer,
+ const unsigned int minlevel = 0,
+ const unsigned int maxlevel = 1000000);
+
+ /**
+ * Transfer from a vector on the
+ * global grid to vectors defined
+ * on each of the levels
+ * separately, i.a. an #MGVector#.
*/
template<typename number>
- void copy_to_mg(const Vector<number>& src);
+ void copy_to_mg (const Vector<number> &src);
/**
- * Transfer from multi-grid vector to
+ * Transfer from multi-level vector to
* normal vector.
*
- * Copies data from active portions
- * of an MGVector into the
- * respective positions of a
- * Vector<double>. All other entries of
- * #src# are zero.
+ * Copies data from active
+ * portions of an MGVector into
+ * the respective positions of a
+ * #Vector<double>#. In order to
+ * keep the result consistent,
+ * constrained degrees of freedom are set to zero.
*/
template<typename number>
- void copy_from_mg(Vector<number>& dst) const;
+ void copy_from_mg (Vector<number> &dst) const;
/**
* Negative #vmult# on a level.
+//TODO: what does this function do?
* @see MGBase.
*/
- virtual void level_vmult(unsigned int,
- Vector<double> &,
- const Vector<double> &,
- const Vector<double> &);
+ virtual void level_vmult (const unsigned int level,
+ Vector<double> &dest,
+ const Vector<double> &src,
+ const Vector<double> &rhs);
/**
* Read-only access to level matrices.
*/
- const SparseMatrix<double>& get_matrix(unsigned int level) const;
+ const SparseMatrix<double>& get_matrix (const unsigned int level) const;
private:
/**
* Associated #MGDoFHandler#.
*/
- SmartPointer<const MGDoFHandler<dim> > dofs;
+ SmartPointer<const MGDoFHandler<dim> > mg_dof_handler;
/**
* Matrices for each level.
* the constructor of #MG# and can
* be accessed for assembling.
*/
- SmartPointer<const MGMatrix<SparseMatrix<double> > > matrices;
- ConstraintMatrix& hanging_nodes;
+ SmartPointer<const MGMatrix<SparseMatrix<double> > > level_matrices;
+
+ /**
+ * Pointer to the object holding
+ * constraints.
+ */
+ SmartPointer<const ConstraintMatrix> constraints;
};
+
+
+
/**
* Implementation of the #MGTransferBase# interface for which the transfer
* operations are prebuilt upon construction of the object of this class as
class MGTransferPrebuilt : public MGTransferBase
{
public:
- /**
- * Preliminary constructor
- */
- MGTransferPrebuilt()
- {}
-
-
/**
* Destructor.
*/
/**
* Prolongate a vector from level
- * #to_level-1# to level #to_level#.
+ * #to_level-1# to level
+ * #to_level#. The previous
+ * content of #dst# is
+ * overwritten.
*
* #src# is assumed to be a vector with
* as many elements as there are degrees
* are degrees of freedom on the finer
* level.
*/
- virtual void prolongate (const unsigned int to_level,
+ virtual void prolongate (const unsigned int to_level,
Vector<double> &dst,
const Vector<double> &src) const;
/**
* Restrict a vector from level
* #from_level# to level
- * #from_level-1#.
+ * #from_level-1# and add this
+ * restriction to
+ * #dst#. Obviously, if the
+ * refined region on level
+ * #from_level# is smaller than
+ * that on level #from_level-1#,
+ * some degrees of freedom in
+ * #dst# are not covered and will
+ * not be altered. For the other
+ * degress of freedom, the result
+ * of the restriction is added.
*
* #src# is assumed to be a vector with
* as many elements as there are degrees
* are degrees of freedom on the coarser
* level.
*/
- virtual void restrict (const unsigned int from_level,
- Vector<double> &dst,
- const Vector<double> &src) const;
+ virtual void restrict_and_add (const unsigned int from_level,
+ Vector<double> &dst,
+ const Vector<double> &src) const;
private:
};
+
+
+/* --------------------------- inline functions --------------------- */
+
+
template<int dim>
inline
const SparseMatrix<double>&
MG<dim>::get_matrix(unsigned int level) const
{
- Assert((level>=minlevel) && (level<maxlevel), ExcIndexRange(level, minlevel, maxlevel));
+ Assert((level>=minlevel) && (level<maxlevel),
+ ExcIndexRange(level, minlevel, maxlevel));
- return (*matrices)[level];
+ return (*level_matrices)[level];
}
// $Id$
+#include <grid/dof_constraints.h>
#include <numerics/multigrid.h>
#include <algorithm>
#include <fstream>
-template <int dim>
-MG<dim>::MG(const MGDoFHandler<dim>& dofs,
- ConstraintMatrix& hanging_nodes,
- const MGMatrix<SparseMatrix<double> >& matrices,
- const MGTransferBase& transfer,
- unsigned int minl, unsigned int maxl)
- :
- MGBase(transfer, minl,
- min(dofs.get_tria().n_levels()-1, maxl)),
- dofs(&dofs),
- matrices(&matrices),
- hanging_nodes(hanging_nodes)
-{
- for (unsigned int level = minlevel; level<=maxlevel ; ++level)
- {
- solution[level].reinit(matrices[level].m());
- defect[level].reinit(matrices[level].m());
- };
-};
-
-
-
template <int dim>
template <typename number>
void
MG<dim>::copy_to_mg(const Vector<number>& src)
{
- const unsigned int fe_dofs = dofs->get_fe().total_dofs;
- const unsigned int face_dofs = dofs->get_fe().dofs_per_face;
+ const unsigned int fe_dofs = mg_dof_handler->get_fe().total_dofs;
+ const unsigned int face_dofs = mg_dof_handler->get_fe().dofs_per_face;
// set the elements of the vectors
// on all levels to zero
defect.clear();
-// hanging_nodes.condense(src);
+// constraints->condense(src);
vector<int> index(fe_dofs);
vector<int> mgindex(fe_dofs);
vector<int> mgfaceindex(face_dofs);
-
-// {
-// ofstream of("CT");
-// DataOut<2> out;
-// out.attach_dof_handler(*dofs);
-// out.add_data_vector(src,"solution");
-// out.write_gnuplot(of,1);
-// }
- for (int level = maxlevel; level >= static_cast<int>(minlevel); --level)
+ for (int level=maxlevel; level>=static_cast<int>(minlevel); --level)
{
- DoFHandler<dim>::active_cell_iterator dc = dofs->DoFHandler<dim>::begin_active(level);
- MGDoFHandler<dim>::active_cell_iterator c = dofs->begin_active(level);
-
+ DoFHandler<dim>::active_cell_iterator
+ dc = mg_dof_handler->DoFHandler<dim>::begin_active(level);
+ MGDoFHandler<dim>::active_cell_iterator
+ level_cell = mg_dof_handler->begin_active(level);
+ const MGDoFHandler<dim>::active_cell_iterator
+ level_end = mg_dof_handler->end_active(level);
+
//TODO: Document on difference between #end# and #end_level#
- for (; c != dofs->end_active(level) ; ++c, ++dc)
+ for (; level_cell!=level_end; ++level_cell, ++dc)
{
dc->get_dof_indices(index);
- c->get_mg_dof_indices(mgindex);
- for (unsigned int i=0;i<fe_dofs;++i)
- defect[c->level()](mgindex[i]) = src(index[i]);
+ level_cell->get_mg_dof_indices (mgindex);
+ for (unsigned int i=0; i<fe_dofs; ++i)
+ defect[level](mgindex[i]) = src(index[i]);
+//TODO: what happens here?
// Delete values on refinement edge
- for (unsigned int face_n = 0; face_n< GeometryInfo<dim>::faces_per_cell; ++face_n)
+ for (unsigned int face_n=0; face_n<GeometryInfo<dim>::faces_per_cell; ++face_n)
{
- MGDoFHandler<dim>::face_iterator face = c->face(face_n);
+ const MGDoFHandler<dim>::face_iterator face = level_cell->face(face_n);
if (face->has_children())
{
face->get_mg_dof_indices(mgfaceindex);
- for (unsigned int i=0;i<face_dofs;++i)
- defect[c->level()](mgfaceindex[i]) = 0.;
- }
- }
- }
-// if (level > (int) minlevel)
-// transfer->restrict(level, d[level-1], d[level]);
- if (level < (int) maxlevel)
- transfer->restrict(level+1, defect[level], defect[level+1]);
- }
-
-// for (unsigned int i=minlevel; i<= maxlevel; ++i)
-// {
-// char name[10];
-// sprintf(name,"CT%d",i);
-// ofstream of(name);
-// write_gnuplot(*dofs, d[i], i, of);
-// }
-}
+ for (unsigned int i=0; i<face_dofs; ++i)
+ defect[level](mgfaceindex[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)
+ transfer->restrict_and_add (level+1, defect[level], defect[level+1]);
+ };
+};
+
+
template <int dim>
template <typename number>
void
MG<dim>::copy_from_mg(Vector<number>& dst) const
{
- const unsigned int fe_dofs = dofs->get_fe().total_dofs;
+ const unsigned int fe_dofs = mg_dof_handler->get_fe().total_dofs;
vector<int> index(fe_dofs);
vector<int> mgindex(fe_dofs);
// write_gnuplot(*dofs, s[i], i, of);
// }
- DoFHandler<dim>::active_cell_iterator dc = dofs->DoFHandler<dim>::begin_active();
- for (MGDoFHandler<dim>::active_cell_iterator c = dofs->begin_active();
- c != dofs->end(); ++c, ++dc)
+ DoFHandler<dim>::active_cell_iterator dc = mg_dof_handler->DoFHandler<dim>::begin_active();
+ for (MGDoFHandler<dim>::active_cell_iterator c = mg_dof_handler->begin_active();
+ c != mg_dof_handler->end(); ++c, ++dc)
{
dc->get_dof_indices(index);
c->get_mg_dof_indices(mgindex);
for (unsigned int i=0;i<fe_dofs;++i)
dst(index[i]) = solution[c->level()](mgindex[i]);
}
- hanging_nodes.set_zero(dst);
+ constraints->set_zero(dst);
// {
// ofstream of("CF");
// DataOut<2> out;
// }
}
-
-
-template <int dim>
-void
-MG<dim>::level_vmult(unsigned int l,
- Vector<double> &r,
- const Vector<double> &u,
- const Vector<double> &)
-{
- (*matrices)[l].vmult(r,u);
- r.scale(-1.);
-}
-
-/*
-template <int dim>
-MG::
-template <int dim>
-MG::
-template <int dim>
-MG::
-
-*/
#include <vector>
+
+
/**
* Implementation of multigrid with matrices.
* While #MGBase# was only an abstract framework for the v-cycle,
* @author Wolfgang Bangerth, Guido Kanschat, 1999
*/
template<int dim>
-class MG
- :
- public MGBase
+class MG : public MGBase
{
public:
- MG(const MGDoFHandler<dim>&,
- ConstraintMatrix& hanging_nodes,
- const MGMatrix<SparseMatrix<double> >&,
- const MGTransferBase& transfer,
- unsigned int minlevel = 0, unsigned int maxlevel = 10000);
-
- /** Transfer from dVector to
- * MGVector.
+ /**
+ * Constructor. Take the object
+ * describing the multilevel
+ * hierarchy of degrees of
+ * freedom, an object describing
+ * constraints to degrees of
+ * freedom on the global grid
+ * (for example for hanging
+ * nodes), a set of matrices
+ * which describe the equation
+ * discretized on each level, and
+ * an object mediating the
+ * transfer of solutions between
+ * different levels.
*
- * This function copies data from a
- * dVector, that is, data on the
- * locally finest level, into the
- * corresponding levels of an
- * MGVector.
+ * This function already
+ * initializes the vectors which
+ * will be used later on in the
+ * course of the
+ * computations. You should
+ * therefore create objects of
+ * this type as late as possible.
+ */
+//TODO: minlevel, maxlevel?
+ MG(const MGDoFHandler<dim> &mg_dof_handler,
+ ConstraintMatrix &constraints,
+ const MGMatrix<SparseMatrix<double> > &level_matrices,
+ const MGTransferBase &transfer,
+ const unsigned int minlevel = 0,
+ const unsigned int maxlevel = 1000000);
+
+ /**
+ * Transfer from a vector on the
+ * global grid to vectors defined
+ * on each of the levels
+ * separately, i.a. an #MGVector#.
*/
template<typename number>
- void copy_to_mg(const Vector<number>& src);
+ void copy_to_mg (const Vector<number> &src);
/**
- * Transfer from multi-grid vector to
+ * Transfer from multi-level vector to
* normal vector.
*
- * Copies data from active portions
- * of an MGVector into the
- * respective positions of a
- * Vector<double>. All other entries of
- * #src# are zero.
+ * Copies data from active
+ * portions of an MGVector into
+ * the respective positions of a
+ * #Vector<double>#. In order to
+ * keep the result consistent,
+ * constrained degrees of freedom are set to zero.
*/
template<typename number>
- void copy_from_mg(Vector<number>& dst) const;
+ void copy_from_mg (Vector<number> &dst) const;
/**
* Negative #vmult# on a level.
+//TODO: what does this function do?
* @see MGBase.
*/
- virtual void level_vmult(unsigned int,
- Vector<double> &,
- const Vector<double> &,
- const Vector<double> &);
+ virtual void level_vmult (const unsigned int level,
+ Vector<double> &dest,
+ const Vector<double> &src,
+ const Vector<double> &rhs);
/**
* Read-only access to level matrices.
*/
- const SparseMatrix<double>& get_matrix(unsigned int level) const;
+ const SparseMatrix<double>& get_matrix (const unsigned int level) const;
private:
/**
* Associated #MGDoFHandler#.
*/
- SmartPointer<const MGDoFHandler<dim> > dofs;
+ SmartPointer<const MGDoFHandler<dim> > mg_dof_handler;
/**
* Matrices for each level.
* the constructor of #MG# and can
* be accessed for assembling.
*/
- SmartPointer<const MGMatrix<SparseMatrix<double> > > matrices;
- ConstraintMatrix& hanging_nodes;
+ SmartPointer<const MGMatrix<SparseMatrix<double> > > level_matrices;
+
+ /**
+ * Pointer to the object holding
+ * constraints.
+ */
+ SmartPointer<const ConstraintMatrix> constraints;
};
+
+
+
/**
* Implementation of the #MGTransferBase# interface for which the transfer
* operations are prebuilt upon construction of the object of this class as
class MGTransferPrebuilt : public MGTransferBase
{
public:
- /**
- * Preliminary constructor
- */
- MGTransferPrebuilt()
- {}
-
-
/**
* Destructor.
*/
/**
* Prolongate a vector from level
- * #to_level-1# to level #to_level#.
+ * #to_level-1# to level
+ * #to_level#. The previous
+ * content of #dst# is
+ * overwritten.
*
* #src# is assumed to be a vector with
* as many elements as there are degrees
* are degrees of freedom on the finer
* level.
*/
- virtual void prolongate (const unsigned int to_level,
+ virtual void prolongate (const unsigned int to_level,
Vector<double> &dst,
const Vector<double> &src) const;
/**
* Restrict a vector from level
* #from_level# to level
- * #from_level-1#.
+ * #from_level-1# and add this
+ * restriction to
+ * #dst#. Obviously, if the
+ * refined region on level
+ * #from_level# is smaller than
+ * that on level #from_level-1#,
+ * some degrees of freedom in
+ * #dst# are not covered and will
+ * not be altered. For the other
+ * degress of freedom, the result
+ * of the restriction is added.
*
* #src# is assumed to be a vector with
* as many elements as there are degrees
* are degrees of freedom on the coarser
* level.
*/
- virtual void restrict (const unsigned int from_level,
- Vector<double> &dst,
- const Vector<double> &src) const;
+ virtual void restrict_and_add (const unsigned int from_level,
+ Vector<double> &dst,
+ const Vector<double> &src) const;
private:
};
+
+
+/* --------------------------- inline functions --------------------- */
+
+
template<int dim>
inline
const SparseMatrix<double>&
MG<dim>::get_matrix(unsigned int level) const
{
- Assert((level>=minlevel) && (level<maxlevel), ExcIndexRange(level, minlevel, maxlevel));
+ Assert((level>=minlevel) && (level<maxlevel),
+ ExcIndexRange(level, minlevel, maxlevel));
- return (*matrices)[level];
+ return (*level_matrices)[level];
}
// $Id$
+#include <grid/dof_constraints.h>
#include <numerics/multigrid.h>
#include <algorithm>
#include <fstream>
-template <int dim>
-MG<dim>::MG(const MGDoFHandler<dim>& dofs,
- ConstraintMatrix& hanging_nodes,
- const MGMatrix<SparseMatrix<double> >& matrices,
- const MGTransferBase& transfer,
- unsigned int minl, unsigned int maxl)
- :
- MGBase(transfer, minl,
- min(dofs.get_tria().n_levels()-1, maxl)),
- dofs(&dofs),
- matrices(&matrices),
- hanging_nodes(hanging_nodes)
-{
- for (unsigned int level = minlevel; level<=maxlevel ; ++level)
- {
- solution[level].reinit(matrices[level].m());
- defect[level].reinit(matrices[level].m());
- };
-};
-
-
-
template <int dim>
template <typename number>
void
MG<dim>::copy_to_mg(const Vector<number>& src)
{
- const unsigned int fe_dofs = dofs->get_fe().total_dofs;
- const unsigned int face_dofs = dofs->get_fe().dofs_per_face;
+ const unsigned int fe_dofs = mg_dof_handler->get_fe().total_dofs;
+ const unsigned int face_dofs = mg_dof_handler->get_fe().dofs_per_face;
// set the elements of the vectors
// on all levels to zero
defect.clear();
-// hanging_nodes.condense(src);
+// constraints->condense(src);
vector<int> index(fe_dofs);
vector<int> mgindex(fe_dofs);
vector<int> mgfaceindex(face_dofs);
-
-// {
-// ofstream of("CT");
-// DataOut<2> out;
-// out.attach_dof_handler(*dofs);
-// out.add_data_vector(src,"solution");
-// out.write_gnuplot(of,1);
-// }
- for (int level = maxlevel; level >= static_cast<int>(minlevel); --level)
+ for (int level=maxlevel; level>=static_cast<int>(minlevel); --level)
{
- DoFHandler<dim>::active_cell_iterator dc = dofs->DoFHandler<dim>::begin_active(level);
- MGDoFHandler<dim>::active_cell_iterator c = dofs->begin_active(level);
-
+ DoFHandler<dim>::active_cell_iterator
+ dc = mg_dof_handler->DoFHandler<dim>::begin_active(level);
+ MGDoFHandler<dim>::active_cell_iterator
+ level_cell = mg_dof_handler->begin_active(level);
+ const MGDoFHandler<dim>::active_cell_iterator
+ level_end = mg_dof_handler->end_active(level);
+
//TODO: Document on difference between #end# and #end_level#
- for (; c != dofs->end_active(level) ; ++c, ++dc)
+ for (; level_cell!=level_end; ++level_cell, ++dc)
{
dc->get_dof_indices(index);
- c->get_mg_dof_indices(mgindex);
- for (unsigned int i=0;i<fe_dofs;++i)
- defect[c->level()](mgindex[i]) = src(index[i]);
+ level_cell->get_mg_dof_indices (mgindex);
+ for (unsigned int i=0; i<fe_dofs; ++i)
+ defect[level](mgindex[i]) = src(index[i]);
+//TODO: what happens here?
// Delete values on refinement edge
- for (unsigned int face_n = 0; face_n< GeometryInfo<dim>::faces_per_cell; ++face_n)
+ for (unsigned int face_n=0; face_n<GeometryInfo<dim>::faces_per_cell; ++face_n)
{
- MGDoFHandler<dim>::face_iterator face = c->face(face_n);
+ const MGDoFHandler<dim>::face_iterator face = level_cell->face(face_n);
if (face->has_children())
{
face->get_mg_dof_indices(mgfaceindex);
- for (unsigned int i=0;i<face_dofs;++i)
- defect[c->level()](mgfaceindex[i]) = 0.;
- }
- }
- }
-// if (level > (int) minlevel)
-// transfer->restrict(level, d[level-1], d[level]);
- if (level < (int) maxlevel)
- transfer->restrict(level+1, defect[level], defect[level+1]);
- }
-
-// for (unsigned int i=minlevel; i<= maxlevel; ++i)
-// {
-// char name[10];
-// sprintf(name,"CT%d",i);
-// ofstream of(name);
-// write_gnuplot(*dofs, d[i], i, of);
-// }
-}
+ for (unsigned int i=0; i<face_dofs; ++i)
+ defect[level](mgfaceindex[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)
+ transfer->restrict_and_add (level+1, defect[level], defect[level+1]);
+ };
+};
+
+
template <int dim>
template <typename number>
void
MG<dim>::copy_from_mg(Vector<number>& dst) const
{
- const unsigned int fe_dofs = dofs->get_fe().total_dofs;
+ const unsigned int fe_dofs = mg_dof_handler->get_fe().total_dofs;
vector<int> index(fe_dofs);
vector<int> mgindex(fe_dofs);
// write_gnuplot(*dofs, s[i], i, of);
// }
- DoFHandler<dim>::active_cell_iterator dc = dofs->DoFHandler<dim>::begin_active();
- for (MGDoFHandler<dim>::active_cell_iterator c = dofs->begin_active();
- c != dofs->end(); ++c, ++dc)
+ DoFHandler<dim>::active_cell_iterator dc = mg_dof_handler->DoFHandler<dim>::begin_active();
+ for (MGDoFHandler<dim>::active_cell_iterator c = mg_dof_handler->begin_active();
+ c != mg_dof_handler->end(); ++c, ++dc)
{
dc->get_dof_indices(index);
c->get_mg_dof_indices(mgindex);
for (unsigned int i=0;i<fe_dofs;++i)
dst(index[i]) = solution[c->level()](mgindex[i]);
}
- hanging_nodes.set_zero(dst);
+ constraints->set_zero(dst);
// {
// ofstream of("CF");
// DataOut<2> out;
// }
}
-
-
-template <int dim>
-void
-MG<dim>::level_vmult(unsigned int l,
- Vector<double> &r,
- const Vector<double> &u,
- const Vector<double> &)
-{
- (*matrices)[l].vmult(r,u);
- r.scale(-1.);
-}
-
-/*
-template <int dim>
-MG::
-template <int dim>
-MG::
-template <int dim>
-MG::
-
-*/
#include <numerics/mg_smoother.h>
#include <lac/vector.h>
-extern void write_gnuplot (const MGDoFHandler<2>& dofs, const Vector<double>& v, unsigned int level,
- ostream &out, unsigned int accuracy);
+
+
+/* --------------------------------- MG ------------------------------------ */
+
+template <int dim>
+MG<dim>::MG(const MGDoFHandler<dim> &mg_dof_handler,
+ ConstraintMatrix &constraints,
+ const MGMatrix<SparseMatrix<double> > &level_matrices,
+ const MGTransferBase &transfer,
+ unsigned int minl, unsigned int maxl)
+ :
+ MGBase(transfer, minl,
+ min(mg_dof_handler.get_tria().n_levels()-1, maxl)),
+ mg_dof_handler(&mg_dof_handler),
+ level_matrices(&level_matrices),
+ constraints(&constraints)
+{
+ // initialize the objects with
+ // their correct size
+ for (unsigned int level=minlevel; level<=maxlevel ; ++level)
+ {
+ solution[level].reinit(level_matrices[level].m());
+ defect[level].reinit(level_matrices[level].m());
+ };
+};
+
+
+
+template <int dim>
+void
+MG<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.);
+}
+
+
+
+
/* ----------------------------- MGTransferPrebuilt ------------------------ */
// cell
prolongation_sparsities.back().reinit (mg_dof.n_dofs(level+1),
mg_dof.n_dofs(level),
-// evil hack, must be corrected!!!
+//TODO: evil hack, must be corrected!!!
dofs_per_cell+1);
for (typename MGDoFHandler<dim>::cell_iterator cell=mg_dof.begin(level);
-void MGTransferPrebuilt::restrict (const unsigned int from_level,
- Vector<double> &dst,
- const Vector<double> &src) const
+void MGTransferPrebuilt::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));
// explicit instatations
+template class MG<deal_II_dimension>;
template
void MGTransferPrebuilt::build_matrices (const MGDoFHandler<deal_II_dimension> &mg_dof);
-template MG<deal_II_dimension>;
#include <numerics/mg_smoother.h>
#include <lac/vector.h>
-extern void write_gnuplot (const MGDoFHandler<2>& dofs, const Vector<double>& v, unsigned int level,
- ostream &out, unsigned int accuracy);
+
+
+/* --------------------------------- MG ------------------------------------ */
+
+template <int dim>
+MG<dim>::MG(const MGDoFHandler<dim> &mg_dof_handler,
+ ConstraintMatrix &constraints,
+ const MGMatrix<SparseMatrix<double> > &level_matrices,
+ const MGTransferBase &transfer,
+ unsigned int minl, unsigned int maxl)
+ :
+ MGBase(transfer, minl,
+ min(mg_dof_handler.get_tria().n_levels()-1, maxl)),
+ mg_dof_handler(&mg_dof_handler),
+ level_matrices(&level_matrices),
+ constraints(&constraints)
+{
+ // initialize the objects with
+ // their correct size
+ for (unsigned int level=minlevel; level<=maxlevel ; ++level)
+ {
+ solution[level].reinit(level_matrices[level].m());
+ defect[level].reinit(level_matrices[level].m());
+ };
+};
+
+
+
+template <int dim>
+void
+MG<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.);
+}
+
+
+
+
/* ----------------------------- MGTransferPrebuilt ------------------------ */
// cell
prolongation_sparsities.back().reinit (mg_dof.n_dofs(level+1),
mg_dof.n_dofs(level),
-// evil hack, must be corrected!!!
+//TODO: evil hack, must be corrected!!!
dofs_per_cell+1);
for (typename MGDoFHandler<dim>::cell_iterator cell=mg_dof.begin(level);
-void MGTransferPrebuilt::restrict (const unsigned int from_level,
- Vector<double> &dst,
- const Vector<double> &src) const
+void MGTransferPrebuilt::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));
// explicit instatations
+template class MG<deal_II_dimension>;
template
void MGTransferPrebuilt::build_matrices (const MGDoFHandler<deal_II_dimension> &mg_dof);
-template MG<deal_II_dimension>;