From 188cde7ac48c8147d829617012e0328ec84f4818 Mon Sep 17 00:00:00 2001 From: Marc Fehling Date: Tue, 16 Oct 2018 17:49:40 -0600 Subject: [PATCH] Refactored FESeries to allow spacedim template for FECollection. --- doc/news/changes/minor/20181017MarcFehling | 5 + include/deal.II/fe/fe_series.h | 36 +- source/fe/CMakeLists.txt | 4 + source/fe/fe_series.cc | 361 +-------------------- source/fe/fe_series_fourier.cc | 221 +++++++++++++ source/fe/fe_series_fourier.inst.in | 24 ++ source/fe/fe_series_legendre.cc | 239 ++++++++++++++ source/fe/fe_series_legendre.inst.in | 24 ++ 8 files changed, 530 insertions(+), 384 deletions(-) create mode 100644 doc/news/changes/minor/20181017MarcFehling create mode 100644 source/fe/fe_series_fourier.cc create mode 100644 source/fe/fe_series_fourier.inst.in create mode 100644 source/fe/fe_series_legendre.cc create mode 100644 source/fe/fe_series_legendre.inst.in diff --git a/doc/news/changes/minor/20181017MarcFehling b/doc/news/changes/minor/20181017MarcFehling new file mode 100644 index 0000000000..61f020bbaa --- /dev/null +++ b/doc/news/changes/minor/20181017MarcFehling @@ -0,0 +1,5 @@ +New: Classes FESeries::Fourier and FESeries::Legendre now take two +template arguments dim and spacedim, which allow them to use corresponding +hp::FECollection objects. +
+(Marc Fehling, 2018/10/17) diff --git a/include/deal.II/fe/fe_series.h b/include/deal.II/fe/fe_series.h index 652ce9c24b..1b76e67e54 100644 --- a/include/deal.II/fe/fe_series.h +++ b/include/deal.II/fe/fe_series.h @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// Copyright (C) 2016 - 2017 by the deal.II authors +// Copyright (C) 2016 - 2018 by the deal.II authors // // This file is part of the deal.II library. // @@ -87,7 +87,7 @@ namespace FESeries * * @author Denis Davydov, 2016. */ - template + template class Fourier : public Subscriptor { public: @@ -98,9 +98,9 @@ namespace FESeries * used to integrate the expansion for each FiniteElement * in @p fe_collection. */ - Fourier(const unsigned int size_in_each_direction, - const hp::FECollection &fe_collection, - const hp::QCollection & q_collection); + Fourier(const unsigned int size_in_each_direction, + const hp::FECollection &fe_collection, + const hp::QCollection & q_collection); /** * Calculate @p fourier_coefficients of the cell vector field given by @@ -116,20 +116,13 @@ namespace FESeries /** * hp::FECollection for which transformation matrices will be calculated. */ - SmartPointer> fe_collection; + SmartPointer> fe_collection; /** * hp::QCollection used in calculation of transformation matrices. */ SmartPointer> q_collection; - /** - * Ensure that the transformation matrix for FiniteElement index - * @p fe_index is calculated. If not, calculate it. - */ - void - ensure_existence(const unsigned int fe_index); - /** * Angular frequencies $ 2 \pi {\bf k} $ . */ @@ -189,7 +182,7 @@ namespace FESeries * * @author Denis Davydov, 2016. */ - template + template class Legendre : public Subscriptor { public: @@ -200,9 +193,9 @@ namespace FESeries * used to integrate the expansion for each FiniteElement * in @p fe_collection. */ - Legendre(const unsigned int size_in_each_direction, - const hp::FECollection &fe_collection, - const hp::QCollection & q_collection); + Legendre(const unsigned int size_in_each_direction, + const hp::FECollection &fe_collection, + const hp::QCollection & q_collection); /** * Calculate @p legendre_coefficients of the cell vector field given by @@ -223,20 +216,13 @@ namespace FESeries /** * hp::FECollection for which transformation matrices will be calculated. */ - SmartPointer> fe_collection; + SmartPointer> fe_collection; /** * hp::QCollection used in calculation of transformation matrices. */ SmartPointer> q_collection; - /** - * Ensure that the transformation matrix for FiniteElement index - * @p fe_index is calculated. If not, calculate it. - */ - void - ensure_existence(const unsigned int fe_index); - /** * Transformation matrices for each FiniteElement. */ diff --git a/source/fe/CMakeLists.txt b/source/fe/CMakeLists.txt index a82f5c9947..79a5bb8466 100644 --- a/source/fe/CMakeLists.txt +++ b/source/fe/CMakeLists.txt @@ -47,6 +47,8 @@ SET(_unity_include_src fe_raviart_thomas_nodal.cc fe_rt_bubbles.cc fe_series.cc + fe_series_fourier.cc + fe_series_legendre.cc fe_system.cc fe_enriched.cc fe_tools.cc @@ -110,6 +112,8 @@ SET(_inst fe_raviart_thomas.inst.in fe_raviart_thomas_nodal.inst.in fe_rt_bubbles.inst.in + fe_series_fourier.inst.in + fe_series_legendre.inst.in fe_system.inst.in fe_enriched.inst.in fe_tools.inst.in diff --git a/source/fe/fe_series.cc b/source/fe/fe_series.cc index f2979b92cf..e1f7ca40f8 100644 --- a/source/fe/fe_series.cc +++ b/source/fe/fe_series.cc @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// Copyright (C) 2016 by the deal.II authors +// Copyright (C) 2016 - 2018 by the deal.II authors // // This file is part of the deal.II library. // @@ -17,358 +17,15 @@ #include -#include - #include -#include #include -#ifdef DEAL_II_WITH_GSL -# include -#endif DEAL_II_NAMESPACE_OPEN namespace FESeries { - /*-------------- Fourier -------------------------------*/ - - void set_k_vectors(Table<1, Tensor<1, 1>> &k_vectors, const unsigned int N) - { - k_vectors.reinit(TableIndices<1>(N)); - for (unsigned int i = 0; i < N; ++i) - k_vectors(i)[0] = 2. * numbers::PI * i; - } - - void set_k_vectors(Table<2, Tensor<1, 2>> &k_vectors, const unsigned int N) - { - k_vectors.reinit(TableIndices<2>(N, N)); - for (unsigned int i = 0; i < N; ++i) - for (unsigned int j = 0; j < N; ++j) - { - k_vectors(i, j)[0] = 2. * numbers::PI * i; - k_vectors(i, j)[1] = 2. * numbers::PI * j; - } - } - - void set_k_vectors(Table<3, Tensor<1, 3>> &k_vectors, const unsigned int N) - { - k_vectors.reinit(TableIndices<3>(N, N, N)); - for (unsigned int i = 0; i < N; ++i) - for (unsigned int j = 0; j < N; ++j) - for (unsigned int k = 0; k < N; ++k) - { - k_vectors(i, j, k)[0] = 2. * numbers::PI * i; - k_vectors(i, j, k)[0] = 2. * numbers::PI * j; - k_vectors(i, j, k)[0] = 2. * numbers::PI * k; - } - } - - - template - Fourier::Fourier(const unsigned int N, - const hp::FECollection &fe_collection, - const hp::QCollection & q_collection) - : fe_collection(&fe_collection) - , q_collection(&q_collection) - , fourier_transform_matrices(fe_collection.size()) - { - set_k_vectors(k_vectors, N); - unrolled_coefficients.resize(k_vectors.n_elements()); - } - - template - void - Fourier::calculate( - const Vector & local_dof_values, - const unsigned int cell_active_fe_index, - Table> &fourier_coefficients) - { - ensure_existence(cell_active_fe_index); - const FullMatrix> &matrix = - fourier_transform_matrices[cell_active_fe_index]; - - std::fill(unrolled_coefficients.begin(), - unrolled_coefficients.end(), - std::complex(0.)); - - Assert(unrolled_coefficients.size() == matrix.m(), ExcInternalError()); - - Assert(local_dof_values.size() == matrix.n(), - ExcDimensionMismatch(local_dof_values.size(), matrix.n())); - - for (unsigned int i = 0; i < unrolled_coefficients.size(); i++) - for (unsigned int j = 0; j < local_dof_values.size(); j++) - unrolled_coefficients[i] += matrix[i][j] * local_dof_values[j]; - - fourier_coefficients.fill(unrolled_coefficients.begin()); - } - - template - std::complex - integrate(const FiniteElement &fe, - const Quadrature & quadrature, - const Tensor<1, dim> & k_vector, - const unsigned int j) - { - 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); - } - return sum; - } - - - template <> - void - Fourier<2>::ensure_existence(const unsigned int fe) - { - Assert(fe < fe_collection->size(), - ExcIndexRange(fe, 0, fe_collection->size())) - - if (fourier_transform_matrices[fe].m() == 0) - { - fourier_transform_matrices[fe].reinit(k_vectors.n_elements(), - (*fe_collection)[fe].dofs_per_cell); - - unsigned int k = 0; - for (unsigned int k1 = 0; k1 < k_vectors.size(0); ++k1) - for (unsigned int k2 = 0; k2 < k_vectors.size(1); ++k2, k++) - for (unsigned int j = 0; j < (*fe_collection)[fe].dofs_per_cell; ++j) - fourier_transform_matrices[fe](k, j) = integrate( - (*fe_collection)[fe], (*q_collection)[fe], k_vectors(k1, k2), j); - } - } - - template <> - void - Fourier<3>::ensure_existence(const unsigned int fe) - { - Assert(fe < fe_collection->size(), - ExcIndexRange(fe, 0, fe_collection->size())) - - if (fourier_transform_matrices[fe].m() == 0) - { - fourier_transform_matrices[fe].reinit(k_vectors.n_elements(), - (*fe_collection)[fe].dofs_per_cell); - - unsigned int k = 0; - for (unsigned int k1 = 0; k1 < k_vectors.size(0); ++k1) - for (unsigned int k2 = 0; k2 < k_vectors.size(1); ++k2) - for (unsigned int k3 = 0; k3 < k_vectors.size(2); ++k3, k++) - for (unsigned int j = 0; j < (*fe_collection)[fe].dofs_per_cell; - ++j) - fourier_transform_matrices[fe](k, j) = - integrate((*fe_collection)[fe], - (*q_collection)[fe], - k_vectors(k1, k2, k3), - j); - } - } - - template <> - void - Fourier<1>::ensure_existence(const unsigned int fe) - { - Assert(fe < fe_collection->size(), - ExcIndexRange(fe, 0, fe_collection->size())) - - if (fourier_transform_matrices[fe].m() == 0) - { - fourier_transform_matrices[fe].reinit(k_vectors.n_elements(), - (*fe_collection)[fe].dofs_per_cell); - - for (unsigned int k = 0; k < k_vectors.size(0); ++k) - for (unsigned int j = 0; j < (*fe_collection)[fe].dofs_per_cell; ++j) - fourier_transform_matrices[fe](k, j) = integrate((*fe_collection)[fe], - (*q_collection)[fe], - k_vectors(k), - j); - } - } - - - /*-------------- Legendre -------------------------------*/ - DeclException2(ExcLegendre, - int, - double, - << "x[" << arg1 << "] = " << arg2 << " is not in [0,1]"); - - /* dim dimensional Legendre function with indices @p indices - * evaluated at @p x_q in [0,1]^dim. - */ - template - double - Lh(const Point &x_q, const TableIndices &indices) - { -#ifdef DEAL_II_WITH_GSL - double res = 1.0; - for (unsigned int d = 0; d < dim; d++) - { - const double x = 2.0 * (x_q[d] - 0.5); - Assert((x_q[d] <= 1.0) && (x_q[d] >= 0.), ExcLegendre(d, x_q[d])); - const int ind = indices[d]; - res *= std::sqrt(2.0) * gsl_sf_legendre_Pl(ind, x); - } - return res; - -#else - - (void)x_q; - (void)indices; - AssertThrow(false, - ExcMessage("deal.II has to be configured with GSL" - "in order to use Legendre transformation.")); - return 0; -#endif - } - - /* - * Multiplier in Legendre coefficients - */ - template - double - multiplier(const TableIndices &indices) - { - double res = 1.0; - for (unsigned int d = 0; d < dim; d++) - res *= (0.5 + indices[d]); - - return res; - } - - - template - Legendre::Legendre(const unsigned int size_in_each_direction, - const hp::FECollection &fe_collection, - const hp::QCollection & q_collection) - : N(size_in_each_direction) - , fe_collection(&fe_collection) - , q_collection(&q_collection) - , legendre_transform_matrices(fe_collection.size()) - , unrolled_coefficients(Utilities::fixed_power(N), 0.) - {} - - template - void - Legendre::calculate(const dealii::Vector &local_dof_values, - const unsigned int cell_active_fe_index, - Table & legendre_coefficients) - { - ensure_existence(cell_active_fe_index); - const FullMatrix &matrix = - legendre_transform_matrices[cell_active_fe_index]; - - std::fill(unrolled_coefficients.begin(), unrolled_coefficients.end(), 0.); - - Assert(unrolled_coefficients.size() == matrix.m(), ExcInternalError()); - - Assert(local_dof_values.size() == matrix.n(), - ExcDimensionMismatch(local_dof_values.size(), matrix.n())); - - for (unsigned int i = 0; i < unrolled_coefficients.size(); i++) - for (unsigned int j = 0; j < local_dof_values.size(); j++) - unrolled_coefficients[i] += matrix[i][j] * local_dof_values[j]; - - legendre_coefficients.fill(unrolled_coefficients.begin()); - } - - - template - double - integrate_Legendre(const FiniteElement &fe, - const Quadrature & quadrature, - const TableIndices & indices, - const unsigned int dof) - { - 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); - } - return sum * multiplier(indices); - } - - template <> - void - Legendre<1>::ensure_existence(const unsigned int fe) - { - Assert(fe < fe_collection->size(), - ExcIndexRange( - fe, - 0, - fe_collection->size())) if (legendre_transform_matrices[fe].m() == - 0) - { - legendre_transform_matrices[fe].reinit( - N, (*fe_collection)[fe].dofs_per_cell); - - for (unsigned int k = 0; k < N; ++k) - for (unsigned int j = 0; j < (*fe_collection)[fe].dofs_per_cell; ++j) - legendre_transform_matrices[fe](k, j) = integrate_Legendre( - (*fe_collection)[fe], (*q_collection)[fe], TableIndices<1>(k), j); - } - } - - - - template <> - void - Legendre<2>::ensure_existence(const unsigned int fe) - { - Assert(fe < fe_collection->size(), - ExcIndexRange(fe, 0, fe_collection->size())) - - if (legendre_transform_matrices[fe].m() == 0) - { - legendre_transform_matrices[fe].reinit( - N * N, (*fe_collection)[fe].dofs_per_cell); - - unsigned int k = 0; - for (unsigned int k1 = 0; k1 < N; ++k1) - for (unsigned int k2 = 0; k2 < N; ++k2, k++) - for (unsigned int j = 0; j < (*fe_collection)[fe].dofs_per_cell; ++j) - legendre_transform_matrices[fe](k, j) = - integrate_Legendre((*fe_collection)[fe], - (*q_collection)[fe], - TableIndices<2>(k1, k2), - j); - } - } - - template <> - void - Legendre<3>::ensure_existence(const unsigned int fe) - { - Assert(fe < fe_collection->size(), - ExcIndexRange(fe, 0, fe_collection->size())) - - if (legendre_transform_matrices[fe].m() == 0) - { - legendre_transform_matrices[fe].reinit( - N * N * N, (*fe_collection)[fe].dofs_per_cell); - - unsigned int k = 0; - for (unsigned int k1 = 0; k1 < N; ++k1) - for (unsigned int k2 = 0; k2 < N; ++k2) - for (unsigned int k3 = 0; k3 < N; ++k3, k++) - for (unsigned int j = 0; j < (*fe_collection)[fe].dofs_per_cell; - ++j) - legendre_transform_matrices[fe](k, j) = - integrate_Legendre((*fe_collection)[fe], - (*q_collection)[fe], - TableIndices<3>(k1, k2, k3), - j); - } - } - - /*-------------- linear_regression -------------------------------*/ std::pair linear_regression(const std::vector &x, const std::vector &y) { @@ -406,20 +63,6 @@ namespace FESeries return std::make_pair(X(1), X(0)); } - - -} // end of namespace FESeries - - - -/*-------------- Explicit Instantiations -------------------------------*/ -template class FESeries::Fourier<1>; -template class FESeries::Fourier<2>; -template class FESeries::Fourier<3>; -template class FESeries::Legendre<1>; -template class FESeries::Legendre<2>; -template class FESeries::Legendre<3>; - - +} // namespace FESeries DEAL_II_NAMESPACE_CLOSE diff --git a/source/fe/fe_series_fourier.cc b/source/fe/fe_series_fourier.cc new file mode 100644 index 0000000000..7763754986 --- /dev/null +++ b/source/fe/fe_series_fourier.cc @@ -0,0 +1,221 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 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. +// +// --------------------------------------------------------------------- + + + +#include + +#include + +#include + +#include + + +DEAL_II_NAMESPACE_OPEN + +namespace +{ + void set_k_vectors(Table<1, Tensor<1, 1>> &k_vectors, const unsigned int N) + { + k_vectors.reinit(TableIndices<1>(N)); + for (unsigned int i = 0; i < N; ++i) + k_vectors(i)[0] = 2. * numbers::PI * i; + } + + void set_k_vectors(Table<2, Tensor<1, 2>> &k_vectors, const unsigned int N) + { + k_vectors.reinit(TableIndices<2>(N, N)); + for (unsigned int i = 0; i < N; ++i) + for (unsigned int j = 0; j < N; ++j) + { + k_vectors(i, j)[0] = 2. * numbers::PI * i; + k_vectors(i, j)[1] = 2. * numbers::PI * j; + } + } + + void set_k_vectors(Table<3, Tensor<1, 3>> &k_vectors, const unsigned int N) + { + k_vectors.reinit(TableIndices<3>(N, N, N)); + for (unsigned int i = 0; i < N; ++i) + for (unsigned int j = 0; j < N; ++j) + for (unsigned int k = 0; k < N; ++k) + { + k_vectors(i, j, k)[0] = 2. * numbers::PI * i; + k_vectors(i, j, k)[0] = 2. * numbers::PI * j; + k_vectors(i, j, k)[0] = 2. * numbers::PI * k; + } + } + + + + template + std::complex + integrate(const FiniteElement &fe, + const Quadrature & quadrature, + const Tensor<1, dim> & k_vector, + const unsigned int j) + { + 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); + } + return sum; + } + + + + /* + * Ensure that the transformation matrix for FiniteElement index + * @p fe_index is calculated. If not, calculate it. + */ + template + void + ensure_existence( + const hp::FECollection<1, spacedim> & fe_collection, + const hp::QCollection<1> & q_collection, + const Table<1, Tensor<1, 1>> & k_vectors, + const unsigned int fe, + std::vector>> &fourier_transform_matrices) + { + AssertIndexRange(fe, fe_collection.size()); + + if (fourier_transform_matrices[fe].m() == 0) + { + fourier_transform_matrices[fe].reinit(k_vectors.n_elements(), + fe_collection[fe].dofs_per_cell); + + for (unsigned int k = 0; k < k_vectors.size(0); ++k) + for (unsigned int j = 0; j < fe_collection[fe].dofs_per_cell; ++j) + fourier_transform_matrices[fe](k, j) = + integrate(fe_collection[fe], q_collection[fe], k_vectors(k), j); + } + } + + template + void + ensure_existence( + const hp::FECollection<2, spacedim> & fe_collection, + const hp::QCollection<2> & q_collection, + const Table<2, Tensor<1, 2>> & k_vectors, + const unsigned int fe, + std::vector>> &fourier_transform_matrices) + { + AssertIndexRange(fe, fe_collection.size()); + + if (fourier_transform_matrices[fe].m() == 0) + { + fourier_transform_matrices[fe].reinit(k_vectors.n_elements(), + fe_collection[fe].dofs_per_cell); + + unsigned int k = 0; + for (unsigned int k1 = 0; k1 < k_vectors.size(0); ++k1) + for (unsigned int k2 = 0; k2 < k_vectors.size(1); ++k2, k++) + for (unsigned int j = 0; j < fe_collection[fe].dofs_per_cell; ++j) + fourier_transform_matrices[fe](k, j) = integrate( + fe_collection[fe], q_collection[fe], k_vectors(k1, k2), j); + } + } + + template + void + ensure_existence( + const hp::FECollection<3, spacedim> & fe_collection, + const hp::QCollection<3> & q_collection, + const Table<3, Tensor<1, 3>> & k_vectors, + const unsigned int fe, + std::vector>> &fourier_transform_matrices) + { + AssertIndexRange(fe, fe_collection.size()); + + if (fourier_transform_matrices[fe].m() == 0) + { + fourier_transform_matrices[fe].reinit(k_vectors.n_elements(), + fe_collection[fe].dofs_per_cell); + + unsigned int k = 0; + for (unsigned int k1 = 0; k1 < k_vectors.size(0); ++k1) + for (unsigned int k2 = 0; k2 < k_vectors.size(1); ++k2) + for (unsigned int k3 = 0; k3 < k_vectors.size(2); ++k3, k++) + for (unsigned int j = 0; j < fe_collection[fe].dofs_per_cell; ++j) + fourier_transform_matrices[fe](k, j) = + integrate(fe_collection[fe], + q_collection[fe], + k_vectors(k1, k2, k3), + j); + } + } +} // namespace + + + +namespace FESeries +{ + template + Fourier::Fourier( + const unsigned int N, + const hp::FECollection &fe_collection, + const hp::QCollection & q_collection) + : fe_collection(&fe_collection) + , q_collection(&q_collection) + , fourier_transform_matrices(fe_collection.size()) + { + set_k_vectors(k_vectors, N); + unrolled_coefficients.resize(k_vectors.n_elements()); + } + + + + template + void + Fourier::calculate( + const Vector & local_dof_values, + const unsigned int cell_active_fe_index, + Table> &fourier_coefficients) + { + ensure_existence(*fe_collection, + *q_collection, + k_vectors, + cell_active_fe_index, + fourier_transform_matrices); + + const FullMatrix> &matrix = + fourier_transform_matrices[cell_active_fe_index]; + + std::fill(unrolled_coefficients.begin(), + unrolled_coefficients.end(), + std::complex(0.)); + + Assert(unrolled_coefficients.size() == matrix.m(), ExcInternalError()); + + Assert(local_dof_values.size() == matrix.n(), + ExcDimensionMismatch(local_dof_values.size(), matrix.n())); + + for (unsigned int i = 0; i < unrolled_coefficients.size(); i++) + for (unsigned int j = 0; j < local_dof_values.size(); j++) + unrolled_coefficients[i] += matrix[i][j] * local_dof_values[j]; + + fourier_coefficients.fill(unrolled_coefficients.begin()); + } +} // namespace FESeries + + +// explicit instantiations +#include "fe_series_fourier.inst" + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/fe/fe_series_fourier.inst.in b/source/fe/fe_series_fourier.inst.in new file mode 100644 index 0000000000..d07d7c1b96 --- /dev/null +++ b/source/fe/fe_series_fourier.inst.in @@ -0,0 +1,24 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) + { +#if deal_II_dimension <= deal_II_space_dimension + template class FESeries::Fourier; +#endif + } diff --git a/source/fe/fe_series_legendre.cc b/source/fe/fe_series_legendre.cc new file mode 100644 index 0000000000..566c104863 --- /dev/null +++ b/source/fe/fe_series_legendre.cc @@ -0,0 +1,239 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 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. +// +// --------------------------------------------------------------------- + + + +#include + +#include + +#include +#ifdef DEAL_II_WITH_GSL +# include +#endif + + +DEAL_II_NAMESPACE_OPEN + +namespace +{ + DeclException2(ExcLegendre, + int, + double, + << "x[" << arg1 << "] = " << arg2 << " is not in [0,1]"); + + /* + * dim dimensional Legendre function with indices @p indices + * evaluated at @p x_q in [0,1]^dim. + */ + template + double + Lh(const Point &x_q, const TableIndices &indices) + { +#ifdef DEAL_II_WITH_GSL + double res = 1.0; + for (unsigned int d = 0; d < dim; d++) + { + const double x = 2.0 * (x_q[d] - 0.5); + Assert((x_q[d] <= 1.0) && (x_q[d] >= 0.), ExcLegendre(d, x_q[d])); + const int ind = indices[d]; + res *= std::sqrt(2.0) * gsl_sf_legendre_Pl(ind, x); + } + return res; + +#else + + (void)x_q; + (void)indices; + AssertThrow(false, + ExcMessage("deal.II has to be configured with GSL" + "in order to use Legendre transformation.")); + return 0; +#endif + } + + + + /* + * Multiplier in Legendre coefficients + */ + template + double + multiplier(const TableIndices &indices) + { + double res = 1.0; + for (unsigned int d = 0; d < dim; d++) + res *= (0.5 + indices[d]); + + return res; + } + + + + template + double + integrate(const FiniteElement &fe, + const Quadrature & quadrature, + const TableIndices & indices, + const unsigned int dof) + { + 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); + } + return sum * multiplier(indices); + } + + + + /** + * Ensure that the transformation matrix for FiniteElement index + * @p fe_index is calculated. If not, calculate it. + */ + template + void + ensure_existence(const hp::FECollection<1, spacedim> &fe_collection, + const hp::QCollection<1> & q_collection, + const unsigned int N, + const unsigned int fe, + std::vector> &legendre_transform_matrices) + { + AssertIndexRange(fe, fe_collection.size()); + + if (legendre_transform_matrices[fe].m() == 0) + { + legendre_transform_matrices[fe].reinit(N, + fe_collection[fe].dofs_per_cell); + + for (unsigned int k = 0; k < N; ++k) + for (unsigned int j = 0; j < fe_collection[fe].dofs_per_cell; ++j) + legendre_transform_matrices[fe](k, j) = integrate( + fe_collection[fe], q_collection[fe], TableIndices<1>(k), j); + } + } + + template + void + ensure_existence(const hp::FECollection<2, spacedim> &fe_collection, + const hp::QCollection<2> & q_collection, + const unsigned int N, + const unsigned int fe, + std::vector> &legendre_transform_matrices) + { + AssertIndexRange(fe, fe_collection.size()); + + if (legendre_transform_matrices[fe].m() == 0) + { + legendre_transform_matrices[fe].reinit(N * N, + fe_collection[fe].dofs_per_cell); + + unsigned int k = 0; + for (unsigned int k1 = 0; k1 < N; ++k1) + for (unsigned int k2 = 0; k2 < N; ++k2, k++) + for (unsigned int j = 0; j < fe_collection[fe].dofs_per_cell; ++j) + legendre_transform_matrices[fe](k, j) = + integrate(fe_collection[fe], + q_collection[fe], + TableIndices<2>(k1, k2), + j); + } + } + + template + void + ensure_existence(const hp::FECollection<3, spacedim> &fe_collection, + const hp::QCollection<3> & q_collection, + const unsigned int N, + const unsigned int fe, + std::vector> &legendre_transform_matrices) + { + AssertIndexRange(fe, fe_collection.size()); + + if (legendre_transform_matrices[fe].m() == 0) + { + legendre_transform_matrices[fe].reinit(N * N * N, + fe_collection[fe].dofs_per_cell); + + unsigned int k = 0; + for (unsigned int k1 = 0; k1 < N; ++k1) + for (unsigned int k2 = 0; k2 < N; ++k2) + for (unsigned int k3 = 0; k3 < N; ++k3, k++) + for (unsigned int j = 0; j < fe_collection[fe].dofs_per_cell; ++j) + legendre_transform_matrices[fe](k, j) = + integrate(fe_collection[fe], + q_collection[fe], + TableIndices<3>(k1, k2, k3), + j); + } + } +} // namespace + + + +namespace FESeries +{ + template + Legendre::Legendre( + const unsigned int size_in_each_direction, + const hp::FECollection &fe_collection, + const hp::QCollection & q_collection) + : N(size_in_each_direction) + , fe_collection(&fe_collection) + , q_collection(&q_collection) + , legendre_transform_matrices(fe_collection.size()) + , unrolled_coefficients(Utilities::fixed_power(N), 0.) + {} + + + + template + void + Legendre::calculate( + const dealii::Vector &local_dof_values, + const unsigned int cell_active_fe_index, + Table & legendre_coefficients) + { + ensure_existence(*fe_collection, + *q_collection, + N, + cell_active_fe_index, + legendre_transform_matrices); + + const FullMatrix &matrix = + legendre_transform_matrices[cell_active_fe_index]; + + std::fill(unrolled_coefficients.begin(), unrolled_coefficients.end(), 0.); + + Assert(unrolled_coefficients.size() == matrix.m(), ExcInternalError()); + + Assert(local_dof_values.size() == matrix.n(), + ExcDimensionMismatch(local_dof_values.size(), matrix.n())); + + for (unsigned int i = 0; i < unrolled_coefficients.size(); i++) + for (unsigned int j = 0; j < local_dof_values.size(); j++) + unrolled_coefficients[i] += matrix[i][j] * local_dof_values[j]; + + legendre_coefficients.fill(unrolled_coefficients.begin()); + } +} // namespace FESeries + + +// explicit instantiations +#include "fe_series_legendre.inst" + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/fe/fe_series_legendre.inst.in b/source/fe/fe_series_legendre.inst.in new file mode 100644 index 0000000000..587120ed8d --- /dev/null +++ b/source/fe/fe_series_legendre.inst.in @@ -0,0 +1,24 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) + { +#if deal_II_dimension <= deal_II_space_dimension + template class FESeries::Legendre; +#endif + } -- 2.39.5