]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement integrate sum into values option 16043/head
authorMaximilian Bergbauer <bergbauer@lnm.mw.tum.de>
Mon, 25 Sep 2023 08:20:57 +0000 (10:20 +0200)
committerMaximilian Bergbauer <maximilian.bergbauer@tum.de>
Wed, 4 Oct 2023 09:24:31 +0000 (11:24 +0200)
doc/news/changes/minor/20231004Bergbauer_2 [new file with mode: 0644]
include/deal.II/matrix_free/fe_point_evaluation.h
tests/matrix_free/point_evaluation_28.cc [new file with mode: 0644]
tests/matrix_free/point_evaluation_28.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20231004Bergbauer_2 b/doc/news/changes/minor/20231004Bergbauer_2
new file mode 100644 (file)
index 0000000..a4e0440
--- /dev/null
@@ -0,0 +1,4 @@
+Added: An option to sum into the solution values vector was added
+to FEPointEvaluation integrate() and test_and_sum().
+<br>
+(Maximilian Bergbauer, 2023/10/04)
index c1e2c28216b3e0e11fbb536e4f55ffa4cbeb9c03..15baefe1f0fa3774f2771e4470744b0eac40dfbe 100644 (file)
@@ -885,12 +885,15 @@ public:
    * @param[in] integration_flags Flags specifying which quantities should be
    * integrated at the points.
    *
+   * @param[in] sum_into_values Flag specifying if the integrated values
+   * should be summed into the solution values. Defaults to false.
+   *
    */
   template <std::size_t stride_view>
   void
   integrate(const StridedArrayView<ScalarNumber, stride_view> &solution_values,
-            const EvaluationFlags::EvaluationFlags &integration_flags);
-
+            const EvaluationFlags::EvaluationFlags &integration_flags,
+            const bool                              sum_into_values = false);
 
   /**
    * This function multiplies the quantities passed in by previous
@@ -909,15 +912,14 @@ public:
    * @param[in] integration_flags Flags specifying which quantities should be
    * integrated at the points.
    *
+   * @param[in] sum_into_values Flag specifying if the integrated values
+   * should be summed into the solution values. Defaults to false.
+   *
    */
   void
   integrate(const ArrayView<ScalarNumber>          &solution_values,
-            const EvaluationFlags::EvaluationFlags &integration_flags)
-  {
-    integrate(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
-                                                solution_values.size()),
-              integration_flags);
-  }
+            const EvaluationFlags::EvaluationFlags &integration_flags,
+            const bool                              sum_into_values = false);
 
   /**
    * This function multiplies the quantities passed in by previous
@@ -940,12 +942,16 @@ public:
    * @param[in] integration_flags Flags specifying which quantities should be
    * integrated at the points.
    *
+   * @param[in] sum_into_values Flag specifying if the integrated values
+   * should be summed into the solution values. Defaults to false.
+   *
    */
   template <std::size_t stride_view>
   void
   test_and_sum(
     const StridedArrayView<ScalarNumber, stride_view> &solution_values,
-    const EvaluationFlags::EvaluationFlags            &integration_flags);
+    const EvaluationFlags::EvaluationFlags            &integration_flags,
+    const bool                                         sum_into_values = false);
 
   /**
    * This function multiplies the quantities passed in by previous
@@ -968,15 +974,14 @@ public:
    * @param[in] integration_flags Flags specifying which quantities should be
    * integrated at the points.
    *
+   * @param[in] sum_into_values Flag specifying if the integrated values
+   * should be summed into the solution values. Defaults to false.
+   *
    */
   void
   test_and_sum(const ArrayView<ScalarNumber>          &solution_values,
-               const EvaluationFlags::EvaluationFlags &integration_flags)
-  {
-    test_and_sum(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
-                                                   solution_values.size()),
-                 integration_flags);
-  }
+               const EvaluationFlags::EvaluationFlags &integration_flags,
+               const bool                              sum_into_values = false);
 
   /**
    * Return the value at quadrature point number @p point_index after a call to
@@ -1172,7 +1177,8 @@ private:
   finish_integrate_fast(
     const StridedArrayView<ScalarNumber, stride_view> &solution_values,
     const EvaluationFlags::EvaluationFlags            &integration_flags,
-    vectorized_value_type *solution_values_vectorized_linear);
+    vectorized_value_type *solution_values_vectorized_linear,
+    const bool             sum_into_values);
 
   /**
    * Fast path of the integrate function.
@@ -1184,7 +1190,8 @@ private:
   void
   integrate_fast(
     const StridedArrayView<ScalarNumber, stride_view> &solution_values,
-    const EvaluationFlags::EvaluationFlags            &integration_flags);
+    const EvaluationFlags::EvaluationFlags            &integration_flags,
+    const bool                                         sum_into_values);
 
   /**
    * Slow path of the integrate function using FEValues.
@@ -1193,7 +1200,8 @@ private:
   void
   integrate_slow(
     const StridedArrayView<ScalarNumber, stride_view> &solution_values,
-    const EvaluationFlags::EvaluationFlags            &integration_flags);
+    const EvaluationFlags::EvaluationFlags            &integration_flags,
+    const bool                                         sum_into_values);
 
   /**
    * Implementation of the integrate/test_and_sum function.
@@ -1202,7 +1210,8 @@ private:
   void
   do_integrate(
     const StridedArrayView<ScalarNumber, stride_view> &solution_values,
-    const EvaluationFlags::EvaluationFlags            &integration_flags);
+    const EvaluationFlags::EvaluationFlags            &integration_flags,
+    const bool                                         sum_into_values);
 
   /**
    * Number of quadrature batches of the current cell/face.
@@ -2393,9 +2402,10 @@ inline void
 FEPointEvaluation<n_components_, dim, spacedim, Number>::finish_integrate_fast(
   const StridedArrayView<ScalarNumber, stride_view> &solution_values,
   const EvaluationFlags::EvaluationFlags            &integration_flags,
-  vectorized_value_type *solution_values_vectorized_linear)
+  vectorized_value_type *solution_values_vectorized_linear,
+  const bool             sum_into_values)
 {
-  if (fe->n_components() > n_components)
+  if (!sum_into_values && fe->n_components() > n_components)
     for (unsigned int i = 0; i < solution_values.size(); ++i)
       solution_values[i] = 0;
 
@@ -2434,13 +2444,19 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::finish_integrate_fast(
           if (is_linear || renumber.empty())
             {
               for (unsigned int i = 0; i < dofs_per_comp; ++i)
-                solution_values[i + offset] = output[i];
+                if (sum_into_values)
+                  solution_values[i + offset] += output[i];
+                else
+                  solution_values[i + offset] = output[i];
             }
           else
             {
               const unsigned int *renumber_ptr = renumber.data() + offset;
               for (unsigned int i = 0; i < dofs_per_comp; ++i)
-                solution_values[renumber_ptr[i]] = output[i];
+                if (sum_into_values)
+                  solution_values[renumber_ptr[i]] += output[i];
+                else
+                  solution_values[renumber_ptr[i]] = output[i];
             }
         }
       else
@@ -2448,18 +2464,29 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::finish_integrate_fast(
           if (is_linear || renumber.empty())
             {
               for (unsigned int i = 0; i < dofs_per_comp; ++i)
-                solution_values[i + offset] =
-                  ETT::sum_value(comp,
-                                 is_linear ?
-                                   *(solution_values_vectorized_linear + i) :
-                                   solution_renumbered_vectorized[i]);
+                if (sum_into_values)
+                  solution_values[i + offset] +=
+                    ETT::sum_value(comp,
+                                   is_linear ?
+                                     *(solution_values_vectorized_linear + i) :
+                                     solution_renumbered_vectorized[i]);
+                else
+                  solution_values[i + offset] =
+                    ETT::sum_value(comp,
+                                   is_linear ?
+                                     *(solution_values_vectorized_linear + i) :
+                                     solution_renumbered_vectorized[i]);
             }
           else
             {
               const unsigned int *renumber_ptr = renumber.data() + offset;
               for (unsigned int i = 0; i < dofs_per_comp; ++i)
-                solution_values[renumber_ptr[i]] =
-                  ETT::sum_value(comp, solution_renumbered_vectorized[i]);
+                if (sum_into_values)
+                  solution_values[renumber_ptr[i]] +=
+                    ETT::sum_value(comp, solution_renumbered_vectorized[i]);
+                else
+                  solution_values[renumber_ptr[i]] =
+                    ETT::sum_value(comp, solution_renumbered_vectorized[i]);
             }
         }
     }
@@ -2475,7 +2502,8 @@ template <bool        do_JxW,
 inline void
 FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate_fast(
   const StridedArrayView<ScalarNumber, stride_view> &solution_values,
-  const EvaluationFlags::EvaluationFlags            &integration_flags)
+  const EvaluationFlags::EvaluationFlags            &integration_flags,
+  const bool                                         sum_into_values)
 {
   // zero out lanes of incomplete last quadrature point batch
   if constexpr (stride == 1)
@@ -2547,7 +2575,8 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate_fast(
   finish_integrate_fast<is_face_path, is_linear>(
     solution_values,
     integration_flags,
-    solution_values_vectorized_linear.data());
+    solution_values_vectorized_linear.data(),
+    sum_into_values);
 }
 
 
@@ -2557,14 +2586,16 @@ template <bool do_JxW, std::size_t stride_view>
 inline void
 FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate_slow(
   const StridedArrayView<ScalarNumber, stride_view> &solution_values,
-  const EvaluationFlags::EvaluationFlags            &integration_flags)
+  const EvaluationFlags::EvaluationFlags            &integration_flags,
+  const bool                                         sum_into_values)
 {
   // slow path with FEValues
   Assert(fe_values.get() != nullptr,
          ExcMessage(
            "Not initialized. Please call FEPointEvaluation::reinit()!"));
-  for (unsigned int i = 0; i < solution_values.size(); ++i)
-    solution_values[i] = 0;
+  if (!sum_into_values)
+    for (unsigned int i = 0; i < solution_values.size(); ++i)
+      solution_values[i] = 0;
 
   const std::size_t n_points = fe_values->get_quadrature().size();
 
@@ -2634,7 +2665,8 @@ template <bool do_JxW, std::size_t stride_view>
 void
 FEPointEvaluation<n_components_, dim, spacedim, Number>::do_integrate(
   const StridedArrayView<ScalarNumber, stride_view> &solution_values,
-  const EvaluationFlags::EvaluationFlags            &integration_flags)
+  const EvaluationFlags::EvaluationFlags            &integration_flags,
+  const bool                                         sum_into_values)
 {
   if (!is_reinitialized)
     reinit();
@@ -2646,8 +2678,9 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::do_integrate(
         (integration_flags &
          EvaluationFlags::gradients))) // no integration flags
     {
-      for (unsigned int i = 0; i < solution_values.size(); ++i)
-        solution_values[i] = 0;
+      if (!sum_into_values)
+        for (unsigned int i = 0; i < solution_values.size(); ++i)
+          solution_values[i] = 0;
       return;
     }
 
@@ -2663,23 +2696,27 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::do_integrate(
         {
           if (use_linear_path)
             integrate_fast<do_JxW, true, true>(solution_values,
-                                               integration_flags);
+                                               integration_flags,
+                                               sum_into_values);
           else
             integrate_fast<do_JxW, true, false>(solution_values,
-                                                integration_flags);
+                                                integration_flags,
+                                                sum_into_values);
         }
       else
         {
           if (use_linear_path)
             integrate_fast<do_JxW, false, true>(solution_values,
-                                                integration_flags);
+                                                integration_flags,
+                                                sum_into_values);
           else
             integrate_fast<do_JxW, false, false>(solution_values,
-                                                 integration_flags);
+                                                 integration_flags,
+                                                 sum_into_values);
         }
     }
   else
-    integrate_slow<do_JxW>(solution_values, integration_flags);
+    integrate_slow<do_JxW>(solution_values, integration_flags, sum_into_values);
 }
 
 
@@ -2689,9 +2726,25 @@ template <std::size_t stride_view>
 void
 FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate(
   const StridedArrayView<ScalarNumber, stride_view> &solution_values,
-  const EvaluationFlags::EvaluationFlags            &integration_flags)
+  const EvaluationFlags::EvaluationFlags            &integration_flags,
+  const bool                                         sum_into_values)
+{
+  do_integrate<true>(solution_values, integration_flags, sum_into_values);
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate(
+  const ArrayView<ScalarNumber>          &solution_values,
+  const EvaluationFlags::EvaluationFlags &integration_flags,
+  const bool                              sum_into_values)
 {
-  do_integrate<true>(solution_values, integration_flags);
+  integrate(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
+                                              solution_values.size()),
+            integration_flags,
+            sum_into_values);
 }
 
 
@@ -2701,9 +2754,25 @@ template <std::size_t stride_view>
 void
 FEPointEvaluation<n_components_, dim, spacedim, Number>::test_and_sum(
   const StridedArrayView<ScalarNumber, stride_view> &solution_values,
-  const EvaluationFlags::EvaluationFlags            &integration_flags)
+  const EvaluationFlags::EvaluationFlags            &integration_flags,
+  const bool                                         sum_into_values)
+{
+  do_integrate<false>(solution_values, integration_flags, sum_into_values);
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::test_and_sum(
+  const ArrayView<ScalarNumber>          &solution_values,
+  const EvaluationFlags::EvaluationFlags &integration_flags,
+  const bool                              sum_into_values)
 {
-  do_integrate<false>(solution_values, integration_flags);
+  test_and_sum(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
+                                                 solution_values.size()),
+               integration_flags,
+               sum_into_values);
 }
 
 
diff --git a/tests/matrix_free/point_evaluation_28.cc b/tests/matrix_free/point_evaluation_28.cc
new file mode 100644 (file)
index 0000000..e91a8b3
--- /dev/null
@@ -0,0 +1,142 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 - 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.
+//
+// ---------------------------------------------------------------------
+
+
+// check FEPointEvaluation integrate_add()
+
+#include <deal.II/base/function_lib.h>
+
+#include <deal.II/dofs/dof_handler.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/tria.h>
+
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/matrix_free/fe_point_evaluation.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+
+
+template <int dim, typename Number = double>
+void
+test(const unsigned int degree)
+{
+  using namespace dealii;
+  Triangulation<dim> tria;
+
+  if (dim > 1)
+    GridGenerator::hyper_shell(tria, Point<dim>(), 0.5, 1, 6);
+  else
+    GridGenerator::subdivided_hyper_cube(tria, 2, 0, 1);
+
+  MappingQ<dim> mapping(degree);
+  deallog << "Mapping of degree " << degree << std::endl;
+
+  std::vector<Point<dim>> unit_points;
+  for (unsigned int i = 0; i < Utilities::pow(2, dim); ++i)
+    {
+      Point<dim> p;
+      for (unsigned int d = 0; d < dim; ++d)
+        p[d] = static_cast<Number>(i) / 17. + 0.015625 * d;
+      unit_points.push_back(p);
+    }
+
+  QGauss<dim> q_gauss(2);
+
+  Quadrature<dim> quad(unit_points, q_gauss.get_weights());
+
+  FE_Q<dim> fe(degree);
+
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+  Vector<Number> vector(dof_handler.n_dofs());
+
+  NonMatching::MappingInfo mapping_info(mapping,
+                                        update_values | update_gradients |
+                                          update_JxW_values);
+
+  std::vector<Quadrature<dim>> quad_vec_cell;
+  quad_vec_cell.reserve(tria.n_cells());
+
+  for (const auto &cell : tria.active_cell_iterators())
+    quad_vec_cell.push_back(quad);
+
+  mapping_info.reinit_cells(tria.active_cell_iterators(), quad_vec_cell);
+
+  FEPointEvaluation<1, dim, dim, Number> evaluator(mapping_info, fe);
+
+  Tensor<1, dim, Number> exponents;
+  exponents[0] = 1.;
+  VectorTools::interpolate(mapping,
+                           dof_handler,
+                           Functions::Monomial<dim, Number>(exponents),
+                           vector);
+
+  std::vector<Number> solution_values(fe.dofs_per_cell);
+
+  for (const auto &cell : dof_handler.active_cell_iterators())
+    {
+      cell->get_dof_values(vector,
+                           solution_values.begin(),
+                           solution_values.end());
+
+      evaluator.reinit(cell->active_cell_index());
+      evaluator.evaluate(solution_values,
+                         EvaluationFlags::values | EvaluationFlags::gradients);
+
+      for (const unsigned int i : evaluator.quadrature_point_indices())
+        {
+          evaluator.submit_value(evaluator.get_value(i), i);
+          evaluator.submit_gradient(evaluator.get_gradient(i), i);
+        }
+
+      evaluator.integrate(solution_values,
+                          EvaluationFlags::values | EvaluationFlags::gradients);
+
+      deallog << "Cell with center " << cell->center(true) << std::endl;
+      for (const auto i : solution_values)
+        deallog << i << ' ';
+      deallog << std::endl;
+
+      evaluator.integrate(solution_values,
+                          EvaluationFlags::values | EvaluationFlags::gradients,
+                          true);
+
+      for (const auto i : solution_values)
+        deallog << i << ' ';
+      deallog << std::endl;
+    }
+}
+
+
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(10);
+
+  test<2>(1);
+  test<2>(2);
+}
diff --git a/tests/matrix_free/point_evaluation_28.output b/tests/matrix_free/point_evaluation_28.output
new file mode 100644 (file)
index 0000000..ff41259
--- /dev/null
@@ -0,0 +1,39 @@
+
+DEAL::Mapping of degree 1
+DEAL::Cell with center 0.6495190528 0.3750000000
+DEAL::1.094128055 0.06405871628 -0.7129529194 -0.07186179469 
+DEAL::2.188256110 0.1281174326 -1.425905839 -0.1437235894 
+DEAL::Cell with center 4.592425497e-17 0.7500000000
+DEAL::0.5256716676 -0.3778726922 0.05767485942 -0.04324345223 
+DEAL::1.051343335 -0.7557453844 0.1153497188 -0.08648690447 
+DEAL::Cell with center -0.6495190528 0.3750000000
+DEAL::-0.5684563874 -0.4419314085 0.7706277789 0.02861834246 
+DEAL::-1.136912775 -0.8838628170 1.541255558 0.05723668491 
+DEAL::Cell with center -0.6495190528 -0.3750000000
+DEAL::-1.094128055 -0.06405871628 0.7129529194 0.07186179469 
+DEAL::-2.188256110 -0.1281174326 1.425905839 0.1437235894 
+DEAL::Cell with center -1.377727649e-16 -0.7500000000
+DEAL::-0.5256716676 0.3778726922 -0.05767485942 0.04324345223 
+DEAL::-1.051343335 0.7557453844 -0.1153497188 0.08648690447 
+DEAL::Cell with center 0.6495190528 -0.3750000000
+DEAL::0.5684563874 0.4419314085 -0.7706277789 -0.02861834246 
+DEAL::1.136912775 0.8838628170 -1.541255558 -0.05723668491 
+DEAL::Mapping of degree 2
+DEAL::Cell with center 0.6495190528 0.3750000000
+DEAL::2.482631076 -0.1564725804 0.4910617067 -0.02558423728 -2.474880528 0.1748480233 0.7176813542 0.1215956751 -0.8420572937 
+DEAL::4.965262152 -0.3129451607 0.9821234133 -0.05116847457 -4.949761056 0.3496960465 1.435362708 0.2431913503 -1.684114587 
+DEAL::Cell with center 4.592425497e-17 0.7500000000
+DEAL::1.903653500 0.1627020517 0.1365654459 -0.02631613009 -0.6782922627 0.1422251499 -0.8006550044 0.1439386250 -0.7778032374 
+DEAL::3.807306999 0.3254041034 0.2731308918 -0.05263226019 -1.356584525 0.2844502999 -1.601310009 0.2878772501 -1.555606475 
+DEAL::Cell with center -0.6495190528 0.3750000000
+DEAL::-0.5789775763 0.3191746321 -0.3544962608 -0.0007318928105 1.796588265 -0.03262287333 -1.518336359 0.02234294989 0.06425405628 
+DEAL::-1.157955153 0.6383492641 -0.7089925215 -0.001463785621 3.593176530 -0.06524574667 -3.036672717 0.04468589979 0.1285081126 
+DEAL::Cell with center -0.6495190528 -0.3750000000
+DEAL::-2.482631076 0.1564725804 -0.4910617067 0.02558423728 2.474880528 -0.1748480233 -0.7176813542 -0.1215956751 0.8420572937 
+DEAL::-4.965262152 0.3129451607 -0.9821234133 0.05116847457 4.949761056 -0.3496960465 -1.435362708 -0.2431913503 1.684114587 
+DEAL::Cell with center -1.377727649e-16 -0.7500000000
+DEAL::-1.903653500 -0.1627020517 -0.1365654459 0.02631613009 0.6782922627 -0.1422251499 0.8006550044 -0.1439386250 0.7778032374 
+DEAL::-3.807306999 -0.3254041034 -0.2731308918 0.05263226019 1.356584525 -0.2844502999 1.601310009 -0.2878772501 1.555606475 
+DEAL::Cell with center 0.6495190528 -0.3750000000
+DEAL::0.5789775763 -0.3191746321 0.3544962608 0.0007318928105 -1.796588265 0.03262287333 1.518336359 -0.02234294989 -0.06425405628 
+DEAL::1.157955153 -0.6383492641 0.7089925215 0.001463785621 -3.593176530 0.06524574667 3.036672717 -0.04468589979 -0.1285081126 

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.