template<typename number>
void distribute (Vector<number> &vec) const;
+ /**
+ * Delete hanging nodes in a vector.
+ * Sets all hanging node values to zero.
+ */
+ template<typename number>
+ void set_zero (Vector<number> &vec) const;
+
/**
* Print the constraint lines. Mainly for
+template<typename number>
+void
+ConstraintMatrix::set_zero (Vector<number> &vec) const
+{
+ Assert (sorted == true, ExcMatrixNotClosed());
+
+ if (lines.size() == 0)
+ // nothing to do
+ return;
+
+ vector<ConstraintLine>::const_iterator next_constraint = lines.begin();
+ for (unsigned int row=0; row<vec.size(); ++row)
+ if (row == next_constraint->line)
+ {
+ // set entry to zero
+ vec(row) = 0.;
+
+ ++next_constraint;
+ if (next_constraint == lines.end())
+ // nothing more to do
+ break;
+ };
+};
+
+
+
template<typename number>
void
ConstraintMatrix::distribute (const Vector<number> &condensed,
* about the matrix is removed to
* that class.
*/
- virtual void operator() (unsigned int level, Vector<float>& dst,
- const Vector<float>& src) const = 0;
+ virtual void operator() (unsigned int level, Vector<double>& dst,
+ const Vector<double>& src) const = 0;
};
* things.
*/
virtual void smooth (const unsigned int level,
- Vector<float> &u,
- const Vector<float> &rhs) const = 0;
+ Vector<double> &u,
+ const Vector<double> &rhs) const = 0;
};
* This function does nothing.
*/
virtual void smooth (const unsigned int level,
- Vector<float> &u,
- const Vector<float> &rhs) const;
+ Vector<double> &u,
+ const Vector<double> &rhs) const;
};
/**
* level.
*/
virtual void prolongate (const unsigned int to_level,
- Vector<float> &dst,
- const Vector<float> &src) const = 0;
+ Vector<double> &dst,
+ const Vector<double> &src) const = 0;
/**
* Restrict a vector from level
* level.
*/
virtual void restrict (const unsigned int from_level,
- Vector<float> &dst,
- const Vector<float> &src) const = 0;
+ Vector<double> &dst,
+ const Vector<double> &src) const = 0;
};
/**
* Safe access operator.
*/
const VECTOR& operator[](unsigned int) const;
+ /**
+ * Delete all entries.
+ */
+ void clear();
private:
/**
* Level of first component.
* matrix, vectors and
* preconditioner.
*/
- virtual void operator() (unsigned int level, Vector<float>& dst,
- const Vector<float>& src) const;
+ virtual void operator() (unsigned int level, Vector<double>& dst,
+ const Vector<double>& src) const;
private:
/**
* Reference to the solver.
/**
* Auxiliary vector, defect.
*/
- MGVector<Vector<float> > d;
+ MGVector<Vector<double> > d;
/**
* Auxiliary vector, solution.
*/
- MGVector<Vector<float> > s;
+ MGVector<Vector<double> > s;
+ /**
+ * Prolongation and restriction object.
+ */
+ SmartPointer<const MGTransferBase> transfer;
private:
/**
* Auxiliary vector.
*/
- Vector<float> t;
+ Vector<double> t;
- /**
- * Prolongation and restriction object.
- */
- SmartPointer<const MGTransferBase> transfer;
/**
* Exception.
*/
const MGCoarseGridSolver& cgs);
/**
- * Apply residual operator on all
+ * Apply negative #vmult# on all
* cells of a level.
* This is implemented in a
* derived class.
*/
- virtual void level_residual(unsigned int level,
- Vector<float>& dst,
- const Vector<float>& src,
- const Vector<float>& rhs) = 0;
+ virtual void level_vmult(unsigned int level,
+ Vector<double>& dst,
+ const Vector<double>& src,
+ const Vector<double>& rhs) = 0;
public:
template<class SOLVER, class MATRIX, class PRECOND>
void
MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>::operator()
-(unsigned int, Vector<float>& dst, const Vector<float>& src) const
+(unsigned int, Vector<double>& dst, const Vector<double>& src) const
{
solver.solve(matrix, dst, src, precondition);
}
return vector<VECTOR>::operator[](i-minlevel);
}
+template<class VECTOR>
+void
+MGVector<VECTOR>::clear()
+{
+ vector<VECTOR>::iterator v;
+ for (v = begin(); v != end(); ++v)
+ v->clear();
+
+}
/* ------------------------------------------------------------------- */
*
* @author Wolfgang Bangerth, Guido Kanschat, 1999
*/
-class MGSmoother
- : public MGSmootherBase
+class MGSmoother : public MGSmootherBase
{
private:
/**
* function does nothing in this case.
*/
void set_zero_interior_boundary (const unsigned int level,
- Vector<float> &u) const;
+ Vector<double> &u) const;
/**
* Modify the number of smoothing steps.
*/
};
/**
- * Implementation of an SOR smoother.
- *
- * This class is an implementation of a smootin algorithm for multigrid. It
- * knows the actual instance of the multigrid matrix and uses the SOR method of the level
- * matrices.
+ * Implementation of a smoother using matrix builtin relaxation methods.
+ *
* @author Wolfgang Bangerth, Guido Kanschat, 1999
*/
-class MGSmootherSOR
- :
- public MGSmoother
+template<typename number>
+class MGSmootherRelaxation : public MGSmoother
{
public:
+ /**
+ * Type of the smoothing
+ * function of the matrix.
+ */
+ typedef void
+ (SparseMatrix<number>::* function_ptr)(Vector<double>&, const Vector<double>&,
+ typename SparseMatrix<number>::value_type) const;
+
/**
* Constructor.
* The constructor uses an #MGDoFHandler# to initialize
* data structures for interior level boundary handling
- * in #MGSmoother#. Furthermore, it takes a pointer to the
- * matrices and the number of SOR steps required in each
- * smoothing operation.
+ * in #MGSmoother#.
+ *
+ * Furthermore, it takes a pointer to the
+ * level matrices and their smoothing function.
+ * This function must perform one relaxation step
+ * like #SparseMatrix<number>::SOR_step# does. Do not
+ * use the preconditioning methods because they apply
+ * a preconditioning matrix to the residual.
+ *
+ * The final two parameters are the number of relaxation
+ * steps and the relaxation parameter.
*/
template<int dim>
- MGSmootherSOR(const MGDoFHandler<dim> &mg_dof,
- const MGMatrix<SparseMatrix<float> >& matrix,
- unsigned int steps);
+ MGSmootherRelaxation(const MGDoFHandler<dim> &mg_dof,
+ const MGMatrix<SparseMatrix<number> >& matrix,
+ function_ptr function,
+ unsigned int steps,
+ double omega = 1.);
/**
* Implementation of the interface in #MGSmootherBase#.
* and find a way to keep the boundary values.
*/
virtual void smooth (const unsigned int level,
- Vector<float> &u,
- const Vector<float> &rhs) const;
+ Vector<double> &u,
+ const Vector<double> &rhs) const;
private:
/**
* Pointer to the matrices.
*/
- SmartPointer< const MGMatrix< SparseMatrix<float> > > matrix;
+ SmartPointer< const MGMatrix< SparseMatrix<number> > > matrix;
+ /**
+ * Pointer to the relaxation function.
+ */
+ function_ptr relaxation;
+ /**
+ * Relaxation parameter.
+ */
+ double omega;
+ /**
+ * Auxiliary vector.
+ */
+ mutable Vector<double> h1, h2;
};
inline
--- /dev/null
+// $Id$
+
+template<typename number>
+template<int dim>
+MGSmootherRelaxation<number>::
+MGSmootherRelaxation(const MGDoFHandler<dim> &mg_dof,
+ const MGMatrix<SparseMatrix<number> >& matrix,
+ function_ptr relaxation,
+ unsigned int steps,
+ double omega)
+ :
+ MGSmoother(mg_dof, steps),
+ matrix(&matrix),
+ relaxation(relaxation),
+ omega(omega)
+{}
+
+template<typename number>
+void
+MGSmootherRelaxation<number>::smooth (const unsigned int level,
+ Vector<double> &u,
+ const Vector<double> &rhs) const
+{
+ h1.reinit(u);
+ h2.reinit(u);
+ for(unsigned i=0;i<get_steps();++i)
+ {
+ (*matrix)[level].residual(h1, u, rhs);
+ ((*matrix)[level].*relaxation)(h2,h1,omega);
+ set_zero_interior_boundary(level,h2);
+ u.add(h2);
+ }
+}
+
{
public:
MG(const MGDoFHandler<dim>&,
- const MGMatrix<SparseMatrix<float> >&,
- const MGTransferBase& transfer);
+ ConstraintMatrix& hanging_nodes,
+ const MGMatrix<SparseMatrix<double> >&,
+ const MGTransferBase& transfer,
+ unsigned int minlevel = 0, unsigned int maxlevel = 10000);
/** Transfer from dVector to
* MGVector.
void copy_from_mg(Vector<number>& dst) const;
/**
- * Residual on a level.
+ * Negative #vmult# on a level.
+ * @see MGBase.
*/
- virtual void level_residual(unsigned int,
- Vector<float> &,
- const Vector<float> &,
- const Vector<float> &);
+ virtual void level_vmult(unsigned int,
+ Vector<double> &,
+ const Vector<double> &,
+ const Vector<double> &);
/**
* Read-only access to level matrices.
*/
- const SparseMatrix<float>& get_matrix(unsigned int level) const;
+ const SparseMatrix<double>& get_matrix(unsigned int level) const;
private:
/**
* the constructor of #MG# and can
* be accessed for assembling.
*/
- SmartPointer<const MGMatrix<SparseMatrix<float> > > matrices;
+ SmartPointer<const MGMatrix<SparseMatrix<double> > > matrices;
+ ConstraintMatrix& hanging_nodes;
};
*/
class MGTransferPrebuilt : public MGTransferBase
{
+ const MGSmoother& smoother;
public:
-
+ /**
+ * Preliminary constructor
+ */
+ MGTransferPrebuilt(const MGSmoother& smoother) :
+ smoother(smoother)
+ {}
+
+
/**
* Destructor.
*/
* level.
*/
virtual void prolongate (const unsigned int to_level,
- Vector<float> &dst,
- const Vector<float> &src) const;
+ Vector<double> &dst,
+ const Vector<double> &src) const;
/**
* Restrict a vector from level
* level.
*/
virtual void restrict (const unsigned int from_level,
- Vector<float> &dst,
- const Vector<float> &src) const;
+ Vector<double> &dst,
+ const Vector<double> &src) const;
private:
template<int dim>
inline
-const SparseMatrix<float>&
+const SparseMatrix<double>&
MG<dim>::get_matrix(unsigned int level) const
{
Assert((level>=minlevel) && (level<maxlevel), ExcIndexRange(level, minlevel, maxlevel));
// $Id$
#include <numerics/multigrid.h>
+#include <fstream>
+
+inline int minimum(int a, int b)
+{
+ return ((a<b) ? a : b);
+}
template <int dim>
MG<dim>::MG(const MGDoFHandler<dim>& dofs,
- const MGMatrix<SparseMatrix<float> >& matrices,
- const MGTransferBase& transfer)
+ ConstraintMatrix& hanging_nodes,
+ const MGMatrix<SparseMatrix<double> >& matrices,
+ const MGTransferBase& transfer,
+ unsigned int minl, unsigned int maxl)
:
- MGBase(transfer, 0, dofs.get_tria().n_levels()-1),
+ MGBase(transfer, minl,
+ minimum(dofs.get_tria().n_levels()-1, maxl)),
dofs(&dofs),
- matrices(&matrices)
+ matrices(&matrices),
+ hanging_nodes(hanging_nodes)
{
for (unsigned int level = minlevel; level<=maxlevel ; ++level)
{
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;
+
+ d.clear();
+ hanging_nodes.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","m");
+ out.write_gnuplot(of,1);
+ }
- 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)
+ for (int level = maxlevel ; level >= (int) minlevel ; --level)
{
- 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]);
- }
+ DoFHandler<dim>::active_cell_iterator dc = dofs->DoFHandler<dim>::begin_active(level);
+ MGDoFHandler<dim>::active_cell_iterator c = dofs->begin_active(level);
+ for (; c != dofs->end(level) ; ++c, ++dc)
+ {
+ dc->get_dof_indices(index);
+ c->get_mg_dof_indices(mgindex);
+ for (unsigned int i=0;i<fe_dofs;++i)
+ d[c->level()](mgindex[i]) = src(index[i]);
+
+ // Delete values on refinement edge
+ for (unsigned int face_n = 0; face_n< GeometryInfo<dim>::faces_per_cell; ++face_n)
+ {
+ MGDoFHandler<dim>::face_iterator face = c->face(face_n);
+ if (face->has_children())
+ {
+ face->get_mg_dof_indices(mgfaceindex);
+ for (unsigned int i=0;i<face_dofs;++i)
+ d[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, d[level], d[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);
+ }
}
template <int dim> template <typename number>
vector<int> index(fe_dofs);
vector<int> mgindex(fe_dofs);
+ for (unsigned int i=minlevel; i<= maxlevel; ++i)
+ {
+ char name[10];
+ sprintf(name,"CF%d",i);
+ ofstream of(name);
+ 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)
+ ; 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]);
- }
+ dst(index[i]) = s[c->level()](mgindex[i]);
+ }
+ hanging_nodes.set_zero(dst);
+ {
+ ofstream of("CF");
+ DataOut<2> out;
+ out.attach_dof_handler(*dofs);
+ out.add_data_vector(dst,"solution","m");
+ out.write_gnuplot(of,1);
+ }
}
template <int dim>
void
-MG<dim>::level_residual(unsigned int l,
- Vector<float> &r,
- const Vector<float> &u,
- const Vector<float> &b)
+MG<dim>::level_vmult(unsigned int l,
+ Vector<double> &r,
+ const Vector<double> &u,
+ const Vector<double> &)
{
- (*matrices)[l].residual(r,u,b);
+ (*matrices)[l].vmult(r,u);
+ r.scale(-1.);
}
/*
* used for the height information, within
* the list of DoF data vectors.
*/
- unsigned height_vector;
+ unsigned int height_vector;
/**
* Number of the vector which is to be
* used for the cell shading values, within
* the list of cell data vectors.
*/
- unsigned cell_vector;
+ unsigned int cell_vector;
/**
* Azimuth of the spectators position.
*
* @author Wolfgang Bangerth, Guido Kanschat, 1999
*/
-class MGSmoother
- : public MGSmootherBase
+class MGSmoother : public MGSmootherBase
{
private:
/**
* function does nothing in this case.
*/
void set_zero_interior_boundary (const unsigned int level,
- Vector<float> &u) const;
+ Vector<double> &u) const;
/**
* Modify the number of smoothing steps.
*/
};
/**
- * Implementation of an SOR smoother.
- *
- * This class is an implementation of a smootin algorithm for multigrid. It
- * knows the actual instance of the multigrid matrix and uses the SOR method of the level
- * matrices.
+ * Implementation of a smoother using matrix builtin relaxation methods.
+ *
* @author Wolfgang Bangerth, Guido Kanschat, 1999
*/
-class MGSmootherSOR
- :
- public MGSmoother
+template<typename number>
+class MGSmootherRelaxation : public MGSmoother
{
public:
+ /**
+ * Type of the smoothing
+ * function of the matrix.
+ */
+ typedef void
+ (SparseMatrix<number>::* function_ptr)(Vector<double>&, const Vector<double>&,
+ typename SparseMatrix<number>::value_type) const;
+
/**
* Constructor.
* The constructor uses an #MGDoFHandler# to initialize
* data structures for interior level boundary handling
- * in #MGSmoother#. Furthermore, it takes a pointer to the
- * matrices and the number of SOR steps required in each
- * smoothing operation.
+ * in #MGSmoother#.
+ *
+ * Furthermore, it takes a pointer to the
+ * level matrices and their smoothing function.
+ * This function must perform one relaxation step
+ * like #SparseMatrix<number>::SOR_step# does. Do not
+ * use the preconditioning methods because they apply
+ * a preconditioning matrix to the residual.
+ *
+ * The final two parameters are the number of relaxation
+ * steps and the relaxation parameter.
*/
template<int dim>
- MGSmootherSOR(const MGDoFHandler<dim> &mg_dof,
- const MGMatrix<SparseMatrix<float> >& matrix,
- unsigned int steps);
+ MGSmootherRelaxation(const MGDoFHandler<dim> &mg_dof,
+ const MGMatrix<SparseMatrix<number> >& matrix,
+ function_ptr function,
+ unsigned int steps,
+ double omega = 1.);
/**
* Implementation of the interface in #MGSmootherBase#.
* and find a way to keep the boundary values.
*/
virtual void smooth (const unsigned int level,
- Vector<float> &u,
- const Vector<float> &rhs) const;
+ Vector<double> &u,
+ const Vector<double> &rhs) const;
private:
/**
* Pointer to the matrices.
*/
- SmartPointer< const MGMatrix< SparseMatrix<float> > > matrix;
+ SmartPointer< const MGMatrix< SparseMatrix<number> > > matrix;
+ /**
+ * Pointer to the relaxation function.
+ */
+ function_ptr relaxation;
+ /**
+ * Relaxation parameter.
+ */
+ double omega;
+ /**
+ * Auxiliary vector.
+ */
+ mutable Vector<double> h1, h2;
};
inline
--- /dev/null
+// $Id$
+
+template<typename number>
+template<int dim>
+MGSmootherRelaxation<number>::
+MGSmootherRelaxation(const MGDoFHandler<dim> &mg_dof,
+ const MGMatrix<SparseMatrix<number> >& matrix,
+ function_ptr relaxation,
+ unsigned int steps,
+ double omega)
+ :
+ MGSmoother(mg_dof, steps),
+ matrix(&matrix),
+ relaxation(relaxation),
+ omega(omega)
+{}
+
+template<typename number>
+void
+MGSmootherRelaxation<number>::smooth (const unsigned int level,
+ Vector<double> &u,
+ const Vector<double> &rhs) const
+{
+ h1.reinit(u);
+ h2.reinit(u);
+ for(unsigned i=0;i<get_steps();++i)
+ {
+ (*matrix)[level].residual(h1, u, rhs);
+ ((*matrix)[level].*relaxation)(h2,h1,omega);
+ set_zero_interior_boundary(level,h2);
+ u.add(h2);
+ }
+}
+
{
public:
MG(const MGDoFHandler<dim>&,
- const MGMatrix<SparseMatrix<float> >&,
- const MGTransferBase& transfer);
+ ConstraintMatrix& hanging_nodes,
+ const MGMatrix<SparseMatrix<double> >&,
+ const MGTransferBase& transfer,
+ unsigned int minlevel = 0, unsigned int maxlevel = 10000);
/** Transfer from dVector to
* MGVector.
void copy_from_mg(Vector<number>& dst) const;
/**
- * Residual on a level.
+ * Negative #vmult# on a level.
+ * @see MGBase.
*/
- virtual void level_residual(unsigned int,
- Vector<float> &,
- const Vector<float> &,
- const Vector<float> &);
+ virtual void level_vmult(unsigned int,
+ Vector<double> &,
+ const Vector<double> &,
+ const Vector<double> &);
/**
* Read-only access to level matrices.
*/
- const SparseMatrix<float>& get_matrix(unsigned int level) const;
+ const SparseMatrix<double>& get_matrix(unsigned int level) const;
private:
/**
* the constructor of #MG# and can
* be accessed for assembling.
*/
- SmartPointer<const MGMatrix<SparseMatrix<float> > > matrices;
+ SmartPointer<const MGMatrix<SparseMatrix<double> > > matrices;
+ ConstraintMatrix& hanging_nodes;
};
*/
class MGTransferPrebuilt : public MGTransferBase
{
+ const MGSmoother& smoother;
public:
-
+ /**
+ * Preliminary constructor
+ */
+ MGTransferPrebuilt(const MGSmoother& smoother) :
+ smoother(smoother)
+ {}
+
+
/**
* Destructor.
*/
* level.
*/
virtual void prolongate (const unsigned int to_level,
- Vector<float> &dst,
- const Vector<float> &src) const;
+ Vector<double> &dst,
+ const Vector<double> &src) const;
/**
* Restrict a vector from level
* level.
*/
virtual void restrict (const unsigned int from_level,
- Vector<float> &dst,
- const Vector<float> &src) const;
+ Vector<double> &dst,
+ const Vector<double> &src) const;
private:
template<int dim>
inline
-const SparseMatrix<float>&
+const SparseMatrix<double>&
MG<dim>::get_matrix(unsigned int level) const
{
Assert((level>=minlevel) && (level<maxlevel), ExcIndexRange(level, minlevel, maxlevel));
// $Id$
#include <numerics/multigrid.h>
+#include <fstream>
+
+inline int minimum(int a, int b)
+{
+ return ((a<b) ? a : b);
+}
template <int dim>
MG<dim>::MG(const MGDoFHandler<dim>& dofs,
- const MGMatrix<SparseMatrix<float> >& matrices,
- const MGTransferBase& transfer)
+ ConstraintMatrix& hanging_nodes,
+ const MGMatrix<SparseMatrix<double> >& matrices,
+ const MGTransferBase& transfer,
+ unsigned int minl, unsigned int maxl)
:
- MGBase(transfer, 0, dofs.get_tria().n_levels()-1),
+ MGBase(transfer, minl,
+ minimum(dofs.get_tria().n_levels()-1, maxl)),
dofs(&dofs),
- matrices(&matrices)
+ matrices(&matrices),
+ hanging_nodes(hanging_nodes)
{
for (unsigned int level = minlevel; level<=maxlevel ; ++level)
{
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;
+
+ d.clear();
+ hanging_nodes.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","m");
+ out.write_gnuplot(of,1);
+ }
- 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)
+ for (int level = maxlevel ; level >= (int) minlevel ; --level)
{
- 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]);
- }
+ DoFHandler<dim>::active_cell_iterator dc = dofs->DoFHandler<dim>::begin_active(level);
+ MGDoFHandler<dim>::active_cell_iterator c = dofs->begin_active(level);
+ for (; c != dofs->end(level) ; ++c, ++dc)
+ {
+ dc->get_dof_indices(index);
+ c->get_mg_dof_indices(mgindex);
+ for (unsigned int i=0;i<fe_dofs;++i)
+ d[c->level()](mgindex[i]) = src(index[i]);
+
+ // Delete values on refinement edge
+ for (unsigned int face_n = 0; face_n< GeometryInfo<dim>::faces_per_cell; ++face_n)
+ {
+ MGDoFHandler<dim>::face_iterator face = c->face(face_n);
+ if (face->has_children())
+ {
+ face->get_mg_dof_indices(mgfaceindex);
+ for (unsigned int i=0;i<face_dofs;++i)
+ d[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, d[level], d[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);
+ }
}
template <int dim> template <typename number>
vector<int> index(fe_dofs);
vector<int> mgindex(fe_dofs);
+ for (unsigned int i=minlevel; i<= maxlevel; ++i)
+ {
+ char name[10];
+ sprintf(name,"CF%d",i);
+ ofstream of(name);
+ 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)
+ ; 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]);
- }
+ dst(index[i]) = s[c->level()](mgindex[i]);
+ }
+ hanging_nodes.set_zero(dst);
+ {
+ ofstream of("CF");
+ DataOut<2> out;
+ out.attach_dof_handler(*dofs);
+ out.add_data_vector(dst,"solution","m");
+ out.write_gnuplot(of,1);
+ }
}
template <int dim>
void
-MG<dim>::level_residual(unsigned int l,
- Vector<float> &r,
- const Vector<float> &u,
- const Vector<float> &b)
+MG<dim>::level_vmult(unsigned int l,
+ Vector<double> &r,
+ const Vector<double> &u,
+ const Vector<double> &)
{
- (*matrices)[l].residual(r,u,b);
+ (*matrices)[l].vmult(r,u);
+ r.scale(-1.);
}
/*
template void ConstraintMatrix::condense(const Vector<number> &uncondensed,
Vector<number> &condensed) const;
template void ConstraintMatrix::condense(Vector<number> &vec) const;
+template void ConstraintMatrix::set_zero(Vector<number> &vec) const;
template void ConstraintMatrix::distribute(const Vector<number> &condensed,
Vector<number> &uncondensed) const;
template void ConstraintMatrix::distribute(Vector<number> &vec) const;
template void ConstraintMatrix::condense(const Vector<number> &uncondensed,
Vector<number> &condensed) const;
template void ConstraintMatrix::condense(Vector<number> &vec) const;
+template void ConstraintMatrix::set_zero(Vector<number> &vec) const;
template void ConstraintMatrix::distribute(const Vector<number> &condensed,
Vector<number> &uncondensed) const;
template void ConstraintMatrix::distribute(Vector<number> &vec) const;
+
+#undef number
pre_smooth.smooth(level, s[level], d[level]);
// t = d-As
t.reinit(s[level].size());
- level_residual(level, t, s[level], d[level]);
+ level_vmult(level, t, s[level], d[level]);
// make t rhs of lower level
transfer->restrict(level, d[level-1], t);
void
MGSmootherIdentity::smooth (const unsigned int,
- Vector<float> &,
- const Vector<float> &) const
+ Vector<double> &,
+ const Vector<double> &) const
{}
void
MGSmoother::set_zero_interior_boundary (const unsigned int level,
- Vector<float> &u) const
+ Vector<double> &u) const
{
if (level==0)
return;
//////////////////////////////////////////////////////////////////////
-template<int dim>
-MGSmootherSOR::MGSmootherSOR(const MGDoFHandler<dim> &mg_dof,
- const MGMatrix<SparseMatrix<float> >& matrix,
- unsigned int steps)
- :
- MGSmoother(mg_dof, steps),
- matrix(&matrix)
-{}
-
-void
-MGSmootherSOR::smooth (const unsigned int /*level*/,
- Vector<float> &/*u*/,
- const Vector<float> &/*rhs*/) const
-{
-
-}
+#include <numerics/mg_smoother.templates.h>
// explicit instantiations
// don't do the following instantiation in 1d, since there is a specialized
template MGSmoother::MGSmoother (const MGDoFHandler<deal_II_dimension>&, unsigned int);
#endif
-template MGSmootherSOR::MGSmootherSOR(const MGDoFHandler<deal_II_dimension> &mg_dof,
- const MGMatrix<SparseMatrix<float> >& matrix,
- unsigned int steps);
+template
+MGSmootherRelaxation<float>
+::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension> &mg_dof,
+ const MGMatrix<SparseMatrix<float> >& matrix,
+ function_ptr relaxation,
+ unsigned int steps,
+ double omega);
+template
+MGSmootherRelaxation<double>
+::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension> &mg_dof,
+ const MGMatrix<SparseMatrix<double> >& matrix,
+ function_ptr relaxation,
+ unsigned int steps,
+ double omega);
#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);
/* ----------------------------- MGTransferPrebuilt ------------------------ */
// cell
prolongation_sparsities.back().reinit (mg_dof.n_dofs(level+1),
mg_dof.n_dofs(level),
- dofs_per_cell);
+// evil hack, must be corrected!!!
+ dofs_per_cell+1);
for (typename MGDoFHandler<dim>::cell_iterator cell=mg_dof.begin(level);
cell != mg_dof.end(level); ++cell)
void MGTransferPrebuilt::prolongate (const unsigned int to_level,
- Vector<float> &dst,
- const Vector<float> &src) const
+ Vector<double> &dst,
+ const Vector<double> &src) const
{
Assert ((to_level >= 1) && (to_level<=prolongation_matrices.size()),
ExcIndexRange (to_level, 1, prolongation_matrices.size()+1));
void MGTransferPrebuilt::restrict (const unsigned int from_level,
- Vector<float> &dst,
- const Vector<float> &src) const
+ 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));
-
- prolongation_matrices[from_level-1].Tvmult (dst, src);
+// smoother.set_zero_interior_boundary(from_level, src);
+ prolongation_matrices[from_level-1].Tvmult_add (dst, src);
};
void
MGSmoother::set_zero_interior_boundary (const unsigned int level,
- Vector<float> &u) const
+ Vector<double> &u) const
{
if (level==0)
return;
//////////////////////////////////////////////////////////////////////
-template<int dim>
-MGSmootherSOR::MGSmootherSOR(const MGDoFHandler<dim> &mg_dof,
- const MGMatrix<SparseMatrix<float> >& matrix,
- unsigned int steps)
- :
- MGSmoother(mg_dof, steps),
- matrix(&matrix)
-{}
-
-void
-MGSmootherSOR::smooth (const unsigned int /*level*/,
- Vector<float> &/*u*/,
- const Vector<float> &/*rhs*/) const
-{
-
-}
+#include <numerics/mg_smoother.templates.h>
// explicit instantiations
// don't do the following instantiation in 1d, since there is a specialized
template MGSmoother::MGSmoother (const MGDoFHandler<deal_II_dimension>&, unsigned int);
#endif
-template MGSmootherSOR::MGSmootherSOR(const MGDoFHandler<deal_II_dimension> &mg_dof,
- const MGMatrix<SparseMatrix<float> >& matrix,
- unsigned int steps);
+template
+MGSmootherRelaxation<float>
+::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension> &mg_dof,
+ const MGMatrix<SparseMatrix<float> >& matrix,
+ function_ptr relaxation,
+ unsigned int steps,
+ double omega);
+template
+MGSmootherRelaxation<double>
+::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension> &mg_dof,
+ const MGMatrix<SparseMatrix<double> >& matrix,
+ function_ptr relaxation,
+ unsigned int steps,
+ double omega);
#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);
/* ----------------------------- MGTransferPrebuilt ------------------------ */
// cell
prolongation_sparsities.back().reinit (mg_dof.n_dofs(level+1),
mg_dof.n_dofs(level),
- dofs_per_cell);
+// evil hack, must be corrected!!!
+ dofs_per_cell+1);
for (typename MGDoFHandler<dim>::cell_iterator cell=mg_dof.begin(level);
cell != mg_dof.end(level); ++cell)
void MGTransferPrebuilt::prolongate (const unsigned int to_level,
- Vector<float> &dst,
- const Vector<float> &src) const
+ Vector<double> &dst,
+ const Vector<double> &src) const
{
Assert ((to_level >= 1) && (to_level<=prolongation_matrices.size()),
ExcIndexRange (to_level, 1, prolongation_matrices.size()+1));
void MGTransferPrebuilt::restrict (const unsigned int from_level,
- Vector<float> &dst,
- const Vector<float> &src) const
+ 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));
-
- prolongation_matrices[from_level-1].Tvmult (dst, src);
+// smoother.set_zero_interior_boundary(from_level, src);
+ prolongation_matrices[from_level-1].Tvmult_add (dst, src);
};
* about the matrix is removed to
* that class.
*/
- virtual void operator() (unsigned int level, Vector<float>& dst,
- const Vector<float>& src) const = 0;
+ virtual void operator() (unsigned int level, Vector<double>& dst,
+ const Vector<double>& src) const = 0;
};
* things.
*/
virtual void smooth (const unsigned int level,
- Vector<float> &u,
- const Vector<float> &rhs) const = 0;
+ Vector<double> &u,
+ const Vector<double> &rhs) const = 0;
};
* This function does nothing.
*/
virtual void smooth (const unsigned int level,
- Vector<float> &u,
- const Vector<float> &rhs) const;
+ Vector<double> &u,
+ const Vector<double> &rhs) const;
};
/**
* level.
*/
virtual void prolongate (const unsigned int to_level,
- Vector<float> &dst,
- const Vector<float> &src) const = 0;
+ Vector<double> &dst,
+ const Vector<double> &src) const = 0;
/**
* Restrict a vector from level
* level.
*/
virtual void restrict (const unsigned int from_level,
- Vector<float> &dst,
- const Vector<float> &src) const = 0;
+ Vector<double> &dst,
+ const Vector<double> &src) const = 0;
};
/**
* Safe access operator.
*/
const VECTOR& operator[](unsigned int) const;
+ /**
+ * Delete all entries.
+ */
+ void clear();
private:
/**
* Level of first component.
* matrix, vectors and
* preconditioner.
*/
- virtual void operator() (unsigned int level, Vector<float>& dst,
- const Vector<float>& src) const;
+ virtual void operator() (unsigned int level, Vector<double>& dst,
+ const Vector<double>& src) const;
private:
/**
* Reference to the solver.
/**
* Auxiliary vector, defect.
*/
- MGVector<Vector<float> > d;
+ MGVector<Vector<double> > d;
/**
* Auxiliary vector, solution.
*/
- MGVector<Vector<float> > s;
+ MGVector<Vector<double> > s;
+ /**
+ * Prolongation and restriction object.
+ */
+ SmartPointer<const MGTransferBase> transfer;
private:
/**
* Auxiliary vector.
*/
- Vector<float> t;
+ Vector<double> t;
- /**
- * Prolongation and restriction object.
- */
- SmartPointer<const MGTransferBase> transfer;
/**
* Exception.
*/
const MGCoarseGridSolver& cgs);
/**
- * Apply residual operator on all
+ * Apply negative #vmult# on all
* cells of a level.
* This is implemented in a
* derived class.
*/
- virtual void level_residual(unsigned int level,
- Vector<float>& dst,
- const Vector<float>& src,
- const Vector<float>& rhs) = 0;
+ virtual void level_vmult(unsigned int level,
+ Vector<double>& dst,
+ const Vector<double>& src,
+ const Vector<double>& rhs) = 0;
public:
template<class SOLVER, class MATRIX, class PRECOND>
void
MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND>::operator()
-(unsigned int, Vector<float>& dst, const Vector<float>& src) const
+(unsigned int, Vector<double>& dst, const Vector<double>& src) const
{
solver.solve(matrix, dst, src, precondition);
}
return vector<VECTOR>::operator[](i-minlevel);
}
+template<class VECTOR>
+void
+MGVector<VECTOR>::clear()
+{
+ vector<VECTOR>::iterator v;
+ for (v = begin(); v != end(); ++v)
+ v->clear();
+
+}
/* ------------------------------------------------------------------- */
// and since it makes the backward
// step when applying the decomposition
// significantly faster
- Assert((global_entry(rowstart_indices[row-1]) !=0), ExcDivideByZero());
+ AssertThrow((global_entry(rowstart_indices[row-1]) !=0),
+ ExcDivideByZero());
global_entry (rowstart_indices[row-1])
= 1./global_entry (rowstart_indices[row-1]);
pre_smooth.smooth(level, s[level], d[level]);
// t = d-As
t.reinit(s[level].size());
- level_residual(level, t, s[level], d[level]);
+ level_vmult(level, t, s[level], d[level]);
// make t rhs of lower level
transfer->restrict(level, d[level-1], t);
void
MGSmootherIdentity::smooth (const unsigned int,
- Vector<float> &,
- const Vector<float> &) const
+ Vector<double> &,
+ const Vector<double> &) const
{}