From e759ceee8edfd88f9961d699ae154bc4c136947b Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sat, 28 Mar 2020 13:52:57 +0100 Subject: [PATCH] Merge DoFHandler and hp::DoFHandler --- cmake/config/template-arguments.in | 11 +- .../changes/incompatibilities/20200523Munch | 4 + doc/news/changes/major/20200523Munch | 7 + examples/step-27/step-27.cc | 6 +- examples/step-46/step-46.cc | 16 +- include/deal.II/distributed/cell_weights.h | 28 +- include/deal.II/dofs/block_info.h | 25 - include/deal.II/dofs/dof_accessor.h | 39 +- include/deal.II/dofs/dof_accessor.templates.h | 1998 +++------- include/deal.II/dofs/dof_handler.h | 737 +++- include/deal.II/dofs/dof_objects.h | 3 +- include/deal.II/dofs/dof_renumbering.h | 36 +- include/deal.II/dofs/dof_tools.h | 112 +- include/deal.II/fe/fe_enriched.h | 8 +- include/deal.II/fe/fe_tools.h | 14 +- .../fe/fe_tools_interpolate.templates.h | 16 +- include/deal.II/grid/grid_tools.h | 9 +- include/deal.II/hp/dof_faces.h | 654 ---- include/deal.II/hp/dof_handler.h | 1517 +------- include/deal.II/hp/dof_level.h | 610 --- include/deal.II/hp/refinement.h | 60 +- include/deal.II/matrix_free/matrix_free.h | 19 +- .../matrix_free/matrix_free.templates.h | 253 +- .../numerics/matrix_creator.templates.h | 89 +- include/deal.II/numerics/matrix_tools.h | 40 +- .../deal.II/numerics/smoothness_estimator.h | 41 +- .../deal.II/numerics/vector_tools_boundary.h | 24 +- .../vector_tools_boundary.templates.h | 15 +- .../vector_tools_integrate_difference.h | 22 +- ...tor_tools_integrate_difference.templates.h | 24 +- .../numerics/vector_tools_interpolate.h | 14 +- .../vector_tools_interpolate.templates.h | 16 +- .../numerics/vector_tools_point_gradient.h | 16 +- .../vector_tools_point_gradient.templates.h | 72 +- .../numerics/vector_tools_point_value.h | 30 +- .../vector_tools_point_value.templates.h | 135 +- .../deal.II/numerics/vector_tools_project.h | 4 +- .../numerics/vector_tools_project.templates.h | 5 +- include/deal.II/numerics/vector_tools_rhs.h | 8 +- .../numerics/vector_tools_rhs.templates.h | 12 +- source/distributed/cell_weights.cc | 43 +- source/distributed/solution_transfer.cc | 6 +- source/distributed/solution_transfer.inst.in | 42 +- source/dofs/block_info.cc | 30 +- source/dofs/block_info.inst.in | 14 - source/dofs/dof_accessor.cc | 4 +- source/dofs/dof_accessor.inst.in | 82 - source/dofs/dof_accessor_get.cc | 20 +- source/dofs/dof_accessor_get.inst.in | 28 - source/dofs/dof_accessor_set.cc | 20 +- source/dofs/dof_accessor_set.inst.in | 28 - source/dofs/dof_handler.cc | 3367 +++++++++++++---- source/dofs/dof_handler_policy.cc | 500 ++- source/dofs/dof_handler_policy.inst.in | 6 - source/dofs/dof_renumbering.cc | 81 +- source/dofs/dof_renumbering.inst.in | 128 +- source/dofs/dof_tools.cc | 79 +- source/dofs/dof_tools.inst.in | 208 +- source/dofs/dof_tools_constraints.cc | 156 +- source/dofs/dof_tools_constraints.inst.in | 14 +- source/dofs/dof_tools_sparsity.cc | 852 +++-- source/dofs/dof_tools_sparsity.inst.in | 220 +- source/fe/fe_enriched.cc | 8 +- source/fe/fe_enriched.inst.in | 7 +- source/fe/fe_tools_interpolate.inst.in | 32 - source/fe/fe_values.inst.in | 12 +- source/grid/grid_tools_dof_handlers.cc | 8 +- source/grid/grid_tools_dof_handlers.inst.in | 21 +- source/hp/CMakeLists.txt | 2 - source/hp/dof_faces.cc | 49 - source/hp/dof_handler.cc | 2337 +----------- source/hp/dof_handler.inst.in | 23 - source/hp/dof_level.cc | 301 -- source/hp/fe_values.inst.in | 2 +- source/hp/refinement.cc | 55 +- source/hp/refinement.inst.in | 30 +- source/multigrid/mg_tools.inst.in | 14 +- source/numerics/data_out_faces.inst.in | 2 - source/numerics/data_out_stack.inst.in | 25 - source/numerics/derivative_approximation.cc | 36 +- source/numerics/fe_field_function.inst.in | 7 +- source/numerics/matrix_creator.inst.in | 44 +- source/numerics/smoothness_estimator.cc | 56 +- source/numerics/smoothness_estimator.inst.in | 8 +- source/numerics/solution_transfer.cc | 14 +- source/numerics/vector_tools_boundary.inst.in | 45 +- .../numerics/vector_tools_constraints.inst.in | 38 +- .../vector_tools_integrate_difference.inst.in | 48 +- .../numerics/vector_tools_interpolate.inst.in | 49 +- .../vector_tools_point_gradient.inst.in | 27 - .../numerics/vector_tools_point_value.inst.in | 30 +- .../numerics/vector_tools_project_hp.inst.in | 38 +- source/numerics/vector_tools_rhs.cc | 25 - source/numerics/vector_tools_rhs.inst.in | 8 +- tests/fe/crash_01.cc | 7 +- tests/fe/curl_curl_01.cc | 7 +- tests/fe/fe_nedelec_sz_non_rect_2d.cc | 7 +- tests/fe/fe_nedelec_sz_non_rect_face.cc | 7 +- tests/fe/nedelec_non_rect_2d.cc | 7 +- tests/fe/nedelec_non_rect_face.cc | 7 +- ...of_values_by_interpolation_01.debug.output | 12 +- tests/mpi/project_bv_div_conf.cc | 7 +- tests/numerics/project_bv_curl_conf.cc | 7 +- tests/numerics/project_bv_curl_conf_02.cc | 7 +- .../numerics/project_bv_curl_conf_complex.cc | 7 +- tests/numerics/project_bv_div_conf.cc | 7 +- tests/numerics/project_bv_div_conf_02.cc | 6 +- ...f_handler_01.with_64bit_indices=off.output | 28 +- ...of_handler_01.with_64bit_indices=on.output | 28 +- ...f_handler_01.with_64bit_indices=off.output | 17 +- ...of_handler_01.with_64bit_indices=on.output | 17 +- 111 files changed, 5629 insertions(+), 10634 deletions(-) create mode 100644 doc/news/changes/incompatibilities/20200523Munch create mode 100644 doc/news/changes/major/20200523Munch delete mode 100644 include/deal.II/hp/dof_faces.h delete mode 100644 include/deal.II/hp/dof_level.h delete mode 100644 source/hp/dof_faces.cc delete mode 100644 source/hp/dof_level.cc diff --git a/cmake/config/template-arguments.in b/cmake/config/template-arguments.in index efe2f5111e..15007e7f08 100644 --- a/cmake/config/template-arguments.in +++ b/cmake/config/template-arguments.in @@ -208,9 +208,11 @@ GENERAL_CONTAINER_TYPES := { std::vector; ArrayView } +// DoFHandler (for iterators!) +DOFHANDLER_TEMPLATE := { DoFHandler } + // DoFHandler and hp::DoFHandler templates -DOFHANDLER_TEMPLATES := { DoFHandler; - hp::DoFHandler } +DOFHANDLER_TEMPLATES := { DoFHandler; hp::DoFHandler } // Triangulation and DoFHandler templates TRIANGULATION_AND_DOFHANDLER_TEMPLATES := { Triangulation; @@ -219,6 +221,11 @@ TRIANGULATION_AND_DOFHANDLER_TEMPLATES := { Triangulation; DoFHandler; hp::DoFHandler } +// concrete sequential Triangulation and DoFHandler with hard-coded +// (for iterators!) +SEQUENTIAL_TRIANGULATION_AND_DOFHANDLER := { Triangulation; + DoFHandler } + // concrete sequential Triangulation and DoFHandler types with hard-coded // SEQUENTIAL_TRIANGULATION_AND_DOFHANDLERS := { Triangulation; diff --git a/doc/news/changes/incompatibilities/20200523Munch b/doc/news/changes/incompatibilities/20200523Munch new file mode 100644 index 0000000000..c7be4b598d --- /dev/null +++ b/doc/news/changes/incompatibilities/20200523Munch @@ -0,0 +1,4 @@ +Deprecated: The class hp::DoFHandler has been deprecated, since the DoFHandler +has been extended with its functionalities. +
+(Peter Munch, 2020/05/23) diff --git a/doc/news/changes/major/20200523Munch b/doc/news/changes/major/20200523Munch new file mode 100644 index 0000000000..d75d8ed1c6 --- /dev/null +++ b/doc/news/changes/major/20200523Munch @@ -0,0 +1,7 @@ +Changed: The internal data structures of DoFHandler have been modified to use +compressed row storage, enabling it to also handle hp::DoFHandler functionalities. +Currently, the user can choose between the normal mode and the hp mode during +calling the constructur. Please note that the multigrid functionalities are only +available during normal mode. +
+(Peter Munch, 2020/05/23) diff --git a/examples/step-27/step-27.cc b/examples/step-27/step-27.cc index cd7d451044..0a2d9cf66b 100644 --- a/examples/step-27/step-27.cc +++ b/examples/step-27/step-27.cc @@ -105,7 +105,7 @@ namespace Step27 Triangulation triangulation; - hp::DoFHandler dof_handler; + DoFHandler dof_handler; hp::FECollection fe_collection; hp::QCollection quadrature_collection; hp::QCollection face_quadrature_collection; @@ -186,7 +186,7 @@ namespace Step27 template LaplaceProblem::LaplaceProblem() - : dof_handler(triangulation) + : dof_handler(triangulation, /*enable_hp_capability */ true) , max_degree(dim <= 2 ? 7 : 5) { for (unsigned int degree = 2; degree <= max_degree; ++degree) @@ -447,7 +447,7 @@ namespace Step27 // DoFHandler@, 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> data_out; + DataOut> data_out; data_out.attach_dof_handler(dof_handler); data_out.add_data_vector(solution, "solution"); diff --git a/examples/step-46/step-46.cc b/examples/step-46/step-46.cc index 46bef2efda..165324c736 100644 --- a/examples/step-46/step-46.cc +++ b/examples/step-46/step-46.cc @@ -104,10 +104,10 @@ namespace Step46 }; static bool cell_is_in_fluid_domain( - const typename hp::DoFHandler::cell_iterator &cell); + const typename DoFHandler::cell_iterator &cell); static bool cell_is_in_solid_domain( - const typename hp::DoFHandler::cell_iterator &cell); + const typename DoFHandler::cell_iterator &cell); void make_grid(); @@ -132,7 +132,7 @@ namespace Step46 FESystem stokes_fe; FESystem elasticity_fe; hp::FECollection fe_collection; - hp::DoFHandler dof_handler; + DoFHandler dof_handler; AffineConstraints constraints; @@ -236,7 +236,7 @@ namespace Step46 1, FE_Q(elasticity_degree), dim) - , dof_handler(triangulation) + , dof_handler(triangulation, /*enable_hp_capability */ true) , viscosity(2) , lambda(1) , mu(1) @@ -249,7 +249,7 @@ namespace Step46 template bool FluidStructureProblem::cell_is_in_fluid_domain( - const typename hp::DoFHandler::cell_iterator &cell) + const typename DoFHandler::cell_iterator &cell) { return (cell->material_id() == fluid_domain_id); } @@ -257,7 +257,7 @@ namespace Step46 template bool FluidStructureProblem::cell_is_in_solid_domain( - const typename hp::DoFHandler::cell_iterator &cell) + const typename DoFHandler::cell_iterator &cell) { return (cell->material_id() == solid_domain_id); } @@ -859,12 +859,12 @@ namespace Step46 data_component_interpretation.push_back( DataComponentInterpretation::component_is_part_of_vector); - DataOut> data_out; + DataOut> data_out; data_out.attach_dof_handler(dof_handler); data_out.add_data_vector(solution, solution_names, - DataOut>::type_dof_data, + DataOut>::type_dof_data, data_component_interpretation); data_out.build_patches(); diff --git a/include/deal.II/distributed/cell_weights.h b/include/deal.II/distributed/cell_weights.h index e6556d86f9..6066e30995 100644 --- a/include/deal.II/distributed/cell_weights.h +++ b/include/deal.II/distributed/cell_weights.h @@ -96,9 +96,9 @@ namespace parallel * weight or, in other words, the additional computational load associated * with it. */ - using WeightingFunction = std::function::cell_iterator &, - const FiniteElement &)>; + using WeightingFunction = std::function< + unsigned int(const typename DoFHandler::cell_iterator &, + const FiniteElement &)>; /** * Constructor. @@ -108,8 +108,8 @@ namespace parallel * @param[in] weighting_function The function that determines each * cell's weight during load balancing. */ - CellWeights(const hp::DoFHandler &dof_handler, - const WeightingFunction & weighting_function); + CellWeights(const dealii::DoFHandler &dof_handler, + const WeightingFunction & weighting_function); /** * Destructor. @@ -125,8 +125,8 @@ namespace parallel * Disconnects the function previously connected to the weighting signal. */ void - reinit(const hp::DoFHandler &dof_handler, - const WeightingFunction & weighting_function); + reinit(const DoFHandler &dof_handler, + const WeightingFunction & weighting_function); /** * Converts a @p weighting_function to a different type that qualifies as @@ -139,7 +139,7 @@ namespace parallel static std::function::cell_iterator &cell, const typename Triangulation::CellStatus status)> - make_weighting_callback(const hp::DoFHandler &dof_handler, + make_weighting_callback(const DoFHandler &dof_handler, const WeightingFunction &weighting_function); /** @@ -193,7 +193,7 @@ namespace parallel * determine each cell's finite element. */ DEAL_II_DEPRECATED - CellWeights(const hp::DoFHandler &dof_handler); + CellWeights(const dealii::DoFHandler &dof_handler); /** * Choose a constant weight @p factor on each cell. @@ -235,9 +235,9 @@ namespace parallel */ DEAL_II_DEPRECATED void register_custom_weighting( - const std::function &, - const typename hp::DoFHandler::cell_iterator &)> + const std::function< + unsigned int(const FiniteElement &, + const typename DoFHandler::cell_iterator &)> custom_function); /** @@ -254,7 +254,7 @@ namespace parallel * Pointer to the degree of freedom handler. */ DEAL_II_DEPRECATED - SmartPointer, CellWeights> + SmartPointer, CellWeights> dof_handler; /** @@ -288,7 +288,7 @@ namespace parallel weighting_callback( const typename Triangulation::cell_iterator &cell, const typename Triangulation::CellStatus status, - const hp::DoFHandler & dof_handler, + const DoFHandler & dof_handler, const parallel::TriangulationBase & triangulation, const WeightingFunction &weighting_function); }; diff --git a/include/deal.II/dofs/block_info.h b/include/deal.II/dofs/block_info.h index 1d93b36ff1..f1d134918f 100644 --- a/include/deal.II/dofs/block_info.h +++ b/include/deal.II/dofs/block_info.h @@ -31,11 +31,6 @@ DEAL_II_NAMESPACE_OPEN #ifndef DOXYGEN template class DoFHandler; -namespace hp -{ - template - class DoFHandler; -} #endif @@ -117,17 +112,6 @@ public: bool levels_only = false, bool active_only = false); - /** - * The same as above but for hp::DoFHandler. - * - * @note Not implemented. - */ - template - void - initialize(const hp::DoFHandler &, - bool levels_only = false, - bool active_only = false); - /** * @brief Initialize block structure on cells and compute renumbering * between cell dofs and block cell dofs. @@ -136,15 +120,6 @@ public: void initialize_local(const DoFHandler &); - /** - * The same as above but for hp::DoFHandler. - * - * @note Not implemented. - */ - template - void - initialize_local(const hp::DoFHandler &); - /** * Access the BlockIndices structure of the global system. */ diff --git a/include/deal.II/dofs/dof_accessor.h b/include/deal.II/dofs/dof_accessor.h index cb4b69fb75..e5a7634cdb 100644 --- a/include/deal.II/dofs/dof_accessor.h +++ b/include/deal.II/dofs/dof_accessor.h @@ -20,6 +20,7 @@ #include #include +#include #include @@ -121,6 +122,8 @@ namespace internal */ using BaseClass = dealii::CellAccessor; }; + + struct Implementation; } // namespace DoFAccessorImplementation } // namespace internal @@ -419,7 +422,7 @@ public: void get_dof_indices( std::vector &dof_indices, - const unsigned int fe_index = DoFHandlerType::default_fe_index) const; + const unsigned int fe_index = DoFHandlerType::invalid_fe_index) const; /** * Return the global multilevel indices of the degrees of freedom that live @@ -431,7 +434,7 @@ public: get_mg_dof_indices( const int level, std::vector &dof_indices, - const unsigned int fe_index = DoFHandlerType::default_fe_index) const; + const unsigned int fe_index = DoFHandlerType::invalid_fe_index) const; /** * Set the level DoF indices that are returned by get_mg_dof_indices. @@ -440,7 +443,7 @@ public: set_mg_dof_indices( const int level, const std::vector &dof_indices, - const unsigned int fe_index = DoFHandlerType::default_fe_index); + const unsigned int fe_index = DoFHandlerType::invalid_fe_index); /** * Global DoF index of the i degree associated with the @p vertexth @@ -463,7 +466,7 @@ public: vertex_dof_index( const unsigned int vertex, const unsigned int i, - const unsigned int fe_index = DoFHandlerType::default_fe_index) const; + const unsigned int fe_index = DoFHandlerType::invalid_fe_index) const; /** * Return the global DoF index of the ith degree of freedom @@ -475,7 +478,7 @@ public: const int level, const unsigned int vertex, const unsigned int i, - const unsigned int fe_index = DoFHandlerType::default_fe_index) const; + const unsigned int fe_index = DoFHandlerType::invalid_fe_index) const; /** * Index of the ith degree of freedom of this object. @@ -507,7 +510,7 @@ public: types::global_dof_index dof_index( const unsigned int i, - const unsigned int fe_index = DoFHandlerType::default_fe_index) const; + const unsigned int fe_index = DoFHandlerType::invalid_fe_index) const; /** * Return the dof_index on the given level. Also see dof_index. @@ -686,7 +689,7 @@ protected: set_dof_index( const unsigned int i, const types::global_dof_index index, - const unsigned int fe_index = DoFHandlerType::default_fe_index) const; + const unsigned int fe_index = DoFHandlerType::invalid_fe_index) const; void set_mg_dof_index(const int level, @@ -716,7 +719,7 @@ protected: const unsigned int vertex, const unsigned int i, const types::global_dof_index index, - const unsigned int fe_index = DoFHandlerType::default_fe_index) const; + const unsigned int fe_index = DoFHandlerType::invalid_fe_index) const; void set_mg_vertex_dof_index( @@ -724,7 +727,7 @@ protected: const unsigned int vertex, const unsigned int i, const types::global_dof_index index, - const unsigned int fe_index = DoFHandlerType::default_fe_index) const; + const unsigned int fe_index = DoFHandlerType::invalid_fe_index) const; // Iterator classes need to be friends because they need to access // operator== and operator!=. @@ -973,7 +976,7 @@ public: void get_dof_indices( std::vector &dof_indices, - const unsigned int fe_index = AccessorData::default_fe_index) const; + const unsigned int fe_index = AccessorData::invalid_fe_index) const; /** * Return the global multilevel indices of the degrees of freedom that live @@ -985,7 +988,7 @@ public: get_mg_dof_indices( const int level, std::vector &dof_indices, - const unsigned int fe_index = AccessorData::default_fe_index) const; + const unsigned int fe_index = AccessorData::invalid_fe_index) const; /** * Global DoF index of the i degree associated with the @p vertexth @@ -1009,7 +1012,7 @@ public: vertex_dof_index( const unsigned int vertex, const unsigned int i, - const unsigned int fe_index = AccessorData::default_fe_index) const; + const unsigned int fe_index = AccessorData::invalid_fe_index) const; /** * Index of the ith degree of freedom of this object. @@ -1030,7 +1033,7 @@ public: */ types::global_dof_index dof_index(const unsigned int i, - const unsigned int fe_index = AccessorData::default_fe_index) const; + const unsigned int fe_index = AccessorData::invalid_fe_index) const; /** * @} @@ -1177,7 +1180,7 @@ protected: set_dof_index( const unsigned int i, const types::global_dof_index index, - const unsigned int fe_index = AccessorData::default_fe_index) const; + const unsigned int fe_index = AccessorData::invalid_fe_index) const; /** * Set the global index of the i degree on the @p vertex-th vertex of @@ -1202,7 +1205,7 @@ protected: const unsigned int vertex, const unsigned int i, const types::global_dof_index index, - const unsigned int fe_index = AccessorData::default_fe_index) const; + const unsigned int fe_index = AccessorData::invalid_fe_index) const; // Iterator classes need to be friends because they need to access // operator== and operator!=. @@ -1306,7 +1309,7 @@ public: set_dof_index(const unsigned int i, const types::global_dof_index index, const unsigned int fe_index = - DoFHandler::default_fe_index) const; + DoFHandler::invalid_fe_index) const; }; @@ -1646,7 +1649,7 @@ public: get_interpolated_dof_values( const InputVector &values, Vector & interpolated_values, - const unsigned int fe_index = DoFHandlerType::default_fe_index) const; + const unsigned int fe_index = DoFHandlerType::invalid_fe_index) const; /** * This function is the counterpart to get_interpolated_dof_values(): you @@ -1705,7 +1708,7 @@ public: set_dof_values_by_interpolation( const Vector &local_values, OutputVector & values, - const unsigned int fe_index = DoFHandlerType::default_fe_index) const; + const unsigned int fe_index = DoFHandlerType::invalid_fe_index) const; /** * Distribute a local (cell based) vector to a global one by mapping the diff --git a/include/deal.II/dofs/dof_accessor.templates.h b/include/deal.II/dofs/dof_accessor.templates.h index 4e1e865e3c..4d52fa12bc 100644 --- a/include/deal.II/dofs/dof_accessor.templates.h +++ b/include/deal.II/dofs/dof_accessor.templates.h @@ -28,9 +28,6 @@ #include #include -#include -#include - #include #include @@ -191,8 +188,8 @@ inline TriaIterator> DoFAccessor::child( const unsigned int i) const { - Assert(static_cast(this->level()) < - this->dof_handler->levels.size(), + Assert(static_cast(this->present_level) < + this->dof_handler->object_dof_indices.size(), ExcMessage("DoFHandler not initialized")); TriaIterator - static types::global_dof_index - get_dof_index(const dealii::DoFHandler<1, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const unsigned int local_index, - std::integral_constant) - { - return dof_handler.levels[obj_level]->dof_object.get_dof_index( - dof_handler, obj_index, fe_index, local_index); - } - - - template - static void - set_dof_index(const dealii::DoFHandler<1, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const unsigned int local_index, - std::integral_constant, - const types::global_dof_index global_index) - { - dof_handler.levels[obj_level]->dof_object.set_dof_index( - dof_handler, obj_index, fe_index, local_index, global_index); - } - - - template - static types::global_dof_index - get_dof_index(const dealii::DoFHandler<2, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const unsigned int local_index, - std::integral_constant) - { - (void)obj_level; - // faces have no levels - Assert(obj_level == 0, ExcInternalError()); - return dof_handler.faces->lines.get_dof_index(dof_handler, - obj_index, - fe_index, - local_index); - } - - - template - static void - set_dof_index(const dealii::DoFHandler<2, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const unsigned int local_index, - std::integral_constant, - const types::global_dof_index global_index) - { - (void)obj_level; - // faces have no levels - Assert(obj_level == 0, ExcInternalError()); - dof_handler.faces->lines.set_dof_index( - dof_handler, obj_index, fe_index, local_index, global_index); - } - - - template - static types::global_dof_index - get_dof_index(const dealii::DoFHandler<2, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const unsigned int local_index, - std::integral_constant) - { - return dof_handler.levels[obj_level]->dof_object.get_dof_index( - dof_handler, obj_index, fe_index, local_index); - } - - - template - static void - set_dof_index(const dealii::DoFHandler<2, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const unsigned int local_index, - std::integral_constant, - const types::global_dof_index global_index) - { - dof_handler.levels[obj_level]->dof_object.set_dof_index( - dof_handler, obj_index, fe_index, local_index, global_index); - } - - - template - static types::global_dof_index - get_dof_index(const dealii::DoFHandler<3, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const unsigned int local_index, - std::integral_constant) - { - (void)obj_level; - // faces have no levels - Assert(obj_level == 0, ExcInternalError()); - return dof_handler.faces->lines.get_dof_index(dof_handler, - obj_index, - fe_index, - local_index); - } - - - template - static void - set_dof_index(const dealii::DoFHandler<3, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const unsigned int local_index, - std::integral_constant, - const types::global_dof_index global_index) - { - (void)obj_level; - // faces have no levels - Assert(obj_level == 0, ExcInternalError()); - dof_handler.faces->lines.set_dof_index( - dof_handler, obj_index, fe_index, local_index, global_index); - } - - - - template - static types::global_dof_index - get_dof_index(const dealii::DoFHandler<3, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const unsigned int local_index, - std::integral_constant) - { - (void)obj_level; - // faces have no levels - Assert(obj_level == 0, ExcInternalError()); - return dof_handler.faces->quads.get_dof_index(dof_handler, - obj_index, - fe_index, - local_index); - } - - - template - static void - set_dof_index(const dealii::DoFHandler<3, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const unsigned int local_index, - std::integral_constant, - const types::global_dof_index global_index) + template + static const types::global_dof_index * + get_cache_ptr(DoFHandlerType * dof_handler, + const unsigned int present_level, + const unsigned int present_index, + const unsigned int dofs_per_cell) { - (void)obj_level; - // faces have no levels - Assert(obj_level == 0, ExcInternalError()); - dof_handler.faces->quads.set_dof_index( - dof_handler, obj_index, fe_index, local_index, global_index); - } - + (void)dofs_per_cell; - - template - static types::global_dof_index - get_dof_index(const dealii::DoFHandler<3, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const unsigned int local_index, - std::integral_constant) - { - return dof_handler.levels[obj_level]->dof_object.get_dof_index( - dof_handler, obj_index, fe_index, local_index); + return &dof_handler + ->cell_dof_cache_indices[present_level] + [dof_handler->cell_dof_cache_ptr + [present_level][present_index]]; } - template - static void - set_dof_index(const dealii::DoFHandler<3, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const unsigned int local_index, - std::integral_constant, - const types::global_dof_index global_index) - { - dof_handler.levels[obj_level]->dof_object.set_dof_index( - dof_handler, obj_index, fe_index, local_index, global_index); - } - - template - static types::global_dof_index - get_dof_index(const dealii::hp::DoFHandler<1, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const unsigned int local_index, - const std::integral_constant &) - { - return dof_handler.levels[obj_level]->get_dof_index(obj_index, - fe_index, - local_index); - } - - - template + /** + * Set the @p local_index-th degree of freedom corresponding to the + * finite element specified by @p fe_index on the vertex with global + * number @p vertex_index to @p global_index. + */ + template static void - set_dof_index(const dealii::hp::DoFHandler<1, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const unsigned int local_index, - const std::integral_constant &, + set_dof_index(const DoFHandlerType &dof_handler, + const unsigned int obj_level, + const unsigned int obj_index, + const unsigned int fe_index, + const unsigned int local_index, + const std::integral_constant &, const types::global_dof_index global_index) { - dof_handler.levels[obj_level]->set_dof_index(obj_index, - fe_index, - local_index, - global_index); - } + Assert(d == DoFHandlerType::dimension || obj_level == 0, + ExcNotImplemented()); + // 1) no hp used -> fe_index == 0 + if (dof_handler.hp_capability_enabled == false) + { + AssertDimension(fe_index, DoFHandlerType::default_fe_index); - template - static types::global_dof_index - get_dof_index(const dealii::hp::DoFHandler<2, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const unsigned int local_index, - const std::integral_constant &) - { - return dof_handler.faces->lines.get_dof_index( - dof_handler, obj_index, fe_index, local_index, obj_level); - } + if (d == 0) + dof_handler.object_dof_indices + [0][0][obj_index * dof_handler.get_fe().dofs_per_vertex + + local_index] = global_index; + else + dof_handler.object_dof_indices + [obj_level][d] + [dof_handler.object_dof_ptr[obj_level][d][obj_index] + + local_index] = global_index; + return; + } - template - static void - set_dof_index(const dealii::hp::DoFHandler<2, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const unsigned int local_index, - const std::integral_constant &, - const types::global_dof_index global_index) - { - dof_handler.faces->lines.set_dof_index(dof_handler, - obj_index, - fe_index, - local_index, - global_index, - obj_level); - } + // 2) cell and hp is used -> there is only one fe_index + if (d == DoFHandlerType::dimension) + { + dof_handler.object_dof_indices + [obj_level][d] + [dof_handler.object_dof_ptr[obj_level][d][obj_index] + + local_index] = global_index; + return; + } + // 3) general entity and hp is used + AssertIndexRange(obj_level, dof_handler.object_dof_indices.size()); + AssertIndexRange(d, dof_handler.object_dof_indices[obj_level].size()); - template - static types::global_dof_index - get_dof_index(const dealii::hp::DoFHandler<2, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const unsigned int local_index, - const std::integral_constant &) - { - return dof_handler.levels[obj_level]->get_dof_index(obj_index, - fe_index, - local_index); - } + unsigned int fe_index_; + if (dof_handler.hp_capability_enabled) + { + AssertIndexRange(d, dof_handler.hp_object_fe_ptr.size()); + AssertIndexRange(obj_index, dof_handler.hp_object_fe_ptr[d].size()); + + const auto ptr = + std::find(dof_handler.hp_object_fe_indices[d].begin() + + dof_handler.hp_object_fe_ptr[d][obj_index], + dof_handler.hp_object_fe_indices[d].begin() + + dof_handler.hp_object_fe_ptr[d][obj_index + 1], + fe_index); + + Assert(ptr != dof_handler.hp_object_fe_indices[d].begin() + + dof_handler.hp_object_fe_ptr[d][obj_index + 1], + ExcNotImplemented()); + + fe_index_ = + std::distance(dof_handler.hp_object_fe_indices[d].begin() + + dof_handler.hp_object_fe_ptr[d][obj_index], + ptr); + } - template - static void - set_dof_index(const dealii::hp::DoFHandler<2, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const unsigned int local_index, - const std::integral_constant &, - const types::global_dof_index global_index) - { - dof_handler.levels[obj_level]->set_dof_index(obj_index, - fe_index, - local_index, - global_index); + AssertIndexRange(dof_handler.hp_capability_enabled ? + (dof_handler.hp_object_fe_ptr[d][obj_index] + + fe_index_) : + obj_index, + dof_handler.object_dof_ptr[obj_level][d].size()); + + AssertIndexRange( + dof_handler + .object_dof_ptr[obj_level][d] + [dof_handler.hp_capability_enabled ? + (dof_handler.hp_object_fe_ptr[d][obj_index] + + fe_index_) : + obj_index] + + local_index, + dof_handler.object_dof_indices[obj_level][d].size()); + + dof_handler.object_dof_indices + [obj_level][d] + [dof_handler.object_dof_ptr + [obj_level][d] + [dof_handler.hp_capability_enabled ? + (dof_handler.hp_object_fe_ptr[d][obj_index] + fe_index_) : + obj_index] + + local_index] = global_index; } - template + /** + * Get the @p local_index-th degree of freedom corresponding to the + * finite element specified by @p fe_index on the vertex with global + * number @p vertex_index to @p global_index. + */ + template static types::global_dof_index - get_dof_index(const dealii::hp::DoFHandler<3, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const unsigned int local_index, - const std::integral_constant &) - { - return dof_handler.faces->lines.get_dof_index( - dof_handler, obj_index, fe_index, local_index, obj_level); - } - - - template - static void - set_dof_index(const dealii::hp::DoFHandler<3, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const unsigned int local_index, - const std::integral_constant &, - const types::global_dof_index global_index) - { - dof_handler.faces->lines.set_dof_index(dof_handler, - obj_index, - fe_index, - local_index, - global_index, - obj_level); - } + get_dof_index(const DoFHandlerType &dof_handler, + const unsigned int obj_level, + const unsigned int obj_index, + const unsigned int fe_index, + const unsigned int local_index, + const std::integral_constant &) + { + Assert(d == DoFHandlerType::dimension || obj_level == 0, + ExcNotImplemented()); + + // 1) no hp used -> fe_index == 0 + if (dof_handler.hp_capability_enabled == false) + { + AssertDimension(fe_index, DoFHandlerType::default_fe_index); + if (d == 0) + return dof_handler.object_dof_indices + [0][0][obj_index * dof_handler.get_fe().dofs_per_vertex + + local_index]; + else + return dof_handler.object_dof_indices + [obj_level][d] + [dof_handler.object_dof_ptr[obj_level][d][obj_index] + + local_index]; + } - template - static types::global_dof_index - get_dof_index(const dealii::hp::DoFHandler<3, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const unsigned int local_index, - const std::integral_constant &) - { - return dof_handler.faces->quads.get_dof_index( - dof_handler, obj_index, fe_index, local_index, obj_level); - } + // 2) cell and hp is used -> there is only one fe_index + if (d == DoFHandlerType::dimension) + { + return dof_handler.object_dof_indices + [obj_level][d] + [dof_handler.object_dof_ptr[obj_level][d][obj_index] + + local_index]; + } + // 3) general entity and hp is used + AssertIndexRange(obj_level, dof_handler.object_dof_indices.size()); + AssertIndexRange(d, dof_handler.object_dof_indices[obj_level].size()); - template - static void - set_dof_index(const dealii::hp::DoFHandler<3, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const unsigned int local_index, - const std::integral_constant &, - const types::global_dof_index global_index) - { - dof_handler.faces->quads.set_dof_index(dof_handler, - obj_index, - fe_index, - local_index, - global_index, - obj_level); - } + unsigned int fe_index_; + if (dof_handler.hp_capability_enabled) + { + AssertIndexRange(d, dof_handler.hp_object_fe_ptr.size()); + AssertIndexRange(obj_index, dof_handler.hp_object_fe_ptr[d].size()); + + const auto ptr = + std::find(dof_handler.hp_object_fe_indices[d].begin() + + dof_handler.hp_object_fe_ptr[d][obj_index], + dof_handler.hp_object_fe_indices[d].begin() + + dof_handler.hp_object_fe_ptr[d][obj_index + 1], + fe_index); + + Assert(ptr != dof_handler.hp_object_fe_indices[d].begin() + + dof_handler.hp_object_fe_ptr[d][obj_index + 1], + ExcNotImplemented()); + + fe_index_ = + std::distance(dof_handler.hp_object_fe_indices[d].begin() + + dof_handler.hp_object_fe_ptr[d][obj_index], + ptr); + } - template + AssertIndexRange(dof_handler.hp_capability_enabled ? + (dof_handler.hp_object_fe_ptr[d][obj_index] + + fe_index_) : + obj_index, + dof_handler.object_dof_ptr[obj_level][d].size()); + + AssertIndexRange( + dof_handler + .object_dof_ptr[obj_level][d] + [dof_handler.hp_capability_enabled ? + (dof_handler.hp_object_fe_ptr[d][obj_index] + + fe_index_) : + obj_index] + + local_index, + dof_handler.object_dof_indices[obj_level][d].size()); + + return dof_handler.object_dof_indices + [obj_level][d] + [dof_handler.object_dof_ptr + [obj_level][d] + [dof_handler.hp_capability_enabled ? + (dof_handler.hp_object_fe_ptr[d][obj_index] + fe_index_) : + obj_index] + + local_index]; + } + + + template static types::global_dof_index - get_dof_index(const dealii::hp::DoFHandler<3, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const unsigned int local_index, - const std::integral_constant &) - { - return dof_handler.levels[obj_level]->get_dof_index(obj_index, - fe_index, - local_index); - } - - - template - static void - set_dof_index(const dealii::hp::DoFHandler<3, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const unsigned int local_index, - const std::integral_constant &, - const types::global_dof_index global_index) + mg_vertex_dof_index(const DoFHandlerType &dof_handler, + const int level, + const unsigned int vertex_index, + const unsigned int i) { - dof_handler.levels[obj_level]->set_dof_index(obj_index, - fe_index, - local_index, - global_index); - } - + Assert(dof_handler.hp_capability_enabled == false, + ExcMessage( + "hp::DoFHandler does not implement multilevel DoFs.")); - template - static types::global_dof_index - mg_vertex_dof_index(const dealii::DoFHandler &dof_handler, - const int level, - const unsigned int vertex_index, - const unsigned int i) - { return dof_handler.mg_vertex_dofs[vertex_index].get_index( level, i, dof_handler.get_fe().dofs_per_vertex); } - template - static types::global_dof_index - mg_vertex_dof_index(const dealii::hp::DoFHandler &, - const int, - const unsigned int, - const unsigned int) - { - Assert(false, - ExcMessage( - "hp::DoFHandler does not implement multilevel DoFs.")); - return numbers::invalid_dof_index; - } - - - template - static void - set_mg_vertex_dof_index(dealii::DoFHandler &dof_handler, - const int level, - const unsigned int vertex_index, - const unsigned int i, - types::global_dof_index index) - { - return dof_handler.mg_vertex_dofs[vertex_index].set_index( - level, i, dof_handler.get_fe().dofs_per_vertex, index); - } - - template + template static void - set_mg_vertex_dof_index(dealii::hp::DoFHandler &, - const int, - const unsigned int, - const unsigned int, - types::global_dof_index) + set_mg_vertex_dof_index(DoFHandlerType & dof_handler, + const int level, + const unsigned int vertex_index, + const unsigned int i, + types::global_dof_index index) { - Assert(false, + Assert(dof_handler.hp_capability_enabled == false, ExcMessage( "hp::DoFHandler does not implement multilevel DoFs.")); - } - - - - template - static bool - fe_index_is_active(const dealii::DoFHandler &, - const unsigned int, - const unsigned int, - const unsigned int fe_index, - const std::integral_constant &) - { - return (fe_index == 0); - } - - - - template - static unsigned int - n_active_fe_indices(const dealii::DoFHandler &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const std::integral_constant &) - { - (void)dof_handler; - (void)obj_level; - (void)obj_index; - // check that the object we look - // at is in fact active. the - // problem is that we have - // templatized on the - // dimensionality of the object, - // so it may be a cell, a face, - // or a line. we have a bit of - // trouble doing this all in the - // generic case, so only check if - // it is either a cell or a - // line. the only case this - // leaves out is faces in 3d -- - // let's hope that this never is - // a problem - Assert((dim == structdim ? - TriaRawIterator>( - &dof_handler.get_triangulation(), obj_level, obj_index) - ->used() : - (structdim == 1 ? - typename internal::TriangulationImplementation:: - Iterators::raw_line_iterator( - &dof_handler.get_triangulation(), obj_level, obj_index) - ->used() : - true)) == true, - ExcMessage("This cell is not active and therefore can't be " - "queried for its active FE indices")); - return 1; - } - - - - template - static unsigned int - nth_active_fe_index(const dealii::DoFHandler &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int n, - const std::integral_constant &) - { - (void)dof_handler; - (void)obj_level; - (void)obj_index; - (void)n; - // check that the object we look - // at is in fact active. the - // problem is that we have - // templatized on the - // dimensionality of the object, - // so it may be a cell, a face, - // or a line. we have a bit of - // trouble doing this all in the - // generic case, so only check if - // it is either a cell or a - // line. the only case this - // leaves out is faces in 3d -- - // let's hope that this never is - // a problem - Assert((dim == structdim ? - TriaRawIterator>( - &dof_handler.get_triangulation(), obj_level, obj_index) - ->used() : - (structdim == 1 ? - typename internal::TriangulationImplementation:: - Iterators::raw_line_iterator( - &dof_handler.get_triangulation(), obj_level, obj_index) - ->used() : - true)) == true, - ExcMessage("This cell is not active and therefore can't be " - "queried for its active FE indices")); - AssertIndexRange(n, 1); - - return dealii::DoFHandler::default_fe_index; - } - - - template - static bool - fe_index_is_active(const dealii::hp::DoFHandler<1, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const std::integral_constant &) - { - return dof_handler.levels[obj_level]->fe_index_is_active(obj_index, - fe_index); - } - - - template - static unsigned int - n_active_fe_indices(const dealii::hp::DoFHandler<1, spacedim> &, - const unsigned int /*obj_level*/, - const unsigned int /*obj_index*/, - const std::integral_constant &) - { - // on a cell, the number of active elements is one - return 1; - } - - - - template - static unsigned int - nth_active_fe_index( - const dealii::hp::DoFHandler<1, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int n, - const std::integral_constant &) - { - (void)n; - Assert(n == 0, - ExcMessage("On cells, there can only be one active FE index")); - return dof_handler.levels[obj_level]->active_fe_index(obj_index); - } - - - template - static bool - fe_index_is_active(const dealii::hp::DoFHandler<2, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const std::integral_constant &) - { - return dof_handler.faces->lines.fe_index_is_active(dof_handler, - obj_index, - fe_index, - obj_level); - } - - - template - static unsigned int - n_active_fe_indices( - const dealii::hp::DoFHandler<2, spacedim> &dof_handler, - const unsigned int, - const unsigned int obj_index, - const std::integral_constant &) - { - return dof_handler.faces->lines.n_active_fe_indices(dof_handler, - obj_index); - } - - - template - static unsigned int - nth_active_fe_index( - const dealii::hp::DoFHandler<2, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int n, - const std::integral_constant &) - { - return dof_handler.faces->lines.nth_active_fe_index(dof_handler, - obj_level, - obj_index, - n); - } - - - - template - static bool - fe_index_is_active(const dealii::hp::DoFHandler<2, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const std::integral_constant &) - { - return dof_handler.levels[obj_level]->fe_index_is_active(obj_index, - fe_index); - } - - - template - static unsigned int - n_active_fe_indices(const dealii::hp::DoFHandler<2, spacedim> &, - const unsigned int /*obj_level*/, - const unsigned int /*obj_index*/, - const std::integral_constant &) - { - // on a cell, the number of active elements is one - return 1; - } - - - - template - static unsigned int - nth_active_fe_index( - const dealii::hp::DoFHandler<2, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int n, - const std::integral_constant &) - { - (void)n; - Assert(n == 0, - ExcMessage("On cells, there can only be one active FE index")); - return dof_handler.levels[obj_level]->active_fe_index(obj_index); - } - - - - template - static bool - fe_index_is_active(const dealii::hp::DoFHandler<3, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const std::integral_constant &) - { - return dof_handler.faces->lines.fe_index_is_active(dof_handler, - obj_index, - fe_index, - obj_level); - } - - template - static unsigned int - n_active_fe_indices( - const dealii::hp::DoFHandler<3, spacedim> &dof_handler, - const unsigned int, - const unsigned int obj_index, - const std::integral_constant &) - { - return dof_handler.faces->lines.n_active_fe_indices(dof_handler, - obj_index); - } - - - - template - static unsigned int - nth_active_fe_index( - const dealii::hp::DoFHandler<3, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int n, - const std::integral_constant &) - { - return dof_handler.faces->lines.nth_active_fe_index(dof_handler, - obj_level, - obj_index, - n); - } - - - - template - static bool - fe_index_is_active(const dealii::hp::DoFHandler<3, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const std::integral_constant &) - { - return dof_handler.faces->quads.fe_index_is_active(dof_handler, - obj_index, - fe_index, - obj_level); - } - - template - static bool - fe_index_is_active(const dealii::hp::DoFHandler<3, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int fe_index, - const std::integral_constant &) - { - return dof_handler.levels[obj_level]->fe_index_is_active(obj_index, - fe_index); - } - - - template - static unsigned int - n_active_fe_indices( - const dealii::hp::DoFHandler<3, spacedim> &dof_handler, - const unsigned int, - const unsigned int obj_index, - const std::integral_constant &) - { - return dof_handler.faces->quads.n_active_fe_indices(dof_handler, - obj_index); - } - - - - template - static unsigned int - nth_active_fe_index( - const dealii::hp::DoFHandler<3, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int n, - const std::integral_constant &) - { - return dof_handler.faces->quads.nth_active_fe_index(dof_handler, - obj_level, - obj_index, - n); - } - - - - template - static unsigned int - n_active_fe_indices(const dealii::hp::DoFHandler<3, spacedim> &, - const unsigned int /*obj_level*/, - const unsigned int /*obj_index*/, - const std::integral_constant &) - { - // on a cell, the number of active elements is one - return 1; - } - - - - template - static unsigned int - nth_active_fe_index( - const dealii::hp::DoFHandler<3, spacedim> &dof_handler, - const unsigned int obj_level, - const unsigned int obj_index, - const unsigned int n, - const std::integral_constant &) - { - (void)n; - Assert(n == 0, - ExcMessage("On cells, there can only be one active FE index")); - return dof_handler.levels[obj_level]->active_fe_index(obj_index); - } - - /** - * Set the @p local_index-th degree of freedom corresponding to the - * finite element specified by @p fe_index on the vertex with global - * number @p vertex_index to @p global_index. - */ - template - static void - set_vertex_dof_index(dealii::DoFHandler &dof_handler, - const unsigned int vertex_index, - const unsigned int fe_index, - const unsigned int local_index, - const types::global_dof_index global_index) - { - (void)fe_index; - Assert( - (fe_index == dealii::DoFHandler::default_fe_index), - ExcMessage( - "Only the default FE index is allowed for non-hp DoFHandler objects")); - AssertIndexRange(local_index, dof_handler.get_fe().dofs_per_vertex); - - dof_handler - .vertex_dofs[vertex_index * dof_handler.get_fe().dofs_per_vertex + - local_index] = global_index; - } - - - template - static void - set_vertex_dof_index(dealii::hp::DoFHandler &dof_handler, - const unsigned int vertex_index, - const unsigned int fe_index, - const unsigned int local_index, - const types::global_dof_index global_index) - { - Assert((fe_index != - dealii::hp::DoFHandler::default_fe_index), - ExcMessage("You need to specify a FE index when working " - "with hp DoFHandlers")); - Assert(dof_handler.fe_collection.size() > 0, - ExcMessage("No finite element collection is associated with " - "this DoFHandler")); - AssertIndexRange(local_index, - dof_handler.get_fe(fe_index).dofs_per_vertex); - Assert(fe_index < dof_handler.fe_collection.size(), ExcInternalError()); - Assert(dof_handler.vertex_dof_offsets[vertex_index] != - numbers::invalid_unsigned_int, - ExcMessage( - "This vertex is unused and has no DoFs associated with it")); - - // hop along the list of index - // sets until we find the one - // with the correct fe_index, and - // then poke into that - // part. trigger an exception if - // we can't find a set for this - // particular fe_index - const unsigned int starting_offset = - dof_handler.vertex_dof_offsets[vertex_index]; - types::global_dof_index *pointer = - &dof_handler.vertex_dofs[starting_offset]; - while (true) - { - Assert(pointer <= &dof_handler.vertex_dofs.back(), - ExcInternalError()); - - // a fe index is always small - Assert((*pointer) < std::numeric_limits::max(), - ExcInternalError()); - const types::global_dof_index this_fe_index = *pointer; - - Assert(this_fe_index != numbers::invalid_dof_index, - ExcInternalError()); - Assert(this_fe_index < dof_handler.fe_collection.size(), - ExcInternalError()); - - if (this_fe_index == fe_index) - { - *(pointer + 1 + local_index) = global_index; - return; - } - else - pointer += static_cast( - dof_handler.get_fe(this_fe_index).dofs_per_vertex + 1); - } - } - - - /** - * Get the @p local_index-th degree of freedom corresponding to the - * finite element specified by @p fe_index on the vertex with global - * number @p vertex_index to @p global_index. - */ - - template - static types::global_dof_index - get_vertex_dof_index(const dealii::DoFHandler &dof_handler, - const unsigned int vertex_index, - const unsigned int fe_index, - const unsigned int local_index) - { - (void)fe_index; - Assert( - (fe_index == dealii::DoFHandler::default_fe_index), - ExcMessage( - "Only the default FE index is allowed for non-hp DoFHandler objects")); - AssertIndexRange(local_index, dof_handler.get_fe().dofs_per_vertex); - - return dof_handler - .vertex_dofs[vertex_index * dof_handler.get_fe().dofs_per_vertex + - local_index]; - } - - - template - static types::global_dof_index - get_vertex_dof_index( - const dealii::hp::DoFHandler &dof_handler, - const unsigned int vertex_index, - const unsigned int fe_index, - const unsigned int local_index) - { - Assert((fe_index != - dealii::hp::DoFHandler::default_fe_index), - ExcMessage("You need to specify a FE index when working " - "with hp DoFHandlers")); - Assert(dof_handler.fe_collection.size() > 0, - ExcMessage("No finite element collection is associated with " - "this DoFHandler")); - AssertIndexRange(local_index, - dof_handler.get_fe(fe_index).dofs_per_vertex); - AssertIndexRange(vertex_index, dof_handler.vertex_dof_offsets.size()); - Assert(dof_handler.vertex_dof_offsets[vertex_index] != - numbers::invalid_unsigned_int, - ExcMessage( - "This vertex is unused and has no DoFs associated with it")); - - // hop along the list of index - // sets until we find the one - // with the correct fe_index, and - // then poke into that - // part. trigger an exception if - // we can't find a set for this - // particular fe_index - const unsigned int starting_offset = - dof_handler.vertex_dof_offsets[vertex_index]; - const types::global_dof_index *pointer = - &dof_handler.vertex_dofs[starting_offset]; - while (true) - { - Assert(pointer <= &dof_handler.vertex_dofs.back(), - ExcInternalError()); - - Assert((*pointer) < - std::numeric_limits::max(), - ExcInternalError()); - const types::global_dof_index this_fe_index = *pointer; - - Assert(this_fe_index != numbers::invalid_dof_index, - ExcInternalError()); - Assert(this_fe_index < dof_handler.fe_collection.size(), - ExcInternalError()); - - if (this_fe_index == fe_index) - return *(pointer + 1 + local_index); - else - pointer += static_cast( - dof_handler.get_fe(this_fe_index).dofs_per_vertex + 1); - } + return dof_handler.mg_vertex_dofs[vertex_index].set_index( + level, i, dof_handler.get_fe().dofs_per_vertex, index); } @@ -1205,48 +468,32 @@ namespace internal * Return the number of different finite elements that are active on a * given vertex. */ - template + template static unsigned int - n_active_vertex_fe_indices( - const dealii::hp::DoFHandler &dof_handler, - const unsigned int vertex_index) + n_active_fe_indices(const DoFHandlerType &dof_handler, + const unsigned int obj_level, + const unsigned int obj_index, + const std::integral_constant &) { - Assert(dof_handler.fe_collection.size() > 0, - ExcMessage("No finite element collection is associated with " - "this DoFHandler")); - - // if this vertex is unused, return 0 - if (dof_handler.vertex_dof_offsets[vertex_index] == - numbers::invalid_unsigned_int) - return 0; - - // hop along the list of index - // sets and count the number of - // hops - const unsigned int starting_offset = - dof_handler.vertex_dof_offsets[vertex_index]; - const types::global_dof_index *pointer = - &dof_handler.vertex_dofs[starting_offset]; - - Assert(*pointer != numbers::invalid_dof_index, ExcInternalError()); - - unsigned int counter = 0; - while (true) - { - Assert(pointer <= &dof_handler.vertex_dofs.back(), - ExcInternalError()); + (void)obj_level; - const types::global_dof_index this_fe_index = *pointer; + Assert(d == DoFHandlerType::dimension || obj_level == 0, + ExcNotImplemented()); - if (this_fe_index == numbers::invalid_dof_index) - return counter; - else - { - pointer += static_cast( - dof_handler.get_fe(this_fe_index).dofs_per_vertex + 1); - ++counter; - } - } + // 1) no hp used -> fe_index == 0 + if (dof_handler.hp_capability_enabled == false) + return 1; + + // 2) cell and hp is used -> there is only one fe_index + if (d == DoFHandlerType::dimension) + return 1; + + // 3) general entity and hp is used + AssertIndexRange(d, dof_handler.hp_object_fe_ptr.size()); + AssertIndexRange(obj_index + 1, dof_handler.hp_object_fe_ptr[d].size()); + + return dof_handler.hp_object_fe_ptr[d][obj_index + 1] - + dof_handler.hp_object_fe_ptr[d][obj_index]; } @@ -1255,58 +502,41 @@ namespace internal * Return the fe index of the n-th finite element active on a given * vertex. */ - template + template static unsigned int - nth_active_vertex_fe_index( - const dealii::hp::DoFHandler &dof_handler, - const unsigned int vertex_index, - const unsigned int n) - { - Assert(dof_handler.fe_collection.size() > 0, - ExcMessage("No finite element collection is associated with " - "this DoFHandler")); - Assert(n < n_active_vertex_fe_indices(dof_handler, vertex_index), - ExcIndexRange( - n, 0, n_active_vertex_fe_indices(dof_handler, vertex_index))); - // make sure we don't ask on - // unused vertices - Assert(dof_handler.vertex_dof_offsets[vertex_index] != - numbers::invalid_unsigned_int, - ExcInternalError()); - - // hop along the list of index - // sets and count the number of - // hops - const unsigned int starting_offset = - dof_handler.vertex_dof_offsets[vertex_index]; - const types::global_dof_index *pointer = - &dof_handler.vertex_dofs[starting_offset]; - - Assert(*pointer != numbers::invalid_dof_index, ExcInternalError()); - - unsigned int counter = 0; - while (true) - { - Assert(pointer <= &dof_handler.vertex_dofs.back(), - ExcInternalError()); - - Assert((*pointer) < std::numeric_limits::max(), - ExcInternalError()); - const types::global_dof_index this_fe_index = *pointer; - - Assert(this_fe_index < dof_handler.fe_collection.size(), - ExcInternalError()); - - if (counter == n) - return this_fe_index; + nth_active_fe_index(const DoFHandlerType &dof_handler, + const unsigned int obj_level, + const unsigned int obj_index, + const unsigned int local_index, + const std::integral_constant &) + { + Assert(d == DoFHandlerType::dimension || obj_level == 0, + ExcNotImplemented()); + + // for cells only one active fe index available + Assert(((d == DoFHandlerType::dimension) && + (local_index != DoFHandlerType::default_fe_index)) == false, + ExcNotImplemented()); + + // 1) no hp used -> fe_index == 0 + if (dof_handler.hp_capability_enabled == false) + return DoFHandlerType::default_fe_index; + + // 2) cell and hp is used -> there is only one fe_index + if (d == DoFHandlerType::dimension) + return dof_handler.hp_cell_active_fe_indices[obj_level][obj_index]; + + // 3) general entity and hp is used + AssertIndexRange(d, dof_handler.hp_object_fe_indices.size()); + AssertIndexRange(d, dof_handler.hp_object_fe_ptr.size()); + AssertIndexRange(obj_index, dof_handler.hp_object_fe_ptr[d].size()); + AssertIndexRange(dof_handler.hp_object_fe_ptr[d][obj_index] + + local_index, + dof_handler.hp_object_fe_indices[d].size()); - Assert(this_fe_index != numbers::invalid_dof_index, - ExcInternalError()); - - pointer += static_cast( - dof_handler.get_fe(this_fe_index).dofs_per_vertex + 1); - ++counter; - } + return dof_handler + .hp_object_fe_indices[d][dof_handler.hp_object_fe_ptr[d][obj_index] + + local_index]; } @@ -1317,81 +547,68 @@ namespace internal * The size of the returned set equals the number of finite elements that * are active on this vertex. */ - template + template static std::set - get_active_vertex_fe_indices( - const dealii::hp::DoFHandler &dof_handler, - const unsigned int vertex_index) + get_active_fe_indices(const DoFHandlerType & dof_handler, + const unsigned int obj_level, + const unsigned int obj_index, + const std::integral_constant &t) { + Assert(d == DoFHandlerType::dimension || obj_level == 0, + ExcNotImplemented()); + + // 1) no hp used -> fe_index == 0 + if (dof_handler.hp_capability_enabled == false) + return {DoFHandlerType::default_fe_index}; + + // 2) cell and hp is used -> there is only one fe_index + if (d == DoFHandlerType::dimension) + return {dof_handler.hp_cell_active_fe_indices[obj_level][obj_index]}; + + // 3) general entity and hp is used std::set active_fe_indices; for (unsigned int i = 0; - i < n_active_vertex_fe_indices(dof_handler, vertex_index); + i < n_active_fe_indices(dof_handler, obj_level, obj_index, t); ++i) active_fe_indices.insert( - nth_active_vertex_fe_index(dof_handler, vertex_index, i)); + nth_active_fe_index(dof_handler, obj_level, obj_index, i, t)); return active_fe_indices; } - /** - * Return whether a particular finite element index is active on the - * specified vertex. - */ - template + template static bool - fe_is_active_on_vertex( - const dealii::hp::DoFHandler &dof_handler, - const unsigned int vertex_index, - const unsigned int fe_index) + fe_index_is_active(const DoFHandlerType &dof_handler, + const unsigned int obj_level, + const unsigned int obj_index, + const unsigned int fe_index, + const std::integral_constant &) { - Assert((fe_index != - dealii::hp::DoFHandler::default_fe_index), - ExcMessage("You need to specify a FE index when working " - "with hp DoFHandlers")); - Assert(dof_handler.fe_collection.size() > 0, - ExcMessage("No finite element collection is associated with " - "this DoFHandler")); - Assert(fe_index < dof_handler.fe_collection.size(), ExcInternalError()); - - // make sure we don't ask on - // unused vertices - Assert(dof_handler.vertex_dof_offsets[vertex_index] != - numbers::invalid_unsigned_int, - ExcInternalError()); - - // hop along the list of index - // sets and see whether we find - // the given index - const unsigned int starting_offset = - dof_handler.vertex_dof_offsets[vertex_index]; - const types::global_dof_index *pointer = - &dof_handler.vertex_dofs[starting_offset]; - - Assert(*pointer != numbers::invalid_dof_index, ExcInternalError()); - - while (true) - { - Assert(pointer <= &dof_handler.vertex_dofs.back(), - ExcInternalError()); + Assert(d == DoFHandlerType::dimension || obj_level == 0, + ExcNotImplemented()); - Assert((*pointer) < - std::numeric_limits::max(), - ExcInternalError()); - const types::global_dof_index this_fe_index = *pointer; + // 1) no hp used -> fe_index == 0 + if (dof_handler.hp_capability_enabled == false) + return (fe_index == DoFHandlerType::default_fe_index); - Assert(this_fe_index < dof_handler.fe_collection.size(), - ExcInternalError()); + // 2) cell and hp is used -> there is only one fe_index + if (d == DoFHandlerType::dimension) + return dof_handler.hp_cell_active_fe_indices[obj_level][obj_index] == + fe_index; - if (this_fe_index == numbers::invalid_dof_index) - return false; - else if (this_fe_index == fe_index) - return true; - else - pointer += dof_handler.get_fe(this_fe_index).dofs_per_vertex + 1; - } + // 3) general entity and hp is used + return std::find(dof_handler.hp_object_fe_indices[d].begin() + + dof_handler.hp_object_fe_ptr[d][obj_index], + dof_handler.hp_object_fe_indices[d].begin() + + dof_handler.hp_object_fe_ptr[d][obj_index + 1], + fe_index) != + (dof_handler.hp_object_fe_indices[d].begin() + + dof_handler.hp_object_fe_ptr[d][obj_index + 1]); } + + template static void set_mg_dof_indices( @@ -1844,8 +1061,14 @@ template inline types::global_dof_index DoFAccessor::dof_index( const unsigned int i, - const unsigned int fe_index) const + const unsigned int fe_index_) const { + const unsigned int fe_index = + (this->dof_handler->hp_capability_enabled == false && + fe_index_ == DoFHandlerType::invalid_fe_index) ? + DoFHandlerType::default_fe_index : + fe_index_; + // access the respective DoF return dealii::internal::DoFAccessorImplementation::Implementation:: get_dof_index(*this->dof_handler, @@ -1873,8 +1096,14 @@ inline void DoFAccessor::set_dof_index( const unsigned int i, const types::global_dof_index index, - const unsigned int fe_index) const + const unsigned int fe_index_) const { + const unsigned int fe_index = + (this->dof_handler->hp_capability_enabled == false && + fe_index_ == DoFHandlerType::invalid_fe_index) ? + DoFHandlerType::default_fe_index : + fe_index_; + // access the respective DoF dealii::internal::DoFAccessorImplementation::Implementation::set_dof_index( *this->dof_handler, @@ -1952,13 +1181,21 @@ inline types::global_dof_index DoFAccessor::vertex_dof_index( const unsigned int vertex, const unsigned int i, - const unsigned int fe_index) const + const unsigned int fe_index_) const { + const unsigned int fe_index = + (this->dof_handler->hp_capability_enabled == false && + fe_index_ == DoFHandlerType::invalid_fe_index) ? + DoFHandlerType::default_fe_index : + fe_index_; + return dealii::internal::DoFAccessorImplementation::Implementation:: - get_vertex_dof_index(*this->dof_handler, - this->vertex_index(vertex), - fe_index, - i); + get_dof_index(*this->dof_handler, + 0, + this->vertex_index(vertex), + fe_index, + i, + std::integral_constant()); } @@ -1968,18 +1205,23 @@ DoFAccessor::mg_vertex_dof_index( const int level, const unsigned int vertex, const unsigned int i, - const unsigned int fe_index) const + const unsigned int fe_index_) const { + const unsigned int fe_index = + (this->dof_handler->hp_capability_enabled == false && + fe_index_ == DoFHandlerType::invalid_fe_index) ? + DoFHandlerType::default_fe_index : + fe_index_; (void)fe_index; Assert(this->dof_handler != nullptr, ExcInvalidObject()); AssertIndexRange(vertex, GeometryInfo::vertices_per_cell); AssertIndexRange(i, this->dof_handler->get_fe(fe_index).dofs_per_vertex); - return dealii::internal::DoFAccessorImplementation::Implementation:: - mg_vertex_dof_index(*this->dof_handler, - level, - this->vertex_index(vertex), - i); + Assert(dof_handler->hp_capability_enabled == false, + ExcMessage("hp::DoFHandler does not implement multilevel DoFs.")); + + return this->dof_handler->mg_vertex_dofs[this->vertex_index(vertex)] + .get_index(level, i, this->dof_handler->get_fe().dofs_per_vertex); } @@ -1989,11 +1231,22 @@ DoFAccessor::set_vertex_dof_index( const unsigned int vertex, const unsigned int i, const types::global_dof_index index, - const unsigned int fe_index) const + const unsigned int fe_index_) const { - dealii::internal::DoFAccessorImplementation::Implementation:: - set_vertex_dof_index( - *this->dof_handler, this->vertex_index(vertex), fe_index, i, index); + const unsigned int fe_index = + (this->dof_handler->hp_capability_enabled == false && + fe_index_ == DoFHandlerType::invalid_fe_index) ? + DoFHandlerType::default_fe_index : + fe_index_; + + dealii::internal::DoFAccessorImplementation::Implementation::set_dof_index( + *this->dof_handler, + 0, + this->vertex_index(vertex), + fe_index, + i, + std::integral_constant(), + index); } @@ -2004,16 +1257,23 @@ DoFAccessor:: const unsigned int vertex, const unsigned int i, const types::global_dof_index index, - const unsigned int fe_index) const + const unsigned int fe_index_) const { + const unsigned int fe_index = + (this->dof_handler->hp_capability_enabled == false && + fe_index_ == DoFHandlerType::invalid_fe_index) ? + DoFHandlerType::default_fe_index : + fe_index_; (void)fe_index; Assert(this->dof_handler != nullptr, ExcInvalidObject()); AssertIndexRange(vertex, GeometryInfo::vertices_per_cell); AssertIndexRange(i, this->dof_handler->get_fe(fe_index).dofs_per_vertex); - return dealii::internal::DoFAccessorImplementation::Implementation:: - set_mg_vertex_dof_index( - *this->dof_handler, level, this->vertex_index(vertex), i, index); + Assert(dof_handler->hp_capability_enabled == false, + ExcMessage("hp::DoFHandler does not implement multilevel DoFs.")); + + this->dof_handler->mg_vertex_dofs[this->vertex_index(vertex)].set_index( + level, i, this->dof_handler->get_fe().dofs_per_vertex, index); } @@ -2162,11 +1422,18 @@ template inline void DoFAccessor::get_dof_indices( std::vector &dof_indices, - const unsigned int fe_index) const + const unsigned int fe_index_) const { Assert(this->dof_handler != nullptr, ExcInvalidObject()); + + const unsigned int fe_index = + (this->dof_handler->hp_capability_enabled == false && + fe_index_ == DoFHandlerType::invalid_fe_index) ? + DoFHandlerType::default_fe_index : + fe_index_; + Assert(static_cast(this->level()) < - this->dof_handler->levels.size(), + this->dof_handler->object_dof_indices.size(), ExcMessage( "The DoFHandler to which this accessor points has not " "been initialized, i.e., it doesn't appear that DoF indices " @@ -2235,10 +1502,16 @@ inline void DoFAccessor::get_mg_dof_indices( const int level, std::vector &dof_indices, - const unsigned int fe_index) const + const unsigned int fe_index_) const { Assert(this->dof_handler != nullptr, ExcInvalidObject()); + const unsigned int fe_index = + (this->dof_handler->hp_capability_enabled == false && + fe_index_ == DoFHandlerType::invalid_fe_index) ? + DoFHandlerType::default_fe_index : + fe_index_; + switch (structdim) { case 1: @@ -2293,10 +1566,16 @@ inline void DoFAccessor::set_mg_dof_indices( const int level, const std::vector &dof_indices, - const unsigned int fe_index) + const unsigned int fe_index_) { Assert(this->dof_handler != nullptr, ExcInvalidObject()); + const unsigned int fe_index = + (this->dof_handler->hp_capability_enabled == false && + fe_index_ == DoFHandlerType::invalid_fe_index) ? + DoFHandlerType::default_fe_index : + fe_index_; + switch (structdim) { case 1: @@ -2548,14 +1827,22 @@ template