FEValuesViews::ReorderedView

dealii · Aug 7, 2023 · 7a86faf · 7a86faf
1 parent b30869c
commit 7a86faf
Show file tree

Hide file tree

Showing 9 changed files with 1,146 additions and 0 deletions.
diff --git a/doc/news/changes/minor/20230802LucaHeltai b/doc/news/changes/minor/20230802LucaHeltai
@@ -0,0 +1,4 @@
+New: The new class FEValuesViews::ReorderedView allows one to filter an existing
+FEValuesViews object via two renumbering vectors, one acting on the degrees of
+freedom, and the other acting on the quadrature points. 
+<br> (Luca Heltai, 2023/08/02)
diff --git a/include/deal.II/fe/fe_coupling_values.h b/include/deal.II/fe/fe_coupling_values.h
diff --git a/include/deal.II/fe/fe_values_views.h b/include/deal.II/fe/fe_values_views.h
@@ -52,6 +52,11 @@ DEAL_II_NAMESPACE_OPEN
 #ifndef DOXYGEN
 template <int dim, int spacedim = dim>
 class FEValuesBase;
+namespace FEValuesViews
+{
+  template <typename ViewType>
+  class ReorderedView;
+}
 #endif
 
 namespace internal
@@ -575,6 +580,12 @@ namespace FEValuesViews
      * Store the data about shape functions.
      */
     std::vector<ShapeFunctionData> shape_function_data;
+
+    /**
+     * Make sure that ReorderedView can access the private data of this class.
+     */
+    template <typename>
+    friend class ReorderedView;
   };
 
 
@@ -1298,6 +1309,12 @@ namespace FEValuesViews
      * Store the data about shape functions.
      */
     std::vector<ShapeFunctionData> shape_function_data;
+
+    /**
+     * Make sure that ReorderedView can access the private data of this class.
+     */
+    template <typename>
+    friend class ReorderedView;
   };
 
 
@@ -1611,6 +1628,12 @@ namespace FEValuesViews
      * Store the data about shape functions.
      */
     std::vector<ShapeFunctionData> shape_function_data;
+
+    /**
+     * Make sure that ReorderedView can access the private data of this class.
+     */
+    template <typename>
+    friend class ReorderedView;
   };
 
 
@@ -1985,6 +2008,12 @@ namespace FEValuesViews
      * Store the data about shape functions.
      */
     std::vector<ShapeFunctionData> shape_function_data;
+
+    /**
+     * Make sure that ReorderedView can access the private data of this class.
+     */
+    template <typename>
+    friend class ReorderedView;
   };
 
 } // namespace FEValuesViews

diff --git a/tests/fe/fe_values_views_reordered_01.cc b/tests/fe/fe_values_views_reordered_01.cc
@@ -0,0 +1,138 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2023 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Test renumbering of FEValuesViews::Scalar using
+// FEValuesViews::ReorderedView
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_coupling_values.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+#include "../test_grids.h"
+
+template <int dim>
+void
+test()
+{
+  Triangulation<dim> tria;
+  TestGrids::hyper_line(tria, 2);
+
+  DoFHandler<dim> dofh(tria);
+  FE_Q<dim>       fe(1);
+  dofh.distribute_dofs(fe);
+
+  MappingQ<dim> mapping(1);
+  UpdateFlags   update_flags =
+    update_values | update_quadrature_points | update_JxW_values;
+
+  FEValues<dim> fv1(mapping, fe, QGauss<dim>(fe.degree + 1), update_flags);
+  fv1.reinit(dofh.begin_active());
+
+  FEValuesViews::Scalar<dim> fv1_scalar(fv1, 0);
+
+  using iota = std_cxx20::ranges::iota_view<unsigned int, unsigned int>;
+  {
+    // No renumbering
+    for (const auto &i : fv1.dof_indices())
+      {
+        deallog << "i = " << i << ", fv1_scalar = " << fv1_scalar.value(i, 0)
+                << std::endl;
+      }
+  }
+  {
+    // Trivial renumbering
+    auto                           id = iota(0, fv1.dofs_per_cell);
+    std::vector<unsigned int>      renumber(id.begin(), id.end());
+    FEValuesViews::RenumberingData data(renumber);
+    FEValuesViews::ReorderedView<FEValuesViews::Scalar<dim>>
+      fv1_scalar_reordered(fv1_scalar, data);
+
+    for (const auto &i : fv1.dof_indices())
+      {
+        deallog << "i = " << i << ", fv1_scalar_trivial_renumber = "
+                << fv1_scalar_reordered.value(i, 0) << std::endl;
+      }
+  }
+  {
+    // Inverse renumbering
+    auto                      id = iota(0, fv1.dofs_per_cell);
+    std::vector<unsigned int> renumber(id.begin(), id.end());
+    std::reverse(renumber.begin(), renumber.end());
+    FEValuesViews::RenumberingData data(renumber);
+    FEValuesViews::ReorderedView<FEValuesViews::Scalar<dim>>
+      fv1_scalar_reordered(fv1_scalar, data);
+    for (const auto &i : fv1.dof_indices())
+      {
+        deallog << "i = " << i << ", fv1_scalar_inverse_renumber = "
+                << fv1_scalar_reordered.value(i, 0) << std::endl;
+      }
+  }
+  // Now test quadrature renumbering
+  {
+    // No renumbering
+    for (const auto &q : fv1.quadrature_point_indices())
+      {
+        deallog << "q = " << q
+                << ", fv1_scalar_0_q = " << fv1_scalar.value(0, q) << std::endl;
+      }
+  }
+  {
+    // Trivial renumbering
+    auto                           id = iota(0, fv1.n_quadrature_points);
+    std::vector<unsigned int>      renumber(id.begin(), id.end());
+    FEValuesViews::RenumberingData data({}, renumber);
+    FEValuesViews::ReorderedView<FEValuesViews::Scalar<dim>>
+      fv1_scalar_reordered(fv1_scalar, data);
+    for (const auto &q : fv1.quadrature_point_indices())
+      {
+        deallog << "q = " << q << ", fv1_scalar_0_q_trivial_renumbering = "
+                << fv1_scalar_reordered.value(0, q) << std::endl;
+      }
+  }
+  {
+    // Inverse renumbering
+    auto                      id = iota(0, fv1.n_quadrature_points);
+    std::vector<unsigned int> renumber(id.begin(), id.end());
+    std::reverse(renumber.begin(), renumber.end());
+    FEValuesViews::RenumberingData data({}, renumber);
+    FEValuesViews::ReorderedView<FEValuesViews::Scalar<dim>>
+      fv1_scalar_reordered(fv1_scalar, data);
+    for (const auto &q : fv1.quadrature_point_indices())
+      {
+        deallog << "q = " << q << ", fv1_scalar_0_q_inverse_renumbering = "
+                << fv1_scalar_reordered.value(0, q) << std::endl;
+      }
+  }
+}
+
+int
+main()
+{
+  initlog();
+  test<2>();
+}
diff --git a/tests/fe/fe_values_views_reordered_01.output b/tests/fe/fe_values_views_reordered_01.output
@@ -0,0 +1,25 @@
+
+DEAL::i = 0, fv1_scalar = 0.622008
+DEAL::i = 1, fv1_scalar = 0.166667
+DEAL::i = 2, fv1_scalar = 0.166667
+DEAL::i = 3, fv1_scalar = 0.0446582
+DEAL::i = 0, fv1_scalar_trivial_renumber = 0.622008
+DEAL::i = 1, fv1_scalar_trivial_renumber = 0.166667
+DEAL::i = 2, fv1_scalar_trivial_renumber = 0.166667
+DEAL::i = 3, fv1_scalar_trivial_renumber = 0.0446582
+DEAL::i = 0, fv1_scalar_inverse_renumber = 0.0446582
+DEAL::i = 1, fv1_scalar_inverse_renumber = 0.166667
+DEAL::i = 2, fv1_scalar_inverse_renumber = 0.166667
+DEAL::i = 3, fv1_scalar_inverse_renumber = 0.622008
+DEAL::q = 0, fv1_scalar_0_q = 0.622008
+DEAL::q = 1, fv1_scalar_0_q = 0.166667
+DEAL::q = 2, fv1_scalar_0_q = 0.166667
+DEAL::q = 3, fv1_scalar_0_q = 0.0446582
+DEAL::q = 0, fv1_scalar_0_q_trivial_renumbering = 0.622008
+DEAL::q = 1, fv1_scalar_0_q_trivial_renumbering = 0.166667
+DEAL::q = 2, fv1_scalar_0_q_trivial_renumbering = 0.166667
+DEAL::q = 3, fv1_scalar_0_q_trivial_renumbering = 0.0446582
+DEAL::q = 0, fv1_scalar_0_q_inverse_renumbering = 0.0446582
+DEAL::q = 1, fv1_scalar_0_q_inverse_renumbering = 0.166667
+DEAL::q = 2, fv1_scalar_0_q_inverse_renumbering = 0.166667
+DEAL::q = 3, fv1_scalar_0_q_inverse_renumbering = 0.622008
diff --git a/tests/fe/fe_values_views_reordered_02.cc b/tests/fe/fe_values_views_reordered_02.cc
@@ -0,0 +1,145 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2023 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Test renumbering of FEValuesViews::Scalar using
+// FEValuesViews::ReorderedView. Test get_function_values and
+// get_function_values_from_local_dof_values
+
+#include <deal.II/base/function_parser.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_coupling_values.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+#include "../test_grids.h"
+
+template <int dim>
+void
+test()
+{
+  Triangulation<dim> tria;
+  TestGrids::hyper_line(tria, 2);
+
+  DoFHandler<dim> dofh(tria);
+  FE_Q<dim>       fe(1);
+  dofh.distribute_dofs(fe);
+
+  // Create a non trivial function of x and y
+  FunctionParser<dim> f("x+10*y");
+
+  Vector<double> vec(dofh.n_dofs());
+  VectorTools::interpolate(dofh, f, vec);
+
+  MappingQ<dim> mapping(1);
+  UpdateFlags   update_flags =
+    update_values | update_quadrature_points | update_JxW_values;
+
+  FEValues<dim> fe_v(mapping, fe, QGauss<dim>(fe.degree + 1), update_flags);
+  const auto &  cell = dofh.begin_active();
+  fe_v.reinit(cell);
+
+  const auto &fv1_scalar = fe_v[FEValuesExtractors::Scalar(0)];
+
+  // Small print function
+  auto print = [](const auto &str, const auto &v) {
+    deallog << str << " = ";
+    for (const auto &i : v)
+      deallog << i << " ";
+    deallog << std::endl;
+  };
+
+  // Now build an artificial renumbering vector for dofs and one for quadratures
+  // - The first renumbering, simply shifts the dof indices by 1, as if there
+  //   was an additional dof at the beginning, which is ignored by this view
+  // - The second renumbering, repeats each quadrature point twice, so that the
+  //   returned values are actually repeated.
+  using iota = std_cxx20::ranges::iota_view<unsigned int, unsigned int>;
+
+  std::vector<unsigned int> dof_renumbering(fe_v.dofs_per_cell + 1);
+  std::vector<unsigned int> quad_renumbering(fe_v.n_quadrature_points * 2);
+  {
+    auto id            = iota(0, fe_v.dofs_per_cell);
+    dof_renumbering[0] = numbers::invalid_unsigned_int;
+    std::copy(id.begin(), id.end(), dof_renumbering.begin() + 1);
+  }
+  {
+    auto id = iota(0, fe_v.n_quadrature_points);
+    std::copy(id.begin(), id.end(), quad_renumbering.begin());
+    std::copy(id.begin(), id.end(), quad_renumbering.begin() + id.size());
+  }
+
+  deallog << "dof_renumbering = ";
+  for (const auto &i : dof_renumbering)
+    deallog << (int)i << " ";
+  deallog << std::endl;
+  print("quad_renumbering", quad_renumbering);
+
+  std::vector<double> function_values(fe_v.n_quadrature_points);
+  std::vector<double> function_values_renumbered(fe_v.n_quadrature_points * 2);
+
+  fv1_scalar.get_function_values(vec, function_values);
+  print("function_values", function_values);
+
+  FEValuesViews::RenumberingData data(dof_renumbering, quad_renumbering);
+  FEValuesViews::ReorderedView<FEValuesViews::Scalar<dim>>
+    fv1_scalar_renumbered(fv1_scalar, data);
+
+  fv1_scalar_renumbered.get_function_values(vec, function_values_renumbered);
+  print("function_values_renumbered", function_values_renumbered);
+
+  std::vector<double> dof_values(fe_v.dofs_per_cell);
+  std::vector<double> dof_values_renumbered(fe_v.dofs_per_cell + 1);
+
+  cell->get_dof_values(vec, dof_values.begin(), dof_values.end());
+
+  // Now fill the first dof with a dummy value
+  dof_values_renumbered[0] = 42;
+  std::copy(dof_values.begin(),
+            dof_values.end(),
+            dof_values_renumbered.begin() + 1);
+
+  print("dof_values", dof_values);
+  print("dof_values_renumbered", dof_values_renumbered);
+
+  fv1_scalar.get_function_values_from_local_dof_values(dof_values,
+                                                       function_values);
+  print("function_values from dof_values", function_values);
+
+  fv1_scalar_renumbered.get_function_values_from_local_dof_values(
+    dof_values_renumbered, function_values_renumbered);
+
+  print("function_values_renumbered from dof_values_renumbered",
+        function_values_renumbered);
+}
+
+int
+main()
+{
+  initlog();
+  test<2>();
+}
diff --git a/tests/fe/fe_values_views_reordered_02.output b/tests/fe/fe_values_views_reordered_02.output
@@ -0,0 +1,9 @@
+
+DEAL::dof_renumbering = -1 0 1 2 3 
+DEAL::quad_renumbering = 0 1 2 3 0 1 2 3 
+DEAL::function_values = 2.32457 2.90192 8.09808 8.67543 
+DEAL::function_values_renumbered = 2.32457 2.90192 8.09808 8.67543 2.32457 2.90192 8.09808 8.67543 
+DEAL::dof_values = 0.00000 1.00000 10.0000 11.0000 
+DEAL::dof_values_renumbered = 42.0000 0.00000 1.00000 10.0000 11.0000 
+DEAL::function_values from dof_values = 2.32457 2.90192 8.09808 8.67543 
+DEAL::function_values_renumbered from dof_values_renumbered = 2.32457 2.90192 8.09808 8.67543 2.32457 2.90192 8.09808 8.67543