*/
virtual void vmult(unsigned int level, VECTOR& dst,
const VECTOR& src) const = 0;
+
+ /**
+ * Adding matrix-vector-multiplication on
+ * a certain level.
+ */
+ virtual void vmult_add(unsigned int level, VECTOR& dst,
+ const VECTOR& src) const = 0;
+
+ /**
+ * Transpose
+ * matrix-vector-multiplication on
+ * a certain level.
+ */
+ virtual void Tvmult(unsigned int level, VECTOR& dst,
+ const VECTOR& src) const = 0;
+
+ /**
+ * Adding transpose
+ * matrix-vector-multiplication on
+ * a certain level.
+ */
+ virtual void Tvmult_add(unsigned int level, VECTOR& dst,
+ const VECTOR& src) const = 0;
};
virtual void vmult(unsigned int level, VECTOR& dst,
const VECTOR& src) const;
+ /**
+ * Adding matrix-vector-multiplication on
+ * a certain level.
+ */
+ virtual void vmult_add(unsigned int level, VECTOR& dst,
+ const VECTOR& src) const;
+
+ /**
+ * Transpose
+ * matrix-vector-multiplication on
+ * a certain level.
+ */
+ virtual void Tvmult(unsigned int level, VECTOR& dst,
+ const VECTOR& src) const;
+
+ /**
+ * Adding transpose
+ * matrix-vector-multiplication on
+ * a certain level.
+ */
+ virtual void Tvmult_add(unsigned int level, VECTOR& dst,
+ const VECTOR& src) const;
};
}
+template <class MATRIX, class VECTOR>
+void
+MGMatrix<MATRIX, VECTOR>::vmult_add (unsigned int level,
+ VECTOR& dst,
+ const VECTOR& src) const
+{
+ const MGLevelObject<MATRIX>& m = **this;
+ m[level].vmult_add(dst, src);
+}
+
+
+template <class MATRIX, class VECTOR>
+void
+MGMatrix<MATRIX, VECTOR>::Tvmult (unsigned int level,
+ VECTOR& dst,
+ const VECTOR& src) const
+{
+ const MGLevelObject<MATRIX>& m = **this;
+ m[level].Tvmult(dst, src);
+}
+
+
+template <class MATRIX, class VECTOR>
+void
+MGMatrix<MATRIX, VECTOR>::Tvmult_add (unsigned int level,
+ VECTOR& dst,
+ const VECTOR& src) const
+{
+ const MGLevelObject<MATRIX>& m = **this;
+ m[level].Tvmult_add(dst, src);
+}
+
+
/*---------------------------- mgbase.h ---------------------------*/
#endif
* Make a sparsity pattern including fluxes
* of discontinuous Galerkin methods.
* @ref make_sparsity_pattern
- * $see DoFTools
+ * @ref DoFTools
*/
template <int dim, class SparsityPattern>
static void
/**
* Create sparsity pattern for
* the fluxes at refinement
- * edges.
+ * edges. The matrix maps a
+ * function of the fine level
+ * space @p{level} to the coarser
+ * space.
+ *
* @ref{make_flux_sparsity_pattern}
- * @ref{DoFTools}
*/
template <int dim, class SparsityPattern>
static void
#include <dofs/dof_constraints.h>
#include <numerics/data_out.h>
-#include <multigrid/multigrid.h>
+#include <multigrid/mg_base.h>
#include <algorithm>
#include <fstream>
*/
SmartPointer<const MGSmoother<VECTOR> > post_smooth;
-
+ /**
+ * Edge matrix from fine to coarse.
+ */
+ SmartPointer<const MGMatrixBase<VECTOR> > edge_down;
+
+ /**
+ * Transpose edge matrix from coarse to fine.
+ */
+ SmartPointer<const MGMatrixBase<VECTOR> > edge_up;
/**
* Exception.
};
-
-
/* --------------------------- inline functions --------------------- */
-template <class VECTOR>
-template <int dim>
-Multigrid<VECTOR>::Multigrid (const MGDoFHandler<dim>& mg_dof_handler,
- const MGMatrixBase<VECTOR>& matrix,
- const MGCoarseGrid<VECTOR>& coarse,
- const MGTransfer<VECTOR>& transfer,
- const MGSmoother<VECTOR>& pre_smooth,
- const MGSmoother<VECTOR>& post_smooth,
- 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),
- matrix(&matrix),
- coarse(&coarse),
- transfer(&transfer),
- pre_smooth(&pre_smooth),
- post_smooth(&post_smooth)
-{
-};
-
-
-template <class VECTOR>
-void
-Multigrid<VECTOR>::level_mgstep(const unsigned int level)
-{
-#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)(level, solution[level], defect[level]);
-#ifdef MG_DEBUG
- sprintf(name, "MG%d-solution",level);
- print_vector(level, solution[level], name);
-#endif
- return;
- }
-
- // smoothing of the residual by modifying s
- pre_smooth->smooth(level, solution[level], defect[level]);
-
-#ifdef MG_DEBUG
- sprintf(name, "MG%d-pre",level);
- print_vector(level, solution[level], name);
-#endif
-
- // t = -A*solution[level]
- matrix->vmult(level, t[level], solution[level]);
- t[level].scale(-1);
-
- // make t rhs of lower level
- // The non-refined parts of the
- // coarse-level defect already contain
- // the global defect, the refined parts
- // its restriction.
- for (unsigned int l = level;l>minlevel;--l)
- {
- t[l-1] = 0.;
- transfer->restrict_and_add (l, t[l-1], t[l]);
- defect[l-1] += t[l-1];
- }
-
- // add additional DG contribution
-// edge_vmult(level, defect[level-1], defect[level]);
-
- // do recursion
- level_mgstep(level-1);
-
- // 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
-}
-
-
-template <class VECTOR>
-void
-Multigrid<VECTOR>::vcycle()
-{
- // 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);
-// abort ();
-}
-
-
-/* ------------------------------------------------------------------- */
-
-
template<int dim, class VECTOR, class TRANSFER>
PreconditionMG<dim, VECTOR, TRANSFER>
::PreconditionMG(const MGDoFHandler<dim>& mg_dof_handler,
#define __deal2__multigrid_templates_h
-#include <dofs/dof_constraints.h>
-//#include <numerics/data_out.h>
#include <multigrid/multigrid.h>
-#include <algorithm>
-#include <fstream>
-#include <lac/sparse_matrix.h>
+/* --------------------------- inline functions --------------------- */
+template <class VECTOR>
template <int dim>
+Multigrid<VECTOR>::Multigrid (const MGDoFHandler<dim>& mg_dof_handler,
+ const MGMatrixBase<VECTOR>& matrix,
+ const MGCoarseGrid<VECTOR>& coarse,
+ const MGTransfer<VECTOR>& transfer,
+ const MGSmoother<VECTOR>& pre_smooth,
+ const MGSmoother<VECTOR>& post_smooth,
+ 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),
+ matrix(&matrix),
+ coarse(&coarse),
+ transfer(&transfer),
+ pre_smooth(&pre_smooth),
+ post_smooth(&post_smooth)
+{
+};
+
+
+template <class VECTOR>
+void
+Multigrid<VECTOR>::level_mgstep(const unsigned int level)
+{
+#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)(level, solution[level], defect[level]);
+#ifdef MG_DEBUG
+ sprintf(name, "MG%d-solution",level);
+ print_vector(level, solution[level], name);
+#endif
+ return;
+ }
+
+ // smoothing of the residual by modifying s
+ pre_smooth->smooth(level, solution[level], defect[level]);
+
+#ifdef MG_DEBUG
+ sprintf(name, "MG%d-pre",level);
+ print_vector(level, solution[level], name);
+#endif
+
+ // t = -A*solution[level]
+ matrix->vmult(level, t[level], solution[level]);
+ t[level].scale(-1);
+
+ // make t rhs of lower level
+ // The non-refined parts of the
+ // coarse-level defect already contain
+ // the global defect, the refined parts
+ // its restriction.
+ for (unsigned int l = level;l>minlevel;--l)
+ {
+ t[l-1] = 0.;
+ transfer->restrict_and_add (l, t[l-1], t[l]);
+ defect[l-1] += t[l-1];
+ }
+
+ // add additional DG contribution
+// edge_vmult(level, defect[level-1], defect[level]);
+
+ // do recursion
+ level_mgstep(level-1);
+
+ // 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
+}
+
+
+template <class VECTOR>
void
-Multigrid<dim>::print_vector (const unsigned int /*level*/,
- const Vector<double>& /*v*/,
- const char* /*name*/) const
+Multigrid<VECTOR>::vcycle()
{
+ // 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);
+// abort ();
+}
+
+
+
+// template <int dim>
+// void
+// Multigrid<dim>::print_vector (const unsigned int /*level*/,
+// const Vector<double>& /*v*/,
+// const char* /*name*/) const
+// {
// if (level!=maxlevel)
// return;
// out.add_data_vector(out_vector, "v");
// out.build_patches(5);
// out.write_gnuplot(out_file);
-}
+//}
#endif
template<int dim, class SparsityPattern>
void
MGTools::make_flux_sparsity_pattern_edge (const MGDoFHandler<dim> &dof,
- SparsityPattern &sparsity,
- const unsigned int level)
+ SparsityPattern &sparsity,
+ const unsigned int level)
{
Assert ((level>=1) && (level<dof.get_tria().n_levels()),
ExcIndexRange(level, 1, dof.get_tria().n_levels()));
#include <lac/sparse_matrix.h>
#include <lac/block_sparse_matrix.h>
#include <multigrid/mg_transfer.h>
+#include <multigrid/multigrid.templates.h>
template <typename number>
// Explicit instantiations
+template class Multigrid<Vector<float> >;
+template class Multigrid<Vector<double> >;
+template class Multigrid<BlockVector<float> >;
+template class Multigrid<BlockVector<double> >;
+
template class MGTransferPrebuilt<float>;
template class MGTransferPrebuilt<double>;
template class MGTransferBlock<float>;