--- /dev/null
+/*---------------------------- mg_dof_accessor.h ---------------------------*/
+/* $Id$ */
+#ifndef __mg_dof_accessor_H
+#define __mg_dof_accessor_H
+/*---------------------------- mg_dof_accessor.h ---------------------------*/
+
+
+#include <grid/dof_accessor.h>
+
+
+// forward declaration
+template <int dim> 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 <int dim>
+class MGDoFAccessor {
+ public:
+ /**
+ * Constructor
+ */
+ MGDoFAccessor () : mg_dof_handler(0) {
+ Assert (false, ExcInvalidObject());
+ };
+
+ /**
+ * This should be the default constructor.
+ */
+ MGDoFAccessor (MGDoFHandler<dim> *mg_dof_handler) :
+ mg_dof_handler(mg_dof_handler) {};
+
+ /**
+ * Reset the DoF handler pointer.
+ */
+ void set_mg_dof_handler (MGDoFHandler<dim> *dh) {
+ Assert (dh != 0, ExcInvalidObject());
+ mg_dof_handler = dh;
+ };
+
+ /**
+ * Copy operator.
+ */
+ MGDoFAccessor<dim> & operator = (const MGDoFAccessor<dim> &da) {
+ set_dof_handler (da.mg_dof_handler);
+ return *this;
+ };
+
+ protected:
+ /**
+ * Store the address of the #MGDoFHandler# object
+ * to be accessed.
+ */
+ MGDoFHandler<dim> *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<dim># cannot
+ * itself be derived from #CellAccessor<dim># 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<dim>#, 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 <int dim, typename BaseClass>
+class MGDoFLineAccessor : public MGDoFAccessor<dim>,
+ public DoFLineAccessor<dim, BaseClass> {
+ public:
+ /**
+ * Default constructor, unused thus
+ * not implemented.
+ */
+ MGDoFLineAccessor ();
+
+ /**
+ * Constructor. The #local_data#
+ * argument is assumed to be a pointer
+ * to a #MGDoFHandler<dim># object.
+ */
+ MGDoFLineAccessor (Triangulation<dim> *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<int> &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<dim,MGDoFLineAccessor<dim,BaseClass> > child (const unsigned int) const;
+
+ /**
+ * Implement the copy operator needed
+ * for the iterator classes.
+ */
+ void copy_from (const MGDoFLineAccessor<dim,BaseClass> &a);
+};
+
+
+
+
+
+
+/**
+ * Grant access to the multilevel degrees of freedom located on quads.
+ *
+ * @see DoFLineAccessor
+ */
+template <int dim, typename BaseClass>
+class MGDoFQuadAccessor : public MGDoFAccessor<dim>,
+ public DoFQuadAccessor<dim, BaseClass> {
+ public:
+ /**
+ * Default constructor, unused thus
+ * not implemented.
+ */
+ MGDoFQuadAccessor ();
+
+ /**
+ * Constructor. The #local_data#
+ * argument is assumed to be a pointer
+ * to a #DoFHandler<dim># object.
+ */
+ MGDoFQuadAccessor (Triangulation<dim> *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<int> &dof_indices) const;
+
+ /**
+ * Return a pointer to the #i#th line
+ * bounding this #Quad#.
+ */
+ TriaIterator<dim,MGDoFLineAccessor<dim,LineAccessor<dim> > >
+ 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<dim,MGDoFQuadAccessor<dim,BaseClass> > child (const unsigned int) const;
+
+ /**
+ * Implement the copy operator needed
+ * for the iterator classes.
+ */
+ void copy_from (const DoFQuadAccessor<dim,BaseClass> &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 <int dim>
+class MGDoFSubstructAccessor : public MGDoFAccessor<dim>,
+ public TriaAccessor<dim> {
+ 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<dim>#. 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<dim>
+ // 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<dim>
+ // 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<dim># 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<dim># class).
+ *
+ * @author Wolfgang Bangerth, 1998
+ */
+template <int dim>
+class MGDoFCellAccessor : public MGDoFSubstructAccessor<dim> {
+ public:
+ /**
+ * Constructor
+ */
+ MGDoFCellAccessor (Triangulation<dim> *tria,
+ const int level,
+ const int index,
+ const void *local_data) :
+ MGDoFSubstructAccessor<dim> (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<dim,MGDoFCellAccessor<dim> > 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<dim,MGDoFCellAccessor<dim> > 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<dim>::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 ---------------------------*/
--- /dev/null
+/*---------------------------- mg_dof.h ---------------------------*/
+/* $Id$ */
+#ifndef __mg_dof_H
+#define __mg_dof_H
+/*---------------------------- mg_dof.h ---------------------------*/
+
+
+#include <grid/dof.h>
+
+
+
+// forward declarations
+template <int dim, typename BaseClass> class MGDoFLineAccessor;
+template <int dim, typename BaseClass> class MGDoFLineAccessor;
+template <int dim, typename BaseClass> class MGDoFQuadAccessor;
+template <int dim, typename BaseClass> class MGDoFQuadAccessor;
+
+
+
+
+
+template <int dim>
+class MGDoFHandler : public DoFHandler<dim> {
+ public:
+ /**
+ * Constructor. Take #tria# as the
+ * triangulation to work on.
+ */
+ MGDoFHandler (Triangulation<dim> *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<dim> &);
+
+ /**
+ * 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<int> &starting_points = vector<int>());
+
+
+ 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<DoFLevel<dim>*> 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<MGVertexDoFs> mg_vertex_dofs;
+
+ /**
+ * Vectors storing the number of degrees of
+ * freedom on each level.
+ */
+ vector<unsigned int> mg_used_dofs;
+
+ friend class MGDoFLineAccessor<dim, LineAccessor<dim> >;
+ friend class MGDoFLineAccessor<dim, CellAccessor<dim> >;
+ friend class MGDoFQuadAccessor<dim, QuadAccessor<dim> >;
+ friend class MGDoFQuadAccessor<dim, CellAccessor<dim> >;
+};
+
+
+
+
+
+
+
+/* ----------------------- Inline functions of MGVertexDoFs -------------------*/
+
+
+template <int dim>
+inline
+void MGDoFHandler<dim>::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 <int dim>
+inline
+int MGDoFHandler<dim>::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 ---------------------------*/
--- /dev/null
+/*---------------------------- mg_dof_accessor.h ---------------------------*/
+/* $Id$ */
+#ifndef __mg_dof_accessor_H
+#define __mg_dof_accessor_H
+/*---------------------------- mg_dof_accessor.h ---------------------------*/
+
+
+#include <grid/dof_accessor.h>
+
+
+// forward declaration
+template <int dim> 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 <int dim>
+class MGDoFAccessor {
+ public:
+ /**
+ * Constructor
+ */
+ MGDoFAccessor () : mg_dof_handler(0) {
+ Assert (false, ExcInvalidObject());
+ };
+
+ /**
+ * This should be the default constructor.
+ */
+ MGDoFAccessor (MGDoFHandler<dim> *mg_dof_handler) :
+ mg_dof_handler(mg_dof_handler) {};
+
+ /**
+ * Reset the DoF handler pointer.
+ */
+ void set_mg_dof_handler (MGDoFHandler<dim> *dh) {
+ Assert (dh != 0, ExcInvalidObject());
+ mg_dof_handler = dh;
+ };
+
+ /**
+ * Copy operator.
+ */
+ MGDoFAccessor<dim> & operator = (const MGDoFAccessor<dim> &da) {
+ set_dof_handler (da.mg_dof_handler);
+ return *this;
+ };
+
+ protected:
+ /**
+ * Store the address of the #MGDoFHandler# object
+ * to be accessed.
+ */
+ MGDoFHandler<dim> *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<dim># cannot
+ * itself be derived from #CellAccessor<dim># 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<dim>#, 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 <int dim, typename BaseClass>
+class MGDoFLineAccessor : public MGDoFAccessor<dim>,
+ public DoFLineAccessor<dim, BaseClass> {
+ public:
+ /**
+ * Default constructor, unused thus
+ * not implemented.
+ */
+ MGDoFLineAccessor ();
+
+ /**
+ * Constructor. The #local_data#
+ * argument is assumed to be a pointer
+ * to a #MGDoFHandler<dim># object.
+ */
+ MGDoFLineAccessor (Triangulation<dim> *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<int> &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<dim,MGDoFLineAccessor<dim,BaseClass> > child (const unsigned int) const;
+
+ /**
+ * Implement the copy operator needed
+ * for the iterator classes.
+ */
+ void copy_from (const MGDoFLineAccessor<dim,BaseClass> &a);
+};
+
+
+
+
+
+
+/**
+ * Grant access to the multilevel degrees of freedom located on quads.
+ *
+ * @see DoFLineAccessor
+ */
+template <int dim, typename BaseClass>
+class MGDoFQuadAccessor : public MGDoFAccessor<dim>,
+ public DoFQuadAccessor<dim, BaseClass> {
+ public:
+ /**
+ * Default constructor, unused thus
+ * not implemented.
+ */
+ MGDoFQuadAccessor ();
+
+ /**
+ * Constructor. The #local_data#
+ * argument is assumed to be a pointer
+ * to a #DoFHandler<dim># object.
+ */
+ MGDoFQuadAccessor (Triangulation<dim> *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<int> &dof_indices) const;
+
+ /**
+ * Return a pointer to the #i#th line
+ * bounding this #Quad#.
+ */
+ TriaIterator<dim,MGDoFLineAccessor<dim,LineAccessor<dim> > >
+ 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<dim,MGDoFQuadAccessor<dim,BaseClass> > child (const unsigned int) const;
+
+ /**
+ * Implement the copy operator needed
+ * for the iterator classes.
+ */
+ void copy_from (const DoFQuadAccessor<dim,BaseClass> &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 <int dim>
+class MGDoFSubstructAccessor : public MGDoFAccessor<dim>,
+ public TriaAccessor<dim> {
+ 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<dim>#. 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<dim>
+ // 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<dim>
+ // 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<dim># 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<dim># class).
+ *
+ * @author Wolfgang Bangerth, 1998
+ */
+template <int dim>
+class MGDoFCellAccessor : public MGDoFSubstructAccessor<dim> {
+ public:
+ /**
+ * Constructor
+ */
+ MGDoFCellAccessor (Triangulation<dim> *tria,
+ const int level,
+ const int index,
+ const void *local_data) :
+ MGDoFSubstructAccessor<dim> (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<dim,MGDoFCellAccessor<dim> > 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<dim,MGDoFCellAccessor<dim> > 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<dim>::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 ---------------------------*/
--- /dev/null
+/*---------------------------- mg_dof.h ---------------------------*/
+/* $Id$ */
+#ifndef __mg_dof_H
+#define __mg_dof_H
+/*---------------------------- mg_dof.h ---------------------------*/
+
+
+#include <grid/dof.h>
+
+
+
+// forward declarations
+template <int dim, typename BaseClass> class MGDoFLineAccessor;
+template <int dim, typename BaseClass> class MGDoFLineAccessor;
+template <int dim, typename BaseClass> class MGDoFQuadAccessor;
+template <int dim, typename BaseClass> class MGDoFQuadAccessor;
+
+
+
+
+
+template <int dim>
+class MGDoFHandler : public DoFHandler<dim> {
+ public:
+ /**
+ * Constructor. Take #tria# as the
+ * triangulation to work on.
+ */
+ MGDoFHandler (Triangulation<dim> *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<dim> &);
+
+ /**
+ * 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<int> &starting_points = vector<int>());
+
+
+ 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<DoFLevel<dim>*> 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<MGVertexDoFs> mg_vertex_dofs;
+
+ /**
+ * Vectors storing the number of degrees of
+ * freedom on each level.
+ */
+ vector<unsigned int> mg_used_dofs;
+
+ friend class MGDoFLineAccessor<dim, LineAccessor<dim> >;
+ friend class MGDoFLineAccessor<dim, CellAccessor<dim> >;
+ friend class MGDoFQuadAccessor<dim, QuadAccessor<dim> >;
+ friend class MGDoFQuadAccessor<dim, CellAccessor<dim> >;
+};
+
+
+
+
+
+
+
+/* ----------------------- Inline functions of MGVertexDoFs -------------------*/
+
+
+template <int dim>
+inline
+void MGDoFHandler<dim>::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 <int dim>
+inline
+int MGDoFHandler<dim>::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 ---------------------------*/
--- /dev/null
+/* $Id$ */
+
+
+#include <grid/mg_dof_accessor.h>
+#include <grid/mg_dof.h>
+#include <grid/dof_levels.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_iterator.templates.h>
+#include <fe/fe.h>
+
+#include <lac/dvector.h>
+#include <lac/dfmatrix.h>
+#include <lac/dsmatrix.h>
+
+
+
+
+
+
+/* ------------------------ MGDoFLineAccessor --------------------------- */
+
+template <int dim, typename BaseClass>
+MGDoFLineAccessor<dim,BaseClass>::MGDoFLineAccessor (Triangulation<dim> *tria,
+ const int level,
+ const int index,
+ const void *local_data) :
+ MGDoFAccessor<dim> ((MGDoFHandler<dim>*)local_data),
+ DoFLineAccessor(tria,level,index,
+ static_cast<DoFHandler<dim>*>
+ (reinterpret_cast<MGDoFHandler<dim>*>(local_data))) {};
+
+
+
+template <int dim, typename BaseClass>
+int MGDoFLineAccessor<dim,BaseClass>::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 (i<dof_handler->get_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 <int dim, typename BaseClass>
+void MGDoFLineAccessor<dim,BaseClass>::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 (i<dof_handler->get_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 dim, typename BaseClass>
+int MGDoFLineAccessor<dim,BaseClass>::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 (i<dof_handler->get_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 <int dim, typename BaseClass>
+void MGDoFLineAccessor<dim,BaseClass>::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 (i<dof_handler->get_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 <int dim, typename BaseClass>
+void
+MGDoFLineAccessor<dim,BaseClass>::get_mg_dof_indices (vector<int> &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<int>::iterator next = dof_indices.begin();
+ for (unsigned int vertex=0; vertex<2; ++vertex)
+ for (unsigned int d=0; d<dofs_per_vertex; ++d)
+ *next++ = mg_vertex_dof_index(vertex,d);
+ for (unsigned int d=0; d<dofs_per_line; ++d)
+ *next++ = mg_dof_index(d);
+};
+
+
+
+template <int dim, typename BaseClass>
+TriaIterator<dim,MGDoFLineAccessor<dim,BaseClass> >
+MGDoFLineAccessor<dim,BaseClass>::child (const unsigned int i) const {
+ TriaIterator<dim,MGDoFLineAccessor<dim,BaseClass> > q (tria,
+ present_level+1,
+ child_index (i),
+ mg_dof_handler);
+
+#ifdef DEBUG
+ if (q.state() != past_the_end)
+ Assert (q->used(), typename TriaAccessor<dim>::ExcUnusedCellAsChild());
+#endif
+ return q;
+};
+
+
+
+
+/* ------------------------ MGDoFQuadAccessor --------------------------- */
+
+template <int dim, typename BaseClass>
+MGDoFQuadAccessor<dim,BaseClass>::MGDoFQuadAccessor (Triangulation<dim> *tria,
+ const int level,
+ const int index,
+ const void *local_data) :
+ MGDoFAccessor<dim> ((MGDoFHandler<dim>*)local_data),
+ DoFQuadAccessor(tria,level,index,
+ static_cast<DoFHandler<dim>*>
+ (reinterpret_cast<MGDoFHandler<dim>*>(local_data))) {};
+
+
+
+template <int dim, typename BaseClass>
+inline
+int MGDoFQuadAccessor<dim,BaseClass>::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 (i<dof_handler->get_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 <int dim, typename BaseClass>
+void MGDoFQuadAccessor<dim,BaseClass>::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 (i<dof_handler->get_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 <int dim, typename BaseClass>
+inline
+int MGDoFQuadAccessor<dim,BaseClass>::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 (i<dof_handler->get_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 <int dim, typename BaseClass>
+void MGDoFQuadAccessor<dim,BaseClass>::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 (i<dof_handler->get_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 <int dim, typename BaseClass>
+void
+MGDoFQuadAccessor<dim,BaseClass>::get_mg_dof_indices (vector<int> &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<int>::iterator next = dof_indices.begin();
+ for (unsigned int vertex=0; vertex<4; ++vertex)
+ for (unsigned int d=0; d<dofs_per_vertex; ++d)
+ *next++ = mg_vertex_dof_index(vertex,d);
+ for (unsigned int line=0; line<4; ++line)
+ for (unsigned int d=0; d<dofs_per_line; ++d)
+ *next++ = this->line(line)->mg_dof_index(d);
+ for (unsigned int d=0; d<dofs_per_quad; ++d)
+ *next++ = mg_dof_index(d);
+};
+
+
+
+template <int dim, typename BaseClass>
+TriaIterator<dim,MGDoFLineAccessor<dim,LineAccessor<dim> > >
+MGDoFQuadAccessor<dim,BaseClass>::line (const unsigned int i) const {
+ Assert (i<4, ExcInvalidIndex (i, 0, 4));
+
+ return TriaIterator<dim,MGDoFLineAccessor<dim,LineAccessor<dim> > >
+ (
+ tria,
+ present_level,
+ line_index (i),
+ mg_dof_handler
+ );
+};
+
+
+
+template <int dim, typename BaseClass>
+TriaIterator<dim,MGDoFQuadAccessor<dim,BaseClass> >
+MGDoFQuadAccessor<dim,BaseClass>::child (const unsigned int i) const {
+ TriaIterator<dim,MGDoFQuadAccessor<dim,BaseClass> > q (tria,
+ present_level+1,
+ child_index (i),
+ mg_dof_handler);
+
+#ifdef DEBUG
+ if (q.state() != past_the_end)
+ Assert (q->used(), typename TriaAccessor<dim>::ExcUnusedCellAsChild());
+#endif
+ return q;
+};
+
+
+
+
+
+/*------------------------- Functions: MGDoFCellAccessor -----------------------*/
+
+template <int dim>
+TriaIterator<dim,MGDoFCellAccessor<dim> >
+MGDoFCellAccessor<dim>::neighbor (const unsigned int i) const {
+ TriaIterator<dim,MGDoFCellAccessor<dim> > q (tria,
+ neighbor_level (i),
+ neighbor_index (i),
+ mg_dof_handler);
+
+#ifdef DEBUG
+ if (q.state() != past_the_end)
+ Assert (q->used(), typename TriaAccessor<dim>::ExcUnusedCellAsNeighbor());
+#endif
+ return q;
+};
+
+
+
+template <int dim>
+TriaIterator<dim,MGDoFCellAccessor<dim> >
+MGDoFCellAccessor<dim>::child (const unsigned int i) const {
+ TriaIterator<dim,MGDoFCellAccessor<dim> > q (tria,
+ present_level+1,
+ child_index (i),
+ mg_dof_handler);
+
+#ifdef DEBUG
+ if (q.state() != past_the_end)
+ Assert (q->used(), typename TriaAccessor<dim>::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<double>::iterator next_dof_value=dof_values.begin();
+ for (unsigned int vertex=0; vertex<2; ++vertex)
+ for (unsigned int d=0; d<dofs_per_vertex; ++d)
+ *next_dof_value++ = values(mg_vertex_dof_index(vertex,d));
+ for (unsigned int d=0; d<dofs_per_line; ++d)
+ *next_dof_value++ = values(mg_dof_index(d));
+};
+
+#endif
+
+
+
+#if deal_II_dimension == 2
+
+template <>
+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<double>::iterator next_dof_value=dof_values.begin();
+ for (unsigned int vertex=0; vertex<4; ++vertex)
+ for (unsigned int d=0; d<dofs_per_vertex; ++d)
+ *next_dof_value++ = values(mg_vertex_dof_index(vertex,d));
+ for (unsigned int line=0; line<4; ++line)
+ for (unsigned int d=0; d<dofs_per_line; ++d)
+ *next_dof_value++ = values(this->line(line)->mg_dof_index(d));
+ for (unsigned int d=0; d<dofs_per_quad; ++d)
+ *next_dof_value++ = values(mg_dof_index(d));
+};
+
+#endif
+
+
+
+
+
+// explicit instantiations
+#if deal_II_dimension == 1
+template class MGDoFLineAccessor<1,CellAccessor<1> >;
+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
+
+
--- /dev/null
+/* $Id$ */
+
+#include <grid/mg_dof.h>
+#include <grid/dof_levels.h>
+#include <grid/dof_accessor.h>
+#include <grid/dof_constraints.h>
+#include <grid/tria_levels.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria.h>
+#include <grid/geometry_info.h>
+#include <fe/fe.h>
+
+
+
+/* ------------------------ MGVertexDoFs ----------------------------------- */
+
+template <int dim>
+MGDoFHandler<dim>::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<n_levels*dofs_per_vertex; ++i)
+ indices[i] = -1;
+};
+
+
+
+template <int dim>
+MGDoFHandler<dim>::MGVertexDoFs::~MGVertexDoFs () {
+ delete[] indices;
+};
+
+
+
+
+
+
+
+/* ------------------------ MGDoFHandler ------------------------------------- */
+
+template <int dim>
+MGDoFHandler<dim>::MGDoFHandler (Triangulation<dim> *tria) :
+ DoFHandler<dim> (tria)
+{};
+
+
+
+template <int dim>
+MGDoFHandler<dim>::~MGDoFHandler () {};
+
+
+
+
+template <int dim>
+void MGDoFHandler<dim>::distribute_dofs (const FiniteElementBase<dim> &fe) {
+ // first distribute global dofs
+ DoFHandler<dim>::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<deal_II_dimension>;
--- /dev/null
+/* $Id$ */
+
+
+#include <grid/mg_dof_accessor.h>
+#include <grid/mg_dof.h>
+#include <grid/dof_levels.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_iterator.templates.h>
+#include <fe/fe.h>
+
+#include <lac/dvector.h>
+#include <lac/dfmatrix.h>
+#include <lac/dsmatrix.h>
+
+
+
+
+
+
+/* ------------------------ MGDoFLineAccessor --------------------------- */
+
+template <int dim, typename BaseClass>
+MGDoFLineAccessor<dim,BaseClass>::MGDoFLineAccessor (Triangulation<dim> *tria,
+ const int level,
+ const int index,
+ const void *local_data) :
+ MGDoFAccessor<dim> ((MGDoFHandler<dim>*)local_data),
+ DoFLineAccessor(tria,level,index,
+ static_cast<DoFHandler<dim>*>
+ (reinterpret_cast<MGDoFHandler<dim>*>(local_data))) {};
+
+
+
+template <int dim, typename BaseClass>
+int MGDoFLineAccessor<dim,BaseClass>::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 (i<dof_handler->get_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 <int dim, typename BaseClass>
+void MGDoFLineAccessor<dim,BaseClass>::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 (i<dof_handler->get_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 dim, typename BaseClass>
+int MGDoFLineAccessor<dim,BaseClass>::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 (i<dof_handler->get_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 <int dim, typename BaseClass>
+void MGDoFLineAccessor<dim,BaseClass>::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 (i<dof_handler->get_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 <int dim, typename BaseClass>
+void
+MGDoFLineAccessor<dim,BaseClass>::get_mg_dof_indices (vector<int> &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<int>::iterator next = dof_indices.begin();
+ for (unsigned int vertex=0; vertex<2; ++vertex)
+ for (unsigned int d=0; d<dofs_per_vertex; ++d)
+ *next++ = mg_vertex_dof_index(vertex,d);
+ for (unsigned int d=0; d<dofs_per_line; ++d)
+ *next++ = mg_dof_index(d);
+};
+
+
+
+template <int dim, typename BaseClass>
+TriaIterator<dim,MGDoFLineAccessor<dim,BaseClass> >
+MGDoFLineAccessor<dim,BaseClass>::child (const unsigned int i) const {
+ TriaIterator<dim,MGDoFLineAccessor<dim,BaseClass> > q (tria,
+ present_level+1,
+ child_index (i),
+ mg_dof_handler);
+
+#ifdef DEBUG
+ if (q.state() != past_the_end)
+ Assert (q->used(), typename TriaAccessor<dim>::ExcUnusedCellAsChild());
+#endif
+ return q;
+};
+
+
+
+
+/* ------------------------ MGDoFQuadAccessor --------------------------- */
+
+template <int dim, typename BaseClass>
+MGDoFQuadAccessor<dim,BaseClass>::MGDoFQuadAccessor (Triangulation<dim> *tria,
+ const int level,
+ const int index,
+ const void *local_data) :
+ MGDoFAccessor<dim> ((MGDoFHandler<dim>*)local_data),
+ DoFQuadAccessor(tria,level,index,
+ static_cast<DoFHandler<dim>*>
+ (reinterpret_cast<MGDoFHandler<dim>*>(local_data))) {};
+
+
+
+template <int dim, typename BaseClass>
+inline
+int MGDoFQuadAccessor<dim,BaseClass>::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 (i<dof_handler->get_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 <int dim, typename BaseClass>
+void MGDoFQuadAccessor<dim,BaseClass>::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 (i<dof_handler->get_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 <int dim, typename BaseClass>
+inline
+int MGDoFQuadAccessor<dim,BaseClass>::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 (i<dof_handler->get_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 <int dim, typename BaseClass>
+void MGDoFQuadAccessor<dim,BaseClass>::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 (i<dof_handler->get_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 <int dim, typename BaseClass>
+void
+MGDoFQuadAccessor<dim,BaseClass>::get_mg_dof_indices (vector<int> &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<int>::iterator next = dof_indices.begin();
+ for (unsigned int vertex=0; vertex<4; ++vertex)
+ for (unsigned int d=0; d<dofs_per_vertex; ++d)
+ *next++ = mg_vertex_dof_index(vertex,d);
+ for (unsigned int line=0; line<4; ++line)
+ for (unsigned int d=0; d<dofs_per_line; ++d)
+ *next++ = this->line(line)->mg_dof_index(d);
+ for (unsigned int d=0; d<dofs_per_quad; ++d)
+ *next++ = mg_dof_index(d);
+};
+
+
+
+template <int dim, typename BaseClass>
+TriaIterator<dim,MGDoFLineAccessor<dim,LineAccessor<dim> > >
+MGDoFQuadAccessor<dim,BaseClass>::line (const unsigned int i) const {
+ Assert (i<4, ExcInvalidIndex (i, 0, 4));
+
+ return TriaIterator<dim,MGDoFLineAccessor<dim,LineAccessor<dim> > >
+ (
+ tria,
+ present_level,
+ line_index (i),
+ mg_dof_handler
+ );
+};
+
+
+
+template <int dim, typename BaseClass>
+TriaIterator<dim,MGDoFQuadAccessor<dim,BaseClass> >
+MGDoFQuadAccessor<dim,BaseClass>::child (const unsigned int i) const {
+ TriaIterator<dim,MGDoFQuadAccessor<dim,BaseClass> > q (tria,
+ present_level+1,
+ child_index (i),
+ mg_dof_handler);
+
+#ifdef DEBUG
+ if (q.state() != past_the_end)
+ Assert (q->used(), typename TriaAccessor<dim>::ExcUnusedCellAsChild());
+#endif
+ return q;
+};
+
+
+
+
+
+/*------------------------- Functions: MGDoFCellAccessor -----------------------*/
+
+template <int dim>
+TriaIterator<dim,MGDoFCellAccessor<dim> >
+MGDoFCellAccessor<dim>::neighbor (const unsigned int i) const {
+ TriaIterator<dim,MGDoFCellAccessor<dim> > q (tria,
+ neighbor_level (i),
+ neighbor_index (i),
+ mg_dof_handler);
+
+#ifdef DEBUG
+ if (q.state() != past_the_end)
+ Assert (q->used(), typename TriaAccessor<dim>::ExcUnusedCellAsNeighbor());
+#endif
+ return q;
+};
+
+
+
+template <int dim>
+TriaIterator<dim,MGDoFCellAccessor<dim> >
+MGDoFCellAccessor<dim>::child (const unsigned int i) const {
+ TriaIterator<dim,MGDoFCellAccessor<dim> > q (tria,
+ present_level+1,
+ child_index (i),
+ mg_dof_handler);
+
+#ifdef DEBUG
+ if (q.state() != past_the_end)
+ Assert (q->used(), typename TriaAccessor<dim>::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<double>::iterator next_dof_value=dof_values.begin();
+ for (unsigned int vertex=0; vertex<2; ++vertex)
+ for (unsigned int d=0; d<dofs_per_vertex; ++d)
+ *next_dof_value++ = values(mg_vertex_dof_index(vertex,d));
+ for (unsigned int d=0; d<dofs_per_line; ++d)
+ *next_dof_value++ = values(mg_dof_index(d));
+};
+
+#endif
+
+
+
+#if deal_II_dimension == 2
+
+template <>
+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<double>::iterator next_dof_value=dof_values.begin();
+ for (unsigned int vertex=0; vertex<4; ++vertex)
+ for (unsigned int d=0; d<dofs_per_vertex; ++d)
+ *next_dof_value++ = values(mg_vertex_dof_index(vertex,d));
+ for (unsigned int line=0; line<4; ++line)
+ for (unsigned int d=0; d<dofs_per_line; ++d)
+ *next_dof_value++ = values(this->line(line)->mg_dof_index(d));
+ for (unsigned int d=0; d<dofs_per_quad; ++d)
+ *next_dof_value++ = values(mg_dof_index(d));
+};
+
+#endif
+
+
+
+
+
+// explicit instantiations
+#if deal_II_dimension == 1
+template class MGDoFLineAccessor<1,CellAccessor<1> >;
+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
+
+
--- /dev/null
+/* $Id$ */
+
+#include <grid/mg_dof.h>
+#include <grid/dof_levels.h>
+#include <grid/dof_accessor.h>
+#include <grid/dof_constraints.h>
+#include <grid/tria_levels.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria.h>
+#include <grid/geometry_info.h>
+#include <fe/fe.h>
+
+
+
+/* ------------------------ MGVertexDoFs ----------------------------------- */
+
+template <int dim>
+MGDoFHandler<dim>::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<n_levels*dofs_per_vertex; ++i)
+ indices[i] = -1;
+};
+
+
+
+template <int dim>
+MGDoFHandler<dim>::MGVertexDoFs::~MGVertexDoFs () {
+ delete[] indices;
+};
+
+
+
+
+
+
+
+/* ------------------------ MGDoFHandler ------------------------------------- */
+
+template <int dim>
+MGDoFHandler<dim>::MGDoFHandler (Triangulation<dim> *tria) :
+ DoFHandler<dim> (tria)
+{};
+
+
+
+template <int dim>
+MGDoFHandler<dim>::~MGDoFHandler () {};
+
+
+
+
+template <int dim>
+void MGDoFHandler<dim>::distribute_dofs (const FiniteElementBase<dim> &fe) {
+ // first distribute global dofs
+ DoFHandler<dim>::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<deal_II_dimension>;