From: Peter Munch <peterrmuench@gmail.com>
Date: Thu, 21 Mar 2024 09:15:43 +0000 (+0100)
Subject: Introduce FEPointEvaluation::integrate_value()
X-Git-Tag: v9.6.0-rc1~465^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F16425%2Fhead;p=dealii.git

Introduce FEPointEvaluation::integrate_value()
---

diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h
index a8247d1a44..882d8839b1 100644
--- a/include/deal.II/matrix_free/fe_evaluation.h
+++ b/include/deal.II/matrix_free/fe_evaluation.h
@@ -558,8 +558,9 @@ public:
 #endif
 
   /**
-   * Takes values at quadrature points, multiplies by the Jacobian determinant
-   * and quadrature weights (JxW) and sums the values for all quadrature
+   * Take values collected at quadrature points via the submit_value() function,
+   * multiply by the Jacobian determinant
+   * and quadrature weights (JxW) and sum the values for all quadrature
    * points on the cell. The result is a scalar, representing the integral
    * over the function over the cell. If a vector-element is used, the
    * resulting components are still separated. Moreover, if vectorization is
diff --git a/include/deal.II/matrix_free/fe_point_evaluation.h b/include/deal.II/matrix_free/fe_point_evaluation.h
index e7f1426163..cae51b5b7c 100644
--- a/include/deal.II/matrix_free/fe_point_evaluation.h
+++ b/include/deal.II/matrix_free/fe_point_evaluation.h
@@ -85,6 +85,23 @@ namespace internal
         result[component] = vector_entry;
       }
 
+      static scalar_value_type
+      sum_value(const scalar_value_type &result)
+      {
+        return result;
+      }
+
+      static scalar_value_type
+      sum_value(const vectorized_value_type &result)
+      {
+        scalar_value_type result_scalar = {};
+
+        for (unsigned int c = 0; c < n_components; ++c)
+          result_scalar[c] = result[c].sum();
+
+        return result_scalar;
+      }
+
       static ScalarNumber
       sum_value(const unsigned int           component,
                 const vectorized_value_type &result)
@@ -236,6 +253,18 @@ namespace internal
         result = vector_entry;
       }
 
+      static scalar_value_type
+      sum_value(const scalar_value_type &result)
+      {
+        return result;
+      }
+
+      static scalar_value_type
+      sum_value(const vectorized_value_type &result)
+      {
+        return result.sum();
+      }
+
       static ScalarNumber
       sum_value(const unsigned int, const vectorized_value_type &result)
       {
@@ -390,6 +419,23 @@ namespace internal
         result[component] = vector_entry;
       }
 
+      static scalar_value_type
+      sum_value(const scalar_value_type &result)
+      {
+        return result;
+      }
+
+      static scalar_value_type
+      sum_value(const vectorized_value_type &result)
+      {
+        scalar_value_type result_scalar = {};
+
+        for (unsigned int c = 0; c < dim; ++c)
+          result_scalar[c] = result[c].sum();
+
+        return result_scalar;
+      }
+
       static ScalarNumber
       sum_value(const unsigned int           component,
                 const vectorized_value_type &result)
@@ -541,6 +587,18 @@ namespace internal
         result = vector_entry;
       }
 
+      static scalar_value_type
+      sum_value(const scalar_value_type &result)
+      {
+        return result;
+      }
+
+      static scalar_value_type
+      sum_value(const vectorized_value_type &result)
+      {
+        return result.sum();
+      }
+
       static ScalarNumber
       sum_value(const unsigned int, const vectorized_value_type &result)
       {
@@ -859,6 +917,16 @@ public:
   Point<dim, Number>
   unit_point(const unsigned int point_index) const;
 
+  /**
+   * Take values collected at quadrature points via the submit_value()
+    * function, multiply by the Jacobian determinant
+    * and quadrature weights (JxW) and sum the values for all quadrature
+   * points on the cell. The result is a scalar, representing the integral
+   * of the function over the cell.
+   */
+  scalar_value_type
+  integrate_value() const;
+
   /**
    * Return an object that can be thought of as an array containing all indices
    * from zero to n_quadrature_points. This allows to write code using
@@ -2466,6 +2534,22 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate(
 
 
 
+template <int n_components_, int dim, int spacedim, typename Number>
+inline typename FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
+  scalar_value_type
+  FEPointEvaluationBase<n_components_, dim, spacedim, Number>::integrate_value()
+    const
+{
+  value_type return_value = {};
+
+  for (const auto point_index : this->quadrature_point_indices())
+    return_value += values[point_index] * this->JxW(point_index);
+
+  return ETT::sum_value(return_value);
+}
+
+
+
 template <int n_components_, int dim, int spacedim, typename Number>
 template <std::size_t stride_view>
 void
diff --git a/tests/matrix_free/integrate_value_01.cc b/tests/matrix_free/integrate_value_01.cc
new file mode 100644
index 0000000000..0dfee3258c
--- /dev/null
+++ b/tests/matrix_free/integrate_value_01.cc
@@ -0,0 +1,313 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2024 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 FEEvaluation::integrate_value() and
+// FEPointEvaluation::integrate_value().
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+template <int dim, typename Number>
+class RightHandSideFunction : public Function<dim, Number>
+{
+public:
+  RightHandSideFunction(const unsigned int n_components)
+    : Function<dim, Number>(n_components)
+  {}
+
+  Number
+  value(const Point<dim, Number> &p, const unsigned int component) const
+  {
+    if (component == 0)
+      return p[0];
+    else if (component == 1)
+      return p[1] * p[1];
+    else
+      return 0;
+  }
+
+  VectorizedArray<Number>
+  value(const Point<dim, VectorizedArray<Number>> &p,
+        const unsigned int                         component) const
+  {
+    if (component == 0)
+      return p[0];
+    else if (component == 1)
+      return p[1] * p[1];
+    else
+      return 0;
+  }
+};
+
+
+template <typename Number>
+Number
+sum(const Number &value)
+{
+  return value;
+}
+
+
+template <typename Number, int dim>
+Number
+sum(const Tensor<1, dim, Number> &value)
+{
+  Number result = {};
+
+  for (unsigned int i = 0; i < dim; ++i)
+    result += value[i];
+
+  return result;
+}
+
+
+template <typename Number>
+Number &
+get(Number &value, const unsigned int)
+{
+  return value;
+}
+
+
+template <typename Number, int dim>
+Number &
+get(Tensor<1, dim, Number> &value, const unsigned int component)
+{
+  return value[component];
+}
+
+
+
+template <int dim, int n_components>
+void
+test()
+{
+  const unsigned int fe_degree                 = 2;
+  const unsigned int n_global_mesh_refinements = 2;
+
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(n_global_mesh_refinements);
+
+  FESystem<dim> fe(FE_Q<dim>(fe_degree), n_components);
+  QGauss<dim>   quad(fe_degree + 1);
+
+  MappingQ1<dim> mapping;
+
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  Vector<double> vector(dof_handler.n_dofs());
+
+  RightHandSideFunction<dim, double> fu(n_components);
+
+  Vector<double> integrals(tria.n_active_cells());
+
+  // create reference solution with VectorTools::integrate_difference()
+  {
+    Functions::ConstantFunction<dim, double> zero(0.0, n_components);
+    VectorTools::integrate_difference(mapping,
+                                      dof_handler,
+                                      vector,
+                                      fu,
+                                      integrals,
+                                      quad,
+                                      VectorTools::NormType::L2_norm);
+
+    for (unsigned int i = 0; i < integrals.size(); ++i)
+      integrals[i] = integrals[i] * integrals[i];
+  }
+
+  // create solution with FEEvaluation::integrate_value()
+  {
+    Vector<double> my_integrals(tria.n_active_cells());
+
+    typename MatrixFree<dim>::AdditionalData ad;
+    ad.mapping_update_flags =
+      update_JxW_values | update_values | update_quadrature_points;
+    MatrixFree<dim> matrix_free;
+
+    matrix_free.reinit(
+      mapping, dof_handler, AffineConstraints<double>(), quad, ad);
+
+    FEEvaluation<dim, -1, 0, n_components> fe_eval(matrix_free);
+
+    for (unsigned int cell = 0; cell < matrix_free.n_cell_batches(); ++cell)
+      {
+        fe_eval.reinit(cell);
+
+        for (const auto q : fe_eval.quadrature_point_indices())
+          {
+            const auto point = fe_eval.quadrature_point(q);
+
+            typename FEEvaluation<dim, -1, 0, n_components>::value_type value;
+
+            for (unsigned int c = 0; c < n_components; ++c)
+              get(value, c) = fu.value(point, c) * fu.value(point, c);
+
+            fe_eval.submit_value(value, q);
+          }
+
+        const auto integrated_cell_value = sum(fe_eval.integrate_value());
+
+        for (unsigned int v = 0;
+             v < matrix_free.n_active_entries_per_cell_batch(cell);
+             ++v)
+          my_integrals[matrix_free.get_cell_iterator(cell, v)
+                         ->active_cell_index()] = integrated_cell_value[v];
+      }
+
+    for (unsigned int i = 0; i < my_integrals.size(); ++i)
+      Assert(std::abs(my_integrals[i] - integrals[i]) < 1e-12,
+             ExcInternalError());
+  }
+
+  // create solution with FEEvaluation::integrate_value() (non-vectorized)
+  {
+    Vector<double> my_integrals(tria.n_active_cells());
+
+    NonMatching::MappingInfo<dim, dim, double> mapping_info(mapping,
+                                                            update_JxW_values |
+                                                              update_values);
+
+    std::vector<typename Triangulation<dim>::cell_iterator> cells;
+    std::vector<Quadrature<dim>>                            quadratures;
+
+    for (const auto &cell : tria.active_cell_iterators())
+      {
+        cells.push_back(cell);
+        quadratures.push_back(quad);
+      }
+
+    mapping_info.reinit_cells(cells, quadratures);
+
+    FEPointEvaluation<n_components, dim, dim, double> fe_eval(mapping_info, fe);
+
+    for (unsigned int cell = 0; cell < cells.size(); ++cell)
+      {
+        fe_eval.reinit(cell);
+
+        for (const auto q : fe_eval.quadrature_point_indices())
+          {
+            const auto point = fe_eval.real_point(q);
+
+            typename FEPointEvaluation<n_components, dim, dim, double>::
+              value_type value;
+
+            for (unsigned int c = 0; c < n_components; ++c)
+              get(value, c) = fu.value(point, c) * fu.value(point, c);
+
+            fe_eval.submit_value(value, q);
+          }
+
+        my_integrals[cell] = sum(fe_eval.integrate_value());
+      }
+
+    for (unsigned int i = 0; i < my_integrals.size(); ++i)
+      Assert(std::abs(my_integrals[i] - integrals[i]) < 1e-12,
+             ExcInternalError());
+  }
+
+  // create solution with FEEvaluation::integrate_value() (vectorized)
+  {
+    Vector<double> my_integrals(tria.n_active_cells());
+
+    NonMatching::MappingInfo<dim, dim, VectorizedArray<double>> mapping_info(
+      mapping, update_JxW_values | update_values);
+
+    std::vector<typename Triangulation<dim>::cell_iterator> cells;
+    std::vector<Quadrature<dim>>                            quadratures;
+
+    for (const auto &cell : tria.active_cell_iterators())
+      {
+        cells.push_back(cell);
+        quadratures.push_back(quad);
+      }
+
+    mapping_info.reinit_cells(cells, quadratures);
+
+    FEPointEvaluation<n_components, dim, dim, VectorizedArray<double>> fe_eval(
+      mapping_info, fe);
+
+    for (unsigned int cell = 0; cell < cells.size(); ++cell)
+      {
+        fe_eval.reinit(cell);
+
+        for (const auto q : fe_eval.quadrature_point_indices())
+          {
+            const auto point = fe_eval.real_point(q);
+
+            typename FEPointEvaluation<n_components,
+                                       dim,
+                                       dim,
+                                       VectorizedArray<double>>::value_type
+              value;
+
+            for (unsigned int c = 0; c < n_components; ++c)
+              get(value, c) = fu.value(point, c) * fu.value(point, c);
+
+            fe_eval.submit_value(value, q);
+          }
+
+        my_integrals[cell] = sum(fe_eval.integrate_value());
+      }
+
+    for (unsigned int i = 0; i < my_integrals.size(); ++i)
+      Assert(std::abs(my_integrals[i] - integrals[i]) < 1e-12,
+             ExcInternalError());
+  }
+
+  deallog << "OK" << std::endl;
+}
+
+
+int
+main()
+{
+  initlog();
+
+  test<2, 1>();
+  test<2, 2>();
+  test<2, 3>();
+
+  test<3, 1>();
+  test<3, 2>();
+  test<3, 3>();
+  test<3, 4>();
+}
diff --git a/tests/matrix_free/integrate_value_01.output b/tests/matrix_free/integrate_value_01.output
new file mode 100644
index 0000000000..5cf798a265
--- /dev/null
+++ b/tests/matrix_free/integrate_value_01.output
@@ -0,0 +1,8 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK