Only zero out values FEPointEvaluation is working on

doc/news/changes/minor/20231224Heinz

@@ -0,0 +1,6 @@
+Fixed: From now on FEPointEvaluation::integrate() with sum_into_values = false
+only zeroes out the values it is working on (instead of all values) and is now
+consistent to FEEvaluation.
+
+<br>
+(Johannes Heinz, 2023/12/24)

include/deal.II/matrix_free/fe_point_evaluation.h

@@ -2405,10 +2405,6 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::finish_integrate_fast(
   vectorized_value_type *solution_values_vectorized_linear,
   const bool             sum_into_values)
 {
-  if (!sum_into_values && fe->n_components() > n_components)
-    for (unsigned int i = 0; i < solution_values.size(); ++i)
-      solution_values[i] = 0;
-
   const unsigned int dofs_per_comp =
     is_linear ? Utilities::pow(2, dim) : dofs_per_component;
 
@@ -2593,8 +2589,14 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate_slow(
   Assert(fe_values.get() != nullptr,
          ExcMessage(
            "Not initialized. Please call FEPointEvaluation::reinit()!"));
+
+  const unsigned int dofs_per_comp =
+    use_linear_path ? Utilities::pow(2, dim) : dofs_per_component;
+
   if (!sum_into_values)
-    for (unsigned int i = 0; i < solution_values.size(); ++i)
+    for (unsigned int i = component_in_base_element * dofs_per_comp;
+         i < (component_in_base_element + n_components) * dofs_per_comp;
+         ++i)
       solution_values[i] = 0;
 
   const std::size_t n_points = fe_values->get_quadrature().size();
@@ -2678,9 +2680,15 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::do_integrate(
         (integration_flags &
          EvaluationFlags::gradients))) // no integration flags
     {
+      const unsigned int dofs_per_comp =
+        use_linear_path ? Utilities::pow(2, dim) : dofs_per_component;
+
       if (!sum_into_values)
-        for (unsigned int i = 0; i < solution_values.size(); ++i)
+        for (unsigned int i = component_in_base_element * dofs_per_comp;
+             i < (component_in_base_element + n_components) * dofs_per_comp;
+             ++i)
           solution_values[i] = 0;
+
       return;
     }
 

tests/matrix_free/fe_point_evaluation_01.cc

@@ -0,0 +1,101 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2023 - 2022 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 that during FEPointEvaluation::integrate(), only the components in
+// provided buffer are zeroed out that integrate writes to.
+
+#include <deal.II/base/mpi.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/matrix_free/fe_point_evaluation.h>
+
+#include <deal.II/non_matching/mapping_info.h>
+
+#include <deal.II/../../tests/tests.h>
+
+using namespace dealii;
+
+template <int dim, typename Number>
+void
+do_test(unsigned int degree)
+{
+  // construct basic triangulation
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+  GridGenerator::hyper_cube(tria);
+
+  // construct and fill a vector of quadratures for each face of the cells
+  std::vector<std::vector<Quadrature<dim - 1>>> quadrature_vector;
+  for (const auto &cell : tria.active_cell_iterators())
+    quadrature_vector.emplace_back(
+      std::vector<Quadrature<dim - 1>>(cell->n_faces()));
+
+  for (const auto &cell : tria.active_cell_iterators())
+    if (cell->is_locally_owned())
+      for (unsigned int f = 0; f < cell->n_faces(); ++f)
+        if (cell->face(f)->at_boundary())
+          quadrature_vector[cell->active_cell_index()][f] =
+            Quadrature<dim - 1>(degree + 1);
+
+
+  // setup NonMatching::MappingInfo
+  MappingQ1<dim> mapping;
+
+  NonMatching::MappingInfo<dim, dim, Number> nm_mapping_info(
+    mapping, update_values | update_JxW_values);
+
+  nm_mapping_info.reinit_faces(tria.active_cell_iterators(), quadrature_vector);
+
+
+  // initialize FEPointEvaluation objects
+  FESystem<dim>                          fe(FE_DGQ<dim>(degree), dim + 1);
+  FEPointEvaluation<1, dim, dim, Number> fe_eval(nm_mapping_info, fe, 1);
+  AlignedVector<Number>                  buffer(fe.dofs_per_cell);
+  for (auto &b : buffer)
+    b = 1.0;
+
+  const auto  &cell = tria.begin_active();
+  unsigned int f    = 0;
+  {
+    fe_eval.reinit(cell->active_cell_index(), f);
+
+    fe_eval.integrate(buffer,
+                      EvaluationFlags::values,
+                      /*sum_into_values=*/false);
+  }
+
+  for (const auto &b : buffer)
+    deallog << std::abs(b) << ", ";
+  deallog << std::endl;
+}
+
+int
+main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
+  MPILogInitAll                    all;
+
+  do_test<2, double>(2);
+
+  return 0;
+}

tests/matrix_free/fe_point_evaluation_01.with_p4est=true.mpirun=1.output

@@ -0,0 +1,2 @@
+
+DEAL:0::1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000, 1.00000,