From: Wolfgang Bangerth Date: Wed, 26 Aug 1998 08:27:02 +0000 (+0000) Subject: Next round of multigrid implementations. X-Git-Tag: v8.0.0~22732 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d971130e53e3358309cb1ae2eaa026e89e6260c6;p=dealii.git Next round of multigrid implementations. git-svn-id: https://svn.dealii.org/trunk@522 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/dofs/mg_dof_accessor.h b/deal.II/deal.II/include/dofs/mg_dof_accessor.h index 03274755dc..75ba6a5a74 100644 --- a/deal.II/deal.II/include/dofs/mg_dof_accessor.h +++ b/deal.II/deal.II/include/dofs/mg_dof_accessor.h @@ -51,15 +51,22 @@ class MGDoFAccessor { /** * This should be the default constructor. + * We cast away the #const#ness of the + * pointer which clearly is EVIL but + * we can't help without making all + * functions which could somehow use + * iterators (directly or indirectly) make + * non-const, even if they preserve + * constness. */ - MGDoFAccessor (MGDoFHandler *mg_dof_handler) : - mg_dof_handler(mg_dof_handler) {}; + MGDoFAccessor (const MGDoFHandler *mg_dof_handler) : + mg_dof_handler(const_cast*>(mg_dof_handler)) {}; /** * Reset the DoF handler pointer. */ void set_mg_dof_handler (MGDoFHandler *dh) { - Assert (dh != 0, ExcInvalidObject()); + Assert (dh != 0, DoFAccessor::ExcInvalidObject()); mg_dof_handler = dh; }; @@ -133,9 +140,15 @@ class MGDoFAccessor { * @author Wolfgang Bangerth, 1998 */ template -class MGDoFLineAccessor : public MGDoFAccessor, - public DoFLineAccessor { +class MGDoFLineAccessor : public MGDoFAccessor, public DoFLineAccessor { public: + /** + * Declare the data type that this accessor + * class expects to get passed from the + * iterator classes. + */ + typedef MGDoFHandler AccessorData; + /** * Default constructor, unused thus * not implemented. @@ -150,7 +163,7 @@ class MGDoFLineAccessor : public MGDoFAccessor, MGDoFLineAccessor (Triangulation *tria, const int level, const int index, - const void *local_data); + const AccessorData *local_data); /** * Return the index of the #i#th degree @@ -223,9 +236,15 @@ class MGDoFLineAccessor : public MGDoFAccessor, * @see DoFLineAccessor */ template -class MGDoFQuadAccessor : public MGDoFAccessor, - public DoFQuadAccessor { +class MGDoFQuadAccessor : public MGDoFAccessor, public DoFQuadAccessor { public: + /** + * Declare the data type that this accessor + * class expects to get passed from the + * iterator classes. + */ + typedef MGDoFHandler AccessorData; + /** * Default constructor, unused thus * not implemented. @@ -240,7 +259,7 @@ class MGDoFQuadAccessor : public MGDoFAccessor, MGDoFQuadAccessor (Triangulation *tria, const int level, const int index, - const void *local_data); + const AccessorData *local_data); /** * Return the index of the #i#th degree @@ -307,7 +326,7 @@ class MGDoFQuadAccessor : public MGDoFAccessor, * Implement the copy operator needed * for the iterator classes. */ - void copy_from (const DoFQuadAccessor &a); + void copy_from (const MGDoFQuadAccessor &a); }; @@ -329,8 +348,7 @@ class MGDoFQuadAccessor : public MGDoFAccessor, * the compiler) happens to use this class. */ template -class MGDoFSubstructAccessor : public MGDoFAccessor, - public TriaAccessor { +class MGDoFSubstructAccessor : public MGDoFLineAccessor<1,CellAccessor<1> > { public: DoFSubstructAccessor () { Assert (false, ExcInternalError()); @@ -351,8 +369,8 @@ class MGDoFSubstructAccessor : public MGDoFAccessor, * with template parameters. This class and #DoFSubstructAccessor<2># * wrap the following names: * \begin{verbatim} - * DoFSubstructAccessor<1> := DoFLineAccessor<1,CellAccessor<1> >; - * DoFSubstructAccessor<2> := DoFQuadAccessor<2,CellAccessor<2> >; + * MGDoFSubstructAccessor<1> := MGDoFLineAccessor<1,CellAccessor<1> >; + * MGDoFSubstructAccessor<2> := MGDoFQuadAccessor<2,CellAccessor<2> >; * \end{verbatim} * We do this rather complex (and needless, provided C++ the needed constructs!) * class hierarchy manipulation, since this way we can declare and implement @@ -366,13 +384,20 @@ class MGDoFSubstructAccessor : public MGDoFAccessor, template <> class MGDoFSubstructAccessor<1> : public MGDoFLineAccessor<1,CellAccessor<1> > { public: + /** + * Declare the data type that this accessor + * class expects to get passed from the + * iterator classes. + */ + typedef typename MGDoFLineAccessor<1,CellAccessor<1> >::AccessorData AccessorData; + /** * Constructor */ MGDoFSubstructAccessor (Triangulation<1> *tria, const int level, const int index, - const void *local_data) : + const AccessorData *local_data) : MGDoFLineAccessor<1,CellAccessor<1> > (tria,level,index,local_data) {}; // do this here, since this way the // CellAccessor has the possibility to know @@ -393,13 +418,20 @@ class MGDoFSubstructAccessor<1> : public MGDoFLineAccessor<1,CellAccessor<1> > template <> class MGDoFSubstructAccessor<2> : public MGDoFQuadAccessor<2,CellAccessor<2> > { public: + /** + * Declare the data type that this accessor + * class expects to get passed from the + * iterator classes. + */ + typedef MGDoFQuadAccessor<2,CellAccessor<2> >::AccessorData AccessorData; + /** * Constructor */ MGDoFSubstructAccessor (Triangulation<2> *tria, const int level, const int index, - const void *local_data) : + const AccessorData *local_data) : MGDoFQuadAccessor<2,CellAccessor<2> > (tria,level,index,local_data) {}; // do this here, since this way the // CellAccessor has the possibility to know @@ -431,13 +463,20 @@ class MGDoFSubstructAccessor<2> : public MGDoFQuadAccessor<2,CellAccessor<2> > { template class MGDoFCellAccessor : public MGDoFSubstructAccessor { public: + /** + * Declare the data type that this accessor + * class expects to get passed from the + * iterator classes. + */ + typedef typename MGDoFSubstructAccessor::AccessorData AccessorData; + /** * Constructor */ MGDoFCellAccessor (Triangulation *tria, const int level, const int index, - const void *local_data) : + const AccessorData *local_data) : MGDoFSubstructAccessor (tria,level,index,local_data) {}; /** diff --git a/deal.II/deal.II/include/dofs/mg_dof_handler.h b/deal.II/deal.II/include/dofs/mg_dof_handler.h index 21ad59db20..0e52401f37 100644 --- a/deal.II/deal.II/include/dofs/mg_dof_handler.h +++ b/deal.II/deal.II/include/dofs/mg_dof_handler.h @@ -11,9 +11,78 @@ // forward declarations template class MGDoFLineAccessor; -template class MGDoFLineAccessor; -template class MGDoFQuadAccessor; template class MGDoFQuadAccessor; +template class MGDoFCellAccessor; + + + + +/** + * Define some types which differ between the dimensions. This class + * is analogous to the \Ref{TriaDimensionInfo} class hierarchy. + * + * @see MGDoFDimensionInfo<1> + * @see MGDoFDimensionInfo<2> + */ +template +class MGDoFDimensionInfo; + + + + + +/** + * Define some types for the DoF handling in one dimension. + * + * The types have the same meaning as those declared in \Ref{TriaDimensionInfo<2>}. + */ +class MGDoFDimensionInfo<1> { + public: + typedef TriaRawIterator<1,MGDoFCellAccessor<1> > raw_line_iterator; + typedef TriaIterator<1,MGDoFCellAccessor<1> > line_iterator; + typedef TriaActiveIterator<1,MGDoFCellAccessor<1> > active_line_iterator; + + typedef void * raw_quad_iterator; + typedef void * quad_iterator; + typedef void * active_quad_iterator; + + typedef raw_line_iterator raw_cell_iterator; + typedef line_iterator cell_iterator; + typedef active_line_iterator active_cell_iterator; + + typedef void * raw_face_iterator; + typedef void * face_iterator; + typedef void * active_face_iterator; +}; + + + + + + +/** + * Define some types for the DoF handling in two dimensions. + * + * The types have the same meaning as those declared in \Ref{TriaDimensionInfo<2>}. + */ +class MGDoFDimensionInfo<2> { + public: + typedef TriaRawIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > > raw_line_iterator; + typedef TriaIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > > line_iterator; + typedef TriaActiveIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > > active_line_iterator; + + typedef TriaRawIterator<2,MGDoFCellAccessor<2> > raw_quad_iterator; + typedef TriaIterator<2,MGDoFCellAccessor<2> > quad_iterator; + typedef TriaActiveIterator<2,MGDoFCellAccessor<2> > active_quad_iterator; + + typedef raw_quad_iterator raw_cell_iterator; + typedef quad_iterator cell_iterator; + typedef active_quad_iterator active_cell_iterator; + + typedef raw_line_iterator raw_face_iterator; + typedef line_iterator face_iterator; + typedef active_line_iterator active_face_iterator; +}; @@ -22,7 +91,26 @@ template class MGDoFQuadAccessor; template class MGDoFHandler : public DoFHandler { public: - /** + // insert these definitions for gcc2.8, + // since it can't inherit typedefs (I + // believe it should, but it can't) + typedef typename MGDoFDimensionInfo::raw_line_iterator raw_line_iterator; + typedef typename MGDoFDimensionInfo::line_iterator line_iterator; + typedef typename MGDoFDimensionInfo::active_line_iterator active_line_iterator; + + typedef typename MGDoFDimensionInfo::raw_quad_iterator raw_quad_iterator; + typedef typename MGDoFDimensionInfo::quad_iterator quad_iterator; + typedef typename MGDoFDimensionInfo::active_quad_iterator active_quad_iterator; + + typedef typename MGDoFDimensionInfo::raw_cell_iterator raw_cell_iterator; + typedef typename MGDoFDimensionInfo::cell_iterator cell_iterator; + typedef typename MGDoFDimensionInfo::active_cell_iterator active_cell_iterator; + + typedef typename MGDoFDimensionInfo::raw_face_iterator raw_face_iterator; + typedef typename MGDoFDimensionInfo::face_iterator face_iterator; + typedef typename MGDoFDimensionInfo::active_face_iterator active_face_iterator; + + /** * Constructor. Take #tria# as the * triangulation to work on. */ @@ -69,6 +157,473 @@ class MGDoFHandler : public DoFHandler { const vector &starting_points = vector()); + /*--------------------------------------*/ + + /** + * @name Cell iterator functions + */ + /*@{*/ + /** + * Return iterator to the first cell, used + * or not, on level #level#. If a level + * has no cells, a past-the-end iterator + * is returned. + * + * This function calls #begin_raw_line# + * in 1D and #begin_raw_quad# in 2D. + */ + raw_cell_iterator begin_raw (const unsigned int level = 0) const; + + /** + * Return iterator to the first used cell + * on level #level#. + * + * This function calls #begin_line# + * in 1D and #begin_quad# in 2D. + */ + cell_iterator begin (const unsigned int level = 0) const; + + /** + * Return iterator to the first active + * cell on level #level#. + * + * This function calls #begin_active_line# + * in 1D and #begin_active_quad# in 2D. + */ + active_cell_iterator begin_active(const unsigned int level = 0) const; + + /** + * Return iterator past the end; this + * iterator serves for comparisons of + * iterators with past-the-end or + * before-the-beginning states. + * + * This function calls #end_line# + * in 1D and #end_quad# in 2D. + */ + raw_cell_iterator end () const; + + /** + * Return an iterator which is the first + * iterator not on level. If #level# is + * the last level, then this returns + * #end()#. + */ + cell_iterator end (const unsigned int level) const; + + /** + * Return a raw iterator which is the first + * iterator not on level. If #level# is + * the last level, then this returns + * #end()#. + */ + raw_cell_iterator end_raw (const unsigned int level) const; + + /** + * Return an active iterator which is the + * first iterator not on level. If #level# + * is the last level, then this returns + * #end()#. + */ + active_cell_iterator end_active (const unsigned int level) const; + + /** + * Return an iterator pointing to the + * last cell, used or not. + * + * This function calls #last_raw_line# + * in 1D and #last_raw_quad# in 2D. + */ + raw_cell_iterator last_raw () const; + + /** + * Return an iterator pointing to the last + * cell of the level #level#, used or not. + * + * This function calls #last_raw_line# + * in 1D and #last_raw_quad# in 2D. + */ + raw_cell_iterator last_raw (const unsigned int level) const; + + /** + * Return an iterator pointing to the last + * used cell. + * + * This function calls #last_line# + * in 1D and #last_quad# in 2D. + */ + cell_iterator last () const; + + /** + * Return an iterator pointing to the last + * used cell on level #level#. + * + * This function calls #last_line# + * in 1D and #last_quad# in 2D. + */ + cell_iterator last (const unsigned int level) const; + + /** + * Return an iterator pointing to the last + * active cell. + * + * This function calls #last_active_line# + * in 1D and #last_active_quad# in 2D. + */ + active_cell_iterator last_active () const; + + /** + * Return an iterator pointing to the last + * active cell on level #level#. + * + * This function calls #last_active_line# + * in 1D and #last_active_quad# in 2D. + */ + active_cell_iterator last_active (const unsigned int level) const; + //@} + + /*---------------------------------------*/ + + /** + * @name Face iterator functions + */ + /*@{*/ + /** + * Return iterator to the first face, used + * or not, on level #level#. If a level + * has no faces, a past-the-end iterator + * is returned. + * + * This function calls #begin_raw_line# + * in 2D and #begin_raw_quad# in 3D. + */ + raw_face_iterator begin_raw_face (const unsigned int level = 0) const; + + /** + * Return iterator to the first used face + * on level #level#. + * + * This function calls #begin_line# + * in 2D and #begin_quad# in 3D. + */ + face_iterator begin_face (const unsigned int level = 0) const; + + /** + * Return iterator to the first active + * face on level #level#. + * + * This function calls #begin_active_line# + * in 2D and #begin_active_quad# in 3D. + */ + active_face_iterator begin_active_face(const unsigned int level = 0) const; + + /** + * Return iterator past the end; this + * iterator serves for comparisons of + * iterators with past-the-end or + * before-the-beginning states. + * + * This function calls #end_line# + * in 2D and #end_quad# in 3D. + */ + raw_face_iterator end_face () const; + + /** + * Return an iterator which is the first + * iterator not on level. If #level# is + * the last level, then this returns + * #end()#. + */ + face_iterator end_face (const unsigned int level) const; + + /** + * Return a raw iterator which is the first + * iterator not on level. If #level# is + * the last level, then this returns + * #end()#. + */ + raw_face_iterator end_raw_face (const unsigned int level) const; + + /** + * Return an active iterator which is the + * first iterator not on level. If #level# + * is the last level, then this returns + * #end()#. + */ + active_face_iterator end_active_face (const unsigned int level) const; + + /** + * Return an iterator pointing to the + * last face, used or not. + * + * This function calls #last_raw_line# + * in 2D and #last_raw_quad# in 3D. + */ + raw_face_iterator last_raw_face () const; + + /** + * Return an iterator pointing to the last + * face of the level #level#, used or not. + * + * This function calls #last_raw_line# + * in 2D and #last_raw_quad# in 3D. + */ + raw_face_iterator last_raw_face (const unsigned int level) const; + + /** + * Return an iterator pointing to the last + * used face. + * + * This function calls #last_line# + * in 2D and #last_quad# in 3D. + */ + face_iterator last_face () const; + + /** + * Return an iterator pointing to the last + * used face on level #level#. + * + * This function calls #last_line# + * in 2D and #last_quad# in 3D. + */ + face_iterator last_face (const unsigned int level) const; + + /** + * Return an iterator pointing to the last + * active face. + * + * This function calls #last_active_line# + * in 2D and #last_active_quad# in 3D. + */ + active_face_iterator last_active_face () const; + + /** + * Return an iterator pointing to the last + * active face on level #level#. + * + * This function calls #last_active_line# + * in 2D and #last_active_quad# in 3D. + */ + active_face_iterator last_active_face (const unsigned int level) const; + //@} + + + /*---------------------------------------*/ + + /** + * @name Line iterator functions + */ + /*@{*/ + /** + * Return iterator to the first line, used + * or not, on level #level#. If a level + * has no lines, a past-the-end iterator + * is returned. + */ + raw_line_iterator + begin_raw_line (const unsigned int level = 0) const; + + /** + * Return iterator to the first used line + * on level #level#. + */ + line_iterator + begin_line (const unsigned int level = 0) const; + + /** + * Return iterator to the first active + * line on level #level#. + */ + active_line_iterator + begin_active_line(const unsigned int level = 0) const; + + /** + * Return iterator past the end; this + * iterator serves for comparisons of + * iterators with past-the-end or + * before-the-beginning states. + */ + raw_line_iterator + end_line () const; + + /** + * Return an iterator which is the first + * iterator not on level. If #level# is + * the last level, then this returns + * #end()#. + */ + line_iterator end_line (const unsigned int level) const; + + /** + * Return a raw iterator which is the first + * iterator not on level. If #level# is + * the last level, then this returns + * #end()#. + */ + raw_line_iterator end_raw_line (const unsigned int level) const; + + /** + * Return an active iterator which is the + * first iterator not on level. If #level# + * is the last level, then this returns + * #end()#. + */ + active_line_iterator end_active_line (const unsigned int level) const; + + + /** + * Return an iterator pointing to the + * last line, used or not. + */ + raw_line_iterator + last_raw_line () const; + + /** + * Return an iterator pointing to the last + * line of the level #level#, used or not. + + */ + raw_line_iterator + last_raw_line (const unsigned int level) const; + + /** + * Return an iterator pointing to the last + * used line. + */ + line_iterator + last_line () const; + + /** + * Return an iterator pointing to the last + * used line on level #level#. + */ + line_iterator + last_line (const unsigned int level) const; + + /** + * Return an iterator pointing to the last + * active line. + */ + active_line_iterator + last_active_line () const; + + /** + * Return an iterator pointing to the last + * active line on level #level#. + */ + active_line_iterator + last_active_line (const unsigned int level) const; + /*@}*/ + + /*---------------------------------------*/ + + /** + * @name Quad iterator functions*/ + /*@{ + */ + /** + * Return iterator to the first quad, used + * or not, on level #level#. If a level + * has no quads, a past-the-end iterator + * is returned. + */ + raw_quad_iterator + begin_raw_quad (const unsigned int level = 0) const; + + /** + * Return iterator to the first used quad + * on level #level#. + */ + quad_iterator + begin_quad (const unsigned int level = 0) const; + + /** + * Return iterator to the first active + * quad on level #level#. + */ + active_quad_iterator + begin_active_quad(const unsigned int level = 0) const; + + /** + * Return iterator past the end; this + * iterator serves for comparisons of + * iterators with past-the-end or + * before-the-beginning states. + */ + raw_quad_iterator + end_quad () const; + + /** + * Return an iterator which is the first + * iterator not on level. If #level# is + * the last level, then this returns + * #end()#. + */ + quad_iterator end_quad (const unsigned int level) const; + + /** + * Return a raw iterator which is the first + * iterator not on level. If #level# is + * the last level, then this returns + * #end()#. + */ + raw_quad_iterator end_raw_quad (const unsigned int level) const; + + /** + * Return an active iterator which is the + * first iterator not on level. If #level# + * is the last level, then this returns + * #end()#. + */ + active_quad_iterator end_active_quad (const unsigned int level) const; + + + /** + * Return an iterator pointing to the + * last quad, used or not. + */ + raw_quad_iterator + last_raw_quad () const; + + /** + * Return an iterator pointing to the last + * quad of the level #level#, used or not. + + */ + raw_quad_iterator + last_raw_quad (const unsigned int level) const; + + /** + * Return an iterator pointing to the last + * used quad. + */ + quad_iterator + last_quad () const; + + /** + * Return an iterator pointing to the last + * used quad on level #level#. + */ + quad_iterator + last_quad (const unsigned int level) const; + + /** + * Return an iterator pointing to the last + * active quad. + */ + active_quad_iterator + last_active_quad () const; + + /** + * Return an iterator pointing to the last + * active quad on level #level#. + */ + active_quad_iterator + last_active_quad (const unsigned int level) const; + /*@}*/ + + /*---------------------------------------*/ + private: /** We need for each vertex a list of the @@ -93,12 +648,23 @@ class MGDoFHandler : public DoFHandler { public: /** - * Constructor. Allocates memory and - * sets all indices to #-1#. + * Constructor. This one is empty + * because it is difficult to make it + * efficient to use vector<>'s and + * still construct the object using + * the constructor. Use the #init# + * function to really allocate + * memory. */ - MGVertexDoFs (const unsigned int coarsest_level, - const unsigned int n_levels, - const unsigned int dofs_per_vertex); + MGVertexDoFs (); + + /** + * Allocate memory and + * set all indices to #-1#. + */ + void init (const unsigned int coarsest_level, + const unsigned int n_levels, + const unsigned int dofs_per_vertex); /** * Destructor @@ -148,6 +714,7 @@ class MGDoFHandler : public DoFHandler { DeclException0 (ExcNoMemory); /** * Exception. + */ DeclException0 (ExcInvalidIndex); /** @@ -157,6 +724,10 @@ class MGDoFHandler : public DoFHandler { int, << "The given level number " << arg1 << " is below the " << "coarsest level this vertex lives on."); + /** + * Exception. + */ + DeclException0 (ExcInternalError); private: /** @@ -165,12 +736,12 @@ class MGDoFHandler : public DoFHandler { * as an offset when accessing the * dofs of a specific level. */ - const unsigned int coarsest_level; + unsigned int coarsest_level; /** * Array holding the indices. */ - int *const indices; + int *indices; }; @@ -180,6 +751,23 @@ class MGDoFHandler : public DoFHandler { */ void reserve_space (); + /** + * Distribute dofs on the given cell, + * with new dofs starting with index + * #next_free_dof#. Return the next + * unused index number. The finite + * element used is the one given to + * #distribute_dofs#, which is copied + * to #selected_fe#. + * + * This function is excluded from the + * #distribute_dofs# function since + * it can not be implemented dimension + * independent. + */ + unsigned int distribute_dofs_on_cell (active_cell_iterator &cell, + unsigned int next_free_dof); + /** * Space to store the DoF numbers for the * different levels. Unlike the #levels# diff --git a/deal.II/deal.II/include/multigrid/mg_dof_accessor.h b/deal.II/deal.II/include/multigrid/mg_dof_accessor.h index 03274755dc..75ba6a5a74 100644 --- a/deal.II/deal.II/include/multigrid/mg_dof_accessor.h +++ b/deal.II/deal.II/include/multigrid/mg_dof_accessor.h @@ -51,15 +51,22 @@ class MGDoFAccessor { /** * This should be the default constructor. + * We cast away the #const#ness of the + * pointer which clearly is EVIL but + * we can't help without making all + * functions which could somehow use + * iterators (directly or indirectly) make + * non-const, even if they preserve + * constness. */ - MGDoFAccessor (MGDoFHandler *mg_dof_handler) : - mg_dof_handler(mg_dof_handler) {}; + MGDoFAccessor (const MGDoFHandler *mg_dof_handler) : + mg_dof_handler(const_cast*>(mg_dof_handler)) {}; /** * Reset the DoF handler pointer. */ void set_mg_dof_handler (MGDoFHandler *dh) { - Assert (dh != 0, ExcInvalidObject()); + Assert (dh != 0, DoFAccessor::ExcInvalidObject()); mg_dof_handler = dh; }; @@ -133,9 +140,15 @@ class MGDoFAccessor { * @author Wolfgang Bangerth, 1998 */ template -class MGDoFLineAccessor : public MGDoFAccessor, - public DoFLineAccessor { +class MGDoFLineAccessor : public MGDoFAccessor, public DoFLineAccessor { public: + /** + * Declare the data type that this accessor + * class expects to get passed from the + * iterator classes. + */ + typedef MGDoFHandler AccessorData; + /** * Default constructor, unused thus * not implemented. @@ -150,7 +163,7 @@ class MGDoFLineAccessor : public MGDoFAccessor, MGDoFLineAccessor (Triangulation *tria, const int level, const int index, - const void *local_data); + const AccessorData *local_data); /** * Return the index of the #i#th degree @@ -223,9 +236,15 @@ class MGDoFLineAccessor : public MGDoFAccessor, * @see DoFLineAccessor */ template -class MGDoFQuadAccessor : public MGDoFAccessor, - public DoFQuadAccessor { +class MGDoFQuadAccessor : public MGDoFAccessor, public DoFQuadAccessor { public: + /** + * Declare the data type that this accessor + * class expects to get passed from the + * iterator classes. + */ + typedef MGDoFHandler AccessorData; + /** * Default constructor, unused thus * not implemented. @@ -240,7 +259,7 @@ class MGDoFQuadAccessor : public MGDoFAccessor, MGDoFQuadAccessor (Triangulation *tria, const int level, const int index, - const void *local_data); + const AccessorData *local_data); /** * Return the index of the #i#th degree @@ -307,7 +326,7 @@ class MGDoFQuadAccessor : public MGDoFAccessor, * Implement the copy operator needed * for the iterator classes. */ - void copy_from (const DoFQuadAccessor &a); + void copy_from (const MGDoFQuadAccessor &a); }; @@ -329,8 +348,7 @@ class MGDoFQuadAccessor : public MGDoFAccessor, * the compiler) happens to use this class. */ template -class MGDoFSubstructAccessor : public MGDoFAccessor, - public TriaAccessor { +class MGDoFSubstructAccessor : public MGDoFLineAccessor<1,CellAccessor<1> > { public: DoFSubstructAccessor () { Assert (false, ExcInternalError()); @@ -351,8 +369,8 @@ class MGDoFSubstructAccessor : public MGDoFAccessor, * with template parameters. This class and #DoFSubstructAccessor<2># * wrap the following names: * \begin{verbatim} - * DoFSubstructAccessor<1> := DoFLineAccessor<1,CellAccessor<1> >; - * DoFSubstructAccessor<2> := DoFQuadAccessor<2,CellAccessor<2> >; + * MGDoFSubstructAccessor<1> := MGDoFLineAccessor<1,CellAccessor<1> >; + * MGDoFSubstructAccessor<2> := MGDoFQuadAccessor<2,CellAccessor<2> >; * \end{verbatim} * We do this rather complex (and needless, provided C++ the needed constructs!) * class hierarchy manipulation, since this way we can declare and implement @@ -366,13 +384,20 @@ class MGDoFSubstructAccessor : public MGDoFAccessor, template <> class MGDoFSubstructAccessor<1> : public MGDoFLineAccessor<1,CellAccessor<1> > { public: + /** + * Declare the data type that this accessor + * class expects to get passed from the + * iterator classes. + */ + typedef typename MGDoFLineAccessor<1,CellAccessor<1> >::AccessorData AccessorData; + /** * Constructor */ MGDoFSubstructAccessor (Triangulation<1> *tria, const int level, const int index, - const void *local_data) : + const AccessorData *local_data) : MGDoFLineAccessor<1,CellAccessor<1> > (tria,level,index,local_data) {}; // do this here, since this way the // CellAccessor has the possibility to know @@ -393,13 +418,20 @@ class MGDoFSubstructAccessor<1> : public MGDoFLineAccessor<1,CellAccessor<1> > template <> class MGDoFSubstructAccessor<2> : public MGDoFQuadAccessor<2,CellAccessor<2> > { public: + /** + * Declare the data type that this accessor + * class expects to get passed from the + * iterator classes. + */ + typedef MGDoFQuadAccessor<2,CellAccessor<2> >::AccessorData AccessorData; + /** * Constructor */ MGDoFSubstructAccessor (Triangulation<2> *tria, const int level, const int index, - const void *local_data) : + const AccessorData *local_data) : MGDoFQuadAccessor<2,CellAccessor<2> > (tria,level,index,local_data) {}; // do this here, since this way the // CellAccessor has the possibility to know @@ -431,13 +463,20 @@ class MGDoFSubstructAccessor<2> : public MGDoFQuadAccessor<2,CellAccessor<2> > { template class MGDoFCellAccessor : public MGDoFSubstructAccessor { public: + /** + * Declare the data type that this accessor + * class expects to get passed from the + * iterator classes. + */ + typedef typename MGDoFSubstructAccessor::AccessorData AccessorData; + /** * Constructor */ MGDoFCellAccessor (Triangulation *tria, const int level, const int index, - const void *local_data) : + const AccessorData *local_data) : MGDoFSubstructAccessor (tria,level,index,local_data) {}; /** diff --git a/deal.II/deal.II/include/multigrid/mg_dof_handler.h b/deal.II/deal.II/include/multigrid/mg_dof_handler.h index 21ad59db20..0e52401f37 100644 --- a/deal.II/deal.II/include/multigrid/mg_dof_handler.h +++ b/deal.II/deal.II/include/multigrid/mg_dof_handler.h @@ -11,9 +11,78 @@ // forward declarations template class MGDoFLineAccessor; -template class MGDoFLineAccessor; -template class MGDoFQuadAccessor; template class MGDoFQuadAccessor; +template class MGDoFCellAccessor; + + + + +/** + * Define some types which differ between the dimensions. This class + * is analogous to the \Ref{TriaDimensionInfo} class hierarchy. + * + * @see MGDoFDimensionInfo<1> + * @see MGDoFDimensionInfo<2> + */ +template +class MGDoFDimensionInfo; + + + + + +/** + * Define some types for the DoF handling in one dimension. + * + * The types have the same meaning as those declared in \Ref{TriaDimensionInfo<2>}. + */ +class MGDoFDimensionInfo<1> { + public: + typedef TriaRawIterator<1,MGDoFCellAccessor<1> > raw_line_iterator; + typedef TriaIterator<1,MGDoFCellAccessor<1> > line_iterator; + typedef TriaActiveIterator<1,MGDoFCellAccessor<1> > active_line_iterator; + + typedef void * raw_quad_iterator; + typedef void * quad_iterator; + typedef void * active_quad_iterator; + + typedef raw_line_iterator raw_cell_iterator; + typedef line_iterator cell_iterator; + typedef active_line_iterator active_cell_iterator; + + typedef void * raw_face_iterator; + typedef void * face_iterator; + typedef void * active_face_iterator; +}; + + + + + + +/** + * Define some types for the DoF handling in two dimensions. + * + * The types have the same meaning as those declared in \Ref{TriaDimensionInfo<2>}. + */ +class MGDoFDimensionInfo<2> { + public: + typedef TriaRawIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > > raw_line_iterator; + typedef TriaIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > > line_iterator; + typedef TriaActiveIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > > active_line_iterator; + + typedef TriaRawIterator<2,MGDoFCellAccessor<2> > raw_quad_iterator; + typedef TriaIterator<2,MGDoFCellAccessor<2> > quad_iterator; + typedef TriaActiveIterator<2,MGDoFCellAccessor<2> > active_quad_iterator; + + typedef raw_quad_iterator raw_cell_iterator; + typedef quad_iterator cell_iterator; + typedef active_quad_iterator active_cell_iterator; + + typedef raw_line_iterator raw_face_iterator; + typedef line_iterator face_iterator; + typedef active_line_iterator active_face_iterator; +}; @@ -22,7 +91,26 @@ template class MGDoFQuadAccessor; template class MGDoFHandler : public DoFHandler { public: - /** + // insert these definitions for gcc2.8, + // since it can't inherit typedefs (I + // believe it should, but it can't) + typedef typename MGDoFDimensionInfo::raw_line_iterator raw_line_iterator; + typedef typename MGDoFDimensionInfo::line_iterator line_iterator; + typedef typename MGDoFDimensionInfo::active_line_iterator active_line_iterator; + + typedef typename MGDoFDimensionInfo::raw_quad_iterator raw_quad_iterator; + typedef typename MGDoFDimensionInfo::quad_iterator quad_iterator; + typedef typename MGDoFDimensionInfo::active_quad_iterator active_quad_iterator; + + typedef typename MGDoFDimensionInfo::raw_cell_iterator raw_cell_iterator; + typedef typename MGDoFDimensionInfo::cell_iterator cell_iterator; + typedef typename MGDoFDimensionInfo::active_cell_iterator active_cell_iterator; + + typedef typename MGDoFDimensionInfo::raw_face_iterator raw_face_iterator; + typedef typename MGDoFDimensionInfo::face_iterator face_iterator; + typedef typename MGDoFDimensionInfo::active_face_iterator active_face_iterator; + + /** * Constructor. Take #tria# as the * triangulation to work on. */ @@ -69,6 +157,473 @@ class MGDoFHandler : public DoFHandler { const vector &starting_points = vector()); + /*--------------------------------------*/ + + /** + * @name Cell iterator functions + */ + /*@{*/ + /** + * Return iterator to the first cell, used + * or not, on level #level#. If a level + * has no cells, a past-the-end iterator + * is returned. + * + * This function calls #begin_raw_line# + * in 1D and #begin_raw_quad# in 2D. + */ + raw_cell_iterator begin_raw (const unsigned int level = 0) const; + + /** + * Return iterator to the first used cell + * on level #level#. + * + * This function calls #begin_line# + * in 1D and #begin_quad# in 2D. + */ + cell_iterator begin (const unsigned int level = 0) const; + + /** + * Return iterator to the first active + * cell on level #level#. + * + * This function calls #begin_active_line# + * in 1D and #begin_active_quad# in 2D. + */ + active_cell_iterator begin_active(const unsigned int level = 0) const; + + /** + * Return iterator past the end; this + * iterator serves for comparisons of + * iterators with past-the-end or + * before-the-beginning states. + * + * This function calls #end_line# + * in 1D and #end_quad# in 2D. + */ + raw_cell_iterator end () const; + + /** + * Return an iterator which is the first + * iterator not on level. If #level# is + * the last level, then this returns + * #end()#. + */ + cell_iterator end (const unsigned int level) const; + + /** + * Return a raw iterator which is the first + * iterator not on level. If #level# is + * the last level, then this returns + * #end()#. + */ + raw_cell_iterator end_raw (const unsigned int level) const; + + /** + * Return an active iterator which is the + * first iterator not on level. If #level# + * is the last level, then this returns + * #end()#. + */ + active_cell_iterator end_active (const unsigned int level) const; + + /** + * Return an iterator pointing to the + * last cell, used or not. + * + * This function calls #last_raw_line# + * in 1D and #last_raw_quad# in 2D. + */ + raw_cell_iterator last_raw () const; + + /** + * Return an iterator pointing to the last + * cell of the level #level#, used or not. + * + * This function calls #last_raw_line# + * in 1D and #last_raw_quad# in 2D. + */ + raw_cell_iterator last_raw (const unsigned int level) const; + + /** + * Return an iterator pointing to the last + * used cell. + * + * This function calls #last_line# + * in 1D and #last_quad# in 2D. + */ + cell_iterator last () const; + + /** + * Return an iterator pointing to the last + * used cell on level #level#. + * + * This function calls #last_line# + * in 1D and #last_quad# in 2D. + */ + cell_iterator last (const unsigned int level) const; + + /** + * Return an iterator pointing to the last + * active cell. + * + * This function calls #last_active_line# + * in 1D and #last_active_quad# in 2D. + */ + active_cell_iterator last_active () const; + + /** + * Return an iterator pointing to the last + * active cell on level #level#. + * + * This function calls #last_active_line# + * in 1D and #last_active_quad# in 2D. + */ + active_cell_iterator last_active (const unsigned int level) const; + //@} + + /*---------------------------------------*/ + + /** + * @name Face iterator functions + */ + /*@{*/ + /** + * Return iterator to the first face, used + * or not, on level #level#. If a level + * has no faces, a past-the-end iterator + * is returned. + * + * This function calls #begin_raw_line# + * in 2D and #begin_raw_quad# in 3D. + */ + raw_face_iterator begin_raw_face (const unsigned int level = 0) const; + + /** + * Return iterator to the first used face + * on level #level#. + * + * This function calls #begin_line# + * in 2D and #begin_quad# in 3D. + */ + face_iterator begin_face (const unsigned int level = 0) const; + + /** + * Return iterator to the first active + * face on level #level#. + * + * This function calls #begin_active_line# + * in 2D and #begin_active_quad# in 3D. + */ + active_face_iterator begin_active_face(const unsigned int level = 0) const; + + /** + * Return iterator past the end; this + * iterator serves for comparisons of + * iterators with past-the-end or + * before-the-beginning states. + * + * This function calls #end_line# + * in 2D and #end_quad# in 3D. + */ + raw_face_iterator end_face () const; + + /** + * Return an iterator which is the first + * iterator not on level. If #level# is + * the last level, then this returns + * #end()#. + */ + face_iterator end_face (const unsigned int level) const; + + /** + * Return a raw iterator which is the first + * iterator not on level. If #level# is + * the last level, then this returns + * #end()#. + */ + raw_face_iterator end_raw_face (const unsigned int level) const; + + /** + * Return an active iterator which is the + * first iterator not on level. If #level# + * is the last level, then this returns + * #end()#. + */ + active_face_iterator end_active_face (const unsigned int level) const; + + /** + * Return an iterator pointing to the + * last face, used or not. + * + * This function calls #last_raw_line# + * in 2D and #last_raw_quad# in 3D. + */ + raw_face_iterator last_raw_face () const; + + /** + * Return an iterator pointing to the last + * face of the level #level#, used or not. + * + * This function calls #last_raw_line# + * in 2D and #last_raw_quad# in 3D. + */ + raw_face_iterator last_raw_face (const unsigned int level) const; + + /** + * Return an iterator pointing to the last + * used face. + * + * This function calls #last_line# + * in 2D and #last_quad# in 3D. + */ + face_iterator last_face () const; + + /** + * Return an iterator pointing to the last + * used face on level #level#. + * + * This function calls #last_line# + * in 2D and #last_quad# in 3D. + */ + face_iterator last_face (const unsigned int level) const; + + /** + * Return an iterator pointing to the last + * active face. + * + * This function calls #last_active_line# + * in 2D and #last_active_quad# in 3D. + */ + active_face_iterator last_active_face () const; + + /** + * Return an iterator pointing to the last + * active face on level #level#. + * + * This function calls #last_active_line# + * in 2D and #last_active_quad# in 3D. + */ + active_face_iterator last_active_face (const unsigned int level) const; + //@} + + + /*---------------------------------------*/ + + /** + * @name Line iterator functions + */ + /*@{*/ + /** + * Return iterator to the first line, used + * or not, on level #level#. If a level + * has no lines, a past-the-end iterator + * is returned. + */ + raw_line_iterator + begin_raw_line (const unsigned int level = 0) const; + + /** + * Return iterator to the first used line + * on level #level#. + */ + line_iterator + begin_line (const unsigned int level = 0) const; + + /** + * Return iterator to the first active + * line on level #level#. + */ + active_line_iterator + begin_active_line(const unsigned int level = 0) const; + + /** + * Return iterator past the end; this + * iterator serves for comparisons of + * iterators with past-the-end or + * before-the-beginning states. + */ + raw_line_iterator + end_line () const; + + /** + * Return an iterator which is the first + * iterator not on level. If #level# is + * the last level, then this returns + * #end()#. + */ + line_iterator end_line (const unsigned int level) const; + + /** + * Return a raw iterator which is the first + * iterator not on level. If #level# is + * the last level, then this returns + * #end()#. + */ + raw_line_iterator end_raw_line (const unsigned int level) const; + + /** + * Return an active iterator which is the + * first iterator not on level. If #level# + * is the last level, then this returns + * #end()#. + */ + active_line_iterator end_active_line (const unsigned int level) const; + + + /** + * Return an iterator pointing to the + * last line, used or not. + */ + raw_line_iterator + last_raw_line () const; + + /** + * Return an iterator pointing to the last + * line of the level #level#, used or not. + + */ + raw_line_iterator + last_raw_line (const unsigned int level) const; + + /** + * Return an iterator pointing to the last + * used line. + */ + line_iterator + last_line () const; + + /** + * Return an iterator pointing to the last + * used line on level #level#. + */ + line_iterator + last_line (const unsigned int level) const; + + /** + * Return an iterator pointing to the last + * active line. + */ + active_line_iterator + last_active_line () const; + + /** + * Return an iterator pointing to the last + * active line on level #level#. + */ + active_line_iterator + last_active_line (const unsigned int level) const; + /*@}*/ + + /*---------------------------------------*/ + + /** + * @name Quad iterator functions*/ + /*@{ + */ + /** + * Return iterator to the first quad, used + * or not, on level #level#. If a level + * has no quads, a past-the-end iterator + * is returned. + */ + raw_quad_iterator + begin_raw_quad (const unsigned int level = 0) const; + + /** + * Return iterator to the first used quad + * on level #level#. + */ + quad_iterator + begin_quad (const unsigned int level = 0) const; + + /** + * Return iterator to the first active + * quad on level #level#. + */ + active_quad_iterator + begin_active_quad(const unsigned int level = 0) const; + + /** + * Return iterator past the end; this + * iterator serves for comparisons of + * iterators with past-the-end or + * before-the-beginning states. + */ + raw_quad_iterator + end_quad () const; + + /** + * Return an iterator which is the first + * iterator not on level. If #level# is + * the last level, then this returns + * #end()#. + */ + quad_iterator end_quad (const unsigned int level) const; + + /** + * Return a raw iterator which is the first + * iterator not on level. If #level# is + * the last level, then this returns + * #end()#. + */ + raw_quad_iterator end_raw_quad (const unsigned int level) const; + + /** + * Return an active iterator which is the + * first iterator not on level. If #level# + * is the last level, then this returns + * #end()#. + */ + active_quad_iterator end_active_quad (const unsigned int level) const; + + + /** + * Return an iterator pointing to the + * last quad, used or not. + */ + raw_quad_iterator + last_raw_quad () const; + + /** + * Return an iterator pointing to the last + * quad of the level #level#, used or not. + + */ + raw_quad_iterator + last_raw_quad (const unsigned int level) const; + + /** + * Return an iterator pointing to the last + * used quad. + */ + quad_iterator + last_quad () const; + + /** + * Return an iterator pointing to the last + * used quad on level #level#. + */ + quad_iterator + last_quad (const unsigned int level) const; + + /** + * Return an iterator pointing to the last + * active quad. + */ + active_quad_iterator + last_active_quad () const; + + /** + * Return an iterator pointing to the last + * active quad on level #level#. + */ + active_quad_iterator + last_active_quad (const unsigned int level) const; + /*@}*/ + + /*---------------------------------------*/ + private: /** We need for each vertex a list of the @@ -93,12 +648,23 @@ class MGDoFHandler : public DoFHandler { public: /** - * Constructor. Allocates memory and - * sets all indices to #-1#. + * Constructor. This one is empty + * because it is difficult to make it + * efficient to use vector<>'s and + * still construct the object using + * the constructor. Use the #init# + * function to really allocate + * memory. */ - MGVertexDoFs (const unsigned int coarsest_level, - const unsigned int n_levels, - const unsigned int dofs_per_vertex); + MGVertexDoFs (); + + /** + * Allocate memory and + * set all indices to #-1#. + */ + void init (const unsigned int coarsest_level, + const unsigned int n_levels, + const unsigned int dofs_per_vertex); /** * Destructor @@ -148,6 +714,7 @@ class MGDoFHandler : public DoFHandler { DeclException0 (ExcNoMemory); /** * Exception. + */ DeclException0 (ExcInvalidIndex); /** @@ -157,6 +724,10 @@ class MGDoFHandler : public DoFHandler { int, << "The given level number " << arg1 << " is below the " << "coarsest level this vertex lives on."); + /** + * Exception. + */ + DeclException0 (ExcInternalError); private: /** @@ -165,12 +736,12 @@ class MGDoFHandler : public DoFHandler { * as an offset when accessing the * dofs of a specific level. */ - const unsigned int coarsest_level; + unsigned int coarsest_level; /** * Array holding the indices. */ - int *const indices; + int *indices; }; @@ -180,6 +751,23 @@ class MGDoFHandler : public DoFHandler { */ void reserve_space (); + /** + * Distribute dofs on the given cell, + * with new dofs starting with index + * #next_free_dof#. Return the next + * unused index number. The finite + * element used is the one given to + * #distribute_dofs#, which is copied + * to #selected_fe#. + * + * This function is excluded from the + * #distribute_dofs# function since + * it can not be implemented dimension + * independent. + */ + unsigned int distribute_dofs_on_cell (active_cell_iterator &cell, + unsigned int next_free_dof); + /** * Space to store the DoF numbers for the * different levels. Unlike the #levels# diff --git a/deal.II/deal.II/source/dofs/mg_dof_accessor.cc b/deal.II/deal.II/source/dofs/mg_dof_accessor.cc index 322b395cfe..d65ad36ca2 100644 --- a/deal.II/deal.II/source/dofs/mg_dof_accessor.cc +++ b/deal.II/deal.II/source/dofs/mg_dof_accessor.cc @@ -23,11 +23,9 @@ template MGDoFLineAccessor::MGDoFLineAccessor (Triangulation *tria, const int level, const int index, - const void *local_data) : - MGDoFAccessor ((MGDoFHandler*)local_data), - DoFLineAccessor(tria,level,index, - static_cast*> - (reinterpret_cast*>(local_data))) {}; + const AccessorData *local_data) : + MGDoFAccessor (local_data), + DoFLineAccessor(tria,level,index,local_data) {}; @@ -139,6 +137,15 @@ MGDoFLineAccessor::child (const unsigned int i) const { +template +void +MGDoFLineAccessor::copy_from (const MGDoFLineAccessor &a) { + DoFLineAccessor::copy_from (a); + set_mg_dof_handler (a.mg_dof_handler); +}; + + + /* ------------------------ MGDoFQuadAccessor --------------------------- */ @@ -146,11 +153,9 @@ template MGDoFQuadAccessor::MGDoFQuadAccessor (Triangulation *tria, const int level, const int index, - const void *local_data) : - MGDoFAccessor ((MGDoFHandler*)local_data), - DoFQuadAccessor(tria,level,index, - static_cast*> - (reinterpret_cast*>(local_data))) {}; + const AccessorData *local_data) : + MGDoFAccessor (local_data), + DoFQuadAccessor(tria,level,index,local_data) {}; @@ -282,6 +287,15 @@ MGDoFQuadAccessor::child (const unsigned int i) const { +template +void +MGDoFQuadAccessor::copy_from (const MGDoFQuadAccessor &a) { + DoFQuadAccessor::copy_from (a); + set_mg_dof_handler (a.mg_dof_handler); +}; + + + /*------------------------- Functions: MGDoFCellAccessor -----------------------*/ diff --git a/deal.II/deal.II/source/dofs/mg_dof_handler.cc b/deal.II/deal.II/source/dofs/mg_dof_handler.cc index 14b9a55a01..bb3d7d999f 100644 --- a/deal.II/deal.II/source/dofs/mg_dof_handler.cc +++ b/deal.II/deal.II/source/dofs/mg_dof_handler.cc @@ -2,7 +2,7 @@ #include #include -#include +#include #include #include #include @@ -16,12 +16,22 @@ /* ------------------------ MGVertexDoFs ----------------------------------- */ template -MGDoFHandler::MGVertexDoFs::MGVertexDoFs (const unsigned int coarsest_level, - const unsigned int n_levels, - const unsigned int dofs_per_vertex) : - coarsest_level (coarsest_level), - indices (new int[n_levels * dofs_per_vertex]) -{ +MGDoFHandler::MGVertexDoFs::MGVertexDoFs () : + coarsest_level (1<<30), + indices (0) +{}; + + + +template +void MGDoFHandler::MGVertexDoFs::init (const unsigned int cl, + const unsigned int n_levels, + const unsigned int dofs_per_vertex) { + Assert (indices == 0, ExcInternalError()); + + coarsest_level = cl; + + indices = new int[n_levels * dofs_per_vertex]; Assert (indices != 0, ExcNoMemory ()); for (unsigned int i=0; i::MGVertexDoFs::MGVertexDoFs (const unsigned int coarsest_level template MGDoFHandler::MGVertexDoFs::~MGVertexDoFs () { + Assert (indices != 0, ExcInternalError ()); + delete[] indices; }; @@ -56,6 +68,766 @@ MGDoFHandler::~MGDoFHandler () {}; + + +#if deal_II_dimension == 1 + +template <> +MGDoFHandler<1>::raw_cell_iterator +MGDoFHandler<1>::begin_raw (const unsigned int level) const { + return begin_raw_line (level); +}; + + + +template <> +MGDoFHandler<1>::cell_iterator +MGDoFHandler<1>::begin (const unsigned int level) const { + return begin_line (level); +}; + + + +template <> +MGDoFHandler<1>::active_cell_iterator +MGDoFHandler<1>::begin_active (const unsigned int level) const { + return begin_active_line (level); +}; + + + +template <> +MGDoFHandler<1>::raw_cell_iterator +MGDoFHandler<1>::end () const { + return end_line (); +}; + + + +template <> +MGDoFHandler<1>::raw_cell_iterator +MGDoFHandler<1>::last_raw () const { + return last_raw_line (); +}; + + + +template <> +MGDoFHandler<1>::raw_cell_iterator +MGDoFHandler<1>::last_raw (const unsigned int level) const { + return last_raw_line (level); +}; + + + +template <> +MGDoFHandler<1>::cell_iterator +MGDoFHandler<1>::last () const { + return last_line (); +}; + + + +template <> +MGDoFHandler<1>::cell_iterator +MGDoFHandler<1>::last (const unsigned int level) const { + return last_line (level); +}; + + + +template <> +MGDoFHandler<1>::active_cell_iterator +MGDoFHandler<1>::last_active () const { + return last_active_line (); +}; + + + +template <> +MGDoFHandler<1>::active_cell_iterator +MGDoFHandler<1>::last_active (const unsigned int level) const { + return last_active_line (level); +}; + + + +template <> +MGDoFDimensionInfo<1>::raw_face_iterator +MGDoFHandler<1>::begin_raw_face (const unsigned int) const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +template <> +MGDoFDimensionInfo<1>::face_iterator +MGDoFHandler<1>::begin_face (const unsigned int) const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +template <> +MGDoFDimensionInfo<1>::active_face_iterator +MGDoFHandler<1>::begin_active_face (const unsigned int) const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +template <> +MGDoFDimensionInfo<1>::raw_face_iterator +MGDoFHandler<1>::end_face () const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +template <> +MGDoFDimensionInfo<1>::raw_face_iterator +MGDoFHandler<1>::last_raw_face () const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +template <> +MGDoFDimensionInfo<1>::raw_face_iterator +MGDoFHandler<1>::last_raw_face (const unsigned int) const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +template <> +MGDoFDimensionInfo<1>::face_iterator +MGDoFHandler<1>::last_face () const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +template <> +MGDoFDimensionInfo<1>::face_iterator +MGDoFHandler<1>::last_face (const unsigned int) const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +template <> +MGDoFDimensionInfo<1>::active_face_iterator +MGDoFHandler<1>::last_active_face () const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +template <> +MGDoFDimensionInfo<1>::active_face_iterator +MGDoFHandler<1>::last_active_face (const unsigned int) const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + + +template <> +MGDoFHandler<1>::raw_quad_iterator +MGDoFHandler<1>::begin_raw_quad (const unsigned int) const { + Assert (false, ExcNotImplemented()); + return 0; +}; + + + +template <> +MGDoFHandler<1>::quad_iterator +MGDoFHandler<1>::begin_quad (const unsigned int) const { + Assert (false, ExcNotImplemented()); + return 0; +}; + + + +template <> +MGDoFHandler<1>::active_quad_iterator +MGDoFHandler<1>::begin_active_quad (const unsigned int) const { + Assert (false, ExcNotImplemented()); + return 0; +}; + + + +template <> +MGDoFHandler<1>::raw_quad_iterator +MGDoFHandler<1>::end_quad () const { + Assert (false, ExcNotImplemented()); + return 0; +}; + + + +template <> +MGDoFHandler<1>::raw_quad_iterator +MGDoFHandler<1>::last_raw_quad (const unsigned int) const { + Assert (false, ExcNotImplemented()); + return 0; +}; + + +template <> +MGDoFHandler<1>::quad_iterator +MGDoFHandler<1>::last_quad (const unsigned int) const { + Assert (false, ExcNotImplemented()); + return 0; +}; + + + +template <> +MGDoFHandler<1>::active_quad_iterator +MGDoFHandler<1>::last_active_quad (const unsigned int) const { + Assert (false, ExcNotImplemented()); + return 0; +}; + + + +template <> +MGDoFHandler<1>::raw_quad_iterator +MGDoFHandler<1>::last_raw_quad () const { + Assert (false, ExcNotImplemented()); + return 0; +}; + + + +template <> +MGDoFHandler<1>::quad_iterator +MGDoFHandler<1>::last_quad () const { + Assert (false, ExcNotImplemented()); + return 0; +}; + + + +template <> +MGDoFHandler<1>::active_quad_iterator +MGDoFHandler<1>::last_active_quad () const { + Assert (false, ExcNotImplemented()); + return 0; +}; + +#endif + + +#if deal_II_dimension == 2 + +template <> +MGDoFHandler<2>::raw_cell_iterator +MGDoFHandler<2>::begin_raw (const unsigned int level) const { + return begin_raw_quad (level); +}; + + + +template <> +MGDoFHandler<2>::cell_iterator +MGDoFHandler<2>::begin (const unsigned int level) const { + return begin_quad (level); +}; + + + +template <> +MGDoFHandler<2>::active_cell_iterator +MGDoFHandler<2>::begin_active (const unsigned int level) const { + return begin_active_quad (level); +}; + + + +template <> +MGDoFHandler<2>::raw_cell_iterator +MGDoFHandler<2>::end () const { + return end_quad (); +}; + + + +template <> +MGDoFHandler<2>::raw_cell_iterator +MGDoFHandler<2>::last_raw () const { + return last_raw_quad (); +}; + + + +template <> +MGDoFHandler<2>::raw_cell_iterator +MGDoFHandler<2>::last_raw (const unsigned int level) const { + return last_raw_quad (level); +}; + + + +template <> +MGDoFHandler<2>::cell_iterator +MGDoFHandler<2>::last () const { + return last_quad (); +}; + + + +template <> +MGDoFHandler<2>::cell_iterator +MGDoFHandler<2>::last (const unsigned int level) const { + return last_quad (level); +}; + + + +template <> +MGDoFHandler<2>::active_cell_iterator +MGDoFHandler<2>::last_active () const { + return last_active_quad (); +}; + + + +template <> +MGDoFHandler<2>::active_cell_iterator +MGDoFHandler<2>::last_active (const unsigned int level) const { + return last_active_quad (level); +}; + + +template <> +MGDoFDimensionInfo<2>::raw_face_iterator +MGDoFHandler<2>::begin_raw_face (const unsigned int level) const { + return begin_raw_line (level); +}; + + + +template <> +MGDoFDimensionInfo<2>::face_iterator +MGDoFHandler<2>::begin_face (const unsigned int level) const { + return begin_line (level); +}; + + + +template <> +MGDoFDimensionInfo<2>::active_face_iterator +MGDoFHandler<2>::begin_active_face (const unsigned int level) const { + return begin_active_line (level); +}; + + + +template <> +MGDoFDimensionInfo<2>::raw_face_iterator +MGDoFHandler<2>::end_face () const { + return end_line (); +}; + + + +template <> +MGDoFDimensionInfo<2>::raw_face_iterator +MGDoFHandler<2>::last_raw_face () const { + return last_raw_line (); +}; + + + +template <> +MGDoFDimensionInfo<2>::raw_face_iterator +MGDoFHandler<2>::last_raw_face (const unsigned int level) const { + return last_raw_line (level); +}; + + + +template <> +MGDoFDimensionInfo<2>::face_iterator +MGDoFHandler<2>::last_face () const { + return last_line (); +}; + + + +template <> +MGDoFDimensionInfo<2>::face_iterator +MGDoFHandler<2>::last_face (const unsigned int level) const { + return last_line (level); +}; + + + +template <> +MGDoFDimensionInfo<2>::active_face_iterator +MGDoFHandler<2>::last_active_face () const { + return last_active_line (); +}; + + + +template <> +MGDoFDimensionInfo<2>::active_face_iterator +MGDoFHandler<2>::last_active_face (const unsigned int level) const { + return last_active_line (level); +}; + +#endif + + + + +template +typename MGDoFHandler::raw_line_iterator +MGDoFHandler::begin_raw_line (const unsigned int level) const { + return raw_line_iterator (tria, + tria->begin_raw_line(level)->level(), + tria->begin_raw_line(level)->index(), + this); +}; + + + +template +typename MGDoFHandler::line_iterator +MGDoFHandler::begin_line (const unsigned int level) const { + return line_iterator (tria, + tria->begin_line(level)->level(), + tria->begin_line(level)->index(), + this); +}; + + + +template +typename MGDoFHandler::active_line_iterator +MGDoFHandler::begin_active_line (const unsigned int level) const { + return active_line_iterator (tria, + tria->begin_active_line(level)->level(), + tria->begin_active_line(level)->index(), + this); +}; + + + +template +typename MGDoFHandler::raw_quad_iterator +MGDoFHandler::begin_raw_quad (const unsigned int level) const { + return raw_quad_iterator (tria, + tria->begin_raw_quad(level)->level(), + tria->begin_raw_quad(level)->index(), + this); +}; + + + +template +typename MGDoFHandler::quad_iterator +MGDoFHandler::begin_quad (const unsigned int level) const { + return quad_iterator (tria, + tria->begin_quad(level)->level(), + tria->begin_quad(level)->index(), + this); +}; + + + +template +typename MGDoFHandler::active_quad_iterator +MGDoFHandler::begin_active_quad (const unsigned int level) const { + return active_quad_iterator (tria, + tria->begin_active_quad(level)->level(), + tria->begin_active_quad(level)->index(), + this); +}; + + + +template +typename MGDoFHandler::raw_line_iterator +MGDoFHandler::end_line () const { + return raw_line_iterator (tria, -1, -1, this); +}; + + + +template +typename MGDoFHandler::raw_quad_iterator +MGDoFHandler::end_quad () const { + return raw_quad_iterator (tria, -1, -1, this); +}; + + + +template +typename MGDoFDimensionInfo::raw_cell_iterator +MGDoFHandler::end_raw (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + end() : + begin_raw (level+1)); +}; + + + +template +typename MGDoFDimensionInfo::cell_iterator +MGDoFHandler::end (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + cell_iterator(end()) : + begin (level+1)); +}; + + + +template +typename MGDoFDimensionInfo::active_cell_iterator +MGDoFHandler::end_active (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + active_cell_iterator(end()) : + begin_active (level+1)); +}; + + + +template +typename MGDoFDimensionInfo::raw_face_iterator +MGDoFHandler::end_raw_face (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + end_face() : + begin_raw_face (level+1)); +}; + + + +template +typename MGDoFDimensionInfo::face_iterator +MGDoFHandler::end_face (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + face_iterator(end_face()) : + begin_face (level+1)); +}; + + + +template +typename MGDoFDimensionInfo::active_face_iterator +MGDoFHandler::end_active_face (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + active_face_iterator(end_face()) : + begin_active_face (level+1)); +}; + + + +template +typename MGDoFDimensionInfo::raw_line_iterator +MGDoFHandler::end_raw_line (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + end_line() : + begin_raw_line (level+1)); +}; + + + +template +typename MGDoFDimensionInfo::line_iterator +MGDoFHandler::end_line (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + line_iterator(end_line()) : + begin_line (level+1)); +}; + + + +template +typename MGDoFDimensionInfo::active_line_iterator +MGDoFHandler::end_active_line (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + active_line_iterator(end_line()) : + begin_active_line (level+1)); +}; + + + + +template +typename MGDoFDimensionInfo::raw_quad_iterator +MGDoFHandler::end_raw_quad (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + end_quad() : + begin_raw_quad (level+1)); +}; + + + +template +typename MGDoFDimensionInfo::quad_iterator +MGDoFHandler::end_quad (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + quad_iterator(end_quad()) : + begin_quad (level+1)); +}; + + + +template +typename MGDoFDimensionInfo::active_quad_iterator +MGDoFHandler::end_active_quad (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + active_quad_iterator(end_quad()) : + begin_active_quad (level+1)); +}; + + + +template +typename MGDoFHandler::raw_line_iterator +MGDoFHandler::last_raw_line (const unsigned int level) const { + return raw_line_iterator (tria, + tria->last_raw_line(level)->level(), + tria->last_raw_line(level)->index(), + this); +}; + + + +template +typename MGDoFHandler::line_iterator +MGDoFHandler::last_line (const unsigned int level) const { + return line_iterator (tria, + tria->last_line(level)->level(), + tria->last_line(level)->index(), + this); +}; + + + +template +typename MGDoFHandler::active_line_iterator +MGDoFHandler::last_active_line (const unsigned int level) const { + return active_line_iterator (tria, + tria->last_active_line(level)->level(), + tria->last_active_line(level)->index(), + this); +}; + + + +template +typename MGDoFHandler::raw_quad_iterator +MGDoFHandler::last_raw_quad (const unsigned int level) const { + return raw_quad_iterator (tria, + tria->last_raw_quad(level)->level(), + tria->last_raw_quad(level)->index(), + this); +}; + + + + +template +typename MGDoFHandler::quad_iterator +MGDoFHandler::last_quad (const unsigned int level) const { + return quad_iterator (tria, + tria->last_quad(level)->level(), + tria->last_quad(level)->index(), + this); +}; + + + + +template +typename MGDoFHandler::active_quad_iterator +MGDoFHandler::last_active_quad (const unsigned int level) const { + return active_quad_iterator (tria, + tria->last_active_quad(level)->level(), + tria->last_active_quad(level)->index(), + this); +}; + + + + +template +typename MGDoFHandler::raw_line_iterator +MGDoFHandler::last_raw_line () const { + return last_raw_line (mg_levels.size()-1); +}; + + + +template +typename MGDoFHandler::raw_quad_iterator +MGDoFHandler::last_raw_quad () const { + return last_raw_quad (mg_levels.size()-1); +}; + + + +template +typename MGDoFHandler::line_iterator +MGDoFHandler::last_line () const { + return last_line (mg_levels.size()-1); +}; + + + +template +typename MGDoFHandler::quad_iterator +MGDoFHandler::last_quad () const { + return last_quad (mg_levels.size()-1); +}; + + + +template +typename MGDoFHandler::active_line_iterator +MGDoFHandler::last_active_line () const { + return last_active_line (mg_levels.size()-1); +}; + + + + +template +typename MGDoFHandler::active_quad_iterator +MGDoFHandler::last_active_quad () const { + return last_active_quad (mg_levels.size()-1); +}; + + + + + + +//------------------------------------------------------------------ + + + + + + template void MGDoFHandler::distribute_dofs (const FiniteElementBase &fe) { // first distribute global dofs @@ -64,9 +836,131 @@ void MGDoFHandler::distribute_dofs (const FiniteElementBase &fe) { // reserve space for the MG dof numbers reserve_space (); + + // clear user flags because we will + // need them + tria->clear_user_flags (); + + // now distribute indices on each level + // separately + for (unsigned int level=0; leveln_levels(); ++level) + { + unsigned int next_free_dof = 0; + active_cell_iterator cell = begin(level), + endc = end(level); + + for (; cell != endc; ++cell) + next_free_dof = distribute_dofs_on_cell (cell, next_free_dof); + + mg_used_dofs[level] = next_free_dof; + }; +}; + + + +#if deal_II_dimension == 1 + +template <> +unsigned int +MGDoFHandler<1>::distribute_dofs_on_cell (active_cell_iterator &cell, + unsigned int next_free_dof) { + + // distribute dofs of vertices + if (selected_fe->dofs_per_vertex > 0) + for (unsigned int v=0; v::vertices_per_cell; ++v) + { + cell_iterator neighbor = cell->neighbor(v); + + if (neighbor.state() == valid) + { + // has neighbor already been processed? + if (neighbor->user_flag_set() && + (neighbor->level() == cell->level())) + // copy dofs if the neighbor is on + // the same level (only then are + // mg dofs the same) + { + if (v==0) + for (unsigned int d=0; ddofs_per_vertex; ++d) + cell->set_mg_vertex_dof_index (0, d, + neighbor->mg_vertex_dof_index (1, d)); + else + for (unsigned int d=0; ddofs_per_vertex; ++d) + cell->set_mg_vertex_dof_index (1, d, + neighbor->mg_vertex_dof_index (0, d)); + + // next neighbor + continue; + }; + }; + + // otherwise: create dofs newly + for (unsigned int d=0; ddofs_per_vertex; ++d) + cell->set_mg_vertex_dof_index (v, d, next_free_dof++); + }; + + // dofs of line + if (selected_fe->dofs_per_line > 0) + for (unsigned int d=0; ddofs_per_line; ++d) + cell->set_mg_dof_index (d, next_free_dof++); + + // note that this cell has been processed + cell->set_user_flag (); + + return next_free_dof; +}; + +#endif + + +#if deal_II_dimension == 2 + +template <> +unsigned int +MGDoFHandler<2>::distribute_dofs_on_cell (active_cell_iterator &cell, + unsigned int next_free_dof) { + if (selected_fe->dofs_per_vertex > 0) + // number dofs on vertices + for (unsigned int vertex=0; vertex::vertices_per_cell; ++vertex) + // check whether dofs for this + // vertex have been distributed + // (only check the first dof) + if (cell->mg_vertex_dof_index(vertex, 0) == -1) + for (unsigned int d=0; ddofs_per_vertex; ++d) + cell->set_mg_vertex_dof_index (vertex, d, next_free_dof++); + + // for the four sides + if (selected_fe->dofs_per_line > 0) + for (unsigned int side=0; side::faces_per_cell; ++side) + { + line_iterator line = cell->line(side); + + // distribute dofs if necessary: + // check whether line dof is already + // numbered (check only first dof) + if (line->mg_dof_index(0) == -1) + // if not: distribute dofs + for (unsigned int d=0; ddofs_per_line; ++d) + line->set_mg_dof_index (d, next_free_dof++); + }; + + + // dofs of quad + if (selected_fe->dofs_per_quad > 0) + for (unsigned int d=0; ddofs_per_quad; ++d) + cell->set_mg_dof_index (d, next_free_dof++); + + + // note that this cell has been processed + cell->set_user_flag (); + return next_free_dof; }; +#endif + + + @@ -74,6 +968,81 @@ void MGDoFHandler::distribute_dofs (const FiniteElementBase &fe) { template <> void MGDoFHandler<1>::reserve_space () { + const unsigned int dim = 1; + + Assert (selected_fe != 0, ExcNoFESelected()); + Assert (tria->n_levels() > 0, ExcInvalidTriangulation()); + + ////////////////////////// + // DESTRUCTION + + // delete all levels and set them up + // newly, since vectors are + // troublesome if you want to change + // their size + for (unsigned int i=0; in_levels(); ++i) + { + mg_levels.push_back (new DoFLevel<1>); + + mg_levels.back()->line_dofs = vector(tria->levels[i]->lines.lines.size() * + selected_fe->dofs_per_line, + -1); + }; + + // now allocate space for the + // vertices. To this end, we need + // to construct as many objects as + // there are vertices and let them + // allocate enough space for their + // vertex indices on the levels they + // live on. We need therefore to + // count to how many levels a cell + // belongs to, which we do by looping + // over all cells and storing the + // maximum and minimum level each + // vertex we pass by belongs to + mg_vertex_dofs.resize (tria->vertices.size()); + + vector min_level (tria->vertices.size(), tria->n_levels()); + vector max_level (tria->vertices.size(), 0); + + Triangulation::cell_iterator cell = tria->begin(), + endc = tria->end(); + for (; cell!=endc; ++cell) + for (unsigned int vertex=0; vertex::vertices_per_cell; + ++vertex) + { + const unsigned int vertex_index = cell->vertex_index(vertex); + if (min_level[vertex_index] > static_cast(cell->level())) + min_level[vertex_index] = cell->level(); + if (max_level[vertex_index] < static_cast(cell->level())) + max_level[vertex_index] = cell->level(); + }; + + // now allocate the needed space + for (unsigned int vertex=0; vertexvertices.size(); ++vertex) + { + Assert (min_level[vertex] < tria->n_levels(), ExcInternalError()); + Assert (max_level[vertex] >= min_level[vertex], ExcInternalError()); + + mg_vertex_dofs[vertex].init (min_level[vertex], + max_level[vertex]-min_level[vertex]+1, + selected_fe->dofs_per_vertex); + }; }; #endif @@ -83,6 +1052,85 @@ void MGDoFHandler<1>::reserve_space () { template <> void MGDoFHandler<2>::reserve_space () { + const unsigned int dim = 2; + + Assert (selected_fe != 0, ExcNoFESelected()); + Assert (tria->n_levels() > 0, ExcInvalidTriangulation()); + + //////////////////////////// + // DESTRUCTION + + // delete all levels and set them up + // newly, since vectors are + // troublesome if you want to change + // their size + for (unsigned int i=0; in_levels(); ++i) + { + mg_levels.push_back (new DoFLevel<2>); + + mg_levels.back()->line_dofs = vector (tria->levels[i]->lines.lines.size() * + selected_fe->dofs_per_line, + -1); + mg_levels.back()->quad_dofs = vector (tria->levels[i]->quads.quads.size() * + selected_fe->dofs_per_quad, + -1); + }; + + // now allocate space for the + // vertices. To this end, we need + // to construct as many objects as + // there are vertices and let them + // allocate enough space for their + // vertex indices on the levels they + // live on. We need therefore to + // count to how many levels a cell + // belongs to, which we do by looping + // over all cells and storing the + // maximum and minimum level each + // vertex we pass by belongs to + mg_vertex_dofs.resize (tria->vertices.size()); + + vector min_level (tria->vertices.size(), tria->n_levels()); + vector max_level (tria->vertices.size(), 0); + + Triangulation::cell_iterator cell = tria->begin(), + endc = tria->end(); + for (; cell!=endc; ++cell) + for (unsigned int vertex=0; vertex::vertices_per_cell; + ++vertex) + { + const unsigned int vertex_index = cell->vertex_index(vertex); + if (min_level[vertex_index] > static_cast(cell->level())) + min_level[vertex_index] = cell->level(); + if (max_level[vertex_index] < static_cast(cell->level())) + max_level[vertex_index] = cell->level(); + }; + + // now allocate the needed space + for (unsigned int vertex=0; vertexvertices.size(); ++vertex) + { + Assert (min_level[vertex] < tria->n_levels(), ExcInternalError()); + Assert (max_level[vertex] >= min_level[vertex], ExcInternalError()); + + mg_vertex_dofs[vertex].init (min_level[vertex], + max_level[vertex]-min_level[vertex]+1, + selected_fe->dofs_per_vertex); + }; }; #endif diff --git a/deal.II/deal.II/source/multigrid/mg_dof_accessor.cc b/deal.II/deal.II/source/multigrid/mg_dof_accessor.cc index 322b395cfe..d65ad36ca2 100644 --- a/deal.II/deal.II/source/multigrid/mg_dof_accessor.cc +++ b/deal.II/deal.II/source/multigrid/mg_dof_accessor.cc @@ -23,11 +23,9 @@ template MGDoFLineAccessor::MGDoFLineAccessor (Triangulation *tria, const int level, const int index, - const void *local_data) : - MGDoFAccessor ((MGDoFHandler*)local_data), - DoFLineAccessor(tria,level,index, - static_cast*> - (reinterpret_cast*>(local_data))) {}; + const AccessorData *local_data) : + MGDoFAccessor (local_data), + DoFLineAccessor(tria,level,index,local_data) {}; @@ -139,6 +137,15 @@ MGDoFLineAccessor::child (const unsigned int i) const { +template +void +MGDoFLineAccessor::copy_from (const MGDoFLineAccessor &a) { + DoFLineAccessor::copy_from (a); + set_mg_dof_handler (a.mg_dof_handler); +}; + + + /* ------------------------ MGDoFQuadAccessor --------------------------- */ @@ -146,11 +153,9 @@ template MGDoFQuadAccessor::MGDoFQuadAccessor (Triangulation *tria, const int level, const int index, - const void *local_data) : - MGDoFAccessor ((MGDoFHandler*)local_data), - DoFQuadAccessor(tria,level,index, - static_cast*> - (reinterpret_cast*>(local_data))) {}; + const AccessorData *local_data) : + MGDoFAccessor (local_data), + DoFQuadAccessor(tria,level,index,local_data) {}; @@ -282,6 +287,15 @@ MGDoFQuadAccessor::child (const unsigned int i) const { +template +void +MGDoFQuadAccessor::copy_from (const MGDoFQuadAccessor &a) { + DoFQuadAccessor::copy_from (a); + set_mg_dof_handler (a.mg_dof_handler); +}; + + + /*------------------------- Functions: MGDoFCellAccessor -----------------------*/ diff --git a/deal.II/deal.II/source/multigrid/mg_dof_handler.cc b/deal.II/deal.II/source/multigrid/mg_dof_handler.cc index 14b9a55a01..bb3d7d999f 100644 --- a/deal.II/deal.II/source/multigrid/mg_dof_handler.cc +++ b/deal.II/deal.II/source/multigrid/mg_dof_handler.cc @@ -2,7 +2,7 @@ #include #include -#include +#include #include #include #include @@ -16,12 +16,22 @@ /* ------------------------ MGVertexDoFs ----------------------------------- */ template -MGDoFHandler::MGVertexDoFs::MGVertexDoFs (const unsigned int coarsest_level, - const unsigned int n_levels, - const unsigned int dofs_per_vertex) : - coarsest_level (coarsest_level), - indices (new int[n_levels * dofs_per_vertex]) -{ +MGDoFHandler::MGVertexDoFs::MGVertexDoFs () : + coarsest_level (1<<30), + indices (0) +{}; + + + +template +void MGDoFHandler::MGVertexDoFs::init (const unsigned int cl, + const unsigned int n_levels, + const unsigned int dofs_per_vertex) { + Assert (indices == 0, ExcInternalError()); + + coarsest_level = cl; + + indices = new int[n_levels * dofs_per_vertex]; Assert (indices != 0, ExcNoMemory ()); for (unsigned int i=0; i::MGVertexDoFs::MGVertexDoFs (const unsigned int coarsest_level template MGDoFHandler::MGVertexDoFs::~MGVertexDoFs () { + Assert (indices != 0, ExcInternalError ()); + delete[] indices; }; @@ -56,6 +68,766 @@ MGDoFHandler::~MGDoFHandler () {}; + + +#if deal_II_dimension == 1 + +template <> +MGDoFHandler<1>::raw_cell_iterator +MGDoFHandler<1>::begin_raw (const unsigned int level) const { + return begin_raw_line (level); +}; + + + +template <> +MGDoFHandler<1>::cell_iterator +MGDoFHandler<1>::begin (const unsigned int level) const { + return begin_line (level); +}; + + + +template <> +MGDoFHandler<1>::active_cell_iterator +MGDoFHandler<1>::begin_active (const unsigned int level) const { + return begin_active_line (level); +}; + + + +template <> +MGDoFHandler<1>::raw_cell_iterator +MGDoFHandler<1>::end () const { + return end_line (); +}; + + + +template <> +MGDoFHandler<1>::raw_cell_iterator +MGDoFHandler<1>::last_raw () const { + return last_raw_line (); +}; + + + +template <> +MGDoFHandler<1>::raw_cell_iterator +MGDoFHandler<1>::last_raw (const unsigned int level) const { + return last_raw_line (level); +}; + + + +template <> +MGDoFHandler<1>::cell_iterator +MGDoFHandler<1>::last () const { + return last_line (); +}; + + + +template <> +MGDoFHandler<1>::cell_iterator +MGDoFHandler<1>::last (const unsigned int level) const { + return last_line (level); +}; + + + +template <> +MGDoFHandler<1>::active_cell_iterator +MGDoFHandler<1>::last_active () const { + return last_active_line (); +}; + + + +template <> +MGDoFHandler<1>::active_cell_iterator +MGDoFHandler<1>::last_active (const unsigned int level) const { + return last_active_line (level); +}; + + + +template <> +MGDoFDimensionInfo<1>::raw_face_iterator +MGDoFHandler<1>::begin_raw_face (const unsigned int) const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +template <> +MGDoFDimensionInfo<1>::face_iterator +MGDoFHandler<1>::begin_face (const unsigned int) const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +template <> +MGDoFDimensionInfo<1>::active_face_iterator +MGDoFHandler<1>::begin_active_face (const unsigned int) const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +template <> +MGDoFDimensionInfo<1>::raw_face_iterator +MGDoFHandler<1>::end_face () const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +template <> +MGDoFDimensionInfo<1>::raw_face_iterator +MGDoFHandler<1>::last_raw_face () const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +template <> +MGDoFDimensionInfo<1>::raw_face_iterator +MGDoFHandler<1>::last_raw_face (const unsigned int) const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +template <> +MGDoFDimensionInfo<1>::face_iterator +MGDoFHandler<1>::last_face () const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +template <> +MGDoFDimensionInfo<1>::face_iterator +MGDoFHandler<1>::last_face (const unsigned int) const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +template <> +MGDoFDimensionInfo<1>::active_face_iterator +MGDoFHandler<1>::last_active_face () const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + +template <> +MGDoFDimensionInfo<1>::active_face_iterator +MGDoFHandler<1>::last_active_face (const unsigned int) const { + Assert (false, ExcFunctionNotUseful()); + return 0; +}; + + + + +template <> +MGDoFHandler<1>::raw_quad_iterator +MGDoFHandler<1>::begin_raw_quad (const unsigned int) const { + Assert (false, ExcNotImplemented()); + return 0; +}; + + + +template <> +MGDoFHandler<1>::quad_iterator +MGDoFHandler<1>::begin_quad (const unsigned int) const { + Assert (false, ExcNotImplemented()); + return 0; +}; + + + +template <> +MGDoFHandler<1>::active_quad_iterator +MGDoFHandler<1>::begin_active_quad (const unsigned int) const { + Assert (false, ExcNotImplemented()); + return 0; +}; + + + +template <> +MGDoFHandler<1>::raw_quad_iterator +MGDoFHandler<1>::end_quad () const { + Assert (false, ExcNotImplemented()); + return 0; +}; + + + +template <> +MGDoFHandler<1>::raw_quad_iterator +MGDoFHandler<1>::last_raw_quad (const unsigned int) const { + Assert (false, ExcNotImplemented()); + return 0; +}; + + +template <> +MGDoFHandler<1>::quad_iterator +MGDoFHandler<1>::last_quad (const unsigned int) const { + Assert (false, ExcNotImplemented()); + return 0; +}; + + + +template <> +MGDoFHandler<1>::active_quad_iterator +MGDoFHandler<1>::last_active_quad (const unsigned int) const { + Assert (false, ExcNotImplemented()); + return 0; +}; + + + +template <> +MGDoFHandler<1>::raw_quad_iterator +MGDoFHandler<1>::last_raw_quad () const { + Assert (false, ExcNotImplemented()); + return 0; +}; + + + +template <> +MGDoFHandler<1>::quad_iterator +MGDoFHandler<1>::last_quad () const { + Assert (false, ExcNotImplemented()); + return 0; +}; + + + +template <> +MGDoFHandler<1>::active_quad_iterator +MGDoFHandler<1>::last_active_quad () const { + Assert (false, ExcNotImplemented()); + return 0; +}; + +#endif + + +#if deal_II_dimension == 2 + +template <> +MGDoFHandler<2>::raw_cell_iterator +MGDoFHandler<2>::begin_raw (const unsigned int level) const { + return begin_raw_quad (level); +}; + + + +template <> +MGDoFHandler<2>::cell_iterator +MGDoFHandler<2>::begin (const unsigned int level) const { + return begin_quad (level); +}; + + + +template <> +MGDoFHandler<2>::active_cell_iterator +MGDoFHandler<2>::begin_active (const unsigned int level) const { + return begin_active_quad (level); +}; + + + +template <> +MGDoFHandler<2>::raw_cell_iterator +MGDoFHandler<2>::end () const { + return end_quad (); +}; + + + +template <> +MGDoFHandler<2>::raw_cell_iterator +MGDoFHandler<2>::last_raw () const { + return last_raw_quad (); +}; + + + +template <> +MGDoFHandler<2>::raw_cell_iterator +MGDoFHandler<2>::last_raw (const unsigned int level) const { + return last_raw_quad (level); +}; + + + +template <> +MGDoFHandler<2>::cell_iterator +MGDoFHandler<2>::last () const { + return last_quad (); +}; + + + +template <> +MGDoFHandler<2>::cell_iterator +MGDoFHandler<2>::last (const unsigned int level) const { + return last_quad (level); +}; + + + +template <> +MGDoFHandler<2>::active_cell_iterator +MGDoFHandler<2>::last_active () const { + return last_active_quad (); +}; + + + +template <> +MGDoFHandler<2>::active_cell_iterator +MGDoFHandler<2>::last_active (const unsigned int level) const { + return last_active_quad (level); +}; + + +template <> +MGDoFDimensionInfo<2>::raw_face_iterator +MGDoFHandler<2>::begin_raw_face (const unsigned int level) const { + return begin_raw_line (level); +}; + + + +template <> +MGDoFDimensionInfo<2>::face_iterator +MGDoFHandler<2>::begin_face (const unsigned int level) const { + return begin_line (level); +}; + + + +template <> +MGDoFDimensionInfo<2>::active_face_iterator +MGDoFHandler<2>::begin_active_face (const unsigned int level) const { + return begin_active_line (level); +}; + + + +template <> +MGDoFDimensionInfo<2>::raw_face_iterator +MGDoFHandler<2>::end_face () const { + return end_line (); +}; + + + +template <> +MGDoFDimensionInfo<2>::raw_face_iterator +MGDoFHandler<2>::last_raw_face () const { + return last_raw_line (); +}; + + + +template <> +MGDoFDimensionInfo<2>::raw_face_iterator +MGDoFHandler<2>::last_raw_face (const unsigned int level) const { + return last_raw_line (level); +}; + + + +template <> +MGDoFDimensionInfo<2>::face_iterator +MGDoFHandler<2>::last_face () const { + return last_line (); +}; + + + +template <> +MGDoFDimensionInfo<2>::face_iterator +MGDoFHandler<2>::last_face (const unsigned int level) const { + return last_line (level); +}; + + + +template <> +MGDoFDimensionInfo<2>::active_face_iterator +MGDoFHandler<2>::last_active_face () const { + return last_active_line (); +}; + + + +template <> +MGDoFDimensionInfo<2>::active_face_iterator +MGDoFHandler<2>::last_active_face (const unsigned int level) const { + return last_active_line (level); +}; + +#endif + + + + +template +typename MGDoFHandler::raw_line_iterator +MGDoFHandler::begin_raw_line (const unsigned int level) const { + return raw_line_iterator (tria, + tria->begin_raw_line(level)->level(), + tria->begin_raw_line(level)->index(), + this); +}; + + + +template +typename MGDoFHandler::line_iterator +MGDoFHandler::begin_line (const unsigned int level) const { + return line_iterator (tria, + tria->begin_line(level)->level(), + tria->begin_line(level)->index(), + this); +}; + + + +template +typename MGDoFHandler::active_line_iterator +MGDoFHandler::begin_active_line (const unsigned int level) const { + return active_line_iterator (tria, + tria->begin_active_line(level)->level(), + tria->begin_active_line(level)->index(), + this); +}; + + + +template +typename MGDoFHandler::raw_quad_iterator +MGDoFHandler::begin_raw_quad (const unsigned int level) const { + return raw_quad_iterator (tria, + tria->begin_raw_quad(level)->level(), + tria->begin_raw_quad(level)->index(), + this); +}; + + + +template +typename MGDoFHandler::quad_iterator +MGDoFHandler::begin_quad (const unsigned int level) const { + return quad_iterator (tria, + tria->begin_quad(level)->level(), + tria->begin_quad(level)->index(), + this); +}; + + + +template +typename MGDoFHandler::active_quad_iterator +MGDoFHandler::begin_active_quad (const unsigned int level) const { + return active_quad_iterator (tria, + tria->begin_active_quad(level)->level(), + tria->begin_active_quad(level)->index(), + this); +}; + + + +template +typename MGDoFHandler::raw_line_iterator +MGDoFHandler::end_line () const { + return raw_line_iterator (tria, -1, -1, this); +}; + + + +template +typename MGDoFHandler::raw_quad_iterator +MGDoFHandler::end_quad () const { + return raw_quad_iterator (tria, -1, -1, this); +}; + + + +template +typename MGDoFDimensionInfo::raw_cell_iterator +MGDoFHandler::end_raw (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + end() : + begin_raw (level+1)); +}; + + + +template +typename MGDoFDimensionInfo::cell_iterator +MGDoFHandler::end (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + cell_iterator(end()) : + begin (level+1)); +}; + + + +template +typename MGDoFDimensionInfo::active_cell_iterator +MGDoFHandler::end_active (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + active_cell_iterator(end()) : + begin_active (level+1)); +}; + + + +template +typename MGDoFDimensionInfo::raw_face_iterator +MGDoFHandler::end_raw_face (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + end_face() : + begin_raw_face (level+1)); +}; + + + +template +typename MGDoFDimensionInfo::face_iterator +MGDoFHandler::end_face (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + face_iterator(end_face()) : + begin_face (level+1)); +}; + + + +template +typename MGDoFDimensionInfo::active_face_iterator +MGDoFHandler::end_active_face (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + active_face_iterator(end_face()) : + begin_active_face (level+1)); +}; + + + +template +typename MGDoFDimensionInfo::raw_line_iterator +MGDoFHandler::end_raw_line (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + end_line() : + begin_raw_line (level+1)); +}; + + + +template +typename MGDoFDimensionInfo::line_iterator +MGDoFHandler::end_line (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + line_iterator(end_line()) : + begin_line (level+1)); +}; + + + +template +typename MGDoFDimensionInfo::active_line_iterator +MGDoFHandler::end_active_line (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + active_line_iterator(end_line()) : + begin_active_line (level+1)); +}; + + + + +template +typename MGDoFDimensionInfo::raw_quad_iterator +MGDoFHandler::end_raw_quad (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + end_quad() : + begin_raw_quad (level+1)); +}; + + + +template +typename MGDoFDimensionInfo::quad_iterator +MGDoFHandler::end_quad (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + quad_iterator(end_quad()) : + begin_quad (level+1)); +}; + + + +template +typename MGDoFDimensionInfo::active_quad_iterator +MGDoFHandler::end_active_quad (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + active_quad_iterator(end_quad()) : + begin_active_quad (level+1)); +}; + + + +template +typename MGDoFHandler::raw_line_iterator +MGDoFHandler::last_raw_line (const unsigned int level) const { + return raw_line_iterator (tria, + tria->last_raw_line(level)->level(), + tria->last_raw_line(level)->index(), + this); +}; + + + +template +typename MGDoFHandler::line_iterator +MGDoFHandler::last_line (const unsigned int level) const { + return line_iterator (tria, + tria->last_line(level)->level(), + tria->last_line(level)->index(), + this); +}; + + + +template +typename MGDoFHandler::active_line_iterator +MGDoFHandler::last_active_line (const unsigned int level) const { + return active_line_iterator (tria, + tria->last_active_line(level)->level(), + tria->last_active_line(level)->index(), + this); +}; + + + +template +typename MGDoFHandler::raw_quad_iterator +MGDoFHandler::last_raw_quad (const unsigned int level) const { + return raw_quad_iterator (tria, + tria->last_raw_quad(level)->level(), + tria->last_raw_quad(level)->index(), + this); +}; + + + + +template +typename MGDoFHandler::quad_iterator +MGDoFHandler::last_quad (const unsigned int level) const { + return quad_iterator (tria, + tria->last_quad(level)->level(), + tria->last_quad(level)->index(), + this); +}; + + + + +template +typename MGDoFHandler::active_quad_iterator +MGDoFHandler::last_active_quad (const unsigned int level) const { + return active_quad_iterator (tria, + tria->last_active_quad(level)->level(), + tria->last_active_quad(level)->index(), + this); +}; + + + + +template +typename MGDoFHandler::raw_line_iterator +MGDoFHandler::last_raw_line () const { + return last_raw_line (mg_levels.size()-1); +}; + + + +template +typename MGDoFHandler::raw_quad_iterator +MGDoFHandler::last_raw_quad () const { + return last_raw_quad (mg_levels.size()-1); +}; + + + +template +typename MGDoFHandler::line_iterator +MGDoFHandler::last_line () const { + return last_line (mg_levels.size()-1); +}; + + + +template +typename MGDoFHandler::quad_iterator +MGDoFHandler::last_quad () const { + return last_quad (mg_levels.size()-1); +}; + + + +template +typename MGDoFHandler::active_line_iterator +MGDoFHandler::last_active_line () const { + return last_active_line (mg_levels.size()-1); +}; + + + + +template +typename MGDoFHandler::active_quad_iterator +MGDoFHandler::last_active_quad () const { + return last_active_quad (mg_levels.size()-1); +}; + + + + + + +//------------------------------------------------------------------ + + + + + + template void MGDoFHandler::distribute_dofs (const FiniteElementBase &fe) { // first distribute global dofs @@ -64,9 +836,131 @@ void MGDoFHandler::distribute_dofs (const FiniteElementBase &fe) { // reserve space for the MG dof numbers reserve_space (); + + // clear user flags because we will + // need them + tria->clear_user_flags (); + + // now distribute indices on each level + // separately + for (unsigned int level=0; leveln_levels(); ++level) + { + unsigned int next_free_dof = 0; + active_cell_iterator cell = begin(level), + endc = end(level); + + for (; cell != endc; ++cell) + next_free_dof = distribute_dofs_on_cell (cell, next_free_dof); + + mg_used_dofs[level] = next_free_dof; + }; +}; + + + +#if deal_II_dimension == 1 + +template <> +unsigned int +MGDoFHandler<1>::distribute_dofs_on_cell (active_cell_iterator &cell, + unsigned int next_free_dof) { + + // distribute dofs of vertices + if (selected_fe->dofs_per_vertex > 0) + for (unsigned int v=0; v::vertices_per_cell; ++v) + { + cell_iterator neighbor = cell->neighbor(v); + + if (neighbor.state() == valid) + { + // has neighbor already been processed? + if (neighbor->user_flag_set() && + (neighbor->level() == cell->level())) + // copy dofs if the neighbor is on + // the same level (only then are + // mg dofs the same) + { + if (v==0) + for (unsigned int d=0; ddofs_per_vertex; ++d) + cell->set_mg_vertex_dof_index (0, d, + neighbor->mg_vertex_dof_index (1, d)); + else + for (unsigned int d=0; ddofs_per_vertex; ++d) + cell->set_mg_vertex_dof_index (1, d, + neighbor->mg_vertex_dof_index (0, d)); + + // next neighbor + continue; + }; + }; + + // otherwise: create dofs newly + for (unsigned int d=0; ddofs_per_vertex; ++d) + cell->set_mg_vertex_dof_index (v, d, next_free_dof++); + }; + + // dofs of line + if (selected_fe->dofs_per_line > 0) + for (unsigned int d=0; ddofs_per_line; ++d) + cell->set_mg_dof_index (d, next_free_dof++); + + // note that this cell has been processed + cell->set_user_flag (); + + return next_free_dof; +}; + +#endif + + +#if deal_II_dimension == 2 + +template <> +unsigned int +MGDoFHandler<2>::distribute_dofs_on_cell (active_cell_iterator &cell, + unsigned int next_free_dof) { + if (selected_fe->dofs_per_vertex > 0) + // number dofs on vertices + for (unsigned int vertex=0; vertex::vertices_per_cell; ++vertex) + // check whether dofs for this + // vertex have been distributed + // (only check the first dof) + if (cell->mg_vertex_dof_index(vertex, 0) == -1) + for (unsigned int d=0; ddofs_per_vertex; ++d) + cell->set_mg_vertex_dof_index (vertex, d, next_free_dof++); + + // for the four sides + if (selected_fe->dofs_per_line > 0) + for (unsigned int side=0; side::faces_per_cell; ++side) + { + line_iterator line = cell->line(side); + + // distribute dofs if necessary: + // check whether line dof is already + // numbered (check only first dof) + if (line->mg_dof_index(0) == -1) + // if not: distribute dofs + for (unsigned int d=0; ddofs_per_line; ++d) + line->set_mg_dof_index (d, next_free_dof++); + }; + + + // dofs of quad + if (selected_fe->dofs_per_quad > 0) + for (unsigned int d=0; ddofs_per_quad; ++d) + cell->set_mg_dof_index (d, next_free_dof++); + + + // note that this cell has been processed + cell->set_user_flag (); + return next_free_dof; }; +#endif + + + @@ -74,6 +968,81 @@ void MGDoFHandler::distribute_dofs (const FiniteElementBase &fe) { template <> void MGDoFHandler<1>::reserve_space () { + const unsigned int dim = 1; + + Assert (selected_fe != 0, ExcNoFESelected()); + Assert (tria->n_levels() > 0, ExcInvalidTriangulation()); + + ////////////////////////// + // DESTRUCTION + + // delete all levels and set them up + // newly, since vectors are + // troublesome if you want to change + // their size + for (unsigned int i=0; in_levels(); ++i) + { + mg_levels.push_back (new DoFLevel<1>); + + mg_levels.back()->line_dofs = vector(tria->levels[i]->lines.lines.size() * + selected_fe->dofs_per_line, + -1); + }; + + // now allocate space for the + // vertices. To this end, we need + // to construct as many objects as + // there are vertices and let them + // allocate enough space for their + // vertex indices on the levels they + // live on. We need therefore to + // count to how many levels a cell + // belongs to, which we do by looping + // over all cells and storing the + // maximum and minimum level each + // vertex we pass by belongs to + mg_vertex_dofs.resize (tria->vertices.size()); + + vector min_level (tria->vertices.size(), tria->n_levels()); + vector max_level (tria->vertices.size(), 0); + + Triangulation::cell_iterator cell = tria->begin(), + endc = tria->end(); + for (; cell!=endc; ++cell) + for (unsigned int vertex=0; vertex::vertices_per_cell; + ++vertex) + { + const unsigned int vertex_index = cell->vertex_index(vertex); + if (min_level[vertex_index] > static_cast(cell->level())) + min_level[vertex_index] = cell->level(); + if (max_level[vertex_index] < static_cast(cell->level())) + max_level[vertex_index] = cell->level(); + }; + + // now allocate the needed space + for (unsigned int vertex=0; vertexvertices.size(); ++vertex) + { + Assert (min_level[vertex] < tria->n_levels(), ExcInternalError()); + Assert (max_level[vertex] >= min_level[vertex], ExcInternalError()); + + mg_vertex_dofs[vertex].init (min_level[vertex], + max_level[vertex]-min_level[vertex]+1, + selected_fe->dofs_per_vertex); + }; }; #endif @@ -83,6 +1052,85 @@ void MGDoFHandler<1>::reserve_space () { template <> void MGDoFHandler<2>::reserve_space () { + const unsigned int dim = 2; + + Assert (selected_fe != 0, ExcNoFESelected()); + Assert (tria->n_levels() > 0, ExcInvalidTriangulation()); + + //////////////////////////// + // DESTRUCTION + + // delete all levels and set them up + // newly, since vectors are + // troublesome if you want to change + // their size + for (unsigned int i=0; in_levels(); ++i) + { + mg_levels.push_back (new DoFLevel<2>); + + mg_levels.back()->line_dofs = vector (tria->levels[i]->lines.lines.size() * + selected_fe->dofs_per_line, + -1); + mg_levels.back()->quad_dofs = vector (tria->levels[i]->quads.quads.size() * + selected_fe->dofs_per_quad, + -1); + }; + + // now allocate space for the + // vertices. To this end, we need + // to construct as many objects as + // there are vertices and let them + // allocate enough space for their + // vertex indices on the levels they + // live on. We need therefore to + // count to how many levels a cell + // belongs to, which we do by looping + // over all cells and storing the + // maximum and minimum level each + // vertex we pass by belongs to + mg_vertex_dofs.resize (tria->vertices.size()); + + vector min_level (tria->vertices.size(), tria->n_levels()); + vector max_level (tria->vertices.size(), 0); + + Triangulation::cell_iterator cell = tria->begin(), + endc = tria->end(); + for (; cell!=endc; ++cell) + for (unsigned int vertex=0; vertex::vertices_per_cell; + ++vertex) + { + const unsigned int vertex_index = cell->vertex_index(vertex); + if (min_level[vertex_index] > static_cast(cell->level())) + min_level[vertex_index] = cell->level(); + if (max_level[vertex_index] < static_cast(cell->level())) + max_level[vertex_index] = cell->level(); + }; + + // now allocate the needed space + for (unsigned int vertex=0; vertexvertices.size(); ++vertex) + { + Assert (min_level[vertex] < tria->n_levels(), ExcInternalError()); + Assert (max_level[vertex] >= min_level[vertex], ExcInternalError()); + + mg_vertex_dofs[vertex].init (min_level[vertex], + max_level[vertex]-min_level[vertex]+1, + selected_fe->dofs_per_vertex); + }; }; #endif