Fix FEPointEvaluation for non-zero starting component

doc/news/changes/minor/20211001Gassmoeller

@@ -0,0 +1,3 @@
+Fixed: The FEPointEvaluation class always evaluated a solution starting at the component 0 of the selected vector-valued base element, even if a different component was selected by the caller. This has been fixed to correctly start evaluating at the selected component.
+<br>
+(Rene Gassmoeller, 2021/10/01)

include/deal.II/matrix_free/fe_point_evaluation.h

@@ -653,6 +653,11 @@ class FEPointEvaluation
    */
   unsigned int dofs_per_component;
 
+  /**
+   * The first selected component in the active base element.
+   */
+  unsigned int component_in_base_element;
+
   /**
    * For complicated FiniteElement objects this variable informs us about
    * which unknowns actually carry degrees of freedom in the selected
@@ -742,18 +747,28 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::FEPointEvaluation(
 
   bool         same_base_element   = true;
   unsigned int base_element_number = 0;
+  component_in_base_element        = 0;
   unsigned int component           = 0;
   for (; base_element_number < fe.n_base_elements(); ++base_element_number)
-    if (component + fe.element_multiplicity(base_element_number) >
-        first_selected_component)
-      {
-        if (first_selected_component + n_components >
-            component + fe.element_multiplicity(base_element_number))
-          same_base_element = false;
+    {
+      bool found_base_element = false;
+      for (component_in_base_element = 0;
+           component_in_base_element <
+           fe.element_multiplicity(base_element_number);
+           ++component_in_base_element, ++component)
+        if (component == first_selected_component)
+          {
+            found_base_element = true;
+
+            if (component_in_base_element + n_components >
+                fe.element_multiplicity(base_element_number))
+              same_base_element = false;
+            break;
+          }
+      if (found_base_element == true)
         break;
-      }
-    else
-      component += fe.element_multiplicity(base_element_number);
+    }
+
   if (fill_mapping_data_for_generic_points &&
       internal::FEPointEvaluation::is_fast_path_supported(
         fe, base_element_number) &&
@@ -868,7 +883,9 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::evaluate(
         for (unsigned int i = 0; i < dofs_per_component; ++i)
           internal::FEPointEvaluation::
             EvaluatorTypeTraits<dim, n_components, Number>::read_value(
-              solution_values[renumber[comp * dofs_per_component + i]],
+              solution_values[renumber[(component_in_base_element + comp) *
+                                         dofs_per_component +
+                                       i]],
               comp,
               solution_renumbered[i]);
 

tests/matrix_free/point_evaluation_14.cc

@@ -0,0 +1,160 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 - 2021 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 for vector-valued FE_Q and MappingQ
+// with a starting component in the middle of a vector-valued
+// base element.
+
+#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/fe_system.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/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>
+class MyFunction : public Function<dim>
+{
+public:
+  MyFunction()
+    : Function<dim>(dim)
+  {}
+
+  double
+  value(const Point<dim> &p, const unsigned int component) const override
+  {
+    return component + p[component];
+  }
+};
+
+
+
+template <int dim>
+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 < 13; ++i)
+    {
+      Point<dim> p;
+      for (unsigned int d = 0; d < dim; ++d)
+        p[d] = static_cast<double>(i) / 17. + 0.015625 * d;
+      unit_points.push_back(p);
+    }
+
+  FESystem<dim> fe(FE_Q<dim>(degree), dim);
+  FEValues<dim> fe_values(mapping,
+                          fe,
+                          Quadrature<dim>(unit_points),
+                          update_values | update_gradients);
+
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+  Vector<double> vector(dof_handler.n_dofs());
+
+  FEPointEvaluation<1, dim> evaluator(mapping,
+                                      fe,
+                                      update_values | update_gradients,
+                                      1);
+
+  VectorTools::interpolate(mapping, dof_handler, MyFunction<dim>(), vector);
+
+  std::vector<double>         solution_values(fe.dofs_per_cell);
+  std::vector<double>         function_values(unit_points.size());
+  std::vector<Tensor<1, dim>> function_gradients(unit_points.size());
+
+  FEValuesExtractors::Scalar extractor(1);
+
+  for (const auto &cell : dof_handler.active_cell_iterators())
+    {
+      fe_values.reinit(cell);
+      fe_values[extractor].get_function_values(vector, function_values);
+      fe_values[extractor].get_function_gradients(vector, function_gradients);
+
+      cell->get_dof_values(vector,
+                           solution_values.begin(),
+                           solution_values.end());
+
+      evaluator.reinit(cell, unit_points);
+      evaluator.evaluate(solution_values,
+                         EvaluationFlags::values | EvaluationFlags::gradients);
+
+      deallog << "Cell with center " << cell->center(true) << std::endl;
+      for (unsigned int i = 0; i < function_values.size(); ++i)
+        deallog << mapping.transform_unit_to_real_cell(cell, unit_points[i])
+                << ": " << evaluator.get_value(i) << " error value "
+                << (function_values[i] - evaluator.get_value(i))
+                << " error grad "
+                << (evaluator.get_gradient(i) - function_gradients[i]).norm()
+                << std::endl;
+
+      deallog << std::endl;
+
+      for (unsigned int i = 0; i < unit_points.size(); ++i)
+        {
+          evaluator.submit_value(evaluator.get_value(i), i);
+          evaluator.submit_gradient(evaluator.get_gradient(i), i);
+        }
+
+      evaluator.integrate(solution_values,
+                          EvaluationFlags::values | EvaluationFlags::gradients);
+
+      for (const auto i : solution_values)
+        deallog << i << " ";
+      deallog << std::endl;
+    }
+}
+
+
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(10);
+
+  test<2>(2);
+  test<2>(6);
+  test<3>(5);
+}