From 5341745b377cf5447aaf7c72f1c7ef106e1e5487 Mon Sep 17 00:00:00 2001 From: wolf Date: Tue, 20 Apr 1999 11:10:33 +0000 Subject: [PATCH] Separate the dof renumbering code from the dof handlers into a new class. Add several missing functions for 3d to the MGDoFHandler class to allow for testing the first item. git-svn-id: https://svn.dealii.org/trunk@1190 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/Attic/examples/dof/dof_test.cc | 4 +- deal.II/deal.II/include/dofs/dof_handler.h | 152 +-- .../deal.II/include/dofs/mg_dof_accessor.h | 4 +- deal.II/deal.II/include/dofs/mg_dof_handler.h | 141 ++- .../include/grid/tria_accessor.templates.h | 29 + .../include/multigrid/mg_dof_accessor.h | 4 +- .../include/multigrid/mg_dof_handler.h | 141 ++- .../include/numerics/dof_renumbering.h | 205 ++++ deal.II/deal.II/source/dofs/dof_handler.cc | 196 ---- .../deal.II/source/dofs/dof_renumbering.cc | 379 ++++++++ .../deal.II/source/dofs/mg_dof_accessor.cc | 191 ++++ deal.II/deal.II/source/dofs/mg_dof_handler.cc | 883 ++++++++++++++---- .../source/multigrid/mg_dof_accessor.cc | 191 ++++ .../source/multigrid/mg_dof_handler.cc | 883 ++++++++++++++---- tests/big-tests/dof/dof_test.cc | 4 +- tests/deal.II/dof_test.cc | 4 +- 16 files changed, 2628 insertions(+), 783 deletions(-) create mode 100644 deal.II/deal.II/include/numerics/dof_renumbering.h create mode 100644 deal.II/deal.II/source/dofs/dof_renumbering.cc diff --git a/deal.II/deal.II/Attic/examples/dof/dof_test.cc b/deal.II/deal.II/Attic/examples/dof/dof_test.cc index cc8a9f83f5..8ca5178abf 100644 --- a/deal.II/deal.II/Attic/examples/dof/dof_test.cc +++ b/deal.II/deal.II/Attic/examples/dof/dof_test.cc @@ -12,7 +12,7 @@ #include #include #include - +#include #include #include @@ -354,7 +354,7 @@ void TestCases::run (ParameterHandler &prm) { 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()); diff --git a/deal.II/deal.II/include/dofs/dof_handler.h b/deal.II/deal.II/include/dofs/dof_handler.h index b32a35c04a..6252b03017 100644 --- a/deal.II/deal.II/include/dofs/dof_handler.h +++ b/deal.II/deal.II/include/dofs/dof_handler.h @@ -129,29 +129,6 @@ class DoFDimensionInfo<3> { -/** - * 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 @@ -203,117 +180,14 @@ enum RenumberingMethod { * * 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 @@ -550,22 +424,6 @@ class DoFHandler : public DoFDimensionInfo { */ 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 &starting_points = vector()); - /** * Actually do the renumbering based on * a list of new dof numbers for all the @@ -1513,7 +1371,7 @@ class DoFHandler : public DoFDimensionInfo { * Exception */ DeclException0 (ExcNoFESelected); - /** + /** * Exception */ DeclException0 (ExcRenumberingIncomplete); diff --git a/deal.II/deal.II/include/dofs/mg_dof_accessor.h b/deal.II/deal.II/include/dofs/mg_dof_accessor.h index 31e9920f5e..36109d757c 100644 --- a/deal.II/deal.II/include/dofs/mg_dof_accessor.h +++ b/deal.II/deal.II/include/dofs/mg_dof_accessor.h @@ -416,7 +416,7 @@ class MGDoFHexAccessor : public MGDoFAccessor, public DoFHexAccessor > > + TriaIterator > > quad (const unsigned int i) const; /** @@ -678,6 +678,8 @@ class MGDoFCellAccessor : public MGDoFSubstructAccessor { + + /*---------------------------- mg_dof_accessor.h ---------------------------*/ /* end of #ifndef __mg_dof_accessor_H */ #endif diff --git a/deal.II/deal.II/include/dofs/mg_dof_handler.h b/deal.II/deal.II/include/dofs/mg_dof_handler.h index edf50654ff..1d6a7d7fb1 100644 --- a/deal.II/deal.II/include/dofs/mg_dof_handler.h +++ b/deal.II/deal.II/include/dofs/mg_dof_handler.h @@ -131,6 +131,10 @@ class MGDoFHandler : public DoFHandler { typedef typename MGDoFDimensionInfo::quad_iterator quad_iterator; typedef typename MGDoFDimensionInfo::active_quad_iterator active_quad_iterator; + typedef typename MGDoFDimensionInfo::raw_hex_iterator raw_hex_iterator; + typedef typename MGDoFDimensionInfo::hex_iterator hex_iterator; + typedef typename MGDoFDimensionInfo::active_hex_iterator active_hex_iterator; + typedef typename MGDoFDimensionInfo::raw_cell_iterator raw_cell_iterator; typedef typename MGDoFDimensionInfo::cell_iterator cell_iterator; typedef typename MGDoFDimensionInfo::active_cell_iterator active_cell_iterator; @@ -170,20 +174,18 @@ class MGDoFHandler : public DoFHandler { virtual void distribute_dofs (const FiniteElement &); /** - * 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 &starting_points = vector()); + void renumber_dofs (const unsigned int level, + const vector &new_numbers); /** * Write the sparsity structure of the @@ -653,6 +655,103 @@ class MGDoFHandler : public DoFHandler { /*---------------------------------------*/ + /** + * @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. @@ -831,22 +930,6 @@ class MGDoFHandler : public DoFHandler { 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 &new_numbers); - /** * Reserve enough space for the MG dof * indices for a given triangulation. @@ -881,6 +964,8 @@ class MGDoFHandler : public DoFHandler { friend class MGDoFLineAccessor >; friend class MGDoFQuadAccessor >; friend class MGDoFQuadAccessor >; + friend class MGDoFHexAccessor >; + friend class MGDoFHexAccessor >; }; diff --git a/deal.II/deal.II/include/grid/tria_accessor.templates.h b/deal.II/deal.II/include/grid/tria_accessor.templates.h index 6035ebf6bc..c1addddbd2 100644 --- a/deal.II/deal.II/include/grid/tria_accessor.templates.h +++ b/deal.II/deal.II/include/grid/tria_accessor.templates.h @@ -473,6 +473,35 @@ HexAccessor::quad (const unsigned int i) const { +template +inline +unsigned int +HexAccessor::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 inline unsigned int diff --git a/deal.II/deal.II/include/multigrid/mg_dof_accessor.h b/deal.II/deal.II/include/multigrid/mg_dof_accessor.h index 31e9920f5e..36109d757c 100644 --- a/deal.II/deal.II/include/multigrid/mg_dof_accessor.h +++ b/deal.II/deal.II/include/multigrid/mg_dof_accessor.h @@ -416,7 +416,7 @@ class MGDoFHexAccessor : public MGDoFAccessor, public DoFHexAccessor > > + TriaIterator > > quad (const unsigned int i) const; /** @@ -678,6 +678,8 @@ class MGDoFCellAccessor : public MGDoFSubstructAccessor { + + /*---------------------------- mg_dof_accessor.h ---------------------------*/ /* end of #ifndef __mg_dof_accessor_H */ #endif diff --git a/deal.II/deal.II/include/multigrid/mg_dof_handler.h b/deal.II/deal.II/include/multigrid/mg_dof_handler.h index edf50654ff..1d6a7d7fb1 100644 --- a/deal.II/deal.II/include/multigrid/mg_dof_handler.h +++ b/deal.II/deal.II/include/multigrid/mg_dof_handler.h @@ -131,6 +131,10 @@ class MGDoFHandler : public DoFHandler { typedef typename MGDoFDimensionInfo::quad_iterator quad_iterator; typedef typename MGDoFDimensionInfo::active_quad_iterator active_quad_iterator; + typedef typename MGDoFDimensionInfo::raw_hex_iterator raw_hex_iterator; + typedef typename MGDoFDimensionInfo::hex_iterator hex_iterator; + typedef typename MGDoFDimensionInfo::active_hex_iterator active_hex_iterator; + typedef typename MGDoFDimensionInfo::raw_cell_iterator raw_cell_iterator; typedef typename MGDoFDimensionInfo::cell_iterator cell_iterator; typedef typename MGDoFDimensionInfo::active_cell_iterator active_cell_iterator; @@ -170,20 +174,18 @@ class MGDoFHandler : public DoFHandler { virtual void distribute_dofs (const FiniteElement &); /** - * 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 &starting_points = vector()); + void renumber_dofs (const unsigned int level, + const vector &new_numbers); /** * Write the sparsity structure of the @@ -653,6 +655,103 @@ class MGDoFHandler : public DoFHandler { /*---------------------------------------*/ + /** + * @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. @@ -831,22 +930,6 @@ class MGDoFHandler : public DoFHandler { 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 &new_numbers); - /** * Reserve enough space for the MG dof * indices for a given triangulation. @@ -881,6 +964,8 @@ class MGDoFHandler : public DoFHandler { friend class MGDoFLineAccessor >; friend class MGDoFQuadAccessor >; friend class MGDoFQuadAccessor >; + friend class MGDoFHexAccessor >; + friend class MGDoFHexAccessor >; }; diff --git a/deal.II/deal.II/include/numerics/dof_renumbering.h b/deal.II/deal.II/include/numerics/dof_renumbering.h new file mode 100644 index 0000000000..62c748c761 --- /dev/null +++ b/deal.II/deal.II/include/numerics/dof_renumbering.h @@ -0,0 +1,205 @@ +/*---------------------------- dof_renumbering.h ---------------------------*/ +/* $Id$ */ +#ifndef __dof_renumbering_H +#define __dof_renumbering_H +/*---------------------------- dof_renumbering.h ---------------------------*/ + + +#include +#include + + + + +/** + * 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 + static void renumber_Cuthill_McKee (DoFHandler &dof_handler, + const bool reversed_numbering = false, + const bool use_constraints = false, + const vector &starting_indices = vector()); + + /** + * 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 + static void renumber_Cuthill_McKee (MGDoFHandler &dof_handler, + const unsigned int level, + const bool reversed_numbering = false, + const vector &starting_indices = vector ()); + + + /** + * Exception + */ + DeclException0 (ExcRenumberingIncomplete); +}; + + + + + + + + +/*---------------------------- dof_renumbering.h ---------------------------*/ +/* end of #ifndef __dof_renumbering_H */ +#endif +/*---------------------------- dof_renumbering.h ---------------------------*/ diff --git a/deal.II/deal.II/source/dofs/dof_handler.cc b/deal.II/deal.II/source/dofs/dof_handler.cc index 982dae0a86..fdb33a5dff 100644 --- a/deal.II/deal.II/source/dofs/dof_handler.cc +++ b/deal.II/deal.II/source/dofs/dof_handler.cc @@ -16,7 +16,6 @@ #include #include -#include #include #include @@ -1565,201 +1564,6 @@ unsigned int DoFHandler<3>::distribute_dofs_on_cell (active_cell_iterator &cell, -template -void DoFHandler::renumber_dofs (const RenumberingMethod method, - const bool use_constraints, - const vector &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 new_number(sparsity.n_rows(), -1); - - // store the indices of the dofs renumbered - // in the last round. Default to starting - // points - vector last_round_dofs (starting_points); - - // delete disallowed elements - for (unsigned int i=0; i=n_dofs)) - last_round_dofs[i] = -1; - - remove_if (last_round_dofs.begin(), last_round_dofs.end(), - bind2nd(equal_to(), -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 next_round_dofs; - - // find all neighbors of the - // dofs numbered in the last - // round - for (unsigned int i=0; i::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 dofs_by_coordination; - - // find coordination number for - // each of these dofs - for (vector::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 new_entry (coordination, *s); - dofs_by_coordination.insert (new_entry); - }; - - //// - multimap::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::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 <> diff --git a/deal.II/deal.II/source/dofs/dof_renumbering.cc b/deal.II/deal.II/source/dofs/dof_renumbering.cc new file mode 100644 index 0000000000..ca68898d72 --- /dev/null +++ b/deal.II/deal.II/source/dofs/dof_renumbering.cc @@ -0,0 +1,379 @@ +/* $Id$ */ + +#include +#include +#include +#include +#include + +#include +#include +#include + + + +template +void DoFRenumbering::renumber_Cuthill_McKee (DoFHandler &dof_handler, + const bool reversed_numbering, + const bool use_constraints, + const vector &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 new_number(sparsity.n_rows(), -1); + + // store the indices of the dofs renumbered + // in the last round. Default to starting + // points + vector last_round_dofs (starting_indices); + + // delete disallowed elements + for (unsigned int i=0; i=n_dofs)) + last_round_dofs[i] = -1; + + remove_if (last_round_dofs.begin(), last_round_dofs.end(), + bind2nd(equal_to(), -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 next_round_dofs; + + // find all neighbors of the + // dofs numbered in the last + // round + for (unsigned int i=0; i::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 dofs_by_coordination; + + // find coordination number for + // each of these dofs + for (vector::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 new_entry (coordination, *s); + dofs_by_coordination.insert (new_entry); + }; + + //// + multimap::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::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 +void DoFRenumbering::renumber_Cuthill_McKee (MGDoFHandler &dof_handler, + const unsigned int level, + const bool reversed_numbering, + const vector &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 new_number(n_dofs, -1); + + // store the indices of the dofs renumbered + // in the last round. Default to starting + // points + vector last_round_dofs (starting_indices); + + // delete disallowed elements + for (unsigned int i=0; i=n_dofs)) + last_round_dofs[i] = -1; + + remove_if (last_round_dofs.begin(), last_round_dofs.end(), + bind2nd(equal_to(), -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 next_round_dofs; + + // find all neighbors of the + // dofs numbered in the last + // round + for (unsigned int i=0; i::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 dofs_by_coordination; + + // find coordination number for + // each of these dofs + for (vector::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 new_entry (coordination, *s); + dofs_by_coordination.insert (new_entry); + }; + + //// + multimap::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::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 &dof_handler, + const bool reversed_numbering, + const bool use_constraints, + const vector &starting_indices); + +template +void DoFRenumbering::renumber_Cuthill_McKee (MGDoFHandler &dof_handler, + const unsigned int level, + const bool reversed_numbering, + const vector &starting_indices); + diff --git a/deal.II/deal.II/source/dofs/mg_dof_accessor.cc b/deal.II/deal.II/source/dofs/mg_dof_accessor.cc index 101f9344fe..d6fd3b0b24 100644 --- a/deal.II/deal.II/source/dofs/mg_dof_accessor.cc +++ b/deal.II/deal.II/source/dofs/mg_dof_accessor.cc @@ -299,6 +299,177 @@ MGDoFQuadAccessor::copy_from (const MGDoFQuadAccessor +MGDoFHexAccessor::MGDoFHexAccessor (Triangulation *tria, + const int level, + const int index, + const AccessorData *local_data) : + MGDoFAccessor (local_data), + DoFHexAccessor(tria,level,index,local_data) {}; + + + +template +inline +int MGDoFHexAccessor::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 (iget_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 +void MGDoFHexAccessor::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 (iget_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 +inline +int MGDoFHexAccessor::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 (iget_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 +void MGDoFHexAccessor::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 (iget_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 +void +MGDoFHexAccessor::get_mg_dof_indices (vector &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::iterator next = dof_indices.begin(); + for (unsigned int vertex=0; vertex<8; ++vertex) + for (unsigned int d=0; dline(line)->mg_dof_index(d); + for (unsigned int quad=0; quad<12; ++quad) + for (unsigned int d=0; dquad(quad)->mg_dof_index(d); + for (unsigned int d=0; d +TriaIterator > > +MGDoFHexAccessor::line (const unsigned int i) const { + Assert (i<12, ExcInvalidIndex (i, 0, 12)); + + return TriaIterator > > + ( + tria, + present_level, + line_index (i), + mg_dof_handler + ); +}; + + + +template +TriaIterator > > +MGDoFHexAccessor::quad (const unsigned int i) const { + Assert (i<12, ExcInvalidIndex (i, 0, 6)); + + return TriaIterator > > + ( + tria, + present_level, + quad_index (i), + mg_dof_handler + ); +}; + + + +template +TriaIterator > +MGDoFHexAccessor::child (const unsigned int i) const { + TriaIterator > q (tria, + present_level+1, + child_index (i), + mg_dof_handler); + +#ifdef DEBUG + if (q.state() != past_the_end) + Assert (q->used(), typename TriaAccessor::ExcUnusedCellAsChild()); +#endif + return q; +}; + + + +template +void +MGDoFHexAccessor::copy_from (const MGDoFHexAccessor &a) { + DoFHexAccessor::copy_from (a); + set_mg_dof_handler (a.mg_dof_handler); +}; + + + + /*------------------------- Functions: MGDoFCellAccessor -----------------------*/ @@ -440,3 +611,23 @@ template class TriaActiveIterator<2,MGDoFCellAccessor<2> >; #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 + + diff --git a/deal.II/deal.II/source/dofs/mg_dof_handler.cc b/deal.II/deal.II/source/dofs/mg_dof_handler.cc index 2d6f4ffaf4..26c8d83933 100644 --- a/deal.II/deal.II/source/dofs/mg_dof_handler.cc +++ b/deal.II/deal.II/source/dofs/mg_dof_handler.cc @@ -353,6 +353,94 @@ MGDoFHandler<1>::last_active_quad () const { 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 @@ -516,6 +604,261 @@ MGDoFHandler<2>::last_active_face (const unsigned int level) const { 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 @@ -588,6 +931,39 @@ MGDoFHandler::begin_active_quad (const unsigned int level) const { +template +typename MGDoFHandler::raw_hex_iterator +MGDoFHandler::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 +typename MGDoFHandler::hex_iterator +MGDoFHandler::begin_hex (const unsigned int level) const { + return hex_iterator (tria, + tria->begin_hex(level)->level(), + tria->begin_hex(level)->index(), + this); +}; + + + +template +typename MGDoFHandler::active_hex_iterator +MGDoFHandler::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 typename MGDoFHandler::raw_line_iterator MGDoFHandler::end_line () const { @@ -604,6 +980,14 @@ MGDoFHandler::end_quad () const { +template +typename MGDoFHandler::raw_hex_iterator +MGDoFHandler::end_hex () const { + return raw_hex_iterator (tria, -1, -1, this); +}; + + + template typename MGDoFDimensionInfo::raw_cell_iterator MGDoFHandler::end_raw (const unsigned int level) const { @@ -694,34 +1078,65 @@ MGDoFHandler::end_active_line (const unsigned int level) const { - + +template +typename MGDoFDimensionInfo::raw_quad_iterator +MGDoFHandler::end_raw_quad (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + end_quad() : + begin_raw_quad (level+1)); +}; + + + +template +typename MGDoFDimensionInfo::quad_iterator +MGDoFHandler::end_quad (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + quad_iterator(end_quad()) : + begin_quad (level+1)); +}; + + + + +template +typename MGDoFDimensionInfo::active_quad_iterator +MGDoFHandler::end_active_quad (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + active_quad_iterator(end_quad()) : + begin_active_quad (level+1)); +}; + + + template -typename MGDoFDimensionInfo::raw_quad_iterator -MGDoFHandler::end_raw_quad (const unsigned int level) const { +typename MGDoFDimensionInfo::raw_hex_iterator +MGDoFHandler::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 -typename MGDoFDimensionInfo::quad_iterator -MGDoFHandler::end_quad (const unsigned int level) const { +typename MGDoFDimensionInfo::hex_iterator +MGDoFHandler::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 -typename MGDoFDimensionInfo::active_quad_iterator -MGDoFHandler::end_active_quad (const unsigned int level) const { +typename MGDoFDimensionInfo::active_hex_iterator +MGDoFHandler::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)); }; @@ -795,6 +1210,42 @@ MGDoFHandler::last_active_quad (const unsigned int level) const { +template +typename MGDoFHandler::raw_hex_iterator +MGDoFHandler::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 +typename MGDoFHandler::hex_iterator +MGDoFHandler::last_hex (const unsigned int level) const { + return hex_iterator (tria, + tria->last_hex(level)->level(), + tria->last_hex(level)->index(), + this); +}; + + + + +template +typename MGDoFHandler::active_hex_iterator +MGDoFHandler::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 typename MGDoFHandler::raw_line_iterator MGDoFHandler::last_raw_line () const { @@ -811,6 +1262,14 @@ MGDoFHandler::last_raw_quad () const { +template +typename MGDoFHandler::raw_hex_iterator +MGDoFHandler::last_raw_hex () const { + return last_raw_hex (mg_levels.size()-1); +}; + + + template typename MGDoFHandler::line_iterator MGDoFHandler::last_line () const { @@ -827,6 +1286,14 @@ MGDoFHandler::last_quad () const { +template +typename MGDoFHandler::hex_iterator +MGDoFHandler::last_hex () const { + return last_hex (mg_levels.size()-1); +}; + + + template typename MGDoFHandler::active_line_iterator MGDoFHandler::last_active_line () const { @@ -844,6 +1311,14 @@ MGDoFHandler::last_active_quad () const { +template +typename MGDoFHandler::active_hex_iterator +MGDoFHandler::last_active_hex () const { + return last_active_hex (mg_levels.size()-1); +}; + + + @@ -989,6 +1464,69 @@ MGDoFHandler<2>::distribute_dofs_on_cell (cell_iterator &cell, +#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::vertices_per_cell; ++vertex) + // check whether dofs for this + // vertex have been distributed + // (only check the first dof) + if (cell->mg_vertex_dof_index(vertex, 0) == -1) + for (unsigned int d=0; ddofs_per_vertex; ++d) + cell->set_mg_vertex_dof_index (vertex, d, next_free_dof++); + + // for the lines + if (selected_fe->dofs_per_line > 0) + for (unsigned int l=0; l::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; ddofs_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::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; ddofs_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; ddofs_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 unsigned int MGDoFHandler::n_dofs (const unsigned int level) const { @@ -1025,182 +1563,11 @@ void MGDoFHandler::make_sparsity_pattern (const unsigned int level, - -template -void MGDoFHandler::renumber_dofs (const unsigned int level, - const RenumberingMethod method, - const vector &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 new_number(n_dofs, -1); - - // store the indices of the dofs renumbered - // in the last round. Default to starting - // points - vector last_round_dofs (starting_points); - - // delete disallowed elements - for (unsigned int i=0; i=n_dofs)) - last_round_dofs[i] = -1; - - remove_if (last_round_dofs.begin(), last_round_dofs.end(), - bind2nd(equal_to(), -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 next_round_dofs; - - // find all neighbors of the - // dofs numbered in the last - // round - for (unsigned int i=0; i::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 dofs_by_coordination; - - // find coordination number for - // each of these dofs - for (vector::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 new_entry (coordination, *s); - dofs_by_coordination.insert (new_entry); - }; - - //// - multimap::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::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 &new_numbers) { +void MGDoFHandler<1>::renumber_dofs (const unsigned int level, + const vector &new_numbers) { Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete()); // note that we can not use cell iterators @@ -1236,8 +1603,45 @@ void MGDoFHandler<1>::do_renumbering (const unsigned int level, #if deal_II_dimension == 2 template <> -void MGDoFHandler<2>::do_renumbering (const unsigned int level, - const vector &new_numbers) { +void MGDoFHandler<2>::renumber_dofs (const unsigned int level, + const vector &new_numbers) { + Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete()); + + for (vector::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; ddofs_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::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::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 &new_numbers) { Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete()); for (vector::iterator i=mg_vertex_dofs.begin(); @@ -1257,12 +1661,20 @@ void MGDoFHandler<2>::do_renumbering (const unsigned int level, Assert (*i != -1, ExcInternalError()); *i = new_numbers[*i]; }; + for (vector::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::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 @@ -1445,5 +1857,100 @@ void MGDoFHandler<2>::reserve_space () { #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; in_levels(); ++i) + { + mg_levels.push_back (new DoFLevel<3>); + + mg_levels.back()->line_dofs = vector (tria->levels[i]->lines.lines.size() * + selected_fe->dofs_per_line, + -1); + mg_levels.back()->quad_dofs = vector (tria->levels[i]->quads.quads.size() * + selected_fe->dofs_per_quad, + -1); + mg_levels.back()->hex_dofs = vector (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 min_level (tria->vertices.size(), tria->n_levels()); + vector max_level (tria->vertices.size(), 0); + + Triangulation::cell_iterator cell = tria->begin(), + endc = tria->end(); + for (; cell!=endc; ++cell) + for (unsigned int vertex=0; vertex::vertices_per_cell; + ++vertex) + { + const unsigned int vertex_index = cell->vertex_index(vertex); + if (min_level[vertex_index] > static_cast(cell->level())) + min_level[vertex_index] = cell->level(); + if (max_level[vertex_index] < static_cast(cell->level())) + max_level[vertex_index] = cell->level(); + }; + + + // now allocate the needed space + for (unsigned int vertex=0; vertexvertices.size(); ++vertex) + { + Assert (min_level[vertex] < tria->n_levels(), ExcInternalError()); + Assert (max_level[vertex] >= min_level[vertex], ExcInternalError()); + + mg_vertex_dofs[vertex].init (min_level[vertex], + max_level[vertex], + selected_fe->dofs_per_vertex); + }; +}; + +#endif + + // explicite instantiations template class MGDoFHandler; diff --git a/deal.II/deal.II/source/multigrid/mg_dof_accessor.cc b/deal.II/deal.II/source/multigrid/mg_dof_accessor.cc index 101f9344fe..d6fd3b0b24 100644 --- a/deal.II/deal.II/source/multigrid/mg_dof_accessor.cc +++ b/deal.II/deal.II/source/multigrid/mg_dof_accessor.cc @@ -299,6 +299,177 @@ MGDoFQuadAccessor::copy_from (const MGDoFQuadAccessor +MGDoFHexAccessor::MGDoFHexAccessor (Triangulation *tria, + const int level, + const int index, + const AccessorData *local_data) : + MGDoFAccessor (local_data), + DoFHexAccessor(tria,level,index,local_data) {}; + + + +template +inline +int MGDoFHexAccessor::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 (iget_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 +void MGDoFHexAccessor::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 (iget_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 +inline +int MGDoFHexAccessor::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 (iget_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 +void MGDoFHexAccessor::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 (iget_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 +void +MGDoFHexAccessor::get_mg_dof_indices (vector &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::iterator next = dof_indices.begin(); + for (unsigned int vertex=0; vertex<8; ++vertex) + for (unsigned int d=0; dline(line)->mg_dof_index(d); + for (unsigned int quad=0; quad<12; ++quad) + for (unsigned int d=0; dquad(quad)->mg_dof_index(d); + for (unsigned int d=0; d +TriaIterator > > +MGDoFHexAccessor::line (const unsigned int i) const { + Assert (i<12, ExcInvalidIndex (i, 0, 12)); + + return TriaIterator > > + ( + tria, + present_level, + line_index (i), + mg_dof_handler + ); +}; + + + +template +TriaIterator > > +MGDoFHexAccessor::quad (const unsigned int i) const { + Assert (i<12, ExcInvalidIndex (i, 0, 6)); + + return TriaIterator > > + ( + tria, + present_level, + quad_index (i), + mg_dof_handler + ); +}; + + + +template +TriaIterator > +MGDoFHexAccessor::child (const unsigned int i) const { + TriaIterator > q (tria, + present_level+1, + child_index (i), + mg_dof_handler); + +#ifdef DEBUG + if (q.state() != past_the_end) + Assert (q->used(), typename TriaAccessor::ExcUnusedCellAsChild()); +#endif + return q; +}; + + + +template +void +MGDoFHexAccessor::copy_from (const MGDoFHexAccessor &a) { + DoFHexAccessor::copy_from (a); + set_mg_dof_handler (a.mg_dof_handler); +}; + + + + /*------------------------- Functions: MGDoFCellAccessor -----------------------*/ @@ -440,3 +611,23 @@ template class TriaActiveIterator<2,MGDoFCellAccessor<2> >; #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 + + diff --git a/deal.II/deal.II/source/multigrid/mg_dof_handler.cc b/deal.II/deal.II/source/multigrid/mg_dof_handler.cc index 2d6f4ffaf4..26c8d83933 100644 --- a/deal.II/deal.II/source/multigrid/mg_dof_handler.cc +++ b/deal.II/deal.II/source/multigrid/mg_dof_handler.cc @@ -353,6 +353,94 @@ MGDoFHandler<1>::last_active_quad () const { 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 @@ -516,6 +604,261 @@ MGDoFHandler<2>::last_active_face (const unsigned int level) const { 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 @@ -588,6 +931,39 @@ MGDoFHandler::begin_active_quad (const unsigned int level) const { +template +typename MGDoFHandler::raw_hex_iterator +MGDoFHandler::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 +typename MGDoFHandler::hex_iterator +MGDoFHandler::begin_hex (const unsigned int level) const { + return hex_iterator (tria, + tria->begin_hex(level)->level(), + tria->begin_hex(level)->index(), + this); +}; + + + +template +typename MGDoFHandler::active_hex_iterator +MGDoFHandler::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 typename MGDoFHandler::raw_line_iterator MGDoFHandler::end_line () const { @@ -604,6 +980,14 @@ MGDoFHandler::end_quad () const { +template +typename MGDoFHandler::raw_hex_iterator +MGDoFHandler::end_hex () const { + return raw_hex_iterator (tria, -1, -1, this); +}; + + + template typename MGDoFDimensionInfo::raw_cell_iterator MGDoFHandler::end_raw (const unsigned int level) const { @@ -694,34 +1078,65 @@ MGDoFHandler::end_active_line (const unsigned int level) const { - + +template +typename MGDoFDimensionInfo::raw_quad_iterator +MGDoFHandler::end_raw_quad (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + end_quad() : + begin_raw_quad (level+1)); +}; + + + +template +typename MGDoFDimensionInfo::quad_iterator +MGDoFHandler::end_quad (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + quad_iterator(end_quad()) : + begin_quad (level+1)); +}; + + + + +template +typename MGDoFDimensionInfo::active_quad_iterator +MGDoFHandler::end_active_quad (const unsigned int level) const { + return (level == mg_levels.size()-1 ? + active_quad_iterator(end_quad()) : + begin_active_quad (level+1)); +}; + + + template -typename MGDoFDimensionInfo::raw_quad_iterator -MGDoFHandler::end_raw_quad (const unsigned int level) const { +typename MGDoFDimensionInfo::raw_hex_iterator +MGDoFHandler::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 -typename MGDoFDimensionInfo::quad_iterator -MGDoFHandler::end_quad (const unsigned int level) const { +typename MGDoFDimensionInfo::hex_iterator +MGDoFHandler::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 -typename MGDoFDimensionInfo::active_quad_iterator -MGDoFHandler::end_active_quad (const unsigned int level) const { +typename MGDoFDimensionInfo::active_hex_iterator +MGDoFHandler::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)); }; @@ -795,6 +1210,42 @@ MGDoFHandler::last_active_quad (const unsigned int level) const { +template +typename MGDoFHandler::raw_hex_iterator +MGDoFHandler::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 +typename MGDoFHandler::hex_iterator +MGDoFHandler::last_hex (const unsigned int level) const { + return hex_iterator (tria, + tria->last_hex(level)->level(), + tria->last_hex(level)->index(), + this); +}; + + + + +template +typename MGDoFHandler::active_hex_iterator +MGDoFHandler::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 typename MGDoFHandler::raw_line_iterator MGDoFHandler::last_raw_line () const { @@ -811,6 +1262,14 @@ MGDoFHandler::last_raw_quad () const { +template +typename MGDoFHandler::raw_hex_iterator +MGDoFHandler::last_raw_hex () const { + return last_raw_hex (mg_levels.size()-1); +}; + + + template typename MGDoFHandler::line_iterator MGDoFHandler::last_line () const { @@ -827,6 +1286,14 @@ MGDoFHandler::last_quad () const { +template +typename MGDoFHandler::hex_iterator +MGDoFHandler::last_hex () const { + return last_hex (mg_levels.size()-1); +}; + + + template typename MGDoFHandler::active_line_iterator MGDoFHandler::last_active_line () const { @@ -844,6 +1311,14 @@ MGDoFHandler::last_active_quad () const { +template +typename MGDoFHandler::active_hex_iterator +MGDoFHandler::last_active_hex () const { + return last_active_hex (mg_levels.size()-1); +}; + + + @@ -989,6 +1464,69 @@ MGDoFHandler<2>::distribute_dofs_on_cell (cell_iterator &cell, +#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::vertices_per_cell; ++vertex) + // check whether dofs for this + // vertex have been distributed + // (only check the first dof) + if (cell->mg_vertex_dof_index(vertex, 0) == -1) + for (unsigned int d=0; ddofs_per_vertex; ++d) + cell->set_mg_vertex_dof_index (vertex, d, next_free_dof++); + + // for the lines + if (selected_fe->dofs_per_line > 0) + for (unsigned int l=0; l::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; ddofs_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::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; ddofs_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; ddofs_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 unsigned int MGDoFHandler::n_dofs (const unsigned int level) const { @@ -1025,182 +1563,11 @@ void MGDoFHandler::make_sparsity_pattern (const unsigned int level, - -template -void MGDoFHandler::renumber_dofs (const unsigned int level, - const RenumberingMethod method, - const vector &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 new_number(n_dofs, -1); - - // store the indices of the dofs renumbered - // in the last round. Default to starting - // points - vector last_round_dofs (starting_points); - - // delete disallowed elements - for (unsigned int i=0; i=n_dofs)) - last_round_dofs[i] = -1; - - remove_if (last_round_dofs.begin(), last_round_dofs.end(), - bind2nd(equal_to(), -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 next_round_dofs; - - // find all neighbors of the - // dofs numbered in the last - // round - for (unsigned int i=0; i::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 dofs_by_coordination; - - // find coordination number for - // each of these dofs - for (vector::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 new_entry (coordination, *s); - dofs_by_coordination.insert (new_entry); - }; - - //// - multimap::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::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 &new_numbers) { +void MGDoFHandler<1>::renumber_dofs (const unsigned int level, + const vector &new_numbers) { Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete()); // note that we can not use cell iterators @@ -1236,8 +1603,45 @@ void MGDoFHandler<1>::do_renumbering (const unsigned int level, #if deal_II_dimension == 2 template <> -void MGDoFHandler<2>::do_renumbering (const unsigned int level, - const vector &new_numbers) { +void MGDoFHandler<2>::renumber_dofs (const unsigned int level, + const vector &new_numbers) { + Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete()); + + for (vector::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; ddofs_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::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::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 &new_numbers) { Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete()); for (vector::iterator i=mg_vertex_dofs.begin(); @@ -1257,12 +1661,20 @@ void MGDoFHandler<2>::do_renumbering (const unsigned int level, Assert (*i != -1, ExcInternalError()); *i = new_numbers[*i]; }; + for (vector::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::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 @@ -1445,5 +1857,100 @@ void MGDoFHandler<2>::reserve_space () { #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; in_levels(); ++i) + { + mg_levels.push_back (new DoFLevel<3>); + + mg_levels.back()->line_dofs = vector (tria->levels[i]->lines.lines.size() * + selected_fe->dofs_per_line, + -1); + mg_levels.back()->quad_dofs = vector (tria->levels[i]->quads.quads.size() * + selected_fe->dofs_per_quad, + -1); + mg_levels.back()->hex_dofs = vector (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 min_level (tria->vertices.size(), tria->n_levels()); + vector max_level (tria->vertices.size(), 0); + + Triangulation::cell_iterator cell = tria->begin(), + endc = tria->end(); + for (; cell!=endc; ++cell) + for (unsigned int vertex=0; vertex::vertices_per_cell; + ++vertex) + { + const unsigned int vertex_index = cell->vertex_index(vertex); + if (min_level[vertex_index] > static_cast(cell->level())) + min_level[vertex_index] = cell->level(); + if (max_level[vertex_index] < static_cast(cell->level())) + max_level[vertex_index] = cell->level(); + }; + + + // now allocate the needed space + for (unsigned int vertex=0; vertexvertices.size(); ++vertex) + { + Assert (min_level[vertex] < tria->n_levels(), ExcInternalError()); + Assert (max_level[vertex] >= min_level[vertex], ExcInternalError()); + + mg_vertex_dofs[vertex].init (min_level[vertex], + max_level[vertex], + selected_fe->dofs_per_vertex); + }; +}; + +#endif + + // explicite instantiations template class MGDoFHandler; diff --git a/tests/big-tests/dof/dof_test.cc b/tests/big-tests/dof/dof_test.cc index cc8a9f83f5..8ca5178abf 100644 --- a/tests/big-tests/dof/dof_test.cc +++ b/tests/big-tests/dof/dof_test.cc @@ -12,7 +12,7 @@ #include #include #include - +#include #include #include @@ -354,7 +354,7 @@ void TestCases::run (ParameterHandler &prm) { 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()); diff --git a/tests/deal.II/dof_test.cc b/tests/deal.II/dof_test.cc index f73b99509b..e5f74da243 100644 --- a/tests/deal.II/dof_test.cc +++ b/tests/deal.II/dof_test.cc @@ -14,7 +14,7 @@ #include #include #include - +#include #include #include @@ -287,7 +287,7 @@ void TestCases::run (const unsigned int test_case) { 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()); -- 2.39.5