]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix FESeries with non-primitive FEs. 11830/head
authorCe Qin <qince168@gmail.com>
Sun, 21 Feb 2021 02:44:59 +0000 (10:44 +0800)
committerCe Qin <qince168@gmail.com>
Thu, 11 Mar 2021 03:11:06 +0000 (11:11 +0800)
doc/news/changes/minor/20210301CeQin [new file with mode: 0644]
include/deal.II/fe/fe_series.h
include/deal.II/numerics/smoothness_estimator.h
source/fe/fe_series_fourier.cc
source/fe/fe_series_legendre.cc
source/numerics/smoothness_estimator.cc
source/numerics/smoothness_estimator.inst.in
tests/numerics/smoothness_estimator_03.cc [new file with mode: 0644]
tests/numerics/smoothness_estimator_03.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20210301CeQin b/doc/news/changes/minor/20210301CeQin
new file mode 100644 (file)
index 0000000..b7e5050
--- /dev/null
@@ -0,0 +1,4 @@
+New: Add a component parameter to FESeries::Fourier/Legendre to
+make them working with non-primitive elements.
+<br>
+(Ce Qin, 2021/03/01)
index 855720ca5e113928886a0d934783d1c4278e132e..34bf45a52017a733f8643052c901d6577207b2c5 100644 (file)
@@ -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<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
@@ -204,13 +216,20 @@ namespace FESeries
      * 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[
@@ -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<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
@@ -364,6 +394,12 @@ namespace FESeries
      * 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;
   };
 
 
index b0d11b6531ec89342d58280866b3a939f8740b90..77434ff39da449422cf74d9684cfbd893afd03fb 100644 (file)
@@ -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 <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
 
 
@@ -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 <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
 
index 5fb7aede3b151000f21cb1315cec8990f3218bd8..9c3e54ddead31714ea2a70aa67ac94afc81aebf2 100644 (file)
@@ -65,14 +65,16 @@ namespace
   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;
   }
@@ -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<FullMatrix<std::complex<double>>> &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<FullMatrix<std::complex<double>>> &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<FullMatrix<std::complex<double>>> &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<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());
@@ -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<CoefficientType> &matrix =
index f74f217f3b1ff300a2737e9fd3657eaf6c7524d6..ffc4cc4283258b18ba40310d8970e3198a8e340c 100644 (file)
@@ -74,14 +74,16 @@ namespace
   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);
   }
@@ -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<FullMatrix<double>> &    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<FullMatrix<double>> &    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<FullMatrix<double>> &    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<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(),
@@ -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<CoefficientType> &matrix =
index fa7c0c3c41fcac4af9db8a4e56e4e65036395330..91c8c34bcb8914bdffde3ba32265b8958254028b 100644 (file)
@@ -287,7 +287,8 @@ namespace SmoothnessEstimator
 
     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.
       //
@@ -323,7 +324,8 @@ namespace SmoothnessEstimator
 
       return FESeries::Legendre<dim, spacedim>(n_coefficients_per_direction,
                                                fe_collection,
-                                               q_collection);
+                                               q_collection,
+                                               component);
     }
   } // namespace Legendre
 
@@ -574,7 +576,8 @@ namespace SmoothnessEstimator
 
     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.
       //
@@ -613,7 +616,8 @@ namespace SmoothnessEstimator
 
       return FESeries::Fourier<dim, spacedim>(n_coefficients_per_direction,
                                               fe_collection,
-                                              q_collection);
+                                              q_collection,
+                                              component);
     }
   } // namespace Fourier
 } // namespace SmoothnessEstimator
index f33653ce6291f942eee162a5bd6ba1002d1fbed5..9f2b9112a3872ba8a63bbdecd75f138fecd63160 100644 (file)
@@ -67,11 +67,13 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
     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
   }
diff --git a/tests/numerics/smoothness_estimator_03.cc b/tests/numerics/smoothness_estimator_03.cc
new file mode 100644 (file)
index 0000000..7422150
--- /dev/null
@@ -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 <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;
+}
diff --git a/tests/numerics/smoothness_estimator_03.output b/tests/numerics/smoothness_estimator_03.output
new file mode 100644 (file)
index 0000000..67a2ddb
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Ok
+DEAL::Ok

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.