#include <lac/sparsematrix.h>
#include <base/parameter_handler.h>
#include <grid/dof_constraints.h>
-
+#include <numerics/dof_renumbering.h>
#include <fstream>
#include <cmath>
dof->distribute_dofs (fe);
cout << " Renumbering degrees of freedom..." << endl;
- dof->renumber_dofs (Cuthill_McKee, false);
+ DoFRenumbering::renumber_Cuthill_McKee (*dof);
SparseMatrixStruct sparsity (dof->n_dofs(),
dof->max_couplings_between_dofs());
-/**
- * Give names to the different possibilities of renumbering the degrees
- * of freedom.
- *
- * \begin{itemize}
- * \item #Cuthill_McKee# and #reverse_Cuthill_McKee# traverse the triangulation
- * in a diagonal, advancing front like method and produce matrices with an
- * almost minimal bandwidth.
- * \item #reverse_Cuthill_McKey# does the same thing, but numbers the dofs in
- * the reverse order.
- * \end{itemize}
- *
- * For a description of the algorithms see the book of Schwarz (H.R.Scharz:
- * Methode der finiten Elemente).
- */
-enum RenumberingMethod {
- Cuthill_McKee,
- reverse_Cuthill_McKee
-};
-
-
-
-
/**
* Manage the distribution and numbering of the degrees of freedom for
*
* This numbering implies very large bandwiths of the resulting matrices and
* is thus vastly suboptimal for some solution algorithms. For this reason,
- * the #DoFHandler# class offers the function #renumber_dofs# which reorders
- * the dof numbering according to some scheme. Presently available are the
- * Cuthill-McKey (CM) and the Reverse Cuthill-McKey algorithm. These algorithms
- * have one major drawback: they require a good starting point, i.e. the degree
- * of freedom index afterwards to be numbered zero. This can thus be given by
- * the user, e.g. by exploiting knowledge of the actual topology of the
- * domain. It is also possible to give several starting indices, which may
- * be used to simulate a simple upstream numbering (by giving the inflow
- * dofs as starting values) or to make preconditioning faster (by letting
- * the dirichlet boundary indices be starting points).
- *
- * If no starting index is given, one is chosen by the program, namely one
- * with the smallest coordination number (the coordination number is the
- * number of other dofs this dof couples with). This dof is usually located
- * on the boundary of the domain. There is, however, large ambiguity in this
- * when using the hierarchical meshes used in this library, since in most
- * cases the computational domain is not approximated by tilting and deforming
- * elements and by plugging together variable numbers of elements at vertices,
- * but rather by hierarchical refinement. There is therefore a large number
- * of dofs with equal coordination numbers. The renumbering algorithms will
- * therefore not give optimal results.
- *
- * In the book of Schwarz (H.R.Schwarz: Methode der finiten Elemente), it is
- * advised to test many starting points, if possible all with the smallest
- * coordination number and also those with slightly higher numbers. However,
- * this seems only possible for meshes with at most several dozen or a few
- * hundred elements found in small engineering problems of the early 1980s
- * (the second edition was published in 1984), but certainly not with those
- * used in this library, featuring several 10,000 to a few 100,000 elements.
- *
- * On the other hand, the need to reduce the bandwidth has decreased since
- * with the mentioned number of cells, only iterative solution methods are
- * able to solve the resulting matrix systems. These, however, are not so
- * demanding with respect to the bandwidth as direct solvers used for
- * smaller problems. Things like upstream numbering become much more important
- * in recent times, so the suboptimality of the renumbering algorithms is
- * not that important any more.
- *
- *
- * \subsection{Implementation of renumbering schemes}
- *
- * The renumbering algorithms need quite a lot of memory, since they have
- * to store for each dof with which other dofs it couples. This is done
- * using a #SparseMatrixStruct# object used to store the sparsity pattern of
- * matrices. It
- * is not useful for the user to do anything between distributing the dofs
- * and renumbering, i.e. the calls to #DoFHandler::distribute_dofs# and
- * #DoFHandler::renumber_dofs# should follow each other immediately. If
- * you try to create a sparsity pattern or anything else in between, these
- * will be invalid afterwards.
- *
- * The renumbering may take care of dof-to-dof couplings only induced by
- * eliminating constraints. In addition to the memory consumption mentioned
- * above, this also takes quite some computational time, but it may be
- * switched off upon calling the #renumber_dofs# function. This will then
- * give inferior results, since knots in the graph (representing dofs)
- * are not found to be neighbors even if they would be after condensation.
- *
- * The renumbering algorithms work on a purely algebraic basis, due to the
- * isomorphism between the graph theoretical groundwork underlying the
- * algorithms and binary matrices (matrices of which the entries are binary
- * values) represented by the sparsity patterns. In special, the algorithms
- * do not try to exploit topological knowledge (e.g. corner detection) to
- * find appropriate starting points. This way, however, they work in
- * arbitrary space dimension.
- *
- * If you want to give starting points, you may give a list of dof indices
- * which will form the first step of the renumbering. The dofs of the list
- * will be consecutively numbered starting with zero, i.e. this list is not
- * renumbered according to the coordination number of the nodes. Indices not
- * in the allowed range are deleted. If no index is allowed, the algorithm
- * will search for its own starting point.
- *
- *
- * \subsection{Results of renumbering}
- *
- * The renumbering schemes mentioned above do not lead to optimal results.
- * However, after all there is no algorithm that accomplishes this within
- * reasonable time. There are situations where the lack of optimality even
- * leads to worse results than with the original, crude, levelwise numering
- * scheme; one of these examples is a mesh of four cells of which always
- * those cells are refined which are neighbors to the center (you may call
- * this mesh a `zoom in' mesh). In one such example the bandwidth was
- * increased by about 50 per cent.
- *
- * In most other cases, the bandwith is reduced significantly. The reduction
- * is the better the less structured the grid is. With one grid where the
- * cells were refined according to a random driven algorithm, the bandwidth
- * was reduced by a factor of six.
- *
- * Using the constraint information usually leads to reductions in bandwidth
- * of 10 or 20 per cent, but may for some very unstructured grids also lead
- * to an increase. You have to weigh the decrease in your case with the time
- * spent to use the constraint information, which usually is several times
- * longer than the `pure' renumbering algorithm.
- *
- * In almost all cases, the renumbering scheme finds a corner to start with.
- * Since there is more than one corner in most grids and since even an
- * interior degree of freedom may be a better starting point, giving the
- * starting point by the user may be a viable way if you have a simple
- * scheme to derive a suitable point (e.g. by successively taking the
- * third child of the cell top left of the coarsest level, taking its
- * third vertex and the dof index thereof, if you want the top left corner
- * vertex). If you do not know beforehand what your grid will look like
- * (e.g. when using adaptive algorithms), searching a best starting point
- * may be difficult, however, and in many cases will not justify the effort.
+ * the #DoFRenumbering# class offers the function #renumber_dofs# which reorders
+ * the dof numbering according to some scheme. See there for a discussion of
+ * the implemented algorithms.
*
*
* \subsection{User defined renumbering schemes}
*
- * The #renumber_dofs# function offers a fixed number of renumbering
+ * The #DoFRenumbering# class offers a fixed number of renumbering
* schemes like the Cuthill-McKey scheme. Basically, the function sets
* up an array in which for each degree of freedom the index is stored
* which is to be assigned by the renumbering. Using this array, the
*/
virtual void clear ();
- /**
- * 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 this
- * class for more details.
- */
- void renumber_dofs (const RenumberingMethod method,
- const bool use_constraints = false,
- const vector<int> &starting_points = vector<int>());
-
/**
* Actually do the renumbering based on
* a list of new dof numbers for all the
* Exception
*/
DeclException0 (ExcNoFESelected);
- /**
+ /**
* Exception
*/
DeclException0 (ExcRenumberingIncomplete);
* Return a pointer to the #i#th quad
* bounding this #Hex#.
*/
- TriaIterator<dim,MGDoFLineAccessor<dim,QuadAccessor<dim> > >
+ TriaIterator<dim,MGDoFQuadAccessor<dim,QuadAccessor<dim> > >
quad (const unsigned int i) const;
/**
+
+
/*---------------------------- mg_dof_accessor.h ---------------------------*/
/* end of #ifndef __mg_dof_accessor_H */
#endif
typedef typename MGDoFDimensionInfo<dim>::quad_iterator quad_iterator;
typedef typename MGDoFDimensionInfo<dim>::active_quad_iterator active_quad_iterator;
+ typedef typename MGDoFDimensionInfo<dim>::raw_hex_iterator raw_hex_iterator;
+ typedef typename MGDoFDimensionInfo<dim>::hex_iterator hex_iterator;
+ typedef typename MGDoFDimensionInfo<dim>::active_hex_iterator active_hex_iterator;
+
typedef typename MGDoFDimensionInfo<dim>::raw_cell_iterator raw_cell_iterator;
typedef typename MGDoFDimensionInfo<dim>::cell_iterator cell_iterator;
typedef typename MGDoFDimensionInfo<dim>::active_cell_iterator active_cell_iterator;
virtual void distribute_dofs (const FiniteElement<dim> &);
/**
- * Renumber the degrees of freedom according
- * to the given scheme, eventually
- * using the given
- * starting points. The starting points
- * default to an empty list, the use of
- * constraint information defaults to
- * false.
+ * Actually do the renumbering based on
+ * a list of new dof numbers for all the
+ * dofs.
*
- * See the general documentation of the
- * #DoFHandler# class for more details.
+ * #new_numbers# is an array of integers
+ * with size equal to the number of dofs
+ * on the present level. It stores the new
+ * indices after renumbering in the
+ * order of the old indices.
*/
- void renumber_dofs (const unsigned int level,
- const RenumberingMethod method,
- const vector<int> &starting_points = vector<int>());
+ void renumber_dofs (const unsigned int level,
+ const vector<int> &new_numbers);
/**
* Write the sparsity structure of the
/*---------------------------------------*/
+ /**
+ * @name Hex iterator functions*/
+ /*@{
+ */
+ /**
+ * Return iterator to the first hex, used
+ * or not, on level #level#. If a level
+ * has no hexs, a past-the-end iterator
+ * is returned.
+ */
+ raw_hex_iterator begin_raw_hex (const unsigned int level = 0) const;
+
+ /**
+ * Return iterator to the first used hex
+ * on level #level#.
+ */
+ hex_iterator begin_hex (const unsigned int level = 0) const;
+
+ /**
+ * Return iterator to the first active
+ * hex on level #level#.
+ */
+ active_hex_iterator begin_active_hex(const unsigned int level = 0) const;
+
+ /**
+ * Return iterator past the end; this
+ * iterator serves for comparisons of
+ * iterators with past-the-end or
+ * before-the-beginning states.
+ */
+ raw_hex_iterator end_hex () const;
+
+ /**
+ * Return an iterator which is the first
+ * iterator not on level. If #level# is
+ * the last level, then this returns
+ * #end()#.
+ */
+ hex_iterator end_hex (const unsigned int level) const;
+
+ /**
+ * Return a raw iterator which is the first
+ * iterator not on level. If #level# is
+ * the last level, then this returns
+ * #end()#.
+ */
+ raw_hex_iterator end_raw_hex (const unsigned int level) const;
+
+ /**
+ * Return an active iterator which is the
+ * first iterator not on level. If #level#
+ * is the last level, then this returns
+ * #end()#.
+ */
+ active_hex_iterator end_active_hex (const unsigned int level) const;
+
+
+ /**
+ * Return an iterator pointing to the
+ * last hex, used or not.
+ */
+ raw_hex_iterator last_raw_hex () const;
+
+ /**
+ * Return an iterator pointing to the last
+ * hex of the level #level#, used or not.
+
+ */
+ raw_hex_iterator last_raw_hex (const unsigned int level) const;
+
+ /**
+ * Return an iterator pointing to the last
+ * used hex.
+ */
+ hex_iterator last_hex () const;
+
+ /**
+ * Return an iterator pointing to the last
+ * used hex on level #level#.
+ */
+ hex_iterator last_hex (const unsigned int level) const;
+
+ /**
+ * Return an iterator pointing to the last
+ * active hex.
+ */
+ active_hex_iterator last_active_hex () const;
+
+ /**
+ * Return an iterator pointing to the last
+ * active hex on level #level#.
+ */
+ active_hex_iterator last_active_hex (const unsigned int level) const;
+ /*@}*/
+
+ /*---------------------------------------*/
+
/**
* Return the number of degrees of freedom
* on the specified level.
unsigned int distribute_dofs_on_cell (cell_iterator &cell,
unsigned int next_free_dof);
- /**
- * Actually do the renumbering prepared
- * by the #renumber_dofs# function on
- * the given #level#. Since
- * this is dimension specific, we
- * need to have another function.
- *
- * #new_numbers# is an array of integers
- * with size equal to the number of dofs
- * on the present level. It stores the new
- * indices after renumbering in the
- * order of the old indices.
- */
- void do_renumbering (const unsigned int level,
- const vector<int> &new_numbers);
-
/**
* Reserve enough space for the MG dof
* indices for a given triangulation.
friend class MGDoFLineAccessor<dim, CellAccessor<dim> >;
friend class MGDoFQuadAccessor<dim, QuadAccessor<dim> >;
friend class MGDoFQuadAccessor<dim, CellAccessor<dim> >;
+ friend class MGDoFHexAccessor<dim, HexAccessor<dim> >;
+ friend class MGDoFHexAccessor<dim, CellAccessor<dim> >;
};
+template <int dim>
+inline
+unsigned int
+HexAccessor<dim>::line_index (unsigned int i) const {
+ Assert (i<12, ExcInvalidIndex(i,0,11));
+
+ if (i<4)
+ return quad(0)->line_index(i);
+ else
+ if (i<8)
+ return quad(1)->line_index(i-4);
+ else
+ switch (i)
+ {
+ case 8:
+ return quad(2)->line_index(3);
+ case 9:
+ return quad(2)->line_index(1);
+ case 10:
+ return quad(4)->line_index(1);
+ case 11:
+ return quad(4)->line_index(3);
+ };
+ Assert (false, ExcInvalidIndex(i,0,11));
+ return 0;
+};
+
+
+
template <int dim>
inline
unsigned int
* Return a pointer to the #i#th quad
* bounding this #Hex#.
*/
- TriaIterator<dim,MGDoFLineAccessor<dim,QuadAccessor<dim> > >
+ TriaIterator<dim,MGDoFQuadAccessor<dim,QuadAccessor<dim> > >
quad (const unsigned int i) const;
/**
+
+
/*---------------------------- mg_dof_accessor.h ---------------------------*/
/* end of #ifndef __mg_dof_accessor_H */
#endif
typedef typename MGDoFDimensionInfo<dim>::quad_iterator quad_iterator;
typedef typename MGDoFDimensionInfo<dim>::active_quad_iterator active_quad_iterator;
+ typedef typename MGDoFDimensionInfo<dim>::raw_hex_iterator raw_hex_iterator;
+ typedef typename MGDoFDimensionInfo<dim>::hex_iterator hex_iterator;
+ typedef typename MGDoFDimensionInfo<dim>::active_hex_iterator active_hex_iterator;
+
typedef typename MGDoFDimensionInfo<dim>::raw_cell_iterator raw_cell_iterator;
typedef typename MGDoFDimensionInfo<dim>::cell_iterator cell_iterator;
typedef typename MGDoFDimensionInfo<dim>::active_cell_iterator active_cell_iterator;
virtual void distribute_dofs (const FiniteElement<dim> &);
/**
- * Renumber the degrees of freedom according
- * to the given scheme, eventually
- * using the given
- * starting points. The starting points
- * default to an empty list, the use of
- * constraint information defaults to
- * false.
+ * Actually do the renumbering based on
+ * a list of new dof numbers for all the
+ * dofs.
*
- * See the general documentation of the
- * #DoFHandler# class for more details.
+ * #new_numbers# is an array of integers
+ * with size equal to the number of dofs
+ * on the present level. It stores the new
+ * indices after renumbering in the
+ * order of the old indices.
*/
- void renumber_dofs (const unsigned int level,
- const RenumberingMethod method,
- const vector<int> &starting_points = vector<int>());
+ void renumber_dofs (const unsigned int level,
+ const vector<int> &new_numbers);
/**
* Write the sparsity structure of the
/*---------------------------------------*/
+ /**
+ * @name Hex iterator functions*/
+ /*@{
+ */
+ /**
+ * Return iterator to the first hex, used
+ * or not, on level #level#. If a level
+ * has no hexs, a past-the-end iterator
+ * is returned.
+ */
+ raw_hex_iterator begin_raw_hex (const unsigned int level = 0) const;
+
+ /**
+ * Return iterator to the first used hex
+ * on level #level#.
+ */
+ hex_iterator begin_hex (const unsigned int level = 0) const;
+
+ /**
+ * Return iterator to the first active
+ * hex on level #level#.
+ */
+ active_hex_iterator begin_active_hex(const unsigned int level = 0) const;
+
+ /**
+ * Return iterator past the end; this
+ * iterator serves for comparisons of
+ * iterators with past-the-end or
+ * before-the-beginning states.
+ */
+ raw_hex_iterator end_hex () const;
+
+ /**
+ * Return an iterator which is the first
+ * iterator not on level. If #level# is
+ * the last level, then this returns
+ * #end()#.
+ */
+ hex_iterator end_hex (const unsigned int level) const;
+
+ /**
+ * Return a raw iterator which is the first
+ * iterator not on level. If #level# is
+ * the last level, then this returns
+ * #end()#.
+ */
+ raw_hex_iterator end_raw_hex (const unsigned int level) const;
+
+ /**
+ * Return an active iterator which is the
+ * first iterator not on level. If #level#
+ * is the last level, then this returns
+ * #end()#.
+ */
+ active_hex_iterator end_active_hex (const unsigned int level) const;
+
+
+ /**
+ * Return an iterator pointing to the
+ * last hex, used or not.
+ */
+ raw_hex_iterator last_raw_hex () const;
+
+ /**
+ * Return an iterator pointing to the last
+ * hex of the level #level#, used or not.
+
+ */
+ raw_hex_iterator last_raw_hex (const unsigned int level) const;
+
+ /**
+ * Return an iterator pointing to the last
+ * used hex.
+ */
+ hex_iterator last_hex () const;
+
+ /**
+ * Return an iterator pointing to the last
+ * used hex on level #level#.
+ */
+ hex_iterator last_hex (const unsigned int level) const;
+
+ /**
+ * Return an iterator pointing to the last
+ * active hex.
+ */
+ active_hex_iterator last_active_hex () const;
+
+ /**
+ * Return an iterator pointing to the last
+ * active hex on level #level#.
+ */
+ active_hex_iterator last_active_hex (const unsigned int level) const;
+ /*@}*/
+
+ /*---------------------------------------*/
+
/**
* Return the number of degrees of freedom
* on the specified level.
unsigned int distribute_dofs_on_cell (cell_iterator &cell,
unsigned int next_free_dof);
- /**
- * Actually do the renumbering prepared
- * by the #renumber_dofs# function on
- * the given #level#. Since
- * this is dimension specific, we
- * need to have another function.
- *
- * #new_numbers# is an array of integers
- * with size equal to the number of dofs
- * on the present level. It stores the new
- * indices after renumbering in the
- * order of the old indices.
- */
- void do_renumbering (const unsigned int level,
- const vector<int> &new_numbers);
-
/**
* Reserve enough space for the MG dof
* indices for a given triangulation.
friend class MGDoFLineAccessor<dim, CellAccessor<dim> >;
friend class MGDoFQuadAccessor<dim, QuadAccessor<dim> >;
friend class MGDoFQuadAccessor<dim, CellAccessor<dim> >;
+ friend class MGDoFHexAccessor<dim, HexAccessor<dim> >;
+ friend class MGDoFHexAccessor<dim, CellAccessor<dim> >;
};
--- /dev/null
+/*---------------------------- dof_renumbering.h ---------------------------*/
+/* $Id$ */
+#ifndef __dof_renumbering_H
+#define __dof_renumbering_H
+/*---------------------------- dof_renumbering.h ---------------------------*/
+
+
+#include <base/exceptions.h>
+#include <basic/forward-declarations.h>
+
+
+
+
+/**
+ * Implementation of a number of renumbering algorithms for the degrees of
+ * freedom on a triangulation.
+ *
+ * \section{Cuthill-McKee like algorithms}
+ *
+ * Within this class, the Cuthill-McKee algorithm is implemented. It starts
+ * at a degree of freedom, searches the other DoFs for those which are couple
+ * with the one we started with and numbers these in a certain way. It then
+ * finds the second level of DoFs, namely those that couple with those of
+ * the previous level (which were those that coupled with the initial DoF)
+ * and numbers these. And so on. For the details of the algorithm, especially
+ * the numbering within each level, we refer the reader to the book of
+ * Schwarz (H.R.Schwarz: Methode der finiten Elemente). The reverse Cuthill-McKee
+ * algorithm does the same job, but numbers all elements in the reverse order.
+ *
+ * These algorithms
+ * have one major drawback: they require a good starting point, i.e. the degree
+ * of freedom index afterwards to be numbered zero. This can thus be given by
+ * the user, e.g. by exploiting knowledge of the actual topology of the
+ * domain. It is also possible to give several starting indices, which may
+ * be used to simulate a simple upstream numbering (by giving the inflow
+ * dofs as starting values) or to make preconditioning faster (by letting
+ * the dirichlet boundary indices be starting points).
+ *
+ * If no starting index is given, one is chosen by the program, namely one
+ * with the smallest coordination number (the coordination number is the
+ * number of other dofs this dof couples with). This dof is usually located
+ * on the boundary of the domain. There is, however, large ambiguity in this
+ * when using the hierarchical meshes used in this library, since in most
+ * cases the computational domain is not approximated by tilting and deforming
+ * elements and by plugging together variable numbers of elements at vertices,
+ * but rather by hierarchical refinement. There is therefore a large number
+ * of dofs with equal coordination numbers. The renumbering algorithms will
+ * therefore not give optimal results.
+ *
+ * In the book of Schwarz (H.R.Schwarz: Methode der finiten Elemente), it is
+ * advised to test many starting points, if possible all with the smallest
+ * coordination number and also those with slightly higher numbers. However,
+ * this seems only possible for meshes with at most several dozen or a few
+ * hundred elements found in small engineering problems of the early 1980s
+ * (the second edition was published in 1984), but certainly not with those
+ * used in this library, featuring several 10,000 to a few 100,000 elements.
+ *
+ * On the other hand, the need to reduce the bandwidth has decreased since
+ * with the mentioned number of cells, only iterative solution methods are
+ * able to solve the resulting matrix systems. These, however, are not so
+ * demanding with respect to the bandwidth as direct solvers used for
+ * smaller problems. Things like upstream numbering become much more important
+ * in recent times, so the suboptimality of the renumbering algorithms is
+ * not that important any more.
+ *
+ *
+ * \subsection{Implementation of renumbering schemes}
+ *
+ * The renumbering algorithms need quite a lot of memory, since they have
+ * to store for each dof with which other dofs it couples. This is done
+ * using a #SparseMatrixStruct# object used to store the sparsity pattern of
+ * matrices. It
+ * is not useful for the user to do anything between distributing the dofs
+ * and renumbering, i.e. the calls to #DoFHandler::distribute_dofs# and
+ * #DoFHandler::renumber_dofs# should follow each other immediately. If
+ * you try to create a sparsity pattern or anything else in between, these
+ * will be invalid afterwards.
+ *
+ * The renumbering may take care of dof-to-dof couplings only induced by
+ * eliminating constraints. In addition to the memory consumption mentioned
+ * above, this also takes quite some computational time, but it may be
+ * switched off upon calling the #renumber_dofs# function. This will then
+ * give inferior results, since knots in the graph (representing dofs)
+ * are not found to be neighbors even if they would be after condensation.
+ *
+ * The renumbering algorithms work on a purely algebraic basis, due to the
+ * isomorphism between the graph theoretical groundwork underlying the
+ * algorithms and binary matrices (matrices of which the entries are binary
+ * values) represented by the sparsity patterns. In special, the algorithms
+ * do not try to exploit topological knowledge (e.g. corner detection) to
+ * find appropriate starting points. This way, however, they work in
+ * arbitrary space dimension.
+ *
+ * If you want to give starting points, you may give a list of dof indices
+ * which will form the first step of the renumbering. The dofs of the list
+ * will be consecutively numbered starting with zero, i.e. this list is not
+ * renumbered according to the coordination number of the nodes. Indices not
+ * in the allowed range are deleted. If no index is allowed, the algorithm
+ * will search for its own starting point.
+ *
+ *
+ * \subsection{Results of renumbering}
+ *
+ * The renumbering schemes mentioned above do not lead to optimal results.
+ * However, after all there is no algorithm that accomplishes this within
+ * reasonable time. There are situations where the lack of optimality even
+ * leads to worse results than with the original, crude, levelwise numering
+ * scheme; one of these examples is a mesh of four cells of which always
+ * those cells are refined which are neighbors to the center (you may call
+ * this mesh a `zoom in' mesh). In one such example the bandwidth was
+ * increased by about 50 per cent.
+ *
+ * In most other cases, the bandwith is reduced significantly. The reduction
+ * is the better the less structured the grid is. With one grid where the
+ * cells were refined according to a random driven algorithm, the bandwidth
+ * was reduced by a factor of six.
+ *
+ * Using the constraint information usually leads to reductions in bandwidth
+ * of 10 or 20 per cent, but may for some very unstructured grids also lead
+ * to an increase. You have to weigh the decrease in your case with the time
+ * spent to use the constraint information, which usually is several times
+ * longer than the `pure' renumbering algorithm.
+ *
+ * In almost all cases, the renumbering scheme finds a corner to start with.
+ * Since there is more than one corner in most grids and since even an
+ * interior degree of freedom may be a better starting point, giving the
+ * starting point by the user may be a viable way if you have a simple
+ * scheme to derive a suitable point (e.g. by successively taking the
+ * third child of the cell top left of the coarsest level, taking its
+ * third vertex and the dof index thereof, if you want the top left corner
+ * vertex). If you do not know beforehand what your grid will look like
+ * (e.g. when using adaptive algorithms), searching a best starting point
+ * may be difficult, however, and in many cases will not justify the effort.
+ *
+ *
+ * \subsection{Multigrid DoF numbering}
+ *
+ * Most algorithms also work on multigrid degree of freedom numberings. Refer
+ * to the actual function declarations to get more information on this.
+ *
+ *
+ * @author Wolfgang Bangerth, 1998, 1999
+ */
+class DoFRenumbering
+{
+ public:
+ /**
+ * Renumber the degrees of freedom
+ * according to the Cuthill-McKee method,
+ * eventually using the reverse numbering
+ * scheme.
+ *
+ * See the general documentation of
+ * this class for details on the
+ * different methods.
+ */
+ template <int dim>
+ static void renumber_Cuthill_McKee (DoFHandler<dim> &dof_handler,
+ const bool reversed_numbering = false,
+ const bool use_constraints = false,
+ const vector<int> &starting_indices = vector<int>());
+
+ /**
+ * Renumber the degrees of freedom
+ * according to the Cuthill-McKee method,
+ * eventually using the reverse numbering
+ * scheme, in this case for a multigrid
+ * numbering of degrees of freedom.
+ *
+ * You can give a triangulation level to
+ * which this function is to be applied.
+ * Since with a level-wise numbering there
+ * are no hanging nodes, no constraints
+ * can be used, so the respective
+ * parameter of the previous function is
+ * ommitted.
+ *
+ * See the general documentation of
+ * this class for details on the
+ * different methods.
+ */
+ template <int dim>
+ static void renumber_Cuthill_McKee (MGDoFHandler<dim> &dof_handler,
+ const unsigned int level,
+ const bool reversed_numbering = false,
+ const vector<int> &starting_indices = vector<int> ());
+
+
+ /**
+ * Exception
+ */
+ DeclException0 (ExcRenumberingIncomplete);
+};
+
+
+
+
+
+
+
+
+/*---------------------------- dof_renumbering.h ---------------------------*/
+/* end of #ifndef __dof_renumbering_H */
+#endif
+/*---------------------------- dof_renumbering.h ---------------------------*/
#include <lac/vector.h>
#include <lac/vector.h>
-#include <map>
#include <set>
#include <algorithm>
-template <int dim>
-void DoFHandler<dim>::renumber_dofs (const RenumberingMethod method,
- const bool use_constraints,
- const vector<int> &starting_points) {
- // make the connection graph
- SparseMatrixStruct sparsity (n_dofs(), max_couplings_between_dofs());
- make_sparsity_pattern (sparsity);
-
- if (use_constraints)
- {
- ConstraintMatrix constraints;
- make_hanging_node_constraints (constraints);
- constraints.close ();
- constraints.condense (sparsity);
- };
-
- int n_dofs = sparsity.n_rows();
- // store the new dof numbers; -1 means
- // that no new number was chosen yet
- vector<int> new_number(sparsity.n_rows(), -1);
-
- // store the indices of the dofs renumbered
- // in the last round. Default to starting
- // points
- vector<int> last_round_dofs (starting_points);
-
- // delete disallowed elements
- for (unsigned int i=0; i<last_round_dofs.size(); ++i)
- if ((last_round_dofs[i]<0) || (last_round_dofs[i]>=n_dofs))
- last_round_dofs[i] = -1;
-
- remove_if (last_round_dofs.begin(), last_round_dofs.end(),
- bind2nd(equal_to<int>(), -1));
-
- // now if no valid points remain:
- // find dof with lowest coordination
- // number
-
- if (last_round_dofs.size() == 0)
- {
- int starting_point = -1;
- unsigned int min_coordination = n_dofs;
- for (int row=0; row<n_dofs; ++row)
- {
- unsigned int j;
-
- // loop until we hit the end
- // of this row's entries
- for (j=sparsity.get_rowstart_indices()[row];
- j<sparsity.get_rowstart_indices()[row+1]; ++j)
- if (sparsity.get_column_numbers()[j] == -1)
- break;
- // post-condition after loop:
- // coordination, i.e. the number
- // of entries in this row is now
- // j-rowstart[row]
- if (j-sparsity.get_rowstart_indices()[row] < min_coordination)
- {
- min_coordination = j-sparsity.get_rowstart_indices()[row];
- starting_point = row;
- };
- };
-
- // now we still have to care for the
- // case that no dof has a coordination
- // number less than n_dofs. this rather
- // exotic case only happens if we only
- // have one cell, as far as I can see,
- // but there may be others as well.
- //
- // if that should be the case, we can
- // chose an arbitrary dof as starting
- // point, e.g. the one with number zero
- if (starting_point == -1)
- starting_point = 0;
-
- // initialize the first dof
- last_round_dofs.push_back (starting_point);
- };
-
-
- // store next free dof index
- int next_free_number = 0;
-
- // enumerate the first round dofs
- for (unsigned int i=0; i!=last_round_dofs.size(); ++i)
- new_number[last_round_dofs[i]] = next_free_number++;
-
- bool all_dofs_renumbered = false;
-
- // now do as many steps as needed to
- // renumber all dofs
- while (!all_dofs_renumbered)
- {
- // store the indices of the dofs to be
- // renumbered in the next round
- vector<int> next_round_dofs;
-
- // find all neighbors of the
- // dofs numbered in the last
- // round
- for (unsigned int i=0; i<last_round_dofs.size(); ++i)
- for (unsigned int j=sparsity.get_rowstart_indices()[last_round_dofs[i]];
- j<sparsity.get_rowstart_indices()[last_round_dofs[i]+1]; ++j)
- if (sparsity.get_column_numbers()[j] == -1)
- break;
- else
- next_round_dofs.push_back (sparsity.get_column_numbers()[j]);
-
- // sort dof numbers
- sort (next_round_dofs.begin(), next_round_dofs.end());
-
- // delete multiple entries
- vector<int>::iterator end_sorted;
- end_sorted = unique (next_round_dofs.begin(), next_round_dofs.end());
- next_round_dofs.erase (end_sorted, next_round_dofs.end());
-
- // eliminate dofs which are
- // already numbered
- for (int s=next_round_dofs.size()-1; s>=0; --s)
- if (new_number[next_round_dofs[s]] != -1)
- next_round_dofs.erase (&next_round_dofs[s]);
-
- // check whether there are any new
- // dofs
- all_dofs_renumbered = (next_round_dofs.size() == 0);
- if (all_dofs_renumbered)
- // end loop if possible
- continue;
-
-
- // store for each coordination
- // number the dofs with these
- // coordination number
- multimap<unsigned int, int> dofs_by_coordination;
-
- // find coordination number for
- // each of these dofs
- for (vector<int>::iterator s=next_round_dofs.begin();
- s!=next_round_dofs.end(); ++s)
- {
- unsigned int coordination = 0;
- for (unsigned int j=sparsity.get_rowstart_indices()[*s];
- j<sparsity.get_rowstart_indices()[*s+1]; ++j)
- if (sparsity.get_column_numbers()[j] == -1)
- break;
- else
- ++coordination;
-
- // insert this dof at its
- // coordination number
- const pair<const unsigned int, int> new_entry (coordination, *s);
- dofs_by_coordination.insert (new_entry);
- };
-
- ////
- multimap<unsigned int, int>::iterator i;
- for (i = dofs_by_coordination.begin(); i!=dofs_by_coordination.end(); ++i)
- new_number[i->second] = next_free_number++;
-
- // after that: copy this round's
- // dofs for the next round
- last_round_dofs = next_round_dofs;
- };
-
-#ifdef DEBUG
- // test for all indices numbered
- if (find (new_number.begin(), new_number.end(), -1) != new_number.end())
- Assert (false, ExcRenumberingIncomplete());
- Assert (next_free_number == n_dofs,
- ExcRenumberingIncomplete());
-#endif
-
- switch (method)
- {
- case Cuthill_McKee:
- break;
- case reverse_Cuthill_McKee:
- {
- for (vector<int>::iterator i=new_number.begin(); i!=new_number.end(); ++i)
- *i = n_dofs-*i;
- break;
- };
- default:
- Assert (false, ExcNotImplemented());
- };
-
- // actually perform renumbering;
- // this is dimension specific and
- // thus needs an own function
- renumber_dofs (new_number);
-};
-
-
-
#if deal_II_dimension == 1
template <>
--- /dev/null
+/* $Id$ */
+
+#include <grid/dof.h>
+#include <grid/mg_dof.h>
+#include <grid/dof_constraints.h>
+#include <numerics/dof_renumbering.h>
+#include <lac/sparsematrix.h>
+
+#include <vector>
+#include <map>
+#include <algorithm>
+
+
+
+template <int dim>
+void DoFRenumbering::renumber_Cuthill_McKee (DoFHandler<dim> &dof_handler,
+ const bool reversed_numbering,
+ const bool use_constraints,
+ const vector<int> &starting_indices) {
+ // make the connection graph
+ SparseMatrixStruct sparsity (dof_handler.n_dofs(),
+ dof_handler.max_couplings_between_dofs());
+ dof_handler.make_sparsity_pattern (sparsity);
+
+ if (use_constraints)
+ {
+ ConstraintMatrix constraints;
+ dof_handler.make_hanging_node_constraints (constraints);
+ constraints.close ();
+ constraints.condense (sparsity);
+ };
+
+ const int n_dofs = sparsity.n_rows();
+ // store the new dof numbers; -1 means
+ // that no new number was chosen yet
+ vector<int> new_number(sparsity.n_rows(), -1);
+
+ // store the indices of the dofs renumbered
+ // in the last round. Default to starting
+ // points
+ vector<int> last_round_dofs (starting_indices);
+
+ // delete disallowed elements
+ for (unsigned int i=0; i<last_round_dofs.size(); ++i)
+ if ((last_round_dofs[i]<0) || (last_round_dofs[i]>=n_dofs))
+ last_round_dofs[i] = -1;
+
+ remove_if (last_round_dofs.begin(), last_round_dofs.end(),
+ bind2nd(equal_to<int>(), -1));
+
+ // now if no valid points remain:
+ // find dof with lowest coordination
+ // number
+
+ if (last_round_dofs.size() == 0)
+ {
+ int starting_point = -1;
+ unsigned int min_coordination = n_dofs;
+ for (int row=0; row<n_dofs; ++row)
+ {
+ unsigned int j;
+
+ // loop until we hit the end
+ // of this row's entries
+ for (j=sparsity.get_rowstart_indices()[row];
+ j<sparsity.get_rowstart_indices()[row+1]; ++j)
+ if (sparsity.get_column_numbers()[j] == -1)
+ break;
+ // post-condition after loop:
+ // coordination, i.e. the number
+ // of entries in this row is now
+ // j-rowstart[row]
+ if (j-sparsity.get_rowstart_indices()[row] < min_coordination)
+ {
+ min_coordination = j-sparsity.get_rowstart_indices()[row];
+ starting_point = row;
+ };
+ };
+
+ // now we still have to care for the
+ // case that no dof has a coordination
+ // number less than n_dofs. this rather
+ // exotic case only happens if we only
+ // have one cell, as far as I can see,
+ // but there may be others as well.
+ //
+ // if that should be the case, we can
+ // chose an arbitrary dof as starting
+ // point, e.g. the one with number zero
+ if (starting_point == -1)
+ starting_point = 0;
+
+ // initialize the first dof
+ last_round_dofs.push_back (starting_point);
+ };
+
+
+ // store next free dof index
+ int next_free_number = 0;
+
+ // enumerate the first round dofs
+ for (unsigned int i=0; i!=last_round_dofs.size(); ++i)
+ new_number[last_round_dofs[i]] = next_free_number++;
+
+ bool all_dofs_renumbered = false;
+
+ // now do as many steps as needed to
+ // renumber all dofs
+ while (!all_dofs_renumbered)
+ {
+ // store the indices of the dofs to be
+ // renumbered in the next round
+ vector<int> next_round_dofs;
+
+ // find all neighbors of the
+ // dofs numbered in the last
+ // round
+ for (unsigned int i=0; i<last_round_dofs.size(); ++i)
+ for (unsigned int j=sparsity.get_rowstart_indices()[last_round_dofs[i]];
+ j<sparsity.get_rowstart_indices()[last_round_dofs[i]+1]; ++j)
+ if (sparsity.get_column_numbers()[j] == -1)
+ break;
+ else
+ next_round_dofs.push_back (sparsity.get_column_numbers()[j]);
+
+ // sort dof numbers
+ sort (next_round_dofs.begin(), next_round_dofs.end());
+
+ // delete multiple entries
+ vector<int>::iterator end_sorted;
+ end_sorted = unique (next_round_dofs.begin(), next_round_dofs.end());
+ next_round_dofs.erase (end_sorted, next_round_dofs.end());
+
+ // eliminate dofs which are
+ // already numbered
+ for (int s=next_round_dofs.size()-1; s>=0; --s)
+ if (new_number[next_round_dofs[s]] != -1)
+ next_round_dofs.erase (&next_round_dofs[s]);
+
+ // check whether there are any new
+ // dofs
+ all_dofs_renumbered = (next_round_dofs.size() == 0);
+ if (all_dofs_renumbered)
+ // end loop if possible
+ continue;
+
+
+ // store for each coordination
+ // number the dofs with these
+ // coordination number
+ multimap<unsigned int, int> dofs_by_coordination;
+
+ // find coordination number for
+ // each of these dofs
+ for (vector<int>::iterator s=next_round_dofs.begin();
+ s!=next_round_dofs.end(); ++s)
+ {
+ unsigned int coordination = 0;
+ for (unsigned int j=sparsity.get_rowstart_indices()[*s];
+ j<sparsity.get_rowstart_indices()[*s+1]; ++j)
+ if (sparsity.get_column_numbers()[j] == -1)
+ break;
+ else
+ ++coordination;
+
+ // insert this dof at its
+ // coordination number
+ const pair<const unsigned int, int> new_entry (coordination, *s);
+ dofs_by_coordination.insert (new_entry);
+ };
+
+ ////
+ multimap<unsigned int, int>::iterator i;
+ for (i = dofs_by_coordination.begin(); i!=dofs_by_coordination.end(); ++i)
+ new_number[i->second] = next_free_number++;
+
+ // after that: copy this round's
+ // dofs for the next round
+ last_round_dofs = next_round_dofs;
+ };
+
+#ifdef DEBUG
+ // test for all indices numbered
+ if (find (new_number.begin(), new_number.end(), -1) != new_number.end())
+ Assert (false, ExcRenumberingIncomplete());
+ Assert (next_free_number == n_dofs,
+ ExcRenumberingIncomplete());
+#endif
+
+ if (reversed_numbering)
+ for (vector<int>::iterator i=new_number.begin(); i!=new_number.end(); ++i)
+ *i = n_dofs-*i;
+
+ // actually perform renumbering;
+ // this is dimension specific and
+ // thus needs an own function
+ dof_handler.renumber_dofs (new_number);
+};
+
+
+
+
+template <int dim>
+void DoFRenumbering::renumber_Cuthill_McKee (MGDoFHandler<dim> &dof_handler,
+ const unsigned int level,
+ const bool reversed_numbering,
+ const vector<int> &starting_indices) {
+ // make the connection graph
+ SparseMatrixStruct sparsity (dof_handler.n_dofs(level),
+ dof_handler.max_couplings_between_dofs());
+ dof_handler.make_sparsity_pattern (level, sparsity);
+
+ const int n_dofs = sparsity.n_rows();
+ // store the new dof numbers; -1 means
+ // that no new number was chosen yet
+ vector<int> new_number(n_dofs, -1);
+
+ // store the indices of the dofs renumbered
+ // in the last round. Default to starting
+ // points
+ vector<int> last_round_dofs (starting_indices);
+
+ // delete disallowed elements
+ for (unsigned int i=0; i<last_round_dofs.size(); ++i)
+ if ((last_round_dofs[i]<0) || (last_round_dofs[i]>=n_dofs))
+ last_round_dofs[i] = -1;
+
+ remove_if (last_round_dofs.begin(), last_round_dofs.end(),
+ bind2nd(equal_to<int>(), -1));
+
+ // now if no valid points remain:
+ // find dof with lowest coordination
+ // number
+
+ if (last_round_dofs.size() == 0)
+ {
+ int starting_point = -1;
+ unsigned int min_coordination = n_dofs;
+ for (int row=0; row<n_dofs; ++row)
+ {
+ unsigned int j;
+ for (j=sparsity.get_rowstart_indices()[row];
+ j<sparsity.get_rowstart_indices()[row+1]; ++j)
+ if (sparsity.get_column_numbers()[j] == -1)
+ break;
+ // post-condition after loop:
+ // coordination is now
+ // j-rowstart[row]
+ if (j-sparsity.get_rowstart_indices()[row] < min_coordination)
+ {
+ min_coordination = j-sparsity.get_rowstart_indices()[row];
+ starting_point = row;
+ };
+ };
+ // initialize the first dof
+ last_round_dofs.push_back (starting_point);
+ };
+
+
+ // store next free dof index
+ int next_free_number = 0;
+
+ // enumerate the first round dofs
+ for (unsigned int i=0; i!=last_round_dofs.size(); ++i)
+ new_number[last_round_dofs[i]] = next_free_number++;
+
+ bool all_dofs_renumbered = false;
+
+ // now do as many steps as needed to
+ // renumber all dofs
+ while (!all_dofs_renumbered)
+ {
+ // store the indices of the dofs to be
+ // renumbered in the next round
+ vector<int> next_round_dofs;
+
+ // find all neighbors of the
+ // dofs numbered in the last
+ // round
+ for (unsigned int i=0; i<last_round_dofs.size(); ++i)
+ for (unsigned int j=sparsity.get_rowstart_indices()[last_round_dofs[i]];
+ j<sparsity.get_rowstart_indices()[last_round_dofs[i]+1]; ++j)
+ if (sparsity.get_column_numbers()[j] == -1)
+ break;
+ else
+ next_round_dofs.push_back (sparsity.get_column_numbers()[j]);
+
+ // sort dof numbers
+ sort (next_round_dofs.begin(), next_round_dofs.end());
+
+ // delete multiple entries
+ vector<int>::iterator end_sorted;
+ end_sorted = unique (next_round_dofs.begin(), next_round_dofs.end());
+ next_round_dofs.erase (end_sorted, next_round_dofs.end());
+
+ // eliminate dofs which are
+ // already numbered
+ for (int s=next_round_dofs.size()-1; s>=0; --s)
+ if (new_number[next_round_dofs[s]] != -1)
+ next_round_dofs.erase (&next_round_dofs[s]);
+
+ // check whether there are any new
+ // dofs
+ all_dofs_renumbered = (next_round_dofs.size() == 0);
+ if (all_dofs_renumbered)
+ // end loop if possible
+ continue;
+
+
+ // store for each coordination
+ // number the dofs with these
+ // coordination number
+ multimap<unsigned int, int> dofs_by_coordination;
+
+ // find coordination number for
+ // each of these dofs
+ for (vector<int>::iterator s=next_round_dofs.begin();
+ s!=next_round_dofs.end(); ++s)
+ {
+ unsigned int coordination = 0;
+ for (unsigned int j=sparsity.get_rowstart_indices()[*s];
+ j<sparsity.get_rowstart_indices()[*s+1]; ++j)
+ if (sparsity.get_column_numbers()[j] == -1)
+ break;
+ else
+ ++coordination;
+
+ // insert this dof at its
+ // coordination number
+ const pair<const unsigned int, int> new_entry (coordination, *s);
+ dofs_by_coordination.insert (new_entry);
+ };
+
+ ////
+ multimap<unsigned int, int>::iterator i;
+ for (i = dofs_by_coordination.begin(); i!=dofs_by_coordination.end(); ++i)
+ new_number[i->second] = next_free_number++;
+
+ // after that: copy this round's
+ // dofs for the next round
+ last_round_dofs = next_round_dofs;
+ };
+
+#ifdef DEBUG
+ // test for all indices numbered
+ if (find (new_number.begin(), new_number.end(), -1) != new_number.end())
+ Assert (false, ExcRenumberingIncomplete());
+ Assert (next_free_number == n_dofs,
+ ExcRenumberingIncomplete());
+#endif
+
+ if (reversed_numbering)
+ for (vector<int>::iterator i=new_number.begin(); i!=new_number.end(); ++i)
+ *i = n_dofs-*i;
+
+ // actually perform renumbering;
+ // this is dimension specific and
+ // thus needs an own function
+ dof_handler.renumber_dofs (level, new_number);
+};
+
+
+
+
+
+
+// explicit instantiations
+template
+void DoFRenumbering::renumber_Cuthill_McKee (DoFHandler<deal_II_dimension> &dof_handler,
+ const bool reversed_numbering,
+ const bool use_constraints,
+ const vector<int> &starting_indices);
+
+template
+void DoFRenumbering::renumber_Cuthill_McKee (MGDoFHandler<deal_II_dimension> &dof_handler,
+ const unsigned int level,
+ const bool reversed_numbering,
+ const vector<int> &starting_indices);
+
+/* ------------------------ MGDoFHexAccessor --------------------------- */
+
+template <int dim, typename BaseClass>
+MGDoFHexAccessor<dim,BaseClass>::MGDoFHexAccessor (Triangulation<dim> *tria,
+ const int level,
+ const int index,
+ const AccessorData *local_data) :
+ MGDoFAccessor<dim> (local_data),
+ DoFHexAccessor(tria,level,index,local_data) {};
+
+
+
+template <int dim, typename BaseClass>
+inline
+int MGDoFHexAccessor<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_fe() != 0, ExcInvalidObject());
+ Assert (i<dof_handler->get_fe().dofs_per_hex,
+ ExcInvalidIndex (i, 0, dof_handler->get_fe().dofs_per_hex));
+
+ return mg_dof_handler->mg_levels[present_level]
+ ->hex_dofs[present_index*dof_handler->get_fe().dofs_per_hex+i];
+};
+
+
+
+template <int dim, typename BaseClass>
+void MGDoFHexAccessor<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_fe() != 0, ExcInvalidObject());
+ Assert (i<dof_handler->get_fe().dofs_per_hex,
+ ExcInvalidIndex (i, 0, dof_handler->get_fe().dofs_per_hex));
+
+ mg_dof_handler->mg_levels[present_level]
+ ->hex_dofs[present_index*dof_handler->get_fe().dofs_per_hex+i] = index;
+};
+
+
+
+template <int dim, typename BaseClass>
+inline
+int MGDoFHexAccessor<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_fe() != 0, ExcInvalidObject());
+ Assert (vertex<8, ExcInvalidIndex (i,0,8));
+ Assert (i<dof_handler->get_fe().dofs_per_vertex,
+ ExcInvalidIndex (i, 0, dof_handler->get_fe().dofs_per_vertex));
+
+ return (mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)]
+ .get_index (present_level, i, dof_handler->get_fe().dofs_per_vertex));
+};
+
+
+
+template <int dim, typename BaseClass>
+void MGDoFHexAccessor<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_fe() != 0, ExcInvalidObject());
+ Assert (vertex<4, ExcInvalidIndex (i,0,4));
+ Assert (i<dof_handler->get_fe().dofs_per_vertex,
+ ExcInvalidIndex (i, 0, dof_handler->get_fe().dofs_per_vertex));
+
+ mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)]
+ .set_index (present_level, i, dof_handler->get_fe().dofs_per_vertex, index);
+};
+
+
+
+template <int dim, typename BaseClass>
+void
+MGDoFHexAccessor<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_fe() != 0, ExcInvalidObject());
+ Assert (dof_indices.size() == (8*dof_handler->get_fe().dofs_per_vertex +
+ 12*dof_handler->get_fe().dofs_per_line +
+ 6*dof_handler->get_fe().dofs_per_quad +
+ dof_handler->get_fe().dofs_per_hex),
+ ExcVectorDoesNotMatch());
+
+ const unsigned int dofs_per_vertex = dof_handler->get_fe().dofs_per_vertex,
+ dofs_per_line = dof_handler->get_fe().dofs_per_line,
+ dofs_per_quad = dof_handler->get_fe().dofs_per_quad,
+ dofs_per_hex = dof_handler->get_fe().dofs_per_hex;
+ vector<int>::iterator next = dof_indices.begin();
+ for (unsigned int vertex=0; vertex<8; ++vertex)
+ for (unsigned int d=0; d<dofs_per_vertex; ++d)
+ *next++ = mg_vertex_dof_index(vertex,d);
+ for (unsigned int line=0; line<12; ++line)
+ for (unsigned int d=0; d<dofs_per_line; ++d)
+ *next++ = this->line(line)->mg_dof_index(d);
+ for (unsigned int quad=0; quad<12; ++quad)
+ for (unsigned int d=0; d<dofs_per_quad; ++d)
+ *next++ = this->quad(quad)->mg_dof_index(d);
+ for (unsigned int d=0; d<dofs_per_hex; ++d)
+ *next++ = mg_dof_index(d);
+};
+
+
+
+template <int dim, typename BaseClass>
+TriaIterator<dim,MGDoFLineAccessor<dim,LineAccessor<dim> > >
+MGDoFHexAccessor<dim,BaseClass>::line (const unsigned int i) const {
+ Assert (i<12, ExcInvalidIndex (i, 0, 12));
+
+ 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,QuadAccessor<dim> > >
+MGDoFHexAccessor<dim,BaseClass>::quad (const unsigned int i) const {
+ Assert (i<12, ExcInvalidIndex (i, 0, 6));
+
+ return TriaIterator<dim,MGDoFQuadAccessor<dim,QuadAccessor<dim> > >
+ (
+ tria,
+ present_level,
+ quad_index (i),
+ mg_dof_handler
+ );
+};
+
+
+
+template <int dim, typename BaseClass>
+TriaIterator<dim,MGDoFHexAccessor<dim,BaseClass> >
+MGDoFHexAccessor<dim,BaseClass>::child (const unsigned int i) const {
+ TriaIterator<dim,MGDoFHexAccessor<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;
+};
+
+
+
+template <int dim, typename BaseClass>
+void
+MGDoFHexAccessor<dim,BaseClass>::copy_from (const MGDoFHexAccessor<dim,BaseClass> &a) {
+ DoFHexAccessor::copy_from (a);
+ set_mg_dof_handler (a.mg_dof_handler);
+};
+
+
+
+
/*------------------------- Functions: MGDoFCellAccessor -----------------------*/
#endif
+#if deal_II_dimension == 3
+template class MGDoFLineAccessor<3,LineAccessor<3> >;
+template class MGDoFQuadAccessor<3,QuadAccessor<3> >;
+template class MGDoFHexAccessor<3,HexAccessor<3> >;
+template class MGDoFHexAccessor<3,CellAccessor<3> >;
+template class MGDoFCellAccessor<3>;
+
+template class TriaRawIterator<3,MGDoFLineAccessor<3,LineAccessor<3> > >;
+template class TriaRawIterator<3,MGDoFQuadAccessor<3,QuadAccessor<3> > >;
+template class TriaRawIterator<3,MGDoFHexAccessor<3,HexAccessor<3> > >;
+template class TriaRawIterator<3,MGDoFCellAccessor<3> >;
+template class TriaIterator<3,MGDoFLineAccessor<3,LineAccessor<3> > >;
+template class TriaIterator<3,MGDoFQuadAccessor<3,QuadAccessor<3> > >;
+template class TriaIterator<3,MGDoFCellAccessor<3> >;
+template class TriaActiveIterator<3,MGDoFLineAccessor<3,LineAccessor<3> > >;
+template class TriaActiveIterator<3,MGDoFQuadAccessor<3,QuadAccessor<3> > >;
+template class TriaActiveIterator<3,MGDoFCellAccessor<3> >;
+#endif
+
+
return 0;
};
+
+template <>
+MGDoFHandler<1>::raw_hex_iterator
+MGDoFHandler<1>::begin_raw_hex (const unsigned int) const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::hex_iterator
+MGDoFHandler<1>::begin_hex (const unsigned int) const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::active_hex_iterator
+MGDoFHandler<1>::begin_active_hex (const unsigned int) const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::raw_hex_iterator
+MGDoFHandler<1>::end_hex () const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::raw_hex_iterator
+MGDoFHandler<1>::last_raw_hex (const unsigned int) const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+template <>
+MGDoFHandler<1>::hex_iterator
+MGDoFHandler<1>::last_hex (const unsigned int) const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::active_hex_iterator
+MGDoFHandler<1>::last_active_hex (const unsigned int) const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::raw_hex_iterator
+MGDoFHandler<1>::last_raw_hex () const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::hex_iterator
+MGDoFHandler<1>::last_hex () const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::active_hex_iterator
+MGDoFHandler<1>::last_active_hex () const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
#endif
return last_active_line (level);
};
+
+
+template <>
+MGDoFHandler<2>::raw_hex_iterator
+MGDoFHandler<2>::begin_raw_hex (const unsigned int) const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::hex_iterator
+MGDoFHandler<2>::begin_hex (const unsigned int) const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::active_hex_iterator
+MGDoFHandler<2>::begin_active_hex (const unsigned int) const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::raw_hex_iterator
+MGDoFHandler<2>::end_hex () const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::raw_hex_iterator
+MGDoFHandler<2>::last_raw_hex (const unsigned int) const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+template <>
+MGDoFHandler<2>::hex_iterator
+MGDoFHandler<2>::last_hex (const unsigned int) const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::active_hex_iterator
+MGDoFHandler<2>::last_active_hex (const unsigned int) const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::raw_hex_iterator
+MGDoFHandler<2>::last_raw_hex () const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::hex_iterator
+MGDoFHandler<2>::last_hex () const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::active_hex_iterator
+MGDoFHandler<2>::last_active_hex () const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+#endif
+
+
+
+
+#if deal_II_dimension == 3
+
+template <>
+MGDoFHandler<3>::raw_cell_iterator
+MGDoFHandler<3>::begin_raw (const unsigned int level) const {
+ return begin_raw_hex (level);
+};
+
+
+
+template <>
+MGDoFHandler<3>::cell_iterator
+MGDoFHandler<3>::begin (const unsigned int level) const {
+ return begin_hex (level);
+};
+
+
+
+template <>
+MGDoFHandler<3>::active_cell_iterator
+MGDoFHandler<3>::begin_active (const unsigned int level) const {
+ return begin_active_hex (level);
+};
+
+
+
+template <>
+MGDoFHandler<3>::raw_cell_iterator
+MGDoFHandler<3>::end () const {
+ return end_hex ();
+};
+
+
+
+template <>
+MGDoFHandler<3>::raw_cell_iterator
+MGDoFHandler<3>::last_raw () const {
+ return last_raw_hex ();
+};
+
+
+
+template <>
+MGDoFHandler<3>::raw_cell_iterator
+MGDoFHandler<3>::last_raw (const unsigned int level) const {
+ return last_raw_hex (level);
+};
+
+
+
+template <>
+MGDoFHandler<3>::cell_iterator
+MGDoFHandler<3>::last () const {
+ return last_hex ();
+};
+
+
+
+template <>
+MGDoFHandler<3>::cell_iterator
+MGDoFHandler<3>::last (const unsigned int level) const {
+ return last_hex (level);
+};
+
+
+
+template <>
+MGDoFHandler<3>::active_cell_iterator
+MGDoFHandler<3>::last_active () const {
+ return last_active_hex ();
+};
+
+
+
+template <>
+MGDoFHandler<3>::active_cell_iterator
+MGDoFHandler<3>::last_active (const unsigned int level) const {
+ return last_active_hex (level);
+};
+
+
+template <>
+MGDoFDimensionInfo<3>::raw_face_iterator
+MGDoFHandler<3>::begin_raw_face (const unsigned int level) const {
+ return begin_raw_quad (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::face_iterator
+MGDoFHandler<3>::begin_face (const unsigned int level) const {
+ return begin_quad (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::active_face_iterator
+MGDoFHandler<3>::begin_active_face (const unsigned int level) const {
+ return begin_active_quad (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::raw_face_iterator
+MGDoFHandler<3>::end_face () const {
+ return end_quad ();
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::raw_face_iterator
+MGDoFHandler<3>::last_raw_face () const {
+ return last_raw_quad ();
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::raw_face_iterator
+MGDoFHandler<3>::last_raw_face (const unsigned int level) const {
+ return last_raw_quad (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::face_iterator
+MGDoFHandler<3>::last_face () const {
+ return last_quad ();
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::face_iterator
+MGDoFHandler<3>::last_face (const unsigned int level) const {
+ return last_quad (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::active_face_iterator
+MGDoFHandler<3>::last_active_face () const {
+ return last_active_quad ();
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::active_face_iterator
+MGDoFHandler<3>::last_active_face (const unsigned int level) const {
+ return last_active_quad (level);
+};
+
#endif
+template <int dim>
+typename MGDoFHandler<dim>::raw_hex_iterator
+MGDoFHandler<dim>::begin_raw_hex (const unsigned int level) const {
+ return raw_hex_iterator (tria,
+ tria->begin_raw_hex(level)->level(),
+ tria->begin_raw_hex(level)->index(),
+ this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::hex_iterator
+MGDoFHandler<dim>::begin_hex (const unsigned int level) const {
+ return hex_iterator (tria,
+ tria->begin_hex(level)->level(),
+ tria->begin_hex(level)->index(),
+ this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::active_hex_iterator
+MGDoFHandler<dim>::begin_active_hex (const unsigned int level) const {
+ return active_hex_iterator (tria,
+ tria->begin_active_hex(level)->level(),
+ tria->begin_active_hex(level)->index(),
+ this);
+};
+
+
+
template <int dim>
typename MGDoFHandler<dim>::raw_line_iterator
MGDoFHandler<dim>::end_line () const {
+template <int dim>
+typename MGDoFHandler<dim>::raw_hex_iterator
+MGDoFHandler<dim>::end_hex () const {
+ return raw_hex_iterator (tria, -1, -1, this);
+};
+
+
+
template <int dim>
typename MGDoFDimensionInfo<dim>::raw_cell_iterator
MGDoFHandler<dim>::end_raw (const unsigned int level) const {
-
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::raw_quad_iterator
+MGDoFHandler<dim>::end_raw_quad (const unsigned int level) const {
+ return (level == mg_levels.size()-1 ?
+ end_quad() :
+ begin_raw_quad (level+1));
+};
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::quad_iterator
+MGDoFHandler<dim>::end_quad (const unsigned int level) const {
+ return (level == mg_levels.size()-1 ?
+ quad_iterator(end_quad()) :
+ begin_quad (level+1));
+};
+
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::active_quad_iterator
+MGDoFHandler<dim>::end_active_quad (const unsigned int level) const {
+ return (level == mg_levels.size()-1 ?
+ active_quad_iterator(end_quad()) :
+ begin_active_quad (level+1));
+};
+
+
+
template <int dim>
-typename MGDoFDimensionInfo<dim>::raw_quad_iterator
-MGDoFHandler<dim>::end_raw_quad (const unsigned int level) const {
+typename MGDoFDimensionInfo<dim>::raw_hex_iterator
+MGDoFHandler<dim>::end_raw_hex (const unsigned int level) const {
return (level == mg_levels.size()-1 ?
- end_quad() :
- begin_raw_quad (level+1));
+ end_hex() :
+ begin_raw_hex (level+1));
};
template <int dim>
-typename MGDoFDimensionInfo<dim>::quad_iterator
-MGDoFHandler<dim>::end_quad (const unsigned int level) const {
+typename MGDoFDimensionInfo<dim>::hex_iterator
+MGDoFHandler<dim>::end_hex (const unsigned int level) const {
return (level == mg_levels.size()-1 ?
- quad_iterator(end_quad()) :
- begin_quad (level+1));
+ hex_iterator(end_hex()) :
+ begin_hex (level+1));
};
template <int dim>
-typename MGDoFDimensionInfo<dim>::active_quad_iterator
-MGDoFHandler<dim>::end_active_quad (const unsigned int level) const {
+typename MGDoFDimensionInfo<dim>::active_hex_iterator
+MGDoFHandler<dim>::end_active_hex (const unsigned int level) const {
return (level == mg_levels.size()-1 ?
- active_quad_iterator(end_quad()) :
- begin_active_quad (level+1));
+ active_hex_iterator(end_hex()) :
+ begin_active_hex (level+1));
};
+template <int dim>
+typename MGDoFHandler<dim>::raw_hex_iterator
+MGDoFHandler<dim>::last_raw_hex (const unsigned int level) const {
+ return raw_hex_iterator (tria,
+ tria->last_raw_hex(level)->level(),
+ tria->last_raw_hex(level)->index(),
+ this);
+};
+
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::hex_iterator
+MGDoFHandler<dim>::last_hex (const unsigned int level) const {
+ return hex_iterator (tria,
+ tria->last_hex(level)->level(),
+ tria->last_hex(level)->index(),
+ this);
+};
+
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::active_hex_iterator
+MGDoFHandler<dim>::last_active_hex (const unsigned int level) const {
+ return active_hex_iterator (tria,
+ tria->last_active_hex(level)->level(),
+ tria->last_active_hex(level)->index(),
+ this);
+};
+
+
+
+
template <int dim>
typename MGDoFHandler<dim>::raw_line_iterator
MGDoFHandler<dim>::last_raw_line () const {
+template <int dim>
+typename MGDoFHandler<dim>::raw_hex_iterator
+MGDoFHandler<dim>::last_raw_hex () const {
+ return last_raw_hex (mg_levels.size()-1);
+};
+
+
+
template <int dim>
typename MGDoFHandler<dim>::line_iterator
MGDoFHandler<dim>::last_line () const {
+template <int dim>
+typename MGDoFHandler<dim>::hex_iterator
+MGDoFHandler<dim>::last_hex () const {
+ return last_hex (mg_levels.size()-1);
+};
+
+
+
template <int dim>
typename MGDoFHandler<dim>::active_line_iterator
MGDoFHandler<dim>::last_active_line () const {
+template <int dim>
+typename MGDoFHandler<dim>::active_hex_iterator
+MGDoFHandler<dim>::last_active_hex () const {
+ return last_active_hex (mg_levels.size()-1);
+};
+
+
+
+#if deal_II_dimension == 3
+
+template <>
+unsigned int
+MGDoFHandler<3>::distribute_dofs_on_cell (cell_iterator &cell,
+ unsigned int next_free_dof) {
+ if (selected_fe->dofs_per_vertex > 0)
+ // number dofs on vertices
+ for (unsigned int vertex=0; vertex<GeometryInfo<3>::vertices_per_cell; ++vertex)
+ // check whether dofs for this
+ // vertex have been distributed
+ // (only check the first dof)
+ if (cell->mg_vertex_dof_index(vertex, 0) == -1)
+ for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
+ cell->set_mg_vertex_dof_index (vertex, d, next_free_dof++);
+
+ // for the lines
+ if (selected_fe->dofs_per_line > 0)
+ for (unsigned int l=0; l<GeometryInfo<3>::lines_per_cell; ++l)
+ {
+ line_iterator line = cell->line(l);
+
+ // distribute dofs if necessary:
+ // check whether line dof is already
+ // numbered (check only first dof)
+ if (line->mg_dof_index(0) == -1)
+ // if not: distribute dofs
+ for (unsigned int d=0; d<selected_fe->dofs_per_line; ++d)
+ line->set_mg_dof_index (d, next_free_dof++);
+ };
+
+ // for the quads
+ if (selected_fe->dofs_per_quad > 0)
+ for (unsigned int q=0; q<GeometryInfo<3>::lines_per_cell; ++q)
+ {
+ quad_iterator quad = cell->quad(q);
+
+ // distribute dofs if necessary:
+ // check whether line dof is already
+ // numbered (check only first dof)
+ if (quad->mg_dof_index(0) == -1)
+ // if not: distribute dofs
+ for (unsigned int d=0; d<selected_fe->dofs_per_line; ++d)
+ quad->set_mg_dof_index (d, next_free_dof++);
+ };
+
+
+ // dofs of cell
+ if (selected_fe->dofs_per_hex > 0)
+ for (unsigned int d=0; d<selected_fe->dofs_per_hex; ++d)
+ cell->set_mg_dof_index (d, next_free_dof++);
+
+
+ // note that this cell has been processed
+ cell->set_user_flag ();
+
+ return next_free_dof;
+};
+
+#endif
+
+
+
template <int dim>
unsigned int MGDoFHandler<dim>::n_dofs (const unsigned int level) const {
-
-template <int dim>
-void MGDoFHandler<dim>::renumber_dofs (const unsigned int level,
- const RenumberingMethod method,
- const vector<int> &starting_points) {
- // make the connection graph
- SparseMatrixStruct sparsity (n_dofs(level), max_couplings_between_dofs());
- make_sparsity_pattern (level, sparsity);
-
- int n_dofs = sparsity.n_rows();
- // store the new dof numbers; -1 means
- // that no new number was chosen yet
- vector<int> new_number(n_dofs, -1);
-
- // store the indices of the dofs renumbered
- // in the last round. Default to starting
- // points
- vector<int> last_round_dofs (starting_points);
-
- // delete disallowed elements
- for (unsigned int i=0; i<last_round_dofs.size(); ++i)
- if ((last_round_dofs[i]<0) || (last_round_dofs[i]>=n_dofs))
- last_round_dofs[i] = -1;
-
- remove_if (last_round_dofs.begin(), last_round_dofs.end(),
- bind2nd(equal_to<int>(), -1));
-
- // now if no valid points remain:
- // find dof with lowest coordination
- // number
-
- if (last_round_dofs.size() == 0)
- {
- int starting_point = -1;
- unsigned int min_coordination = n_dofs;
- for (int row=0; row<n_dofs; ++row)
- {
- unsigned int j;
- for (j=sparsity.get_rowstart_indices()[row];
- j<sparsity.get_rowstart_indices()[row+1]; ++j)
- if (sparsity.get_column_numbers()[j] == -1)
- break;
- // post-condition after loop:
- // coordination is now
- // j-rowstart[row]
- if (j-sparsity.get_rowstart_indices()[row] < min_coordination)
- {
- min_coordination = j-sparsity.get_rowstart_indices()[row];
- starting_point = row;
- };
- };
- // initialize the first dof
- last_round_dofs.push_back (starting_point);
- };
-
-
- // store next free dof index
- int next_free_number = 0;
-
- // enumerate the first round dofs
- for (unsigned int i=0; i!=last_round_dofs.size(); ++i)
- new_number[last_round_dofs[i]] = next_free_number++;
-
- bool all_dofs_renumbered = false;
-
- // now do as many steps as needed to
- // renumber all dofs
- while (!all_dofs_renumbered)
- {
- // store the indices of the dofs to be
- // renumbered in the next round
- vector<int> next_round_dofs;
-
- // find all neighbors of the
- // dofs numbered in the last
- // round
- for (unsigned int i=0; i<last_round_dofs.size(); ++i)
- for (unsigned int j=sparsity.get_rowstart_indices()[last_round_dofs[i]];
- j<sparsity.get_rowstart_indices()[last_round_dofs[i]+1]; ++j)
- if (sparsity.get_column_numbers()[j] == -1)
- break;
- else
- next_round_dofs.push_back (sparsity.get_column_numbers()[j]);
-
- // sort dof numbers
- sort (next_round_dofs.begin(), next_round_dofs.end());
-
- // delete multiple entries
- vector<int>::iterator end_sorted;
- end_sorted = unique (next_round_dofs.begin(), next_round_dofs.end());
- next_round_dofs.erase (end_sorted, next_round_dofs.end());
-
- // eliminate dofs which are
- // already numbered
- for (int s=next_round_dofs.size()-1; s>=0; --s)
- if (new_number[next_round_dofs[s]] != -1)
- next_round_dofs.erase (&next_round_dofs[s]);
-
- // check whether there are any new
- // dofs
- all_dofs_renumbered = (next_round_dofs.size() == 0);
- if (all_dofs_renumbered)
- // end loop if possible
- continue;
-
-
- // store for each coordination
- // number the dofs with these
- // coordination number
- multimap<unsigned int, int> dofs_by_coordination;
-
- // find coordination number for
- // each of these dofs
- for (vector<int>::iterator s=next_round_dofs.begin();
- s!=next_round_dofs.end(); ++s)
- {
- unsigned int coordination = 0;
- for (unsigned int j=sparsity.get_rowstart_indices()[*s];
- j<sparsity.get_rowstart_indices()[*s+1]; ++j)
- if (sparsity.get_column_numbers()[j] == -1)
- break;
- else
- ++coordination;
-
- // insert this dof at its
- // coordination number
- const pair<const unsigned int, int> new_entry (coordination, *s);
- dofs_by_coordination.insert (new_entry);
- };
-
- ////
- multimap<unsigned int, int>::iterator i;
- for (i = dofs_by_coordination.begin(); i!=dofs_by_coordination.end(); ++i)
- new_number[i->second] = next_free_number++;
-
- // after that: copy this round's
- // dofs for the next round
- last_round_dofs = next_round_dofs;
- };
-
-#ifdef DEBUG
- // test for all indices numbered
- if (find (new_number.begin(), new_number.end(), -1) != new_number.end())
- Assert (false, ExcRenumberingIncomplete());
- Assert (next_free_number == n_dofs,
- ExcRenumberingIncomplete());
-#endif
-
- switch (method)
- {
- case Cuthill_McKee:
- break;
- case reverse_Cuthill_McKee:
- {
- for (vector<int>::iterator i=new_number.begin(); i!=new_number.end(); ++i)
- *i = n_dofs-*i;
- break;
- };
- default:
- Assert (false, ExcNotImplemented());
- };
-
- // actually perform renumbering;
- // this is dimension specific and
- // thus needs an own function
- do_renumbering (level, new_number);
-};
-
-
-
-
#if deal_II_dimension == 1
template <>
-void MGDoFHandler<1>::do_renumbering (const unsigned int level,
- const vector<int> &new_numbers) {
+void MGDoFHandler<1>::renumber_dofs (const unsigned int level,
+ const vector<int> &new_numbers) {
Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete());
// note that we can not use cell iterators
#if deal_II_dimension == 2
template <>
-void MGDoFHandler<2>::do_renumbering (const unsigned int level,
- const vector<int> &new_numbers) {
+void MGDoFHandler<2>::renumber_dofs (const unsigned int level,
+ const vector<int> &new_numbers) {
+ Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete());
+
+ for (vector<MGVertexDoFs>::iterator i=mg_vertex_dofs.begin();
+ i!=mg_vertex_dofs.end(); ++i)
+ // if the present vertex lives on
+ // the present level
+ if ((i->get_coarsest_level() <= level) &&
+ (i->get_finest_level() >= level))
+ for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
+ i->set_index (level, d, selected_fe->dofs_per_vertex,
+ new_numbers[i->get_index (level, d,
+ selected_fe->dofs_per_vertex)]);
+
+ for (vector<int>::iterator i=mg_levels[level]->line_dofs.begin();
+ i!=mg_levels[level]->line_dofs.end(); ++i)
+ {
+ Assert (*i != -1, ExcInternalError());
+ *i = new_numbers[*i];
+ };
+
+ for (vector<int>::iterator i=mg_levels[level]->quad_dofs.begin();
+ i!=mg_levels[level]->quad_dofs.end(); ++i)
+ {
+ Assert (*i != -1, ExcInternalError());
+ *i = new_numbers[*i];
+ };
+};
+
+#endif
+
+
+
+#if deal_II_dimension == 3
+
+template <>
+void MGDoFHandler<3>::renumber_dofs (const unsigned int level,
+ const vector<int> &new_numbers) {
Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete());
for (vector<MGVertexDoFs>::iterator i=mg_vertex_dofs.begin();
Assert (*i != -1, ExcInternalError());
*i = new_numbers[*i];
};
+
for (vector<int>::iterator i=mg_levels[level]->quad_dofs.begin();
i!=mg_levels[level]->quad_dofs.end(); ++i)
{
Assert (*i != -1, ExcInternalError());
*i = new_numbers[*i];
};
+
+ for (vector<int>::iterator i=mg_levels[level]->hex_dofs.begin();
+ i!=mg_levels[level]->hex_dofs.end(); ++i)
+ {
+ Assert (*i != -1, ExcInternalError());
+ *i = new_numbers[*i];
+ };
};
#endif
#endif
+
+
+#if deal_II_dimension == 3
+
+template <>
+void MGDoFHandler<3>::reserve_space () {
+ const unsigned int dim = 3;
+
+ Assert (selected_fe != 0, ExcNoFESelected());
+ Assert (tria->n_levels() > 0, ExcInvalidTriangulation());
+
+ ////////////////////////////
+ // DESTRUCTION
+
+ // delete all levels and set them up
+ // newly, since vectors are
+ // troublesome if you want to change
+ // their size
+ for (unsigned int i=0; i<mg_levels.size(); ++i)
+ delete mg_levels[i];
+ mg_levels.clear ();
+
+ // also delete vector of vertex indices
+ // this calls the destructor which
+ // must free the space
+ mg_vertex_dofs.clear ();
+
+
+ ////////////////////////////
+ // CONSTRUCTION
+
+ // first allocate space for the
+ // lines and quads on each level
+ for (unsigned int i=0; i<tria->n_levels(); ++i)
+ {
+ mg_levels.push_back (new DoFLevel<3>);
+
+ mg_levels.back()->line_dofs = vector<int> (tria->levels[i]->lines.lines.size() *
+ selected_fe->dofs_per_line,
+ -1);
+ mg_levels.back()->quad_dofs = vector<int> (tria->levels[i]->quads.quads.size() *
+ selected_fe->dofs_per_quad,
+ -1);
+ mg_levels.back()->hex_dofs = vector<int> (tria->levels[i]->hexes.hexes.size() *
+ selected_fe->dofs_per_hex,
+ -1);
+ };
+
+
+ // now allocate space for the
+ // vertices. To this end, we need
+ // to construct as many objects as
+ // there are vertices and let them
+ // allocate enough space for their
+ // vertex indices on the levels they
+ // live on. We need therefore to
+ // count to how many levels a cell
+ // belongs to, which we do by looping
+ // over all cells and storing the
+ // maximum and minimum level each
+ // vertex we pass by belongs to
+ mg_vertex_dofs.resize (tria->vertices.size());
+
+ vector<unsigned int> min_level (tria->vertices.size(), tria->n_levels());
+ vector<unsigned int> max_level (tria->vertices.size(), 0);
+
+ Triangulation<dim>::cell_iterator cell = tria->begin(),
+ endc = tria->end();
+ for (; cell!=endc; ++cell)
+ for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
+ ++vertex)
+ {
+ const unsigned int vertex_index = cell->vertex_index(vertex);
+ if (min_level[vertex_index] > static_cast<unsigned int>(cell->level()))
+ min_level[vertex_index] = cell->level();
+ if (max_level[vertex_index] < static_cast<unsigned int>(cell->level()))
+ max_level[vertex_index] = cell->level();
+ };
+
+
+ // now allocate the needed space
+ for (unsigned int vertex=0; vertex<tria->vertices.size(); ++vertex)
+ {
+ Assert (min_level[vertex] < tria->n_levels(), ExcInternalError());
+ Assert (max_level[vertex] >= min_level[vertex], ExcInternalError());
+
+ mg_vertex_dofs[vertex].init (min_level[vertex],
+ max_level[vertex],
+ selected_fe->dofs_per_vertex);
+ };
+};
+
+#endif
+
+
// explicite instantiations
template class MGDoFHandler<deal_II_dimension>;
+/* ------------------------ MGDoFHexAccessor --------------------------- */
+
+template <int dim, typename BaseClass>
+MGDoFHexAccessor<dim,BaseClass>::MGDoFHexAccessor (Triangulation<dim> *tria,
+ const int level,
+ const int index,
+ const AccessorData *local_data) :
+ MGDoFAccessor<dim> (local_data),
+ DoFHexAccessor(tria,level,index,local_data) {};
+
+
+
+template <int dim, typename BaseClass>
+inline
+int MGDoFHexAccessor<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_fe() != 0, ExcInvalidObject());
+ Assert (i<dof_handler->get_fe().dofs_per_hex,
+ ExcInvalidIndex (i, 0, dof_handler->get_fe().dofs_per_hex));
+
+ return mg_dof_handler->mg_levels[present_level]
+ ->hex_dofs[present_index*dof_handler->get_fe().dofs_per_hex+i];
+};
+
+
+
+template <int dim, typename BaseClass>
+void MGDoFHexAccessor<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_fe() != 0, ExcInvalidObject());
+ Assert (i<dof_handler->get_fe().dofs_per_hex,
+ ExcInvalidIndex (i, 0, dof_handler->get_fe().dofs_per_hex));
+
+ mg_dof_handler->mg_levels[present_level]
+ ->hex_dofs[present_index*dof_handler->get_fe().dofs_per_hex+i] = index;
+};
+
+
+
+template <int dim, typename BaseClass>
+inline
+int MGDoFHexAccessor<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_fe() != 0, ExcInvalidObject());
+ Assert (vertex<8, ExcInvalidIndex (i,0,8));
+ Assert (i<dof_handler->get_fe().dofs_per_vertex,
+ ExcInvalidIndex (i, 0, dof_handler->get_fe().dofs_per_vertex));
+
+ return (mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)]
+ .get_index (present_level, i, dof_handler->get_fe().dofs_per_vertex));
+};
+
+
+
+template <int dim, typename BaseClass>
+void MGDoFHexAccessor<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_fe() != 0, ExcInvalidObject());
+ Assert (vertex<4, ExcInvalidIndex (i,0,4));
+ Assert (i<dof_handler->get_fe().dofs_per_vertex,
+ ExcInvalidIndex (i, 0, dof_handler->get_fe().dofs_per_vertex));
+
+ mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)]
+ .set_index (present_level, i, dof_handler->get_fe().dofs_per_vertex, index);
+};
+
+
+
+template <int dim, typename BaseClass>
+void
+MGDoFHexAccessor<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_fe() != 0, ExcInvalidObject());
+ Assert (dof_indices.size() == (8*dof_handler->get_fe().dofs_per_vertex +
+ 12*dof_handler->get_fe().dofs_per_line +
+ 6*dof_handler->get_fe().dofs_per_quad +
+ dof_handler->get_fe().dofs_per_hex),
+ ExcVectorDoesNotMatch());
+
+ const unsigned int dofs_per_vertex = dof_handler->get_fe().dofs_per_vertex,
+ dofs_per_line = dof_handler->get_fe().dofs_per_line,
+ dofs_per_quad = dof_handler->get_fe().dofs_per_quad,
+ dofs_per_hex = dof_handler->get_fe().dofs_per_hex;
+ vector<int>::iterator next = dof_indices.begin();
+ for (unsigned int vertex=0; vertex<8; ++vertex)
+ for (unsigned int d=0; d<dofs_per_vertex; ++d)
+ *next++ = mg_vertex_dof_index(vertex,d);
+ for (unsigned int line=0; line<12; ++line)
+ for (unsigned int d=0; d<dofs_per_line; ++d)
+ *next++ = this->line(line)->mg_dof_index(d);
+ for (unsigned int quad=0; quad<12; ++quad)
+ for (unsigned int d=0; d<dofs_per_quad; ++d)
+ *next++ = this->quad(quad)->mg_dof_index(d);
+ for (unsigned int d=0; d<dofs_per_hex; ++d)
+ *next++ = mg_dof_index(d);
+};
+
+
+
+template <int dim, typename BaseClass>
+TriaIterator<dim,MGDoFLineAccessor<dim,LineAccessor<dim> > >
+MGDoFHexAccessor<dim,BaseClass>::line (const unsigned int i) const {
+ Assert (i<12, ExcInvalidIndex (i, 0, 12));
+
+ 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,QuadAccessor<dim> > >
+MGDoFHexAccessor<dim,BaseClass>::quad (const unsigned int i) const {
+ Assert (i<12, ExcInvalidIndex (i, 0, 6));
+
+ return TriaIterator<dim,MGDoFQuadAccessor<dim,QuadAccessor<dim> > >
+ (
+ tria,
+ present_level,
+ quad_index (i),
+ mg_dof_handler
+ );
+};
+
+
+
+template <int dim, typename BaseClass>
+TriaIterator<dim,MGDoFHexAccessor<dim,BaseClass> >
+MGDoFHexAccessor<dim,BaseClass>::child (const unsigned int i) const {
+ TriaIterator<dim,MGDoFHexAccessor<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;
+};
+
+
+
+template <int dim, typename BaseClass>
+void
+MGDoFHexAccessor<dim,BaseClass>::copy_from (const MGDoFHexAccessor<dim,BaseClass> &a) {
+ DoFHexAccessor::copy_from (a);
+ set_mg_dof_handler (a.mg_dof_handler);
+};
+
+
+
+
/*------------------------- Functions: MGDoFCellAccessor -----------------------*/
#endif
+#if deal_II_dimension == 3
+template class MGDoFLineAccessor<3,LineAccessor<3> >;
+template class MGDoFQuadAccessor<3,QuadAccessor<3> >;
+template class MGDoFHexAccessor<3,HexAccessor<3> >;
+template class MGDoFHexAccessor<3,CellAccessor<3> >;
+template class MGDoFCellAccessor<3>;
+
+template class TriaRawIterator<3,MGDoFLineAccessor<3,LineAccessor<3> > >;
+template class TriaRawIterator<3,MGDoFQuadAccessor<3,QuadAccessor<3> > >;
+template class TriaRawIterator<3,MGDoFHexAccessor<3,HexAccessor<3> > >;
+template class TriaRawIterator<3,MGDoFCellAccessor<3> >;
+template class TriaIterator<3,MGDoFLineAccessor<3,LineAccessor<3> > >;
+template class TriaIterator<3,MGDoFQuadAccessor<3,QuadAccessor<3> > >;
+template class TriaIterator<3,MGDoFCellAccessor<3> >;
+template class TriaActiveIterator<3,MGDoFLineAccessor<3,LineAccessor<3> > >;
+template class TriaActiveIterator<3,MGDoFQuadAccessor<3,QuadAccessor<3> > >;
+template class TriaActiveIterator<3,MGDoFCellAccessor<3> >;
+#endif
+
+
return 0;
};
+
+template <>
+MGDoFHandler<1>::raw_hex_iterator
+MGDoFHandler<1>::begin_raw_hex (const unsigned int) const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::hex_iterator
+MGDoFHandler<1>::begin_hex (const unsigned int) const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::active_hex_iterator
+MGDoFHandler<1>::begin_active_hex (const unsigned int) const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::raw_hex_iterator
+MGDoFHandler<1>::end_hex () const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::raw_hex_iterator
+MGDoFHandler<1>::last_raw_hex (const unsigned int) const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+template <>
+MGDoFHandler<1>::hex_iterator
+MGDoFHandler<1>::last_hex (const unsigned int) const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::active_hex_iterator
+MGDoFHandler<1>::last_active_hex (const unsigned int) const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::raw_hex_iterator
+MGDoFHandler<1>::last_raw_hex () const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::hex_iterator
+MGDoFHandler<1>::last_hex () const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::active_hex_iterator
+MGDoFHandler<1>::last_active_hex () const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
#endif
return last_active_line (level);
};
+
+
+template <>
+MGDoFHandler<2>::raw_hex_iterator
+MGDoFHandler<2>::begin_raw_hex (const unsigned int) const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::hex_iterator
+MGDoFHandler<2>::begin_hex (const unsigned int) const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::active_hex_iterator
+MGDoFHandler<2>::begin_active_hex (const unsigned int) const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::raw_hex_iterator
+MGDoFHandler<2>::end_hex () const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::raw_hex_iterator
+MGDoFHandler<2>::last_raw_hex (const unsigned int) const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+template <>
+MGDoFHandler<2>::hex_iterator
+MGDoFHandler<2>::last_hex (const unsigned int) const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::active_hex_iterator
+MGDoFHandler<2>::last_active_hex (const unsigned int) const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::raw_hex_iterator
+MGDoFHandler<2>::last_raw_hex () const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::hex_iterator
+MGDoFHandler<2>::last_hex () const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::active_hex_iterator
+MGDoFHandler<2>::last_active_hex () const {
+ Assert (false, ExcNotImplemented());
+ return 0;
+};
+
+
+
+#endif
+
+
+
+
+#if deal_II_dimension == 3
+
+template <>
+MGDoFHandler<3>::raw_cell_iterator
+MGDoFHandler<3>::begin_raw (const unsigned int level) const {
+ return begin_raw_hex (level);
+};
+
+
+
+template <>
+MGDoFHandler<3>::cell_iterator
+MGDoFHandler<3>::begin (const unsigned int level) const {
+ return begin_hex (level);
+};
+
+
+
+template <>
+MGDoFHandler<3>::active_cell_iterator
+MGDoFHandler<3>::begin_active (const unsigned int level) const {
+ return begin_active_hex (level);
+};
+
+
+
+template <>
+MGDoFHandler<3>::raw_cell_iterator
+MGDoFHandler<3>::end () const {
+ return end_hex ();
+};
+
+
+
+template <>
+MGDoFHandler<3>::raw_cell_iterator
+MGDoFHandler<3>::last_raw () const {
+ return last_raw_hex ();
+};
+
+
+
+template <>
+MGDoFHandler<3>::raw_cell_iterator
+MGDoFHandler<3>::last_raw (const unsigned int level) const {
+ return last_raw_hex (level);
+};
+
+
+
+template <>
+MGDoFHandler<3>::cell_iterator
+MGDoFHandler<3>::last () const {
+ return last_hex ();
+};
+
+
+
+template <>
+MGDoFHandler<3>::cell_iterator
+MGDoFHandler<3>::last (const unsigned int level) const {
+ return last_hex (level);
+};
+
+
+
+template <>
+MGDoFHandler<3>::active_cell_iterator
+MGDoFHandler<3>::last_active () const {
+ return last_active_hex ();
+};
+
+
+
+template <>
+MGDoFHandler<3>::active_cell_iterator
+MGDoFHandler<3>::last_active (const unsigned int level) const {
+ return last_active_hex (level);
+};
+
+
+template <>
+MGDoFDimensionInfo<3>::raw_face_iterator
+MGDoFHandler<3>::begin_raw_face (const unsigned int level) const {
+ return begin_raw_quad (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::face_iterator
+MGDoFHandler<3>::begin_face (const unsigned int level) const {
+ return begin_quad (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::active_face_iterator
+MGDoFHandler<3>::begin_active_face (const unsigned int level) const {
+ return begin_active_quad (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::raw_face_iterator
+MGDoFHandler<3>::end_face () const {
+ return end_quad ();
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::raw_face_iterator
+MGDoFHandler<3>::last_raw_face () const {
+ return last_raw_quad ();
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::raw_face_iterator
+MGDoFHandler<3>::last_raw_face (const unsigned int level) const {
+ return last_raw_quad (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::face_iterator
+MGDoFHandler<3>::last_face () const {
+ return last_quad ();
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::face_iterator
+MGDoFHandler<3>::last_face (const unsigned int level) const {
+ return last_quad (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::active_face_iterator
+MGDoFHandler<3>::last_active_face () const {
+ return last_active_quad ();
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::active_face_iterator
+MGDoFHandler<3>::last_active_face (const unsigned int level) const {
+ return last_active_quad (level);
+};
+
#endif
+template <int dim>
+typename MGDoFHandler<dim>::raw_hex_iterator
+MGDoFHandler<dim>::begin_raw_hex (const unsigned int level) const {
+ return raw_hex_iterator (tria,
+ tria->begin_raw_hex(level)->level(),
+ tria->begin_raw_hex(level)->index(),
+ this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::hex_iterator
+MGDoFHandler<dim>::begin_hex (const unsigned int level) const {
+ return hex_iterator (tria,
+ tria->begin_hex(level)->level(),
+ tria->begin_hex(level)->index(),
+ this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::active_hex_iterator
+MGDoFHandler<dim>::begin_active_hex (const unsigned int level) const {
+ return active_hex_iterator (tria,
+ tria->begin_active_hex(level)->level(),
+ tria->begin_active_hex(level)->index(),
+ this);
+};
+
+
+
template <int dim>
typename MGDoFHandler<dim>::raw_line_iterator
MGDoFHandler<dim>::end_line () const {
+template <int dim>
+typename MGDoFHandler<dim>::raw_hex_iterator
+MGDoFHandler<dim>::end_hex () const {
+ return raw_hex_iterator (tria, -1, -1, this);
+};
+
+
+
template <int dim>
typename MGDoFDimensionInfo<dim>::raw_cell_iterator
MGDoFHandler<dim>::end_raw (const unsigned int level) const {
-
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::raw_quad_iterator
+MGDoFHandler<dim>::end_raw_quad (const unsigned int level) const {
+ return (level == mg_levels.size()-1 ?
+ end_quad() :
+ begin_raw_quad (level+1));
+};
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::quad_iterator
+MGDoFHandler<dim>::end_quad (const unsigned int level) const {
+ return (level == mg_levels.size()-1 ?
+ quad_iterator(end_quad()) :
+ begin_quad (level+1));
+};
+
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::active_quad_iterator
+MGDoFHandler<dim>::end_active_quad (const unsigned int level) const {
+ return (level == mg_levels.size()-1 ?
+ active_quad_iterator(end_quad()) :
+ begin_active_quad (level+1));
+};
+
+
+
template <int dim>
-typename MGDoFDimensionInfo<dim>::raw_quad_iterator
-MGDoFHandler<dim>::end_raw_quad (const unsigned int level) const {
+typename MGDoFDimensionInfo<dim>::raw_hex_iterator
+MGDoFHandler<dim>::end_raw_hex (const unsigned int level) const {
return (level == mg_levels.size()-1 ?
- end_quad() :
- begin_raw_quad (level+1));
+ end_hex() :
+ begin_raw_hex (level+1));
};
template <int dim>
-typename MGDoFDimensionInfo<dim>::quad_iterator
-MGDoFHandler<dim>::end_quad (const unsigned int level) const {
+typename MGDoFDimensionInfo<dim>::hex_iterator
+MGDoFHandler<dim>::end_hex (const unsigned int level) const {
return (level == mg_levels.size()-1 ?
- quad_iterator(end_quad()) :
- begin_quad (level+1));
+ hex_iterator(end_hex()) :
+ begin_hex (level+1));
};
template <int dim>
-typename MGDoFDimensionInfo<dim>::active_quad_iterator
-MGDoFHandler<dim>::end_active_quad (const unsigned int level) const {
+typename MGDoFDimensionInfo<dim>::active_hex_iterator
+MGDoFHandler<dim>::end_active_hex (const unsigned int level) const {
return (level == mg_levels.size()-1 ?
- active_quad_iterator(end_quad()) :
- begin_active_quad (level+1));
+ active_hex_iterator(end_hex()) :
+ begin_active_hex (level+1));
};
+template <int dim>
+typename MGDoFHandler<dim>::raw_hex_iterator
+MGDoFHandler<dim>::last_raw_hex (const unsigned int level) const {
+ return raw_hex_iterator (tria,
+ tria->last_raw_hex(level)->level(),
+ tria->last_raw_hex(level)->index(),
+ this);
+};
+
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::hex_iterator
+MGDoFHandler<dim>::last_hex (const unsigned int level) const {
+ return hex_iterator (tria,
+ tria->last_hex(level)->level(),
+ tria->last_hex(level)->index(),
+ this);
+};
+
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::active_hex_iterator
+MGDoFHandler<dim>::last_active_hex (const unsigned int level) const {
+ return active_hex_iterator (tria,
+ tria->last_active_hex(level)->level(),
+ tria->last_active_hex(level)->index(),
+ this);
+};
+
+
+
+
template <int dim>
typename MGDoFHandler<dim>::raw_line_iterator
MGDoFHandler<dim>::last_raw_line () const {
+template <int dim>
+typename MGDoFHandler<dim>::raw_hex_iterator
+MGDoFHandler<dim>::last_raw_hex () const {
+ return last_raw_hex (mg_levels.size()-1);
+};
+
+
+
template <int dim>
typename MGDoFHandler<dim>::line_iterator
MGDoFHandler<dim>::last_line () const {
+template <int dim>
+typename MGDoFHandler<dim>::hex_iterator
+MGDoFHandler<dim>::last_hex () const {
+ return last_hex (mg_levels.size()-1);
+};
+
+
+
template <int dim>
typename MGDoFHandler<dim>::active_line_iterator
MGDoFHandler<dim>::last_active_line () const {
+template <int dim>
+typename MGDoFHandler<dim>::active_hex_iterator
+MGDoFHandler<dim>::last_active_hex () const {
+ return last_active_hex (mg_levels.size()-1);
+};
+
+
+
+#if deal_II_dimension == 3
+
+template <>
+unsigned int
+MGDoFHandler<3>::distribute_dofs_on_cell (cell_iterator &cell,
+ unsigned int next_free_dof) {
+ if (selected_fe->dofs_per_vertex > 0)
+ // number dofs on vertices
+ for (unsigned int vertex=0; vertex<GeometryInfo<3>::vertices_per_cell; ++vertex)
+ // check whether dofs for this
+ // vertex have been distributed
+ // (only check the first dof)
+ if (cell->mg_vertex_dof_index(vertex, 0) == -1)
+ for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
+ cell->set_mg_vertex_dof_index (vertex, d, next_free_dof++);
+
+ // for the lines
+ if (selected_fe->dofs_per_line > 0)
+ for (unsigned int l=0; l<GeometryInfo<3>::lines_per_cell; ++l)
+ {
+ line_iterator line = cell->line(l);
+
+ // distribute dofs if necessary:
+ // check whether line dof is already
+ // numbered (check only first dof)
+ if (line->mg_dof_index(0) == -1)
+ // if not: distribute dofs
+ for (unsigned int d=0; d<selected_fe->dofs_per_line; ++d)
+ line->set_mg_dof_index (d, next_free_dof++);
+ };
+
+ // for the quads
+ if (selected_fe->dofs_per_quad > 0)
+ for (unsigned int q=0; q<GeometryInfo<3>::lines_per_cell; ++q)
+ {
+ quad_iterator quad = cell->quad(q);
+
+ // distribute dofs if necessary:
+ // check whether line dof is already
+ // numbered (check only first dof)
+ if (quad->mg_dof_index(0) == -1)
+ // if not: distribute dofs
+ for (unsigned int d=0; d<selected_fe->dofs_per_line; ++d)
+ quad->set_mg_dof_index (d, next_free_dof++);
+ };
+
+
+ // dofs of cell
+ if (selected_fe->dofs_per_hex > 0)
+ for (unsigned int d=0; d<selected_fe->dofs_per_hex; ++d)
+ cell->set_mg_dof_index (d, next_free_dof++);
+
+
+ // note that this cell has been processed
+ cell->set_user_flag ();
+
+ return next_free_dof;
+};
+
+#endif
+
+
+
template <int dim>
unsigned int MGDoFHandler<dim>::n_dofs (const unsigned int level) const {
-
-template <int dim>
-void MGDoFHandler<dim>::renumber_dofs (const unsigned int level,
- const RenumberingMethod method,
- const vector<int> &starting_points) {
- // make the connection graph
- SparseMatrixStruct sparsity (n_dofs(level), max_couplings_between_dofs());
- make_sparsity_pattern (level, sparsity);
-
- int n_dofs = sparsity.n_rows();
- // store the new dof numbers; -1 means
- // that no new number was chosen yet
- vector<int> new_number(n_dofs, -1);
-
- // store the indices of the dofs renumbered
- // in the last round. Default to starting
- // points
- vector<int> last_round_dofs (starting_points);
-
- // delete disallowed elements
- for (unsigned int i=0; i<last_round_dofs.size(); ++i)
- if ((last_round_dofs[i]<0) || (last_round_dofs[i]>=n_dofs))
- last_round_dofs[i] = -1;
-
- remove_if (last_round_dofs.begin(), last_round_dofs.end(),
- bind2nd(equal_to<int>(), -1));
-
- // now if no valid points remain:
- // find dof with lowest coordination
- // number
-
- if (last_round_dofs.size() == 0)
- {
- int starting_point = -1;
- unsigned int min_coordination = n_dofs;
- for (int row=0; row<n_dofs; ++row)
- {
- unsigned int j;
- for (j=sparsity.get_rowstart_indices()[row];
- j<sparsity.get_rowstart_indices()[row+1]; ++j)
- if (sparsity.get_column_numbers()[j] == -1)
- break;
- // post-condition after loop:
- // coordination is now
- // j-rowstart[row]
- if (j-sparsity.get_rowstart_indices()[row] < min_coordination)
- {
- min_coordination = j-sparsity.get_rowstart_indices()[row];
- starting_point = row;
- };
- };
- // initialize the first dof
- last_round_dofs.push_back (starting_point);
- };
-
-
- // store next free dof index
- int next_free_number = 0;
-
- // enumerate the first round dofs
- for (unsigned int i=0; i!=last_round_dofs.size(); ++i)
- new_number[last_round_dofs[i]] = next_free_number++;
-
- bool all_dofs_renumbered = false;
-
- // now do as many steps as needed to
- // renumber all dofs
- while (!all_dofs_renumbered)
- {
- // store the indices of the dofs to be
- // renumbered in the next round
- vector<int> next_round_dofs;
-
- // find all neighbors of the
- // dofs numbered in the last
- // round
- for (unsigned int i=0; i<last_round_dofs.size(); ++i)
- for (unsigned int j=sparsity.get_rowstart_indices()[last_round_dofs[i]];
- j<sparsity.get_rowstart_indices()[last_round_dofs[i]+1]; ++j)
- if (sparsity.get_column_numbers()[j] == -1)
- break;
- else
- next_round_dofs.push_back (sparsity.get_column_numbers()[j]);
-
- // sort dof numbers
- sort (next_round_dofs.begin(), next_round_dofs.end());
-
- // delete multiple entries
- vector<int>::iterator end_sorted;
- end_sorted = unique (next_round_dofs.begin(), next_round_dofs.end());
- next_round_dofs.erase (end_sorted, next_round_dofs.end());
-
- // eliminate dofs which are
- // already numbered
- for (int s=next_round_dofs.size()-1; s>=0; --s)
- if (new_number[next_round_dofs[s]] != -1)
- next_round_dofs.erase (&next_round_dofs[s]);
-
- // check whether there are any new
- // dofs
- all_dofs_renumbered = (next_round_dofs.size() == 0);
- if (all_dofs_renumbered)
- // end loop if possible
- continue;
-
-
- // store for each coordination
- // number the dofs with these
- // coordination number
- multimap<unsigned int, int> dofs_by_coordination;
-
- // find coordination number for
- // each of these dofs
- for (vector<int>::iterator s=next_round_dofs.begin();
- s!=next_round_dofs.end(); ++s)
- {
- unsigned int coordination = 0;
- for (unsigned int j=sparsity.get_rowstart_indices()[*s];
- j<sparsity.get_rowstart_indices()[*s+1]; ++j)
- if (sparsity.get_column_numbers()[j] == -1)
- break;
- else
- ++coordination;
-
- // insert this dof at its
- // coordination number
- const pair<const unsigned int, int> new_entry (coordination, *s);
- dofs_by_coordination.insert (new_entry);
- };
-
- ////
- multimap<unsigned int, int>::iterator i;
- for (i = dofs_by_coordination.begin(); i!=dofs_by_coordination.end(); ++i)
- new_number[i->second] = next_free_number++;
-
- // after that: copy this round's
- // dofs for the next round
- last_round_dofs = next_round_dofs;
- };
-
-#ifdef DEBUG
- // test for all indices numbered
- if (find (new_number.begin(), new_number.end(), -1) != new_number.end())
- Assert (false, ExcRenumberingIncomplete());
- Assert (next_free_number == n_dofs,
- ExcRenumberingIncomplete());
-#endif
-
- switch (method)
- {
- case Cuthill_McKee:
- break;
- case reverse_Cuthill_McKee:
- {
- for (vector<int>::iterator i=new_number.begin(); i!=new_number.end(); ++i)
- *i = n_dofs-*i;
- break;
- };
- default:
- Assert (false, ExcNotImplemented());
- };
-
- // actually perform renumbering;
- // this is dimension specific and
- // thus needs an own function
- do_renumbering (level, new_number);
-};
-
-
-
-
#if deal_II_dimension == 1
template <>
-void MGDoFHandler<1>::do_renumbering (const unsigned int level,
- const vector<int> &new_numbers) {
+void MGDoFHandler<1>::renumber_dofs (const unsigned int level,
+ const vector<int> &new_numbers) {
Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete());
// note that we can not use cell iterators
#if deal_II_dimension == 2
template <>
-void MGDoFHandler<2>::do_renumbering (const unsigned int level,
- const vector<int> &new_numbers) {
+void MGDoFHandler<2>::renumber_dofs (const unsigned int level,
+ const vector<int> &new_numbers) {
+ Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete());
+
+ for (vector<MGVertexDoFs>::iterator i=mg_vertex_dofs.begin();
+ i!=mg_vertex_dofs.end(); ++i)
+ // if the present vertex lives on
+ // the present level
+ if ((i->get_coarsest_level() <= level) &&
+ (i->get_finest_level() >= level))
+ for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
+ i->set_index (level, d, selected_fe->dofs_per_vertex,
+ new_numbers[i->get_index (level, d,
+ selected_fe->dofs_per_vertex)]);
+
+ for (vector<int>::iterator i=mg_levels[level]->line_dofs.begin();
+ i!=mg_levels[level]->line_dofs.end(); ++i)
+ {
+ Assert (*i != -1, ExcInternalError());
+ *i = new_numbers[*i];
+ };
+
+ for (vector<int>::iterator i=mg_levels[level]->quad_dofs.begin();
+ i!=mg_levels[level]->quad_dofs.end(); ++i)
+ {
+ Assert (*i != -1, ExcInternalError());
+ *i = new_numbers[*i];
+ };
+};
+
+#endif
+
+
+
+#if deal_II_dimension == 3
+
+template <>
+void MGDoFHandler<3>::renumber_dofs (const unsigned int level,
+ const vector<int> &new_numbers) {
Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete());
for (vector<MGVertexDoFs>::iterator i=mg_vertex_dofs.begin();
Assert (*i != -1, ExcInternalError());
*i = new_numbers[*i];
};
+
for (vector<int>::iterator i=mg_levels[level]->quad_dofs.begin();
i!=mg_levels[level]->quad_dofs.end(); ++i)
{
Assert (*i != -1, ExcInternalError());
*i = new_numbers[*i];
};
+
+ for (vector<int>::iterator i=mg_levels[level]->hex_dofs.begin();
+ i!=mg_levels[level]->hex_dofs.end(); ++i)
+ {
+ Assert (*i != -1, ExcInternalError());
+ *i = new_numbers[*i];
+ };
};
#endif
#endif
+
+
+#if deal_II_dimension == 3
+
+template <>
+void MGDoFHandler<3>::reserve_space () {
+ const unsigned int dim = 3;
+
+ Assert (selected_fe != 0, ExcNoFESelected());
+ Assert (tria->n_levels() > 0, ExcInvalidTriangulation());
+
+ ////////////////////////////
+ // DESTRUCTION
+
+ // delete all levels and set them up
+ // newly, since vectors are
+ // troublesome if you want to change
+ // their size
+ for (unsigned int i=0; i<mg_levels.size(); ++i)
+ delete mg_levels[i];
+ mg_levels.clear ();
+
+ // also delete vector of vertex indices
+ // this calls the destructor which
+ // must free the space
+ mg_vertex_dofs.clear ();
+
+
+ ////////////////////////////
+ // CONSTRUCTION
+
+ // first allocate space for the
+ // lines and quads on each level
+ for (unsigned int i=0; i<tria->n_levels(); ++i)
+ {
+ mg_levels.push_back (new DoFLevel<3>);
+
+ mg_levels.back()->line_dofs = vector<int> (tria->levels[i]->lines.lines.size() *
+ selected_fe->dofs_per_line,
+ -1);
+ mg_levels.back()->quad_dofs = vector<int> (tria->levels[i]->quads.quads.size() *
+ selected_fe->dofs_per_quad,
+ -1);
+ mg_levels.back()->hex_dofs = vector<int> (tria->levels[i]->hexes.hexes.size() *
+ selected_fe->dofs_per_hex,
+ -1);
+ };
+
+
+ // now allocate space for the
+ // vertices. To this end, we need
+ // to construct as many objects as
+ // there are vertices and let them
+ // allocate enough space for their
+ // vertex indices on the levels they
+ // live on. We need therefore to
+ // count to how many levels a cell
+ // belongs to, which we do by looping
+ // over all cells and storing the
+ // maximum and minimum level each
+ // vertex we pass by belongs to
+ mg_vertex_dofs.resize (tria->vertices.size());
+
+ vector<unsigned int> min_level (tria->vertices.size(), tria->n_levels());
+ vector<unsigned int> max_level (tria->vertices.size(), 0);
+
+ Triangulation<dim>::cell_iterator cell = tria->begin(),
+ endc = tria->end();
+ for (; cell!=endc; ++cell)
+ for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
+ ++vertex)
+ {
+ const unsigned int vertex_index = cell->vertex_index(vertex);
+ if (min_level[vertex_index] > static_cast<unsigned int>(cell->level()))
+ min_level[vertex_index] = cell->level();
+ if (max_level[vertex_index] < static_cast<unsigned int>(cell->level()))
+ max_level[vertex_index] = cell->level();
+ };
+
+
+ // now allocate the needed space
+ for (unsigned int vertex=0; vertex<tria->vertices.size(); ++vertex)
+ {
+ Assert (min_level[vertex] < tria->n_levels(), ExcInternalError());
+ Assert (max_level[vertex] >= min_level[vertex], ExcInternalError());
+
+ mg_vertex_dofs[vertex].init (min_level[vertex],
+ max_level[vertex],
+ selected_fe->dofs_per_vertex);
+ };
+};
+
+#endif
+
+
// explicite instantiations
template class MGDoFHandler<deal_II_dimension>;
#include <lac/sparsematrix.h>
#include <base/parameter_handler.h>
#include <grid/dof_constraints.h>
-
+#include <numerics/dof_renumbering.h>
#include <fstream>
#include <cmath>
dof->distribute_dofs (fe);
cout << " Renumbering degrees of freedom..." << endl;
- dof->renumber_dofs (Cuthill_McKee, false);
+ DoFRenumbering::renumber_Cuthill_McKee (*dof);
SparseMatrixStruct sparsity (dof->n_dofs(),
dof->max_couplings_between_dofs());
#include <grid/tria_accessor.h>
#include <lac/sparsematrix.h>
#include <grid/dof_constraints.h>
-
+#include <numerics/dof_renumbering.h>
#include <cmath>
#include <cstdlib>
dof->distribute_dofs (fe);
cout << " Renumbering degrees of freedom..." << endl;
- dof->renumber_dofs (Cuthill_McKee, false);
+ DoFRenumbering::renumber_Cuthill_McKee (*dof);
SparseMatrixStruct sparsity (dof->n_dofs(),
dof->max_couplings_between_dofs());