From ed6822c9fd87f82835dac1819e5219b43523fcaf Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Thu, 30 Jul 2020 09:20:11 +0200 Subject: [PATCH] Non local dofs within FE. --- include/deal.II/dofs/dof_accessor.templates.h | 24 +++-- include/deal.II/fe/fe.h | 93 +++++++++++++++++-- source/dofs/dof_handler.cc | 22 +++-- tests/non_local/fe_q1_nonlocal.h | 75 +++++++-------- tests/non_local/fe_q1_nonlocal_01.cc | 3 +- tests/non_local/fe_q1_nonlocal_02.cc | 3 +- tests/non_local/fe_q1_nonlocal_03.cc | 4 +- tests/non_local/fe_q1_nonlocal_04.cc | 4 +- 8 files changed, 157 insertions(+), 71 deletions(-) diff --git a/include/deal.II/dofs/dof_accessor.templates.h b/include/deal.II/dofs/dof_accessor.templates.h index 08a9006f02..bef69bba72 100644 --- a/include/deal.II/dofs/dof_accessor.templates.h +++ b/include/deal.II/dofs/dof_accessor.templates.h @@ -822,11 +822,12 @@ namespace internal { const auto &fe = accessor.get_fe(fe_index_); - const unsigned int // - dofs_per_vertex = fe.n_dofs_per_vertex(), // - dofs_per_line = fe.n_dofs_per_line(), // - dofs_per_quad = fe.n_dofs_per_quad(), // - dofs_per_hex = fe.n_dofs_per_hex(); // + const unsigned int // + dofs_per_vertex = fe.n_dofs_per_vertex(), // + dofs_per_line = fe.n_dofs_per_line(), // + dofs_per_quad = fe.n_dofs_per_quad(), // + dofs_per_hex = fe.n_dofs_per_hex(), // + non_local_dofs = fe.n_non_local_dofs_per_cell(); // const unsigned int inner_dofs = structdim == 1 ? dofs_per_line : @@ -860,6 +861,10 @@ namespace internal // 4) INNER dofs index += inner_dofs; + // 5) non local dofs + if (dim == structdim) + index += non_local_dofs; + return index; } else @@ -874,6 +879,9 @@ namespace internal const auto diff = [](const auto &p) { return p.second - p.first; }; + const unsigned int non_local_dofs = + accessor.get_fe(fe_index).n_non_local_dofs_per_cell(); + // 1) VERTEX dofs for (const auto vertex : accessor.vertex_indices()) index += @@ -898,6 +906,10 @@ namespace internal // 4) INNER dofs index += diff(process_object_range(accessor, fe_index)); + // 5) non local dofs + if (dim == structdim) + index += non_local_dofs; + return index; } } @@ -926,7 +938,7 @@ namespace internal AssertDimension(dof_indices.size(), n_dof_indices(accessor, fe_index)); const auto &fe = accessor.get_fe(fe_index); - const auto non_local_dofs_per_cell = fe.n_non_local_dofs(); + const auto non_local_dofs_per_cell = fe.n_non_local_dofs_per_cell(); unsigned int index = 0; diff --git a/include/deal.II/fe/fe.h b/include/deal.II/fe/fe.h index 5ccb1589f9..7882b2d05e 100644 --- a/include/deal.II/fe/fe.h +++ b/include/deal.II/fe/fe.h @@ -44,6 +44,9 @@ template class FESubfaceValues; template class FESystem; +template +class DoFCellAccessor; + /** * This is the base class for finite elements in arbitrary dimensions. It @@ -2267,10 +2270,50 @@ public: * @{ */ /** - * Return an object that knows how to handle non local dof indices. + * Return the non local dof indices associated to the current cell, for + * active cell accessors. + */ + virtual std::vector + get_non_local_dof_indices( + const DoFCellAccessor &accessor) const; + + /** + * Return the non local dof indices associated to the current cell, for + * level cell accessors. + */ + virtual std::vector + get_non_local_dof_indices( + const DoFCellAccessor &accessor) const; + + /** + * Return the global number of non local dof indices that are required in + * addition to the local ones. + */ + virtual types::global_dof_index + n_global_non_local_dofs() const; + + /** + * Return an identification string that uniquely identifies the non local + * behaviour of the finite element space. + * + * For non local finite element spaces, n_non_local_dofs_per_cell() and + * n_global_non_local_dofs() may return a non zero number, meaning that there + * are degrees of freedom that are not associated to vertices, edges, faces, + * or cells (for example, they may be associated to patches of cells, or be + * global non zero basis functions), that are non zero on certain cells. + * + * In these cases, one usually uses the hp support of the library, and defines + * an FECollection where each FiniteElement of the collection has the same non + * local behaviour. This id is checked when calling + * DoFHandler::distribute_dofs() with an FECollection as argument, and an + * assertion is thrown if the ids do not coincide for all FiniteElement + * objects of the collection. + * + * By default, FiniteElement spaces are local, and this function returns the + * string "Local FiniteElement space". */ - virtual std::shared_ptr> - get_non_local_dof_handler() const; + virtual std::string + get_non_local_id() const; //@} /** @@ -3314,14 +3357,52 @@ FiniteElement::get_associated_geometry_primitive( +// Non local dofs support. +template +inline std::vector +FiniteElement::get_non_local_dof_indices( + const DoFCellAccessor &) const +{ + Assert(this->n_non_local_dofs_per_cell() == 0 && + n_global_non_local_dofs() == 0, + ExcPureFunctionCalled()); + return std::vector(); +} + + + +template +inline std::vector +FiniteElement::get_non_local_dof_indices( + const DoFCellAccessor &) const +{ + Assert(this->n_non_local_dofs_per_cell() == 0 && + n_global_non_local_dofs() == 0, + ExcPureFunctionCalled()); + return std::vector(); +} + + + template -inline std::shared_ptr> -FiniteElement::get_non_local_dof_handler() const +inline types::global_dof_index +FiniteElement::n_global_non_local_dofs() const { - return std::make_shared>(); + return 0; } + +template +inline std::string +FiniteElement::get_non_local_id() const +{ + Assert(this->n_non_local_dofs_per_cell() == 0 && + n_global_non_local_dofs() == 0, + ExcPureFunctionCalled()); + return "Local FiniteElement space"; +} + DEAL_II_NAMESPACE_CLOSE #endif diff --git a/source/dofs/dof_handler.cc b/source/dofs/dof_handler.cc index 1485977780..9f5c6a07f4 100644 --- a/source/dofs/dof_handler.cc +++ b/source/dofs/dof_handler.cc @@ -2612,18 +2612,22 @@ DoFHandler::distribute_non_local_dofs() ExcMessage( "Distribute active DoFs using distribute_dofs() before calling distribute_non_local_dofs().")); - auto non_local_dh = begin()->get_fe().get_non_local_dof_handler(); + const auto non_local_id = begin()->get_fe().get_non_local_id(); + const auto n_global_non_local_dofs = + begin()->get_fe().n_global_non_local_dofs(); + for (auto cell : active_cell_iterators()) { - // Make sure all cells use the same NonLocalDoFHandler. - Assert(non_local_dh == cell->get_fe().get_non_local_dof_handler(), - ExcMessage( - "You are trying to use different NonLocalDoFHandler objects")); - cell->set_non_local_dof_indices( - non_local_dh->get_non_local_dof_indices(*cell)); + const auto &fe = cell->get_fe(); + // Make sure all cells handle non local dofs in the same way. + Assert( + non_local_id == fe.get_non_local_id(), + ExcMessage( + "You are trying to use different non local finite element spaces")); + AssertDimension(n_global_non_local_dofs, fe.n_global_non_local_dofs()); + cell->set_non_local_dof_indices(fe.get_non_local_dof_indices(*cell)); } - this->number_cache.n_global_dofs += - non_local_dh->n_additional_non_local_dofs(); + this->number_cache.n_global_dofs += n_global_non_local_dofs; } diff --git a/tests/non_local/fe_q1_nonlocal.h b/tests/non_local/fe_q1_nonlocal.h index 1d713624ac..f9a7c1c111 100644 --- a/tests/non_local/fe_q1_nonlocal.h +++ b/tests/non_local/fe_q1_nonlocal.h @@ -21,6 +21,8 @@ #include #include +#include + #include #include @@ -36,45 +38,11 @@ get_dpo_vector() return dpo; } - -/** - * Associates non local dofs according to the vertex index. - */ -template -class NonLocalQ1DoFHandler : public NonLocalDoFHandler -{ -public: - virtual std::vector - get_non_local_dof_indices( - const DoFCellAccessor &accessor) const override - { - if (!tria) - tria = &(accessor.get_triangulation()); - std::vector dofs(accessor.n_vertices()); - for (unsigned int i = 0; i < dofs.size(); ++i) - dofs[i] = accessor.vertex_index(i); - return dofs; - } - - virtual types::global_dof_index - n_additional_non_local_dofs() const override - { - Assert(tria, ExcInternalError()); - return tria->n_vertices(); - } - -private: - mutable SmartPointer> tria = nullptr; -}; - - - template class FE_Q1_Nonlocal : public FE_Q_Base, dim, dim> { public: - FE_Q1_Nonlocal(const std::shared_ptr> &ptr = - std::shared_ptr>()) + FE_Q1_Nonlocal(const Triangulation &tria) : FE_Q_Base, dim, dim>( TensorProductPolynomials( Polynomials::generate_complete_Lagrange_basis( @@ -84,7 +52,7 @@ public: 1, FiniteElementData::H1), std::vector(1, false)) - , non_local_dh(ptr ? ptr : std::make_shared>()) + , tria(&tria) { this->unit_support_points = QTrapez().get_points(); } @@ -92,7 +60,7 @@ public: virtual std::unique_ptr> clone() const override { - return std::make_unique>(non_local_dh); + return std::make_unique>(*tria); } virtual std::string @@ -101,12 +69,6 @@ public: return "FE_Q_Nonlocal"; } - virtual std::shared_ptr> - get_non_local_dof_handler() const override - { - return non_local_dh; - } - virtual void convert_generalized_support_point_values_to_dof_values( const std::vector> &support_point_values, @@ -124,6 +86,31 @@ public: } } + + virtual std::vector + get_non_local_dof_indices( + const DoFCellAccessor &accessor) const override + { + std::vector dofs(accessor.n_vertices()); + for (unsigned int i = 0; i < dofs.size(); ++i) + dofs[i] = accessor.vertex_index(i); + return dofs; + } + + + virtual types::global_dof_index + n_global_non_local_dofs() const override + { + Assert(tria, ExcInternalError()); + return tria->n_vertices(); + } + + virtual std::string + get_non_local_id() const override + { + return "NonLocal FEQ1"; + } + private: - std::shared_ptr> non_local_dh; + mutable SmartPointer> tria = nullptr; }; diff --git a/tests/non_local/fe_q1_nonlocal_01.cc b/tests/non_local/fe_q1_nonlocal_01.cc index 5633da7015..56f2c73c46 100644 --- a/tests/non_local/fe_q1_nonlocal_01.cc +++ b/tests/non_local/fe_q1_nonlocal_01.cc @@ -31,7 +31,8 @@ template void test() { - FE_Q1_Nonlocal fe; + Triangulation tria; + FE_Q1_Nonlocal fe(tria); if (fe.n_dofs_per_cell() != GeometryInfo::vertices_per_cell) deallog << "Not OK" << std::endl; else diff --git a/tests/non_local/fe_q1_nonlocal_02.cc b/tests/non_local/fe_q1_nonlocal_02.cc index a490caa192..0339d61642 100644 --- a/tests/non_local/fe_q1_nonlocal_02.cc +++ b/tests/non_local/fe_q1_nonlocal_02.cc @@ -31,7 +31,8 @@ void test() { deallog << "DIM: " << dim << std::endl << std::endl; - FE_Q1_Nonlocal fe; + Triangulation tria; + FE_Q1_Nonlocal fe(tria); // Check that the basis functions are what we expect const auto &quad = fe.get_unit_support_points(); diff --git a/tests/non_local/fe_q1_nonlocal_03.cc b/tests/non_local/fe_q1_nonlocal_03.cc index 5746bba088..8a06b270b2 100644 --- a/tests/non_local/fe_q1_nonlocal_03.cc +++ b/tests/non_local/fe_q1_nonlocal_03.cc @@ -31,9 +31,9 @@ template void test() { - FE_Q1_Nonlocal fe; - Triangulation tria; + Triangulation tria; GridGenerator::hyper_cube(tria); + FE_Q1_Nonlocal fe(tria); tria.refine_global(1); DoFHandler dh(tria); dh.distribute_dofs(fe); diff --git a/tests/non_local/fe_q1_nonlocal_04.cc b/tests/non_local/fe_q1_nonlocal_04.cc index cebe016ec8..5150249355 100644 --- a/tests/non_local/fe_q1_nonlocal_04.cc +++ b/tests/non_local/fe_q1_nonlocal_04.cc @@ -35,9 +35,9 @@ template void test() { - FE_Q1_Nonlocal fe; - Triangulation tria; + Triangulation tria; GridGenerator::hyper_cube(tria); + FE_Q1_Nonlocal fe(tria); tria.refine_global(1); DoFHandler dh(tria); dh.distribute_dofs(fe); -- 2.39.5