From: Ce Qin Date: Sun, 21 Feb 2021 02:44:59 +0000 (+0800) Subject: Fix FESeries with non-primitive FEs. X-Git-Tag: v9.3.0-rc1~344^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ca90a6e06fb7cfa7cafd62c8f5eaaa3a2801191d;p=dealii.git Fix FESeries with non-primitive FEs. --- diff --git a/doc/news/changes/minor/20210301CeQin b/doc/news/changes/minor/20210301CeQin new file mode 100644 index 0000000000..b7e5050189 --- /dev/null +++ b/doc/news/changes/minor/20210301CeQin @@ -0,0 +1,4 @@ +New: Add a component parameter to FESeries::Fourier/Legendre to +make them working with non-primitive elements. +
+(Ce Qin, 2021/03/01) diff --git a/include/deal.II/fe/fe_series.h b/include/deal.II/fe/fe_series.h index 855720ca5e..34bf45a520 100644 --- a/include/deal.II/fe/fe_series.h +++ b/include/deal.II/fe/fe_series.h @@ -55,9 +55,10 @@ DEAL_II_NAMESPACE_OPEN namespace FESeries { /** - * A class to calculate expansion of a scalar FE field into Fourier series - * on a reference element. The exponential form of the Fourier series is - * based on completeness and Hermitian orthogonality of the set of exponential + * A class to calculate expansion of a scalar FE (or a single component + * of vector-valued FE) field into Fourier series on a reference element. + * The exponential form of the Fourier series is based on completeness + * and Hermitian orthogonality of the set of exponential * functions $ \phi_{\bf k}({\bf x}) = \exp(2 \pi i\, {\bf k} \cdot {\bf x})$. * For example in 1D the L2-orthogonality condition reads * @f[ @@ -96,10 +97,21 @@ namespace FESeries * each direction, @p fe_collection is the hp::FECollection for which * expansion will be used and @p q_collection is the hp::QCollection used to * integrate the expansion for each FiniteElement in @p fe_collection. + * + * As the Fourier expansion can only be performed on scalar fields, this + * class does not operate on vector-valued finite elements and will + * therefore throw an assertion. However, each component of a finite element + * field can be treated as a scalar field, respectively, on which Fourier + * expansions are again possible. For this purpose, the optional parameter + * @p component defines which component of each FiniteElement will be used. + * The default value of @p component only applies to scalar FEs, in which + * case it indicates that the sole component is to be decomposed. For + * vector-valued FEs, a non-default value must be explicitly provided. */ Fourier(const std::vector & n_coefficients_per_direction, const hp::FECollection &fe_collection, - const hp::QCollection & q_collection); + const hp::QCollection & q_collection, + const unsigned int component = numbers::invalid_unsigned_int); /** * A non-default constructor. The @p n_coefficients_per_direction defines the @@ -204,13 +216,20 @@ namespace FESeries * Auxiliary vector to store unrolled coefficients. */ std::vector unrolled_coefficients; + + /** + * Which component of FiniteElement should be used to calculate the + * expansion. + */ + const unsigned int component; }; /** - * A class to calculate expansion of a scalar FE field into series of Legendre - * functions on a reference element. + * A class to calculate expansion of a scalar FE (or a single component + * of vector-valued FE) field into series of Legendre functions on a + * reference element. * * Legendre functions are solutions to Legendre's differential equation * @f[ @@ -262,10 +281,21 @@ namespace FESeries * each direction, @p fe_collection is the hp::FECollection for which * expansion will be used and @p q_collection is the hp::QCollection used to * integrate the expansion for each FiniteElement in @p fe_collection. + * + * As the Legendre expansion can only be performed on scalar fields, this + * class does not operate on vector-valued finite elements and will + * therefore throw an assertion. However, each component of a finite element + * field can be treated as a scalar field, respectively, on which Legendre + * expansions are again possible. For this purpose, the optional parameter + * @p component defines which component of each FiniteElement will be used. + * The default value of @p component only applies to scalar FEs, in which + * case it indicates that the sole component is to be decomposed. For + * vector-valued FEs, a non-default value must be explicitly provided. */ Legendre(const std::vector &n_coefficients_per_direction, const hp::FECollection &fe_collection, - const hp::QCollection & q_collection); + const hp::QCollection & q_collection, + const unsigned int component = numbers::invalid_unsigned_int); /** * A non-default constructor. The @p size_in_each_direction defines the number @@ -364,6 +394,12 @@ namespace FESeries * Auxiliary vector to store unrolled coefficients. */ std::vector unrolled_coefficients; + + /** + * Which component of FiniteElement should be used to calculate the + * expansion. + */ + const unsigned int component; }; diff --git a/include/deal.II/numerics/smoothness_estimator.h b/include/deal.II/numerics/smoothness_estimator.h index b0d11b6531..77434ff39d 100644 --- a/include/deal.II/numerics/smoothness_estimator.h +++ b/include/deal.II/numerics/smoothness_estimator.h @@ -232,10 +232,22 @@ namespace SmoothnessEstimator * polynomial which is just a constant. Further for each element, we use a * Gaussian quadrature designed to yield exact results for the highest order * Legendre polynomial used. + * + * As the Legendre expansion can only be performed on scalar fields, this + * class does not operate on vector-valued finite elements and will + * therefore throw an assertion. However, each component of a finite element + * field can be treated as a scalar field, respectively, on which Legendre + * expansions are again possible. For this purpose, the optional parameter + * @p component defines which component of each FiniteElement will be used. + * The default value of @p component only applies to scalar FEs, in which + * case it indicates that the sole component is to be decomposed. For + * vector-valued FEs, a non-default value must be explicitly provided. */ template FESeries::Legendre - default_fe_series(const hp::FECollection &fe_collection); + default_fe_series( + const hp::FECollection &fe_collection, + const unsigned int component = numbers::invalid_unsigned_int); } // namespace Legendre @@ -458,10 +470,22 @@ namespace SmoothnessEstimator * a 5-point Gaussian quarature iterated in each dimension by the maximal * wave number, which is the number of modes decreased by one since we start * with $k = 0$. + * + * As the Fourier expansion can only be performed on scalar fields, this + * class does not operate on vector-valued finite elements and will + * therefore throw an assertion. However, each component of a finite element + * field can be treated as a scalar field, respectively, on which Fourier + * expansions are again possible. For this purpose, the optional parameter + * @p component defines which component of each FiniteElement will be used. + * The default value of @p component only applies to scalar FEs, in which + * case it indicates that the sole component is to be decomposed. For + * vector-valued FEs, a non-default value must be explicitly provided. */ template FESeries::Fourier - default_fe_series(const hp::FECollection &fe_collection); + default_fe_series( + const hp::FECollection &fe_collection, + const unsigned int component = numbers::invalid_unsigned_int); } // namespace Fourier } // namespace SmoothnessEstimator diff --git a/source/fe/fe_series_fourier.cc b/source/fe/fe_series_fourier.cc index 5fb7aede3b..9c3e54ddea 100644 --- a/source/fe/fe_series_fourier.cc +++ b/source/fe/fe_series_fourier.cc @@ -65,14 +65,16 @@ namespace integrate(const FiniteElement &fe, const Quadrature & quadrature, const Tensor<1, dim> & k_vector, - const unsigned int j) + const unsigned int j, + const unsigned int component) { std::complex sum = 0; for (unsigned int q = 0; q < quadrature.size(); ++q) { const Point &x_q = quadrature.point(q); sum += std::exp(std::complex(0, 1) * (k_vector * x_q)) * - fe.shape_value(j, x_q) * quadrature.weight(q); + fe.shape_value_component(j, x_q, component) * + quadrature.weight(q); } return sum; } @@ -91,6 +93,7 @@ namespace const hp::QCollection<1> & q_collection, const Table<1, Tensor<1, 1>> & k_vectors, const unsigned int fe, + const unsigned int component, std::vector>> &fourier_transform_matrices) { AssertIndexRange(fe, fe_collection.size()); @@ -103,8 +106,8 @@ namespace for (unsigned int k = 0; k < n_coefficients_per_direction[fe]; ++k) for (unsigned int j = 0; j < fe_collection[fe].n_dofs_per_cell(); ++j) - fourier_transform_matrices[fe](k, j) = - integrate(fe_collection[fe], q_collection[fe], k_vectors(k), j); + fourier_transform_matrices[fe](k, j) = integrate( + fe_collection[fe], q_collection[fe], k_vectors(k), j, component); } } @@ -116,6 +119,7 @@ namespace const hp::QCollection<2> & q_collection, const Table<2, Tensor<1, 2>> & k_vectors, const unsigned int fe, + const unsigned int component, std::vector>> &fourier_transform_matrices) { AssertIndexRange(fe, fe_collection.size()); @@ -132,8 +136,12 @@ namespace ++k2, ++k) for (unsigned int j = 0; j < fe_collection[fe].n_dofs_per_cell(); ++j) - fourier_transform_matrices[fe](k, j) = integrate( - fe_collection[fe], q_collection[fe], k_vectors(k1, k2), j); + fourier_transform_matrices[fe](k, j) = + integrate(fe_collection[fe], + q_collection[fe], + k_vectors(k1, k2), + j, + component); } } @@ -145,6 +153,7 @@ namespace const hp::QCollection<3> & q_collection, const Table<3, Tensor<1, 3>> & k_vectors, const unsigned int fe, + const unsigned int component, std::vector>> &fourier_transform_matrices) { AssertIndexRange(fe, fe_collection.size()); @@ -166,7 +175,8 @@ namespace integrate(fe_collection[fe], q_collection[fe], k_vectors(k1, k2, k3), - j); + j, + component); } } } // namespace @@ -179,16 +189,28 @@ namespace FESeries Fourier::Fourier( const std::vector & n_coefficients_per_direction, const hp::FECollection &fe_collection, - const hp::QCollection & q_collection) + const hp::QCollection & q_collection, + const unsigned int component_) : n_coefficients_per_direction(n_coefficients_per_direction) , fe_collection(&fe_collection) , q_collection(q_collection) , fourier_transform_matrices(fe_collection.size()) + , component(component_ != numbers::invalid_unsigned_int ? component_ : 0) { Assert(n_coefficients_per_direction.size() == fe_collection.size() && n_coefficients_per_direction.size() == q_collection.size(), ExcMessage("All parameters are supposed to have the same size.")); + if (fe_collection[0].n_components() > 1) + Assert( + component_ != numbers::invalid_unsigned_int, + ExcMessage( + "For vector-valued problems, you need to explicitly specify for " + "which vector component you will want to do a Fourier decomposition " + "by setting the 'component' argument of this constructor.")); + + AssertIndexRange(component, fe_collection[0].n_components()); + const unsigned int max_n_coefficients_per_direction = *std::max_element(n_coefficients_per_direction.cbegin(), n_coefficients_per_direction.cend()); @@ -224,7 +246,8 @@ namespace FESeries (*fe_collection == *(fourier.fe_collection)) && (q_collection == fourier.q_collection) && (k_vectors == fourier.k_vectors) && - (fourier_transform_matrices == fourier.fourier_transform_matrices)); + (fourier_transform_matrices == fourier.fourier_transform_matrices) && + (component == fourier.component)); } @@ -241,6 +264,7 @@ namespace FESeries q_collection, k_vectors, fe, + component, fourier_transform_matrices); }); @@ -276,6 +300,7 @@ namespace FESeries q_collection, k_vectors, cell_active_fe_index, + component, fourier_transform_matrices); const FullMatrix &matrix = diff --git a/source/fe/fe_series_legendre.cc b/source/fe/fe_series_legendre.cc index f74f217f3b..ffc4cc4283 100644 --- a/source/fe/fe_series_legendre.cc +++ b/source/fe/fe_series_legendre.cc @@ -74,14 +74,16 @@ namespace integrate(const FiniteElement &fe, const Quadrature & quadrature, const TableIndices & indices, - const unsigned int dof) + const unsigned int dof, + const unsigned int component) { double sum = 0; for (unsigned int q = 0; q < quadrature.size(); ++q) { const Point &x_q = quadrature.point(q); - sum += - Lh(x_q, indices) * fe.shape_value(dof, x_q) * quadrature.weight(q); + sum += Lh(x_q, indices) * + fe.shape_value_component(dof, x_q, component) * + quadrature.weight(q); } return sum * multiplier(indices); } @@ -99,6 +101,7 @@ namespace const hp::FECollection<1, spacedim> &fe_collection, const hp::QCollection<1> & q_collection, const unsigned int fe, + const unsigned int component, std::vector> & legendre_transform_matrices) { AssertIndexRange(fe, fe_collection.size()); @@ -111,8 +114,12 @@ namespace for (unsigned int k = 0; k < n_coefficients_per_direction[fe]; ++k) for (unsigned int j = 0; j < fe_collection[fe].n_dofs_per_cell(); ++j) - legendre_transform_matrices[fe](k, j) = integrate( - fe_collection[fe], q_collection[fe], TableIndices<1>(k), j); + legendre_transform_matrices[fe](k, j) = + integrate(fe_collection[fe], + q_collection[fe], + TableIndices<1>(k), + j, + component); } } @@ -123,6 +130,7 @@ namespace const hp::FECollection<2, spacedim> &fe_collection, const hp::QCollection<2> & q_collection, const unsigned int fe, + const unsigned int component, std::vector> & legendre_transform_matrices) { AssertIndexRange(fe, fe_collection.size()); @@ -143,7 +151,8 @@ namespace integrate(fe_collection[fe], q_collection[fe], TableIndices<2>(k1, k2), - j); + j, + component); } } @@ -154,6 +163,7 @@ namespace const hp::FECollection<3, spacedim> &fe_collection, const hp::QCollection<3> & q_collection, const unsigned int fe, + const unsigned int component, std::vector> & legendre_transform_matrices) { AssertIndexRange(fe, fe_collection.size()); @@ -175,7 +185,8 @@ namespace integrate(fe_collection[fe], q_collection[fe], TableIndices<3>(k1, k2, k3), - j); + j, + component); } } } // namespace @@ -188,16 +199,28 @@ namespace FESeries Legendre::Legendre( const std::vector & n_coefficients_per_direction, const hp::FECollection &fe_collection, - const hp::QCollection & q_collection) + const hp::QCollection & q_collection, + const unsigned int component_) : n_coefficients_per_direction(n_coefficients_per_direction) , fe_collection(&fe_collection) , q_collection(q_collection) , legendre_transform_matrices(fe_collection.size()) + , component(component_ != numbers::invalid_unsigned_int ? component_ : 0) { Assert(n_coefficients_per_direction.size() == fe_collection.size() && n_coefficients_per_direction.size() == q_collection.size(), ExcMessage("All parameters are supposed to have the same size.")); + if (fe_collection[0].n_components() > 1) + Assert( + component_ != numbers::invalid_unsigned_int, + ExcMessage( + "For vector-valued problems, you need to explicitly specify for " + "which vector component you will want to do a Fourier decomposition " + "by setting the 'component' argument of this constructor.")); + + AssertIndexRange(component, fe_collection[0].n_components()); + // reserve sufficient memory const unsigned int max_n_coefficients_per_direction = *std::max_element(n_coefficients_per_direction.cbegin(), @@ -231,7 +254,8 @@ namespace FESeries (n_coefficients_per_direction == legendre.n_coefficients_per_direction) && (*fe_collection == *(legendre.fe_collection)) && (q_collection == legendre.q_collection) && - (legendre_transform_matrices == legendre.legendre_transform_matrices)); + (legendre_transform_matrices == legendre.legendre_transform_matrices) && + (component == legendre.component)); } @@ -247,6 +271,7 @@ namespace FESeries *fe_collection, q_collection, fe, + component, legendre_transform_matrices); }); @@ -281,6 +306,7 @@ namespace FESeries *fe_collection, q_collection, cell_active_fe_index, + component, legendre_transform_matrices); const FullMatrix &matrix = diff --git a/source/numerics/smoothness_estimator.cc b/source/numerics/smoothness_estimator.cc index fa7c0c3c41..91c8c34bcb 100644 --- a/source/numerics/smoothness_estimator.cc +++ b/source/numerics/smoothness_estimator.cc @@ -287,7 +287,8 @@ namespace SmoothnessEstimator template FESeries::Legendre - default_fe_series(const hp::FECollection &fe_collection) + default_fe_series(const hp::FECollection &fe_collection, + const unsigned int component) { // Default number of coefficients per direction. // @@ -323,7 +324,8 @@ namespace SmoothnessEstimator return FESeries::Legendre(n_coefficients_per_direction, fe_collection, - q_collection); + q_collection, + component); } } // namespace Legendre @@ -574,7 +576,8 @@ namespace SmoothnessEstimator template FESeries::Fourier - default_fe_series(const hp::FECollection &fe_collection) + default_fe_series(const hp::FECollection &fe_collection, + const unsigned int component) { // Default number of coefficients per direction. // @@ -613,7 +616,8 @@ namespace SmoothnessEstimator return FESeries::Fourier(n_coefficients_per_direction, fe_collection, - q_collection); + q_collection, + component); } } // namespace Fourier } // namespace SmoothnessEstimator diff --git a/source/numerics/smoothness_estimator.inst.in b/source/numerics/smoothness_estimator.inst.in index f33653ce62..9f2b9112a3 100644 --- a/source/numerics/smoothness_estimator.inst.in +++ b/source/numerics/smoothness_estimator.inst.in @@ -67,11 +67,13 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) template FESeries::Legendre SmoothnessEstimator::Legendre::default_fe_series( - const hp::FECollection &); + const hp::FECollection &, + const unsigned int); template FESeries::Fourier SmoothnessEstimator::Fourier::default_fe_series( - const hp::FECollection &); + const hp::FECollection &, + const unsigned int); #endif } diff --git a/tests/numerics/smoothness_estimator_03.cc b/tests/numerics/smoothness_estimator_03.cc new file mode 100644 index 0000000000..7422150e1a --- /dev/null +++ b/tests/numerics/smoothness_estimator_03.cc @@ -0,0 +1,90 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 - 2021 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. +// +// --------------------------------------------------------------------- + + + +// This test is used to make sure that FESeries::Fourier/Legendre +// work with non-primitive elements + + +#include + +#include +#include + +#include +#include + +#include + +#include + +#include + +#include "../tests.h" + +template +void +test() +{ + hp::FECollection fe_collection; + fe_collection.push_back(FE_Nedelec(0)); + fe_collection.push_back(FE_Nedelec(1)); + + Triangulation tria; + GridGenerator::subdivided_hyper_cube(tria, 2); + + DoFHandler dh(tria); + dh.distribute_dofs(fe_collection); + + Vector solution(dh.n_dofs()); + Vector smoothness(tria.n_active_cells()); + + for (unsigned int i = 0; i < dim; ++i) + { + smoothness = 0.0; + FESeries::Fourier fourier = + SmoothnessEstimator::Fourier::default_fe_series(fe_collection, i); + SmoothnessEstimator::Fourier::coefficient_decay(fourier, + dh, + solution, + smoothness); + } + + + for (unsigned int i = 0; i < dim; ++i) + { + smoothness = 0.0; + FESeries::Legendre legendre = + SmoothnessEstimator::Legendre::default_fe_series(fe_collection, i); + SmoothnessEstimator::Legendre::coefficient_decay(legendre, + dh, + solution, + smoothness); + } + + deallog << "Ok" << std::endl; +} + +int +main() +{ + initlog(); + + test<2>(); + test<3>(); + + return 0; +} diff --git a/tests/numerics/smoothness_estimator_03.output b/tests/numerics/smoothness_estimator_03.output new file mode 100644 index 0000000000..67a2ddb88c --- /dev/null +++ b/tests/numerics/smoothness_estimator_03.output @@ -0,0 +1,3 @@ + +DEAL::Ok +DEAL::Ok