From 1c5c8ee0c04d65d9e3a620b4f2059c1f2315d39e Mon Sep 17 00:00:00 2001
From: Peter Munch <peterrmuench@gmail.com>
Date: Fri, 15 Jul 2022 09:19:57 +0200
Subject: [PATCH] Generalized VT::point_gradients() for multiple components

Co-authored-by: Peter Munch <peterrmuench@gmail.com>
Co-authored-by: Martin Kronbichler <martin.kronbichler@uni-a.de>

rename file
---
 .../deal.II/numerics/vector_tools_evaluate.h  |  63 ++-
 .../vector_tools_evaluate_at_points_07.cc     | 176 ++++++++
 ...aluate_at_points_07.with_p4est=true.output | 404 ++++++++++++++++++
 3 files changed, 640 insertions(+), 3 deletions(-)
 create mode 100644 tests/remote_point_evaluation/vector_tools_evaluate_at_points_07.cc
 create mode 100644 tests/remote_point_evaluation/vector_tools_evaluate_at_points_07.with_p4est=true.output

diff --git a/include/deal.II/numerics/vector_tools_evaluate.h b/include/deal.II/numerics/vector_tools_evaluate.h
index 8660498930..fbb44a6ed7 100644
--- a/include/deal.II/numerics/vector_tools_evaluate.h
+++ b/include/deal.II/numerics/vector_tools_evaluate.h
@@ -343,6 +343,42 @@ namespace VectorTools
 
 
 
+    /**
+     * Perform reduction for tensors of tensors (e.g., gradient of
+     * vectorial quantities).
+     */
+    template <int n_components, int rank, int dim, typename Number>
+    Tensor<1, n_components, Tensor<rank, dim, Number>>
+    reduce(
+      const EvaluationFlags::EvaluationFlags &flags,
+      const ArrayView<const Tensor<1, n_components, Tensor<rank, dim, Number>>>
+        &values)
+    {
+      switch (flags)
+        {
+          case EvaluationFlags::avg:
+            {
+              Tensor<1, n_components, Tensor<rank, dim, Number>> temp;
+
+              for (unsigned int j = 0; j < values.size(); ++j)
+                for (unsigned int i = 0; i < n_components; ++i)
+                  temp[i] = temp[i] + values[j][i];
+
+              for (unsigned int i = 0; i < n_components; ++i)
+                temp[i] /= Number(values.size());
+
+              return temp;
+            }
+          case EvaluationFlags::insert:
+            return values[0];
+          default:
+            Assert(false, ExcNotImplemented());
+            return values[0];
+        }
+    }
+
+
+
     template <int n_components,
               int dim,
               int spacedim,
@@ -459,6 +495,26 @@ namespace VectorTools
 
 
 
+    template <typename Number, typename Number2>
+    void
+    set_value(Number &dst, const Number2 &src)
+    {
+      dst = src;
+    }
+
+
+
+    template <typename Number, int rank, int dim, typename Number2>
+    void
+    set_value(Tensor<rank, dim, Number> &, const Number2 &)
+    {
+      Assert(false,
+             ExcMessage(
+               "A cell-data vector can only have a single component."));
+    }
+
+
+
     template <int n_components,
               int dim,
               int spacedim,
@@ -491,8 +547,9 @@ namespace VectorTools
     {
       (void)evaluation_flags;
 
-      static_assert(n_components == 1,
-                    "A cell-data vector can only have a single component.");
+      Assert(n_components == 1,
+             ExcMessage(
+               "A cell-data vector can only have a single component."));
 
       Assert(evaluation_flags ==
                dealii::EvaluationFlags::EvaluationFlags::values,
@@ -506,7 +563,7 @@ namespace VectorTools
       for (unsigned int q = cell_data.reference_point_ptrs[i];
            q < cell_data.reference_point_ptrs[i + 1];
            ++q)
-        values[q] = value;
+        set_value(values[q], value);
     }
 
 
diff --git a/tests/remote_point_evaluation/vector_tools_evaluate_at_points_07.cc b/tests/remote_point_evaluation/vector_tools_evaluate_at_points_07.cc
new file mode 100644
index 0000000000..7a95765752
--- /dev/null
+++ b/tests/remote_point_evaluation/vector_tools_evaluate_at_points_07.cc
@@ -0,0 +1,176 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+// Evaluate solution vector along a line.
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_control.h>
+
+#include <deal.II/matrix_free/fe_point_evaluation.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/vector_tools_evaluate.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <typename MeshType>
+MPI_Comm
+get_mpi_comm(const MeshType &mesh)
+{
+  const auto *tria_parallel = dynamic_cast<
+    const parallel::TriangulationBase<MeshType::dimension,
+                                      MeshType::space_dimension> *>(
+    &(mesh.get_triangulation()));
+
+  return tria_parallel != nullptr ? tria_parallel->get_communicator() :
+                                    MPI_COMM_SELF;
+}
+
+template <int dim, int spacedim>
+std::shared_ptr<const Utilities::MPI::Partitioner>
+create_partitioner(const DoFHandler<dim, spacedim> &dof_handler)
+{
+  IndexSet locally_relevant_dofs;
+
+  DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
+
+  return std::make_shared<const Utilities::MPI::Partitioner>(
+    dof_handler.locally_owned_dofs(),
+    locally_relevant_dofs,
+    get_mpi_comm(dof_handler));
+}
+
+template <int dim>
+void
+print(const Mapping<dim> &                              mapping,
+      const DoFHandler<dim> &                           dof_handler,
+      const LinearAlgebra::distributed::Vector<double> &result,
+      const unsigned int                                counter)
+{
+  result.update_ghost_values();
+
+  DataOutBase::VtkFlags flags;
+  flags.write_higher_order_cells = true;
+
+  DataOut<dim> data_out;
+  data_out.set_flags(flags);
+  data_out.attach_dof_handler(dof_handler);
+
+  const auto &tria = dof_handler.get_triangulation();
+
+  Vector<double> ranks(tria.n_active_cells());
+  for (const auto &cell :
+       tria.active_cell_iterators() | IteratorFilters::LocallyOwnedCell())
+    ranks(cell->active_cell_index()) = cell->subdomain_id();
+  data_out.add_data_vector(ranks, "rank");
+  data_out.add_data_vector(result, "result");
+
+  data_out.build_patches(mapping, dof_handler.get_fe().tensor_degree() + 1);
+  data_out.write_vtu_with_pvtu_record(
+    "./", "example-7", counter, MPI_COMM_WORLD, 1, 1);
+
+  result.zero_out_ghost_values();
+}
+
+template <int dim>
+class AnalyticalFunction : public Function<dim>
+{
+public:
+  AnalyticalFunction(const unsigned int n_components)
+    : Function<dim>(n_components)
+  {}
+
+  virtual double
+  value(const Point<dim> &p, const unsigned int component) const
+  {
+    return p[component] * p[component];
+  }
+};
+
+void
+test()
+{
+  const unsigned int dim             = 2;
+  const unsigned int degree          = 2;
+  const unsigned int n_components    = dim;
+  const unsigned int n_refinements_1 = 3;
+
+  parallel::distributed::Triangulation<dim> tria_1(MPI_COMM_WORLD);
+  GridGenerator::hyper_cube(tria_1, -1.0, +1.0);
+  tria_1.refine_global(n_refinements_1);
+  DoFHandler<dim>       dof_handler_1(tria_1);
+  dealii::FESystem<dim> fe_1(dealii::FE_Q<dim>(degree), n_components);
+  dof_handler_1.distribute_dofs(fe_1);
+
+  LinearAlgebra::distributed::Vector<double> vector_1(
+    create_partitioner(dof_handler_1));
+
+  const MappingQ1<dim> mapping_1;
+  VectorTools::interpolate(mapping_1,
+                           dof_handler_1,
+                           AnalyticalFunction<dim>(n_components),
+                           vector_1);
+
+  std::vector<Point<dim>> evaluation_points;
+
+  const unsigned int n_intervals = 100;
+
+  if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+    for (unsigned int i = 0; i <= n_intervals; ++i)
+      evaluation_points.emplace_back(-1.0 + 2.0 / n_intervals * i, 0.0);
+
+  vector_1.update_ghost_values();
+  Utilities::MPI::RemotePointEvaluation<dim> evaluation_cache;
+  const auto                                 evaluation_point_gradient_results =
+    VectorTools::point_gradients<n_components>(
+      mapping_1, dof_handler_1, vector_1, evaluation_points, evaluation_cache);
+  vector_1.zero_out_ghost_values();
+
+  if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
+    for (unsigned int i = 0; i <= n_intervals; ++i)
+      for (unsigned int c1 = 0; c1 < n_components; ++c1)
+        for (unsigned int c2 = 0; c2 < n_components; ++c2)
+          deallog << evaluation_point_gradient_results[i][c1][c2] << std::endl;
+}
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
+  initlog();
+  deallog << std::setprecision(8) << std::fixed;
+
+  test();
+}
diff --git a/tests/remote_point_evaluation/vector_tools_evaluate_at_points_07.with_p4est=true.output b/tests/remote_point_evaluation/vector_tools_evaluate_at_points_07.with_p4est=true.output
new file mode 100644
index 0000000000..b7bd8d724a
--- /dev/null
+++ b/tests/remote_point_evaluation/vector_tools_evaluate_at_points_07.with_p4est=true.output
@@ -0,0 +1,404 @@
+-2.00000000
+0.00000000
+0.00000000
+0.00000000
+-1.96000000
+0.00000000
+0.00000000
+0.00000000
+-1.92000000
+-0.00000000
+0.00000000
+0.00000000
+-1.88000000
+0.00000000
+0.00000000
+0.00000000
+-1.84000000
+0.00000000
+0.00000000
+0.00000000
+-1.80000000
+0.00000000
+0.00000000
+0.00000000
+-1.76000000
+0.00000000
+0.00000000
+0.00000000
+-1.72000000
+-0.00000000
+0.00000000
+0.00000000
+-1.68000000
+0.00000000
+0.00000000
+0.00000000
+-1.64000000
+0.00000000
+0.00000000
+0.00000000
+-1.60000000
+0.00000000
+0.00000000
+0.00000000
+-1.56000000
+0.00000000
+0.00000000
+0.00000000
+-1.52000000
+-0.00000000
+0.00000000
+0.00000000
+-1.48000000
+0.00000000
+0.00000000
+0.00000000
+-1.44000000
+0.00000000
+0.00000000
+0.00000000
+-1.40000000
+0.00000000
+0.00000000
+0.00000000
+-1.36000000
+0.00000000
+0.00000000
+0.00000000
+-1.32000000
+0.00000000
+0.00000000
+0.00000000
+-1.28000000
+-0.00000000
+0.00000000
+0.00000000
+-1.24000000
+-0.00000000
+0.00000000
+0.00000000
+-1.20000000
+-0.00000000
+0.00000000
+0.00000000
+-1.16000000
+0.00000000
+0.00000000
+0.00000000
+-1.12000000
+0.00000000
+0.00000000
+0.00000000
+-1.08000000
+0.00000000
+0.00000000
+0.00000000
+-1.04000000
+-0.00000000
+0.00000000
+0.00000000
+-1.00000000
+0.00000000
+0.00000000
+0.00000000
+-0.96000000
+-0.00000000
+0.00000000
+0.00000000
+-0.92000000
+0.00000000
+0.00000000
+0.00000000
+-0.88000000
+0.00000000
+0.00000000
+0.00000000
+-0.84000000
+-0.00000000
+0.00000000
+0.00000000
+-0.80000000
+0.00000000
+0.00000000
+0.00000000
+-0.76000000
+-0.00000000
+0.00000000
+0.00000000
+-0.72000000
+0.00000000
+0.00000000
+0.00000000
+-0.68000000
+-0.00000000
+0.00000000
+0.00000000
+-0.64000000
+0.00000000
+0.00000000
+0.00000000
+-0.60000000
+0.00000000
+0.00000000
+0.00000000
+-0.56000000
+0.00000000
+0.00000000
+0.00000000
+-0.52000000
+-0.00000000
+0.00000000
+0.00000000
+-0.48000000
+-0.00000000
+0.00000000
+0.00000000
+-0.44000000
+0.00000000
+0.00000000
+0.00000000
+-0.40000000
+0.00000000
+0.00000000
+0.00000000
+-0.36000000
+0.00000000
+0.00000000
+0.00000000
+-0.32000000
+0.00000000
+0.00000000
+0.00000000
+-0.28000000
+0.00000000
+0.00000000
+0.00000000
+-0.24000000
+-0.00000000
+0.00000000
+0.00000000
+-0.20000000
+0.00000000
+0.00000000
+0.00000000
+-0.16000000
+-0.00000000
+0.00000000
+0.00000000
+-0.12000000
+0.00000000
+0.00000000
+0.00000000
+-0.08000000
+0.00000000
+0.00000000
+0.00000000
+-0.04000000
+0.00000000
+0.00000000
+0.00000000
+0.00000000
+0.00000000
+0.00000000
+0.00000000
+0.04000000
+0.00000000
+0.00000000
+0.00000000
+0.08000000
+0.00000000
+0.00000000
+0.00000000
+0.12000000
+-0.00000000
+0.00000000
+0.00000000
+0.16000000
+-0.00000000
+0.00000000
+0.00000000
+0.20000000
+-0.00000000
+0.00000000
+0.00000000
+0.24000000
+-0.00000000
+0.00000000
+0.00000000
+0.28000000
+0.00000000
+0.00000000
+0.00000000
+0.32000000
+-0.00000000
+0.00000000
+0.00000000
+0.36000000
+0.00000000
+0.00000000
+0.00000000
+0.40000000
+0.00000000
+0.00000000
+0.00000000
+0.44000000
+0.00000000
+0.00000000
+0.00000000
+0.48000000
+-0.00000000
+0.00000000
+0.00000000
+0.52000000
+0.00000000
+0.00000000
+0.00000000
+0.56000000
+0.00000000
+0.00000000
+0.00000000
+0.60000000
+0.00000000
+0.00000000
+0.00000000
+0.64000000
+0.00000000
+0.00000000
+0.00000000
+0.68000000
+0.00000000
+0.00000000
+0.00000000
+0.72000000
+-0.00000000
+0.00000000
+0.00000000
+0.76000000
+0.00000000
+0.00000000
+0.00000000
+0.80000000
+0.00000000
+0.00000000
+0.00000000
+0.84000000
+-0.00000000
+0.00000000
+0.00000000
+0.88000000
+0.00000000
+0.00000000
+0.00000000
+0.92000000
+0.00000000
+0.00000000
+0.00000000
+0.96000000
+-0.00000000
+0.00000000
+0.00000000
+1.00000000
+0.00000000
+0.00000000
+0.00000000
+1.04000000
+0.00000000
+0.00000000
+0.00000000
+1.08000000
+0.00000000
+0.00000000
+0.00000000
+1.12000000
+0.00000000
+0.00000000
+0.00000000
+1.16000000
+0.00000000
+0.00000000
+0.00000000
+1.20000000
+0.00000000
+0.00000000
+0.00000000
+1.24000000
+-0.00000000
+0.00000000
+0.00000000
+1.28000000
+0.00000000
+0.00000000
+0.00000000
+1.32000000
+0.00000000
+0.00000000
+0.00000000
+1.36000000
+0.00000000
+0.00000000
+0.00000000
+1.40000000
+0.00000000
+0.00000000
+0.00000000
+1.44000000
+0.00000000
+0.00000000
+0.00000000
+1.48000000
+0.00000000
+0.00000000
+0.00000000
+1.52000000
+-0.00000000
+0.00000000
+0.00000000
+1.56000000
+0.00000000
+0.00000000
+0.00000000
+1.60000000
+0.00000000
+0.00000000
+0.00000000
+1.64000000
+0.00000000
+0.00000000
+0.00000000
+1.68000000
+-0.00000000
+0.00000000
+0.00000000
+1.72000000
+0.00000000
+0.00000000
+0.00000000
+1.76000000
+-0.00000000
+0.00000000
+0.00000000
+1.80000000
+0.00000000
+0.00000000
+0.00000000
+1.84000000
+0.00000000
+0.00000000
+0.00000000
+1.88000000
+0.00000000
+0.00000000
+0.00000000
+1.92000000
+-0.00000000
+0.00000000
+0.00000000
+1.96000000
+0.00000000
+0.00000000
+0.00000000
+2.00000000
+0.00000000
+0.00000000
+0.00000000
-- 
2.39.5