From: guido Date: Mon, 28 Jun 1999 15:19:09 +0000 (+0000) Subject: MG works on locally refined meshes, needs improvement and cleaning X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7d902ad83fdfbb874b4c482fa45600dc543416b4;p=dealii-svn.git MG works on locally refined meshes, needs improvement and cleaning git-svn-id: https://svn.dealii.org/trunk@1486 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/dofs/dof_constraints.h b/deal.II/deal.II/include/dofs/dof_constraints.h index 92a14f7bc8..34fb1ab7d8 100644 --- a/deal.II/deal.II/include/dofs/dof_constraints.h +++ b/deal.II/deal.II/include/dofs/dof_constraints.h @@ -247,6 +247,13 @@ class ConstraintMatrix { template void distribute (Vector &vec) const; + /** + * Delete hanging nodes in a vector. + * Sets all hanging node values to zero. + */ + template + void set_zero (Vector &vec) const; + /** * Print the constraint lines. Mainly for diff --git a/deal.II/deal.II/include/dofs/dof_constraints.templates.h b/deal.II/deal.II/include/dofs/dof_constraints.templates.h index 4e964898b8..5c2f975c32 100644 --- a/deal.II/deal.II/include/dofs/dof_constraints.templates.h +++ b/deal.II/deal.II/include/dofs/dof_constraints.templates.h @@ -348,6 +348,32 @@ ConstraintMatrix::condense (Vector &vec) const +template +void +ConstraintMatrix::set_zero (Vector &vec) const +{ + Assert (sorted == true, ExcMatrixNotClosed()); + + if (lines.size() == 0) + // nothing to do + return; + + vector::const_iterator next_constraint = lines.begin(); + for (unsigned int row=0; rowline) + { + // set entry to zero + vec(row) = 0.; + + ++next_constraint; + if (next_constraint == lines.end()) + // nothing more to do + break; + }; +}; + + + template void ConstraintMatrix::distribute (const Vector &condensed, diff --git a/deal.II/deal.II/include/multigrid/mg_base.h b/deal.II/deal.II/include/multigrid/mg_base.h index 18704a4799..fd6a230c29 100644 --- a/deal.II/deal.II/include/multigrid/mg_base.h +++ b/deal.II/deal.II/include/multigrid/mg_base.h @@ -40,8 +40,8 @@ class MGCoarseGridSolver * about the matrix is removed to * that class. */ - virtual void operator() (unsigned int level, Vector& dst, - const Vector& src) const = 0; + virtual void operator() (unsigned int level, Vector& dst, + const Vector& src) const = 0; }; @@ -74,8 +74,8 @@ class MGSmootherBase * things. */ virtual void smooth (const unsigned int level, - Vector &u, - const Vector &rhs) const = 0; + Vector &u, + const Vector &rhs) const = 0; }; @@ -93,8 +93,8 @@ class MGSmootherIdentity * This function does nothing. */ virtual void smooth (const unsigned int level, - Vector &u, - const Vector &rhs) const; + Vector &u, + const Vector &rhs) const; }; /** @@ -129,8 +129,8 @@ class MGTransferBase * level. */ virtual void prolongate (const unsigned int to_level, - Vector &dst, - const Vector &src) const = 0; + Vector &dst, + const Vector &src) const = 0; /** * Restrict a vector from level @@ -146,8 +146,8 @@ class MGTransferBase * level. */ virtual void restrict (const unsigned int from_level, - Vector &dst, - const Vector &src) const = 0; + Vector &dst, + const Vector &src) const = 0; }; /** @@ -175,6 +175,10 @@ class MGVector * Safe access operator. */ const VECTOR& operator[](unsigned int) const; + /** + * Delete all entries. + */ + void clear(); private: /** * Level of first component. @@ -245,8 +249,8 @@ class MGCoarseGridLACIteration * matrix, vectors and * preconditioner. */ - virtual void operator() (unsigned int level, Vector& dst, - const Vector& src) const; + virtual void operator() (unsigned int level, Vector& dst, + const Vector& src) const; private: /** * Reference to the solver. @@ -303,23 +307,23 @@ class MGBase /** * Auxiliary vector, defect. */ - MGVector > d; + MGVector > d; /** * Auxiliary vector, solution. */ - MGVector > s; + MGVector > s; + /** + * Prolongation and restriction object. + */ + SmartPointer transfer; private: /** * Auxiliary vector. */ - Vector t; + Vector t; - /** - * Prolongation and restriction object. - */ - SmartPointer transfer; /** * Exception. */ @@ -344,15 +348,15 @@ class MGBase 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& dst, - const Vector& src, - const Vector& rhs) = 0; + virtual void level_vmult(unsigned int level, + Vector& dst, + const Vector& src, + const Vector& rhs) = 0; public: @@ -450,7 +454,7 @@ MGCoarseGridLACIteration(SOLVER& s, const MATRIX& m, const PRECOND& p) template void MGCoarseGridLACIteration::operator() -(unsigned int, Vector& dst, const Vector& src) const +(unsigned int, Vector& dst, const Vector& src) const { solver.solve(matrix, dst, src, precondition); } @@ -485,6 +489,15 @@ MGVector::operator[](unsigned int i) const return vector::operator[](i-minlevel); } +template +void +MGVector::clear() +{ + vector::iterator v; + for (v = begin(); v != end(); ++v) + v->clear(); + +} /* ------------------------------------------------------------------- */ diff --git a/deal.II/deal.II/include/multigrid/mg_smoother.h b/deal.II/deal.II/include/multigrid/mg_smoother.h index 6db42e8039..110cca5eb1 100644 --- a/deal.II/deal.II/include/multigrid/mg_smoother.h +++ b/deal.II/deal.II/include/multigrid/mg_smoother.h @@ -22,8 +22,7 @@ * * @author Wolfgang Bangerth, Guido Kanschat, 1999 */ -class MGSmoother - : public MGSmootherBase +class MGSmoother : public MGSmootherBase { private: /** @@ -82,7 +81,7 @@ class MGSmoother * function does nothing in this case. */ void set_zero_interior_boundary (const unsigned int level, - Vector &u) const; + Vector &u) const; /** * Modify the number of smoothing steps. */ @@ -113,31 +112,45 @@ class MGSmoother }; /** - * 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 +class MGSmootherRelaxation : public MGSmoother { public: + /** + * Type of the smoothing + * function of the matrix. + */ + typedef void + (SparseMatrix::* function_ptr)(Vector&, const Vector&, + typename SparseMatrix::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::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 - MGSmootherSOR(const MGDoFHandler &mg_dof, - const MGMatrix >& matrix, - unsigned int steps); + MGSmootherRelaxation(const MGDoFHandler &mg_dof, + const MGMatrix >& matrix, + function_ptr function, + unsigned int steps, + double omega = 1.); /** * Implementation of the interface in #MGSmootherBase#. @@ -145,13 +158,25 @@ class MGSmootherSOR * and find a way to keep the boundary values. */ virtual void smooth (const unsigned int level, - Vector &u, - const Vector &rhs) const; + Vector &u, + const Vector &rhs) const; private: /** * Pointer to the matrices. */ - SmartPointer< const MGMatrix< SparseMatrix > > matrix; + SmartPointer< const MGMatrix< SparseMatrix > > matrix; + /** + * Pointer to the relaxation function. + */ + function_ptr relaxation; + /** + * Relaxation parameter. + */ + double omega; + /** + * Auxiliary vector. + */ + mutable Vector h1, h2; }; inline diff --git a/deal.II/deal.II/include/multigrid/mg_smoother.templates.h b/deal.II/deal.II/include/multigrid/mg_smoother.templates.h new file mode 100644 index 0000000000..f72af10b02 --- /dev/null +++ b/deal.II/deal.II/include/multigrid/mg_smoother.templates.h @@ -0,0 +1,34 @@ +// $Id$ + +template +template +MGSmootherRelaxation:: +MGSmootherRelaxation(const MGDoFHandler &mg_dof, + const MGMatrix >& matrix, + function_ptr relaxation, + unsigned int steps, + double omega) + : + MGSmoother(mg_dof, steps), + matrix(&matrix), + relaxation(relaxation), + omega(omega) +{} + +template +void +MGSmootherRelaxation::smooth (const unsigned int level, + Vector &u, + const Vector &rhs) const +{ + h1.reinit(u); + h2.reinit(u); + for(unsigned i=0;i&, - const MGMatrix >&, - const MGTransferBase& transfer); + ConstraintMatrix& hanging_nodes, + const MGMatrix >&, + const MGTransferBase& transfer, + unsigned int minlevel = 0, unsigned int maxlevel = 10000); /** Transfer from dVector to * MGVector. @@ -58,17 +60,18 @@ class MG void copy_from_mg(Vector& dst) const; /** - * Residual on a level. + * Negative #vmult# on a level. + * @see MGBase. */ - virtual void level_residual(unsigned int, - Vector &, - const Vector &, - const Vector &); + virtual void level_vmult(unsigned int, + Vector &, + const Vector &, + const Vector &); /** * Read-only access to level matrices. */ - const SparseMatrix& get_matrix(unsigned int level) const; + const SparseMatrix& get_matrix(unsigned int level) const; private: /** @@ -82,7 +85,8 @@ class MG * the constructor of #MG# and can * be accessed for assembling. */ - SmartPointer > > matrices; + SmartPointer > > matrices; + ConstraintMatrix& hanging_nodes; }; @@ -97,8 +101,16 @@ class MG */ class MGTransferPrebuilt : public MGTransferBase { + const MGSmoother& smoother; public: - + /** + * Preliminary constructor + */ + MGTransferPrebuilt(const MGSmoother& smoother) : + smoother(smoother) + {} + + /** * Destructor. */ @@ -124,8 +136,8 @@ class MGTransferPrebuilt : public MGTransferBase * level. */ virtual void prolongate (const unsigned int to_level, - Vector &dst, - const Vector &src) const; + Vector &dst, + const Vector &src) const; /** * Restrict a vector from level @@ -141,8 +153,8 @@ class MGTransferPrebuilt : public MGTransferBase * level. */ virtual void restrict (const unsigned int from_level, - Vector &dst, - const Vector &src) const; + Vector &dst, + const Vector &src) const; private: @@ -162,7 +174,7 @@ class MGTransferPrebuilt : public MGTransferBase template inline -const SparseMatrix& +const SparseMatrix& MG::get_matrix(unsigned int level) const { Assert((level>=minlevel) && (level +#include + +inline int minimum(int a, int b) +{ + return ((a MG::MG(const MGDoFHandler& dofs, - const MGMatrix >& matrices, - const MGTransferBase& transfer) + ConstraintMatrix& hanging_nodes, + const MGMatrix >& 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) { @@ -24,19 +34,59 @@ void MG::copy_to_mg(const Vector& 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 index(fe_dofs); vector mgindex(fe_dofs); + vector 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::active_cell_iterator dc = dofs->DoFHandler::begin_active(); - for (MGDoFHandler::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::active_cell_iterator dc = dofs->DoFHandler::begin_active(level); + MGDoFHandler::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;ilevel()](mgindex[i]) = src(index[i]); + + // Delete values on refinement edge + for (unsigned int face_n = 0; face_n< GeometryInfo::faces_per_cell; ++face_n) + { + MGDoFHandler::face_iterator face = c->face(face_n); + if (face->has_children()) + { + face->get_mg_dof_indices(mgfaceindex); + for (unsigned int i=0;ilevel()](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 template @@ -48,25 +98,42 @@ MG::copy_from_mg(Vector& dst) const vector index(fe_dofs); vector 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::active_cell_iterator dc = dofs->DoFHandler::begin_active(); for (MGDoFHandler::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;ilevel()](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 void -MG::level_residual(unsigned int l, - Vector &r, - const Vector &u, - const Vector &b) +MG::level_vmult(unsigned int l, + Vector &r, + const Vector &u, + const Vector &) { - (*matrices)[l].residual(r,u,b); + (*matrices)[l].vmult(r,u); + r.scale(-1.); } /* diff --git a/deal.II/deal.II/include/numerics/data_io.h b/deal.II/deal.II/include/numerics/data_io.h index 538f4428af..75b6ee8c13 100644 --- a/deal.II/deal.II/include/numerics/data_io.h +++ b/deal.II/deal.II/include/numerics/data_io.h @@ -809,14 +809,14 @@ class EpsOutputData{ * 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. diff --git a/deal.II/deal.II/include/numerics/mg_smoother.h b/deal.II/deal.II/include/numerics/mg_smoother.h index 6db42e8039..110cca5eb1 100644 --- a/deal.II/deal.II/include/numerics/mg_smoother.h +++ b/deal.II/deal.II/include/numerics/mg_smoother.h @@ -22,8 +22,7 @@ * * @author Wolfgang Bangerth, Guido Kanschat, 1999 */ -class MGSmoother - : public MGSmootherBase +class MGSmoother : public MGSmootherBase { private: /** @@ -82,7 +81,7 @@ class MGSmoother * function does nothing in this case. */ void set_zero_interior_boundary (const unsigned int level, - Vector &u) const; + Vector &u) const; /** * Modify the number of smoothing steps. */ @@ -113,31 +112,45 @@ class MGSmoother }; /** - * 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 +class MGSmootherRelaxation : public MGSmoother { public: + /** + * Type of the smoothing + * function of the matrix. + */ + typedef void + (SparseMatrix::* function_ptr)(Vector&, const Vector&, + typename SparseMatrix::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::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 - MGSmootherSOR(const MGDoFHandler &mg_dof, - const MGMatrix >& matrix, - unsigned int steps); + MGSmootherRelaxation(const MGDoFHandler &mg_dof, + const MGMatrix >& matrix, + function_ptr function, + unsigned int steps, + double omega = 1.); /** * Implementation of the interface in #MGSmootherBase#. @@ -145,13 +158,25 @@ class MGSmootherSOR * and find a way to keep the boundary values. */ virtual void smooth (const unsigned int level, - Vector &u, - const Vector &rhs) const; + Vector &u, + const Vector &rhs) const; private: /** * Pointer to the matrices. */ - SmartPointer< const MGMatrix< SparseMatrix > > matrix; + SmartPointer< const MGMatrix< SparseMatrix > > matrix; + /** + * Pointer to the relaxation function. + */ + function_ptr relaxation; + /** + * Relaxation parameter. + */ + double omega; + /** + * Auxiliary vector. + */ + mutable Vector h1, h2; }; inline diff --git a/deal.II/deal.II/include/numerics/mg_smoother.templates.h b/deal.II/deal.II/include/numerics/mg_smoother.templates.h new file mode 100644 index 0000000000..f72af10b02 --- /dev/null +++ b/deal.II/deal.II/include/numerics/mg_smoother.templates.h @@ -0,0 +1,34 @@ +// $Id$ + +template +template +MGSmootherRelaxation:: +MGSmootherRelaxation(const MGDoFHandler &mg_dof, + const MGMatrix >& matrix, + function_ptr relaxation, + unsigned int steps, + double omega) + : + MGSmoother(mg_dof, steps), + matrix(&matrix), + relaxation(relaxation), + omega(omega) +{} + +template +void +MGSmootherRelaxation::smooth (const unsigned int level, + Vector &u, + const Vector &rhs) const +{ + h1.reinit(u); + h2.reinit(u); + for(unsigned i=0;i&, - const MGMatrix >&, - const MGTransferBase& transfer); + ConstraintMatrix& hanging_nodes, + const MGMatrix >&, + const MGTransferBase& transfer, + unsigned int minlevel = 0, unsigned int maxlevel = 10000); /** Transfer from dVector to * MGVector. @@ -58,17 +60,18 @@ class MG void copy_from_mg(Vector& dst) const; /** - * Residual on a level. + * Negative #vmult# on a level. + * @see MGBase. */ - virtual void level_residual(unsigned int, - Vector &, - const Vector &, - const Vector &); + virtual void level_vmult(unsigned int, + Vector &, + const Vector &, + const Vector &); /** * Read-only access to level matrices. */ - const SparseMatrix& get_matrix(unsigned int level) const; + const SparseMatrix& get_matrix(unsigned int level) const; private: /** @@ -82,7 +85,8 @@ class MG * the constructor of #MG# and can * be accessed for assembling. */ - SmartPointer > > matrices; + SmartPointer > > matrices; + ConstraintMatrix& hanging_nodes; }; @@ -97,8 +101,16 @@ class MG */ class MGTransferPrebuilt : public MGTransferBase { + const MGSmoother& smoother; public: - + /** + * Preliminary constructor + */ + MGTransferPrebuilt(const MGSmoother& smoother) : + smoother(smoother) + {} + + /** * Destructor. */ @@ -124,8 +136,8 @@ class MGTransferPrebuilt : public MGTransferBase * level. */ virtual void prolongate (const unsigned int to_level, - Vector &dst, - const Vector &src) const; + Vector &dst, + const Vector &src) const; /** * Restrict a vector from level @@ -141,8 +153,8 @@ class MGTransferPrebuilt : public MGTransferBase * level. */ virtual void restrict (const unsigned int from_level, - Vector &dst, - const Vector &src) const; + Vector &dst, + const Vector &src) const; private: @@ -162,7 +174,7 @@ class MGTransferPrebuilt : public MGTransferBase template inline -const SparseMatrix& +const SparseMatrix& MG::get_matrix(unsigned int level) const { Assert((level>=minlevel) && (level +#include + +inline int minimum(int a, int b) +{ + return ((a MG::MG(const MGDoFHandler& dofs, - const MGMatrix >& matrices, - const MGTransferBase& transfer) + ConstraintMatrix& hanging_nodes, + const MGMatrix >& 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) { @@ -24,19 +34,59 @@ void MG::copy_to_mg(const Vector& 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 index(fe_dofs); vector mgindex(fe_dofs); + vector 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::active_cell_iterator dc = dofs->DoFHandler::begin_active(); - for (MGDoFHandler::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::active_cell_iterator dc = dofs->DoFHandler::begin_active(level); + MGDoFHandler::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;ilevel()](mgindex[i]) = src(index[i]); + + // Delete values on refinement edge + for (unsigned int face_n = 0; face_n< GeometryInfo::faces_per_cell; ++face_n) + { + MGDoFHandler::face_iterator face = c->face(face_n); + if (face->has_children()) + { + face->get_mg_dof_indices(mgfaceindex); + for (unsigned int i=0;ilevel()](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 template @@ -48,25 +98,42 @@ MG::copy_from_mg(Vector& dst) const vector index(fe_dofs); vector 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::active_cell_iterator dc = dofs->DoFHandler::begin_active(); for (MGDoFHandler::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;ilevel()](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 void -MG::level_residual(unsigned int l, - Vector &r, - const Vector &u, - const Vector &b) +MG::level_vmult(unsigned int l, + Vector &r, + const Vector &u, + const Vector &) { - (*matrices)[l].residual(r,u,b); + (*matrices)[l].vmult(r,u); + r.scale(-1.); } /* diff --git a/deal.II/deal.II/source/dofs/dof_constraints.cc b/deal.II/deal.II/source/dofs/dof_constraints.cc index 1af3e7944b..e66b41e8cf 100644 --- a/deal.II/deal.II/source/dofs/dof_constraints.cc +++ b/deal.II/deal.II/source/dofs/dof_constraints.cc @@ -332,6 +332,7 @@ template void ConstraintMatrix::condense(SparseMatrix &uncondensed) cons template void ConstraintMatrix::condense(const Vector &uncondensed, Vector &condensed) const; template void ConstraintMatrix::condense(Vector &vec) const; +template void ConstraintMatrix::set_zero(Vector &vec) const; template void ConstraintMatrix::distribute(const Vector &condensed, Vector &uncondensed) const; template void ConstraintMatrix::distribute(Vector &vec) const; @@ -345,6 +346,9 @@ template void ConstraintMatrix::condense(SparseMatrix &uncondensed) cons template void ConstraintMatrix::condense(const Vector &uncondensed, Vector &condensed) const; template void ConstraintMatrix::condense(Vector &vec) const; +template void ConstraintMatrix::set_zero(Vector &vec) const; template void ConstraintMatrix::distribute(const Vector &condensed, Vector &uncondensed) const; template void ConstraintMatrix::distribute(Vector &vec) const; + +#undef number diff --git a/deal.II/deal.II/source/multigrid/mg_base.cc b/deal.II/deal.II/source/multigrid/mg_base.cc index e6340762e8..9df47ff594 100644 --- a/deal.II/deal.II/source/multigrid/mg_base.cc +++ b/deal.II/deal.II/source/multigrid/mg_base.cc @@ -63,7 +63,7 @@ MGBase::level_mgstep(unsigned level, 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); @@ -100,8 +100,8 @@ MGSmootherBase::~MGSmootherBase() void MGSmootherIdentity::smooth (const unsigned int, - Vector &, - const Vector &) const + Vector &, + const Vector &) const {} diff --git a/deal.II/deal.II/source/multigrid/mg_smoother.cc b/deal.II/deal.II/source/multigrid/mg_smoother.cc index 23bea7acf9..9fad33f9c8 100644 --- a/deal.II/deal.II/source/multigrid/mg_smoother.cc +++ b/deal.II/deal.II/source/multigrid/mg_smoother.cc @@ -92,7 +92,7 @@ MGSmoother::MGSmoother (const MGDoFHandler &mg_dof, unsigned int steps) void MGSmoother::set_zero_interior_boundary (const unsigned int level, - Vector &u) const + Vector &u) const { if (level==0) return; @@ -104,22 +104,7 @@ MGSmoother::set_zero_interior_boundary (const unsigned int level, ////////////////////////////////////////////////////////////////////// -template -MGSmootherSOR::MGSmootherSOR(const MGDoFHandler &mg_dof, - const MGMatrix >& matrix, - unsigned int steps) - : - MGSmoother(mg_dof, steps), - matrix(&matrix) -{} - -void -MGSmootherSOR::smooth (const unsigned int /*level*/, - Vector &/*u*/, - const Vector &/*rhs*/) const -{ - -} +#include // explicit instantiations // don't do the following instantiation in 1d, since there is a specialized @@ -128,7 +113,18 @@ MGSmootherSOR::smooth (const unsigned int /*level*/, template MGSmoother::MGSmoother (const MGDoFHandler&, unsigned int); #endif -template MGSmootherSOR::MGSmootherSOR(const MGDoFHandler &mg_dof, - const MGMatrix >& matrix, - unsigned int steps); +template +MGSmootherRelaxation +::MGSmootherRelaxation(const MGDoFHandler &mg_dof, + const MGMatrix >& matrix, + function_ptr relaxation, + unsigned int steps, + double omega); +template +MGSmootherRelaxation +::MGSmootherRelaxation(const MGDoFHandler &mg_dof, + const MGMatrix >& matrix, + function_ptr relaxation, + unsigned int steps, + double omega); diff --git a/deal.II/deal.II/source/multigrid/multigrid.cc b/deal.II/deal.II/source/multigrid/multigrid.cc index 429f929615..2785d302c6 100644 --- a/deal.II/deal.II/source/multigrid/multigrid.cc +++ b/deal.II/deal.II/source/multigrid/multigrid.cc @@ -11,6 +11,8 @@ #include #include +extern void write_gnuplot (const MGDoFHandler<2>& dofs, const Vector& v, unsigned int level, + ostream &out, unsigned int accuracy); /* ----------------------------- MGTransferPrebuilt ------------------------ */ @@ -59,7 +61,8 @@ void MGTransferPrebuilt::build_matrices (const MGDoFHandler &mg_dof) // 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::cell_iterator cell=mg_dof.begin(level); cell != mg_dof.end(level); ++cell) @@ -134,8 +137,8 @@ void MGTransferPrebuilt::build_matrices (const MGDoFHandler &mg_dof) void MGTransferPrebuilt::prolongate (const unsigned int to_level, - Vector &dst, - const Vector &src) const + Vector &dst, + const Vector &src) const { Assert ((to_level >= 1) && (to_level<=prolongation_matrices.size()), ExcIndexRange (to_level, 1, prolongation_matrices.size()+1)); @@ -146,13 +149,13 @@ void MGTransferPrebuilt::prolongate (const unsigned int to_level, void MGTransferPrebuilt::restrict (const unsigned int from_level, - Vector &dst, - const Vector &src) const + Vector &dst, + const Vector &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); }; diff --git a/deal.II/deal.II/source/numerics/mg_smoother.cc b/deal.II/deal.II/source/numerics/mg_smoother.cc index 23bea7acf9..9fad33f9c8 100644 --- a/deal.II/deal.II/source/numerics/mg_smoother.cc +++ b/deal.II/deal.II/source/numerics/mg_smoother.cc @@ -92,7 +92,7 @@ MGSmoother::MGSmoother (const MGDoFHandler &mg_dof, unsigned int steps) void MGSmoother::set_zero_interior_boundary (const unsigned int level, - Vector &u) const + Vector &u) const { if (level==0) return; @@ -104,22 +104,7 @@ MGSmoother::set_zero_interior_boundary (const unsigned int level, ////////////////////////////////////////////////////////////////////// -template -MGSmootherSOR::MGSmootherSOR(const MGDoFHandler &mg_dof, - const MGMatrix >& matrix, - unsigned int steps) - : - MGSmoother(mg_dof, steps), - matrix(&matrix) -{} - -void -MGSmootherSOR::smooth (const unsigned int /*level*/, - Vector &/*u*/, - const Vector &/*rhs*/) const -{ - -} +#include // explicit instantiations // don't do the following instantiation in 1d, since there is a specialized @@ -128,7 +113,18 @@ MGSmootherSOR::smooth (const unsigned int /*level*/, template MGSmoother::MGSmoother (const MGDoFHandler&, unsigned int); #endif -template MGSmootherSOR::MGSmootherSOR(const MGDoFHandler &mg_dof, - const MGMatrix >& matrix, - unsigned int steps); +template +MGSmootherRelaxation +::MGSmootherRelaxation(const MGDoFHandler &mg_dof, + const MGMatrix >& matrix, + function_ptr relaxation, + unsigned int steps, + double omega); +template +MGSmootherRelaxation +::MGSmootherRelaxation(const MGDoFHandler &mg_dof, + const MGMatrix >& matrix, + function_ptr relaxation, + unsigned int steps, + double omega); diff --git a/deal.II/deal.II/source/numerics/multigrid.cc b/deal.II/deal.II/source/numerics/multigrid.cc index 429f929615..2785d302c6 100644 --- a/deal.II/deal.II/source/numerics/multigrid.cc +++ b/deal.II/deal.II/source/numerics/multigrid.cc @@ -11,6 +11,8 @@ #include #include +extern void write_gnuplot (const MGDoFHandler<2>& dofs, const Vector& v, unsigned int level, + ostream &out, unsigned int accuracy); /* ----------------------------- MGTransferPrebuilt ------------------------ */ @@ -59,7 +61,8 @@ void MGTransferPrebuilt::build_matrices (const MGDoFHandler &mg_dof) // 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::cell_iterator cell=mg_dof.begin(level); cell != mg_dof.end(level); ++cell) @@ -134,8 +137,8 @@ void MGTransferPrebuilt::build_matrices (const MGDoFHandler &mg_dof) void MGTransferPrebuilt::prolongate (const unsigned int to_level, - Vector &dst, - const Vector &src) const + Vector &dst, + const Vector &src) const { Assert ((to_level >= 1) && (to_level<=prolongation_matrices.size()), ExcIndexRange (to_level, 1, prolongation_matrices.size()+1)); @@ -146,13 +149,13 @@ void MGTransferPrebuilt::prolongate (const unsigned int to_level, void MGTransferPrebuilt::restrict (const unsigned int from_level, - Vector &dst, - const Vector &src) const + Vector &dst, + const Vector &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); }; diff --git a/deal.II/lac/include/lac/mgbase.h b/deal.II/lac/include/lac/mgbase.h index d464b68f08..8e860139e7 100644 --- a/deal.II/lac/include/lac/mgbase.h +++ b/deal.II/lac/include/lac/mgbase.h @@ -40,8 +40,8 @@ class MGCoarseGridSolver * about the matrix is removed to * that class. */ - virtual void operator() (unsigned int level, Vector& dst, - const Vector& src) const = 0; + virtual void operator() (unsigned int level, Vector& dst, + const Vector& src) const = 0; }; @@ -74,8 +74,8 @@ class MGSmootherBase * things. */ virtual void smooth (const unsigned int level, - Vector &u, - const Vector &rhs) const = 0; + Vector &u, + const Vector &rhs) const = 0; }; @@ -93,8 +93,8 @@ class MGSmootherIdentity * This function does nothing. */ virtual void smooth (const unsigned int level, - Vector &u, - const Vector &rhs) const; + Vector &u, + const Vector &rhs) const; }; /** @@ -129,8 +129,8 @@ class MGTransferBase * level. */ virtual void prolongate (const unsigned int to_level, - Vector &dst, - const Vector &src) const = 0; + Vector &dst, + const Vector &src) const = 0; /** * Restrict a vector from level @@ -146,8 +146,8 @@ class MGTransferBase * level. */ virtual void restrict (const unsigned int from_level, - Vector &dst, - const Vector &src) const = 0; + Vector &dst, + const Vector &src) const = 0; }; /** @@ -175,6 +175,10 @@ class MGVector * Safe access operator. */ const VECTOR& operator[](unsigned int) const; + /** + * Delete all entries. + */ + void clear(); private: /** * Level of first component. @@ -245,8 +249,8 @@ class MGCoarseGridLACIteration * matrix, vectors and * preconditioner. */ - virtual void operator() (unsigned int level, Vector& dst, - const Vector& src) const; + virtual void operator() (unsigned int level, Vector& dst, + const Vector& src) const; private: /** * Reference to the solver. @@ -303,23 +307,23 @@ class MGBase /** * Auxiliary vector, defect. */ - MGVector > d; + MGVector > d; /** * Auxiliary vector, solution. */ - MGVector > s; + MGVector > s; + /** + * Prolongation and restriction object. + */ + SmartPointer transfer; private: /** * Auxiliary vector. */ - Vector t; + Vector t; - /** - * Prolongation and restriction object. - */ - SmartPointer transfer; /** * Exception. */ @@ -344,15 +348,15 @@ class MGBase 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& dst, - const Vector& src, - const Vector& rhs) = 0; + virtual void level_vmult(unsigned int level, + Vector& dst, + const Vector& src, + const Vector& rhs) = 0; public: @@ -450,7 +454,7 @@ MGCoarseGridLACIteration(SOLVER& s, const MATRIX& m, const PRECOND& p) template void MGCoarseGridLACIteration::operator() -(unsigned int, Vector& dst, const Vector& src) const +(unsigned int, Vector& dst, const Vector& src) const { solver.solve(matrix, dst, src, precondition); } @@ -485,6 +489,15 @@ MGVector::operator[](unsigned int i) const return vector::operator[](i-minlevel); } +template +void +MGVector::clear() +{ + vector::iterator v; + for (v = begin(); v != end(); ++v) + v->clear(); + +} /* ------------------------------------------------------------------- */ diff --git a/deal.II/lac/include/lac/sparse_ilu.templates.h b/deal.II/lac/include/lac/sparse_ilu.templates.h index 198697afb1..4f0acff54a 100644 --- a/deal.II/lac/include/lac/sparse_ilu.templates.h +++ b/deal.II/lac/include/lac/sparse_ilu.templates.h @@ -139,7 +139,8 @@ void SparseILU::decompose (const SparseMatrix &matrix, // 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]); diff --git a/deal.II/lac/source/mgbase.cc b/deal.II/lac/source/mgbase.cc index e6340762e8..9df47ff594 100644 --- a/deal.II/lac/source/mgbase.cc +++ b/deal.II/lac/source/mgbase.cc @@ -63,7 +63,7 @@ MGBase::level_mgstep(unsigned level, 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); @@ -100,8 +100,8 @@ MGSmootherBase::~MGSmootherBase() void MGSmootherIdentity::smooth (const unsigned int, - Vector &, - const Vector &) const + Vector &, + const Vector &) const {}