--- /dev/null
+New: Add a component parameter to FESeries::Fourier/Legendre to
+make them working with non-primitive elements.
+<br>
+(Ce Qin, 2021/03/01)
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[
* 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<unsigned int> & n_coefficients_per_direction,
const hp::FECollection<dim, spacedim> &fe_collection,
- const hp::QCollection<dim> & q_collection);
+ const hp::QCollection<dim> & q_collection,
+ const unsigned int component = numbers::invalid_unsigned_int);
/**
* A non-default constructor. The @p n_coefficients_per_direction defines the
* Auxiliary vector to store unrolled coefficients.
*/
std::vector<CoefficientType> 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[
* 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<unsigned int> &n_coefficients_per_direction,
const hp::FECollection<dim, spacedim> &fe_collection,
- const hp::QCollection<dim> & q_collection);
+ const hp::QCollection<dim> & q_collection,
+ const unsigned int component = numbers::invalid_unsigned_int);
/**
* A non-default constructor. The @p size_in_each_direction defines the number
* Auxiliary vector to store unrolled coefficients.
*/
std::vector<CoefficientType> unrolled_coefficients;
+
+ /**
+ * Which component of FiniteElement should be used to calculate the
+ * expansion.
+ */
+ const unsigned int component;
};
* 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 <int dim, int spacedim>
FESeries::Legendre<dim, spacedim>
- default_fe_series(const hp::FECollection<dim, spacedim> &fe_collection);
+ default_fe_series(
+ const hp::FECollection<dim, spacedim> &fe_collection,
+ const unsigned int component = numbers::invalid_unsigned_int);
} // namespace Legendre
* 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 <int dim, int spacedim>
FESeries::Fourier<dim, spacedim>
- default_fe_series(const hp::FECollection<dim, spacedim> &fe_collection);
+ default_fe_series(
+ const hp::FECollection<dim, spacedim> &fe_collection,
+ const unsigned int component = numbers::invalid_unsigned_int);
} // namespace Fourier
} // namespace SmoothnessEstimator
integrate(const FiniteElement<dim, spacedim> &fe,
const Quadrature<dim> & quadrature,
const Tensor<1, dim> & k_vector,
- const unsigned int j)
+ const unsigned int j,
+ const unsigned int component)
{
std::complex<double> sum = 0;
for (unsigned int q = 0; q < quadrature.size(); ++q)
{
const Point<dim> &x_q = quadrature.point(q);
sum += std::exp(std::complex<double>(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;
}
const hp::QCollection<1> & q_collection,
const Table<1, Tensor<1, 1>> & k_vectors,
const unsigned int fe,
+ const unsigned int component,
std::vector<FullMatrix<std::complex<double>>> &fourier_transform_matrices)
{
AssertIndexRange(fe, fe_collection.size());
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);
}
}
const hp::QCollection<2> & q_collection,
const Table<2, Tensor<1, 2>> & k_vectors,
const unsigned int fe,
+ const unsigned int component,
std::vector<FullMatrix<std::complex<double>>> &fourier_transform_matrices)
{
AssertIndexRange(fe, fe_collection.size());
++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);
}
}
const hp::QCollection<3> & q_collection,
const Table<3, Tensor<1, 3>> & k_vectors,
const unsigned int fe,
+ const unsigned int component,
std::vector<FullMatrix<std::complex<double>>> &fourier_transform_matrices)
{
AssertIndexRange(fe, fe_collection.size());
integrate(fe_collection[fe],
q_collection[fe],
k_vectors(k1, k2, k3),
- j);
+ j,
+ component);
}
}
} // namespace
Fourier<dim, spacedim>::Fourier(
const std::vector<unsigned int> & n_coefficients_per_direction,
const hp::FECollection<dim, spacedim> &fe_collection,
- const hp::QCollection<dim> & q_collection)
+ const hp::QCollection<dim> & 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());
(*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));
}
q_collection,
k_vectors,
fe,
+ component,
fourier_transform_matrices);
});
q_collection,
k_vectors,
cell_active_fe_index,
+ component,
fourier_transform_matrices);
const FullMatrix<CoefficientType> &matrix =
integrate(const FiniteElement<dim, spacedim> &fe,
const Quadrature<dim> & quadrature,
const TableIndices<dim> & 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<dim> &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);
}
const hp::FECollection<1, spacedim> &fe_collection,
const hp::QCollection<1> & q_collection,
const unsigned int fe,
+ const unsigned int component,
std::vector<FullMatrix<double>> & legendre_transform_matrices)
{
AssertIndexRange(fe, fe_collection.size());
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);
}
}
const hp::FECollection<2, spacedim> &fe_collection,
const hp::QCollection<2> & q_collection,
const unsigned int fe,
+ const unsigned int component,
std::vector<FullMatrix<double>> & legendre_transform_matrices)
{
AssertIndexRange(fe, fe_collection.size());
integrate(fe_collection[fe],
q_collection[fe],
TableIndices<2>(k1, k2),
- j);
+ j,
+ component);
}
}
const hp::FECollection<3, spacedim> &fe_collection,
const hp::QCollection<3> & q_collection,
const unsigned int fe,
+ const unsigned int component,
std::vector<FullMatrix<double>> & legendre_transform_matrices)
{
AssertIndexRange(fe, fe_collection.size());
integrate(fe_collection[fe],
q_collection[fe],
TableIndices<3>(k1, k2, k3),
- j);
+ j,
+ component);
}
}
} // namespace
Legendre<dim, spacedim>::Legendre(
const std::vector<unsigned int> & n_coefficients_per_direction,
const hp::FECollection<dim, spacedim> &fe_collection,
- const hp::QCollection<dim> & q_collection)
+ const hp::QCollection<dim> & 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(),
(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));
}
*fe_collection,
q_collection,
fe,
+ component,
legendre_transform_matrices);
});
*fe_collection,
q_collection,
cell_active_fe_index,
+ component,
legendre_transform_matrices);
const FullMatrix<CoefficientType> &matrix =
template <int dim, int spacedim>
FESeries::Legendre<dim, spacedim>
- default_fe_series(const hp::FECollection<dim, spacedim> &fe_collection)
+ default_fe_series(const hp::FECollection<dim, spacedim> &fe_collection,
+ const unsigned int component)
{
// Default number of coefficients per direction.
//
return FESeries::Legendre<dim, spacedim>(n_coefficients_per_direction,
fe_collection,
- q_collection);
+ q_collection,
+ component);
}
} // namespace Legendre
template <int dim, int spacedim>
FESeries::Fourier<dim, spacedim>
- default_fe_series(const hp::FECollection<dim, spacedim> &fe_collection)
+ default_fe_series(const hp::FECollection<dim, spacedim> &fe_collection,
+ const unsigned int component)
{
// Default number of coefficients per direction.
//
return FESeries::Fourier<dim, spacedim>(n_coefficients_per_direction,
fe_collection,
- q_collection);
+ q_collection,
+ component);
}
} // namespace Fourier
} // namespace SmoothnessEstimator
template FESeries::Legendre<deal_II_dimension, deal_II_space_dimension>
SmoothnessEstimator::Legendre::default_fe_series<deal_II_dimension,
deal_II_space_dimension>(
- const hp::FECollection<deal_II_dimension, deal_II_space_dimension> &);
+ const hp::FECollection<deal_II_dimension, deal_II_space_dimension> &,
+ const unsigned int);
template FESeries::Fourier<deal_II_dimension, deal_II_space_dimension>
SmoothnessEstimator::Fourier::default_fe_series<deal_II_dimension,
deal_II_space_dimension>(
- const hp::FECollection<deal_II_dimension, deal_II_space_dimension> &);
+ const hp::FECollection<deal_II_dimension, deal_II_space_dimension> &,
+ const unsigned int);
#endif
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_nedelec.h>
+#include <deal.II/fe/fe_series.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/hp/fe_collection.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/smoothness_estimator.h>
+
+#include "../tests.h"
+
+template <int dim>
+void
+test()
+{
+ hp::FECollection<dim> fe_collection;
+ fe_collection.push_back(FE_Nedelec<dim>(0));
+ fe_collection.push_back(FE_Nedelec<dim>(1));
+
+ Triangulation<dim> tria;
+ GridGenerator::subdivided_hyper_cube(tria, 2);
+
+ DoFHandler<dim> dh(tria);
+ dh.distribute_dofs(fe_collection);
+
+ Vector<double> solution(dh.n_dofs());
+ Vector<float> smoothness(tria.n_active_cells());
+
+ for (unsigned int i = 0; i < dim; ++i)
+ {
+ smoothness = 0.0;
+ FESeries::Fourier<dim> 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<dim> 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;
+}
--- /dev/null
+
+DEAL::Ok
+DEAL::Ok
projects / dealii.git / commitdiff