From 1b7aa2a8651c25100fa8b6ff13cdc0d1b9526069 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 25 Aug 1998 16:19:05 +0000 Subject: [PATCH] Initial version of the multigrid components. git-svn-id: https://svn.dealii.org/trunk@518 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/include/dofs/mg_dof_accessor.h | 505 ++++++++++++++++++ deal.II/deal.II/include/dofs/mg_dof_handler.h | 253 +++++++++ .../include/multigrid/mg_dof_accessor.h | 505 ++++++++++++++++++ .../include/multigrid/mg_dof_handler.h | 253 +++++++++ .../deal.II/source/dofs/mg_dof_accessor.cc | 426 +++++++++++++++ deal.II/deal.II/source/dofs/mg_dof_handler.cc | 92 ++++ .../source/multigrid/mg_dof_accessor.cc | 426 +++++++++++++++ .../source/multigrid/mg_dof_handler.cc | 92 ++++ 8 files changed, 2552 insertions(+) create mode 100644 deal.II/deal.II/include/dofs/mg_dof_accessor.h create mode 100644 deal.II/deal.II/include/dofs/mg_dof_handler.h create mode 100644 deal.II/deal.II/include/multigrid/mg_dof_accessor.h create mode 100644 deal.II/deal.II/include/multigrid/mg_dof_handler.h create mode 100644 deal.II/deal.II/source/dofs/mg_dof_accessor.cc create mode 100644 deal.II/deal.II/source/dofs/mg_dof_handler.cc create mode 100644 deal.II/deal.II/source/multigrid/mg_dof_accessor.cc create mode 100644 deal.II/deal.II/source/multigrid/mg_dof_handler.cc diff --git a/deal.II/deal.II/include/dofs/mg_dof_accessor.h b/deal.II/deal.II/include/dofs/mg_dof_accessor.h new file mode 100644 index 0000000000..03274755dc --- /dev/null +++ b/deal.II/deal.II/include/dofs/mg_dof_accessor.h @@ -0,0 +1,505 @@ +/*---------------------------- mg_dof_accessor.h ---------------------------*/ +/* $Id$ */ +#ifndef __mg_dof_accessor_H +#define __mg_dof_accessor_H +/*---------------------------- mg_dof_accessor.h ---------------------------*/ + + +#include + + +// forward declaration +template class MGDoFHandler; + + + + +/** + * Define the basis for accessors to the degrees of freedom for + * a multigrid DoF handler object. + * + * Note that it is allowed to construct an object of which the + * #mg_dof_handler# pointer is a Null pointer. Such an object would + * result in a strange kind of behaviour, though every reasonable + * operating system should disallow access through that pointer. + * The reason we do not check for the null pointer in the + * constructor which gets passed the #DoFHandler# pointer is that + * if we did we could not make dof iterators member of other classes + * (like in the #FEValues# class) if we did not know about the + * #DoFHandler# object to be used upon construction of that object. + * Through the way this class is implemented here, we allow the + * creation of a kind of virgin object which only gets useful if + * assigned to from another object before first usage. + * + * Opposite to construction, it is not possible to copy an object + * which has an invalid dof handler pointer. This is to guarantee + * that every iterator which is once assigned to is a valid + * object. However, this assertion only holds in debug mode, when + * the #Assert# macro is switched on. + * + * @author Wolfgang Bangerth, 1998 + */ +template +class MGDoFAccessor { + public: + /** + * Constructor + */ + MGDoFAccessor () : mg_dof_handler(0) { + Assert (false, ExcInvalidObject()); + }; + + /** + * This should be the default constructor. + */ + MGDoFAccessor (MGDoFHandler *mg_dof_handler) : + mg_dof_handler(mg_dof_handler) {}; + + /** + * Reset the DoF handler pointer. + */ + void set_mg_dof_handler (MGDoFHandler *dh) { + Assert (dh != 0, ExcInvalidObject()); + mg_dof_handler = dh; + }; + + /** + * Copy operator. + */ + MGDoFAccessor & operator = (const MGDoFAccessor &da) { + set_dof_handler (da.mg_dof_handler); + return *this; + }; + + protected: + /** + * Store the address of the #MGDoFHandler# object + * to be accessed. + */ + MGDoFHandler *mg_dof_handler; +}; + + + + + +/** + * Grant access to the degrees of freedom located on lines. + * This class follows mainly the route laid out by the accessor library + * declared in the triangulation library (\Ref{TriaAccessor}). It enables + * the user to access the degrees of freedom on the lines (there are similar + * versions for the DoFs on quads, etc), where the dimension of the underlying + * triangulation does not really matter (i.e. this accessor works with the + * lines in 1D-, 2D-, etc dimensions). + * + * + * \subsection{Usage} + * + * The \Ref{DoFDimensionInfo} classes inherited by the \Ref{DoFHandler} classes + * declare typedefs to iterators using the accessors declared in this class + * hierarchy tree. Usage is best to happens through these typedefs, since they + * are more secure to changes in the class naming and template interface as well + * as they provide easier typing (much less complicated names!). + * + * + * \subsection{Notes about the class hierarchy structure} + * + * The class hierarchy seems to be a bit confused here. The reason for this is + * that we would really like to derive a #DoFLineAccessor# from a #LineAccessor#. + * Unfortunately, we would run into problems, if we wanted a #DoFLineAccessor# + * in one spatial dimension, in which case a line is also a cell. The traditional + * solution would be to declare a #DoFCellAccessor<1># which is derived from + * #DoFLineAccessor<1># and #CellAccessor<1># (the #DoFLineAccessor# cannot + * itself be derived from #CellAccessor# since a line is not a cell + * unless in one space dimension), but since a #DoFLineAccessor# and a + * #CellAccessor# are both derived from #TriaAccessor#, we would have to make + * the last derivation virtual. + * + * Since we want to avoid virtual inheritance since this involves another + * indirection in every member variable access, we chose another way: we + * pass a second template parameter to a #DoFLineAccessor# which tells it + * which class to be derived from: if we are in one spatial dimension, the + * base class is to be #CellAccessor<1>#, in two or more dimensions it + * is a #LineAccessor#, i.e. am accessor to lines without the missing + * functionality needed for cells (neighbors, etc.). + * + * This way we can declare a #DoFCellAccessor# in one dimension by deriving + * from #DoFLineAccessor<1,CellAccessor<1> >#, thus getting the cell + * functionality through the #DoFLineAccessor# instead of through a virtual + * multiple inheritance of #DoFLineAccessor# and #CellAccessor<1>#. + * + * The same concept is used with #DoFQuadAccessor# classes etc. + * + * @author Wolfgang Bangerth, 1998 + */ +template +class MGDoFLineAccessor : public MGDoFAccessor, + public DoFLineAccessor { + public: + /** + * Default constructor, unused thus + * not implemented. + */ + MGDoFLineAccessor (); + + /** + * Constructor. The #local_data# + * argument is assumed to be a pointer + * to a #MGDoFHandler# object. + */ + MGDoFLineAccessor (Triangulation *tria, + const int level, + const int index, + const void *local_data); + + /** + * Return the index of the #i#th degree + * of freedom of this line on the level + * this line lives on. + */ + int mg_dof_index (const unsigned int i) const; + + /** + * Set the index of the #i#th degree + * of freedom of this line on the + * level this line lives on to #index#. + */ + void set_mg_dof_index (const unsigned int i, const int index) const; + + /** + * Return the index of the #i#th degree + * on the #vertex#th vertex for the + * level this line lives on. + */ + int mg_vertex_dof_index (const unsigned int vertex, + const unsigned int i) const; + + /** + * Set the index of the #i#th degree + * on the #vertex#th vertex to #index# + * for the level this line lives on. + */ + void set_mg_vertex_dof_index (const unsigned int vertex, + const unsigned int i, + const int index) const; + + /** + * Return the indices of the dofs of this + * line in the standard ordering: dofs + * on vertex 0, dofs on vertex 1, + * dofs on line. + * + * It is assumed that the vector already + * has the right size beforehand. The + * indices refer to the local numbering + * for the level this line lives on. + */ + void get_mg_dof_indices (vector &dof_indices) const; + + /** + * Return the #i#th child as a MGDoF line + * iterator. This function is needed since + * the child function of the base + * class returns a line accessor without + * access to the MG data. + */ + TriaIterator > child (const unsigned int) const; + + /** + * Implement the copy operator needed + * for the iterator classes. + */ + void copy_from (const MGDoFLineAccessor &a); +}; + + + + + + +/** + * Grant access to the multilevel degrees of freedom located on quads. + * + * @see DoFLineAccessor + */ +template +class MGDoFQuadAccessor : public MGDoFAccessor, + public DoFQuadAccessor { + public: + /** + * Default constructor, unused thus + * not implemented. + */ + MGDoFQuadAccessor (); + + /** + * Constructor. The #local_data# + * argument is assumed to be a pointer + * to a #DoFHandler# object. + */ + MGDoFQuadAccessor (Triangulation *tria, + const int level, + const int index, + const void *local_data); + + /** + * Return the index of the #i#th degree + * of freedom of this quad on the level + * this quad lives on. + */ + int mg_dof_index (const unsigned int i) const; + + /** + * Set the index of the #i#th degree + * of freedom of this quad on the + * level this quad lives on to #index#. + */ + void set_mg_dof_index (const unsigned int i, const int index) const; + + /** + * Return the index of the #i#th degree + * on the #vertex#th vertex for the level + * this quad lives on. + */ + int mg_vertex_dof_index (const unsigned int vertex, + const unsigned int i) const; + + /** + * Set the index of the #i#th degree + * on the #vertex#th vertex for the + * level this quad lives on to #index#. + */ + void set_mg_vertex_dof_index (const unsigned int vertex, + const unsigned int i, + const int index) const; + + /** + * Return the indices of the dofs of this + * quad in the standard ordering: dofs + * on vertex 0, dofs on vertex 1, etc, + * dofs on line 0, dofs on line 1, etc, + * dofs on quad 0, etc. + * + * It is assumed that the vector already + * has the right size beforehand. The + * indices refer to the local numbering + * for the level this quad lives on. + */ + void get_mg_dof_indices (vector &dof_indices) const; + + /** + * Return a pointer to the #i#th line + * bounding this #Quad#. + */ + TriaIterator > > + line (const unsigned int i) const; + + /** + * Return the #i#th child as a DoF quad + * iterator. This function is needed since + * the child function of the base + * class returns a quad accessor without + * access to the DoF data. + */ + TriaIterator > child (const unsigned int) const; + + /** + * Implement the copy operator needed + * for the iterator classes. + */ + void copy_from (const DoFQuadAccessor &a); +}; + + + + + +/** + * Intermediate, "typedef"-class, not for public use. + * + * Rationale for the declaration of members for this class: gcc 2.8 has a bug + * when deriving from explicitely specialized classes which materializes in + * the calculation of wrong addresses of member variables. By declaring the + * general template of #DoFSubstructAccessor# to have the same object layout as + * the specialized versions (using the same base classes), we fool the compiler, + * which still looks in the wrong place for the addresses but finds the + * right information. This way, at least ot works. + * + * Insert a guard, however, in the constructor to avoid that anyone (including + * the compiler) happens to use this class. + */ +template +class MGDoFSubstructAccessor : public MGDoFAccessor, + public TriaAccessor { + public: + DoFSubstructAccessor () { + Assert (false, ExcInternalError()); + }; + + DeclException0 (ExcInternalError); +}; + + + + +/** + * Intermediate, "typedef"-class, not for public use. + * + * \subsection{Rationale} + * + * This class is only a wrapper class used to do kind of a typedef + * 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> >; + * \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 + * the \Ref{DoFCellAccessor} dimension independent as an inheritance from + * #DoFSubstructAccessor#. If we had not declared these + * types, we would have to write two class declarations, one for + * #DoFCellAccessor<1>#, derived from #DoFLineAccessor<1,CellAccessor<1> ># + * and one for #DoFCellAccessor<2>#, derived from + * #DoFQuadAccessor<2,CellAccessor<2> >#. + */ +template <> +class MGDoFSubstructAccessor<1> : public MGDoFLineAccessor<1,CellAccessor<1> > { + public: + /** + * Constructor + */ + MGDoFSubstructAccessor (Triangulation<1> *tria, + const int level, + const int index, + const void *local_data) : + MGDoFLineAccessor<1,CellAccessor<1> > (tria,level,index,local_data) {}; + // do this here, since this way the + // CellAccessor has the possibility to know + // what a face_iterator is. Otherwise + // it would have to ask the MGDoFHandler + // which would need another big include + // file and maybe cyclic interdependence + typedef void * face_iterator; +}; + + + +/** + * Intermediate, "typedef"-class, not for public use. + * + * @see MGDoFSubstructAccessor<1> + */ +template <> +class MGDoFSubstructAccessor<2> : public MGDoFQuadAccessor<2,CellAccessor<2> > { + public: + /** + * Constructor + */ + MGDoFSubstructAccessor (Triangulation<2> *tria, + const int level, + const int index, + const void *local_data) : + MGDoFQuadAccessor<2,CellAccessor<2> > (tria,level,index,local_data) {}; + // do this here, since this way the + // CellAccessor has the possibility to know + // what a face_iterator is. Otherwise + // it would have to ask the DoFHandler + // which would need another big include + // file and maybe cyclic interdependence + typedef TriaIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > > face_iterator; +}; + + + + + + +/** + * Grant access to the degrees of freedom on a cell. In fact, since all + * access to the degrees of freedom has been enabled by the classes + * #DoFLineAccessor<1># and #DoFQuadAccessor<2># for the space dimension + * one and two, respectively, this class only collects the pieces + * together by deriving from the appropriate #DoF*Accessor# and the + * right #CellAccessor# and finally adding two functions which give + * access to the neighbors and children as #DoFCellAccessor# objects + * rather than #CellAccessor# objects (the latter function was inherited + * from the #CellAccessor# class). + * + * @author Wolfgang Bangerth, 1998 + */ +template +class MGDoFCellAccessor : public MGDoFSubstructAccessor { + public: + /** + * Constructor + */ + MGDoFCellAccessor (Triangulation *tria, + const int level, + const int index, + const void *local_data) : + MGDoFSubstructAccessor (tria,level,index,local_data) {}; + + /** + * Return the #i#th neighbor as a MGDoF cell + * iterator. This function is needed since + * the neighbor function of the base + * class returns a cell accessor without + * access to the MGDoF data. + */ + TriaIterator > neighbor (const unsigned int) const; + + /** + * Return the #i#th child as a MGDoF cell + * iterator. This function is needed since + * the child function of the base + * class returns a cell accessor without + * access to the DoF data. + */ + TriaIterator > child (const unsigned int) const; + + /** + * Return an iterator to the #i#th face + * of this cell. + * + * This function is not implemented in 1D, + * and maps to MGDoFQuadAccessor::line in 2D. + */ + typename MGDoFSubstructAccessor::face_iterator + face (const unsigned int i) const; + + /** + * Return the value of the given vector + * restricted to the dofs of this + * cell in the standard ordering: dofs + * on vertex 0, dofs on vertex 1, etc, + * dofs on line 0, dofs on line 1, etc, + * dofs on quad 0, etc. + * + * It is assumed that the vector already + * has the right size beforehand. The + * indices refer to the multilevel + * numbering local to the present + * level of this cell. The vector shall + * therefore have the same number of + * entries as there are degrees of + * freedom on this level. + */ + void get_mg_dof_values (const dVector &values, + dVector &dof_values) const; + + + /** + * Exception + */ + DeclException0 (ExcNotUsefulForThisDimension); +}; + + + + + +/*---------------------------- mg_dof_accessor.h ---------------------------*/ +/* end of #ifndef __mg_dof_accessor_H */ +#endif +/*---------------------------- mg_dof_accessor.h ---------------------------*/ diff --git a/deal.II/deal.II/include/dofs/mg_dof_handler.h b/deal.II/deal.II/include/dofs/mg_dof_handler.h new file mode 100644 index 0000000000..21ad59db20 --- /dev/null +++ b/deal.II/deal.II/include/dofs/mg_dof_handler.h @@ -0,0 +1,253 @@ +/*---------------------------- mg_dof.h ---------------------------*/ +/* $Id$ */ +#ifndef __mg_dof_H +#define __mg_dof_H +/*---------------------------- mg_dof.h ---------------------------*/ + + +#include + + + +// forward declarations +template class MGDoFLineAccessor; +template class MGDoFLineAccessor; +template class MGDoFQuadAccessor; +template class MGDoFQuadAccessor; + + + + + +template +class MGDoFHandler : public DoFHandler { + public: + /** + * Constructor. Take #tria# as the + * triangulation to work on. + */ + MGDoFHandler (Triangulation *tria); + + /** + * Destructor + */ + virtual ~MGDoFHandler (); + + /** + * Go through the triangulation and + * distribute the degrees of freedoms + * needed for the given finite element + * according to the given distribution + * method. We first call the #DoFHandler#'s + * function and then distribute the + * levelwise numbers. + * + * A copy of the transferred finite + * element is stored. + * + * This function uses the user flags of the + * triangulation object, so make sure you + * don't need them after calling this + * function, or if so store them. + */ + virtual void distribute_dofs (const FiniteElementBase &); + + /** + * Renumber the degrees of freedom according + * to the given scheme, eventually using + * constraint information and the given + * starting points. The starting points + * default to an empty list, the use of + * constraint information defaults to + * false. + * + * See the general documentation of the + * #DoFHandler# class for more details. + */ + virtual void renumber_dofs (const RenumberingMethod method, + const bool use_constraints = false, + const vector &starting_points = vector()); + + + private: + + /** We need for each vertex a list of the + * degree of freedom indices on each of the + * levels this vertex lives on. Since most + * vertices live only on a few levels, it + * is not economical to reserve indices for + * all the levels there are; rather, we + * create an object which holds the indices + * on those levels only where the vertex + * lives. To construct such an array, it + * is necessary to know beforehand which + * is the coarsest level the vertex lives + * on, how many levels it lives on and + * how many dofs there are on each vertex. + * If we have this information, we can + * allocate exactly the amount of memory + * which is needed and need not handle + * growing arrays and the like. + */ + class MGVertexDoFs { + public: + + /** + * Constructor. Allocates memory and + * sets all indices to #-1#. + */ + MGVertexDoFs (const unsigned int coarsest_level, + const unsigned int n_levels, + const unsigned int dofs_per_vertex); + + /** + * Destructor + */ + ~MGVertexDoFs (); + + /** + * Set the index with number + * #dof_number# of this vertex on + * #level# to the given index. To + * compute the position in the array, + * one has to specify how many + * dofs per vertex there are. It is + * not checked that the level number + * is below the number of the finest + * level this vertex lives on. + * + * The function is inline, so should + * be reasonably fast. + */ + void set_index (const unsigned int level, + const unsigned int dof_number, + const unsigned int dofs_per_vertex, + const unsigned int index); + + /** + * Return the index with number + * #dof_number# of this vertex on + * #level#. To + * compute the position in the array, + * one has to specify how many + * dofs per vertex there are. It is + * not checked that the level number + * is below the number of the finest + * level this vertex lives on. + * + * The function is inline, so should + * be reasonably fast. + */ + int get_index (const unsigned int level, + const unsigned int dof_number, + const unsigned int dofs_per_vertex) const; + + /** + * Exception. + */ + DeclException0 (ExcNoMemory); + /** + * Exception. + */ + DeclException0 (ExcInvalidIndex); + /** + * Exception. + */ + DeclException1 (ExcInvalidLevel, + int, + << "The given level number " << arg1 << " is below the " + << "coarsest level this vertex lives on."); + + private: + /** + * Store the coarsest level this + * vertex lives on. This is used + * as an offset when accessing the + * dofs of a specific level. + */ + const unsigned int coarsest_level; + + /** + * Array holding the indices. + */ + int *const indices; + }; + + + /** + * Reserve enough space for the MG dof + * indices for a given triangulation. + */ + void reserve_space (); + + /** + * Space to store the DoF numbers for the + * different levels. Unlike the #levels# + * object in the #DoFHandler#, these are + * not global numbers but rather are + * numbers which start from zero on each + * level. + */ + vector*> mg_levels; + + /** + * For each vertex there is a list of + * indices of the degrees of freedom indices + * on the different levels it lives on and + * which are these levels. + */ + vector mg_vertex_dofs; + + /** + * Vectors storing the number of degrees of + * freedom on each level. + */ + vector mg_used_dofs; + + friend class MGDoFLineAccessor >; + friend class MGDoFLineAccessor >; + friend class MGDoFQuadAccessor >; + friend class MGDoFQuadAccessor >; +}; + + + + + + + +/* ----------------------- Inline functions of MGVertexDoFs -------------------*/ + + +template +inline +void MGDoFHandler::MGVertexDoFs::set_index (const unsigned int level, + const unsigned int dof_number, + const unsigned int dofs_per_vertex, + const unsigned int index) { + Assert (level >= coarsest_level, ExcInvalidLevel(level)); + Assert (dof_number < dofs_per_vertex, ExcInvalidIndex ()); + + indices[(level-coarsest_level)*dofs_per_vertex + dof_number] = index; +}; + + + + +template +inline +int MGDoFHandler::MGVertexDoFs::get_index (const unsigned int level, + const unsigned int dof_number, + const unsigned int dofs_per_vertex) const { + Assert (level >= coarsest_level, ExcInvalidLevel(level)); + Assert (dof_number < dofs_per_vertex, ExcInvalidIndex ()); + + return indices[(level-coarsest_level)*dofs_per_vertex + dof_number]; +}; + + + +/*---------------------------- mg_dof.h ---------------------------*/ +/* end of #ifndef __mg_dof_H */ +#endif +/*---------------------------- mg_dof.h ---------------------------*/ diff --git a/deal.II/deal.II/include/multigrid/mg_dof_accessor.h b/deal.II/deal.II/include/multigrid/mg_dof_accessor.h new file mode 100644 index 0000000000..03274755dc --- /dev/null +++ b/deal.II/deal.II/include/multigrid/mg_dof_accessor.h @@ -0,0 +1,505 @@ +/*---------------------------- mg_dof_accessor.h ---------------------------*/ +/* $Id$ */ +#ifndef __mg_dof_accessor_H +#define __mg_dof_accessor_H +/*---------------------------- mg_dof_accessor.h ---------------------------*/ + + +#include + + +// forward declaration +template class MGDoFHandler; + + + + +/** + * Define the basis for accessors to the degrees of freedom for + * a multigrid DoF handler object. + * + * Note that it is allowed to construct an object of which the + * #mg_dof_handler# pointer is a Null pointer. Such an object would + * result in a strange kind of behaviour, though every reasonable + * operating system should disallow access through that pointer. + * The reason we do not check for the null pointer in the + * constructor which gets passed the #DoFHandler# pointer is that + * if we did we could not make dof iterators member of other classes + * (like in the #FEValues# class) if we did not know about the + * #DoFHandler# object to be used upon construction of that object. + * Through the way this class is implemented here, we allow the + * creation of a kind of virgin object which only gets useful if + * assigned to from another object before first usage. + * + * Opposite to construction, it is not possible to copy an object + * which has an invalid dof handler pointer. This is to guarantee + * that every iterator which is once assigned to is a valid + * object. However, this assertion only holds in debug mode, when + * the #Assert# macro is switched on. + * + * @author Wolfgang Bangerth, 1998 + */ +template +class MGDoFAccessor { + public: + /** + * Constructor + */ + MGDoFAccessor () : mg_dof_handler(0) { + Assert (false, ExcInvalidObject()); + }; + + /** + * This should be the default constructor. + */ + MGDoFAccessor (MGDoFHandler *mg_dof_handler) : + mg_dof_handler(mg_dof_handler) {}; + + /** + * Reset the DoF handler pointer. + */ + void set_mg_dof_handler (MGDoFHandler *dh) { + Assert (dh != 0, ExcInvalidObject()); + mg_dof_handler = dh; + }; + + /** + * Copy operator. + */ + MGDoFAccessor & operator = (const MGDoFAccessor &da) { + set_dof_handler (da.mg_dof_handler); + return *this; + }; + + protected: + /** + * Store the address of the #MGDoFHandler# object + * to be accessed. + */ + MGDoFHandler *mg_dof_handler; +}; + + + + + +/** + * Grant access to the degrees of freedom located on lines. + * This class follows mainly the route laid out by the accessor library + * declared in the triangulation library (\Ref{TriaAccessor}). It enables + * the user to access the degrees of freedom on the lines (there are similar + * versions for the DoFs on quads, etc), where the dimension of the underlying + * triangulation does not really matter (i.e. this accessor works with the + * lines in 1D-, 2D-, etc dimensions). + * + * + * \subsection{Usage} + * + * The \Ref{DoFDimensionInfo} classes inherited by the \Ref{DoFHandler} classes + * declare typedefs to iterators using the accessors declared in this class + * hierarchy tree. Usage is best to happens through these typedefs, since they + * are more secure to changes in the class naming and template interface as well + * as they provide easier typing (much less complicated names!). + * + * + * \subsection{Notes about the class hierarchy structure} + * + * The class hierarchy seems to be a bit confused here. The reason for this is + * that we would really like to derive a #DoFLineAccessor# from a #LineAccessor#. + * Unfortunately, we would run into problems, if we wanted a #DoFLineAccessor# + * in one spatial dimension, in which case a line is also a cell. The traditional + * solution would be to declare a #DoFCellAccessor<1># which is derived from + * #DoFLineAccessor<1># and #CellAccessor<1># (the #DoFLineAccessor# cannot + * itself be derived from #CellAccessor# since a line is not a cell + * unless in one space dimension), but since a #DoFLineAccessor# and a + * #CellAccessor# are both derived from #TriaAccessor#, we would have to make + * the last derivation virtual. + * + * Since we want to avoid virtual inheritance since this involves another + * indirection in every member variable access, we chose another way: we + * pass a second template parameter to a #DoFLineAccessor# which tells it + * which class to be derived from: if we are in one spatial dimension, the + * base class is to be #CellAccessor<1>#, in two or more dimensions it + * is a #LineAccessor#, i.e. am accessor to lines without the missing + * functionality needed for cells (neighbors, etc.). + * + * This way we can declare a #DoFCellAccessor# in one dimension by deriving + * from #DoFLineAccessor<1,CellAccessor<1> >#, thus getting the cell + * functionality through the #DoFLineAccessor# instead of through a virtual + * multiple inheritance of #DoFLineAccessor# and #CellAccessor<1>#. + * + * The same concept is used with #DoFQuadAccessor# classes etc. + * + * @author Wolfgang Bangerth, 1998 + */ +template +class MGDoFLineAccessor : public MGDoFAccessor, + public DoFLineAccessor { + public: + /** + * Default constructor, unused thus + * not implemented. + */ + MGDoFLineAccessor (); + + /** + * Constructor. The #local_data# + * argument is assumed to be a pointer + * to a #MGDoFHandler# object. + */ + MGDoFLineAccessor (Triangulation *tria, + const int level, + const int index, + const void *local_data); + + /** + * Return the index of the #i#th degree + * of freedom of this line on the level + * this line lives on. + */ + int mg_dof_index (const unsigned int i) const; + + /** + * Set the index of the #i#th degree + * of freedom of this line on the + * level this line lives on to #index#. + */ + void set_mg_dof_index (const unsigned int i, const int index) const; + + /** + * Return the index of the #i#th degree + * on the #vertex#th vertex for the + * level this line lives on. + */ + int mg_vertex_dof_index (const unsigned int vertex, + const unsigned int i) const; + + /** + * Set the index of the #i#th degree + * on the #vertex#th vertex to #index# + * for the level this line lives on. + */ + void set_mg_vertex_dof_index (const unsigned int vertex, + const unsigned int i, + const int index) const; + + /** + * Return the indices of the dofs of this + * line in the standard ordering: dofs + * on vertex 0, dofs on vertex 1, + * dofs on line. + * + * It is assumed that the vector already + * has the right size beforehand. The + * indices refer to the local numbering + * for the level this line lives on. + */ + void get_mg_dof_indices (vector &dof_indices) const; + + /** + * Return the #i#th child as a MGDoF line + * iterator. This function is needed since + * the child function of the base + * class returns a line accessor without + * access to the MG data. + */ + TriaIterator > child (const unsigned int) const; + + /** + * Implement the copy operator needed + * for the iterator classes. + */ + void copy_from (const MGDoFLineAccessor &a); +}; + + + + + + +/** + * Grant access to the multilevel degrees of freedom located on quads. + * + * @see DoFLineAccessor + */ +template +class MGDoFQuadAccessor : public MGDoFAccessor, + public DoFQuadAccessor { + public: + /** + * Default constructor, unused thus + * not implemented. + */ + MGDoFQuadAccessor (); + + /** + * Constructor. The #local_data# + * argument is assumed to be a pointer + * to a #DoFHandler# object. + */ + MGDoFQuadAccessor (Triangulation *tria, + const int level, + const int index, + const void *local_data); + + /** + * Return the index of the #i#th degree + * of freedom of this quad on the level + * this quad lives on. + */ + int mg_dof_index (const unsigned int i) const; + + /** + * Set the index of the #i#th degree + * of freedom of this quad on the + * level this quad lives on to #index#. + */ + void set_mg_dof_index (const unsigned int i, const int index) const; + + /** + * Return the index of the #i#th degree + * on the #vertex#th vertex for the level + * this quad lives on. + */ + int mg_vertex_dof_index (const unsigned int vertex, + const unsigned int i) const; + + /** + * Set the index of the #i#th degree + * on the #vertex#th vertex for the + * level this quad lives on to #index#. + */ + void set_mg_vertex_dof_index (const unsigned int vertex, + const unsigned int i, + const int index) const; + + /** + * Return the indices of the dofs of this + * quad in the standard ordering: dofs + * on vertex 0, dofs on vertex 1, etc, + * dofs on line 0, dofs on line 1, etc, + * dofs on quad 0, etc. + * + * It is assumed that the vector already + * has the right size beforehand. The + * indices refer to the local numbering + * for the level this quad lives on. + */ + void get_mg_dof_indices (vector &dof_indices) const; + + /** + * Return a pointer to the #i#th line + * bounding this #Quad#. + */ + TriaIterator > > + line (const unsigned int i) const; + + /** + * Return the #i#th child as a DoF quad + * iterator. This function is needed since + * the child function of the base + * class returns a quad accessor without + * access to the DoF data. + */ + TriaIterator > child (const unsigned int) const; + + /** + * Implement the copy operator needed + * for the iterator classes. + */ + void copy_from (const DoFQuadAccessor &a); +}; + + + + + +/** + * Intermediate, "typedef"-class, not for public use. + * + * Rationale for the declaration of members for this class: gcc 2.8 has a bug + * when deriving from explicitely specialized classes which materializes in + * the calculation of wrong addresses of member variables. By declaring the + * general template of #DoFSubstructAccessor# to have the same object layout as + * the specialized versions (using the same base classes), we fool the compiler, + * which still looks in the wrong place for the addresses but finds the + * right information. This way, at least ot works. + * + * Insert a guard, however, in the constructor to avoid that anyone (including + * the compiler) happens to use this class. + */ +template +class MGDoFSubstructAccessor : public MGDoFAccessor, + public TriaAccessor { + public: + DoFSubstructAccessor () { + Assert (false, ExcInternalError()); + }; + + DeclException0 (ExcInternalError); +}; + + + + +/** + * Intermediate, "typedef"-class, not for public use. + * + * \subsection{Rationale} + * + * This class is only a wrapper class used to do kind of a typedef + * 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> >; + * \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 + * the \Ref{DoFCellAccessor} dimension independent as an inheritance from + * #DoFSubstructAccessor#. If we had not declared these + * types, we would have to write two class declarations, one for + * #DoFCellAccessor<1>#, derived from #DoFLineAccessor<1,CellAccessor<1> ># + * and one for #DoFCellAccessor<2>#, derived from + * #DoFQuadAccessor<2,CellAccessor<2> >#. + */ +template <> +class MGDoFSubstructAccessor<1> : public MGDoFLineAccessor<1,CellAccessor<1> > { + public: + /** + * Constructor + */ + MGDoFSubstructAccessor (Triangulation<1> *tria, + const int level, + const int index, + const void *local_data) : + MGDoFLineAccessor<1,CellAccessor<1> > (tria,level,index,local_data) {}; + // do this here, since this way the + // CellAccessor has the possibility to know + // what a face_iterator is. Otherwise + // it would have to ask the MGDoFHandler + // which would need another big include + // file and maybe cyclic interdependence + typedef void * face_iterator; +}; + + + +/** + * Intermediate, "typedef"-class, not for public use. + * + * @see MGDoFSubstructAccessor<1> + */ +template <> +class MGDoFSubstructAccessor<2> : public MGDoFQuadAccessor<2,CellAccessor<2> > { + public: + /** + * Constructor + */ + MGDoFSubstructAccessor (Triangulation<2> *tria, + const int level, + const int index, + const void *local_data) : + MGDoFQuadAccessor<2,CellAccessor<2> > (tria,level,index,local_data) {}; + // do this here, since this way the + // CellAccessor has the possibility to know + // what a face_iterator is. Otherwise + // it would have to ask the DoFHandler + // which would need another big include + // file and maybe cyclic interdependence + typedef TriaIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > > face_iterator; +}; + + + + + + +/** + * Grant access to the degrees of freedom on a cell. In fact, since all + * access to the degrees of freedom has been enabled by the classes + * #DoFLineAccessor<1># and #DoFQuadAccessor<2># for the space dimension + * one and two, respectively, this class only collects the pieces + * together by deriving from the appropriate #DoF*Accessor# and the + * right #CellAccessor# and finally adding two functions which give + * access to the neighbors and children as #DoFCellAccessor# objects + * rather than #CellAccessor# objects (the latter function was inherited + * from the #CellAccessor# class). + * + * @author Wolfgang Bangerth, 1998 + */ +template +class MGDoFCellAccessor : public MGDoFSubstructAccessor { + public: + /** + * Constructor + */ + MGDoFCellAccessor (Triangulation *tria, + const int level, + const int index, + const void *local_data) : + MGDoFSubstructAccessor (tria,level,index,local_data) {}; + + /** + * Return the #i#th neighbor as a MGDoF cell + * iterator. This function is needed since + * the neighbor function of the base + * class returns a cell accessor without + * access to the MGDoF data. + */ + TriaIterator > neighbor (const unsigned int) const; + + /** + * Return the #i#th child as a MGDoF cell + * iterator. This function is needed since + * the child function of the base + * class returns a cell accessor without + * access to the DoF data. + */ + TriaIterator > child (const unsigned int) const; + + /** + * Return an iterator to the #i#th face + * of this cell. + * + * This function is not implemented in 1D, + * and maps to MGDoFQuadAccessor::line in 2D. + */ + typename MGDoFSubstructAccessor::face_iterator + face (const unsigned int i) const; + + /** + * Return the value of the given vector + * restricted to the dofs of this + * cell in the standard ordering: dofs + * on vertex 0, dofs on vertex 1, etc, + * dofs on line 0, dofs on line 1, etc, + * dofs on quad 0, etc. + * + * It is assumed that the vector already + * has the right size beforehand. The + * indices refer to the multilevel + * numbering local to the present + * level of this cell. The vector shall + * therefore have the same number of + * entries as there are degrees of + * freedom on this level. + */ + void get_mg_dof_values (const dVector &values, + dVector &dof_values) const; + + + /** + * Exception + */ + DeclException0 (ExcNotUsefulForThisDimension); +}; + + + + + +/*---------------------------- mg_dof_accessor.h ---------------------------*/ +/* end of #ifndef __mg_dof_accessor_H */ +#endif +/*---------------------------- mg_dof_accessor.h ---------------------------*/ diff --git a/deal.II/deal.II/include/multigrid/mg_dof_handler.h b/deal.II/deal.II/include/multigrid/mg_dof_handler.h new file mode 100644 index 0000000000..21ad59db20 --- /dev/null +++ b/deal.II/deal.II/include/multigrid/mg_dof_handler.h @@ -0,0 +1,253 @@ +/*---------------------------- mg_dof.h ---------------------------*/ +/* $Id$ */ +#ifndef __mg_dof_H +#define __mg_dof_H +/*---------------------------- mg_dof.h ---------------------------*/ + + +#include + + + +// forward declarations +template class MGDoFLineAccessor; +template class MGDoFLineAccessor; +template class MGDoFQuadAccessor; +template class MGDoFQuadAccessor; + + + + + +template +class MGDoFHandler : public DoFHandler { + public: + /** + * Constructor. Take #tria# as the + * triangulation to work on. + */ + MGDoFHandler (Triangulation *tria); + + /** + * Destructor + */ + virtual ~MGDoFHandler (); + + /** + * Go through the triangulation and + * distribute the degrees of freedoms + * needed for the given finite element + * according to the given distribution + * method. We first call the #DoFHandler#'s + * function and then distribute the + * levelwise numbers. + * + * A copy of the transferred finite + * element is stored. + * + * This function uses the user flags of the + * triangulation object, so make sure you + * don't need them after calling this + * function, or if so store them. + */ + virtual void distribute_dofs (const FiniteElementBase &); + + /** + * Renumber the degrees of freedom according + * to the given scheme, eventually using + * constraint information and the given + * starting points. The starting points + * default to an empty list, the use of + * constraint information defaults to + * false. + * + * See the general documentation of the + * #DoFHandler# class for more details. + */ + virtual void renumber_dofs (const RenumberingMethod method, + const bool use_constraints = false, + const vector &starting_points = vector()); + + + private: + + /** We need for each vertex a list of the + * degree of freedom indices on each of the + * levels this vertex lives on. Since most + * vertices live only on a few levels, it + * is not economical to reserve indices for + * all the levels there are; rather, we + * create an object which holds the indices + * on those levels only where the vertex + * lives. To construct such an array, it + * is necessary to know beforehand which + * is the coarsest level the vertex lives + * on, how many levels it lives on and + * how many dofs there are on each vertex. + * If we have this information, we can + * allocate exactly the amount of memory + * which is needed and need not handle + * growing arrays and the like. + */ + class MGVertexDoFs { + public: + + /** + * Constructor. Allocates memory and + * sets all indices to #-1#. + */ + MGVertexDoFs (const unsigned int coarsest_level, + const unsigned int n_levels, + const unsigned int dofs_per_vertex); + + /** + * Destructor + */ + ~MGVertexDoFs (); + + /** + * Set the index with number + * #dof_number# of this vertex on + * #level# to the given index. To + * compute the position in the array, + * one has to specify how many + * dofs per vertex there are. It is + * not checked that the level number + * is below the number of the finest + * level this vertex lives on. + * + * The function is inline, so should + * be reasonably fast. + */ + void set_index (const unsigned int level, + const unsigned int dof_number, + const unsigned int dofs_per_vertex, + const unsigned int index); + + /** + * Return the index with number + * #dof_number# of this vertex on + * #level#. To + * compute the position in the array, + * one has to specify how many + * dofs per vertex there are. It is + * not checked that the level number + * is below the number of the finest + * level this vertex lives on. + * + * The function is inline, so should + * be reasonably fast. + */ + int get_index (const unsigned int level, + const unsigned int dof_number, + const unsigned int dofs_per_vertex) const; + + /** + * Exception. + */ + DeclException0 (ExcNoMemory); + /** + * Exception. + */ + DeclException0 (ExcInvalidIndex); + /** + * Exception. + */ + DeclException1 (ExcInvalidLevel, + int, + << "The given level number " << arg1 << " is below the " + << "coarsest level this vertex lives on."); + + private: + /** + * Store the coarsest level this + * vertex lives on. This is used + * as an offset when accessing the + * dofs of a specific level. + */ + const unsigned int coarsest_level; + + /** + * Array holding the indices. + */ + int *const indices; + }; + + + /** + * Reserve enough space for the MG dof + * indices for a given triangulation. + */ + void reserve_space (); + + /** + * Space to store the DoF numbers for the + * different levels. Unlike the #levels# + * object in the #DoFHandler#, these are + * not global numbers but rather are + * numbers which start from zero on each + * level. + */ + vector*> mg_levels; + + /** + * For each vertex there is a list of + * indices of the degrees of freedom indices + * on the different levels it lives on and + * which are these levels. + */ + vector mg_vertex_dofs; + + /** + * Vectors storing the number of degrees of + * freedom on each level. + */ + vector mg_used_dofs; + + friend class MGDoFLineAccessor >; + friend class MGDoFLineAccessor >; + friend class MGDoFQuadAccessor >; + friend class MGDoFQuadAccessor >; +}; + + + + + + + +/* ----------------------- Inline functions of MGVertexDoFs -------------------*/ + + +template +inline +void MGDoFHandler::MGVertexDoFs::set_index (const unsigned int level, + const unsigned int dof_number, + const unsigned int dofs_per_vertex, + const unsigned int index) { + Assert (level >= coarsest_level, ExcInvalidLevel(level)); + Assert (dof_number < dofs_per_vertex, ExcInvalidIndex ()); + + indices[(level-coarsest_level)*dofs_per_vertex + dof_number] = index; +}; + + + + +template +inline +int MGDoFHandler::MGVertexDoFs::get_index (const unsigned int level, + const unsigned int dof_number, + const unsigned int dofs_per_vertex) const { + Assert (level >= coarsest_level, ExcInvalidLevel(level)); + Assert (dof_number < dofs_per_vertex, ExcInvalidIndex ()); + + return indices[(level-coarsest_level)*dofs_per_vertex + dof_number]; +}; + + + +/*---------------------------- mg_dof.h ---------------------------*/ +/* end of #ifndef __mg_dof_H */ +#endif +/*---------------------------- mg_dof.h ---------------------------*/ diff --git a/deal.II/deal.II/source/dofs/mg_dof_accessor.cc b/deal.II/deal.II/source/dofs/mg_dof_accessor.cc new file mode 100644 index 0000000000..322b395cfe --- /dev/null +++ b/deal.II/deal.II/source/dofs/mg_dof_accessor.cc @@ -0,0 +1,426 @@ +/* $Id$ */ + + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + + + + + + +/* ------------------------ MGDoFLineAccessor --------------------------- */ + +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))) {}; + + + +template +int MGDoFLineAccessor::mg_dof_index (const unsigned int i) const { + Assert (dof_handler != 0, ExcInvalidObject()); + Assert (mg_dof_handler != 0, ExcInvalidObject()); + // make sure a FE has been selected + // and enough room was reserved + Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (iget_selected_fe().dofs_per_line, + ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_line)); + + return mg_dof_handler->mg_levels[present_level] + ->line_dofs[present_index*dof_handler->get_selected_fe().dofs_per_line+i]; +}; + + + + + +template +void MGDoFLineAccessor::set_mg_dof_index (const unsigned int i, + const int index) const { + Assert (dof_handler != 0, ExcInvalidObject()); + Assert (mg_dof_handler != 0, ExcInvalidObject()); + // make sure a FE has been selected + // and enough room was reserved + Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (iget_selected_fe().dofs_per_line, + ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_line)); + + mg_dof_handler->mg_levels[present_level] + ->line_dofs[present_index*dof_handler->get_selected_fe().dofs_per_line+i] = index; +}; + + + + +template +int MGDoFLineAccessor::mg_vertex_dof_index (const unsigned int vertex, + const unsigned int i) const { + Assert (dof_handler != 0, ExcInvalidObject()); + Assert (mg_dof_handler != 0, ExcInvalidObject()); + Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (vertex<2, ExcInvalidIndex (i,0,2)); + Assert (iget_selected_fe().dofs_per_vertex, + ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_vertex)); + + return (mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)] + .get_index (present_level, i, dof_handler->get_selected_fe().dofs_per_vertex)); +}; + + + +template +void MGDoFLineAccessor::set_mg_vertex_dof_index (const unsigned int vertex, + const unsigned int i, + const int index) const { + Assert (dof_handler != 0, ExcInvalidObject()); + Assert (mg_dof_handler != 0, ExcInvalidObject()); + Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (vertex<2, ExcInvalidIndex (i,0,2)); + Assert (iget_selected_fe().dofs_per_vertex, + ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_vertex)); + + mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)] + .set_index (present_level, i, dof_handler->get_selected_fe().dofs_per_vertex, index); +}; + + + +template +void +MGDoFLineAccessor::get_mg_dof_indices (vector &dof_indices) const { + Assert (dof_handler != 0, ExcInvalidObject()); + Assert (mg_dof_handler != 0, ExcInvalidObject()); + Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (dof_indices.size() == (2*dof_handler->get_selected_fe().dofs_per_vertex + + dof_handler->get_selected_fe().dofs_per_line), + ExcVectorDoesNotMatch()); + + const unsigned int dofs_per_vertex = dof_handler->get_selected_fe().dofs_per_vertex, + dofs_per_line = dof_handler->get_selected_fe().dofs_per_line; + vector::iterator next = dof_indices.begin(); + for (unsigned int vertex=0; vertex<2; ++vertex) + for (unsigned int d=0; d +TriaIterator > +MGDoFLineAccessor::child (const unsigned int i) const { + TriaIterator > q (tria, + present_level+1, + child_index (i), + mg_dof_handler); + +#ifdef DEBUG + if (q.state() != past_the_end) + Assert (q->used(), typename TriaAccessor::ExcUnusedCellAsChild()); +#endif + return q; +}; + + + + +/* ------------------------ MGDoFQuadAccessor --------------------------- */ + +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))) {}; + + + +template +inline +int MGDoFQuadAccessor::mg_dof_index (const unsigned int i) const { + Assert (dof_handler != 0, ExcInvalidObject()); + Assert (mg_dof_handler != 0, ExcInvalidObject()); + // make sure a FE has been selected + // and enough room was reserved + Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (iget_selected_fe().dofs_per_quad, + ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_quad)); + + return mg_dof_handler->mg_levels[present_level] + ->quad_dofs[present_index*dof_handler->get_selected_fe().dofs_per_quad+i]; +}; + + + +template +void MGDoFQuadAccessor::set_mg_dof_index (const unsigned int i, + const int index) const { + Assert (dof_handler != 0, ExcInvalidObject()); + Assert (mg_dof_handler != 0, ExcInvalidObject()); + // make sure a FE has been selected + // and enough room was reserved + Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (iget_selected_fe().dofs_per_quad, + ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_quad)); + + mg_dof_handler->mg_levels[present_level] + ->quad_dofs[present_index*dof_handler->get_selected_fe().dofs_per_quad+i] = index; +}; + + + +template +inline +int MGDoFQuadAccessor::mg_vertex_dof_index (const unsigned int vertex, + const unsigned int i) const { + Assert (dof_handler != 0, ExcInvalidObject()); + Assert (mg_dof_handler != 0, ExcInvalidObject()); + Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (vertex<4, ExcInvalidIndex (i,0,4)); + Assert (iget_selected_fe().dofs_per_vertex, + ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_vertex)); + + return (mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)] + .get_index (present_level, i, dof_handler->get_selected_fe().dofs_per_vertex)); +}; + + + +template +void MGDoFQuadAccessor::set_mg_vertex_dof_index (const unsigned int vertex, + const unsigned int i, + const int index) const { + Assert (dof_handler != 0, ExcInvalidObject()); + Assert (mg_dof_handler != 0, ExcInvalidObject()); + Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (vertex<4, ExcInvalidIndex (i,0,4)); + Assert (iget_selected_fe().dofs_per_vertex, + ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_vertex)); + + mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)] + .set_index (present_level, i, dof_handler->get_selected_fe().dofs_per_vertex, index); +}; + + + +template +void +MGDoFQuadAccessor::get_mg_dof_indices (vector &dof_indices) const { + Assert (dof_handler != 0, ExcInvalidObject()); + Assert (mg_dof_handler != 0, ExcInvalidObject()); + Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (dof_indices.size() == (4*dof_handler->get_selected_fe().dofs_per_vertex + + 4*dof_handler->get_selected_fe().dofs_per_line + + dof_handler->get_selected_fe().dofs_per_quad), + ExcVectorDoesNotMatch()); + + const unsigned int dofs_per_vertex = dof_handler->get_selected_fe().dofs_per_vertex, + dofs_per_line = dof_handler->get_selected_fe().dofs_per_line, + dofs_per_quad = dof_handler->get_selected_fe().dofs_per_quad; + vector::iterator next = dof_indices.begin(); + for (unsigned int vertex=0; vertex<4; ++vertex) + for (unsigned int d=0; dline(line)->mg_dof_index(d); + for (unsigned int d=0; d +TriaIterator > > +MGDoFQuadAccessor::line (const unsigned int i) const { + Assert (i<4, ExcInvalidIndex (i, 0, 4)); + + return TriaIterator > > + ( + tria, + present_level, + line_index (i), + mg_dof_handler + ); +}; + + + +template +TriaIterator > +MGDoFQuadAccessor::child (const unsigned int i) const { + TriaIterator > q (tria, + present_level+1, + child_index (i), + mg_dof_handler); + +#ifdef DEBUG + if (q.state() != past_the_end) + Assert (q->used(), typename TriaAccessor::ExcUnusedCellAsChild()); +#endif + return q; +}; + + + + + +/*------------------------- Functions: MGDoFCellAccessor -----------------------*/ + +template +TriaIterator > +MGDoFCellAccessor::neighbor (const unsigned int i) const { + TriaIterator > q (tria, + neighbor_level (i), + neighbor_index (i), + mg_dof_handler); + +#ifdef DEBUG + if (q.state() != past_the_end) + Assert (q->used(), typename TriaAccessor::ExcUnusedCellAsNeighbor()); +#endif + return q; +}; + + + +template +TriaIterator > +MGDoFCellAccessor::child (const unsigned int i) const { + TriaIterator > q (tria, + present_level+1, + child_index (i), + mg_dof_handler); + +#ifdef DEBUG + if (q.state() != past_the_end) + Assert (q->used(), typename TriaAccessor::ExcUnusedCellAsChild()); +#endif + return q; +}; + + + +#if deal_II_dimension == 1 + +template <> +MGDoFSubstructAccessor<1>::face_iterator +MGDoFCellAccessor<1>::face (const unsigned int) const { + Assert (false, ExcNotUsefulForThisDimension()); + return 0; +}; + + + +template <> +void +MGDoFCellAccessor<1>::get_mg_dof_values (const dVector &values, + dVector &dof_values) const { + Assert (dof_handler != 0, ExcInvalidObject()); + Assert (mg_dof_handler != 0, ExcInvalidObject()); + Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (dof_values.size() == dof_handler->get_selected_fe().total_dofs, + ExcVectorDoesNotMatch()); + Assert (values.size() == dof_handler->n_dofs(), + ExcVectorDoesNotMatch()); + + const unsigned int dofs_per_vertex = dof_handler->get_selected_fe().dofs_per_vertex, + dofs_per_line = dof_handler->get_selected_fe().dofs_per_line; + vector::iterator next_dof_value=dof_values.begin(); + for (unsigned int vertex=0; vertex<2; ++vertex) + for (unsigned int d=0; d +MGDoFSubstructAccessor<2>::face_iterator +MGDoFCellAccessor<2>::face (const unsigned int i) const { + return line(i); +}; + + + +template <> +void +MGDoFCellAccessor<2>::get_mg_dof_values (const dVector &values, + dVector &dof_values) const { + Assert (dof_handler != 0, ExcInvalidObject()); + Assert (mg_dof_handler != 0, ExcInvalidObject()); + Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (dof_values.size() == dof_handler->get_selected_fe().total_dofs, + ExcVectorDoesNotMatch()); + Assert (values.size() == dof_handler->n_dofs(), + ExcVectorDoesNotMatch()); + + const unsigned int dofs_per_vertex = dof_handler->get_selected_fe().dofs_per_vertex, + dofs_per_line = dof_handler->get_selected_fe().dofs_per_line, + dofs_per_quad = dof_handler->get_selected_fe().dofs_per_quad; + vector::iterator next_dof_value=dof_values.begin(); + for (unsigned int vertex=0; vertex<4; ++vertex) + for (unsigned int d=0; dline(line)->mg_dof_index(d)); + for (unsigned int d=0; d >; +template class MGDoFCellAccessor<1>; + +template class TriaRawIterator<1,MGDoFCellAccessor<1> >; +template class TriaIterator<1,MGDoFCellAccessor<1> >; +template class TriaActiveIterator<1,MGDoFCellAccessor<1> >; +#endif + +#if deal_II_dimension == 2 +template class MGDoFLineAccessor<2,LineAccessor<2> >; +template class MGDoFQuadAccessor<2,QuadAccessor<2> >; +template class MGDoFQuadAccessor<2,CellAccessor<2> >; +template class MGDoFCellAccessor<2>; + +template class TriaRawIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > >; +template class TriaRawIterator<2,MGDoFCellAccessor<2> >; +template class TriaIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > >; +template class TriaIterator<2,MGDoFCellAccessor<2> >; +template class TriaActiveIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > >; +template class TriaActiveIterator<2,MGDoFCellAccessor<2> >; +#endif + + diff --git a/deal.II/deal.II/source/dofs/mg_dof_handler.cc b/deal.II/deal.II/source/dofs/mg_dof_handler.cc new file mode 100644 index 0000000000..14b9a55a01 --- /dev/null +++ b/deal.II/deal.II/source/dofs/mg_dof_handler.cc @@ -0,0 +1,92 @@ +/* $Id$ */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + + +/* ------------------------ 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]) +{ + Assert (indices != 0, ExcNoMemory ()); + + for (unsigned int i=0; i +MGDoFHandler::MGVertexDoFs::~MGVertexDoFs () { + delete[] indices; +}; + + + + + + + +/* ------------------------ MGDoFHandler ------------------------------------- */ + +template +MGDoFHandler::MGDoFHandler (Triangulation *tria) : + DoFHandler (tria) +{}; + + + +template +MGDoFHandler::~MGDoFHandler () {}; + + + + +template +void MGDoFHandler::distribute_dofs (const FiniteElementBase &fe) { + // first distribute global dofs + DoFHandler::distribute_dofs (fe); + + + // reserve space for the MG dof numbers + reserve_space (); + +}; + + + + +#if deal_II_dimension == 1 + +template <> +void MGDoFHandler<1>::reserve_space () { +}; + +#endif + + +#if deal_II_dimension == 2 + +template <> +void MGDoFHandler<2>::reserve_space () { +}; + +#endif + + +// explicite instantiations +template class MGDoFHandler; diff --git a/deal.II/deal.II/source/multigrid/mg_dof_accessor.cc b/deal.II/deal.II/source/multigrid/mg_dof_accessor.cc new file mode 100644 index 0000000000..322b395cfe --- /dev/null +++ b/deal.II/deal.II/source/multigrid/mg_dof_accessor.cc @@ -0,0 +1,426 @@ +/* $Id$ */ + + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + + + + + + +/* ------------------------ MGDoFLineAccessor --------------------------- */ + +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))) {}; + + + +template +int MGDoFLineAccessor::mg_dof_index (const unsigned int i) const { + Assert (dof_handler != 0, ExcInvalidObject()); + Assert (mg_dof_handler != 0, ExcInvalidObject()); + // make sure a FE has been selected + // and enough room was reserved + Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (iget_selected_fe().dofs_per_line, + ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_line)); + + return mg_dof_handler->mg_levels[present_level] + ->line_dofs[present_index*dof_handler->get_selected_fe().dofs_per_line+i]; +}; + + + + + +template +void MGDoFLineAccessor::set_mg_dof_index (const unsigned int i, + const int index) const { + Assert (dof_handler != 0, ExcInvalidObject()); + Assert (mg_dof_handler != 0, ExcInvalidObject()); + // make sure a FE has been selected + // and enough room was reserved + Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (iget_selected_fe().dofs_per_line, + ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_line)); + + mg_dof_handler->mg_levels[present_level] + ->line_dofs[present_index*dof_handler->get_selected_fe().dofs_per_line+i] = index; +}; + + + + +template +int MGDoFLineAccessor::mg_vertex_dof_index (const unsigned int vertex, + const unsigned int i) const { + Assert (dof_handler != 0, ExcInvalidObject()); + Assert (mg_dof_handler != 0, ExcInvalidObject()); + Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (vertex<2, ExcInvalidIndex (i,0,2)); + Assert (iget_selected_fe().dofs_per_vertex, + ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_vertex)); + + return (mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)] + .get_index (present_level, i, dof_handler->get_selected_fe().dofs_per_vertex)); +}; + + + +template +void MGDoFLineAccessor::set_mg_vertex_dof_index (const unsigned int vertex, + const unsigned int i, + const int index) const { + Assert (dof_handler != 0, ExcInvalidObject()); + Assert (mg_dof_handler != 0, ExcInvalidObject()); + Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (vertex<2, ExcInvalidIndex (i,0,2)); + Assert (iget_selected_fe().dofs_per_vertex, + ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_vertex)); + + mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)] + .set_index (present_level, i, dof_handler->get_selected_fe().dofs_per_vertex, index); +}; + + + +template +void +MGDoFLineAccessor::get_mg_dof_indices (vector &dof_indices) const { + Assert (dof_handler != 0, ExcInvalidObject()); + Assert (mg_dof_handler != 0, ExcInvalidObject()); + Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (dof_indices.size() == (2*dof_handler->get_selected_fe().dofs_per_vertex + + dof_handler->get_selected_fe().dofs_per_line), + ExcVectorDoesNotMatch()); + + const unsigned int dofs_per_vertex = dof_handler->get_selected_fe().dofs_per_vertex, + dofs_per_line = dof_handler->get_selected_fe().dofs_per_line; + vector::iterator next = dof_indices.begin(); + for (unsigned int vertex=0; vertex<2; ++vertex) + for (unsigned int d=0; d +TriaIterator > +MGDoFLineAccessor::child (const unsigned int i) const { + TriaIterator > q (tria, + present_level+1, + child_index (i), + mg_dof_handler); + +#ifdef DEBUG + if (q.state() != past_the_end) + Assert (q->used(), typename TriaAccessor::ExcUnusedCellAsChild()); +#endif + return q; +}; + + + + +/* ------------------------ MGDoFQuadAccessor --------------------------- */ + +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))) {}; + + + +template +inline +int MGDoFQuadAccessor::mg_dof_index (const unsigned int i) const { + Assert (dof_handler != 0, ExcInvalidObject()); + Assert (mg_dof_handler != 0, ExcInvalidObject()); + // make sure a FE has been selected + // and enough room was reserved + Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (iget_selected_fe().dofs_per_quad, + ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_quad)); + + return mg_dof_handler->mg_levels[present_level] + ->quad_dofs[present_index*dof_handler->get_selected_fe().dofs_per_quad+i]; +}; + + + +template +void MGDoFQuadAccessor::set_mg_dof_index (const unsigned int i, + const int index) const { + Assert (dof_handler != 0, ExcInvalidObject()); + Assert (mg_dof_handler != 0, ExcInvalidObject()); + // make sure a FE has been selected + // and enough room was reserved + Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (iget_selected_fe().dofs_per_quad, + ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_quad)); + + mg_dof_handler->mg_levels[present_level] + ->quad_dofs[present_index*dof_handler->get_selected_fe().dofs_per_quad+i] = index; +}; + + + +template +inline +int MGDoFQuadAccessor::mg_vertex_dof_index (const unsigned int vertex, + const unsigned int i) const { + Assert (dof_handler != 0, ExcInvalidObject()); + Assert (mg_dof_handler != 0, ExcInvalidObject()); + Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (vertex<4, ExcInvalidIndex (i,0,4)); + Assert (iget_selected_fe().dofs_per_vertex, + ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_vertex)); + + return (mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)] + .get_index (present_level, i, dof_handler->get_selected_fe().dofs_per_vertex)); +}; + + + +template +void MGDoFQuadAccessor::set_mg_vertex_dof_index (const unsigned int vertex, + const unsigned int i, + const int index) const { + Assert (dof_handler != 0, ExcInvalidObject()); + Assert (mg_dof_handler != 0, ExcInvalidObject()); + Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (vertex<4, ExcInvalidIndex (i,0,4)); + Assert (iget_selected_fe().dofs_per_vertex, + ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_vertex)); + + mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)] + .set_index (present_level, i, dof_handler->get_selected_fe().dofs_per_vertex, index); +}; + + + +template +void +MGDoFQuadAccessor::get_mg_dof_indices (vector &dof_indices) const { + Assert (dof_handler != 0, ExcInvalidObject()); + Assert (mg_dof_handler != 0, ExcInvalidObject()); + Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (dof_indices.size() == (4*dof_handler->get_selected_fe().dofs_per_vertex + + 4*dof_handler->get_selected_fe().dofs_per_line + + dof_handler->get_selected_fe().dofs_per_quad), + ExcVectorDoesNotMatch()); + + const unsigned int dofs_per_vertex = dof_handler->get_selected_fe().dofs_per_vertex, + dofs_per_line = dof_handler->get_selected_fe().dofs_per_line, + dofs_per_quad = dof_handler->get_selected_fe().dofs_per_quad; + vector::iterator next = dof_indices.begin(); + for (unsigned int vertex=0; vertex<4; ++vertex) + for (unsigned int d=0; dline(line)->mg_dof_index(d); + for (unsigned int d=0; d +TriaIterator > > +MGDoFQuadAccessor::line (const unsigned int i) const { + Assert (i<4, ExcInvalidIndex (i, 0, 4)); + + return TriaIterator > > + ( + tria, + present_level, + line_index (i), + mg_dof_handler + ); +}; + + + +template +TriaIterator > +MGDoFQuadAccessor::child (const unsigned int i) const { + TriaIterator > q (tria, + present_level+1, + child_index (i), + mg_dof_handler); + +#ifdef DEBUG + if (q.state() != past_the_end) + Assert (q->used(), typename TriaAccessor::ExcUnusedCellAsChild()); +#endif + return q; +}; + + + + + +/*------------------------- Functions: MGDoFCellAccessor -----------------------*/ + +template +TriaIterator > +MGDoFCellAccessor::neighbor (const unsigned int i) const { + TriaIterator > q (tria, + neighbor_level (i), + neighbor_index (i), + mg_dof_handler); + +#ifdef DEBUG + if (q.state() != past_the_end) + Assert (q->used(), typename TriaAccessor::ExcUnusedCellAsNeighbor()); +#endif + return q; +}; + + + +template +TriaIterator > +MGDoFCellAccessor::child (const unsigned int i) const { + TriaIterator > q (tria, + present_level+1, + child_index (i), + mg_dof_handler); + +#ifdef DEBUG + if (q.state() != past_the_end) + Assert (q->used(), typename TriaAccessor::ExcUnusedCellAsChild()); +#endif + return q; +}; + + + +#if deal_II_dimension == 1 + +template <> +MGDoFSubstructAccessor<1>::face_iterator +MGDoFCellAccessor<1>::face (const unsigned int) const { + Assert (false, ExcNotUsefulForThisDimension()); + return 0; +}; + + + +template <> +void +MGDoFCellAccessor<1>::get_mg_dof_values (const dVector &values, + dVector &dof_values) const { + Assert (dof_handler != 0, ExcInvalidObject()); + Assert (mg_dof_handler != 0, ExcInvalidObject()); + Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (dof_values.size() == dof_handler->get_selected_fe().total_dofs, + ExcVectorDoesNotMatch()); + Assert (values.size() == dof_handler->n_dofs(), + ExcVectorDoesNotMatch()); + + const unsigned int dofs_per_vertex = dof_handler->get_selected_fe().dofs_per_vertex, + dofs_per_line = dof_handler->get_selected_fe().dofs_per_line; + vector::iterator next_dof_value=dof_values.begin(); + for (unsigned int vertex=0; vertex<2; ++vertex) + for (unsigned int d=0; d +MGDoFSubstructAccessor<2>::face_iterator +MGDoFCellAccessor<2>::face (const unsigned int i) const { + return line(i); +}; + + + +template <> +void +MGDoFCellAccessor<2>::get_mg_dof_values (const dVector &values, + dVector &dof_values) const { + Assert (dof_handler != 0, ExcInvalidObject()); + Assert (mg_dof_handler != 0, ExcInvalidObject()); + Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject()); + Assert (dof_values.size() == dof_handler->get_selected_fe().total_dofs, + ExcVectorDoesNotMatch()); + Assert (values.size() == dof_handler->n_dofs(), + ExcVectorDoesNotMatch()); + + const unsigned int dofs_per_vertex = dof_handler->get_selected_fe().dofs_per_vertex, + dofs_per_line = dof_handler->get_selected_fe().dofs_per_line, + dofs_per_quad = dof_handler->get_selected_fe().dofs_per_quad; + vector::iterator next_dof_value=dof_values.begin(); + for (unsigned int vertex=0; vertex<4; ++vertex) + for (unsigned int d=0; dline(line)->mg_dof_index(d)); + for (unsigned int d=0; d >; +template class MGDoFCellAccessor<1>; + +template class TriaRawIterator<1,MGDoFCellAccessor<1> >; +template class TriaIterator<1,MGDoFCellAccessor<1> >; +template class TriaActiveIterator<1,MGDoFCellAccessor<1> >; +#endif + +#if deal_II_dimension == 2 +template class MGDoFLineAccessor<2,LineAccessor<2> >; +template class MGDoFQuadAccessor<2,QuadAccessor<2> >; +template class MGDoFQuadAccessor<2,CellAccessor<2> >; +template class MGDoFCellAccessor<2>; + +template class TriaRawIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > >; +template class TriaRawIterator<2,MGDoFCellAccessor<2> >; +template class TriaIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > >; +template class TriaIterator<2,MGDoFCellAccessor<2> >; +template class TriaActiveIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > >; +template class TriaActiveIterator<2,MGDoFCellAccessor<2> >; +#endif + + diff --git a/deal.II/deal.II/source/multigrid/mg_dof_handler.cc b/deal.II/deal.II/source/multigrid/mg_dof_handler.cc new file mode 100644 index 0000000000..14b9a55a01 --- /dev/null +++ b/deal.II/deal.II/source/multigrid/mg_dof_handler.cc @@ -0,0 +1,92 @@ +/* $Id$ */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + + +/* ------------------------ 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]) +{ + Assert (indices != 0, ExcNoMemory ()); + + for (unsigned int i=0; i +MGDoFHandler::MGVertexDoFs::~MGVertexDoFs () { + delete[] indices; +}; + + + + + + + +/* ------------------------ MGDoFHandler ------------------------------------- */ + +template +MGDoFHandler::MGDoFHandler (Triangulation *tria) : + DoFHandler (tria) +{}; + + + +template +MGDoFHandler::~MGDoFHandler () {}; + + + + +template +void MGDoFHandler::distribute_dofs (const FiniteElementBase &fe) { + // first distribute global dofs + DoFHandler::distribute_dofs (fe); + + + // reserve space for the MG dof numbers + reserve_space (); + +}; + + + + +#if deal_II_dimension == 1 + +template <> +void MGDoFHandler<1>::reserve_space () { +}; + +#endif + + +#if deal_II_dimension == 2 + +template <> +void MGDoFHandler<2>::reserve_space () { +}; + +#endif + + +// explicite instantiations +template class MGDoFHandler; -- 2.39.5