From: Luca Heltai Date: Thu, 23 Jul 2020 13:15:05 +0000 (+0200) Subject: NonLocalDoFHandler. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b7c91bb6603dc881b34cffbeb7f52089effee616;p=dealii.git NonLocalDoFHandler. --- diff --git a/include/deal.II/dofs/dof_accessor.templates.h b/include/deal.II/dofs/dof_accessor.templates.h index 9f2d6dd98c..08a9006f02 100644 --- a/include/deal.II/dofs/dof_accessor.templates.h +++ b/include/deal.II/dofs/dof_accessor.templates.h @@ -2252,88 +2252,6 @@ namespace internal - /** - * Implement setting dof indices on a cell. TO CHECK ZHAOWEI + LH - */ - template - static void - set_dof_indices( - const DoFCellAccessor &accessor, - const std::vector &local_dof_indices) - { - Assert(accessor.has_children() == false, ExcInternalError()); - - const unsigned int dofs_per_vertex = - accessor.get_fe().n_dofs_per_vertex(), - dofs_per_line = accessor.get_fe().n_dofs_per_line(), - dofs_per_quad = accessor.get_fe().n_dofs_per_quad(), - dofs_per_hex = accessor.get_fe().n_dofs_per_hex(); - - Assert(local_dof_indices.size() == accessor.get_fe().dofs_per_cell, - ExcInternalError()); - - unsigned int index = 0; - - for (unsigned int vertex = 0; - vertex < GeometryInfo::vertices_per_cell; - ++vertex) - for (unsigned int d = 0; d < dofs_per_vertex; ++d, ++index) - accessor.set_vertex_dof_index(vertex, - d, - local_dof_indices[index], - accessor.active_fe_index()); - // now copy dof numbers into the line. for lines in 3d with the - // wrong orientation, we have already made sure that we're ok - // by picking the correct vertices (this happens automatically - // in the vertex() function). however, if the line is in wrong - // orientation, we look at it in flipped orientation and we - // will have to adjust the shape function indices that we see - // to correspond to the correct (cell-local) ordering. - // - // of course, if dim<3, then there is nothing to adjust - for (unsigned int line = 0; line < GeometryInfo::lines_per_cell; - ++line) - for (unsigned int d = 0; d < dofs_per_line; ++d, ++index) - accessor.line(line)->set_dof_index( - dim < 3 ? - d : - accessor.get_fe().adjust_line_dof_index_for_line_orientation( - d, accessor.line_orientation(line)), - local_dof_indices[index], - accessor.active_fe_index()); - // now copy dof numbers into the face. for faces in 3d with the - // wrong orientation, we have already made sure that we're ok - // by picking the correct lines and vertices (this happens - // automatically in the line() and vertex() - // functions). however, if the face is in wrong orientation, - // we look at it in flipped orientation and we will have to - // adjust the shape function indices that we see to correspond - // to the correct (cell-local) ordering. The same applies, if - // the face_rotation or face_orientation is non-standard - // - // again, if dim<3, then there is nothing to adjust - for (unsigned int quad = 0; quad < GeometryInfo::quads_per_cell; - ++quad) - for (unsigned int d = 0; d < dofs_per_quad; ++d, ++index) - accessor.quad(quad)->set_dof_index( - dim < 3 ? - d : - accessor.get_fe().adjust_quad_dof_index_for_face_orientation( - d, - accessor.face_orientation(quad), - accessor.face_flip(quad), - accessor.face_rotation(quad)), - local_dof_indices[index], - accessor.active_fe_index()); - for (unsigned int d = 0; d < dofs_per_hex; ++d, ++index) - accessor.set_dof_index(d, - local_dof_indices[index], - accessor.active_fe_index()); - Assert(index == accessor.get_fe().dofs_per_cell, ExcInternalError()); - } - - - /** * Implement setting non-local dof indices on a cell. */ @@ -2345,7 +2263,16 @@ namespace internal { Assert(accessor.has_children() == false, ExcInternalError()); - const unsigned int dofs_per_cell = accessor.get_fe().dofs_per_cell; + const unsigned int dofs_per_cell = accessor.get_fe().n_dofs_per_cell(); + const unsigned int n_non_local_dofs = + accessor.get_fe().n_non_local_dofs_per_cell(); + const unsigned int n_local_dofs = dofs_per_cell - n_non_local_dofs; + + AssertDimension(local_non_local_dof_indices.size(), n_non_local_dofs); + + // First the easy case. + if (n_non_local_dofs == 0) + return; types::global_dof_index *next_dof_index = const_cast( @@ -2353,22 +2280,11 @@ namespace internal get_cache_ptr(accessor.dof_handler, accessor.present_level, accessor.present_index, - dofs_per_cell)); - - const unsigned int non_local_dofs = - accessor.get_fe().n_non_local_dofs_per_cell(), - n_dofs = accessor.get_fe().n_dofs_per_cell(); - - unsigned int index = 0; + dofs_per_cell)) + + n_local_dofs; - for (unsigned int d = 0; d < n_dofs; ++d, ++next_dof_index) - if (d >= n_dofs - non_local_dofs) - { - *next_dof_index = local_non_local_dof_indices[index]; - ++index; - } - Assert(index == accessor.get_fe().n_non_local_dofs_per_cell(), - ExcInternalError()); + for (unsigned int d = 0; d < n_non_local_dofs; ++d, ++next_dof_index) + *next_dof_index = local_non_local_dof_indices[d]; } diff --git a/include/deal.II/dofs/dof_handler.h b/include/deal.II/dofs/dof_handler.h index f727ef1e27..a0d8a2f79c 100644 --- a/include/deal.II/dofs/dof_handler.h +++ b/include/deal.II/dofs/dof_handler.h @@ -592,6 +592,13 @@ public: void distribute_mg_dofs(); + /** + * Distribute non local degrees of freedom. The local DoFs need to be + * distributed using distribute_dofs() before calling this function. + */ + void + distribute_non_local_dofs(); + /** * This function returns whether this DoFHandler has DoFs distributed on * each multigrid level or in other words if distribute_mg_dofs() has been diff --git a/include/deal.II/dofs/non_local_dof_handler.h b/include/deal.II/dofs/non_local_dof_handler.h new file mode 100644 index 0000000000..7cfe068e5d --- /dev/null +++ b/include/deal.II/dofs/non_local_dof_handler.h @@ -0,0 +1,149 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii_non_local_dof_handler_h +#define dealii_non_local_dof_handler_h + + +#include + +#include + +DEAL_II_NAMESPACE_OPEN + +// Forward declarations +#ifndef DOXYGEN +template +class DoFCellAccessor; +#endif + +/** + * This class is used to enumerate non local degrees of freedom. Its default + * implementation does nothing, since in general FiniteElement spaces only + * define degrees of freedom on vertices, edges, faces, or cells. There are + * cases, however, in which a FiniteElement space may define some non-local + * basis functions which are non-zero on a given cell, even if the basis + * function cannot be attributed to local subobjects of the given cell. + * + * In these cases, the DoFHandler class needs to know how to distribute these + * extra degrees of freedom, and since it cannot do so on its own, it asks + * the FiniteElement to return a NonLocalDoFHandler class that knows how to + * enumerate and distribute these non local degrees of freedom. + * + * This class is intended as a base class for all FiniteElement spaces that + * need to enumerate non local degrees of freedom. + * + * @ingroup dofs + */ +template +class NonLocalDoFHandler : public Subscriptor +{ +public: + /** + * Make the dimension available in function templates. + */ + static const unsigned int dimension = dim; + + /** + * Make the space dimension available in function templates. + */ + static const unsigned int space_dimension = spacedim; + + /** + * Standard constructor, not initializing any data. + */ + NonLocalDoFHandler() = default; + + /** + * Destructor. + */ + virtual ~NonLocalDoFHandler() = default; + + /** + * Copy operator. NonLocalDoFHandler objects may be large and expensive. + * They should not be copied, in particular not by accident, but + * rather deliberately constructed. As a consequence, this operator + * is explicitly removed from the interface of this class. + */ + NonLocalDoFHandler & + operator=(const NonLocalDoFHandler &) = delete; + + /** + * Return the non local dof indices associated to the current cell, for + * active cell iterators. + */ + virtual std::vector + get_non_local_dof_indices( + const DoFCellAccessor &) const; + + /** + * Return the non local dof indices associated to the current cell, for + * level cell iterators. + */ + virtual std::vector + get_non_local_dof_indices(const DoFCellAccessor &) const; + + /** + * Return the number of non local dof indices that are required in addition to + * the local ones. + */ + virtual types::global_dof_index + n_additional_non_local_dofs() const; + + /** + * At any given moment, only one NonLocalDoFHandler can be used. Throw an + * exception if two cells define an active FiniteElement for which + * get_non_local_dof_handler() return a different object. + * + * @ingroup Exceptions + */ + DeclException0(ExcDifferentNonLocalDoFHandlers); +}; + + + +// ---------------------------------------------------------------------------- +template +inline std::vector +NonLocalDoFHandler::get_non_local_dof_indices( + const DoFCellAccessor &) const +{ + // By default, we return an empty vector. + return std::vector(); +} + + + +template +inline std::vector +NonLocalDoFHandler::get_non_local_dof_indices( + const DoFCellAccessor &) const +{ + // By default, we return an empty vector. + return std::vector(); +} + + + +template +inline types::global_dof_index +NonLocalDoFHandler::n_additional_non_local_dofs() const +{ + return 0; +} + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/include/deal.II/fe/fe.h b/include/deal.II/fe/fe.h index 586a3784ff..5ccb1589f9 100644 --- a/include/deal.II/fe/fe.h +++ b/include/deal.II/fe/fe.h @@ -18,6 +18,8 @@ #include +#include + #include #include #include @@ -2260,6 +2262,17 @@ public: //@} + /** + * @name Non local dofs support + * @{ + */ + /** + * Return an object that knows how to handle non local dof indices. + */ + virtual std::shared_ptr> + get_non_local_dof_handler() const; + //@} + /** * Determine an estimate for the memory consumption (in bytes) of this * object. @@ -3301,6 +3314,14 @@ FiniteElement::get_associated_geometry_primitive( +template +inline std::shared_ptr> +FiniteElement::get_non_local_dof_handler() const +{ + return std::make_shared>(); +} + + DEAL_II_NAMESPACE_CLOSE #endif diff --git a/source/dofs/dof_handler.cc b/source/dofs/dof_handler.cc index 70a4c40cb1..1485977780 100644 --- a/source/dofs/dof_handler.cc +++ b/source/dofs/dof_handler.cc @@ -2565,6 +2565,8 @@ DoFHandler::distribute_dofs( &*this->tria) == nullptr) this->block_info_object.initialize(*this, false, true); } + // now distribute non_local_dofs. + distribute_non_local_dofs(); } @@ -2601,6 +2603,31 @@ DoFHandler::distribute_mg_dofs() +template +void +DoFHandler::distribute_non_local_dofs() +{ + Assert( + this->object_dof_indices.size() > 0, + 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(); + 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)); + } + this->number_cache.n_global_dofs += + non_local_dh->n_additional_non_local_dofs(); +} + + + template void DoFHandler::initialize_local_block_info() diff --git a/tests/non_local/fe_q1_nonlocal.h b/tests/non_local/fe_q1_nonlocal.h index 0e1a518af1..1d713624ac 100644 --- a/tests/non_local/fe_q1_nonlocal.h +++ b/tests/non_local/fe_q1_nonlocal.h @@ -36,11 +36,45 @@ 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() + FE_Q1_Nonlocal(const std::shared_ptr> &ptr = + std::shared_ptr>()) : FE_Q_Base, dim, dim>( TensorProductPolynomials( Polynomials::generate_complete_Lagrange_basis( @@ -50,6 +84,7 @@ public: 1, FiniteElementData::H1), std::vector(1, false)) + , non_local_dh(ptr ? ptr : std::make_shared>()) { this->unit_support_points = QTrapez().get_points(); } @@ -57,7 +92,7 @@ public: virtual std::unique_ptr> clone() const override { - return std::make_unique>(); + return std::make_unique>(non_local_dh); } virtual std::string @@ -65,4 +100,30 @@ 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, + std::vector & nodal_values) const override + { + AssertDimension(support_point_values.size(), + this->get_unit_support_points().size()); + AssertDimension(support_point_values.size(), nodal_values.size()); + AssertDimension(this->n_dofs_per_cell(), nodal_values.size()); + + for (unsigned int i = 0; i < this->n_dofs_per_cell(); ++i) + { + AssertDimension(support_point_values[i].size(), 1); + nodal_values[i] = support_point_values[i](0); + } + } + +private: + std::shared_ptr> non_local_dh; }; diff --git a/tests/non_local/fe_q1_nonlocal_03.output b/tests/non_local/fe_q1_nonlocal_03.output index e69de29bb2..fb71de2867 100644 --- a/tests/non_local/fe_q1_nonlocal_03.output +++ b/tests/non_local/fe_q1_nonlocal_03.output @@ -0,0 +1,4 @@ + +DEAL::OK +DEAL::OK +DEAL::OK diff --git a/tests/non_local/fe_q1_nonlocal_04.cc b/tests/non_local/fe_q1_nonlocal_04.cc new file mode 100644 index 0000000000..cebe016ec8 --- /dev/null +++ b/tests/non_local/fe_q1_nonlocal_04.cc @@ -0,0 +1,64 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 2018 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// build an FE_Q1_Nonlocal finite element, distribute dofs on a simple +// Triangulation, and interpolate a scalar function. + +#include + +#include + +#include +#include + +#include + +#include + +#include "../tests.h" + +#include "../non_local/fe_q1_nonlocal.h" + +template +void +test() +{ + FE_Q1_Nonlocal fe; + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(1); + DoFHandler dh(tria); + dh.distribute_dofs(fe); + + Tensor<1, dim> ones; + for (unsigned int d = 0; d < dim; ++d) + ones[d] = 1; + + Vector solution(dh.n_dofs()); + + VectorTools::interpolate(dh, Functions::Monomial(ones), solution); + + deallog << solution << std::endl; +} + +int +main() +{ + initlog(); + + test<1>(); + test<2>(); + test<3>(); +} diff --git a/tests/non_local/fe_q1_nonlocal_04.output b/tests/non_local/fe_q1_nonlocal_04.output new file mode 100644 index 0000000000..e25d0a7fbc --- /dev/null +++ b/tests/non_local/fe_q1_nonlocal_04.output @@ -0,0 +1,4 @@ + +DEAL::0.00000 1.00000 0.500000 +DEAL::0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.500000 0.500000 0.250000 +DEAL::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.500000 0.00000 0.00000 0.500000 0.500000 0.00000 0.00000 0.00000 0.250000 0.250000 0.250000 0.125000