]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add factory functions for FERemoteEvaluationCommunicator
authorJohannes Heinz <43043310+jh66637@users.noreply.github.com>
Mon, 22 Jan 2024 08:52:34 +0000 (09:52 +0100)
committerJohannes Heinz <43043310+jh66637@users.noreply.github.com>
Sun, 28 Jan 2024 17:30:24 +0000 (18:30 +0100)
doc/news/changes/minor/20231231Heinz [new file with mode: 0644]
include/deal.II/matrix_free/fe_remote_evaluation.h
tests/matrix_free/fe_remote_evaluation_03.cc [new file with mode: 0644]
tests/matrix_free/fe_remote_evaluation_03.with_p4est=true.mpirun=1.with_cgal=true.output [new file with mode: 0644]
tests/matrix_free/fe_remote_evaluation_03.with_p4est=true.mpirun=2.with_cgal=true.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20231231Heinz b/doc/news/changes/minor/20231231Heinz
new file mode 100644 (file)
index 0000000..306b95c
--- /dev/null
@@ -0,0 +1,4 @@
+New: Add utility functions that setup FERemoteEvaluationCommunicator for point-to-point
+interpolation and Nitsche-type mortaring.
+<br>
+(Johannes Heinz, 2023/12/31)
index 11bf9f874a81cde786ca2aa1f5200cedefad68a7..7bf81d25f66d267ec007e41b6a2aa725a57cf20b 100644 (file)
 
 #include <deal.II/base/mpi_remote_point_evaluation.h>
 
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_tools_cache.h>
+
 #include <deal.II/matrix_free/fe_evaluation.h>
 #include <deal.II/matrix_free/fe_point_evaluation.h>
 #include <deal.II/matrix_free/matrix_free.h>
 
 #include <deal.II/numerics/vector_tools.h>
 
+#include <algorithm>
 #include <variant>
 
 DEAL_II_NAMESPACE_OPEN
@@ -446,6 +450,76 @@ private:
 
 
 
+/**
+ * The namespace collects convenience functions to setup
+ * FERemoteEvaluationCommunicator for typical use cases.
+ */
+namespace Utilities
+{
+  /**
+   * A factory function for the FERemoteEvaluationCommunicator in the case of
+   * point-to-point interpolation.
+   *
+   * @param[in] matrix_free MatrixFree object that is used to distribute
+   * quadrature points at non-matching faces.
+   * @param[in] non_matching_faces_marked_vertices A vector of boundary face IDs
+   * that relate to non-matching interfaces. Each boundary face ID is
+   * accompanied by a lambda that marks the vertices of cells which are
+   * considered during the search of remote points (quadrature points of faces
+   * with given boundary face ID).
+   * @param[in] quad_no Quadrature number in @p matrix_free.
+   * @param[in] dof_no DoFHandler number in @p matrix_free.
+   * @param[in] tolerance Tolerance to find remote points.
+   */
+  template <int dim,
+            typename Number,
+            typename VectorizedArrayType = VectorizedArray<Number>>
+  FERemoteEvaluationCommunicator<dim>
+  setup_remote_communicator_faces_point_to_point_interpolation(
+    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+    const std::vector<
+      std::pair<types::boundary_id, std::function<std::vector<bool>()>>>
+                      &non_matching_faces_marked_vertices,
+    const unsigned int quad_no   = 0,
+    const unsigned int dof_no    = 0,
+    const double       tolerance = 1e-9);
+
+
+
+  /**
+   * A factory function for the FERemoteEvaluationCommunicator in the case of
+   * Nitsche-type mortaring.
+   *
+   * @param[in] matrix_free MatrixFree object that is used to distribute
+   * @param[in] non_matching_faces_marked_vertices A vector of boundary face IDs
+   * that relate to non-matching interfaces. Each boundary face ID is
+   * accompanied by a lambda that marks the vertices of cells which are
+   * considered during the process of computing intersections.
+   * @param[in] n_q_pnts_1D The number of 1D quadrature points per intersection.
+   * @param[in] dof_no DoFHandler number in @p matrix_free.
+   * distributed to each intersection (given for one coordinate direction).
+   * @param[in] nm_mapping_info In case nm_mapping_info is not a `nullptr` it is
+   * setup for cell based loops, such that it can be used in combination with
+   * FERemoteEvaluation.
+   * @param[in] tolerance Tolerance to find intersections.
+   */
+  template <int dim,
+            typename Number,
+            typename VectorizedArrayType = VectorizedArray<Number>>
+  FERemoteEvaluationCommunicator<dim>
+  setup_remote_communicator_faces_nitsche_type_mortaring(
+    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+    const std::vector<
+      std::pair<types::boundary_id, std::function<std::vector<bool>()>>>
+                      &non_matching_faces_marked_vertices,
+    const unsigned int n_q_pnts_1D,
+    const unsigned int dof_no                                   = 0,
+    NonMatching::MappingInfo<dim, dim, Number> *nm_mapping_info = nullptr,
+    const double                                tolerance       = 1e-9);
+} // namespace Utilities
+
+
+
 /**
  * Class to access data in a matrix-free loop for non-matching discretizations.
  * Interfaces are named with FEEvaluation in mind.
@@ -999,6 +1073,373 @@ FERemoteEvaluationCommunicator<dim>::CopyInstructions::copy_data_entries(
 
 
 
+namespace Utilities
+{
+  template <int dim, typename Number, typename VectorizedArrayType>
+  FERemoteEvaluationCommunicator<dim>
+  setup_remote_communicator_faces_point_to_point_interpolation(
+    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+    const std::vector<
+      std::pair<types::boundary_id, std::function<std::vector<bool>()>>>
+                      &non_matching_faces_marked_vertices,
+    const unsigned int quad_no,
+    const unsigned int dof_no,
+    const double       tolerance)
+  {
+    const auto &dof_handler = matrix_free.get_dof_handler(dof_no);
+    const auto &tria        = dof_handler.get_triangulation();
+    const auto &mapping     = *matrix_free.get_mapping_info().mapping;
+
+    // Communication objects know about the communication pattern. I.e.,
+    // they know about the cells and quadrature points that have to be
+    // evaluated at remote faces. This information is given via
+    // RemotePointEvaluation. Additionally, the communication objects
+    // have to be able to match the quadrature points of the remote
+    // points (that provide exterior information) to the quadrature points
+    // defined at the interior cell. In case of point-to-point interpolation
+    // a vector of pairs with face batch Ids and the number of faces in the
+    // batch is needed. @c FERemoteCommunicationObjectEntityBatches
+    // is a container to store this information.
+    //
+    // We need multiple communication objects (one for each non-matching face
+    // ID).
+    std::vector<FERemoteCommunicationObjectEntityBatches<dim>> comm_objects;
+
+    // Additionally to the communication objects we need a vector
+    // that stores quadrature rules for every face batch.
+    // The quadrature can be empty in case of non non-matching faces,
+    // i.e. boundary faces. Internally this information is needed to correctly
+    // access values over multiple communication objects.
+    std::vector<Quadrature<dim>> global_quadrature_vector(
+      matrix_free.n_boundary_face_batches());
+
+    // Get the range of face batches we have to look at during construction of
+    // the communication objects. We only have to look at boundary faces.
+    const auto face_batch_range =
+      std::make_pair(matrix_free.n_inner_face_batches(),
+                     matrix_free.n_inner_face_batches() +
+                       matrix_free.n_boundary_face_batches());
+
+    // Iterate over all non-matching face IDs.
+    for (const auto &[nm_face, marked_vertices] :
+         non_matching_faces_marked_vertices)
+      {
+        // Construct the communication object for every face ID:
+        // 1) RemotePointEvaluation with user specified function for marked
+        // vertices.
+        auto rpe = std::make_shared<Utilities::MPI::RemotePointEvaluation<dim>>(
+          tolerance, false, 0, marked_vertices);
+
+        // 2) Face batch IDs and number of faces in batch.
+        std::vector<std::pair<unsigned int, unsigned int>>
+          face_batch_id_n_faces;
+
+        // Points that are searched by rpe.
+        std::vector<Point<dim>> points;
+
+        // Temporarily setup FEFaceEvaluation to access the quadrature points at
+        // the faces on the non-matching interface.
+        FEFaceEvaluation<dim, -1, 0, 1, Number> phi(matrix_free,
+                                                    true,
+                                                    dof_no,
+                                                    quad_no);
+
+        // Iterate over the boundary faces.
+        for (unsigned int bface = 0;
+             bface < face_batch_range.second - face_batch_range.first;
+             ++bface)
+          {
+            const unsigned int face = face_batch_range.first + bface;
+
+            if (matrix_free.get_boundary_id(face) == nm_face)
+              {
+                phi.reinit(face);
+
+                // If @c face is on the current side of the non-matching
+                // interface. Add the face batch ID and the number of faces in
+                // the batch to the corresponding data structure.
+                const unsigned int n_faces =
+                  matrix_free.n_active_entries_per_face_batch(face);
+                face_batch_id_n_faces.emplace_back(
+                  std::make_pair(face, n_faces));
+
+                // Append the quadrature points to the points we need to search
+                // for.
+                for (unsigned int v = 0; v < n_faces; ++v)
+                  {
+                    for (unsigned int q : phi.quadrature_point_indices())
+                      {
+                        const auto point = phi.quadrature_point(q);
+                        Point<dim> temp;
+                        for (unsigned int i = 0; i < dim; ++i)
+                          temp[i] = point[i][v];
+
+                        points.push_back(temp);
+                      }
+                  }
+
+                // Insert a quadrature rule of correct size into the global
+                // quadrature vector. First check that each face is only
+                // considered once.
+                Assert(global_quadrature_vector[bface].size() == 0,
+                       ExcMessage(
+                         "Quadrature for given face already provided."));
+
+                global_quadrature_vector[bface] =
+                  Quadrature<dim>(phi.n_q_points);
+              }
+          }
+
+        // Reinit RPE and ensure all points are found.
+        rpe->reinit(points, tria, mapping);
+        Assert(rpe->all_points_found(),
+               ExcMessage("Not all remote points found."));
+
+        // Add communication object to the list of objects.
+        FERemoteCommunicationObjectEntityBatches<dim> co;
+        co.batch_id_n_entities = face_batch_id_n_faces;
+        co.rpe                 = rpe;
+        comm_objects.push_back(co);
+      }
+
+    // Reinit the communicator `FERemoteEvaluationCommunicator`
+    // with the communication objects.
+    FERemoteEvaluationCommunicator<dim> remote_communicator;
+
+    remote_communicator.reinit_faces(comm_objects,
+                                     face_batch_range,
+                                     global_quadrature_vector);
+
+    return remote_communicator;
+  }
+
+
+
+  template <int dim, typename Number, typename VectorizedArrayType>
+  FERemoteEvaluationCommunicator<dim>
+  setup_remote_communicator_faces_nitsche_type_mortaring(
+    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+    const std::vector<
+      std::pair<types::boundary_id, std::function<std::vector<bool>()>>>
+                      &non_matching_faces_marked_vertices,
+    const unsigned int n_q_pnts_1D,
+    const unsigned int dof_no,
+    NonMatching::MappingInfo<dim, dim, Number> *nm_mapping_info,
+    const double                                tolerance)
+  {
+    const auto &dof_handler = matrix_free.get_dof_handler(dof_no);
+    const auto &tria        = dof_handler.get_triangulation();
+    const auto &mapping     = *matrix_free.get_mapping_info().mapping;
+
+    constexpr unsigned int n_lanes = VectorizedArray<Number>::size();
+
+    std::pair<unsigned int, unsigned int> face_range =
+      std::make_pair(matrix_free.n_inner_face_batches(),
+                     matrix_free.n_inner_face_batches() +
+                       matrix_free.n_boundary_face_batches());
+
+    std::vector<Quadrature<dim - 1>> global_quadrature_vector(
+      (matrix_free.n_inner_face_batches() +
+       matrix_free.n_boundary_face_batches()) *
+      n_lanes);
+
+    // In case of Nitsche-type mortaring a vector of face indices is
+    // needed as communication object.
+    // @c FERemoteCommunicationObjectFaces is a container to store this
+    // information.
+    //
+    // We need multiple communication objects (one for each non-matching face
+    // ID).
+    std::vector<FERemoteCommunicationObject<dim>> comm_objects;
+
+    // Create bounding boxes and GridTools::Cache which is needed in
+    // the following loop.
+    std::vector<BoundingBox<dim>> local_boxes;
+    for (const auto &cell : tria.active_cell_iterators())
+      if (cell->is_locally_owned())
+        local_boxes.emplace_back(mapping.get_bounding_box(cell));
+
+    // Create r-tree of bounding boxes
+    const auto local_tree = pack_rtree(local_boxes);
+
+    // Compress r-tree to a minimal set of bounding boxes
+    std::vector<std::vector<BoundingBox<dim>>> global_bboxes(1);
+    global_bboxes[0] = extract_rtree_level(local_tree, 0);
+
+    const GridTools::Cache<dim, dim> cache(tria, mapping);
+
+    // Iterate over all sides of the non-matching interface.
+    for (const auto &[nm_face, marked_vertices] :
+         non_matching_faces_marked_vertices)
+      {
+        // 1) compute cell face pairs
+        std::vector<
+          std::pair<typename Triangulation<dim>::cell_iterator, unsigned int>>
+          cell_face_pairs;
+
+        std::vector<unsigned int> indices;
+
+        for (unsigned int face = face_range.first; face < face_range.second;
+             ++face)
+          {
+            if (matrix_free.get_boundary_id(face) == nm_face)
+              {
+                for (unsigned int v = 0;
+                     v < matrix_free.n_active_entries_per_face_batch(face);
+                     ++v)
+                  {
+                    const auto &[c, f] = matrix_free.get_face_iterator(face, v);
+
+                    cell_face_pairs.emplace_back(std::make_pair(c, f));
+                    indices.push_back(face * n_lanes + v);
+                  }
+              }
+          }
+
+        // 2) Create RPE.
+        // In the Nitsche-type case we do not collect points for the setup
+        // of RemotePointEvaluation. Instead we compute intersections between
+        // the faces and setup RemotePointEvaluation with the computed
+        // intersections.
+
+        // Build intersection requests. Intersection requests
+        // correspond to vertices at faces.
+        std::vector<std::vector<Point<dim>>> intersection_requests;
+        for (const auto &[cell, f] : cell_face_pairs)
+          {
+            std::vector<Point<dim>> vertices(cell->face(f)->n_vertices());
+            std::copy_n(mapping.get_vertices(cell, f).begin(),
+                        cell->face(f)->n_vertices(),
+                        vertices.begin());
+            intersection_requests.emplace_back(vertices);
+          }
+
+        // Compute intersection data with user specified function for marked
+        // vertices.
+        auto intersection_data =
+          GridTools::internal::distributed_compute_intersection_locations<
+            dim - 1>(cache,
+                     intersection_requests,
+                     global_bboxes,
+                     marked_vertices(),
+                     tolerance);
+
+        // Convert to RPE.
+        std::vector<Quadrature<dim>> mapped_quadratures_recv_comp;
+
+        auto rpe =
+          std::make_shared<Utilities::MPI::RemotePointEvaluation<dim>>();
+        rpe->reinit(
+          intersection_data
+            .template convert_to_distributed_compute_point_locations_internal<
+              dim>(n_q_pnts_1D, tria, mapping, &mapped_quadratures_recv_comp),
+          tria,
+          mapping);
+
+        // 3) Fill global quadrature vector.
+        for (unsigned int i = 0; i < intersection_requests.size(); ++i)
+          {
+            const auto idx = indices[i];
+
+            // We do not use a structural binding here, since with
+            // C++17 caputuring structural bindings in lambdas leads
+            // to an ill formed program.
+            const auto &cell = std::get<0>(cell_face_pairs[i]);
+            const auto &f    = std::get<1>(cell_face_pairs[i]);
+
+            std::vector<Point<dim - 1>> q_points;
+            std::vector<double>         weights;
+            for (unsigned int ptr = intersection_data.recv_ptrs[i];
+                 ptr < intersection_data.recv_ptrs[i + 1];
+                 ++ptr)
+              {
+                const auto &quad = mapped_quadratures_recv_comp[ptr];
+
+                const auto &ps = quad.get_points();
+                std::transform(
+                  ps.begin(),
+                  ps.end(),
+                  std::back_inserter(q_points),
+                  [&](const Point<dim> &p) {
+                    return mapping.project_real_point_to_unit_point_on_face(
+                      cell, f, p);
+                  });
+
+                const auto &ws = quad.get_weights();
+                weights.insert(weights.end(), ws.begin(), ws.end());
+              }
+            Quadrature<dim - 1> quad(q_points, weights);
+
+            Assert(global_quadrature_vector[idx].size() == 0,
+                   ExcMessage("Quadrature for given face already provided."));
+
+            global_quadrature_vector[idx] = quad;
+          }
+
+        // Add communication object.
+        FERemoteCommunicationObject<dim> co;
+        co.indices = indices;
+        co.rpe     = rpe;
+        comm_objects.push_back(co);
+      }
+
+    // Reinit the communicator with the communication objects.
+    FERemoteEvaluationCommunicator<dim> remote_communicator;
+
+    remote_communicator.reinit_faces(
+      comm_objects,
+      std::make_pair(0, global_quadrature_vector.size()),
+      global_quadrature_vector);
+
+    if (nm_mapping_info != nullptr)
+      {
+        std::vector<
+          std::pair<typename DoFHandler<dim>::cell_iterator, unsigned int>>
+          vector_face_accessors;
+        vector_face_accessors.reserve((matrix_free.n_inner_face_batches() +
+                                       matrix_free.n_boundary_face_batches()) *
+                                      n_lanes);
+
+        // fill container for inner face batches
+        unsigned int face_batch = 0;
+        for (; face_batch < matrix_free.n_inner_face_batches(); ++face_batch)
+          {
+            for (unsigned int v = 0; v < n_lanes; ++v)
+              {
+                if (v < matrix_free.n_active_entries_per_face_batch(face_batch))
+                  vector_face_accessors.push_back(
+                    matrix_free.get_face_iterator(face_batch, v));
+                else
+                  vector_face_accessors.push_back(
+                    matrix_free.get_face_iterator(face_batch, 0));
+              }
+          }
+        // and boundary face batches
+        for (; face_batch < (matrix_free.n_inner_face_batches() +
+                             matrix_free.n_boundary_face_batches());
+             ++face_batch)
+          {
+            for (unsigned int v = 0; v < n_lanes; ++v)
+              {
+                if (v < matrix_free.n_active_entries_per_face_batch(face_batch))
+                  vector_face_accessors.push_back(
+                    matrix_free.get_face_iterator(face_batch, v));
+                else
+                  vector_face_accessors.push_back(
+                    matrix_free.get_face_iterator(face_batch, 0));
+              }
+          }
+
+        nm_mapping_info->reinit_faces(vector_face_accessors,
+                                      global_quadrature_vector);
+      }
+
+    return remote_communicator;
+  }
+} // namespace Utilities
+
+
+
 template <int dim, int n_components, typename value_type>
 template <typename MeshType>
 FERemoteEvaluation<dim, n_components, value_type>::FERemoteEvaluation(
diff --git a/tests/matrix_free/fe_remote_evaluation_03.cc b/tests/matrix_free/fe_remote_evaluation_03.cc
new file mode 100644 (file)
index 0000000..9d12be1
--- /dev/null
@@ -0,0 +1,308 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+// This test checks the correct setup of FERemoteEvaluationCommunicator
+// via the factory functions
+// setup_remote_communicator_faces_point_to_point_interpolation()
+// and setup_remote_communicator_faces_nitsche_type_mortaring().
+
+#include <deal.II/base/mpi.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_system.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/la_parallel_block_vector.h>
+
+#include <deal.II/matrix_free/fe_remote_evaluation.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+using Number = double;
+
+using VectorType = LinearAlgebra::distributed::Vector<Number>;
+
+constexpr unsigned int dim = 2;
+
+template <int dim>
+struct Coordinates : public Function<dim>
+{
+  Coordinates()
+    : Function<dim>(dim, 0.0)
+  {}
+
+  double
+  value(const Point<dim> &p, const unsigned int component) const final
+  {
+    if (component == 0)
+      return p[component];
+    return p[component];
+  }
+};
+
+
+void
+test_point_to_point(const MatrixFree<dim, Number> &matrix_free)
+{
+  const Triangulation<dim> &tria =
+    matrix_free.get_dof_handler().get_triangulation();
+
+  std::vector<std::pair<types::boundary_id, std::function<std::vector<bool>()>>>
+    non_matching_faces_marked_vertices;
+
+  non_matching_faces_marked_vertices.emplace_back(std::make_pair(1, [&]() {
+    std::vector<bool> mask(tria.n_vertices(), true);
+
+    for (const auto &cell : tria.active_cell_iterators())
+      for (auto const &f : cell->face_indices())
+        if (cell->face(f)->at_boundary() && cell->face(f)->boundary_id() == 1)
+          for (const auto v : cell->vertex_indices())
+            mask[cell->vertex_index(v)] = false;
+
+    return mask;
+  }));
+
+  auto remote_communicator =
+    Utilities::setup_remote_communicator_faces_point_to_point_interpolation(
+      matrix_free, non_matching_faces_marked_vertices);
+
+  // 1) Allocate vectors and interpolate grid coordinates in DoFs.
+  VectorType src;
+  VectorType dst;
+  matrix_free.initialize_dof_vector(src);
+  matrix_free.initialize_dof_vector(dst);
+
+  VectorTools::interpolate(*matrix_free.get_mapping_info().mapping,
+                           matrix_free.get_dof_handler(),
+                           Coordinates<dim>(),
+                           src);
+
+  // 2) Construct remote evaluator...
+  FERemoteEvaluation<dim, dim, VectorizedArray<Number>> phi_remote_eval(
+    remote_communicator, matrix_free.get_dof_handler());
+
+  // ...and precompute remote values.
+  phi_remote_eval.gather_evaluate(src, EvaluationFlags::values);
+
+  // 3) Access values at quadrature points
+  const auto boundary_function =
+    [&](const auto &data, auto &dst, const auto &src, const auto face_range) {
+      FEFaceEvaluation<dim, -1, 0, dim, Number> phi(matrix_free, true);
+      auto phi_remote = phi_remote_eval.get_data_accessor();
+
+      for (unsigned int face = face_range.first; face < face_range.second;
+           ++face)
+        {
+          if (matrix_free.get_boundary_id(face) == 1)
+            {
+              phi.reinit(face);
+              phi_remote.reinit(face);
+
+              for (unsigned int q : phi.quadrature_point_indices())
+                {
+                  const auto error =
+                    phi.quadrature_point(q) - phi_remote.get_value(q);
+
+                  // only consider active entries in global error
+                  for (unsigned int d = 0; d < dim; ++d)
+                    for (unsigned int v = 0;
+                         v < matrix_free.n_active_entries_per_face_batch(face);
+                         ++v)
+                      AssertThrow(std::abs(error[d][v]) < 1e-13,
+                                  ExcMessage("Error too large."));
+                }
+            }
+        }
+    };
+
+  matrix_free.template loop<VectorType, VectorType>(
+    {}, {}, boundary_function, dst, src, true);
+}
+
+void
+test_nitsche_type_mortaring(const MatrixFree<dim, Number> &matrix_free)
+{
+  const Triangulation<dim> &tria =
+    matrix_free.get_dof_handler().get_triangulation();
+
+  std::vector<std::pair<types::boundary_id, std::function<std::vector<bool>()>>>
+    non_matching_faces_marked_vertices;
+
+  non_matching_faces_marked_vertices.emplace_back(std::make_pair(1, [&]() {
+    std::vector<bool> mask(tria.n_vertices(), true);
+
+    for (const auto &cell : tria.active_cell_iterators())
+      for (auto const &f : cell->face_indices())
+        if (cell->face(f)->at_boundary() && cell->face(f)->boundary_id() == 1)
+          for (const auto v : cell->vertex_indices())
+            mask[cell->vertex_index(v)] = false;
+
+    return mask;
+  }));
+
+  typename NonMatching::MappingInfo<dim, dim, Number>::AdditionalData
+    additional_data;
+  additional_data.use_global_weights = true;
+
+  NonMatching::MappingInfo<dim, dim, Number> nm_mapping_info(
+    *matrix_free.get_mapping_info().mapping,
+    update_quadrature_points,
+    additional_data);
+
+  auto remote_communicator =
+    Utilities::setup_remote_communicator_faces_nitsche_type_mortaring(
+      matrix_free, non_matching_faces_marked_vertices, 4, 0, &nm_mapping_info);
+
+  // 1) Allocate vectors and interpolate grid coordinates in DoFs.
+  VectorType src;
+  VectorType dst;
+  matrix_free.initialize_dof_vector(src);
+  matrix_free.initialize_dof_vector(dst);
+
+  VectorTools::interpolate(*matrix_free.get_mapping_info().mapping,
+                           matrix_free.get_dof_handler(),
+                           Coordinates<dim>(),
+                           src);
+
+  // 2) Construct remote evaluator...
+  FERemoteEvaluation<dim, dim, Number> phi_remote_eval(
+    remote_communicator, matrix_free.get_dof_handler());
+
+  // ...and precompute remote values.
+  phi_remote_eval.gather_evaluate(src, EvaluationFlags::values);
+
+  // 3) Access values at quadrature points
+  const auto boundary_function =
+    [&](const auto &data, auto &dst, const auto &src, const auto face_range) {
+      FEFacePointEvaluation<dim, dim, dim, Number> phi(
+        nm_mapping_info, matrix_free.get_dof_handler().get_fe());
+      auto phi_remote = phi_remote_eval.get_data_accessor();
+
+
+      for (unsigned int face = face_range.first; face < face_range.second;
+           ++face)
+        {
+          if (matrix_free.get_boundary_id(face) == 1)
+            {
+              for (unsigned int v = 0;
+                   v < matrix_free.n_active_entries_per_face_batch(face);
+                   ++v)
+                {
+                  constexpr unsigned int n_lanes =
+                    VectorizedArray<Number>::size();
+                  phi.reinit(face * n_lanes + v);
+                  phi_remote.reinit(face * n_lanes + v);
+
+                  for (unsigned int q : phi.quadrature_point_indices())
+                    {
+                      const auto             point = phi.real_point(q);
+                      Tensor<1, dim, Number> temp;
+                      for (unsigned int i = 0; i < dim; ++i)
+                        temp[i] = point[i];
+
+                      const auto error = temp - phi_remote.get_value(q);
+                      for (unsigned int d = 0; d < dim; ++d)
+                        AssertThrow(std::abs(error[d]) < 1e-13,
+                                    ExcMessage("Error too large."));
+                    }
+                }
+            }
+        }
+    };
+
+  matrix_free.template loop<VectorType, VectorType>(
+    {}, {}, boundary_function, dst, src, true);
+}
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
+  MPILogInitAll                    all;
+
+  // Construct non-matching triangulation.
+  const double       length = 1.0;
+  Triangulation<dim> tria_left;
+  const unsigned int subdiv_left = 11;
+  GridGenerator::subdivided_hyper_rectangle(tria_left,
+                                            {subdiv_left, 2 * subdiv_left},
+                                            {0.0, 0.0},
+                                            {0.525 * length, length});
+
+  for (const auto &face : tria_left.active_face_iterators())
+    if (face->at_boundary())
+      {
+        face->set_boundary_id(0);
+        if (face->center()[0] > 0.525 * length - 1e-6)
+          face->set_boundary_id(1);
+      }
+
+  Triangulation<dim> tria_right;
+  const unsigned int subdiv_right = 4;
+  GridGenerator::subdivided_hyper_rectangle(tria_right,
+                                            {subdiv_right, 2 * subdiv_right},
+                                            {0.525 * length, 0.0},
+                                            {length, length});
+
+  for (const auto &face : tria_right.active_face_iterators())
+    if (face->at_boundary())
+      face->set_boundary_id(0);
+
+  // Merge triangulations with tolerance 0 to ensure no vertices
+  // are merged.
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+  GridGenerator::merge_triangulations(tria_left,
+                                      tria_right,
+                                      tria,
+                                      /*tolerance*/ 0.,
+                                      /*copy_manifold_ids*/ false,
+                                      /*copy_boundary_ids*/ true);
+
+  // Create DoFHandler
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(FESystem<dim>(FE_DGQ<dim>(3), dim));
+
+  // Constraints are not considered in this test
+  AffineConstraints<Number> constraints;
+  constraints.close();
+
+  // Setup MatrixFree
+  typename MatrixFree<dim, Number>::AdditionalData data;
+  data.mapping_update_flags = update_values;
+  // We are not considering inner faces in this test. Therefore,
+  // data.mapping_update_flags_inner_faces is not set.
+  data.mapping_update_flags_boundary_faces =
+    update_quadrature_points | update_values;
+
+
+  MatrixFree<dim, Number> matrix_free;
+  matrix_free.reinit(
+    MappingQ1<dim>(), dof_handler, constraints, QGauss<dim>(4), data);
+
+  // 2) test point-to-point
+  test_point_to_point(matrix_free);
+
+  deallog << "OK" << std::endl;
+
+  // 3) test Nitsche-type
+  test_nitsche_type_mortaring(matrix_free);
+
+  deallog << "OK" << std::endl;
+
+  return 0;
+}
diff --git a/tests/matrix_free/fe_remote_evaluation_03.with_p4est=true.mpirun=1.with_cgal=true.output b/tests/matrix_free/fe_remote_evaluation_03.with_p4est=true.mpirun=1.with_cgal=true.output
new file mode 100644 (file)
index 0000000..2d64160
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL:0::OK
+DEAL:0::OK
diff --git a/tests/matrix_free/fe_remote_evaluation_03.with_p4est=true.mpirun=2.with_cgal=true.output b/tests/matrix_free/fe_remote_evaluation_03.with_p4est=true.mpirun=2.with_cgal=true.output
new file mode 100644 (file)
index 0000000..abeeb87
--- /dev/null
@@ -0,0 +1,7 @@
+
+DEAL:0::OK
+DEAL:0::OK
+
+DEAL:1::OK
+DEAL:1::OK
+

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.