From be4fcc0a2569e74269124835fdff6e14c3464d1f Mon Sep 17 00:00:00 2001 From: Marc Fehling Date: Tue, 26 Feb 2019 16:34:09 +0100 Subject: [PATCH] Added hierarchy to hp::FECollection. --- doc/doxygen/headers/hp.h | 15 ++++ doc/news/changes/minor/20190228Fehling | 7 ++ include/deal.II/hp/fe_collection.h | 118 ++++++++++++++++++++++++- source/hp/fe_collection.cc | 68 ++++++++++++++ tests/hp/fe_hierarchy.cc | 61 +++++++++++++ tests/hp/fe_hierarchy.output | 11 +++ 6 files changed, 277 insertions(+), 3 deletions(-) create mode 100644 doc/news/changes/minor/20190228Fehling create mode 100644 tests/hp/fe_hierarchy.cc create mode 100644 tests/hp/fe_hierarchy.output diff --git a/doc/doxygen/headers/hp.h b/doc/doxygen/headers/hp.h index 5fceed9e24..b9c6953168 100644 --- a/doc/doxygen/headers/hp.h +++ b/doc/doxygen/headers/hp.h @@ -89,6 +89,21 @@ * collection of mappings is used, the same holds for hp::MappingCollection * objects as well. * + * Whenever p adaptivity is considered in an hp finite element program, + * a hierarchy of finite elements needs to be established to determine + * succeeding finite elements for refinement and preceding ones for coarsening. + * Typically, this hierarchy considers how finite element spaces are nested: + * for example, a $Q_1$ element describes a sub-space of a $Q_2$ element, + * and so doing $p$ refinement usually means using a larger (more accurate) + * finite element space. In other words, the hierarchy of finite elements is built + * by considering whether some elements of the collection are sub- or + * super-spaces of others. + * + * By default, we assume that finite elements are stored in an ascending order + * based on their polynomial degree. If the order of elements differs, + * a corresponding hierarchy needs to be supplied to the collection via the + * hp::FECollection::set_hierarchy() member function. + * * @ingroup hp */ diff --git a/doc/news/changes/minor/20190228Fehling b/doc/news/changes/minor/20190228Fehling new file mode 100644 index 0000000000..85c0089d88 --- /dev/null +++ b/doc/news/changes/minor/20190228Fehling @@ -0,0 +1,7 @@ +New: Hierarchy of finite elements in hp::FECollection objects. Get succeeding +and preceding indices via hp::FECollection::next_in_hierarchy() and +hp::FECollection::prev_in_hierarchy(). By default, a hierarchy corresponding +to indices is established. Hierarchy can be overridden via +hp::FECollection::set_hierarchy(). +
+(Marc Fehling, 2019/02/28) diff --git a/include/deal.II/hp/fe_collection.h b/include/deal.II/hp/fe_collection.h index a58bef7ea2..c4668b89dc 100644 --- a/include/deal.II/hp/fe_collection.h +++ b/include/deal.II/hp/fe_collection.h @@ -54,11 +54,53 @@ namespace hp class FECollection : public Subscriptor { public: + /** + * Whenever p adaptivity is considered in an hp finite element program, + * a hierarchy of finite elements needs to be established to determine + * succeeding finite elements for refinement and preceding ones for + * coarsening. + * + * In this struct, we supply a hierachy that is imposed on all FECollection + * objects by default. + */ + struct DefaultHierarchy + { + /** + * Return the index succeeding @p fe_index in the @p fe_collection. + * + * Once the last element of the @p fe_collection is reached, there is no element on a higher level in + * the hierarchy and thus we return the last value. + */ + static unsigned int + next_index(const typename hp::FECollection &fe_collection, + const unsigned int fe_index) + { + return ((fe_index + 1) < fe_collection.size()) ? fe_index + 1 : + fe_index; + } + + /** + * Return the index preceding @p fe_index in the @p fe_collection. + * + * Once the first element of the @p fe_collection is reached, there is no element on a lower level in + * the hierarchy and thus we return the first value. + */ + static unsigned int + previous_index( + const typename hp::FECollection &fe_collection, + const unsigned int fe_index) + { + (void)fe_collection; + return (fe_index > 0) ? fe_index - 1 : fe_index; + } + }; + /** * Default constructor. Leads to an empty collection that can later be - * filled using push_back(). + * filled using push_back(). Establishes a hierarchy of finite elements + * corresponding to their index in the collection. */ - FECollection() = default; + FECollection(); /** * Conversion constructor. This constructor creates a FECollection from a @@ -94,7 +136,7 @@ namespace hp /** * Move constructor. */ - FECollection(FECollection &&) noexcept = default; + FECollection(FECollection &&) = default; /** * Move assignment operator. @@ -371,6 +413,60 @@ namespace hp find_dominating_fe_in_subset(const std::set &fes, const unsigned int codim = 0) const; + /** + * Set functions determining the hierarchy of finite elements, i.e. a + * function @p next that returns the index of the finite element following + * the given one, and a function @p prev returning the preceding one. + * + * Both functions expect an hp::FECollection to be passed along with a + * finite element index, on whose basis the new index will be found and + * returned. + * + * @note Both passed and returned indices have to be valid within the index + * range of this collection, i.e. within [0, size()). + */ + void + set_hierarchy(const std::function &, + const unsigned int)> &next, + const std::function &, + const unsigned int)> &prev); + + /** + * Set the default hierarchy corresponding to the index of each finite + * element in the collection. + * + * This default hierarchy is established with functions + * DefaultHierarchy::next_index() and DefaultHierarchy::previous_index(). + */ + void + set_default_hierarchy(); + + /** + * Function returning the index of the finite element following the given + * @p fe_index in hierarchy. + * + * By default, the index succeeding @p fe_index will be returned. If @p fe_index + * already corresponds to the last index, the last index will be returned. + * A custom hierarchy can be supplied via the member function + * set_hierachy(). + */ + unsigned int + next_in_hierarchy(const unsigned int fe_index) const; + + /** + * Function returning the index of the finite element preceding the given + * @p fe_index in hierarchy. + * + * By default, the index preceding @p fe_index will be returned. If @p fe_index + * already corresponds to the first index, the first index will be returned. + * A custom hierarchy can be supplied via the member function + * set_hierachy(). + */ + unsigned int + previous_in_hierarchy(const unsigned int fe_index) const; + /** * Return a component mask with as many elements as this object has vector * components and of which exactly the one component is true that @@ -576,6 +672,22 @@ namespace hp */ std::vector>> finite_elements; + + /** + * Function returning the index of the finite element following the given + * one in hierarchy. + */ + std::function &, + const unsigned int)> + hierarchy_next; + + /** + * Function returning the index of the finite element preceding the given + * one in hierarchy. + */ + std::function &, + const unsigned int)> + hierarchy_prev; }; diff --git a/source/hp/fe_collection.cc b/source/hp/fe_collection.cc index 616359767d..cc3c26634a 100644 --- a/source/hp/fe_collection.cc +++ b/source/hp/fe_collection.cc @@ -160,9 +160,18 @@ namespace hp + template + FECollection::FECollection() + { + set_default_hierarchy(); + } + + + template FECollection::FECollection( const FiniteElement &fe) + : FECollection() { push_back(fe); } @@ -172,6 +181,7 @@ namespace hp template FECollection::FECollection( const std::vector *> &fes) + : FECollection() { Assert(fes.size() > 0, ExcMessage("Need to pass at least one finite element.")); @@ -202,6 +212,64 @@ namespace hp + template + void + FECollection::set_hierarchy( + const std::function< + unsigned int(const typename hp::FECollection &, + const unsigned int)> &next, + const std::function< + unsigned int(const typename hp::FECollection &, + const unsigned int)> &prev) + { + // copy hierarchy functions + hierarchy_next = next; + hierarchy_prev = prev; + } + + + + template + void + FECollection::set_default_hierarchy() + { + // establish hierarchy corresponding to order of indices + set_hierarchy(&DefaultHierarchy::next_index, + &DefaultHierarchy::previous_index); + } + + + + template + unsigned int + FECollection::next_in_hierarchy( + const unsigned int fe_index) const + { + Assert(fe_index < size(), ExcIndexRange(fe_index, 0, size())); + + const unsigned int new_fe_index = hierarchy_next(*this, fe_index); + Assert(new_fe_index < size(), ExcIndexRange(new_fe_index, 0, size())); + + return new_fe_index; + } + + + + template + unsigned int + FECollection::previous_in_hierarchy( + const unsigned int fe_index) const + { + Assert(fe_index < size(), ExcIndexRange(fe_index, 0, size())); + + const unsigned int new_fe_index = hierarchy_prev(*this, fe_index); + Assert(new_fe_index < size(), ExcIndexRange(new_fe_index, 0, size())); + + return new_fe_index; + } + + + template ComponentMask FECollection::component_mask( diff --git a/tests/hp/fe_hierarchy.cc b/tests/hp/fe_hierarchy.cc new file mode 100644 index 0000000000..185421afa4 --- /dev/null +++ b/tests/hp/fe_hierarchy.cc @@ -0,0 +1,61 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2005 - 2019 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. +// +// --------------------------------------------------------------------- + + + +// Test default hierarchy of any hp::FECollection object. +// In the default implementation, the new indices are either the +// succeeding or preceding indices until limits are reached. + + +#include + +#include + +#include "../tests.h" + + +template +void +test() +{ + hp::FECollection fe_collection; + + while (fe_collection.size() < 3) + { + // add dummy fe to collection + fe_collection.push_back(FE_Q(1)); + deallog << "size:" << fe_collection.size() << std::endl; + + for (unsigned int fe_index = 0; fe_index < fe_collection.size(); + ++fe_index) + deallog << " idx:" << fe_index + << " next:" << fe_collection.next_in_hierarchy(fe_index) + << " prev:" << fe_collection.previous_in_hierarchy(fe_index) + << std::endl; + } +} + + + +int +main() +{ + initlog(); + + test<1>(); + + deallog << "OK" << std::endl; +} diff --git a/tests/hp/fe_hierarchy.output b/tests/hp/fe_hierarchy.output new file mode 100644 index 0000000000..64e52b61e7 --- /dev/null +++ b/tests/hp/fe_hierarchy.output @@ -0,0 +1,11 @@ + +DEAL::size:1 +DEAL:: idx:0 next:0 prev:0 +DEAL::size:2 +DEAL:: idx:0 next:1 prev:0 +DEAL:: idx:1 next:1 prev:0 +DEAL::size:3 +DEAL:: idx:0 next:1 prev:0 +DEAL:: idx:1 next:2 prev:0 +DEAL:: idx:2 next:2 prev:1 +DEAL::OK -- 2.39.5