* Prolongation and restriction object.
*/
SmartPointer<const MGTransferBase> transfer;
+ /**
+ * Exception.
+ */
+ DeclException2(ExcSwitchedLevels, int, int,
+ << "minlevel and maxlevel switched, should be: "
+ << arg1 << "<=" << arg2);
protected:
* Constructor, subject to change.
*/
MGBase(const MGTransferBase& transfer,
- unsigned int maxlevel, unsigned int minlevel);
+ unsigned int minlevel, unsigned int maxlevel);
/**
* Virtual destructor.
*/
public MGBase
{
public:
- MG(const MGDoFHandler<dim>&, const MGTransferBase& transfer);
+ MG(const MGDoFHandler<dim>&,
+ const MGMatrix<SparseMatrix<float> >&,
+ const MGTransferBase& transfer);
/** Transfer from dVector to
* MGVector.
* MGVector.
*/
template<typename number>
- void copy_to_mg(vector<Vector<float> >& dst,
- const Vector<number>& src) const;
+ void copy_to_mg(const Vector<number>& src);
/**
* Transfer from multi-grid vector to
* #src# are zero.
*/
template<typename number>
- void copy_from_mg(Vector<number>& dst,
- const vector<Vector<float> >& src) const;
+ void copy_from_mg(Vector<number>& dst) const;
/**
- * Access to matrix structure.
+ * Residual on a level.
*/
- SparseMatrixStruct& get_sparsity_pattern(unsigned int level);
-
- /**
- * Read-only access to matrix structure.
- */
- const SparseMatrixStruct& get_sparsity_pattern(unsigned int level) const;
-
- /**
- * Access to level matrices.
- */
- SparseMatrix<float>& get_matrix(unsigned int level);
+ virtual void level_residual(unsigned int,
+ Vector<float> &,
+ const Vector<float> &,
+ const Vector<float> &);
/**
* Read-only access to level matrices.
*/
SmartPointer<const MGDoFHandler<dim> > dofs;
- /**
- * Matrix structures for each level.
- * These are generated by the constructor
- * of #MG# and can be accessed
- * later.
- */
- vector<SparseMatrixStruct> matrix_structures;
-
/**
* Matrices for each level.
* The matrices are prepared by
* the constructor of #MG# and can
* be accessed for assembling.
*/
- vector<SparseMatrix<float> > matrices;
+ SmartPointer<const MGMatrix<SparseMatrix<float> > > matrices;
};
};
-template<int dim>
-inline
-SparseMatrix<float>&
-MG<dim>::get_matrix(unsigned int level)
-{
- Assert((level>=minlevel) && (level<maxlevel), ExcIndexRange(level, minlevel, maxlevel));
-
- return matrices[level];
-}
-
template<int dim>
inline
const SparseMatrix<float>&
{
Assert((level>=minlevel) && (level<maxlevel), ExcIndexRange(level, minlevel, maxlevel));
- return matrices[level];
+ return (*matrices)[level];
}
#include <numerics/multigrid.h>
-/*
template <int dim>
-MultiGrid::
-*/
-
-template <int dim>
-MG<dim>::MG(const MGDoFHandler<dim>& dofs, const MGTransferBase& transfer)
+MG<dim>::MG(const MGDoFHandler<dim>& dofs,
+ const MGMatrix<SparseMatrix<float> >& matrices,
+ const MGTransferBase& transfer)
:
MGBase(transfer, 0, dofs.get_tria().n_levels()-1),
dofs(&dofs),
- matrix_structures(maxlevel-minlevel),
- matrices(maxlevel-minlevel)
+ matrices(&matrices)
{
- for(unsigned l = minlevel; l < maxlevel; ++l)
+ for (unsigned int level = minlevel; level<=maxlevel ; ++level)
{
- dofs.make_sparsity_pattern(l, matrix_structures[l-minlevel]);
- matrices[l-minlevel].reinit(matrix_structures[l-minlevel]);
+ s[level].reinit(matrices[level].m());
+ d[level].reinit(matrices[level].m());
}
}
-template <int dim> template <typename number>
+template <int dim>
+template <typename number>
void
-MG<dim>::copy_to_mg(vector<Vector<float> >& dst,
- const Vector<number>& src) const
+MG<dim>::copy_to_mg(const Vector<number>& src)
{
+ const unsigned int fe_dofs = dofs->get_fe().total_dofs;
+
+ vector<int> index(fe_dofs);
+ vector<int> mgindex(fe_dofs);
+ 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)
+ {
+ dc->get_dof_indices(index);
+ c->get_mg_dof_indices(mgindex);
+ for (unsigned int i=0;i<fe_dofs;++i)
+ d[maxlevel](mgindex[i]) = src(index[i]);
+ }
}
template <int dim> template <typename number>
void
-MG<dim>::copy_from_mg(Vector<number>& dst,
- const vector<Vector<float> >& src) const
+MG<dim>::copy_from_mg(Vector<number>& dst) const
{
-// active iterator
+ const unsigned int fe_dofs = dofs->get_fe().total_dofs;
+
+ vector<int> index(fe_dofs);
+ vector<int> mgindex(fe_dofs);
+
+ 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)
+ {
+ dc->get_dof_indices(index);
+ c->get_mg_dof_indices(mgindex);
+ for (unsigned int i=0;i<fe_dofs;++i)
+ dst(index[i]) = s[maxlevel](mgindex[i]);
+ }
}
+template <int dim>
+void
+MG<dim>::level_residual(unsigned int l,
+ Vector<float> &r,
+ const Vector<float> &u,
+ const Vector<float> &b)
+{
+ (*matrices)[l].residual(r,u,b);
+}
/*
template <int dim>
MG::
template <int dim>
MG::
-template <int dim>
-MG::
*/
public MGBase
{
public:
- MG(const MGDoFHandler<dim>&, const MGTransferBase& transfer);
+ MG(const MGDoFHandler<dim>&,
+ const MGMatrix<SparseMatrix<float> >&,
+ const MGTransferBase& transfer);
/** Transfer from dVector to
* MGVector.
* MGVector.
*/
template<typename number>
- void copy_to_mg(vector<Vector<float> >& dst,
- const Vector<number>& src) const;
+ void copy_to_mg(const Vector<number>& src);
/**
* Transfer from multi-grid vector to
* #src# are zero.
*/
template<typename number>
- void copy_from_mg(Vector<number>& dst,
- const vector<Vector<float> >& src) const;
+ void copy_from_mg(Vector<number>& dst) const;
/**
- * Access to matrix structure.
+ * Residual on a level.
*/
- SparseMatrixStruct& get_sparsity_pattern(unsigned int level);
-
- /**
- * Read-only access to matrix structure.
- */
- const SparseMatrixStruct& get_sparsity_pattern(unsigned int level) const;
-
- /**
- * Access to level matrices.
- */
- SparseMatrix<float>& get_matrix(unsigned int level);
+ virtual void level_residual(unsigned int,
+ Vector<float> &,
+ const Vector<float> &,
+ const Vector<float> &);
/**
* Read-only access to level matrices.
*/
SmartPointer<const MGDoFHandler<dim> > dofs;
- /**
- * Matrix structures for each level.
- * These are generated by the constructor
- * of #MG# and can be accessed
- * later.
- */
- vector<SparseMatrixStruct> matrix_structures;
-
/**
* Matrices for each level.
* The matrices are prepared by
* the constructor of #MG# and can
* be accessed for assembling.
*/
- vector<SparseMatrix<float> > matrices;
+ SmartPointer<const MGMatrix<SparseMatrix<float> > > matrices;
};
};
-template<int dim>
-inline
-SparseMatrix<float>&
-MG<dim>::get_matrix(unsigned int level)
-{
- Assert((level>=minlevel) && (level<maxlevel), ExcIndexRange(level, minlevel, maxlevel));
-
- return matrices[level];
-}
-
template<int dim>
inline
const SparseMatrix<float>&
{
Assert((level>=minlevel) && (level<maxlevel), ExcIndexRange(level, minlevel, maxlevel));
- return matrices[level];
+ return (*matrices)[level];
}
#include <numerics/multigrid.h>
-/*
template <int dim>
-MultiGrid::
-*/
-
-template <int dim>
-MG<dim>::MG(const MGDoFHandler<dim>& dofs, const MGTransferBase& transfer)
+MG<dim>::MG(const MGDoFHandler<dim>& dofs,
+ const MGMatrix<SparseMatrix<float> >& matrices,
+ const MGTransferBase& transfer)
:
MGBase(transfer, 0, dofs.get_tria().n_levels()-1),
dofs(&dofs),
- matrix_structures(maxlevel-minlevel),
- matrices(maxlevel-minlevel)
+ matrices(&matrices)
{
- for(unsigned l = minlevel; l < maxlevel; ++l)
+ for (unsigned int level = minlevel; level<=maxlevel ; ++level)
{
- dofs.make_sparsity_pattern(l, matrix_structures[l-minlevel]);
- matrices[l-minlevel].reinit(matrix_structures[l-minlevel]);
+ s[level].reinit(matrices[level].m());
+ d[level].reinit(matrices[level].m());
}
}
-template <int dim> template <typename number>
+template <int dim>
+template <typename number>
void
-MG<dim>::copy_to_mg(vector<Vector<float> >& dst,
- const Vector<number>& src) const
+MG<dim>::copy_to_mg(const Vector<number>& src)
{
+ const unsigned int fe_dofs = dofs->get_fe().total_dofs;
+
+ vector<int> index(fe_dofs);
+ vector<int> mgindex(fe_dofs);
+ 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)
+ {
+ dc->get_dof_indices(index);
+ c->get_mg_dof_indices(mgindex);
+ for (unsigned int i=0;i<fe_dofs;++i)
+ d[maxlevel](mgindex[i]) = src(index[i]);
+ }
}
template <int dim> template <typename number>
void
-MG<dim>::copy_from_mg(Vector<number>& dst,
- const vector<Vector<float> >& src) const
+MG<dim>::copy_from_mg(Vector<number>& dst) const
{
-// active iterator
+ const unsigned int fe_dofs = dofs->get_fe().total_dofs;
+
+ vector<int> index(fe_dofs);
+ vector<int> mgindex(fe_dofs);
+
+ 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)
+ {
+ dc->get_dof_indices(index);
+ c->get_mg_dof_indices(mgindex);
+ for (unsigned int i=0;i<fe_dofs;++i)
+ dst(index[i]) = s[maxlevel](mgindex[i]);
+ }
}
+template <int dim>
+void
+MG<dim>::level_residual(unsigned int l,
+ Vector<float> &r,
+ const Vector<float> &u,
+ const Vector<float> &b)
+{
+ (*matrices)[l].residual(r,u,b);
+}
/*
template <int dim>
MG::
template <int dim>
MG::
-template <int dim>
-MG::
*/
MGBase::MGBase(const MGTransferBase& transfer,
- unsigned maxlevel, unsigned minlevel)
+ unsigned minlevel, unsigned maxlevel)
:
maxlevel(maxlevel), minlevel(minlevel),
d(minlevel,maxlevel),
s(minlevel,maxlevel),
transfer(&transfer)
{
+ Assert(minlevel <= maxlevel, ExcSwitchedLevels(minlevel, maxlevel));
}
void
//////////////////////////////////////////////////////////////////////
-MGSmootherBase::~MGSmootherBase ()
-{};
-
-//////////////////////////////////////////////////////////////////////
-
void
MGSmoother::set_zero_interior_boundary (const unsigned int level,
Vector<float> &u) const
Vector<float> &dst,
const Vector<float> &src) const
{
- Assert ((to_level >= 1) && (to_level<prolongation_matrices.size()-1),
- ExcIndexRange (to_level, 1, prolongation_matrices.size()-1));
+ Assert ((to_level >= 1) && (to_level<=prolongation_matrices.size()),
+ ExcIndexRange (to_level, 1, prolongation_matrices.size()+1));
prolongation_matrices[to_level-1].vmult (dst, src);
};
Vector<float> &dst,
const Vector<float> &src) const
{
- Assert ((from_level >= 1) && (from_level<prolongation_matrices.size()-1),
- ExcIndexRange (from_level, 1, prolongation_matrices.size()-1));
+ Assert ((from_level >= 1) && (from_level<=prolongation_matrices.size()),
+ ExcIndexRange (from_level, 1, prolongation_matrices.size()+1));
prolongation_matrices[from_level-1].Tvmult (dst, src);
};
//////////////////////////////////////////////////////////////////////
-MGSmootherBase::~MGSmootherBase ()
-{};
-
-//////////////////////////////////////////////////////////////////////
-
void
MGSmoother::set_zero_interior_boundary (const unsigned int level,
Vector<float> &u) const
Vector<float> &dst,
const Vector<float> &src) const
{
- Assert ((to_level >= 1) && (to_level<prolongation_matrices.size()-1),
- ExcIndexRange (to_level, 1, prolongation_matrices.size()-1));
+ Assert ((to_level >= 1) && (to_level<=prolongation_matrices.size()),
+ ExcIndexRange (to_level, 1, prolongation_matrices.size()+1));
prolongation_matrices[to_level-1].vmult (dst, src);
};
Vector<float> &dst,
const Vector<float> &src) const
{
- Assert ((from_level >= 1) && (from_level<prolongation_matrices.size()-1),
- ExcIndexRange (from_level, 1, prolongation_matrices.size()-1));
+ Assert ((from_level >= 1) && (from_level<=prolongation_matrices.size()),
+ ExcIndexRange (from_level, 1, prolongation_matrices.size()+1));
prolongation_matrices[from_level-1].Tvmult (dst, src);
};
* Prolongation and restriction object.
*/
SmartPointer<const MGTransferBase> transfer;
+ /**
+ * Exception.
+ */
+ DeclException2(ExcSwitchedLevels, int, int,
+ << "minlevel and maxlevel switched, should be: "
+ << arg1 << "<=" << arg2);
protected:
* Constructor, subject to change.
*/
MGBase(const MGTransferBase& transfer,
- unsigned int maxlevel, unsigned int minlevel);
+ unsigned int minlevel, unsigned int maxlevel);
/**
* Virtual destructor.
*/
template <typename number>
-const SparseMatrixStruct & SparseMatrix<number>::get_sparsity_pattern () const {
+const SparseMatrixStruct &
+SparseMatrix<number>::get_sparsity_pattern () const
+{
Assert (cols != 0, ExcMatrixNotInitialized());
return *cols;
};
MGBase::MGBase(const MGTransferBase& transfer,
- unsigned maxlevel, unsigned minlevel)
+ unsigned minlevel, unsigned maxlevel)
:
maxlevel(maxlevel), minlevel(minlevel),
d(minlevel,maxlevel),
s(minlevel,maxlevel),
transfer(&transfer)
{
+ Assert(minlevel <= maxlevel, ExcSwitchedLevels(minlevel, maxlevel));
}
void
dof.distribute_dofs(fe);
QTrapez<2> q;
- FEValues<2> fevalues(fe,q,update_gradients,tria.get_boundary());
+ FEValues<2> fevalues(fe,q,update_gradients);
fevalues.reinit (dof.begin_active());
--- /dev/null
+// $Id$
+
+#include <lac/mgbase.h>
+
+/**
+ * Generator for system matrix and right hand side.
+ */
+class Helmholtz
+{
+ public:
+ template<int dim, class MATRIX, class VECTOR>
+ void build_all(MATRIX& A,
+ VECTOR& f,
+ const DoFHandler<dim>& dof,
+ const Quadrature<dim>& quadrature,
+ const Function<dim>& rhs);
+
+ template<int dim, class MATRIX>
+ void build_mgmatrix(MGMatrix<MATRIX>& A,
+ const MGDoFHandler<dim>& dof,
+ const Quadrature<dim>& quadrature);
+};
+
+#include "helmholtz1.th"
+#include "helmholtz1mg.th"
+
--- /dev/null
+// $Id$
+
+template<int dim, class MATRIX, class VECTOR>
+void
+Helmholtz::build_all(MATRIX& A,
+ VECTOR& f,
+ const DoFHandler<dim>& dof,
+ const Quadrature<dim>& quadrature,
+ const Function<dim>& rhs)
+{
+ const unsigned int fe_dofs = dof.get_fe().total_dofs;
+ FEValues<dim> fe(dof.get_fe(), quadrature,
+ UpdateFlags(update_gradients | update_JxW_values |
+ update_q_points));
+
+ Vector<double> elvec(fe_dofs);
+ FullMatrix<double> elmat(fe_dofs);
+ vector<int> indices(fe_dofs);
+
+////////////////////////////////////////////////////
+// Loop for computing the element matrices
+////////////////////////////////////////////////////
+
+ for (DoFHandler<dim>::active_cell_iterator c = dof.begin_active()
+ ; c != dof.end() ; ++c)
+ {
+ elmat.clear();
+ elvec.clear();
+ fe.reinit(c);
+
+ // Quadrature loop
+
+ for (unsigned k=0;k<quadrature.n_quadrature_points;++k)
+ {
+ Point<dim> loc = fe.quadrature_point(k);
+
+ // Test function loop
+ for (unsigned i=0;i<fe_dofs;++i)
+ {
+ const Point<dim> dv = fe.shape_grad(i,k);
+ double v = fe.shape_value(i,k);
+
+ elvec(i) += fe.JxW(k) *
+ rhs(loc) * v;
+
+ //Trial function loop
+ for (unsigned j=0;j<fe_dofs;++j)
+ {
+ const Point<dim> du = fe.shape_grad(j,k);
+ double u = fe.shape_value(j,k);
+
+ elmat(i,j) += fe.JxW(k) *
+ (1.*u*v + du*dv)
+ ;
+ }
+ }
+ }
+
+ c->get_dof_indices(indices);
+ for (unsigned i=0;i<fe_dofs;++i)
+ {
+ f(indices[i]) += elvec(i);
+
+ for (unsigned j=0;j<fe_dofs;++j)
+ {
+ A.add(indices[i], indices[j], elmat(i,j));
+ }
+ }
+ }
+}
--- /dev/null
+// $Id$
+
+template<int dim, class MATRIX>
+void
+Helmholtz::build_mgmatrix(MGMatrix<MATRIX>& A,
+ const MGDoFHandler<dim>& dof,
+ const Quadrature<dim>& quadrature)
+{
+ const unsigned int fe_dofs = dof.get_fe().total_dofs;
+ FEValues<dim> fe(dof.get_fe(), quadrature,
+ UpdateFlags(update_gradients | update_JxW_values |
+ update_q_points));
+
+ Vector<double> elvec(fe_dofs);
+ FullMatrix<double> elmat(fe_dofs);
+ vector<int> indices(fe_dofs);
+
+////////////////////////////////////////////////////
+// Loop for computing the element matrices
+////////////////////////////////////////////////////
+
+ DoFHandler<dim>::cell_iterator dc = dof.DoFHandler<dim>::begin();
+ for (MGDoFHandler<dim>::cell_iterator c = dof.begin()
+ ; c != dof.end() ; ++c, ++dc)
+ {
+ elmat.clear();
+ elvec.clear();
+ fe.reinit(dc);
+
+ // Quadrature loop
+
+ for (unsigned k=0;k<quadrature.n_quadrature_points;++k)
+ {
+// Point<dim> loc = fe.quadrature_point(k);
+
+ // Test function loop
+ for (unsigned i=0;i<fe_dofs;++i)
+ {
+ const Point<dim> dv = fe.shape_grad(i,k);
+ double v = fe.shape_value(i,k);
+
+ //Trial function loop
+ for (unsigned j=0;j<fe_dofs;++j)
+ {
+ const Point<dim> du = fe.shape_grad(j,k);
+ double u = fe.shape_value(j,k);
+
+ elmat(i,j) += fe.JxW(k) *
+ (1.*u*v + du*dv)
+ ;
+ }
+ }
+ }
+
+ const unsigned int level = c->level();
+
+ c->get_mg_dof_indices(indices);
+ for (unsigned i=0;i<fe_dofs;++i)
+ {
+ for (unsigned j=0;j<fe_dofs;++j)
+ {
+ A[level].add(indices[i], indices[j], elmat(i,j));
+ }
+ }
+ }
+}
--- /dev/null
+// $Id$
+
+// deal_II_libraries.g=-ldeal_II_2d.g
+// deal_II_libraries=-ldeal_II_2d
+
+#include <base/function.h>
+#include <base/quadrature_lib.h>
+#include <lac/vector.h>
+#include <lac/sparsematrix.h>
+#include <lac/solver_cg.h>
+#include <lac/solver_richardson.h>
+#include <lac/vector_memory.h>
+#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <grid/tria_iterator.h>
+#include <grid/mg_dof.h>
+#include <grid/dof_accessor.h>
+#include <grid/mg_dof_accessor.h>
+#include <grid/grid_generator.h>
+#include <fe/fe_lib.lagrange.h>
+#include <fe/fe_values.h>
+#include <numerics/multigrid.h>
+#include <numerics/multigrid.templates.h>
+#include <numerics/mg_smoother.h>
+
+#include "helmholtz.h"
+
+template<int dim>
+class RHSFunction
+ :
+ public Function<dim>
+{
+ public:
+ virtual double operator() (const Point<dim>&) const;
+// virtual Tensor<1,dim> gradient (const Point<dim> &p) const;
+};
+
+class MGSmootherLAC
+ :
+ public MGSmootherBase
+{
+ private:
+ SmartPointer<MGMatrix<SparseMatrix<float> > >matrices;
+ public:
+ MGSmootherLAC(MGMatrix<SparseMatrix<float> >&);
+
+ virtual void smooth (const unsigned int level,
+ Vector<float> &u,
+ const Vector<float> &rhs) const;
+
+};
+
+main(int argc, const char** argv)
+{
+ unsigned int refine;
+ if (argc<=1)
+ {
+ cerr << "Usage:" << argv[0] << " <refinement>" << endl;
+ refine = 0;
+ }
+ else
+ refine = atoi(argv[1]);
+
+
+ Helmholtz equation;
+ RHSFunction<2> rhs;
+ QGauss3<2> quadrature;
+
+ FELinear<2> fe;
+ Triangulation<2> tr;
+ MGDoFHandler<2> mgdof(&tr);
+ DoFHandler<2>& dof(mgdof);
+
+ GridGenerator::hyper_cube(tr,-1.,1.);
+ tr.refine_global(refine);
+ dof.distribute_dofs(fe);
+
+ const unsigned int size = dof.n_dofs();
+
+ SparseMatrixStruct structure(size, dof.max_couplings_between_dofs());
+ dof.make_sparsity_pattern(structure);
+ structure.compress();
+
+ SparseMatrix<double> A(structure);
+ Vector<double> f(size);
+
+ equation.build_all(A,f,dof, quadrature, rhs);
+
+ Vector<double> u;
+ u.reinit(f);
+ PrimitiveVectorMemory<Vector<double> > mem;
+ SolverControl control(100, 1.e-15);
+ SolverCG<SparseMatrix<double>, Vector<double> >solver(control, mem);
+
+ PreconditionIdentity<Vector<double> > precondition;
+ solver.solve(A,u,f,precondition);
+
+ u = 0.;
+ vector<SparseMatrixStruct> mgstruct(tr.n_levels());
+ MGMatrix<SparseMatrix<float> > mgA(0,tr.n_levels()-1);
+ for (unsigned int i=0;i<tr.n_levels();++i)
+ {
+ mgstruct[i].reinit(mgdof.n_dofs(i), mgdof.n_dofs(i),
+ mgdof.max_couplings_between_dofs());
+ mgdof.make_sparsity_pattern(i, mgstruct[i]);
+ mgstruct[i].compress();
+
+ mgA[i].reinit(mgstruct[i]);
+ }
+ equation.build_mgmatrix(mgA, mgdof, quadrature);
+
+ SolverControl cgcontrol(10,0., false, false);
+
+ PrimitiveVectorMemory<Vector<float> > cgmem;
+ SolverCG<SparseMatrix<float>, Vector<float> > cgsolver(cgcontrol, cgmem);
+ PreconditionIdentity<Vector<float> > cgprec;
+ MGCoarseGridLACIteration<SolverCG<SparseMatrix<float>, Vector<float> >,
+ SparseMatrix<float>, PreconditionIdentity<Vector<float> > >
+ coarse(cgsolver, mgA[0], cgprec);
+
+ MGSmootherLAC smoother(mgA);
+ MGTransferPrebuilt transfer;
+ transfer.build_matrices(mgdof);
+
+ deallog << "Levels: " << tr.n_levels() << endl;
+
+ MG<2> multigrid(mgdof, mgA, transfer);
+ PreconditionMG<MG<2>, Vector<double> >
+ mgprecondition(multigrid, smoother, smoother, coarse);
+
+ solver.solve(A, u, f, mgprecondition);
+
+ cerr << "OK" << endl;
+}
+
+template<int dim>
+double
+RHSFunction<dim>::operator() (const Point<dim>&) const
+{
+ return 1.;
+}
+
+MGSmootherLAC::MGSmootherLAC(MGMatrix<SparseMatrix<float> >& matrix)
+ :
+ matrices(&matrix)
+{}
+
+void
+MGSmootherLAC::smooth (const unsigned int level,
+ Vector<float> &u,
+ const Vector<float> &rhs) const
+{
+ SolverControl control(1,1.e-300,false,false);
+ PrimitiveVectorMemory<Vector<float> > mem;
+ SolverRichardson<SparseMatrix<float> , Vector<float> > rich(control, mem);
+ PreconditionRelaxation<SparseMatrix<float> , Vector<float> >
+ prec((*matrices)[level], &SparseMatrix<float> ::template precondition_SSOR<float>, 1.);
+
+ rich.solve((*matrices)[level], u, rhs, prec);
+}
+++ /dev/null
-// $Id$
-
-// deal_II_libraries.g=-ldeal_II_2d.g
-// deal_II_libraries=-ldeal_II_2d
-
-#include <numerics/multigrid.h>
-#include <numerics/mg_smoother.h>
-#include <grid/tria.h>
-#include <grid/mg_dof.h>
-#include <grid/grid_generator.h>
-#include <fe/fe_lib.lagrange.h>
-
-main()
-{
- Triangulation<2> tr;
- MGDoFHandler<2> dof(&tr);
- FELinear<2> fe;
-
- GridGenerator::hyper_cube(tr,-1.,1.);
- tr.refine_global(3);
- dof.distribute_dofs(fe);
-}
-
StraightBoundary<2> b;
QTrapez<2> q;
- FEValues<2> fevalues(fe,q,update_second_derivatives,b);
+ FEValues<2> fevalues(fe,q,update_second_derivatives);
Vector<double> val(4);
FDMG::FDMG(unsigned int maxlevel, MGMatrix<SparseMatrix<float> >& matrices,
FDMGTransfer& transfer)
:
- MGBase(transfer, maxlevel, 0),
+ MGBase(transfer, 0, maxlevel),
matrices(&matrices)
{
for (unsigned int level = minlevel; level<=maxlevel ; ++level)