// $Id$
// Version: $Name$
//
-// Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors
+// Copyright (C) 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
* published in 1984), but certainly not with those used in this
* library, featuring several 10,000 to a few 100,000 elements.
*
- *
+ *
* <h4>Implementation of renumbering schemes</h4>
*
* The renumbering algorithms need quite a lot of memory, since they
* 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
* deleted. If no index is allowed, the algorithm will search for its
* own starting point.
*
- *
+ *
* <h4>Results of renumbering</h4>
*
* The renumbering schemes mentioned above do not lead to optimal
* @ingroup dofs
* @author Wolfgang Bangerth, Guido Kanschat, 1998, 1999, 2000, 2004, 2007, 2008
*/
-namespace DoFRenumbering
+namespace DoFRenumbering
{
/**
* Direction based comparator for
*/
CompareDownstream (const Point<dim> &dir)
:
- dir(dir)
+ dir(dir)
{}
/**
* Return true if c1 less c2.
const Point<dim> diff = c2->center() - c1->center();
return (diff*dir > 0);
}
-
+
+ private:
+ /**
+ * Flow direction.
+ */
+ const Point<dim> dir;
+ };
+
+
+ /**
+ * Point based comparator for downstream
+ * directions: it returns @p true if the
+ * second point is downstream of the first
+ * one with respect to the direction given
+ * to the constructor. If the points are
+ * the same with respect to the downstream
+ * direction, the point with the lower DoF
+ * number is considered smaller.
+ */
+ template <int dim>
+ struct ComparePointwiseDownstream
+ {
+ /**
+ * Constructor.
+ */
+ ComparePointwiseDownstream (const Point<dim> &dir)
+ :
+ dir(dir)
+ {}
+ /**
+ * Return true if c1 less c2.
+ */
+ bool operator () (const std::pair<Point<dim>,unsigned int> &c1,
+ const std::pair<Point<dim>,unsigned int> &c2) const
+ {
+ const Point<dim> diff = c2.first-c1.first;
+ return (diff*dir > 0 || (diff*dir==0 && c1.second<c2.second));
+ }
+
private:
/**
* Flow direction.
*/
const Point<dim> dir;
};
-
+
/**
* A namespace for the
* implementation of some
* DoFRenumbering namespace.
*/
template <class DH>
- void
+ void
Cuthill_McKee (DH& dof_handler,
const bool reversed_numbering = false,
const bool use_constraints = false);
* the DoFHandler dofs but
* returns the renumbering
* vector.
- */
+ */
template <class DH>
void
compute_Cuthill_McKee (std::vector<unsigned int>& new_dof_indices,
* step-22.
*/
template <class DH>
- void
+ void
king_ordering (DH& dof_handler,
const bool reversed_numbering = false,
const bool use_constraints = false);
* DoFRenumbering namespace.
*/
template <class DH>
- void
+ void
minimum_degree (DH& dof_handler,
const bool reversed_numbering = false,
const bool use_constraints = false);
const bool reversed_numbering = false,
const bool use_constraints = false);
}
-
+
/**
* Renumber the degrees of
* freedom according to the
* DoFRenumbering namespace.
*/
template <class DH>
- void
+ void
Cuthill_McKee (DH& dof_handler,
const bool reversed_numbering = false,
const bool use_constraints = false,
* the DoFHandler dofs but
* returns the renumbering
* vector.
- */
+ */
template <class DH>
void
compute_Cuthill_McKee (std::vector<unsigned int>& new_dof_indices,
* the different methods.
*/
template <int dim>
- void
+ void
Cuthill_McKee (MGDoFHandler<dim> &dof_handler,
const unsigned int level,
const bool reversed_numbering = false,
* are associated with the first
* vector component to which they
* belong.
- *
+ *
* For finite elements with only
* one component, or a single
* non-primitive base element,
void
component_wise (MGDoFHandler<dim>& dof_handler,
const std::vector<unsigned int>& target_component = std::vector<unsigned int>());
-
+
/**
* Computes the renumbering
* vector needed by the
* the renumbering on the
* DoFHandler dofs but returns
* the renumbering vector.
- */
+ */
template <int dim, class ITERATOR, class ENDITERATOR>
static unsigned int
compute_component_wise (std::vector<unsigned int>& new_dof_indices,
cell_wise (MGDoFHandler<dim> &dof_handler,
const unsigned int level,
const std::vector<typename MGDoFHandler<dim>::cell_iterator> &cell_order);
-
+
/**
* @deprecated Use cell_wise() instead.
*
cell_wise_dg (MGDoFHandler<dim> &dof_handler,
const unsigned int level,
const std::vector<typename MGDoFHandler<dim>::cell_iterator> &cell_order);
-
+
/**
* Computes the renumbering
* vector needed by the
/**
- * Cell-wise downstream numbering
- * with respect to a constant
- * flow direction.
+ * Downstream numbering with respect to a
+ * constant flow direction. If the
+ * additional argument @p
+ * dof_wise_renumbering is set to @p false,
+ * the numbering is performed cell-wise. If
+ * it is set to @p true, the dofs are
+ * renumbered for finite elements that are
+ * based on support points. Note that in
+ * this case, the numbering of points will
+ * be unaffected in case two points are
+ * located the same with respect to the
+ * downstream direction.
*
* The cells are sorted such that
* the centers of higher numbers
template <class DH, int dim>
void
downstream (DH& dof_handler,
- const Point<dim>& direction);
+ const Point<dim>& direction,
+ const bool dof_wise_renumbering = false);
/**
* @deprecated Use downstream() instead.
*/
template <int dim>
void
- downstream (MGDoFHandler<dim> &dof_handler,
+ downstream (MGDoFHandler<dim> &dof_handler,
const unsigned int level,
- const Point<dim> &direction);
+ const Point<dim> &direction,
+ const bool dof_wise_renumbering = false);
+
/**
* @deprecated Use downstream()
* instead.
*/
template <int dim>
void
- downstream_dg (MGDoFHandler<dim> &dof_handler,
+ downstream_dg (MGDoFHandler<dim> &dof_handler,
const unsigned int level,
- const Point<dim> &direction);
+ const Point<dim> &direction);
/**
* @deprecated The new function
compute_downstream (std::vector<unsigned int>& new_dof_indices,
std::vector<unsigned int>& reverse,
const DH& dof_handler,
- const Point<dim>& direction);
-
+ const Point<dim>& direction,
+ const bool dof_wise_renumbering);
+
/**
* @deprecated Use
* compute_downstream() instead
std::vector<unsigned int>& reverse,
const MGDoFHandler<dim>& dof_handler,
const unsigned int level,
- const Point<dim>& direction);
-
+ const Point<dim>& direction,
+ const bool dof_wise_renumbering);
+
/**
* @deprecated Use
* compute_downstream() instead
* renumbering on the DoFHandler
* dofs but returns the
* renumbering vector.
- */
+ */
template <class DH>
void
compute_random (std::vector<unsigned int> &new_dof_indices,
* renumbering on the @p
* DoFHandler dofs but returns
* the renumbering vector.
- */
+ */
template <class DH>
void
compute_subdomain_wise (std::vector<unsigned int> &new_dof_indices,
const DH &dof_handler);
-
+
/**
* Exception
*
// $id$
// Version: $name$
//
-// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors
+// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
#include <base/thread_management.h>
#include <base/utilities.h>
+#include <base/quadrature_lib.h>
+
#include <lac/sparsity_pattern.h>
#include <lac/sparsity_tools.h>
#include <lac/compressed_simple_sparsity_pattern.h>
#include <lac/constraint_matrix.h>
+
#include <dofs/dof_accessor.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_tools.h>
+#include <dofs/dof_renumbering.h>
+
#include <grid/tria_iterator.h>
#include <grid/tria.h>
-#include <dofs/dof_handler.h>
+
+#include <fe/fe.h>
#include <hp/dof_handler.h>
#include <hp/fe_collection.h>
-#include <dofs/dof_tools.h>
-#include <fe/fe.h>
-#include <dofs/dof_renumbering.h>
+#include <hp/fe_values.h>
#include <multigrid/mg_dof_handler.h>
#include <multigrid/mg_dof_accessor.h>
template <class DH, int dim>
void
- downstream (DH& dof, const Point<dim>& direction)
+ downstream (DH& dof, const Point<dim>& direction,
+ const bool dof_wise_renumbering)
{
std::vector<unsigned int> renumbering(dof.n_dofs());
std::vector<unsigned int> reverse(dof.n_dofs());
- compute_downstream(renumbering, reverse, dof, direction);
+ compute_downstream(renumbering, reverse, dof, direction,
+ dof_wise_renumbering);
dof.renumber_dofs(renumbering);
}
std::vector<unsigned int>& new_indices,
std::vector<unsigned int>& reverse,
const DH& dof,
- const Point<dim>& direction)
+ const Point<dim>& direction,
+ const bool dof_wise_renumbering)
{
- std::vector<typename DH::cell_iterator>
- ordered_cells(dof.get_tria().n_active_cells());
- const CompareDownstream<typename DH::cell_iterator, dim> comparator(direction);
+ if (dof_wise_renumbering == false)
+ {
+ std::vector<typename DH::cell_iterator>
+ ordered_cells(dof.get_tria().n_active_cells());
+ const CompareDownstream<typename DH::cell_iterator, dim> comparator(direction);
- typename DH::active_cell_iterator begin = dof.begin_active();
- typename DH::active_cell_iterator end = dof.end();
+ typename DH::active_cell_iterator begin = dof.begin_active();
+ typename DH::active_cell_iterator end = dof.end();
- copy (begin, end, ordered_cells.begin());
- std::sort (ordered_cells.begin(), ordered_cells.end(), comparator);
+ copy (begin, end, ordered_cells.begin());
+ std::sort (ordered_cells.begin(), ordered_cells.end(), comparator);
+
+ compute_cell_wise(new_indices, reverse, dof, ordered_cells);
+ }
+ else
+ {
+ // similar code as for
+ // DoFTools::map_dofs_to_support_points, but
+ // need to do this for general DH classes and
+ // want to be able to sort the result
+ // (otherwise, could use something like
+ // DoFTools::map_support_points_to_dofs)
+ const unsigned int n_dofs = dof.n_dofs();
+ std::vector<std::pair<Point<dim>,unsigned int> > support_point_list
+ (n_dofs);
+
+ const hp::FECollection<dim> fe_collection (dof.get_fe ());
+ Assert (fe_collection[0].has_support_points(),
+ DoFTools::ExcFEHasNoSupportPoints());
+ hp::QCollection<dim> quadrature_collection;
+ for (unsigned int comp=0; comp<fe_collection.size(); ++comp)
+ {
+ Assert (fe_collection[comp].has_support_points(),
+ DoFTools::ExcFEHasNoSupportPoints());
+ quadrature_collection.push_back
+ (Quadrature<DH::dimension> (fe_collection[comp].
+ get_unit_support_points()));
+ }
+ hp::FEValues<DH::dimension,DH::space_dimension>
+ hp_fe_values (fe_collection, quadrature_collection,
+ update_quadrature_points);
+
+ std::vector<bool> already_touched (n_dofs, false);
+
+ std::vector<unsigned int> local_dof_indices;
+ typename DH::active_cell_iterator begin = dof.begin_active();
+ typename DH::active_cell_iterator end = dof.end();
+ for ( ; begin != end; ++begin)
+ {
+ const unsigned int dofs_per_cell = begin->get_fe().dofs_per_cell;
+ local_dof_indices.resize (dofs_per_cell);
+ hp_fe_values.reinit (begin);
+ const FEValues<dim> &fe_values =
+ hp_fe_values.get_present_fe_values ();
+ begin->get_dof_indices(local_dof_indices);
+ const std::vector<Point<DH::space_dimension> > & points
+ = fe_values.get_quadrature_points ();
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ if (!already_touched[local_dof_indices[i]])
+ {
+ support_point_list[local_dof_indices[i]].first = points[i];
+ support_point_list[local_dof_indices[i]].second =
+ local_dof_indices[i];
+ already_touched[local_dof_indices[i]] = true;
+ }
+ }
- compute_cell_wise(new_indices, reverse, dof, ordered_cells);
+ ComparePointwiseDownstream<dim> comparator (direction);
+ std::sort (support_point_list.begin(), support_point_list.end(),
+ comparator);
+ for (unsigned int i=0; i<n_dofs; ++i)
+ new_indices[support_point_list[i].second] = i;
+ }
}
template <int dim>
void downstream_dg (MGDoFHandler<dim>& dof,
const unsigned int level,
- const Point<dim>& direction)
+ const Point<dim>& direction)
{
std::vector<unsigned int> renumbering(dof.n_dofs(level));
std::vector<unsigned int> reverse(dof.n_dofs(level));
template <int dim>
void downstream (MGDoFHandler<dim>& dof,
const unsigned int level,
- const Point<dim>& direction)
+ const Point<dim>& direction,
+ const bool dof_wise_renumbering)
{
std::vector<unsigned int> renumbering(dof.n_dofs(level));
std::vector<unsigned int> reverse(dof.n_dofs(level));
- compute_downstream(renumbering, reverse, dof, level, direction);
+ compute_downstream(renumbering, reverse, dof, level, direction,
+ dof_wise_renumbering);
dof.renumber_dofs(level, renumbering);
}
{
std::vector<typename MGDoFHandler<dim>::cell_iterator>
ordered_cells(dof.get_tria().n_cells(level));
- const CompareDownstream<typename MGDoFHandler<dim>::cell_iterator, dim> comparator(direction);
+ const CompareDownstream<typename MGDoFHandler<dim>::cell_iterator, dim>
+ comparator(direction);
typename MGDoFHandler<dim>::cell_iterator begin = dof.begin(level);
typename MGDoFHandler<dim>::cell_iterator end = dof.end(level);
std::vector<unsigned int>& reverse,
const MGDoFHandler<dim>& dof,
const unsigned int level,
- const Point<dim>& direction)
+ const Point<dim>& direction,
+ const bool dof_wise_renumbering)
{
- std::vector<typename MGDoFHandler<dim>::cell_iterator>
- ordered_cells(dof.get_tria().n_cells(level));
- const CompareDownstream<typename MGDoFHandler<dim>::cell_iterator, dim> comparator(direction);
+ if (dof_wise_renumbering == false)
+ {
+ std::vector<typename MGDoFHandler<dim>::cell_iterator>
+ ordered_cells(dof.get_tria().n_cells(level));
+ const CompareDownstream<typename MGDoFHandler<dim>::cell_iterator, dim> comparator(direction);
- typename MGDoFHandler<dim>::cell_iterator begin = dof.begin(level);
- typename MGDoFHandler<dim>::cell_iterator end = dof.end(level);
+ typename MGDoFHandler<dim>::cell_iterator begin = dof.begin(level);
+ typename MGDoFHandler<dim>::cell_iterator end = dof.end(level);
- std::copy (begin, end, ordered_cells.begin());
- std::sort (ordered_cells.begin(), ordered_cells.end(), comparator);
+ std::copy (begin, end, ordered_cells.begin());
+ std::sort (ordered_cells.begin(), ordered_cells.end(), comparator);
- compute_cell_wise(new_indices, reverse, dof, level, ordered_cells);
+ compute_cell_wise(new_indices, reverse, dof, level, ordered_cells);
+ }
+ else
+ {
+ Assert (dof.get_fe().has_support_points(),
+ DoFTools::ExcFEHasNoSupportPoints());
+ const unsigned int n_dofs = dof.n_dofs(level);
+ std::vector<std::pair<Point<dim>,unsigned int> > support_point_list
+ (n_dofs);
+
+ Quadrature<dim> q_dummy(dof.get_fe().get_unit_support_points());
+ FEValues<dim,dim> fe_values (dof.get_fe(), q_dummy,
+ update_quadrature_points);
+
+ std::vector<bool> already_touched (dof.n_dofs(), false);
+
+ const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell;
+ std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+ typename MGDoFHandler<dim>::cell_iterator begin = dof.begin(level);
+ typename MGDoFHandler<dim>::cell_iterator end = dof.end(level);
+ for ( ; begin != end; ++begin)
+ {
+ begin->get_mg_dof_indices(local_dof_indices);
+ fe_values.reinit (begin);
+ const std::vector<Point<dim> > & points
+ = fe_values.get_quadrature_points ();
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ if (!already_touched[local_dof_indices[i]])
+ {
+ support_point_list[local_dof_indices[i]].first = points[i];
+ support_point_list[local_dof_indices[i]].second =
+ local_dof_indices[i];
+ already_touched[local_dof_indices[i]] = true;
+ }
+ }
+
+ ComparePointwiseDownstream<dim> comparator (direction);
+ std::sort (support_point_list.begin(), support_point_list.end(),
+ comparator);
+ for (unsigned int i=0; i<n_dofs; ++i)
+ new_indices[support_point_list[i].second] = i;
+ }
}
template void
downstream<DoFHandler<deal_II_dimension> >
- (DoFHandler<deal_II_dimension>&, const Point<deal_II_dimension>&);
+ (DoFHandler<deal_II_dimension>&, const Point<deal_II_dimension>&,
+ const bool);
template void
compute_downstream<DoFHandler<deal_II_dimension> >
(std::vector<unsigned int>&,std::vector<unsigned int>&,
- const DoFHandler<deal_II_dimension>&, const Point<deal_II_dimension>&);
+ const DoFHandler<deal_II_dimension>&, const Point<deal_II_dimension>&,
+ const bool);
template
void
template void
downstream<hp::DoFHandler<deal_II_dimension> >
(hp::DoFHandler<deal_II_dimension>&,
- const Point<deal_II_dimension>&);
+ const Point<deal_II_dimension>&,
+ const bool);
template void
compute_downstream<hp::DoFHandler<deal_II_dimension> >
(std::vector<unsigned int>&,std::vector<unsigned int>&,
const hp::DoFHandler<deal_II_dimension>&,
- const Point<deal_II_dimension>&);
+ const Point<deal_II_dimension>&,
+ const bool);
template
void
template
void downstream<deal_II_dimension>
(MGDoFHandler<deal_II_dimension>&, const unsigned int,
- const Point<deal_II_dimension>&);
+ const Point<deal_II_dimension>&,
+ const bool);
template
void clockwise_dg<deal_II_dimension>