/*---------------------------- multigrid.h ---------------------------*/
#include <vector>
-#include <base/subscriptor.h>
+#include <base/smartpointer.h>
#include <lac/forward-declarations.h>
+#include <lac/vector.h>
+
+class MGTransferBase;
+
/**
* Basic matrix class for multigrid preconditioning.
MultiGridBase(const MultiGridBase&);
const MultiGridBase& operator=(const MultiGridBase&);
- /** Auxiliary vectors.
+ /**
+ * Auxiliary vector.
*/
vector<Vector<float> > d;
+
+ /**
+ * Auxiliary vector.
+ */
vector<Vector<float> > s;
+
+ /**
+ * Auxiliary vector.
+ */
+ Vector<float> t;
/**
* Highest level of cells.
* Level for coarse grid solution.
*/
unsigned minlevel;
-
+
+ /**
+ * Prolongation and restriction object.
+ */
+ SmartPointer<MGTransferBase> transfer;
+
/** Tranfer from dVector to
* MGVector.
*
virtual void smooth(unsigned level,
Vector<float>& x,
const Vector<float>& b,
- unsigned steps);
+ unsigned steps) = 0;
/**
* The post-smoothing algorithm.
unsigned steps);
/**
- * Apply operator on all
+ * Apply residual operator on all
* cells of a level.
- *
+ * This is implemented in a
+ * derived class.
*/
- virtual void level_vmult(unsigned level,
+ virtual void level_residual(unsigned level,
Vector<float>& dst,
- const Vector<float>& src);
+ const Vector<float>& src,
+ const Vector<float>& rhs) = 0;
/**
* Solve exactly on coarsest grid.
+ * Usually, this function should
+ * be overloaded by a more
+ * sophisticated derived
+ * class. Still, there is a
+ * standard implementation doing
+ * #10 * (n_pre_smooth +
+ * n_post_smooth)#
+ * smoothing steps.
*/
virtual void coarse_grid_solution(unsigned l,
Vector<float>& dst,
/**
* Constructor, subject to change.
*/
- MultiGridBase();
+ MultiGridBase(MGTransferBase& transfer,
+ unsigned maxlevel, unsigned minlevel,
+ unsigned pre_smooth, unsigned post_smooth);
virtual ~MultiGridBase();
};
/*---------------------------- multigrid.h ---------------------------*/
#include <vector>
-#include <base/subscriptor.h>
+#include <base/smartpointer.h>
#include <lac/forward-declarations.h>
+#include <lac/vector.h>
+
+class MGTransferBase;
+
/**
* Basic matrix class for multigrid preconditioning.
MultiGridBase(const MultiGridBase&);
const MultiGridBase& operator=(const MultiGridBase&);
- /** Auxiliary vectors.
+ /**
+ * Auxiliary vector.
*/
vector<Vector<float> > d;
+
+ /**
+ * Auxiliary vector.
+ */
vector<Vector<float> > s;
+
+ /**
+ * Auxiliary vector.
+ */
+ Vector<float> t;
/**
* Highest level of cells.
* Level for coarse grid solution.
*/
unsigned minlevel;
-
+
+ /**
+ * Prolongation and restriction object.
+ */
+ SmartPointer<MGTransferBase> transfer;
+
/** Tranfer from dVector to
* MGVector.
*
virtual void smooth(unsigned level,
Vector<float>& x,
const Vector<float>& b,
- unsigned steps);
+ unsigned steps) = 0;
/**
* The post-smoothing algorithm.
unsigned steps);
/**
- * Apply operator on all
+ * Apply residual operator on all
* cells of a level.
- *
+ * This is implemented in a
+ * derived class.
*/
- virtual void level_vmult(unsigned level,
+ virtual void level_residual(unsigned level,
Vector<float>& dst,
- const Vector<float>& src);
+ const Vector<float>& src,
+ const Vector<float>& rhs) = 0;
/**
* Solve exactly on coarsest grid.
+ * Usually, this function should
+ * be overloaded by a more
+ * sophisticated derived
+ * class. Still, there is a
+ * standard implementation doing
+ * #10 * (n_pre_smooth +
+ * n_post_smooth)#
+ * smoothing steps.
*/
virtual void coarse_grid_solution(unsigned l,
Vector<float>& dst,
/**
* Constructor, subject to change.
*/
- MultiGridBase();
+ MultiGridBase(MGTransferBase& transfer,
+ unsigned maxlevel, unsigned minlevel,
+ unsigned pre_smooth, unsigned post_smooth);
virtual ~MultiGridBase();
};
#include <numerics/multigrid.h>
#include <lac/vector.h>
+MultiGridBase::MultiGridBase(MGTransferBase& transfer,
+ unsigned maxlevel, unsigned minlevel,
+ unsigned pre_smooth, unsigned post_smooth)
+ :
+ maxlevel(maxlevel), minlevel(minlevel),
+ transfer(&transfer),
+ n_pre_smooth(pre_smooth), n_post_smooth(post_smooth)
+{}
+
void
MultiGridBase::level_mgstep(unsigned level)
{
return;
}
+ // smoothing of the residual by modifying s
smooth(level, d[level], s[level], n_pre_smooth);
+ // t = d-As
+ level_residual(level, t, s[level], d[level]);
+
+ // make t rhs of lower level
+ transfer->restrict(level, t, d[level-1]);
+ // do recursion
+ level_mgstep(level-1);
+ // do coarse grid correction
+ transfer->prolongate(level, s[level], s[level-1]);
+
+ // smoothing (modify s again)
post_smooth(level, d[level], s[level], n_post_smooth);
}
-void MultiGridBase::post_smooth(unsigned level,
- Vector<float>& dst, const Vector<float>& src,
- unsigned steps)
+void
+MultiGridBase::post_smooth(unsigned level,
+ Vector<float>& dst, const Vector<float>& src,
+ unsigned steps)
{
smooth(level, dst, src, steps);
}
+void
+MultiGridBase::coarse_grid_solution(unsigned level,
+ Vector<float>& dst,
+ const Vector<float>& src)
+{
+ smooth(level, dst, src, 10 * (n_pre_smooth + n_post_smooth));
+}
+
// ab hier Wolfgang
#include <numerics/multigrid.h>
#include <lac/vector.h>
+MultiGridBase::MultiGridBase(MGTransferBase& transfer,
+ unsigned maxlevel, unsigned minlevel,
+ unsigned pre_smooth, unsigned post_smooth)
+ :
+ maxlevel(maxlevel), minlevel(minlevel),
+ transfer(&transfer),
+ n_pre_smooth(pre_smooth), n_post_smooth(post_smooth)
+{}
+
void
MultiGridBase::level_mgstep(unsigned level)
{
return;
}
+ // smoothing of the residual by modifying s
smooth(level, d[level], s[level], n_pre_smooth);
+ // t = d-As
+ level_residual(level, t, s[level], d[level]);
+
+ // make t rhs of lower level
+ transfer->restrict(level, t, d[level-1]);
+ // do recursion
+ level_mgstep(level-1);
+ // do coarse grid correction
+ transfer->prolongate(level, s[level], s[level-1]);
+
+ // smoothing (modify s again)
post_smooth(level, d[level], s[level], n_post_smooth);
}
-void MultiGridBase::post_smooth(unsigned level,
- Vector<float>& dst, const Vector<float>& src,
- unsigned steps)
+void
+MultiGridBase::post_smooth(unsigned level,
+ Vector<float>& dst, const Vector<float>& src,
+ unsigned steps)
{
smooth(level, dst, src, steps);
}
+void
+MultiGridBase::coarse_grid_solution(unsigned level,
+ Vector<float>& dst,
+ const Vector<float>& src)
+{
+ smooth(level, dst, src, 10 * (n_pre_smooth + n_post_smooth));
+}
+
// ab hier Wolfgang