]> https://gitweb.dealii.org/ - dealii.git/commitdiff
implement VectorTools::add_constant 14813/head
authorTimo Heister <timo.heister@gmail.com>
Fri, 17 Feb 2023 20:35:29 +0000 (15:35 -0500)
committerTimo Heister <timo.heister@gmail.com>
Sat, 18 Feb 2023 13:07:55 +0000 (08:07 -0500)
doc/news/changes/minor/20230217tjhei [new file with mode: 0644]
include/deal.II/numerics/vector_tools_mean_value.h
include/deal.II/numerics/vector_tools_mean_value.templates.h
source/numerics/vector_tools_mean_value.inst.in
tests/vector_tools/add_constant_01.cc [new file with mode: 0644]
tests/vector_tools/add_constant_01.output [new file with mode: 0644]
tests/vector_tools/add_constant_02.cc [new file with mode: 0644]
tests/vector_tools/add_constant_02.mpirun=2.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20230217tjhei b/doc/news/changes/minor/20230217tjhei
new file mode 100644 (file)
index 0000000..9ad1a0b
--- /dev/null
@@ -0,0 +1,4 @@
+New: Implement `VectorTools::add_constant` that adds a constant
+to a given component of a finite element solution.
+<br>
+(Timo Heister, 2023/02/17)
index 2a0646be2c96483135f791c65ece50cfe426c15a..6e0246832c8dde499c1bd07ba51e8c3f6d1693ea 100644 (file)
@@ -100,6 +100,38 @@ namespace VectorTools
   void
   subtract_mean_value(VectorType &v, const std::vector<bool> &p_select = {});
 
+  /**
+   * Add the constant @p constant_adjustment to the specified
+   * component of the finite element function given by the coefficient
+   * vector @p solution defined by the given DoFHandler.
+   *
+   * This operation is a common operation to compute a solution with
+   * mean pressure zero for a Stokes flow problem. Here, one can
+   * use VectorTools::compute_mean_value() to compute the value to
+   * subtract.
+   *
+   * For a nodal finite element like FE_Q, this function will simply add
+   * the value @p constant_adjustment to each coefficient of the corresponding
+   * component. If you have the component in a separate block @p b, you
+   * could directly use
+   * <code>solution.block(b) += constant_adjustment</code>
+   * instead of calling this function (and this would be more efficient).
+   * For other finite element spaces like FE_DGP, the logic is more
+   * complicated and handled correctly by this function.
+   *
+   * @note Not all kinds of finite elements are supported and the selected
+   * component must not be part of a non-primitive element for this
+   * implementation to work.
+   *
+   * @note In contrast to subtract_mean_value(), this function can
+   * adjust a single component of a distributed vector.
+   */
+  template <class VectorType, int dim, int spacedim = dim>
+  void
+  add_constant(VectorType &                          solution,
+               const DoFHandler<dim, spacedim> &     dof_handler,
+               const unsigned int                    component,
+               const typename VectorType::value_type constant_adjustment);
 
   /**
    * Compute the mean value of one component of the solution.
index 5d360f60394e101d576599ace3dac8d06a0cbb8a..67d451e768c16ca72fb6980cdbd743fc123ff249 100644 (file)
 
 #include <deal.II/distributed/tria_base.h>
 
+#include <deal.II/fe/fe_dgp.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
 #include <deal.II/fe/fe_values.h>
 
 #include <deal.II/grid/filtered_iterator.h>
@@ -124,6 +128,151 @@ namespace VectorTools
     }
   } // namespace internal
 
+
+
+  template <class VectorType, int dim, int spacedim>
+  void
+  add_constant(VectorType &                          solution,
+               const DoFHandler<dim, spacedim> &     dof_handler,
+               const unsigned int                    component,
+               const typename VectorType::value_type constant_adjustment)
+  {
+    Assert(dof_handler.has_hp_capabilities() == false, ExcNotImplemented());
+
+    AssertDimension(solution.size(), dof_handler.n_dofs());
+    AssertIndexRange(component, dof_handler.get_fe().n_components());
+
+    const FiniteElement<dim, spacedim> &fe_system = dof_handler.get_fe();
+    const FiniteElement<dim, spacedim> &fe = fe_system.get_sub_fe(component, 1);
+
+    if ((dynamic_cast<const FE_DGP<dim, spacedim> *>(&fe) != nullptr))
+      {
+        // The FE to modify is an FE_DGP, which is not a nodal
+        // element. The first shape function of a DGP element happens
+        // to be the constant function, so we just have to adjust the
+        // corresponding DoF on each cell:
+        std::vector<types::global_dof_index> local_dof_indices(
+          dof_handler.get_fe().dofs_per_cell);
+
+        for (const auto &cell : dof_handler.active_cell_iterators())
+          if (cell->is_locally_owned())
+            {
+              cell->get_dof_indices(local_dof_indices);
+              const unsigned int first_pressure_dof =
+                fe_system.component_to_system_index(component, 0);
+
+              // Make sure that this DoF is really owned by the
+              // current processor:
+              Assert(dof_handler.locally_owned_dofs().is_element(
+                       local_dof_indices[first_pressure_dof]),
+                     ExcInternalError());
+
+              // Then adjust its value:
+              solution(local_dof_indices[first_pressure_dof]) +=
+                constant_adjustment;
+            }
+
+        solution.compress(VectorOperation::add);
+      }
+    else if ((dynamic_cast<const FE_Q<dim, spacedim> *>(&fe) != nullptr) ||
+             (dynamic_cast<const FE_SimplexP<dim, spacedim> *>(&fe) != nullptr))
+      {
+        // We need to make sure to not touch DoFs shared between cells more
+        // than once. Instead of counting or limiting the number of times
+        // we touch an individual vector entry, we make a copy of the vector
+        // and use that as the input and add the constant to it. This way
+        // it does not matter how often an individual entry is touched.
+
+        VectorType copy(solution);
+        copy = solution;
+
+        std::vector<types::global_dof_index> local_dof_indices(
+          dof_handler.get_fe().dofs_per_cell);
+
+        for (const auto &cell : dof_handler.active_cell_iterators())
+          if (cell->is_locally_owned())
+            {
+              cell->get_dof_indices(local_dof_indices);
+              for (unsigned i = 0; i < dof_handler.get_fe().dofs_per_cell; ++i)
+                {
+                  if (!fe_system.is_primitive(i))
+                    continue;
+
+                  const auto component_and_index =
+                    fe_system.system_to_component_index(i);
+                  if (component_and_index.first == component)
+                    {
+                      const types::global_dof_index idx = local_dof_indices[i];
+                      // Make sure that this DoF is really owned by the
+                      // current processor:
+                      if (!dof_handler.locally_owned_dofs().is_element(idx))
+                        continue;
+
+                      // Then adjust its value:
+                      solution(idx) = copy(idx) + constant_adjustment;
+                    }
+                }
+            }
+
+        solution.compress(VectorOperation::insert);
+      }
+    else if ((dynamic_cast<const FE_DGQ<dim, spacedim> *>(&fe) != nullptr) ||
+             (dynamic_cast<const FE_SimplexDGP<dim, spacedim> *>(&fe) !=
+              nullptr))
+      {
+        // Add the constant to every single shape function per cell
+
+        std::vector<types::global_dof_index> local_dof_indices(
+          dof_handler.get_fe().dofs_per_cell);
+
+        for (const auto &cell : dof_handler.active_cell_iterators())
+          if (cell->is_locally_owned())
+            {
+              cell->get_dof_indices(local_dof_indices);
+              for (unsigned i = 0; i < dof_handler.get_fe().dofs_per_cell; ++i)
+                {
+                  if (!fe_system.is_primitive(i))
+                    continue;
+
+                  const auto component_and_index =
+                    fe_system.system_to_component_index(i);
+                  if (component_and_index.first == component)
+                    {
+                      // Make sure that this DoF is really owned by the
+                      // current processor:
+                      Assert(dof_handler.locally_owned_dofs().is_element(
+                               local_dof_indices[i]),
+                             ExcInternalError());
+
+                      // Then adjust its value:
+                      solution(local_dof_indices[i]) += constant_adjustment;
+                    }
+                }
+            }
+
+        solution.compress(VectorOperation::add);
+      }
+    else
+      AssertThrow(false, ExcNotImplemented());
+  }
+
+
+
+#ifdef DEAL_II_WITH_TRILINOS
+  template <int dim, int spacedim>
+  void
+  add_constant(LinearAlgebra::EpetraWrappers::Vector &,
+               const DoFHandler<dim, spacedim> &,
+               const unsigned int,
+               const double)
+  {
+    // TODO: no vector access using operator()
+    AssertThrow(false, ExcNotImplemented());
+  }
+#endif
+
+
+
   template <int dim, typename VectorType, int spacedim>
   typename VectorType::value_type
   compute_mean_value(
index 246a32d601ee26439cb11fe2b881f5e6dd1a14b0..0085f615a9245fc4e07d297df6eb390fb42dbeaa 100644 (file)
@@ -36,6 +36,12 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS;
         const VEC &,
         const unsigned int);
 
+      template void
+      add_constant<>(VEC &                                      solution,
+                     const DoFHandler<deal_II_dimension,
+                                      deal_II_space_dimension> &dof_handler,
+                     const unsigned int                         component,
+                     const typename VEC::value_type constant_adjustment);
     \}
 #endif
   }
diff --git a/tests/vector_tools/add_constant_01.cc b/tests/vector_tools/add_constant_01.cc
new file mode 100644 (file)
index 0000000..8519215
--- /dev/null
@@ -0,0 +1,176 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2003 - 2020 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 VectorTools::add_constant()
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_dgp.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_raviart_thomas.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+
+// 0 (dim components)
+// 1
+// 2
+// 3
+template <int dim>
+class Reference : public Function<dim>
+{
+public:
+  Reference()
+    : Function<dim>(dim + 3)
+  {}
+
+  double
+  value(const Point<dim> &p, const unsigned int c) const
+  {
+    if (c == dim)
+      return 1.0;
+    if (c == dim + 1)
+      return 2.0;
+    if (c == dim + 2)
+      return 3.0;
+    return 0.0;
+  }
+};
+
+
+
+template <int dim, class VectorType>
+void
+test()
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(2);
+
+  FESystem<dim>   fe(FE_RaviartThomas<dim>(1),
+                   1,
+                   FE_Q<dim>(2),
+                   1,
+                   FE_DGQ<dim>(2),
+                   1,
+                   FE_DGP<dim>(1),
+                   1);
+  DoFHandler<dim> dofhandler(tria);
+  dofhandler.distribute_dofs(fe);
+
+  VectorType vec(dofhandler.n_dofs());
+  vec = 0.0;
+  for (unsigned int c = dim; c < dim + 3; ++c)
+    VectorTools::add_constant(vec, dofhandler, c, 1.0 * (c - dim + 1));
+
+  Vector<double> cellwise_errors(tria.n_active_cells());
+  QIterated<dim> quadrature(QTrapezoid<1>(), 5);
+
+  const dealii::Function<dim, double> *w = nullptr;
+  VectorTools::integrate_difference(dofhandler,
+                                    vec,
+                                    Reference<dim>(),
+                                    cellwise_errors,
+                                    quadrature,
+                                    VectorTools::L2_norm);
+
+  const double error = VectorTools::compute_global_error(tria,
+                                                         cellwise_errors,
+                                                         VectorTools::L2_norm);
+
+  AssertThrow(error < 1e-10, ExcMessage("Error in integrate_difference"));
+}
+
+template <int dim, class VectorType>
+void
+test_simplex()
+{
+  Triangulation<dim> tria;
+  GridGenerator::subdivided_hyper_cube_with_simplices(tria, 4);
+
+  FESystem<dim>   fe(FE_SimplexP<dim>(1),
+                   dim,
+                   FE_SimplexP<dim>(2),
+                   1,
+                   FE_SimplexP<dim>(1),
+                   1,
+                   FE_SimplexDGP<dim>(1),
+                   1);
+  DoFHandler<dim> dofhandler(tria);
+  dofhandler.distribute_dofs(fe);
+
+  VectorType vec(dofhandler.n_dofs());
+  vec = 0.0;
+  for (unsigned int c = dim; c < dim + 3; ++c)
+    VectorTools::add_constant(vec, dofhandler, c, 1.0 * (c - dim + 1));
+
+  Vector<double>     cellwise_errors(tria.n_active_cells());
+  QGaussSimplex<dim> quadrature(2);
+
+  const dealii::Function<dim, double> *w = nullptr;
+  MappingFE<dim>                       mapping(FE_SimplexP<dim>(1));
+  VectorTools::integrate_difference(mapping,
+                                    dofhandler,
+                                    vec,
+                                    Reference<dim>(),
+                                    cellwise_errors,
+                                    quadrature,
+                                    VectorTools::L2_norm);
+
+  const double error = VectorTools::compute_global_error(tria,
+                                                         cellwise_errors,
+                                                         VectorTools::L2_norm);
+
+  AssertThrow(error < 1e-10, ExcMessage("Error in integrate_difference"));
+}
+
+
+template <int dim>
+void
+big()
+{
+  test<dim, Vector<double>>();
+  test_simplex<dim, Vector<double>>();
+}
+
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+
+  big<2>();
+  big<3>();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/vector_tools/add_constant_01.output b/tests/vector_tools/add_constant_01.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/vector_tools/add_constant_02.cc b/tests/vector_tools/add_constant_02.cc
new file mode 100644 (file)
index 0000000..40b2e11
--- /dev/null
@@ -0,0 +1,203 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2003 - 2020 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 VectorTools::add_constant() in parallel
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/distributed/shared_tria.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_dgp.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_raviart_thomas.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+
+#include <deal.II/lac/petsc_vector.h>
+#include <deal.II/lac/trilinos_vector.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+
+// 0 (dim components)
+// 1
+// 2
+// 3
+template <int dim>
+class Reference : public Function<dim>
+{
+public:
+  Reference()
+    : Function<dim>(dim + 3)
+  {}
+
+  double
+  value(const Point<dim> &p, const unsigned int c) const
+  {
+    if (c == dim)
+      return 1.0;
+    if (c == dim + 1)
+      return 2.0;
+    if (c == dim + 2)
+      return 3.0;
+    return 0.0;
+  }
+};
+
+
+
+template <int dim, class VectorType>
+void
+test()
+{
+  parallel::shared::Triangulation<dim> tria(MPI_COMM_WORLD);
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(2);
+
+  FESystem<dim>   fe(FE_RaviartThomas<dim>(1),
+                   1,
+                   FE_Q<dim>(2),
+                   1,
+                   FE_DGQ<dim>(2),
+                   1,
+                   FE_DGP<dim>(1),
+                   1);
+  DoFHandler<dim> dofhandler(tria);
+  dofhandler.distribute_dofs(fe);
+
+  VectorType vec(dofhandler.locally_owned_dofs(), MPI_COMM_WORLD);
+  vec = 0.0;
+
+  for (unsigned int c = dim; c < dim + 3; ++c)
+    VectorTools::add_constant(vec, dofhandler, c, 1.0 * (c - dim + 1));
+
+  Vector<double> cellwise_errors(tria.n_active_cells());
+  QIterated<dim> quadrature(QTrapezoid<1>(), 5);
+
+  const dealii::Function<dim, double> *w = nullptr;
+  IndexSet                             relevant_dofs;
+  DoFTools::extract_locally_relevant_dofs(dofhandler, relevant_dofs);
+  VectorType ghosted(dofhandler.locally_owned_dofs(),
+                     relevant_dofs,
+                     MPI_COMM_WORLD);
+  ghosted = vec;
+  VectorTools::integrate_difference(dofhandler,
+                                    ghosted,
+                                    Reference<dim>(),
+                                    cellwise_errors,
+                                    quadrature,
+                                    VectorTools::L2_norm);
+
+  const double error = VectorTools::compute_global_error(tria,
+                                                         cellwise_errors,
+                                                         VectorTools::L2_norm);
+
+  AssertThrow(error < 1e-10, ExcMessage("Error in integrate_difference"));
+}
+
+template <int dim, class VectorType>
+void
+test_simplex()
+{
+  parallel::shared::Triangulation<dim> tria(MPI_COMM_WORLD);
+  GridGenerator::subdivided_hyper_cube_with_simplices(tria, 4);
+
+  FESystem<dim>   fe(FE_SimplexP<dim>(1),
+                   dim,
+                   FE_SimplexP<dim>(1),
+                   1,
+                   FE_SimplexP<dim>(1),
+                   1,
+                   FE_SimplexDGP<dim>(1),
+                   1);
+  DoFHandler<dim> dofhandler(tria);
+  dofhandler.distribute_dofs(fe);
+
+  VectorType vec(dofhandler.locally_owned_dofs(), MPI_COMM_WORLD);
+  vec = 0.0;
+  for (unsigned int c = dim; c < dim + 3; ++c)
+    VectorTools::add_constant(vec, dofhandler, c, 1.0 * (c - dim + 1));
+
+  Vector<double>     cellwise_errors(tria.n_active_cells());
+  QGaussSimplex<dim> quadrature(2);
+
+  const dealii::Function<dim, double> *w = nullptr;
+  MappingFE<dim>                       mapping(FE_SimplexP<dim>(1));
+
+  IndexSet relevant_dofs;
+  DoFTools::extract_locally_relevant_dofs(dofhandler, relevant_dofs);
+  VectorType ghosted(dofhandler.locally_owned_dofs(),
+                     relevant_dofs,
+                     MPI_COMM_WORLD);
+  ghosted = vec;
+  VectorTools::integrate_difference(mapping,
+                                    dofhandler,
+                                    ghosted,
+                                    Reference<dim>(),
+                                    cellwise_errors,
+                                    quadrature,
+                                    VectorTools::L2_norm);
+
+  const double error = VectorTools::compute_global_error(tria,
+                                                         cellwise_errors,
+                                                         VectorTools::L2_norm);
+
+  AssertThrow(error < 1e-10, ExcMessage("Error in integrate_difference"));
+}
+
+
+template <int dim>
+void
+big()
+{
+  test<dim, LinearAlgebra::distributed::Vector<double>>();
+  test_simplex<dim, LinearAlgebra::distributed::Vector<double>>();
+
+#ifdef DEAL_II_WITH_TRILINOS
+  test<dim, TrilinosWrappers::MPI::Vector>();
+  test_simplex<dim, TrilinosWrappers::MPI::Vector>();
+#endif
+#ifdef DEAL_II_WITH_PETSC
+  test<dim, PETScWrappers::MPI::Vector>();
+  test_simplex<dim, PETScWrappers::MPI::Vector>();
+#endif
+}
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  MPILogInitAll                    log;
+
+  big<2>();
+  big<3>();
+
+  deallog << "OK" << std::endl;
+}
diff --git a/tests/vector_tools/add_constant_02.mpirun=2.output b/tests/vector_tools/add_constant_02.mpirun=2.output
new file mode 100644 (file)
index 0000000..5210393
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL:0::OK
+
+DEAL:1::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.