]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Introduce VectorTools::point_gradients() 12333/head
authorPeter Munch <peterrmuench@gmail.com>
Sun, 30 May 2021 06:52:56 +0000 (08:52 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Sun, 30 May 2021 06:52:56 +0000 (08:52 +0200)
include/deal.II/numerics/vector_tools_evaluate.h
tests/remote_point_evaluation/vector_tools_evaluate_at_points_04.cc

index cde9f9508f719045ccd82ab07724ec5d0f58539a..ec8aec1776e3738439fbdc6db37fb28347fedad0 100644 (file)
@@ -32,12 +32,12 @@ DEAL_II_NAMESPACE_OPEN
 namespace VectorTools
 {
   /**
-   * Namespace for the flags for point_values().
+   * Namespace for the flags for point_values() and point_gradients().
    */
   namespace EvaluationFlags
   {
     /**
-     * Flags for point_values().
+     * Flags for point_values() and point_gradients().
      */
     enum EvaluationFlags
     {
@@ -100,6 +100,36 @@ namespace VectorTools
     const VectorType &                                          vector,
     const EvaluationFlags::EvaluationFlags flags = EvaluationFlags::avg);
 
+  /**
+   * Given a (distributed) solution vector @p vector, evaluate the gradients at
+   * the (arbitrary and even remote) points specified by @p evaluation_points.
+   */
+  template <int n_components, int dim, int spacedim, typename VectorType>
+  std::vector<typename FEPointEvaluation<n_components, dim>::gradient_type>
+  point_gradients(
+    const Mapping<dim> &                                  mapping,
+    const DoFHandler<dim, spacedim> &                     dof_handler,
+    const VectorType &                                    vector,
+    const std::vector<Point<spacedim>> &                  evaluation_points,
+    Utilities::MPI::RemotePointEvaluation<dim, spacedim> &cache,
+    const EvaluationFlags::EvaluationFlags flags = EvaluationFlags::avg);
+
+  /**
+   * Given a (distributed) solution vector @p vector, evaluate the gradients at
+   * the points specified by @p cache which might have been set up by the
+   * above function.
+   *
+   * @note Refinement/coarsening/repartitioning leads to the invalidation of the
+   *   cache so that the above function has to be called again.
+   */
+  template <int n_components, int dim, int spacedim, typename VectorType>
+  std::vector<typename FEPointEvaluation<n_components, dim>::gradient_type>
+  point_gradients(
+    const Utilities::MPI::RemotePointEvaluation<dim, spacedim> &cache,
+    const DoFHandler<dim, spacedim> &                           dof_handler,
+    const VectorType &                                          vector,
+    const EvaluationFlags::EvaluationFlags flags = EvaluationFlags::avg);
+
 
 
   // inlined functions
@@ -122,6 +152,23 @@ namespace VectorTools
 
 
 
+  template <int n_components, int dim, int spacedim, typename VectorType>
+  inline std::vector<
+    typename FEPointEvaluation<n_components, dim>::gradient_type>
+  point_gradients(const Mapping<dim> &                mapping,
+                  const DoFHandler<dim, spacedim> &   dof_handler,
+                  const VectorType &                  vector,
+                  const std::vector<Point<spacedim>> &evaluation_points,
+                  Utilities::MPI::RemotePointEvaluation<dim, spacedim> &cache,
+                  const EvaluationFlags::EvaluationFlags                flags)
+  {
+    cache.reinit(evaluation_points, dof_handler.get_triangulation(), mapping);
+
+    return point_gradients<n_components>(cache, dof_handler, vector, flags);
+  }
+
+
+
   namespace internal
   {
     /**
@@ -177,10 +224,133 @@ namespace VectorTools
             return values[0];
         }
     }
-  } // namespace internal
 
 
 
+    template <int n_components,
+              int dim,
+              int spacedim,
+              typename VectorType,
+              typename value_type>
+    inline std::vector<value_type>
+    evaluate_at_points(
+      const Utilities::MPI::RemotePointEvaluation<dim, spacedim> &cache,
+      const DoFHandler<dim, spacedim> &                           dof_handler,
+      const VectorType &                                          vector,
+      const EvaluationFlags::EvaluationFlags                      flags,
+      const UpdateFlags                                           update_flags,
+      const dealii::EvaluationFlags::EvaluationFlags evaluation_flags,
+      const std::function<
+        value_type(const FEPointEvaluation<n_components, dim> &,
+                   const unsigned int &)> process_quadrature_point)
+    {
+      Assert(cache.is_ready(),
+             ExcMessage(
+               "Utilities::MPI::RemotePointEvaluation is not ready yet! "
+               "Please call Utilities::MPI::RemotePointEvaluation::reinit() "
+               "yourself or another function that does this for you."));
+
+      Assert(
+        &dof_handler.get_triangulation() == &cache.get_triangulation(),
+        ExcMessage(
+          "The provided Utilities::MPI::RemotePointEvaluation and DoFHandler "
+          "object have been set up with different Triangulation objects, "
+          "a scenario not supported!"));
+
+      // evaluate values at points if possible
+      const auto evaluation_point_results = [&]() {
+        // helper function for accessing the global vector and interpolating
+        // the results onto the points
+        const auto evaluation_function = [&](auto &      values,
+                                             const auto &cell_data) {
+          std::vector<typename VectorType::value_type> solution_values;
+
+          std::vector<std::unique_ptr<FEPointEvaluation<n_components, dim>>>
+            evaluators(dof_handler.get_fe_collection().size());
+
+          const auto get_evaluator = [&](const unsigned int active_fe_index)
+            -> FEPointEvaluation<n_components, dim> & {
+            if (evaluators[active_fe_index] == nullptr)
+              evaluators[active_fe_index] =
+                std::make_unique<FEPointEvaluation<n_components, dim>>(
+                  cache.get_mapping(),
+                  dof_handler.get_fe(active_fe_index),
+                  update_flags);
+
+            return *evaluators[active_fe_index];
+          };
+
+          for (unsigned int i = 0; i < cell_data.cells.size(); ++i)
+            {
+              typename DoFHandler<dim>::active_cell_iterator cell = {
+                &cache.get_triangulation(),
+                cell_data.cells[i].first,
+                cell_data.cells[i].second,
+                &dof_handler};
+
+              const ArrayView<const Point<dim>> unit_points(
+                cell_data.reference_point_values.data() +
+                  cell_data.reference_point_ptrs[i],
+                cell_data.reference_point_ptrs[i + 1] -
+                  cell_data.reference_point_ptrs[i]);
+
+              solution_values.resize(
+                dof_handler.get_fe(cell->active_fe_index()).n_dofs_per_cell());
+              cell->get_dof_values(vector,
+                                   solution_values.begin(),
+                                   solution_values.end());
+
+              auto &evaluator = get_evaluator(cell->active_fe_index());
+
+              evaluator.reinit(cell, unit_points);
+              evaluator.evaluate(solution_values, evaluation_flags);
+
+              for (unsigned int q = 0; q < unit_points.size(); ++q)
+                values[q + cell_data.reference_point_ptrs[i]] =
+                  process_quadrature_point(evaluator, q);
+            }
+        };
+
+        std::vector<value_type> evaluation_point_results;
+        std::vector<value_type> buffer;
+
+        cache.template evaluate_and_process<value_type>(
+          evaluation_point_results, buffer, evaluation_function);
+
+        return evaluation_point_results;
+      }();
+
+      if (cache.is_map_unique())
+        {
+          // each point has exactly one result (unique map)
+          return evaluation_point_results;
+        }
+      else
+        {
+          // map is not unique (multiple or no results): postprocessing is
+          // needed
+          std::vector<value_type> unique_evaluation_point_results(
+            cache.get_point_ptrs().size() - 1);
+
+          const auto &ptr = cache.get_point_ptrs();
+
+          for (unsigned int i = 0; i < ptr.size() - 1; ++i)
+            {
+              const auto n_entries = ptr[i + 1] - ptr[i];
+              if (n_entries == 0)
+                continue;
+
+              unique_evaluation_point_results[i] =
+                reduce(flags,
+                       ArrayView<const value_type>(
+                         evaluation_point_results.data() + ptr[i], n_entries));
+            }
+
+          return unique_evaluation_point_results;
+        }
+    }
+  } // namespace internal
+
   template <int n_components, int dim, int spacedim, typename VectorType>
   inline std::vector<typename FEPointEvaluation<n_components, dim>::value_type>
   point_values(
@@ -189,117 +359,49 @@ namespace VectorTools
     const VectorType &                                          vector,
     const EvaluationFlags::EvaluationFlags                      flags)
   {
-    using value_type =
-      typename FEPointEvaluation<n_components, dim>::value_type;
-
-    Assert(cache.is_ready(),
-           ExcMessage(
-             "Utilities::MPI::RemotePointEvaluation is not ready yet! "
-             "Please call Utilities::MPI::RemotePointEvaluation::reinit() "
-             "yourself or the other point_values(), which does this for"
-             "you."));
-
-    Assert(
-      &dof_handler.get_triangulation() == &cache.get_triangulation(),
-      ExcMessage(
-        "The provided Utilities::MPI::RemotePointEvaluation and DoFHandler "
-        "object have been set up with different Triangulation objects, "
-        "a scenario not supported!"));
-
-    // evaluate values at points if possible
-    const auto evaluation_point_results = [&]() {
-      // helper function for accessing the global vector and interpolating
-      // the results onto the points
-      const auto evaluation_function = [&](auto &      values,
-                                           const auto &cell_data) {
-        std::vector<typename VectorType::value_type> solution_values;
-
-        std::vector<std::unique_ptr<FEPointEvaluation<n_components, dim>>>
-          evaluators(dof_handler.get_fe_collection().size());
-
-        const auto get_evaluator = [&](const unsigned int active_fe_index)
-          -> FEPointEvaluation<n_components, dim> & {
-          if (evaluators[active_fe_index] == nullptr)
-            evaluators[active_fe_index] =
-              std::make_unique<FEPointEvaluation<n_components, dim>>(
-                cache.get_mapping(),
-                dof_handler.get_fe(active_fe_index),
-                update_values);
-
-          return *evaluators[active_fe_index];
-        };
+    return internal::evaluate_at_points<
+      n_components,
+      dim,
+      spacedim,
+      VectorType,
+      typename FEPointEvaluation<n_components, dim>::value_type>(
+      cache,
+      dof_handler,
+      vector,
+      flags,
+      update_values,
+      dealii::EvaluationFlags::values,
+      [](const auto &evaluator, const auto &q) {
+        return evaluator.get_value(q);
+      });
+  }
 
-        for (unsigned int i = 0; i < cell_data.cells.size(); ++i)
-          {
-            typename DoFHandler<dim>::active_cell_iterator cell = {
-              &cache.get_triangulation(),
-              cell_data.cells[i].first,
-              cell_data.cells[i].second,
-              &dof_handler};
-
-            const ArrayView<const Point<dim>> unit_points(
-              cell_data.reference_point_values.data() +
-                cell_data.reference_point_ptrs[i],
-              cell_data.reference_point_ptrs[i + 1] -
-                cell_data.reference_point_ptrs[i]);
-
-            solution_values.resize(
-              dof_handler.get_fe(cell->active_fe_index()).n_dofs_per_cell());
-            cell->get_dof_values(vector,
-                                 solution_values.begin(),
-                                 solution_values.end());
-
-            auto &evaluator = get_evaluator(cell->active_fe_index());
-
-            evaluator.reinit(cell, unit_points);
-            evaluator.evaluate(solution_values,
-                               dealii::EvaluationFlags::values);
-
-            for (unsigned int q = 0; q < unit_points.size(); ++q)
-              values[q + cell_data.reference_point_ptrs[i]] =
-                evaluator.get_value(q);
-          }
-      };
-
-      std::vector<value_type> evaluation_point_results;
-      std::vector<value_type> buffer;
-
-      cache.template evaluate_and_process<value_type>(evaluation_point_results,
-                                                      buffer,
-                                                      evaluation_function);
-
-      return evaluation_point_results;
-    }();
-
-    if (cache.is_map_unique())
-      {
-        // each point has exactly one result (unique map)
-        return evaluation_point_results;
-      }
-    else
-      {
-        // map is not unique (multiple or no results): postprocessing is needed
-        std::vector<value_type> unique_evaluation_point_results(
-          cache.get_point_ptrs().size() - 1);
-
-        const auto &ptr = cache.get_point_ptrs();
-
-        for (unsigned int i = 0; i < ptr.size() - 1; ++i)
-          {
-            const auto n_entries = ptr[i + 1] - ptr[i];
-            if (n_entries == 0)
-              continue;
-
-            unique_evaluation_point_results[i] =
-              internal::reduce(flags,
-                               ArrayView<const value_type>(
-                                 evaluation_point_results.data() + ptr[i],
-                                 n_entries));
-          }
-
-        return unique_evaluation_point_results;
-      }
+  template <int n_components, int dim, int spacedim, typename VectorType>
+  inline std::vector<
+    typename FEPointEvaluation<n_components, dim>::gradient_type>
+  point_gradients(
+    const Utilities::MPI::RemotePointEvaluation<dim, spacedim> &cache,
+    const DoFHandler<dim, spacedim> &                           dof_handler,
+    const VectorType &                                          vector,
+    const EvaluationFlags::EvaluationFlags                      flags)
+  {
+    return internal::evaluate_at_points<
+      n_components,
+      dim,
+      spacedim,
+      VectorType,
+      typename FEPointEvaluation<n_components, dim>::gradient_type>(
+      cache,
+      dof_handler,
+      vector,
+      flags,
+      update_gradients,
+      dealii::EvaluationFlags::gradients,
+      [](const auto &evaluator, const auto &q) {
+        return evaluator.get_gradient(q);
+      });
   }
+
 #endif
 } // namespace VectorTools
 
index 1e15ae272751873a62e6b0c217b828fca26c81c1..c30f0ec9bc43c05fa373f9a07a2c5775f18c569c 100644 (file)
@@ -155,13 +155,17 @@ test()
   Utilities::MPI::RemotePointEvaluation<dim> evaluation_cache;
   const auto evaluation_point_results = VectorTools::point_values<1>(
     mapping_1, dof_handler_1, vector_1, evaluation_points, evaluation_cache);
+  const auto evaluation_point_gradient_results =
+    VectorTools::point_gradients<1>(
+      mapping_1, dof_handler_1, vector_1, evaluation_points, evaluation_cache);
   vector_1.zero_out_ghosts();
 
   if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
     for (unsigned int i = 0; i <= n_intervals; ++i)
       {
         if (std::abs(evaluation_points[i][0] - evaluation_point_results[i]) >
-            1e-10)
+              1e-10 ||
+            std::abs(evaluation_point_gradient_results[i][0] - 1.0) > 1e-10)
           Assert(false, ExcNotImplemented());
       }
 

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.