From 566e08c975c95bdf9a65026898dde7da1bec9f25 Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Fri, 6 Oct 2017 14:46:06 -0400 Subject: [PATCH] add FiniteElement::get_sub_fe --- doc/news/changes/minor/20171007Heister | 3 + include/deal.II/fe/fe.h | 87 ++++++++++++ include/deal.II/fe/fe_system.h | 11 ++ source/fe/fe.cc | 51 +++++++ source/fe/fe_system.cc | 26 +++- tests/fe/get_sub_fe_01.cc | 180 +++++++++++++++++++++++++ tests/fe/get_sub_fe_01.debug.output | 53 ++++++++ tests/fe/get_sub_fe_02.cc | 81 +++++++++++ tests/fe/get_sub_fe_02.output | 49 +++++++ 9 files changed, 540 insertions(+), 1 deletion(-) create mode 100644 doc/news/changes/minor/20171007Heister create mode 100644 tests/fe/get_sub_fe_01.cc create mode 100644 tests/fe/get_sub_fe_01.debug.output create mode 100644 tests/fe/get_sub_fe_02.cc create mode 100644 tests/fe/get_sub_fe_02.output diff --git a/doc/news/changes/minor/20171007Heister b/doc/news/changes/minor/20171007Heister new file mode 100644 index 0000000000..fc8b210ded --- /dev/null +++ b/doc/news/changes/minor/20171007Heister @@ -0,0 +1,3 @@ +New: FiniteElement::get_sub_fe() allows you to return a contained FiniteElement based on a ComponentMask. +
+(Timo Heister, 2017/10/07) diff --git a/include/deal.II/fe/fe.h b/include/deal.II/fe/fe.h index e4ae7e555f..cbe84c9faa 100644 --- a/include/deal.II/fe/fe.h +++ b/include/deal.II/fe/fe.h @@ -1573,6 +1573,92 @@ public: unsigned int element_multiplicity (const unsigned int index) const; + /** + * Return a reference to a contained finite element that matches the components + * selected by the given ComponentMask @p mask. + * + * For an arbitrarily nested FESystem, this function returns the inner-most + * FiniteElement that matches the given mask. The method fails if the @p mask + * does not exactly match one of the contained finite elements. It is most + * useful if the current object is an FESystem, as the return value can + * only be @p this in all other cases. + * + * Note that the returned object can be an FESystem if the + * mask matches it but not any of the contained objects. + * + * Let us illustrate the function with the an FESystem @p fe with 7 components: + * @code + * FESystem<2> fe_velocity(FE_Q<2>(2), 2); + * FE_Q<2> fe_pressure(1); + * FE_DGP<2> fe_dg(0); + * FE_BDM<2> fe_nonprim(1); + * FESystem<2> fe(fe_velocity, 1, fe_pressure, 1, fe_dg, 2, fe_nonprim, 1); + * @endcode + * + * The following table lists all possible component masks you can use: + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + *
ComponentMaskResultDescription
[true,true,true,true,true,true,true]FESystem<2>[FESystem<2>[FE_Q<2>(2)^2]-FE_Q<2>(1)-FE_DGP<2>(0)^2-FE_BDM<2>(1)]@p fe itself, the whole @p FESystem
[true,true,false,false,false,false,false]FESystem<2>[FE_Q<2>(2)^2]just the @p fe_velocity
[true,false,false,false,false,false,false]FE_Q<2>(2)The first component in @p fe_velocity
[false,true,false,false,false,false,false]FE_Q<2>(2)The second component in @p fe_velocity
[false,false,true,false,false,false,false]FE_Q<2>(1)@p fe_pressure
[false,false,false,true,false,false,false]FE_DGP<2>(0)first copy of @p fe_dg
[false,false,false,false,true,false,false]FE_DGP<2>(0)second copy of @p fe_dg
[false,false,false,false,false,true,true]FE_BDM<2>(1)both components of @p fe_nonprim
+ */ + const FiniteElement & + get_sub_fe (const ComponentMask &mask) const; + + /** + * Return a reference to a contained finite element that matches the components + * @p n_selected_components components starting at component with index + * @p first_component. + * + * See the other get_sub_fe() function above for more details. + */ + virtual + const FiniteElement & + get_sub_fe (const unsigned int first_component, + const unsigned int n_selected_components) const; + /** * Return for shape function @p index the base element it belongs to, the * number of the copy of this base element (which is between zero and the @@ -2989,6 +3075,7 @@ FiniteElement::face_system_to_component_index (const unsigned int + template inline std::pair,unsigned int> diff --git a/include/deal.II/fe/fe_system.h b/include/deal.II/fe/fe_system.h index f9218330df..9801df566c 100644 --- a/include/deal.II/fe/fe_system.h +++ b/include/deal.II/fe/fe_system.h @@ -519,6 +519,17 @@ public: UpdateFlags requires_update_flags (const UpdateFlags update_flags) const; + // make variant with ComponentMask also available: + using FiniteElement::get_sub_fe; + + /** + * @copydoc FiniteElement::get_sub_fe() + */ + virtual + const FiniteElement & + get_sub_fe (const unsigned int first_component, + const unsigned int n_selected_components) const; + /** * Return the value of the @p ith shape function at the point @p p. @p p is * a point on the reference element. Since this finite element is always diff --git a/source/fe/fe.cc b/source/fe/fe.cc index 27538dbe65..4a5beb7e76 100644 --- a/source/fe/fe.cc +++ b/source/fe/fe.cc @@ -1079,6 +1079,57 @@ FiniteElement::has_support_on_face ( +template +const FiniteElement & +FiniteElement:: +get_sub_fe (const ComponentMask &mask) const +{ + // Translate the ComponentMask into first_selected and n_components after + // some error checking: + const unsigned int n_total_components = this->n_components(); + Assert((n_total_components == mask.size()) || (mask.size()==0), + ExcMessage("The given ComponentMask has the wrong size.")); + + const unsigned int n_selected = mask.n_selected_components(n_total_components); + Assert(n_selected>0, + ExcMessage("You need at least one selected component.")); + + const unsigned int first_selected = mask.first_selected_component(n_total_components); + +#ifdef DEBUG + // check that it is contiguous: + for (unsigned int c=0; c=first_selected && c=first_selected+n_selected && !mask[c]), + ExcMessage("Error: the given ComponentMask is not contiguous!")); +#endif + + return get_sub_fe(first_selected, n_selected); +} + + + +template +const FiniteElement & +FiniteElement:: +get_sub_fe (const unsigned int first_component, + const unsigned int n_selected_components) const +{ + // No complicated logic is needed here, because it is overridden in + // FESystem. Just make sure that what the user chose is valid: + Assert(first_component == 0 && n_selected_components == this->n_components(), + ExcMessage("You can only select a whole FiniteElement, not a part of one.")); + + (void)first_component; + (void)n_selected_components; + return *this; +} + + + template std::pair, std::vector > FiniteElement::get_constant_modes () const diff --git a/source/fe/fe_system.cc b/source/fe/fe_system.cc index c3b5785de0..22e03f11f7 100644 --- a/source/fe/fe_system.cc +++ b/source/fe/fe_system.cc @@ -71,7 +71,6 @@ FESystem::InternalData::~InternalData() } - template typename FiniteElement::InternalDataBase & FESystem:: @@ -333,6 +332,31 @@ FESystem::clone() const +template +const FiniteElement & +FESystem:: +get_sub_fe (const unsigned int first_component, + const unsigned int n_selected_components) const +{ + Assert(first_component+n_selected_components <= this->n_components(), + ExcMessage("Invalid arguments (not a part of this FiniteElement).")); + + const unsigned int base_index = this->component_to_base_table[first_component].first.first; + const unsigned int component_in_base = this->component_to_base_table[first_component].first.second; + const unsigned int base_components = this->base_element(base_index).n_components(); + + // Only select our child base_index if that is all the user wanted. Error + // handling will be done inside the recursion. + if (n_selected_components<=base_components) + return this->base_element(base_index).get_sub_fe(component_in_base, n_selected_components); + + Assert(n_selected_components == this->n_components(), + ExcMessage("You can not select a part of a FiniteElement.")); + return *this; +} + + + template double FESystem::shape_value (const unsigned int i, diff --git a/tests/fe/get_sub_fe_01.cc b/tests/fe/get_sub_fe_01.cc new file mode 100644 index 0000000000..8fedc690ef --- /dev/null +++ b/tests/fe/get_sub_fe_01.cc @@ -0,0 +1,180 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2002 - 2015 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// Test FiniteElement::get_sub_fe() + +#include "../tests.h" +#include + +#include +#include +#include +#include +#include +#include +#include + +template +void works(const FiniteElement &fe, const ComponentMask &m) +{ + deallog << "FE: " << fe.get_name() + << " mask: " << m << std::endl; + + const FiniteElement &child = fe.get_sub_fe(m); + + deallog << " worked: " << child.get_name() << std::endl; +} + +template +void fails(const FiniteElement &fe, const ComponentMask &m) +{ + deallog << "FE: " << fe.get_name() + << " mask: " << m << std::endl; + + try + { + const FiniteElement &child = fe.get_sub_fe(m); + deallog << " ERROR: we succeeded and got " << child.get_name() + << " but we should have failed!" << std::endl; + } + catch (...) + { + deallog << " failed as expected" << std::endl; + } +} + + +template +void check () +{ + auto mask_none = [] (const unsigned int n_components) -> ComponentMask + { + return ComponentMask(n_components, false); + }; + auto mask_all = [] (const unsigned int n_components) -> ComponentMask + { + return ComponentMask(n_components, true); + }; + auto mask_single = [] (const unsigned int n_components, + const unsigned int single) -> ComponentMask + { + ComponentMask c(n_components, false); + c.set(single, true); + return c; + }; + auto mask = [] (const unsigned int n_components, + const unsigned int first, + const unsigned int last) -> ComponentMask + { + ComponentMask c(n_components, false); + for (unsigned int i=first; i<=last; ++i) + c.set(i, true); + return c; + }; + + FE_Q fe_q(2); + FE_DGP fe_dg(2); + FE_BDM fe_nonprim(1); + + // simple FE: + { + works(fe_q, mask_all(1)); + fails(fe_q, mask_none(1)); + } + + // simple non-primitive: + { + works(fe_nonprim, ComponentMask()); // un-sized "select all" + works(fe_nonprim, mask_all(dim)); + fails(fe_nonprim, mask_none(dim)); + fails(fe_nonprim, mask_single(dim, 0)); // can not select part of it + } + + // remove system: + { + FESystem fe(fe_q, 1); + FESystem fe2(fe, 1); + works(fe, mask_all(1)); + works(fe2, mask_all(1)); + } + + // part of system: + { + FESystem fe(fe_q, 3); + works(fe, mask_all(3)); // keep system + works(fe, mask_single(3, 1)); // select one component + fails(fe, mask(3, 1, 2)); // select more than one component but not the whole FESystem + } + + // more systems: + { + FESystem fe(fe_nonprim,1,fe_dg,1); + works(fe, mask(dim+1, 0, dim)); // nonprimitive + works(fe, mask_single(dim+1, dim)); // select one component + } + { + FESystem fe(fe_nonprim,1,fe_dg,2); + // non-contiguous is a fail! + auto m = mask(dim+2, 0, dim); + m.set(1, false); + m.set(dim+1,true); + fails(fe, m); + } + { + FESystem fe(fe_q,dim,fe_q,1); + works(fe, mask_single(dim+1, dim)); + fails(fe, mask(dim+1, 0, dim-1)); // can not select the dim fe_q + } + { + FESystem qs(fe_q,dim); + FESystem fe(qs,1,fe_q,1); + works(fe, mask(dim+1, 0, dim-1)); // but works if nested + } + + { + FESystem fe(fe_q,2,fe_nonprim,2,fe_dg,1); + works(fe, mask_single(2*dim+3, 0)); + works(fe, mask_single(2*dim+3, 1)); + works(fe, mask_single(2*dim+3, 2*dim+3-1)); + works(fe, mask(2*dim+3, 2, 2+dim-1)); // 1st nonprimitive + works(fe, mask(2*dim+3, 2+dim, 2+2*dim-1)); // 2nd nonprimitive + } + + // nesting doll: + { + FESystem inner(fe_dg,1); + FESystem fe(fe_nonprim,1,inner,2); + FESystem outer(fe,2); + FESystem outer_outer(outer,1); + + works(fe, mask_single(dim+2, dim)); + works(fe, mask_single(dim+2, dim+1)); + works(outer, mask_single(2*(dim+2), dim)); + works(outer_outer, mask_single(2*(dim+2), dim)); + } + +} + + + + +int main () +{ + initlog(); + deal_II_exceptions::disable_abort_on_exception(); + + check<2> (); +} + diff --git a/tests/fe/get_sub_fe_01.debug.output b/tests/fe/get_sub_fe_01.debug.output new file mode 100644 index 0000000000..72d8177bf0 --- /dev/null +++ b/tests/fe/get_sub_fe_01.debug.output @@ -0,0 +1,53 @@ + +DEAL::FE: FE_Q<2>(2) mask: [true] +DEAL:: worked: FE_Q<2>(2) +DEAL::FE: FE_Q<2>(2) mask: [false] +DEAL:: failed as expected +DEAL::FE: FE_BDM<2>(1) mask: [all components selected] +DEAL:: worked: FE_BDM<2>(1) +DEAL::FE: FE_BDM<2>(1) mask: [true,true] +DEAL:: worked: FE_BDM<2>(1) +DEAL::FE: FE_BDM<2>(1) mask: [false,false] +DEAL:: failed as expected +DEAL::FE: FE_BDM<2>(1) mask: [true,false] +DEAL:: failed as expected +DEAL::FE: FESystem<2>[FE_Q<2>(2)] mask: [true] +DEAL:: worked: FE_Q<2>(2) +DEAL::FE: FESystem<2>[FESystem<2>[FE_Q<2>(2)]] mask: [true] +DEAL:: worked: FE_Q<2>(2) +DEAL::FE: FESystem<2>[FE_Q<2>(2)^3] mask: [true,true,true] +DEAL:: worked: FESystem<2>[FE_Q<2>(2)^3] +DEAL::FE: FESystem<2>[FE_Q<2>(2)^3] mask: [false,true,false] +DEAL:: worked: FE_Q<2>(2) +DEAL::FE: FESystem<2>[FE_Q<2>(2)^3] mask: [false,true,true] +DEAL:: failed as expected +DEAL::FE: FESystem<2>[FE_BDM<2>(1)-FE_DGP<2>(2)] mask: [true,true,true] +DEAL:: worked: FESystem<2>[FE_BDM<2>(1)-FE_DGP<2>(2)] +DEAL::FE: FESystem<2>[FE_BDM<2>(1)-FE_DGP<2>(2)] mask: [false,false,true] +DEAL:: worked: FE_DGP<2>(2) +DEAL::FE: FESystem<2>[FE_BDM<2>(1)-FE_DGP<2>(2)^2] mask: [true,false,true,true] +DEAL:: failed as expected +DEAL::FE: FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(2)] mask: [false,false,true] +DEAL:: worked: FE_Q<2>(2) +DEAL::FE: FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(2)] mask: [true,true,false] +DEAL:: failed as expected +DEAL::FE: FESystem<2>[FESystem<2>[FE_Q<2>(2)^2]-FE_Q<2>(2)] mask: [true,true,false] +DEAL:: worked: FESystem<2>[FE_Q<2>(2)^2] +DEAL::FE: FESystem<2>[FE_Q<2>(2)^2-FE_BDM<2>(1)^2-FE_DGP<2>(2)] mask: [true,false,false,false,false,false,false] +DEAL:: worked: FE_Q<2>(2) +DEAL::FE: FESystem<2>[FE_Q<2>(2)^2-FE_BDM<2>(1)^2-FE_DGP<2>(2)] mask: [false,true,false,false,false,false,false] +DEAL:: worked: FE_Q<2>(2) +DEAL::FE: FESystem<2>[FE_Q<2>(2)^2-FE_BDM<2>(1)^2-FE_DGP<2>(2)] mask: [false,false,false,false,false,false,true] +DEAL:: worked: FE_DGP<2>(2) +DEAL::FE: FESystem<2>[FE_Q<2>(2)^2-FE_BDM<2>(1)^2-FE_DGP<2>(2)] mask: [false,false,true,true,false,false,false] +DEAL:: worked: FE_BDM<2>(1) +DEAL::FE: FESystem<2>[FE_Q<2>(2)^2-FE_BDM<2>(1)^2-FE_DGP<2>(2)] mask: [false,false,false,false,true,true,false] +DEAL:: worked: FE_BDM<2>(1) +DEAL::FE: FESystem<2>[FE_BDM<2>(1)-FESystem<2>[FE_DGP<2>(2)]^2] mask: [false,false,true,false] +DEAL:: worked: FE_DGP<2>(2) +DEAL::FE: FESystem<2>[FE_BDM<2>(1)-FESystem<2>[FE_DGP<2>(2)]^2] mask: [false,false,false,true] +DEAL:: worked: FE_DGP<2>(2) +DEAL::FE: FESystem<2>[FESystem<2>[FE_BDM<2>(1)-FESystem<2>[FE_DGP<2>(2)]^2]^2] mask: [false,false,true,false,false,false,false,false] +DEAL:: worked: FE_DGP<2>(2) +DEAL::FE: FESystem<2>[FESystem<2>[FESystem<2>[FE_BDM<2>(1)-FESystem<2>[FE_DGP<2>(2)]^2]^2]] mask: [false,false,true,false,false,false,false,false] +DEAL:: worked: FE_DGP<2>(2) diff --git a/tests/fe/get_sub_fe_02.cc b/tests/fe/get_sub_fe_02.cc new file mode 100644 index 0000000000..051f3ae8bc --- /dev/null +++ b/tests/fe/get_sub_fe_02.cc @@ -0,0 +1,81 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2002 - 2015 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// Test FiniteElement::get_sub_fe(), example used in documentation + +#include "../tests.h" +#include + +#include +#include +#include +#include +#include +#include +#include + +template +void check () +{ + FESystem<2> fe_velocity(FE_Q<2>(2), 2); + FE_Q<2> fe_pressure(1); + FE_DGP<2> fe_dg(0); + FE_BDM<2> fe_nonprim(1); + + FESystem<2> fe(fe_velocity, 1, fe_pressure, 1, fe_dg, 2, fe_nonprim, 1); + + // same using component masks to copy over: + auto run = [&](const unsigned int first, + const unsigned int n, + const std::string &desc) + { + const unsigned int n_components = fe.n_components(); + + ComponentMask mask(n_components, false); + for (unsigned int i=first; i\n" + << "" << mask << "\n" + << "" << fe.get_sub_fe(mask).get_name() << "\n" + << "" << desc << "\n" + << "\n"; + + // we should be able to use ComponentMask or indices: + Assert(fe.get_sub_fe(mask) == fe.get_sub_fe(first, n), + ExcInternalError()); + }; + + deallog << "\n\n" + << "\n\n\n\n\n"; + run(0, 7, "@p fe itself, the whole @p FESystem"); + run(0, 2, "just the @p fe_velocity"); + run(0, 1, "The first component in @p fe_velocity"); + run(1, 1, "The second component in @p fe_velocity"); + run(2, 1, "@p fe_pressure"); + run(3, 1, "first copy of @p fe_dg"); + run(4, 1, "second copy of @p fe_dg"); + run(5, 2, "both components of @p fe_nonprim"); + deallog << "
ComponentMaskResultDescription
" << std::endl; +} + +int main () +{ + initlog(); + + check<2> (); +} + diff --git a/tests/fe/get_sub_fe_02.output b/tests/fe/get_sub_fe_02.output new file mode 100644 index 0000000000..b0fb89c175 --- /dev/null +++ b/tests/fe/get_sub_fe_02.output @@ -0,0 +1,49 @@ + +DEAL:: + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
ComponentMaskResultDescription
[true,true,true,true,true,true,true]FESystem<2>[FESystem<2>[FE_Q<2>(2)^2]-FE_Q<2>(1)-FE_DGP<2>(0)^2-FE_BDM<2>(1)]@p fe itself, the whole @p FESystem
[true,true,false,false,false,false,false]FESystem<2>[FE_Q<2>(2)^2]just the @p fe_velocity
[true,false,false,false,false,false,false]FE_Q<2>(2)The first component in @p fe_velocity
[false,true,false,false,false,false,false]FE_Q<2>(2)The second component in @p fe_velocity
[false,false,true,false,false,false,false]FE_Q<2>(1)@p fe_pressure
[false,false,false,true,false,false,false]FE_DGP<2>(0)first copy of @p fe_dg
[false,false,false,false,true,false,false]FE_DGP<2>(0)second copy of @p fe_dg
[false,false,false,false,false,true,true]FE_BDM<2>(1)both components of @p fe_nonprim
-- 2.39.5