#include <grid/tria.h>
#include <dofs/dof_handler.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>
T::operator==;
};
}
-
+
namespace boost
{
namespace types
using namespace ::boost;
using namespace std;
- typedef adjacency_list<vecS, vecS, undirectedS,
+ typedef adjacency_list<vecS, vecS, undirectedS,
property<vertex_color_t, default_color_type,
property<vertex_degree_t,int> > > Graph;
typedef graph_traits<Graph>::vertex_descriptor Vertex;
typedef std::pair<size_type, size_type> Pair;
}
-
+
namespace internal
{
// (faster than directly submitting
// indices)
ConstraintMatrix constraints;
- if (use_constraints)
+ if (use_constraints)
DoFTools::make_hanging_node_constraints (dof_handler, constraints);
constraints.close ();
CompressedSimpleSparsityPattern csp (dof_handler.n_dofs(),
}
}
-
+
template <class DH>
void
Cuthill_McKee (DH& dof_handler,
graph_degree;
internal::create_graph (dof_handler, use_constraints, graph, graph_degree);
-
+
types::property_map<types::Graph, types::vertex_index_t>::type
index_map = get(::boost::vertex_index, graph);
::boost::cuthill_mckee_ordering(graph, inv_perm.begin(),
get(::boost::vertex_color, graph),
make_degree_map(graph));
-
+
for (types::size_type c = 0; c != inv_perm.size(); ++c)
new_dof_indices[index_map[inv_perm[c]]] = c;
graph_degree;
internal::create_graph (dof_handler, use_constraints, graph, graph_degree);
-
+
types::property_map<types::Graph, types::vertex_index_t>::type
index_map = get(::boost::vertex_index, graph);
::boost::king_ordering(graph, inv_perm.rbegin());
else
::boost::king_ordering(graph, inv_perm.begin());
-
+
for (types::size_type c = 0; c != inv_perm.size(); ++c)
new_dof_indices[index_map[inv_perm[c]]] = c;
add_edge (dofs_on_this_cell[j], dofs_on_this_cell[i], G);
}
}
-
-
+
+
typedef std::vector<int> Vector;
}
} // namespace boost
-
+
template <class DH>
// empty, and calling condense with
// it is a no-op
ConstraintMatrix constraints;
- if (use_constraints)
+ if (use_constraints)
DoFTools::make_hanging_node_constraints (dof_handler, constraints);
constraints.close ();
-
+
SparsityPattern sparsity;
if (DH::dimension < 2)
{
// constraints are not needed anymore
constraints.clear ();
-
+
Assert(new_indices.size() == sparsity.n_rows(),
ExcDimensionMismatch(new_indices.size(),
sparsity.n_rows()));
const std::vector<unsigned int> &starting_indices)
{
//TODO: we should be doing the same here as in the other compute_CMK function to preserve some memory
-
+
// make the connection graph
SparsityPattern sparsity (dof_handler.n_dofs(level),
dof_handler.max_couplings_between_dofs());
MGTools::make_sparsity_pattern (dof_handler, sparsity, level);
-
+
std::vector<unsigned int> new_indices(sparsity.n_rows());
SparsityTools::reorder_Cuthill_McKee (sparsity, new_indices,
starting_indices);
template <int dim>
void
-
component_wise (DoFHandler<dim> &dof_handler,
const std::vector<unsigned int> &component_order_arg)
{
const typename DoFHandler<dim>::cell_iterator
end = dof_handler.end();
- unsigned int result =
+ const unsigned int result =
compute_component_wise<dim, ITERATOR,
typename DoFHandler<dim>::cell_iterator>(renumbering,
start, end,
component_order_arg);
if (result == 0) return;
-
+
Assert (result == dof_handler.n_dofs(),
ExcRenumberingIncomplete());
-
+
dof_handler.renumber_dofs (renumbering);
}
template <int dim>
- void
+ void
+ component_wise (hp::DoFHandler<dim> &dof_handler,
+ const std::vector<unsigned int> &component_order_arg)
+ {
+ std::vector<unsigned int> renumbering (dof_handler.n_dofs(),
+ hp::DoFHandler<dim>::invalid_dof_index);
+
+ typedef
+ internal::WrapDoFIterator<typename hp::DoFHandler<dim>::active_cell_iterator> ITERATOR;
+
+ typename hp::DoFHandler<dim>::active_cell_iterator
+ istart = dof_handler.begin_active();
+ ITERATOR start = istart;
+ const typename hp::DoFHandler<dim>::cell_iterator
+ end = dof_handler.end();
+
+ const unsigned int result =
+ compute_component_wise<dim, ITERATOR,
+ typename hp::DoFHandler<dim>::cell_iterator>(renumbering,
+ start, end,
+ component_order_arg);
+
+ if (result == 0) return;
+
+ Assert (result == dof_handler.n_dofs(),
+ ExcRenumberingIncomplete());
+
+ dof_handler.renumber_dofs (renumbering);
+ }
+
+
+ template <int dim>
+ void
component_wise (MGDoFHandler<dim> &dof_handler,
const unsigned int level,
const std::vector<unsigned int> &component_order_arg)
typedef
internal::WrapMGDoFIterator<typename MGDoFHandler<dim>::cell_iterator> ITERATOR;
-
+
typename MGDoFHandler<dim>::cell_iterator
istart =dof_handler.begin(level);
ITERATOR start = istart;
renumbering, start, end, component_order_arg);
if (result == 0) return;
-
+
Assert (result == dof_handler.n_dofs(level),
ExcRenumberingIncomplete());
}
-
template <int dim>
- void
+ void
component_wise (MGDoFHandler<dim> &dof_handler,
const std::vector<unsigned int> &component_order_arg)
{
= &component_wise<dim>;
Threads::Task<>
task = Threads::new_task (non_mg_part, dof_handler, component_order_arg);
-
+
for (unsigned int level=0; level<dof_handler.get_tria().n_levels(); ++level)
component_wise (dof_handler, level, component_order_arg);
const ENDITERATOR& end,
const std::vector<unsigned int> &component_order_arg)
{
-//TODO: Modify for hp
- const FiniteElement<dim>& fe = start->get_fe();
- const unsigned int dofs_per_cell = fe.dofs_per_cell;
-
-// Assert (new_indices.size() == start->dof_handler().n_dofs(),
-// ExcDimensionMismatch(new_indices.size(),
-// start->dof_handler.n_dofs()));
+ const hp::FECollection<dim> fe_collection (start->get_fe ());
// do nothing if the FE has only
// one component
- if (fe.n_components() == 1)
+ if (fe_collection.n_components() == 1)
{
new_indices.resize(0);
return 0;
}
+ const FiniteElement<dim>& fe = fe_collection[0];
+ const unsigned int dofs_per_cell = fe.dofs_per_cell;
+
+// Assert (new_indices.size() == start->dof_handler().n_dofs(),
+// ExcDimensionMismatch(new_indices.size(),
+// start->dof_handler.n_dofs()));
+
// Copy last argument into a
// writable vector.
std::vector<unsigned int> component_order (component_order_arg);
Assert (component_order.size() == fe.n_components(),
ExcDimensionMismatch(component_order.size(), fe.n_components()));
-
+
for (unsigned int i=0; i< component_order.size(); ++i)
Assert(component_order[i] < fe.n_components(),
ExcIndexRange(component_order[i], 0, fe.n_components()));
-
+
// vector to hold the dof indices on
// the cell we visit at a time
std::vector<unsigned int> local_dof_indices(dofs_per_cell);
// component
component_list[i] = component_order[c];
};
-
+
// set up a map where for each
// component the respective degrees
// of freedom are collected.
for (unsigned int i=0; i<dofs_per_cell; ++i)
component_to_dof_map[component_list[i]].push_back (local_dof_indices[i]);
};
-
+
// now we've got all indices sorted
// into buckets labelled with their
// target component number. we've
for (unsigned int c=0; c<fe.n_components(); ++c)
{
const unsigned int component = c;
-
+
const typename std::vector<unsigned int>::const_iterator
begin_of_component = component_to_dof_map[component].begin(),
end_of_component = component_to_dof_map[component].end();
-
+
for (typename std::vector<unsigned int>::const_iterator
dof_index = begin_of_component;
dof_index != end_of_component; ++dof_index)
// their selection state
Assert (new_indices.size() == n_dofs,
ExcDimensionMismatch(new_indices.size(), n_dofs));
-
+
const unsigned int n_selected_dofs = std::count (selected_dofs.begin(),
selected_dofs.end(),
false);
-
+
unsigned int next_unselected = 0;
unsigned int next_selected = n_selected_dofs;
for (unsigned int i=0; i<n_dofs; ++i)
std::vector<unsigned int> renumbering(dof.n_dofs());
std::vector<unsigned int> reverse(dof.n_dofs());
compute_cell_wise_dg(renumbering, reverse, dof, cells);
-
+
dof.renumber_dofs(renumbering);
}
-
+
template <class DH>
void
std::vector<unsigned int> renumbering(dof.n_dofs());
std::vector<unsigned int> reverse(dof.n_dofs());
compute_cell_wise(renumbering, reverse, dof, cells);
-
+
dof.renumber_dofs(renumbering);
}
// Later, irs inverse is computed
// into new_indices, which is the
// return argument.
-
+
Assert(new_indices.size() == n_global_dofs,
ExcDimensionMismatch(new_indices.size(), n_global_dofs));
Assert(reverse.size() == n_global_dofs,
ExcDimensionMismatch(reverse.size(), n_global_dofs));
-
+
std::vector<unsigned int> cell_dofs;
unsigned int global_index = 0;
-
+
typename std::vector<typename DH::cell_iterator>::const_iterator cell;
for(cell = cells.begin(); cell != cells.end(); ++cell)
// numbers.
unsigned int n_cell_dofs = (*cell)->get_fe().n_dofs_per_cell();
cell_dofs.resize(n_cell_dofs);
-
+
(*cell)->get_dof_indices(cell_dofs);
// Sort here to make sure that
dof.get_tria().n_active_cells()));
unsigned int n_global_dofs = dof.n_dofs();
-
+
// Actually, we compute the
// inverse of the reordering
// vector, called reverse here.
// Later, irs inverse is computed
// into new_indices, which is the
// return argument.
-
+
Assert(new_indices.size() == n_global_dofs,
ExcDimensionMismatch(new_indices.size(), n_global_dofs));
Assert(reverse.size() == n_global_dofs,
std::vector<unsigned int> cell_dofs;
unsigned int global_index = 0;
-
+
typename std::vector<typename DH::cell_iterator>::const_iterator cell;
for(cell = cells.begin(); cell != cells.end(); ++cell)
// numbers.
unsigned int n_cell_dofs = (*cell)->get_fe().n_dofs_per_cell();
cell_dofs.resize(n_cell_dofs);
-
+
(*cell)->get_dof_indices(cell_dofs);
// Sort here to make sure that
new_indices[reverse[i]] = i;
}
-
+
template <int dim>
void cell_wise_dg (
{
std::vector<unsigned int> renumbering(dof.n_dofs(level));
std::vector<unsigned int> reverse(dof.n_dofs(level));
-
+
compute_cell_wise_dg(renumbering, reverse, dof, level, cells);
dof.renumber_dofs(level, renumbering);
}
-
+
template <int dim>
void cell_wise (
{
std::vector<unsigned int> renumbering(dof.n_dofs(level));
std::vector<unsigned int> reverse(dof.n_dofs(level));
-
+
compute_cell_wise(renumbering, reverse, dof, level, cells);
dof.renumber_dofs(level, renumbering);
}
-
+
template <int dim>
void compute_cell_wise_dg (
Assert(dof.get_fe().n_dofs_per_vertex()==0,
ExcNotDGFEM());
}
-
+
Assert (new_order.size() == dof.n_dofs(level),
ExcDimensionMismatch(new_order.size(), dof.n_dofs(level)));
Assert (reverse.size() == dof.n_dofs(level),
ExcDimensionMismatch(reverse.size(), dof.n_dofs(level)));
-
+
unsigned int n_global_dofs = dof.n_dofs(level);
unsigned int n_cell_dofs = dof.get_fe().n_dofs_per_cell();
-
+
std::vector<unsigned int> cell_dofs(n_cell_dofs);
-
+
unsigned int global_index = 0;
-
+
typename std::vector<typename MGDoFHandler<dim>::cell_iterator>::const_iterator cell;
for(cell = cells.begin(); cell != cells.end(); ++cell)
Assert(global_index == n_global_dofs, ExcRenumberingIncomplete());
for (unsigned int i=0;i<new_order.size(); ++i)
- new_order[reverse[i]] = i;
+ new_order[reverse[i]] = i;
}
ExcDimensionMismatch(new_order.size(), dof.n_dofs(level)));
Assert (reverse.size() == dof.n_dofs(level),
ExcDimensionMismatch(reverse.size(), dof.n_dofs(level)));
-
+
unsigned int n_global_dofs = dof.n_dofs(level);
unsigned int n_cell_dofs = dof.get_fe().n_dofs_per_cell();
-
+
std::vector<bool> already_sorted(n_global_dofs, false);
std::vector<unsigned int> cell_dofs(n_cell_dofs);
-
+
unsigned int global_index = 0;
-
+
typename std::vector<typename MGDoFHandler<dim>::cell_iterator>::const_iterator cell;
for(cell = cells.begin(); cell != cells.end(); ++cell)
(*cell)->get_mg_dof_indices(cell_dofs);
std::sort(cell_dofs.begin(), cell_dofs.end());
-
+
for (unsigned int i=0;i<n_cell_dofs;++i)
{
if (!already_sorted[cell_dofs[i]])
Assert(global_index == n_global_dofs, ExcRenumberingIncomplete());
for (unsigned int i=0;i<new_order.size(); ++i)
- new_order[reverse[i]] = i;
+ new_order[reverse[i]] = i;
}
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();
-
+
copy (begin, end, ordered_cells.begin());
std::sort (ordered_cells.begin(), ordered_cells.end(), comparator);
compute_cell_wise_dg(new_indices, reverse, dof, ordered_cells);
}
-
+
template <class DH, int dim>
void
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();
-
+
copy (begin, end, ordered_cells.begin());
std::sort (ordered_cells.begin(), ordered_cells.end(), comparator);
-
+
compute_cell_wise_dg(new_indices, reverse, dof, ordered_cells);
}
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();
-
+
copy (begin, end, ordered_cells.begin());
std::sort (ordered_cells.begin(), ordered_cells.end(), comparator);
-
+
compute_cell_wise(new_indices, reverse, dof, ordered_cells);
}
-
+
template <int dim>
void downstream_dg (MGDoFHandler<dim>& dof,
const unsigned int level,
std::vector<unsigned int> renumbering(dof.n_dofs(level));
std::vector<unsigned int> reverse(dof.n_dofs(level));
compute_downstream_dg(renumbering, reverse, dof, level, direction);
-
+
dof.renumber_dofs(level, renumbering);
}
-
+
template <int dim>
void downstream (MGDoFHandler<dim>& dof,
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);
-
+
dof.renumber_dofs(level, renumbering);
}
-
+
template <int dim>
void
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);
-
+
std::copy (begin, end, ordered_cells.begin());
std::sort (ordered_cells.begin(), ordered_cells.end(), comparator);
-
+
compute_cell_wise_dg(new_indices, reverse, dof, level, ordered_cells);
}
-
+
template <int dim>
void
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);
-
+
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);
}
* Revert sorting order.
*/
bool counter;
-
+
/**
* Constructor.
*/
bool operator () (const DHCellIterator& c1,
const DHCellIterator& c2) const
{
-
+
const Point<dim> v1 = c1->center() - center;
const Point<dim> v2 = c2->center() - center;
const double s1 = std::atan2(v1(0), v1(1));
};
}
-
+
template <class DH, int dim>
void
std::vector<typename DH::cell_iterator>
ordered_cells(dof.get_tria().n_active_cells());
internal::ClockCells<dim> comparator(center, counter);
-
+
typename DH::active_cell_iterator begin = dof.begin_active();
typename DH::active_cell_iterator end = dof.end();
-
+
std::copy (begin, end, ordered_cells.begin());
std::sort (ordered_cells.begin(), ordered_cells.end(), comparator);
compute_cell_wise_dg(new_indices, reverse, dof, ordered_cells);
}
-
+
template <int dim>
void clockwise_dg (MGDoFHandler<dim>& dof,
std::vector<typename MGDoFHandler<dim>::cell_iterator>
ordered_cells(dof.get_tria().n_cells(level));
internal::ClockCells<dim> comparator(center, counter);
-
+
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);
dof_handler.renumber_dofs(renumbering);
}
-
+
template <class DH>
void
const unsigned int n_dofs = dof_handler.n_dofs();
Assert(new_indices.size() == n_dofs,
ExcDimensionMismatch(new_indices.size(), n_dofs));
-
+
for (unsigned i=0; i<n_dofs; ++i)
new_indices[i] = i;
-
+
std::random_shuffle (new_indices.begin(), new_indices.end());
}
dof_handler.renumber_dofs(renumbering);
}
-
+
template <class DH>
void
const unsigned int n_subdomains
= *std::max_element (subdomain_association.begin(),
subdomain_association.end()) + 1;
-
+
// then renumber the subdomains by first
// looking at those belonging to subdomain
// 0, then those of subdomain 1, etc. note
compute_minimum_degree (std::vector<unsigned int> &,
const DoFHandler<deal_II_dimension> &, bool, bool);
}
-
+
// non-boost functions:
(DoFHandler<deal_II_dimension>&,
const std::vector<unsigned int>&);
+ template
+ void component_wise<deal_II_dimension>
+ (hp::DoFHandler<deal_II_dimension>&,
+ const std::vector<unsigned int>&);
+
// template
// void
// compute_component_wise<deal_II_dimension>