-<h3>The hp::DoFHandler class, associating cells with finite elements, and constraints</h3>
+<h3>The DoFHandler class in <i>hp</i>-mode, associating cells with finite elements, and constraints</h3>
The next task we have to consider is what to do with the list of
finite element objects we want to use. In previous tutorial programs,
version of the DoFHandler class that can deal with this situation, and
(ii) a way to tell the DoF handler which element to use on which cell.
-The first of these two things is implemented in the hp::DoFHandler
-class: rather than associating it with a triangulation and a single
-finite element object, it is associated with a triangulation and a
-finite element collection. The second part is achieved by a loop over
-all cells of this hp::DoFHandler and for each cell setting the index
+The first of these two things is implemented in the <i>hp</i>-mode of
+the DoFHandler class: rather than associating it with a triangulation
+and a single finite element object, it is associated with a triangulation
+and a finite element collection. The second part is achieved by a loop
+over all cells of this DoFHandler and for each cell setting the index
of the finite element within the collection that shall be used on this
cell. We call the index of the finite element object within the
collection that shall be used on a cell the cell's <i>active FE
inactive on it. The general outline of this reads like this:
@code
- hp::DoFHandler<dim> dof_handler(triangulation);
+ DoFHandler<dim> dof_handler(triangulation);
for (auto &cell: dof_handler.active_cell_iterators())
cell->set_active_fe_index(...);
dof_handler.distribute_dofs(fe_collection);
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
-// These are the new files we need. The first and second provide <i>hp</i>
-// versions of the DoFHandler and FEValues classes as described in the
-// introduction of this program. The last one provides Fourier transformation
-// class on the unit cell.
-#include <deal.II/hp/dof_handler.h>
+// These are the new files we need. The first and second provide the
+// FECollection and the <i>hp</i> version of the FEValues class as described in
+// the introduction of this program. The last one provides Fourier
+// transformation class on the unit cell.
+#include <deal.II/hp/fe_collection.h>
#include <deal.II/hp/fe_values.h>
#include <deal.II/fe/fe_series.h>
// computes this estimated smoothness, as discussed in the introduction.
//
// As far as member variables are concerned, we use the same structure as
- // already used in step-6, but instead of a regular DoFHandler we use an
- // object of type hp::DoFHandler, and we need collections instead of
+ // already used in step-6, but we need collections instead of
// individual finite element, quadrature, and face quadrature objects. We
// will fill these collections in the constructor of the class. The last
// variable, <code>max_degree</code>, indicates the maximal polynomial
// @sect4{LaplaceProblem::LaplaceProblem constructor}
// The constructor of this class is fairly straightforward. It associates
- // the hp::DoFHandler object with the triangulation, and then sets the
+ // the DoFHandler object with the triangulation, and then sets the
// maximal polynomial degree to 7 (in 1d and 2d) or 5 (in 3d and higher). We
// do so because using higher order polynomial degrees becomes prohibitively
// expensive, especially in higher space dimensions.
template <int dim>
LaplaceProblem<dim>::LaplaceProblem()
- : dof_handler(triangulation, /*enable_hp_capability */ true)
+ : dof_handler(triangulation)
, max_degree(dim <= 2 ? 7 : 5)
{
for (unsigned int degree = 2; degree <= max_degree; ++degree)
// With now all data vectors available -- solution, estimated errors and
// smoothness indicators, and finite element degrees --, we create a
- // DataOut object for graphical output and attach all data. Note that
- // the DataOut class has a second template argument (which defaults to
- // DoFHandler@<dim@>, which is why we have never seen it in previous
- // tutorial programs) that indicates the type of DoF handler to be
- // used. Here, we have to use the hp::DoFHandler class:
- DataOut<dim, DoFHandler<dim>> data_out;
+ // DataOut object for graphical output and attach all data:
+ DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution, "solution");
and facilities of the hp namespace to assign different elements to different
cells. In other words, we will use collect the two finite elements in an
hp::FECollection, will integrate with an appropriate hp::QCollection using an
-hp::FEValues object, and our DoF handler will be of type hp::DoFHandler. You
+hp::FEValues object, and our DoFHandler will be in <i>hp</i>-mode. You
may wish to take a look at step-27 for an overview of all of these concepts.
Before going on describing the testcase, let us clarify a bit <i>why</i> this
material_id numbers and will use them to identify which part of the domain a
cell is on.
-Secondly, we use an object of type hp::DoFHandler. This class needs to know
-which cells will use the Stokes and which the elasticity finite element. At
-the beginning of each refinement cycle we will therefore have to walk over
-all cells and set the (in hp parlance) active FE index to whatever is
-appropriate in the current situation. While we can use symbolic names for the
-material id, the active FE index is in fact a number that will frequently be
-used to index into collections of objects (e.g. of type hp::FECollection and
-hp::QCollection); that means that the active FE index actually has to have
-value zero for the fluid and one for the elastic part of the domain.
+Secondly, we use an object of type DoFHandler operating in <i>hp</i>-mode. This
+class needs to know which cells will use the Stokes and which the elasticity
+finite element. At the beginning of each refinement cycle we will therefore
+have to walk over all cells and set the (in hp parlance) active FE index to
+whatever is appropriate in the current situation. While we can use symbolic
+names for the material id, the active FE index is in fact a number that will
+frequently be used to index into collections of objects (e.g. of type
+hp::FECollection and hp::QCollection); that means that the active FE index
+actually has to have value zero for the fluid and one for the elastic part of
+the domain.
<h4>Linear solvers</h4>
is likely to be very large at the locations where the solution is
discontinuous and extended by zero; it also doesn't become smaller as the mesh
is refined. The KellyErrorEstimator class can't just ignore the interface
-because it essentially only sees an hp::DoFHandler where the element type
-changes from one cell to another — precisely the thing that the
-hp::DoFHandler was designed for, the interface in the current program looks no
+because it essentially only sees a DoFHandler in <i>hp</i>-mode where the element
+type changes from one cell to another — precisely the thing that the
+<i>hp</i>-mode was designed for, the interface in the current program looks no
different than the interfaces in step-27, for example, and certainly no less
legitimate. Be that as it may, the end results is that there is a layer of
cells on both sides of the interface between the two subdomains where error
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
-#include <deal.II/hp/dof_handler.h>
#include <deal.II/hp/fe_collection.h>
#include <deal.II/hp/fe_values.h>
// This is the main class. It is, if you want, a combination of step-8 and
// step-22 in that it has member variables that either address the global
- // problem (the Triangulation and hp::DoFHandler objects, as well as the
+ // problem (the Triangulation and DoFHandler objects, as well as the
// hp::FECollection and various linear algebra objects) or that pertain to
// either the elasticity or Stokes sub-problems. The general structure of
// the class, however, is like that of most of the other programs
1,
FE_Q<dim>(elasticity_degree),
dim)
- , dof_handler(triangulation, /*enable_hp_capability */ true)
+ , dof_handler(triangulation)
, viscosity(2)
, lambda(1)
, mu(1)
// Generating graphical output is rather trivial here: all we have to do is
// identify which components of the solution vector belong to scalars and/or
// vectors (see, for example, step-22 for a previous example), and then pass
- // it all on to the DataOut class (with the second template argument equal
- // to hp::DoFHandler instead of the usual default DoFHandler):
+ // it all on to the DataOut class:
template <int dim>
void FluidStructureProblem<dim>::output_results(
const unsigned int refinement_cycle) const
data_component_interpretation.push_back(
DataComponentInterpretation::component_is_part_of_vector);
- DataOut<dim, DoFHandler<dim>> data_out;
+ DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution,
solution_names,
- DataOut<dim, DoFHandler<dim>>::type_dof_data,
+ DataOut<dim>::type_dof_data,
data_component_interpretation);
data_out.build_patches();
*/
static const unsigned int space_dimension = spacedim;
- /**
- * Boolean indicating whether or not the current DoFHander has hp
- * capabilities.
- */
- const bool hp_capability_enabled;
-
/**
* The default index of the finite element to be used on a given cell.
*/
* object with this constructor, use initialize() to make a valid
* DoFHandler.
*/
- DoFHandler(const bool hp_capability_enabled = false);
+ DoFHandler();
/**
* Constructor. Take @p tria as the triangulation to work on.
*/
- DoFHandler(const Triangulation<dim, spacedim> &tria,
- const bool hp_capability_enabled = false);
+ DoFHandler(const Triangulation<dim, spacedim> &tria);
/**
* Copy constructor. DoFHandler objects are large and expensive.
void
distribute_mg_dofs();
+ /**
+ * Returns whether this DoFHandler has hp capabilities.
+ */
+ bool
+ has_hp_capabilities() const;
+
/**
* This function returns whether this DoFHandler has DoFs distributed on
* each multigrid level or in other words if distribute_mg_dofs() has been
*/
BlockInfo block_info_object;
+ /**
+ * Boolean indicating whether or not the current DoFHandler has hp
+ * capabilities.
+ */
+ bool hp_capability_enabled;
+
/**
* Address of the triangulation to work on.
*/
setup_policy();
/**
- * Setup DoFHandler policy and listeners (in the hp-context).
+ * Setup connections to refinement signals of the underlying triangulation.
+ * Necessary for the hp-mode.
*/
void
- setup_policy_and_listeners();
+ connect_to_triangulation_signals();
/**
* Create default tables for the active_fe_indices.
*/
+template <int dim, int spacedim>
+inline bool
+DoFHandler<dim, spacedim>::has_hp_capabilities() const
+{
+ return hp_capability_enabled;
+}
+
+
+
template <int dim, int spacedim>
inline bool
DoFHandler<dim, spacedim>::has_level_dofs() const
dof_handler[0]->get_triangulation(),
cell_level_index,
face_info,
- dof_handler[0]->hp_capability_enabled ?
+ dof_handler[0]->has_hp_capabilities() ?
dof_info[0].cell_active_fe_index :
std::vector<unsigned int>(),
mapping,
std::vector<Tensor<1, spacedim, typename VectorType::value_type>>
&gradients)
{
- if (dof.hp_capability_enabled == false)
+ if (dof.has_hp_capabilities() == false)
point_gradient(StaticMappingQ1<dim, spacedim>::mapping,
dof,
fe_function,
const VectorType & fe_function,
const Point<spacedim> & point)
{
- if (dof.hp_capability_enabled == false)
+ if (dof.has_hp_capabilities() == false)
return point_gradient(StaticMappingQ1<dim, spacedim>::mapping,
dof,
fe_function,
const Point<spacedim> & point,
Vector<typename VectorType::value_type> &value)
{
- if (dof.hp_capability_enabled == false)
+ if (dof.has_hp_capabilities() == false)
point_value(StaticMappingQ1<dim, spacedim>::mapping,
dof,
fe_function,
const VectorType & fe_function,
const Point<spacedim> & point)
{
- if (dof.hp_capability_enabled == false)
+ if (dof.has_hp_capabilities() == false)
return point_value(StaticMappingQ1<dim, spacedim>::mapping,
dof,
fe_function,
const Point<spacedim> & p,
Vector<double> & rhs_vector)
{
- if (dof_handler.hp_capability_enabled)
+ if (dof_handler.has_hp_capabilities())
create_point_source_vector(hp::StaticMappingQ1<dim>::mapping_collection,
dof_handler,
p,
const Point<dim> & orientation,
Vector<double> & rhs_vector)
{
- if (dof_handler.hp_capability_enabled)
+ if (dof_handler.has_hp_capabilities())
create_point_source_vector(hp::StaticMappingQ1<dim>::mapping_collection,
dof_handler,
p,
cell_iterator &cell_,
const typename Triangulation<dim, DoFHandlerType::space_dimension>::
CellStatus status) { return this->pack_callback(cell_, status); },
- /*returns_variable_size_data=*/dof_handler->hp_capability_enabled);
+ /*returns_variable_size_data=*/dof_handler->has_hp_capabilities());
}
input_vectors.size());
unsigned int fe_index = 0;
- if (dof_handler->hp_capability_enabled)
+ if (dof_handler->has_hp_capabilities())
{
switch (status)
{
typename DoFHandlerType::cell_iterator cell(*cell_, dof_handler);
unsigned int fe_index = 0;
- if (dof_handler->hp_capability_enabled)
+ if (dof_handler->has_hp_capabilities())
{
switch (status)
{
bool levels_only,
bool active_only)
{
- Assert(dof.hp_capability_enabled == false,
+ Assert(dof.has_hp_capabilities() == false,
(typename DoFHandler<dim, spacedim>::ExcNotImplementedWithHP()));
if (!levels_only && dof.has_active_dofs())
void
BlockInfo::initialize_local(const DoFHandler<dim, spacedim> &dof)
{
- Assert(dof.hp_capability_enabled == false,
+ Assert(dof.has_hp_capabilities() == false,
(typename DoFHandler<dim, spacedim>::ExcNotImplementedWithHP()));
const FiniteElement<dim, spacedim> & fe = dof.get_fe();
template <int dim, int spacedim>
-DoFHandler<dim, spacedim>::DoFHandler(const bool hp_capability_enabled)
- : hp_capability_enabled(hp_capability_enabled)
+DoFHandler<dim, spacedim>::DoFHandler()
+ : hp_capability_enabled(true)
, tria(nullptr, typeid(*this).name())
, mg_faces(nullptr)
{}
template <int dim, int spacedim>
-DoFHandler<dim, spacedim>::DoFHandler(const Triangulation<dim, spacedim> &tria,
- const bool hp_capability_enabled)
- : hp_capability_enabled(hp_capability_enabled)
+DoFHandler<dim, spacedim>::DoFHandler(const Triangulation<dim, spacedim> &tria)
+ : hp_capability_enabled(true)
, tria(&tria, typeid(*this).name())
, mg_faces(nullptr)
{
- if (hp_capability_enabled)
- {
- this->setup_policy_and_listeners();
- this->create_active_fe_table();
- }
- else
- {
- this->setup_policy();
- }
+ this->setup_policy();
+
+ // provide all hp functionalities before finite elements are registered
+ this->connect_to_triangulation_signals();
+ this->create_active_fe_table();
}
+
+
template <int dim, int spacedim>
DoFHandler<dim, spacedim>::~DoFHandler()
{
- if (hp_capability_enabled)
- {
- // unsubscribe as a listener to refinement of the underlying
- // triangulation
- for (auto &connection : this->tria_listeners)
- connection.disconnect();
- this->tria_listeners.clear();
-
- // ...and release allocated memory
- // virtual functions called in constructors and destructors never use the
- // override in a derived class
- // for clarity be explicit on which function is called
- DoFHandler<dim, spacedim>::clear();
- }
- else
- {
- // release allocated memory
- // virtual functions called in constructors and destructors never use the
- // override in a derived class
- // for clarity be explicit on which function is called
- DoFHandler<dim, spacedim>::clear();
-
- // also release the policy. this needs to happen before the
- // current object disappears because the policy objects
- // store references to the DoFhandler object they work on
- this->policy.reset();
- }
+ // unsubscribe all attachments to refinement signals of the underlying
+ // triangulation
+ for (auto &connection : this->tria_listeners)
+ connection.disconnect();
+ this->tria_listeners.clear();
+
+ // release allocated memory
+ // virtual functions called in constructors and destructors never use the
+ // override in a derived class
+ // for clarity be explicit on which function is called
+ DoFHandler<dim, spacedim>::clear();
+
+ // also release the policy. this needs to happen before the
+ // current object disappears because the policy objects
+ // store references to the DoFhandler object they work on
+ this->policy.reset();
}
if (this->tria != &tria)
{
+ // remove association with old triangulation
for (auto &connection : this->tria_listeners)
connection.disconnect();
this->tria_listeners.clear();
-
- this->tria = &tria;
-
- this->setup_policy_and_listeners();
}
-
- this->create_active_fe_table();
-
- this->distribute_dofs(fe);
}
else
{
- this->tria = &tria;
// this->faces = nullptr;
this->number_cache.n_global_dofs = 0;
+ }
+ if (this->tria != &tria)
+ {
+ // establish connection to new triangulation
+ this->tria = &tria;
this->setup_policy();
- this->distribute_dofs(fe);
+ // start in hp-mode and let distribute_dofs toggle it if necessary
+ hp_capability_enabled = true;
+ this->connect_to_triangulation_signals();
+ this->create_active_fe_table();
}
+
+ this->distribute_dofs(fe);
}
// don't create a new object if the one we have is already appropriate
if (this->fe_collection != ff)
- this->fe_collection = hp::FECollection<dim, spacedim>(ff);
+ {
+ this->fe_collection = hp::FECollection<dim, spacedim>(ff);
+
+ const bool contains_multiple_fes = (this->fe_collection.size() > 1);
+
+ // disable hp-mode if only a single finite element has been registered
+ if (hp_capability_enabled && !contains_multiple_fes)
+ {
+ hp_capability_enabled = false;
+
+ // unsubscribe connections to signals
+ for (auto &connection : this->tria_listeners)
+ connection.disconnect();
+ this->tria_listeners.clear();
+
+ // release active and future finite element tables
+ this->hp_cell_active_fe_indices.clear();
+ this->hp_cell_active_fe_indices.shrink_to_fit();
+ this->hp_cell_future_fe_indices.clear();
+ this->hp_cell_future_fe_indices.shrink_to_fit();
+ }
+
+ // re-enabling hp-mode is not permitted since the active and future fe
+ // tables are no longer available
+ AssertThrow(
+ hp_capability_enabled || !contains_multiple_fes,
+ ExcMessage(
+ "You cannot re-enable hp capabilities after you registered a single "
+ "finite element. Please create a new DoFHandler object instead."));
+ }
if (hp_capability_enabled)
{
- // ensure that the active_fe_indices vectors are initialized correctly
- this->create_active_fe_table();
-
// make sure every processor knows the active_fe_indices
// on both its own cells and all ghost cells
dealii::internal::hp::DoFHandlerImplementation::Implementation::
DoFHandler<dim, spacedim>::distribute_dofs(
const hp::FECollection<dim, spacedim> &ff)
{
+ this->set_fe(ff);
+
if (hp_capability_enabled)
{
object_dof_indices.resize(this->tria->n_levels());
object_dof_ptr.resize(this->tria->n_levels());
cell_dof_cache_indices.resize(this->tria->n_levels());
cell_dof_cache_ptr.resize(this->tria->n_levels());
- hp_cell_active_fe_indices.resize(this->tria->n_levels());
- hp_cell_future_fe_indices.resize(this->tria->n_levels());
- // assign the fe_collection and initialize all active_fe_indices
- this->set_fe(ff);
// If an underlying shared::Tria allows artificial cells,
// then save the current set of subdomain ids, and set
}
else
{
- // first, assign the finite_element
- this->set_fe(ff);
-
// delete all levels and set them up newly. note that we still have to
// allocate space for all degrees of freedom on this mesh (including ghost
// and cells that are entirely stored on different processors), though we
if (hp_capability_enabled)
{
this->hp_cell_active_fe_indices.clear();
+ this->hp_cell_active_fe_indices.shrink_to_fit();
this->hp_cell_future_fe_indices.clear();
-
- object_dof_indices.clear();
+ this->hp_cell_future_fe_indices.shrink_to_fit();
}
else
{
template <int dim, int spacedim>
void
-DoFHandler<dim, spacedim>::setup_policy_and_listeners()
+DoFHandler<dim, spacedim>::connect_to_triangulation_signals()
{
// connect functions to signals of the underlying triangulation
this->tria_listeners.push_back(this->tria->signals.pre_refinement.connect(
const dealii::parallel::DistributedTriangulationBase<dim, spacedim> *>(
&this->get_triangulation()))
{
- this->policy =
- std::make_unique<internal::DoFHandlerImplementation::Policy::
- ParallelDistributed<dim, spacedim>>(*this);
-
// repartitioning signals
this->tria_listeners.push_back(
this->tria->signals.pre_distributed_repartition.connect([this]() {
const dealii::parallel::shared::Triangulation<dim, spacedim> *>(
&this->get_triangulation()) != nullptr)
{
- this->policy =
- std::make_unique<internal::DoFHandlerImplementation::Policy::
- ParallelShared<dim, spacedim>>(*this);
-
// partitioning signals
this->tria_listeners.push_back(
this->tria->signals.pre_partition.connect([this]() {
}
else
{
- this->policy = std::make_unique<
- internal::DoFHandlerImplementation::Policy::Sequential<dim, spacedim>>(
- *this);
-
// refinement signals
this->tria_listeners.push_back(this->tria->signals.pre_refinement.connect(
[this] { this->pre_active_fe_index_transfer(); }));
const Table<2, Coupling> & table,
std::vector<Table<2, Coupling>> &tables_by_block)
{
- if (dof_handler.hp_capability_enabled == false)
+ if (dof_handler.has_hp_capabilities() == false)
{
const FiniteElement<dim, spacedim> &fe = dof_handler.get_fe();
const unsigned int nb = fe.n_blocks();
unsigned int
n_finite_elements(const DoFHandler<dim, spacedim> &dof_handler)
{
- if (dof_handler.hp_capability_enabled == true)
+ if (dof_handler.has_hp_capabilities() == true)
return dof_handler.get_fe_collection().size();
else
return 1;
std::set<unsigned int> fe_ind_face_subface;
fe_ind_face_subface.insert(cell->active_fe_index());
- if (dof_handler.hp_capability_enabled)
+ if (dof_handler.has_hp_capabilities())
for (unsigned int c = 0;
c < cell->face(face)->number_of_children();
++c)
//
// since this is something that can only happen for hp
// dof handlers, add a check here...
- Assert(dof_handler.hp_capability_enabled == true,
+ Assert(dof_handler.has_hp_capabilities() == true,
ExcInternalError());
const dealii::hp::FECollection<dim, spacedim>
// Only if there is a neighbor with a different active_fe_index
// and the same h-level, some action has to be taken.
- if ((dof_handler.hp_capability_enabled) &&
+ if ((dof_handler.has_hp_capabilities()) &&
!cell->face(face)->at_boundary() &&
(cell->neighbor(face)->active_fe_index() !=
cell->active_fe_index()) &&
bool(const typename DoFHandler<dim, spacedim>::active_cell_iterator &,
const unsigned int)> &face_has_flux_coupling)
{
- if (dof.hp_capability_enabled == false)
+ if (dof.has_hp_capabilities() == false)
{
const FiniteElement<dim, spacedim> &fe = dof.get_fe();
{
template <int dim, int spacedim>
DoFHandler<dim, spacedim>::DoFHandler()
- : dealii::DoFHandler<dim, spacedim>(true)
+ : dealii::DoFHandler<dim, spacedim>()
{}
template <int dim, int spacedim>
DoFHandler<dim, spacedim>::DoFHandler(
const Triangulation<dim, spacedim> &tria)
- : dealii::DoFHandler<dim, spacedim>(tria, true)
+ : dealii::DoFHandler<dim, spacedim>(tria)
{}
} // namespace hp
extract_interpolation_matrices(const dealii::DoFHandler<dim, spacedim> &dof,
dealii::Table<2, FullMatrix<double>> &matrices)
{
- if (dof.hp_capability_enabled == false)
+ if (dof.has_hp_capabilities() == false)
return;
const dealii::hp::FECollection<dim, spacedim> &fe = dof.get_fe_collection();
fes.push_back(FE_Q<dim>(1));
fes.push_back(FE_Q<dim>(2));
- DoFHandler<dim> dofh(tria, /*hp_capability_enabled=*/true);
+ DoFHandler<dim> dofh(tria);
dofh.set_fe(fes);
dofh.begin_active()->set_active_fe_index(1);
// create a DoFHandler on which we have both cells with FE_Q as well as
// FE_Nothing
- DoFHandler<dim> dof_handler(tria, true);
+ DoFHandler<dim> dof_handler(tria);
dof_handler.begin(0)->child(0)->set_active_fe_index(1);
IndexSet locally_relevant_dofs;
Triangulation<2> tria;
- MGLevelObject<DoFHandler<2>> dof_handlers(0, 3, tria, true);
+ MGLevelObject<DoFHandler<2>> dof_handlers(0, 3, tria);
- dof_handlers.resize(0, 5, tria, true);
+ dof_handlers.resize(0, 5, tria);
deallog << "OK!" << std::endl;
}
projects / dealii.git / commitdiff