]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement communication-free owner rank search (p4est interface). 12341/head
authorSimon, Dr. Konrad <konrad.simon@uni-hamburg.de>
Sun, 23 May 2021 08:40:37 +0000 (10:40 +0200)
committerKonrad Simon <konrad.simon@uni-hamburg.de>
Wed, 13 Oct 2021 19:12:10 +0000 (21:12 +0200)
26 files changed:
cmake/configure/configure_p4est.cmake
cmake/modules/FindP4EST.cmake
doc/news/changes/major/20210924Simon [new file with mode: 0644]
include/deal.II/distributed/p4est_wrappers.h
include/deal.II/distributed/tria.h
source/distributed/p4est_wrappers.cc
source/distributed/tria.cc
tests/distributed_grids/coarse_grid_common.h
tests/mpi/codim_01.cc
tests/mpi/p4est_2d_find_point_owner_rank_01.cc [new file with mode: 0644]
tests/mpi/p4est_2d_find_point_owner_rank_01.mpirun=1.with_p4est=true.with_p4est_with_search_local=true.output [new file with mode: 0644]
tests/mpi/p4est_2d_find_point_owner_rank_01.mpirun=3.with_p4est=true.with_p4est_with_search_local=true.output [new file with mode: 0644]
tests/mpi/p4est_2d_find_point_owner_rank_01.mpirun=5.with_p4est=true.with_p4est_with_search_local=true.output [new file with mode: 0644]
tests/mpi/p4est_2d_find_point_owner_rank_02.cc [new file with mode: 0644]
tests/mpi/p4est_2d_find_point_owner_rank_02.mpirun=1.with_p4est=true.with_p4est_with_search_local=true.output [new file with mode: 0644]
tests/mpi/p4est_2d_find_point_owner_rank_02.mpirun=3.with_p4est=true.with_p4est_with_search_local=true.output [new file with mode: 0644]
tests/mpi/p4est_2d_find_point_owner_rank_02.mpirun=5.with_p4est=true.with_p4est_with_search_local=true.output [new file with mode: 0644]
tests/mpi/p4est_3d_find_point_owner_rank_01.cc [new file with mode: 0644]
tests/mpi/p4est_3d_find_point_owner_rank_01.mpirun=1.with_p4est=true.with_p4est_with_search_local=true.output [new file with mode: 0644]
tests/mpi/p4est_3d_find_point_owner_rank_01.mpirun=3.with_p4est=true.with_p4est_with_search_local=true.output [new file with mode: 0644]
tests/mpi/p4est_3d_find_point_owner_rank_01.mpirun=5.with_p4est=true.with_p4est_with_search_local=true.output [new file with mode: 0644]
tests/mpi/p4est_3d_find_point_owner_rank_02.cc [new file with mode: 0644]
tests/mpi/p4est_3d_find_point_owner_rank_02.mpirun=1.with_p4est=true.with_p4est_with_search_local=true.output [new file with mode: 0644]
tests/mpi/p4est_3d_find_point_owner_rank_02.mpirun=3.with_p4est=true.with_p4est_with_search_local=true.output [new file with mode: 0644]
tests/mpi/p4est_3d_find_point_owner_rank_02.mpirun=5.with_p4est=true.with_p4est_with_search_local=true.output [new file with mode: 0644]
tests/mpi/p4est_find_point_owner_rank.h [new file with mode: 0644]

index cc9c3c00164011997a95d0af18bff578661d3ef8..1f0465aca18f84a3880a97ce93a47d893dba7279 100644 (file)
@@ -77,6 +77,7 @@ ENDMACRO()
 
 MACRO(FEATURE_P4EST_CONFIGURE_EXTERNAL)
   SET(DEAL_II_P4EST_WITH_VTK_BINARY ${P4EST_WITH_VTK_BINARY})
+  SET(DEAL_II_P4EST_WITH_SEARCH_LOCAL ${P4EST_WITH_SEARCH_LOCAL})
 ENDMACRO()
 
 
index eb87d2f47e15af1b21aa9e69f186a8c66d39d0ad..a91ef26876d180163a6a1ac201379103ac6fd669 100644 (file)
@@ -22,6 +22,7 @@
 #   P4EST_WITH_MPI
 #   P4EST_WITH_ZLIB
 #   P4EST_WITH_VTK_BINARY
+#   P4EST_WITH_SEARCH_LOCAL
 #   P4EST_VERSION
 #   P4EST_VERSION_MAJOR
 #   P4EST_VERSION_MINOR
@@ -139,6 +140,17 @@ IF(EXISTS ${P4EST_INCLUDE_DIR}/p4est_config.h)
   ELSE()
     SET(P4EST_WITH_VTK_BINARY TRUE)
   ENDIF()
+  
+  #
+  # Does p4est have search local?
+  #
+  FILE(STRINGS "${P4EST_INCLUDE_DIR}/p4est_base.h" P4EST_SEARCH_LOCAL_STRING
+    REGEX "#define.*P4EST_SEARCH_LOCAL")
+  IF("${P4EST_SEARCH_LOCAL_STRING}" STREQUAL "")
+    SET(P4EST_WITH_SEARCH_LOCAL FALSE)
+  ELSE()
+    SET(P4EST_WITH_SEARCH_LOCAL TRUE)
+  ENDIF()
 
   #
   # Extract version numbers:
diff --git a/doc/news/changes/major/20210924Simon b/doc/news/changes/major/20210924Simon
new file mode 100644 (file)
index 0000000..04418f7
--- /dev/null
@@ -0,0 +1,7 @@
+New: Added an interface to p4est_search.h (>=v2.2) to find the MPI
+rank of cells owning remote points without communication. This
+works in the special (but not uncommon) case that no nontrivial 
+manifold is attached to the triangulation, i.e., the triangulation
+has only FlatManifold attached to it.
+<br>
+(Konrad Simon, 2021/09/24)
index 21410aa792d8549e8590cc3fba05d99c6b7d0dd6..bf63f4369778650863496f370be663939d4aaaff 100644 (file)
 #  include <p4est_extended.h>
 #  include <p4est_ghost.h>
 #  include <p4est_iterate.h>
+#  include <p4est_search.h>
 #  include <p4est_vtk.h>
 #  include <p8est_bits.h>
 #  include <p8est_communication.h>
 #  include <p8est_extended.h>
 #  include <p8est_ghost.h>
 #  include <p8est_iterate.h>
+#  include <p8est_search.h>
 #  include <p8est_vtk.h>
 
 #  include <map>
@@ -87,12 +89,16 @@ namespace internal
       using forest           = p4est_t;
       using tree             = p4est_tree_t;
       using quadrant         = p4est_quadrant_t;
+      using quadrant_coord   = p4est_qcoord_t;
       using topidx           = p4est_topidx_t;
       using locidx           = p4est_locidx_t;
       using gloidx           = p4est_gloidx_t;
       using balance_type     = p4est_connect_type_t;
       using ghost            = p4est_ghost_t;
       using transfer_context = p4est_transfer_context_t;
+#  ifdef P4EST_SEARCH_LOCAL
+      using search_partition_callback = p4est_search_partition_t;
+#  endif
     };
 
     template <>
@@ -102,12 +108,16 @@ namespace internal
       using forest           = p8est_t;
       using tree             = p8est_tree_t;
       using quadrant         = p8est_quadrant_t;
+      using quadrant_coord   = p4est_qcoord_t;
       using topidx           = p4est_topidx_t;
       using locidx           = p4est_locidx_t;
       using gloidx           = p4est_gloidx_t;
       using balance_type     = p8est_connect_type_t;
       using ghost            = p8est_ghost_t;
       using transfer_context = p8est_transfer_context_t;
+#  ifdef P4EST_SEARCH_LOCAL
+      using search_partition_callback = p8est_search_partition_t;
+#  endif
     };
 
 
@@ -306,6 +316,22 @@ namespace internal
         const int *             src_sizes);
 
       static void (&transfer_custom_end)(types<2>::transfer_context *tc);
+
+#  ifdef P4EST_SEARCH_LOCAL
+      static void (&search_partition)(
+        types<2>::forest *                  forest,
+        int                                 call_post,
+        types<2>::search_partition_callback quadrant_fn,
+        types<2>::search_partition_callback point_fn,
+        sc_array_t *                        points);
+#  endif
+
+      static void (&quadrant_coord_to_vertex)(
+        types<2>::connectivity * connectivity,
+        types<2>::topidx         treeid,
+        types<2>::quadrant_coord x,
+        types<2>::quadrant_coord y,
+        double                   vxyz[3]);
     };
 
 
@@ -490,6 +516,23 @@ namespace internal
         const int *             src_sizes);
 
       static void (&transfer_custom_end)(types<3>::transfer_context *tc);
+
+#  ifdef P4EST_SEARCH_LOCAL
+      static void (&search_partition)(
+        types<3>::forest *                  forest,
+        int                                 call_post,
+        types<3>::search_partition_callback quadrant_fn,
+        types<3>::search_partition_callback point_fn,
+        sc_array_t *                        points);
+#  endif
+
+      static void (&quadrant_coord_to_vertex)(
+        types<3>::connectivity * connectivity,
+        types<3>::topidx         treeid,
+        types<3>::quadrant_coord x,
+        types<3>::quadrant_coord y,
+        types<3>::quadrant_coord z,
+        double                   vxyz[3]);
     };
 
 
index 51f22b9712dbc569eca2fbf3139d420f61870efe..dd4f34e903b11611979b57325add715a4da2f525 100644 (file)
@@ -50,7 +50,6 @@
 #  include <p8est_ghost.h>
 #endif
 
-
 DEAL_II_NAMESPACE_OPEN
 
 #ifdef DEAL_II_WITH_P4EST
@@ -324,7 +323,15 @@ namespace parallel
          * after a refinement cycle. It can be executed manually by calling
          * repartition().
          */
-        no_automatic_repartitioning = 0x4
+        no_automatic_repartitioning = 0x4,
+        /**
+         * Setting this flag will communicate vertices to p4est. In debug mode,
+         * the vertices will always be communicated. This way one can
+         * use the 'find_point_owner_rank()' to find the MPI rank of the active
+         * cell that owns an arbitrary point in case all attached manifolds are
+         * flat.
+         */
+        communicate_vertices_to_p4est = 0x8
       };
 
 
@@ -386,6 +393,12 @@ namespace parallel
       bool
       is_multilevel_hierarchy_constructed() const override;
 
+      /**
+       * Return if vertices will be communicated to p4est in release mode.
+       */
+      bool
+      are_vertices_communicated_to_p4est() const;
+
       /**
        * Transfer data across forests.
        *
@@ -441,6 +454,36 @@ namespace parallel
         const TriangulationDescription::Description<dim, spacedim>
           &construction_data) override;
 
+      /**
+       * Find the MPI rank of the cell that contains this point in a distributed
+       * mesh.
+       *
+       * @note This function calls `find_point_owner_rank(const std::vector<Point<dim>> &points)`
+       * (requires p4est v2.2 and higher). Please see the documentation of
+       * `find_point_owner_rank(const std::vector<Point<dim>> &points)`.
+       */
+      types::subdomain_id
+      find_point_owner_rank(const Point<dim> &p);
+
+      /**
+       * Find the MPI rank of the cells that contain the input points in a
+       * distributed mesh. If any point is not owned by any mesh cell its return
+       * value will be `numbers::invalid_subdomain_id`.
+       *
+       * @note The query points do not need to be owned locally or in the ghost layer.
+       *
+       * @note This function can only be used with p4est v2.2 and higher, flat manifolds
+       * and requires the settings flag
+       * `Settings::communicate_vertices_to_p4est` to be set.
+       *
+       * @note The algorithm is free of communication.
+       *
+       * @param[in] points a list of query points
+       * @return list of owner ranks
+       */
+      std::vector<types::subdomain_id>
+      find_point_owner_rank(const std::vector<Point<dim>> &points);
+
       /**
        * Coarsen and refine the mesh according to refinement and coarsening
        * flags set.
@@ -564,6 +607,9 @@ namespace parallel
        *
        * More than anything else, this function is useful for debugging the
        * interface between deal.II and p4est.
+       *
+       * @note To use the function the flag
+       * `Settings::communicate_vertices_to_p4est` must be set.
        */
       void
       write_mesh_vtk(const std::string &file_basename) const;
@@ -896,6 +942,13 @@ namespace parallel
       virtual bool
       is_multilevel_hierarchy_constructed() const override;
 
+      /**
+       * This function is not implemented, but needs to be present for the
+       * compiler.
+       */
+      bool
+      are_vertices_communicated_to_p4est() const;
+
       /**
        * This function is not implemented, but needs to be present for the
        * compiler.
@@ -972,7 +1025,8 @@ namespace parallel
         default_setting                          = 0x0,
         mesh_reconstruction_after_repartitioning = 0x1,
         construct_multigrid_hierarchy            = 0x2,
-        no_automatic_repartitioning              = 0x4
+        no_automatic_repartitioning              = 0x4,
+        communicate_vertices_to_p4est            = 0x8
       };
 
       /**
@@ -996,6 +1050,16 @@ namespace parallel
         return false;
       }
 
+      /**
+       * Dummy replacement to allow for better error messages when compiling
+       * this class.
+       */
+      bool
+      are_vertices_communicated_to_p4est() const
+      {
+        return false;
+      }
+
       /**
        * Dummy replacement to allow for better error messages when compiling
        * this class.
index 996a0fcc3256d1eabf4378dd16c673873ec88774..85540da99b764f54979f1a30611a564404ed430e 100644 (file)
@@ -535,7 +535,21 @@ namespace internal
     void (&functions<2>::transfer_custom_end)(types<2>::transfer_context *tc) =
       p4est_transfer_custom_end;
 
-
+#  ifdef P4EST_SEARCH_LOCAL
+    void (&functions<2>::search_partition)(
+      types<2>::forest *                  p4est,
+      int                                 call_post,
+      types<2>::search_partition_callback quadrant_fn,
+      types<2>::search_partition_callback point_fn,
+      sc_array_t *                        points) = p4est_search_partition;
+#  endif
+
+    void (&functions<2>::quadrant_coord_to_vertex)(
+      types<2>::connectivity * connectivity,
+      types<2>::topidx         treeid,
+      types<2>::quadrant_coord x,
+      types<2>::quadrant_coord y,
+      double                   vxyz[3]) = p4est_qcoord_to_vertex;
 
     int (&functions<3>::quadrant_compare)(const void *v1, const void *v2) =
       p8est_quadrant_compare;
@@ -743,7 +757,22 @@ namespace internal
     void (&functions<3>::transfer_custom_end)(types<3>::transfer_context *tc) =
       p8est_transfer_custom_end;
 
-
+#  ifdef P4EST_SEARCH_LOCAL
+    void (&functions<3>::search_partition)(
+      types<3>::forest *                  p4est,
+      int                                 call_post,
+      types<3>::search_partition_callback quadrant_fn,
+      types<3>::search_partition_callback point_fn,
+      sc_array_t *                        points) = p8est_search_partition;
+#  endif
+
+    void (&functions<3>::quadrant_coord_to_vertex)(
+      types<3>::connectivity * connectivity,
+      types<3>::topidx         treeid,
+      types<3>::quadrant_coord x,
+      types<3>::quadrant_coord y,
+      types<3>::quadrant_coord z,
+      double                   vxyz[3]) = p8est_qcoord_to_vertex;
 
     template <int dim>
     void
index 648aab7efbcda4345725f2d61c81e440c817e8a6..e53164f288be976c69c9e54d7389fa8cc3265110 100644 (file)
@@ -16,6 +16,7 @@
 
 #include <deal.II/base/logstream.h>
 #include <deal.II/base/memory_consumption.h>
+#include <deal.II/base/point.h>
 #include <deal.II/base/utilities.h>
 
 #include <deal.II/distributed/p4est_wrappers.h>
@@ -567,6 +568,583 @@ namespace
   }
 
 
+#  ifdef P4EST_SEARCH_LOCAL
+  template <int dim>
+  class PartitionSearch
+  {
+  public:
+    PartitionSearch()
+    {
+      Assert(dim > 1, ExcNotImplemented());
+    }
+
+    PartitionSearch(const PartitionSearch<dim> &other) = delete;
+
+    PartitionSearch<dim> &
+    operator=(const PartitionSearch<dim> &other) = delete;
+
+  public:
+    /**
+     * Callback exectuted before point function. Last argument is always
+     * nullptr.
+     *
+     * @return `int` interpreted as a C "bool". Zero means "stop the recursion".
+     *
+     * @note We never stop the recursion in this callback since we search for
+     * each point individually.
+     */
+    static int
+    local_quadrant_fn(typename internal::p4est::types<dim>::forest *forest,
+                      typename internal::p4est::types<dim>::topidx  which_tree,
+                      typename internal::p4est::types<dim>::quadrant *quadrant,
+                      int   rank_begin,
+                      int   rank_end,
+                      void *point);
+
+    /**
+     * Callback for point function. Check whether a point is in a (physical)
+     * quadrant.
+     *
+     * @note We can handle a quadrant that is mapped by bi-linear or tri-linear
+     * mappings. Checking for a point in a cell of a curved domain required
+     * knowledge of the attached manifold.
+     *
+     * @return `int` interpreted as a C "bool". Zero means "stop the recursion".
+     * This can happen once we know the owner rank or if we know that a point
+     * does not belong to a quadrant.
+     */
+    static int
+    local_point_fn(typename internal::p4est::types<dim>::forest *  forest,
+                   typename internal::p4est::types<dim>::topidx    which_tree,
+                   typename internal::p4est::types<dim>::quadrant *quadrant,
+                   int                                             rank_begin,
+                   int                                             rank_end,
+                   void *                                          point);
+
+  private:
+    /**
+     * Simple struct to keep relavant data. Can be accessed though p4est's user
+     * pointer.
+     */
+    class QuadrantData
+    {
+    public:
+      QuadrantData();
+
+      void
+      set_cell_vertices(
+        typename internal::p4est::types<dim>::forest *  forest,
+        typename internal::p4est::types<dim>::topidx    which_tree,
+        typename internal::p4est::types<dim>::quadrant *quadrant,
+        const typename internal::p4est::types<dim>::quadrant_coord
+          quad_length_on_level);
+
+      void
+      initialize_mapping();
+
+      Point<dim>
+      map_real_to_unit_cell(const Point<dim> &p) const;
+
+      bool
+      is_in_this_quadrant(const Point<dim> &p) const;
+
+    private:
+      std::vector<Point<dim>> cell_vertices;
+
+      /**
+       * Matrix holds coefficients mapping from this physical cell to unit
+       * cell.
+       */
+      FullMatrix<double> quadrant_mapping_matrix;
+
+      bool are_vertices_initialized;
+
+      bool is_reference_mapping_initialized;
+    };
+
+    /**
+     * Quadrant data to be filled upon call of `local_quadrant_fn`.
+     */
+    QuadrantData quadrant_data;
+  }; // class PartitionSearch
+
+
+
+  template <int dim>
+  int
+  PartitionSearch<dim>::local_quadrant_fn(
+    typename internal::p4est::types<dim>::forest *  forest,
+    typename internal::p4est::types<dim>::topidx    which_tree,
+    typename internal::p4est::types<dim>::quadrant *quadrant,
+    int /* rank_begin */,
+    int /* rank_end */,
+    void * /* this is always nullptr */ point)
+  {
+    // point must be nullptr here
+    Assert(point == nullptr, dealii::ExcInternalError());
+
+    // we need the user pointer
+    // note that this is not available since function is static
+    PartitionSearch<dim> *this_object =
+      reinterpret_cast<PartitionSearch<dim> *>(forest->user_pointer);
+
+    // Avoid p4est macros, instead do bitshifts manually with fixed size types
+    const typename internal::p4est::types<dim>::quadrant_coord
+      quad_length_on_level =
+        1 << (static_cast<typename internal::p4est::types<dim>::quadrant_coord>(
+                (dim == 2 ? P4EST_MAXLEVEL : P8EST_MAXLEVEL)) -
+              static_cast<typename internal::p4est::types<dim>::quadrant_coord>(
+                quadrant->level));
+
+    this_object->quadrant_data.set_cell_vertices(forest,
+                                                 which_tree,
+                                                 quadrant,
+                                                 quad_length_on_level);
+
+    // from cell vertices we can initialize the mapping
+    this_object->quadrant_data.initialize_mapping();
+
+    // always return true since we must decide by point
+    return /* true */ 1;
+  }
+
+
+
+  template <int dim>
+  int
+  PartitionSearch<dim>::local_point_fn(
+    typename internal::p4est::types<dim>::forest *forest,
+    typename internal::p4est::types<dim>::topidx /* which_tree */,
+    typename internal::p4est::types<dim>::quadrant * /* quadrant */,
+    int   rank_begin,
+    int   rank_end,
+    void *point)
+  {
+    // point must NOT be be nullptr here
+    Assert(point != nullptr, dealii::ExcInternalError());
+
+    // we need the user pointer
+    // note that this is not available since function is static
+    PartitionSearch<dim> *this_object =
+      reinterpret_cast<PartitionSearch<dim> *>(forest->user_pointer);
+
+    // point with rank as double pointer
+    double *this_point_dptr = static_cast<double *>(point);
+
+    Point<dim> this_point =
+      (dim == 2 ? Point<dim>(this_point_dptr[0], this_point_dptr[1]) :
+                  Point<dim>(this_point_dptr[0],
+                             this_point_dptr[1],
+                             this_point_dptr[2]));
+
+    // use reference mapping to decide whether this point is in this quadrant
+    const bool is_in_this_quadrant =
+      this_object->quadrant_data.is_in_this_quadrant(this_point);
+
+
+
+    if (!is_in_this_quadrant)
+      {
+        // no need to search further, stop recursion
+        return /* false */ 0;
+      }
+
+
+
+    // From here we have a candidate
+    if (rank_begin < rank_end)
+      {
+        // continue recursion
+        return /* true */ 1;
+      }
+
+    // Now, we know that the point is found (rank_begin==rank_end) and we have
+    // the MPI rank, so no need to search further.
+    this_point_dptr[dim] = static_cast<double>(rank_begin);
+
+    // stop recursion.
+    return /* false */ 0;
+  }
+
+
+
+  template <int dim>
+  bool
+  PartitionSearch<dim>::QuadrantData::is_in_this_quadrant(
+    const Point<dim> &p) const
+  {
+    const Point<dim> p_ref = map_real_to_unit_cell(p);
+
+    return GeometryInfo<dim>::is_inside_unit_cell(p_ref);
+  }
+
+
+
+  template <int dim>
+  Point<dim>
+  PartitionSearch<dim>::QuadrantData::map_real_to_unit_cell(
+    const Point<dim> &p) const
+  {
+    Assert(is_reference_mapping_initialized,
+           dealii::ExcMessage(
+             "Cell vertices and mapping coefficients must be fully "
+             "initialized before transforming a point to the unit cell."));
+
+    Point<dim> p_out;
+
+    if (dim == 2)
+      {
+        for (unsigned int alpha = 0;
+             alpha < GeometryInfo<dim>::vertices_per_cell;
+             ++alpha)
+          {
+            const Point<dim> &p_ref =
+              GeometryInfo<dim>::unit_cell_vertex(alpha);
+
+            p_out += (quadrant_mapping_matrix(alpha, 0) +
+                      quadrant_mapping_matrix(alpha, 1) * p(0) +
+                      quadrant_mapping_matrix(alpha, 2) * p(1) +
+                      quadrant_mapping_matrix(alpha, 3) * p(0) * p(1)) *
+                     p_ref;
+          }
+      }
+    else
+      {
+        for (unsigned int alpha = 0;
+             alpha < GeometryInfo<dim>::vertices_per_cell;
+             ++alpha)
+          {
+            const Point<dim> &p_ref =
+              GeometryInfo<dim>::unit_cell_vertex(alpha);
+
+            p_out += (quadrant_mapping_matrix(alpha, 0) +
+                      quadrant_mapping_matrix(alpha, 1) * p(0) +
+                      quadrant_mapping_matrix(alpha, 2) * p(1) +
+                      quadrant_mapping_matrix(alpha, 3) * p(2) +
+                      quadrant_mapping_matrix(alpha, 4) * p(0) * p(1) +
+                      quadrant_mapping_matrix(alpha, 5) * p(1) * p(2) +
+                      quadrant_mapping_matrix(alpha, 6) * p(0) * p(2) +
+                      quadrant_mapping_matrix(alpha, 7) * p(0) * p(1) * p(2)) *
+                     p_ref;
+          }
+      }
+
+    return p_out;
+  }
+
+
+  template <int dim>
+  PartitionSearch<dim>::QuadrantData::QuadrantData()
+    : cell_vertices(GeometryInfo<dim>::vertices_per_cell)
+    , quadrant_mapping_matrix(GeometryInfo<dim>::vertices_per_cell,
+                              GeometryInfo<dim>::vertices_per_cell)
+    , are_vertices_initialized(false)
+    , is_reference_mapping_initialized(false)
+  {}
+
+
+
+  template <int dim>
+  void
+  PartitionSearch<dim>::QuadrantData::initialize_mapping()
+  {
+    Assert(
+      are_vertices_initialized,
+      dealii::ExcMessage(
+        "Cell vertices must be initialized before the cell mapping can be filled."));
+
+    FullMatrix<double> point_matrix(GeometryInfo<dim>::vertices_per_cell,
+                                    GeometryInfo<dim>::vertices_per_cell);
+
+    if (dim == 2)
+      {
+        for (unsigned int alpha = 0;
+             alpha < GeometryInfo<dim>::vertices_per_cell;
+             ++alpha)
+          {
+            // point matrix to be inverted
+            point_matrix(0, alpha) = 1;
+            point_matrix(1, alpha) = cell_vertices[alpha](0);
+            point_matrix(2, alpha) = cell_vertices[alpha](1);
+            point_matrix(3, alpha) =
+              cell_vertices[alpha](0) * cell_vertices[alpha](1);
+          }
+
+        /*
+         * Rows of quadrant_mapping_matrix are the coefficients of the basis
+         * on the physical cell
+         */
+        quadrant_mapping_matrix.invert(point_matrix);
+      }
+    else
+      {
+        for (unsigned int alpha = 0;
+             alpha < GeometryInfo<dim>::vertices_per_cell;
+             ++alpha)
+          {
+            // point matrix to be inverted
+            point_matrix(0, alpha) = 1;
+            point_matrix(1, alpha) = cell_vertices[alpha](0);
+            point_matrix(2, alpha) = cell_vertices[alpha](1);
+            point_matrix(3, alpha) = cell_vertices[alpha](2);
+            point_matrix(4, alpha) =
+              cell_vertices[alpha](0) * cell_vertices[alpha](1);
+            point_matrix(5, alpha) =
+              cell_vertices[alpha](1) * cell_vertices[alpha](2);
+            point_matrix(6, alpha) =
+              cell_vertices[alpha](0) * cell_vertices[alpha](2);
+            point_matrix(7, alpha) = cell_vertices[alpha](0) *
+                                     cell_vertices[alpha](1) *
+                                     cell_vertices[alpha](2);
+          }
+
+        /*
+         * Rows of quadrant_mapping_matrix are the coefficients of the basis
+         * on the physical cell
+         */
+        quadrant_mapping_matrix.invert(point_matrix);
+      }
+
+    is_reference_mapping_initialized = true;
+  }
+
+
+
+  template <>
+  void
+  PartitionSearch<2>::QuadrantData::set_cell_vertices(
+    typename internal::p4est::types<2>::forest *  forest,
+    typename internal::p4est::types<2>::topidx    which_tree,
+    typename internal::p4est::types<2>::quadrant *quadrant,
+    const typename internal::p4est::types<2>::quadrant_coord
+      quad_length_on_level)
+  {
+    constexpr int dim = 2;
+
+    // p4est for some reason always needs double vxyz[3] as last argument to
+    // quadrant_coord_to_vertex
+    double corner_point[dim + 1] = {0};
+
+    // A lambda to avoid code duplication.
+    const auto copy_vertex = [&](unsigned int vertex_index) -> void {
+      // copy into local struct
+      for (size_t d = 0; d < dim; ++d)
+        {
+          cell_vertices[vertex_index](d) = corner_point[d];
+          // reset
+          corner_point[d] = 0;
+        }
+    };
+
+    // Fill points of QuadrantData in lexicographic order
+    /*
+     * Corner #0
+     */
+    unsigned int vertex_index = 0;
+    internal::p4est::functions<dim>::quadrant_coord_to_vertex(
+      forest->connectivity, which_tree, quadrant->x, quadrant->y, corner_point);
+
+    // copy into local struct
+    copy_vertex(vertex_index);
+
+    /*
+     * Corner #1
+     */
+    vertex_index = 1;
+    internal::p4est::functions<dim>::quadrant_coord_to_vertex(
+      forest->connectivity,
+      which_tree,
+      quadrant->x + quad_length_on_level,
+      quadrant->y,
+      corner_point);
+
+    // copy into local struct
+    copy_vertex(vertex_index);
+
+    /*
+     * Corner #2
+     */
+    vertex_index = 2;
+    internal::p4est::functions<dim>::quadrant_coord_to_vertex(
+      forest->connectivity,
+      which_tree,
+      quadrant->x,
+      quadrant->y + quad_length_on_level,
+      corner_point);
+
+    // copy into local struct
+    copy_vertex(vertex_index);
+
+    /*
+     * Corner #3
+     */
+    vertex_index = 3;
+    internal::p4est::functions<dim>::quadrant_coord_to_vertex(
+      forest->connectivity,
+      which_tree,
+      quadrant->x + quad_length_on_level,
+      quadrant->y + quad_length_on_level,
+      corner_point);
+
+    // copy into local struct
+    copy_vertex(vertex_index);
+
+    are_vertices_initialized = true;
+  }
+
+
+
+  template <>
+  void
+  PartitionSearch<3>::QuadrantData::set_cell_vertices(
+    typename internal::p4est::types<3>::forest *  forest,
+    typename internal::p4est::types<3>::topidx    which_tree,
+    typename internal::p4est::types<3>::quadrant *quadrant,
+    const typename internal::p4est::types<3>::quadrant_coord
+      quad_length_on_level)
+  {
+    constexpr int dim = 3;
+
+    double corner_point[dim] = {0};
+
+    // A lambda to avoid code duplication.
+    auto copy_vertex = [&](unsigned int vertex_index) -> void {
+      // copy into local struct
+      for (size_t d = 0; d < dim; ++d)
+        {
+          cell_vertices[vertex_index](d) = corner_point[d];
+          // reset
+          corner_point[d] = 0;
+        }
+    };
+
+    // Fill points of QuadrantData in lexicographic order
+    /*
+     * Corner #0
+     */
+    unsigned int vertex_index = 0;
+    internal::p4est::functions<dim>::quadrant_coord_to_vertex(
+      forest->connectivity,
+      which_tree,
+      quadrant->x,
+      quadrant->y,
+      quadrant->z,
+      corner_point);
+
+    // copy into local struct
+    copy_vertex(vertex_index);
+
+
+    /*
+     * Corner #1
+     */
+    vertex_index = 1;
+    internal::p4est::functions<dim>::quadrant_coord_to_vertex(
+      forest->connectivity,
+      which_tree,
+      quadrant->x + quad_length_on_level,
+      quadrant->y,
+      quadrant->z,
+      corner_point);
+
+    // copy into local struct
+    copy_vertex(vertex_index);
+
+    /*
+     * Corner #2
+     */
+    vertex_index = 2;
+    internal::p4est::functions<dim>::quadrant_coord_to_vertex(
+      forest->connectivity,
+      which_tree,
+      quadrant->x,
+      quadrant->y + quad_length_on_level,
+      quadrant->z,
+      corner_point);
+
+    // copy into local struct
+    copy_vertex(vertex_index);
+
+    /*
+     * Corner #3
+     */
+    vertex_index = 3;
+    internal::p4est::functions<dim>::quadrant_coord_to_vertex(
+      forest->connectivity,
+      which_tree,
+      quadrant->x + quad_length_on_level,
+      quadrant->y + quad_length_on_level,
+      quadrant->z,
+      corner_point);
+
+    // copy into local struct
+    copy_vertex(vertex_index);
+
+    /*
+     * Corner #4
+     */
+    vertex_index = 4;
+    internal::p4est::functions<dim>::quadrant_coord_to_vertex(
+      forest->connectivity,
+      which_tree,
+      quadrant->x,
+      quadrant->y,
+      quadrant->z + quad_length_on_level,
+      corner_point);
+
+    // copy into local struct
+    copy_vertex(vertex_index);
+
+    /*
+     * Corner #5
+     */
+    vertex_index = 5;
+    internal::p4est::functions<dim>::quadrant_coord_to_vertex(
+      forest->connectivity,
+      which_tree,
+      quadrant->x + quad_length_on_level,
+      quadrant->y,
+      quadrant->z + quad_length_on_level,
+      corner_point);
+
+    // copy into local struct
+    copy_vertex(vertex_index);
+
+    /*
+     * Corner #6
+     */
+    vertex_index = 6;
+    internal::p4est::functions<dim>::quadrant_coord_to_vertex(
+      forest->connectivity,
+      which_tree,
+      quadrant->x,
+      quadrant->y + quad_length_on_level,
+      quadrant->z + quad_length_on_level,
+      corner_point);
+
+    // copy into local struct
+    copy_vertex(vertex_index);
+
+    /*
+     * Corner #7
+     */
+    vertex_index = 7;
+    internal::p4est::functions<dim>::quadrant_coord_to_vertex(
+      forest->connectivity,
+      which_tree,
+      quadrant->x + quad_length_on_level,
+      quadrant->y + quad_length_on_level,
+      quadrant->z + quad_length_on_level,
+      corner_point);
+
+    // copy into local struct
+    copy_vertex(vertex_index);
+
+
+    are_vertices_initialized = true;
+  }
+#  endif //  P4EST_SEARCH_LOCAL defined
+
 
   /**
    * A data structure that we use to store which cells (indicated by
@@ -1127,7 +1705,14 @@ namespace parallel
             Triangulation<dim, spacedim>::limit_level_difference_at_vertices) :
           smooth_grid,
         false)
+#  ifdef DEBUG
+      // Always communicate vertices to p4est in debug mode
+      , settings(
+          static_cast<Settings>(settings | communicate_vertices_to_p4est))
+#  else
       , settings(settings)
+#  endif
+
       , triangulation_has_content(false)
       , connectivity(nullptr)
       , parallel_forest(nullptr)
@@ -1269,6 +1854,16 @@ namespace parallel
 
 
 
+    template <int dim, int spacedim>
+    bool
+    Triangulation<dim, spacedim>::are_vertices_communicated_to_p4est() const
+    {
+      return settings &
+             Triangulation<dim, spacedim>::communicate_vertices_to_p4est;
+    }
+
+
+
     template <int dim, int spacedim>
     void
     Triangulation<dim, spacedim>::execute_transfer(
@@ -1414,6 +2009,12 @@ namespace parallel
     {
       Assert(parallel_forest != nullptr,
              ExcMessage("Can't produce output when no forest is created yet."));
+
+      AssertThrow(are_vertices_communicated_to_p4est(),
+                  ExcMessage(
+                    "To use this function the triangulation's flag "
+                    "Settings::communicate_vertices_to_p4est must be set."));
+
       dealii::internal::p4est::functions<dim>::vtk_write_file(
         parallel_forest, nullptr, file_basename.c_str());
     }
@@ -1708,13 +2309,7 @@ namespace parallel
       // now create a connectivity object with the right sizes for all
       // arrays. set vertex information only in debug mode (saves a few bytes
       // in optimized mode)
-      const bool set_vertex_info
-#    ifdef DEBUG
-        = true
-#    else
-        = false
-#    endif
-        ;
+      const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
 
       connectivity = dealii::internal::p4est::functions<2>::connectivity_new(
         (set_vertex_info == true ? this->n_vertices() : 0),
@@ -1774,13 +2369,7 @@ namespace parallel
       // now create a connectivity object with the right sizes for all
       // arrays. set vertex information only in debug mode (saves a few bytes
       // in optimized mode)
-      const bool set_vertex_info
-#    ifdef DEBUG
-        = true
-#    else
-        = false
-#    endif
-        ;
+      const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
 
       connectivity = dealii::internal::p4est::functions<2>::connectivity_new(
         (set_vertex_info == true ? this->n_vertices() : 0),
@@ -1843,13 +2432,7 @@ namespace parallel
         std::accumulate(edge_touch_count.begin(), edge_touch_count.end(), 0u);
 
       // now create a connectivity object with the right sizes for all arrays
-      const bool set_vertex_info
-#    ifdef DEBUG
-        = true
-#    else
-        = false
-#    endif
-        ;
+      const bool set_vertex_info = this->are_vertices_communicated_to_p4est();
 
       connectivity = dealii::internal::p4est::functions<3>::connectivity_new(
         (set_vertex_info == true ? this->n_vertices() : 0),
@@ -1974,20 +2557,22 @@ namespace parallel
         for (unsigned int i = 0; i < topological_vertex_numbering.size(); ++i)
           topological_vertex_numbering[i] = i;
         // combine vertices that have different locations (and thus, different
-        // vertex_index) but represent the same topological entity over periodic
-        // boundaries. The vector topological_vertex_numbering contains a linear
-        // map from 0 to n_vertices at input and at output relates periodic
-        // vertices with only one vertex index. The output is used to always
-        // identify the same vertex according to the periodicity, e.g. when
-        // finding the maximum cell level around a vertex.
+        // vertex_index) but represent the same topological entity over
+        // periodic boundaries. The vector topological_vertex_numbering
+        // contains a linear map from 0 to n_vertices at input and at output
+        // relates periodic vertices with only one vertex index. The output is
+        // used to always identify the same vertex according to the
+        // periodicity, e.g. when finding the maximum cell level around a
+        // vertex.
         //
-        // Example: On a 3D cell with vertices numbered from 0 to 7 and periodic
-        // boundary conditions in x direction, the vector
+        // Example: On a 3D cell with vertices numbered from 0 to 7 and
+        // periodic boundary conditions in x direction, the vector
         // topological_vertex_numbering will contain the numbers
         // {0,0,2,2,4,4,6,6} (because the vertex pairs {0,1}, {2,3}, {4,5},
-        // {6,7} belong together, respectively). If periodicity is set in x and
-        // z direction, the output is {0,0,2,2,0,0,2,2}, and if periodicity is
-        // in all directions, the output is simply {0,0,0,0,0,0,0,0}.
+        // {6,7} belong together, respectively). If periodicity is set in x
+        // and z direction, the output is {0,0,2,2,0,0,2,2}, and if
+        // periodicity is in all directions, the output is simply
+        // {0,0,0,0,0,0,0,0}.
         using cell_iterator =
           typename Triangulation<dim, spacedim>::cell_iterator;
         typename std::map<std::pair<cell_iterator, unsigned int>,
@@ -2009,8 +2594,8 @@ namespace parallel
                      v < GeometryInfo<dim - 1>::vertices_per_cell;
                      ++v)
                   {
-                    // take possible non-standard orientation of face on cell[0]
-                    // into account
+                    // take possible non-standard orientation of face on
+                    // cell[0] into account
                     const unsigned int vface0 =
                       GeometryInfo<dim>::standard_to_real_face_vertex(
                         v,
@@ -2191,8 +2776,9 @@ namespace parallel
           mesh_changed = enforce_mesh_balance_over_periodic_boundaries(*this);
 
           // We can't be sure that we won't run into a situation where we can
-          // not reconcile mesh smoothing and balancing of periodic faces. As we
-          // don't know what else to do, at least abort with an error message.
+          // not reconcile mesh smoothing and balancing of periodic faces. As
+          // we don't know what else to do, at least abort with an error
+          // message.
           ++loop_counter;
           AssertThrow(
             loop_counter < 32,
@@ -2417,8 +3003,9 @@ namespace parallel
 
       // fill level_subdomain_ids for geometric multigrid
       // the level ownership of a cell is defined as the owner if the cell is
-      // active or as the owner of child(0) we need this information for all our
-      // ancestors and the same-level neighbors of our own cells (=level ghosts)
+      // active or as the owner of child(0) we need this information for all
+      // our ancestors and the same-level neighbors of our own cells (=level
+      // ghosts)
       if (settings & construct_multigrid_hierarchy)
         {
           // step 1: We set our own ids all the way down and all the others to
@@ -2447,8 +3034,8 @@ namespace parallel
                 }
             }
 
-          // step 2: make sure all the neighbors to our level_cells exist. Need
-          // to look up in p4est...
+          // step 2: make sure all the neighbors to our level_cells exist.
+          // Need to look up in p4est...
           std::vector<std::vector<bool>> marked_vertices(this->n_levels());
           for (unsigned int lvl = 0; lvl < this->n_levels(); ++lvl)
             marked_vertices[lvl] = mark_locally_active_vertices_on_level(lvl);
@@ -2552,6 +3139,119 @@ namespace parallel
 
 
 
+    template <int dim, int spacedim>
+    types::subdomain_id
+    Triangulation<dim, spacedim>::find_point_owner_rank(const Point<dim> &p)
+    {
+      // Call the other function
+      std::vector<Point<dim>>          point{p};
+      std::vector<types::subdomain_id> owner = find_point_owner_rank(point);
+
+      return owner[0];
+    }
+
+
+
+    template <int dim, int spacedim>
+    std::vector<types::subdomain_id>
+    Triangulation<dim, spacedim>::find_point_owner_rank(
+      const std::vector<Point<dim>> &points)
+    {
+#  ifndef P4EST_SEARCH_LOCAL
+      (void)points;
+      AssertThrow(
+        false,
+        ExcMessage(
+          "This function is only available if p4est is version 2.2 and higher."));
+      // Just return to satisfy compiler
+      return std::vector<unsigned int>(1,
+                                       dealii::numbers::invalid_subdomain_id);
+#  else
+      // We can only use this function if vertices are communicated to p4est
+      AssertThrow(this->are_vertices_communicated_to_p4est(),
+                  ExcMessage(
+                    "Vertices need to be communicated to p4est to use this "
+                    "function. This must explicitly be turned on in the "
+                    "settings of the triangulation's constructor."));
+
+      // We can only use this function if all manifolds are flat
+      for (const auto &manifold_id : this->get_manifold_ids())
+        {
+          AssertThrow(
+            manifold_id == numbers::flat_manifold_id,
+            ExcMessage(
+              "This function can only be used if the triangulation "
+              "has no other manifold than a Cartesian (flat) manifold attached."));
+        }
+
+      // Create object for callback
+      PartitionSearch<dim> partition_search;
+
+      // Pointer should be this triangulation before we set it to something else
+      Assert(parallel_forest->user_pointer == this, ExcInternalError());
+
+      // re-assign p4est's user pointer
+      parallel_forest->user_pointer = &partition_search;
+
+      //
+      // Copy points into p4est internal array data struct
+      //
+      // pointer to an array of points.
+      sc_array_t *point_sc_array;
+      // allocate memory for a number of dim-dimensional points including their
+      // MPI rank, i.e., dim + 1 fields
+      point_sc_array =
+        sc_array_new_count(sizeof(double[dim + 1]), points.size());
+
+      // now assign the actual value
+      for (size_t i = 0; i < points.size(); ++i)
+        {
+          // alias
+          const Point<dim> &p = points[i];
+          // get a non-const view of the array
+          double *this_sc_point =
+            static_cast<double *>(sc_array_index_ssize_t(point_sc_array, i));
+          // fill this with the point data
+          for (unsigned int d = 0; d < dim; ++d)
+            {
+              this_sc_point[d] = p(d);
+            }
+          this_sc_point[dim] = -1.0; // owner rank
+        }
+
+      dealii::internal::p4est::functions<dim>::search_partition(
+        parallel_forest,
+        /* execute quadrant function when leaving quadrant */
+        static_cast<int>(false),
+        &PartitionSearch<dim>::local_quadrant_fn,
+        &PartitionSearch<dim>::local_point_fn,
+        point_sc_array);
+
+      // copy the points found to an std::array
+      std::vector<types::subdomain_id> owner_rank(
+        points.size(), numbers::invalid_subdomain_id);
+
+      // fill the array
+      for (size_t i = 0; i < points.size(); ++i)
+        {
+          // get a non-const view of the array
+          double *this_sc_point =
+            static_cast<double *>(sc_array_index_ssize_t(point_sc_array, i));
+          owner_rank[i] = static_cast<types::subdomain_id>(this_sc_point[dim]);
+        }
+
+      // reset the internal pointer to this triangulation
+      parallel_forest->user_pointer = this;
+
+      // release the memory (otherwise p4est will complain)
+      sc_array_destroy_null(&point_sc_array);
+
+      return owner_rank;
+#  endif // P4EST_SEARCH_LOCAL defined
+    }
+
+
+
     template <int dim, int spacedim>
     void
     Triangulation<dim, spacedim>::execute_coarsening_and_refinement()
@@ -2655,9 +3355,8 @@ namespace parallel
       // been updated
       this->signals.post_p4est_refinement();
 
-      // before repartitioning the mesh, save a copy of the current positions of
-      // quadrants
-      // only if data needs to be transferred later
+      // before repartitioning the mesh, save a copy of the current positions
+      // of quadrants only if data needs to be transferred later
       std::vector<typename dealii::internal::p4est::types<dim>::gloidx>
         previous_global_first_quadrant;
 
@@ -2673,8 +3372,8 @@ namespace parallel
 
       if (!(settings & no_automatic_repartitioning))
         {
-          // partition the new mesh between all processors. If cell weights have
-          // not been given balance the number of cells.
+          // partition the new mesh between all processors. If cell weights
+          // have not been given balance the number of cells.
           if (this->signals.cell_weight.num_slots() == 0)
             dealii::internal::p4est::functions<dim>::partition(
               parallel_forest,
@@ -2824,9 +3523,8 @@ namespace parallel
       // signal that repartitioning is going to happen
       this->signals.pre_distributed_repartition();
 
-      // before repartitioning the mesh, save a copy of the current positions of
-      // quadrants
-      // only if data needs to be transferred later
+      // before repartitioning the mesh, save a copy of the current positions
+      // of quadrants only if data needs to be transferred later
       std::vector<typename dealii::internal::p4est::types<dim>::gloidx>
         previous_global_first_quadrant;
 
@@ -2856,8 +3554,8 @@ namespace parallel
 
           PartitionWeights<dim, spacedim> partition_weights(cell_weights);
 
-          // attach (temporarily) a pointer to the cell weights through p4est's
-          // user_pointer object
+          // attach (temporarily) a pointer to the cell weights through
+          // p4est's user_pointer object
           Assert(parallel_forest->user_pointer == this, ExcInternalError());
           parallel_forest->user_pointer = &partition_weights;
 
@@ -3173,7 +3871,8 @@ namespace parallel
           Assert(false, ExcInternalError());
         }
 
-      // The range of ghost_owners might have changed so update that information
+      // The range of ghost_owners might have changed so update that
+      // information
       this->update_number_cache();
     }
 
@@ -3531,6 +4230,16 @@ namespace parallel
 
 
 
+    template <int spacedim>
+    bool
+    Triangulation<1, spacedim>::are_vertices_communicated_to_p4est() const
+    {
+      Assert(false, ExcNotImplemented());
+      return false;
+    }
+
+
+
     template <int spacedim>
     void
     Triangulation<1, spacedim>::update_cell_relations()
index a72d7ca901b92eaad722df527ad2bcd33dc64726..3deab18dd85ed9a8941e00a3b0136b6ffade7363 100644 (file)
@@ -35,6 +35,11 @@ void
 write_vtk(const parallel::distributed::Triangulation<dim, spacedim> &tria,
           const char *                                               filename)
 {
+  AssertThrow(tria.are_vertices_communicated_to_p4est(),
+              ExcMessage("To use this function the flag "
+                         "Settings::communicate_vertices_to_p4est "
+                         "must be set in the triangulation."));
+
   deallog << "Checksum: " << tria.get_checksum() << std::endl;
 
   tria.write_mesh_vtk(filename);
index 65b8e95962743bbeb611a99246146e470c7aceae..ac15dffe55103a710bd0bb443aad85681add31e0 100644 (file)
@@ -38,6 +38,11 @@ void
 write_vtk(const parallel::distributed::Triangulation<dim, spacedim> &tria,
           const char *                                               filename)
 {
+  AssertThrow(tria.are_vertices_communicated_to_p4est(),
+              ExcMessage("To use this function the flag "
+                         "Settings::communicate_vertices_to_p4est "
+                         "must be set in the triangulation."));
+
   deallog << "Checksum: " << tria.get_checksum() << std::endl;
 
   tria.write_mesh_vtk(filename);
diff --git a/tests/mpi/p4est_2d_find_point_owner_rank_01.cc b/tests/mpi/p4est_2d_find_point_owner_rank_01.cc
new file mode 100644 (file)
index 0000000..4da881e
--- /dev/null
@@ -0,0 +1,42 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/utilities.h>
+
+#include "../tests.h"
+
+// contains the actual templated tests
+#include "p4est_find_point_owner_rank.h"
+
+int
+main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  MPILogInitAll                    log;
+
+
+  /*
+   * We craete a distrubuted mesh with three cells (tree roots), e.g., a
+   * hyper_L. We want to find the mpi rank of a single fixed point. On all
+   * processes we must find the same owner.
+   */
+  const std::vector<Point<2>> point(1, Point<2>(0.23, 0.77));
+
+  deallog.push("2D-Test");
+  test<2>(point);
+  deallog.pop();
+}
diff --git a/tests/mpi/p4est_2d_find_point_owner_rank_01.mpirun=1.with_p4est=true.with_p4est_with_search_local=true.output b/tests/mpi/p4est_2d_find_point_owner_rank_01.mpirun=1.with_p4est=true.with_p4est_with_search_local=true.output
new file mode 100644 (file)
index 0000000..13b37b6
--- /dev/null
@@ -0,0 +1,13 @@
+
+DEAL:0:2D-Test:: Point search on subdivided hyper cube...
+DEAL:0:2D-Test::   Number of MPI processes = 1
+DEAL:0:2D-Test::   Number of this MPI processes = 0
+DEAL:0:2D-Test::   Number of cells = 256
+DEAL:0:2D-Test::   Number of levels = 4
+DEAL:0:2D-Test::   Number of global levels = 4
+DEAL:0:2D-Test::   Triangulation checksum = 1754413825
+DEAL:0:2D-Test::   points = 0.230000 0.770000    
+DEAL:0:2D-Test::
+DEAL:0:2D-Test::   All point owner ranks found on all processes equal each other. TEST PASSED.
+DEAL:0:2D-Test::
+DEAL:0:2D-Test::   >>> Reached end of test <<<
diff --git a/tests/mpi/p4est_2d_find_point_owner_rank_01.mpirun=3.with_p4est=true.with_p4est_with_search_local=true.output b/tests/mpi/p4est_2d_find_point_owner_rank_01.mpirun=3.with_p4est=true.with_p4est_with_search_local=true.output
new file mode 100644 (file)
index 0000000..44f9155
--- /dev/null
@@ -0,0 +1,35 @@
+
+DEAL:0:2D-Test:: Point search on subdivided hyper cube...
+DEAL:0:2D-Test::   Number of MPI processes = 3
+DEAL:0:2D-Test::   Number of this MPI processes = 0
+DEAL:0:2D-Test::   Number of cells = 256
+DEAL:0:2D-Test::   Number of levels = 4
+DEAL:0:2D-Test::   Number of global levels = 4
+DEAL:0:2D-Test::   Triangulation checksum = 1754413825
+DEAL:0:2D-Test::   points = 0.230000 0.770000    
+DEAL:0:2D-Test::
+DEAL:0:2D-Test::   All point owner ranks found on all processes equal each other. TEST PASSED.
+DEAL:0:2D-Test::
+DEAL:0:2D-Test::   >>> Reached end of test <<<
+
+DEAL:1:2D-Test:: Point search on subdivided hyper cube...
+DEAL:1:2D-Test::   Number of MPI processes = 3
+DEAL:1:2D-Test::   Number of this MPI processes = 1
+DEAL:1:2D-Test::   Number of cells = 256
+DEAL:1:2D-Test::   Number of levels = 4
+DEAL:1:2D-Test::   Number of global levels = 4
+DEAL:1:2D-Test::   Triangulation checksum = 0
+DEAL:1:2D-Test::   points = 0.230000 0.770000    
+DEAL:1:2D-Test::   >>> Reached end of test <<<
+
+
+DEAL:2:2D-Test:: Point search on subdivided hyper cube...
+DEAL:2:2D-Test::   Number of MPI processes = 3
+DEAL:2:2D-Test::   Number of this MPI processes = 2
+DEAL:2:2D-Test::   Number of cells = 256
+DEAL:2:2D-Test::   Number of levels = 4
+DEAL:2:2D-Test::   Number of global levels = 4
+DEAL:2:2D-Test::   Triangulation checksum = 0
+DEAL:2:2D-Test::   points = 0.230000 0.770000    
+DEAL:2:2D-Test::   >>> Reached end of test <<<
+
diff --git a/tests/mpi/p4est_2d_find_point_owner_rank_01.mpirun=5.with_p4est=true.with_p4est_with_search_local=true.output b/tests/mpi/p4est_2d_find_point_owner_rank_01.mpirun=5.with_p4est=true.with_p4est_with_search_local=true.output
new file mode 100644 (file)
index 0000000..3373add
--- /dev/null
@@ -0,0 +1,57 @@
+
+DEAL:0:2D-Test:: Point search on subdivided hyper cube...
+DEAL:0:2D-Test::   Number of MPI processes = 5
+DEAL:0:2D-Test::   Number of this MPI processes = 0
+DEAL:0:2D-Test::   Number of cells = 256
+DEAL:0:2D-Test::   Number of levels = 4
+DEAL:0:2D-Test::   Number of global levels = 4
+DEAL:0:2D-Test::   Triangulation checksum = 1754413825
+DEAL:0:2D-Test::   points = 0.230000 0.770000    
+DEAL:0:2D-Test::
+DEAL:0:2D-Test::   All point owner ranks found on all processes equal each other. TEST PASSED.
+DEAL:0:2D-Test::
+DEAL:0:2D-Test::   >>> Reached end of test <<<
+
+DEAL:1:2D-Test:: Point search on subdivided hyper cube...
+DEAL:1:2D-Test::   Number of MPI processes = 5
+DEAL:1:2D-Test::   Number of this MPI processes = 1
+DEAL:1:2D-Test::   Number of cells = 256
+DEAL:1:2D-Test::   Number of levels = 4
+DEAL:1:2D-Test::   Number of global levels = 4
+DEAL:1:2D-Test::   Triangulation checksum = 0
+DEAL:1:2D-Test::   points = 0.230000 0.770000    
+DEAL:1:2D-Test::   >>> Reached end of test <<<
+
+
+DEAL:2:2D-Test:: Point search on subdivided hyper cube...
+DEAL:2:2D-Test::   Number of MPI processes = 5
+DEAL:2:2D-Test::   Number of this MPI processes = 2
+DEAL:2:2D-Test::   Number of cells = 256
+DEAL:2:2D-Test::   Number of levels = 4
+DEAL:2:2D-Test::   Number of global levels = 4
+DEAL:2:2D-Test::   Triangulation checksum = 0
+DEAL:2:2D-Test::   points = 0.230000 0.770000    
+DEAL:2:2D-Test::   >>> Reached end of test <<<
+
+
+DEAL:3:2D-Test:: Point search on subdivided hyper cube...
+DEAL:3:2D-Test::   Number of MPI processes = 5
+DEAL:3:2D-Test::   Number of this MPI processes = 3
+DEAL:3:2D-Test::   Number of cells = 256
+DEAL:3:2D-Test::   Number of levels = 4
+DEAL:3:2D-Test::   Number of global levels = 4
+DEAL:3:2D-Test::   Triangulation checksum = 0
+DEAL:3:2D-Test::   points = 0.230000 0.770000    
+DEAL:3:2D-Test::   >>> Reached end of test <<<
+
+
+DEAL:4:2D-Test:: Point search on subdivided hyper cube...
+DEAL:4:2D-Test::   Number of MPI processes = 5
+DEAL:4:2D-Test::   Number of this MPI processes = 4
+DEAL:4:2D-Test::   Number of cells = 256
+DEAL:4:2D-Test::   Number of levels = 4
+DEAL:4:2D-Test::   Number of global levels = 4
+DEAL:4:2D-Test::   Triangulation checksum = 0
+DEAL:4:2D-Test::   points = 0.230000 0.770000    
+DEAL:4:2D-Test::   >>> Reached end of test <<<
+
diff --git a/tests/mpi/p4est_2d_find_point_owner_rank_02.cc b/tests/mpi/p4est_2d_find_point_owner_rank_02.cc
new file mode 100644 (file)
index 0000000..19b82ee
--- /dev/null
@@ -0,0 +1,50 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/utilities.h>
+
+#include "../tests.h"
+
+// contains the actual templated tests
+#include "p4est_find_point_owner_rank.h"
+
+
+int
+main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  MPILogInitAll                    log;
+
+
+  /*
+   * We craete a distrubuted mesh with three cells (tree roots), e.g., a
+   * hyper_L. We want to find the mpi rank of a set of fixed points. On all
+   * processes we must find the same owner ranks.
+   */
+  std::vector<Point<2>> points;
+  points.emplace_back(Point<2>(0.23, 0.23));
+  points.emplace_back(Point<2>(0.23, 0.77));
+  points.emplace_back(Point<2>(0.77, 0.23));
+  points.emplace_back(Point<2>(0.77, 0.77));
+
+  // this point is outside the mesh
+  points.emplace_back(Point<2>(1.1, 1.1));
+
+  deallog.push("2D-Test");
+  test<2>(points);
+  deallog.pop();
+}
diff --git a/tests/mpi/p4est_2d_find_point_owner_rank_02.mpirun=1.with_p4est=true.with_p4est_with_search_local=true.output b/tests/mpi/p4est_2d_find_point_owner_rank_02.mpirun=1.with_p4est=true.with_p4est_with_search_local=true.output
new file mode 100644 (file)
index 0000000..6c9cd9e
--- /dev/null
@@ -0,0 +1,13 @@
+
+DEAL:0:2D-Test:: Point search on subdivided hyper cube...
+DEAL:0:2D-Test::   Number of MPI processes = 1
+DEAL:0:2D-Test::   Number of this MPI processes = 0
+DEAL:0:2D-Test::   Number of cells = 256
+DEAL:0:2D-Test::   Number of levels = 4
+DEAL:0:2D-Test::   Number of global levels = 4
+DEAL:0:2D-Test::   Triangulation checksum = 1754413825
+DEAL:0:2D-Test::   points = 0.230000 0.230000    0.230000 0.770000    0.770000 0.230000    0.770000 0.770000    1.10000 1.10000    
+DEAL:0:2D-Test::
+DEAL:0:2D-Test::   All point owner ranks found on all processes equal each other. TEST PASSED.
+DEAL:0:2D-Test::
+DEAL:0:2D-Test::   >>> Reached end of test <<<
diff --git a/tests/mpi/p4est_2d_find_point_owner_rank_02.mpirun=3.with_p4est=true.with_p4est_with_search_local=true.output b/tests/mpi/p4est_2d_find_point_owner_rank_02.mpirun=3.with_p4est=true.with_p4est_with_search_local=true.output
new file mode 100644 (file)
index 0000000..718bd08
--- /dev/null
@@ -0,0 +1,35 @@
+
+DEAL:0:2D-Test:: Point search on subdivided hyper cube...
+DEAL:0:2D-Test::   Number of MPI processes = 3
+DEAL:0:2D-Test::   Number of this MPI processes = 0
+DEAL:0:2D-Test::   Number of cells = 256
+DEAL:0:2D-Test::   Number of levels = 4
+DEAL:0:2D-Test::   Number of global levels = 4
+DEAL:0:2D-Test::   Triangulation checksum = 1754413825
+DEAL:0:2D-Test::   points = 0.230000 0.230000    0.230000 0.770000    0.770000 0.230000    0.770000 0.770000    1.10000 1.10000    
+DEAL:0:2D-Test::
+DEAL:0:2D-Test::   All point owner ranks found on all processes equal each other. TEST PASSED.
+DEAL:0:2D-Test::
+DEAL:0:2D-Test::   >>> Reached end of test <<<
+
+DEAL:1:2D-Test:: Point search on subdivided hyper cube...
+DEAL:1:2D-Test::   Number of MPI processes = 3
+DEAL:1:2D-Test::   Number of this MPI processes = 1
+DEAL:1:2D-Test::   Number of cells = 256
+DEAL:1:2D-Test::   Number of levels = 4
+DEAL:1:2D-Test::   Number of global levels = 4
+DEAL:1:2D-Test::   Triangulation checksum = 0
+DEAL:1:2D-Test::   points = 0.230000 0.230000    0.230000 0.770000    0.770000 0.230000    0.770000 0.770000    1.10000 1.10000    
+DEAL:1:2D-Test::   >>> Reached end of test <<<
+
+
+DEAL:2:2D-Test:: Point search on subdivided hyper cube...
+DEAL:2:2D-Test::   Number of MPI processes = 3
+DEAL:2:2D-Test::   Number of this MPI processes = 2
+DEAL:2:2D-Test::   Number of cells = 256
+DEAL:2:2D-Test::   Number of levels = 4
+DEAL:2:2D-Test::   Number of global levels = 4
+DEAL:2:2D-Test::   Triangulation checksum = 0
+DEAL:2:2D-Test::   points = 0.230000 0.230000    0.230000 0.770000    0.770000 0.230000    0.770000 0.770000    1.10000 1.10000    
+DEAL:2:2D-Test::   >>> Reached end of test <<<
+
diff --git a/tests/mpi/p4est_2d_find_point_owner_rank_02.mpirun=5.with_p4est=true.with_p4est_with_search_local=true.output b/tests/mpi/p4est_2d_find_point_owner_rank_02.mpirun=5.with_p4est=true.with_p4est_with_search_local=true.output
new file mode 100644 (file)
index 0000000..ed2fd7c
--- /dev/null
@@ -0,0 +1,57 @@
+
+DEAL:0:2D-Test:: Point search on subdivided hyper cube...
+DEAL:0:2D-Test::   Number of MPI processes = 5
+DEAL:0:2D-Test::   Number of this MPI processes = 0
+DEAL:0:2D-Test::   Number of cells = 256
+DEAL:0:2D-Test::   Number of levels = 4
+DEAL:0:2D-Test::   Number of global levels = 4
+DEAL:0:2D-Test::   Triangulation checksum = 1754413825
+DEAL:0:2D-Test::   points = 0.230000 0.230000    0.230000 0.770000    0.770000 0.230000    0.770000 0.770000    1.10000 1.10000    
+DEAL:0:2D-Test::
+DEAL:0:2D-Test::   All point owner ranks found on all processes equal each other. TEST PASSED.
+DEAL:0:2D-Test::
+DEAL:0:2D-Test::   >>> Reached end of test <<<
+
+DEAL:1:2D-Test:: Point search on subdivided hyper cube...
+DEAL:1:2D-Test::   Number of MPI processes = 5
+DEAL:1:2D-Test::   Number of this MPI processes = 1
+DEAL:1:2D-Test::   Number of cells = 256
+DEAL:1:2D-Test::   Number of levels = 4
+DEAL:1:2D-Test::   Number of global levels = 4
+DEAL:1:2D-Test::   Triangulation checksum = 0
+DEAL:1:2D-Test::   points = 0.230000 0.230000    0.230000 0.770000    0.770000 0.230000    0.770000 0.770000    1.10000 1.10000    
+DEAL:1:2D-Test::   >>> Reached end of test <<<
+
+
+DEAL:2:2D-Test:: Point search on subdivided hyper cube...
+DEAL:2:2D-Test::   Number of MPI processes = 5
+DEAL:2:2D-Test::   Number of this MPI processes = 2
+DEAL:2:2D-Test::   Number of cells = 256
+DEAL:2:2D-Test::   Number of levels = 4
+DEAL:2:2D-Test::   Number of global levels = 4
+DEAL:2:2D-Test::   Triangulation checksum = 0
+DEAL:2:2D-Test::   points = 0.230000 0.230000    0.230000 0.770000    0.770000 0.230000    0.770000 0.770000    1.10000 1.10000    
+DEAL:2:2D-Test::   >>> Reached end of test <<<
+
+
+DEAL:3:2D-Test:: Point search on subdivided hyper cube...
+DEAL:3:2D-Test::   Number of MPI processes = 5
+DEAL:3:2D-Test::   Number of this MPI processes = 3
+DEAL:3:2D-Test::   Number of cells = 256
+DEAL:3:2D-Test::   Number of levels = 4
+DEAL:3:2D-Test::   Number of global levels = 4
+DEAL:3:2D-Test::   Triangulation checksum = 0
+DEAL:3:2D-Test::   points = 0.230000 0.230000    0.230000 0.770000    0.770000 0.230000    0.770000 0.770000    1.10000 1.10000    
+DEAL:3:2D-Test::   >>> Reached end of test <<<
+
+
+DEAL:4:2D-Test:: Point search on subdivided hyper cube...
+DEAL:4:2D-Test::   Number of MPI processes = 5
+DEAL:4:2D-Test::   Number of this MPI processes = 4
+DEAL:4:2D-Test::   Number of cells = 256
+DEAL:4:2D-Test::   Number of levels = 4
+DEAL:4:2D-Test::   Number of global levels = 4
+DEAL:4:2D-Test::   Triangulation checksum = 0
+DEAL:4:2D-Test::   points = 0.230000 0.230000    0.230000 0.770000    0.770000 0.230000    0.770000 0.770000    1.10000 1.10000    
+DEAL:4:2D-Test::   >>> Reached end of test <<<
+
diff --git a/tests/mpi/p4est_3d_find_point_owner_rank_01.cc b/tests/mpi/p4est_3d_find_point_owner_rank_01.cc
new file mode 100644 (file)
index 0000000..1387b18
--- /dev/null
@@ -0,0 +1,42 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/utilities.h>
+
+#include "../tests.h"
+
+// contains the actual templated tests
+#include "p4est_find_point_owner_rank.h"
+
+int
+main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  MPILogInitAll                    log;
+
+
+  /*
+   * We craete a distrubuted mesh with three cells (tree roots), e.g., a
+   * hyper_L. We want to find the mpi rank of a single fixed point. On all
+   * processes we must find the same owner.
+   */
+  const std::vector<Point<3>> point(1, Point<3>(0.23, 0.77, 0.77));
+
+  deallog.push("3D-Test");
+  test<3>(point);
+  deallog.pop();
+}
diff --git a/tests/mpi/p4est_3d_find_point_owner_rank_01.mpirun=1.with_p4est=true.with_p4est_with_search_local=true.output b/tests/mpi/p4est_3d_find_point_owner_rank_01.mpirun=1.with_p4est=true.with_p4est_with_search_local=true.output
new file mode 100644 (file)
index 0000000..82161b9
--- /dev/null
@@ -0,0 +1,13 @@
+
+DEAL:0:3D-Test:: Point search on subdivided hyper cube...
+DEAL:0:3D-Test::   Number of MPI processes = 1
+DEAL:0:3D-Test::   Number of this MPI processes = 0
+DEAL:0:3D-Test::   Number of cells = 4096
+DEAL:0:3D-Test::   Number of levels = 4
+DEAL:0:3D-Test::   Number of global levels = 4
+DEAL:0:3D-Test::   Triangulation checksum = 544659457
+DEAL:0:3D-Test::   points = 0.230000 0.770000 0.770000
+DEAL:0:3D-Test::
+DEAL:0:3D-Test::   All point owner ranks found on all processes equal each other. TEST PASSED.
+DEAL:0:3D-Test::
+DEAL:0:3D-Test::   >>> Reached end of test <<<
diff --git a/tests/mpi/p4est_3d_find_point_owner_rank_01.mpirun=3.with_p4est=true.with_p4est_with_search_local=true.output b/tests/mpi/p4est_3d_find_point_owner_rank_01.mpirun=3.with_p4est=true.with_p4est_with_search_local=true.output
new file mode 100644 (file)
index 0000000..298f2e2
--- /dev/null
@@ -0,0 +1,35 @@
+
+DEAL:0:3D-Test:: Point search on subdivided hyper cube...
+DEAL:0:3D-Test::   Number of MPI processes = 3
+DEAL:0:3D-Test::   Number of this MPI processes = 0
+DEAL:0:3D-Test::   Number of cells = 4096
+DEAL:0:3D-Test::   Number of levels = 4
+DEAL:0:3D-Test::   Number of global levels = 4
+DEAL:0:3D-Test::   Triangulation checksum = 544659457
+DEAL:0:3D-Test::   points = 0.230000 0.770000 0.770000
+DEAL:0:3D-Test::
+DEAL:0:3D-Test::   All point owner ranks found on all processes equal each other. TEST PASSED.
+DEAL:0:3D-Test::
+DEAL:0:3D-Test::   >>> Reached end of test <<<
+
+DEAL:1:3D-Test:: Point search on subdivided hyper cube...
+DEAL:1:3D-Test::   Number of MPI processes = 3
+DEAL:1:3D-Test::   Number of this MPI processes = 1
+DEAL:1:3D-Test::   Number of cells = 4096
+DEAL:1:3D-Test::   Number of levels = 4
+DEAL:1:3D-Test::   Number of global levels = 4
+DEAL:1:3D-Test::   Triangulation checksum = 0
+DEAL:1:3D-Test::   points = 0.230000 0.770000 0.770000
+DEAL:1:3D-Test::   >>> Reached end of test <<<
+
+
+DEAL:2:3D-Test:: Point search on subdivided hyper cube...
+DEAL:2:3D-Test::   Number of MPI processes = 3
+DEAL:2:3D-Test::   Number of this MPI processes = 2
+DEAL:2:3D-Test::   Number of cells = 4096
+DEAL:2:3D-Test::   Number of levels = 4
+DEAL:2:3D-Test::   Number of global levels = 4
+DEAL:2:3D-Test::   Triangulation checksum = 0
+DEAL:2:3D-Test::   points = 0.230000 0.770000 0.770000
+DEAL:2:3D-Test::   >>> Reached end of test <<<
+
diff --git a/tests/mpi/p4est_3d_find_point_owner_rank_01.mpirun=5.with_p4est=true.with_p4est_with_search_local=true.output b/tests/mpi/p4est_3d_find_point_owner_rank_01.mpirun=5.with_p4est=true.with_p4est_with_search_local=true.output
new file mode 100644 (file)
index 0000000..3aa3096
--- /dev/null
@@ -0,0 +1,57 @@
+
+DEAL:0:3D-Test:: Point search on subdivided hyper cube...
+DEAL:0:3D-Test::   Number of MPI processes = 5
+DEAL:0:3D-Test::   Number of this MPI processes = 0
+DEAL:0:3D-Test::   Number of cells = 4096
+DEAL:0:3D-Test::   Number of levels = 4
+DEAL:0:3D-Test::   Number of global levels = 4
+DEAL:0:3D-Test::   Triangulation checksum = 544659457
+DEAL:0:3D-Test::   points = 0.230000 0.770000 0.770000    
+DEAL:0:3D-Test::
+DEAL:0:3D-Test::   All point owner ranks found on all processes equal each other. TEST PASSED.
+DEAL:0:3D-Test::
+DEAL:0:3D-Test::   >>> Reached end of test <<<
+
+DEAL:1:3D-Test:: Point search on subdivided hyper cube...
+DEAL:1:3D-Test::   Number of MPI processes = 5
+DEAL:1:3D-Test::   Number of this MPI processes = 1
+DEAL:1:3D-Test::   Number of cells = 4096
+DEAL:1:3D-Test::   Number of levels = 4
+DEAL:1:3D-Test::   Number of global levels = 4
+DEAL:1:3D-Test::   Triangulation checksum = 0
+DEAL:1:3D-Test::   points = 0.230000 0.770000 0.770000    
+DEAL:1:3D-Test::   >>> Reached end of test <<<
+
+
+DEAL:2:3D-Test:: Point search on subdivided hyper cube...
+DEAL:2:3D-Test::   Number of MPI processes = 5
+DEAL:2:3D-Test::   Number of this MPI processes = 2
+DEAL:2:3D-Test::   Number of cells = 4096
+DEAL:2:3D-Test::   Number of levels = 4
+DEAL:2:3D-Test::   Number of global levels = 4
+DEAL:2:3D-Test::   Triangulation checksum = 0
+DEAL:2:3D-Test::   points = 0.230000 0.770000 0.770000    
+DEAL:2:3D-Test::   >>> Reached end of test <<<
+
+
+DEAL:3:3D-Test:: Point search on subdivided hyper cube...
+DEAL:3:3D-Test::   Number of MPI processes = 5
+DEAL:3:3D-Test::   Number of this MPI processes = 3
+DEAL:3:3D-Test::   Number of cells = 4096
+DEAL:3:3D-Test::   Number of levels = 4
+DEAL:3:3D-Test::   Number of global levels = 4
+DEAL:3:3D-Test::   Triangulation checksum = 0
+DEAL:3:3D-Test::   points = 0.230000 0.770000 0.770000    
+DEAL:3:3D-Test::   >>> Reached end of test <<<
+
+
+DEAL:4:3D-Test:: Point search on subdivided hyper cube...
+DEAL:4:3D-Test::   Number of MPI processes = 5
+DEAL:4:3D-Test::   Number of this MPI processes = 4
+DEAL:4:3D-Test::   Number of cells = 4096
+DEAL:4:3D-Test::   Number of levels = 4
+DEAL:4:3D-Test::   Number of global levels = 4
+DEAL:4:3D-Test::   Triangulation checksum = 0
+DEAL:4:3D-Test::   points = 0.230000 0.770000 0.770000    
+DEAL:4:3D-Test::   >>> Reached end of test <<<
+
diff --git a/tests/mpi/p4est_3d_find_point_owner_rank_02.cc b/tests/mpi/p4est_3d_find_point_owner_rank_02.cc
new file mode 100644 (file)
index 0000000..b5cb75d
--- /dev/null
@@ -0,0 +1,54 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/utilities.h>
+
+#include "../tests.h"
+
+// contains the actual templated tests
+#include "p4est_find_point_owner_rank.h"
+
+
+int
+main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  MPILogInitAll                    log;
+
+
+  /*
+   * We craete a distrubuted mesh with three cells (tree roots), e.g., a
+   * hyper_L. We want to find the mpi rank of a set of fixed points. On all
+   * processes we must find the same owner ranks.
+   */
+  std::vector<Point<3>> points;
+  points.emplace_back(Point<3>(0.23, 0.23, 0.23));
+  points.emplace_back(Point<3>(0.23, 0.77, 0.23));
+  points.emplace_back(Point<3>(0.77, 0.23, 0.23));
+  points.emplace_back(Point<3>(0.77, 0.77, 0.23));
+  points.emplace_back(Point<3>(0.23, 0.23, 0.77));
+  points.emplace_back(Point<3>(0.23, 0.77, 0.77));
+  points.emplace_back(Point<3>(0.77, 0.23, 0.77));
+  points.emplace_back(Point<3>(0.77, 0.77, 0.77));
+
+  // this point is outside the mesh
+  points.emplace_back(Point<3>(1.1, 1.1, 1.1));
+
+  deallog.push("3D-Test");
+  test<3>(points);
+  deallog.pop();
+}
diff --git a/tests/mpi/p4est_3d_find_point_owner_rank_02.mpirun=1.with_p4est=true.with_p4est_with_search_local=true.output b/tests/mpi/p4est_3d_find_point_owner_rank_02.mpirun=1.with_p4est=true.with_p4est_with_search_local=true.output
new file mode 100644 (file)
index 0000000..32d2da3
--- /dev/null
@@ -0,0 +1,13 @@
+
+DEAL:0:3D-Test:: Point search on subdivided hyper cube...
+DEAL:0:3D-Test::   Number of MPI processes = 1
+DEAL:0:3D-Test::   Number of this MPI processes = 0
+DEAL:0:3D-Test::   Number of cells = 4096
+DEAL:0:3D-Test::   Number of levels = 4
+DEAL:0:3D-Test::   Number of global levels = 4
+DEAL:0:3D-Test::   Triangulation checksum = 544659457
+DEAL:0:3D-Test::   points = 0.230000 0.230000 0.230000    0.230000 0.770000 0.230000    0.770000 0.230000 0.230000    0.770000 0.770000 0.230000    0.230000 0.230000 0.770000    0.230000 0.770000 0.770000    0.770000 0.230000 0.770000    0.770000 0.770000 0.770000    1.10000 1.10000 1.10000    
+DEAL:0:3D-Test::
+DEAL:0:3D-Test::   All point owner ranks found on all processes equal each other. TEST PASSED.
+DEAL:0:3D-Test::
+DEAL:0:3D-Test::   >>> Reached end of test <<<
diff --git a/tests/mpi/p4est_3d_find_point_owner_rank_02.mpirun=3.with_p4est=true.with_p4est_with_search_local=true.output b/tests/mpi/p4est_3d_find_point_owner_rank_02.mpirun=3.with_p4est=true.with_p4est_with_search_local=true.output
new file mode 100644 (file)
index 0000000..19f8b54
--- /dev/null
@@ -0,0 +1,35 @@
+
+DEAL:0:3D-Test:: Point search on subdivided hyper cube...
+DEAL:0:3D-Test::   Number of MPI processes = 3
+DEAL:0:3D-Test::   Number of this MPI processes = 0
+DEAL:0:3D-Test::   Number of cells = 4096
+DEAL:0:3D-Test::   Number of levels = 4
+DEAL:0:3D-Test::   Number of global levels = 4
+DEAL:0:3D-Test::   Triangulation checksum = 544659457
+DEAL:0:3D-Test::   points = 0.230000 0.230000 0.230000    0.230000 0.770000 0.230000    0.770000 0.230000 0.230000    0.770000 0.770000 0.230000    0.230000 0.230000 0.770000    0.230000 0.770000 0.770000    0.770000 0.230000 0.770000    0.770000 0.770000 0.770000    1.10000 1.10000 1.10000    
+DEAL:0:3D-Test::
+DEAL:0:3D-Test::   All point owner ranks found on all processes equal each other. TEST PASSED.
+DEAL:0:3D-Test::
+DEAL:0:3D-Test::   >>> Reached end of test <<<
+
+DEAL:1:3D-Test:: Point search on subdivided hyper cube...
+DEAL:1:3D-Test::   Number of MPI processes = 3
+DEAL:1:3D-Test::   Number of this MPI processes = 1
+DEAL:1:3D-Test::   Number of cells = 4096
+DEAL:1:3D-Test::   Number of levels = 4
+DEAL:1:3D-Test::   Number of global levels = 4
+DEAL:1:3D-Test::   Triangulation checksum = 0
+DEAL:1:3D-Test::   points = 0.230000 0.230000 0.230000    0.230000 0.770000 0.230000    0.770000 0.230000 0.230000    0.770000 0.770000 0.230000    0.230000 0.230000 0.770000    0.230000 0.770000 0.770000    0.770000 0.230000 0.770000    0.770000 0.770000 0.770000    1.10000 1.10000 1.10000    
+DEAL:1:3D-Test::   >>> Reached end of test <<<
+
+
+DEAL:2:3D-Test:: Point search on subdivided hyper cube...
+DEAL:2:3D-Test::   Number of MPI processes = 3
+DEAL:2:3D-Test::   Number of this MPI processes = 2
+DEAL:2:3D-Test::   Number of cells = 4096
+DEAL:2:3D-Test::   Number of levels = 4
+DEAL:2:3D-Test::   Number of global levels = 4
+DEAL:2:3D-Test::   Triangulation checksum = 0
+DEAL:2:3D-Test::   points = 0.230000 0.230000 0.230000    0.230000 0.770000 0.230000    0.770000 0.230000 0.230000    0.770000 0.770000 0.230000    0.230000 0.230000 0.770000    0.230000 0.770000 0.770000    0.770000 0.230000 0.770000    0.770000 0.770000 0.770000    1.10000 1.10000 1.10000    
+DEAL:2:3D-Test::   >>> Reached end of test <<<
+
diff --git a/tests/mpi/p4est_3d_find_point_owner_rank_02.mpirun=5.with_p4est=true.with_p4est_with_search_local=true.output b/tests/mpi/p4est_3d_find_point_owner_rank_02.mpirun=5.with_p4est=true.with_p4est_with_search_local=true.output
new file mode 100644 (file)
index 0000000..b934189
--- /dev/null
@@ -0,0 +1,57 @@
+
+DEAL:0:3D-Test:: Point search on subdivided hyper cube...
+DEAL:0:3D-Test::   Number of MPI processes = 5
+DEAL:0:3D-Test::   Number of this MPI processes = 0
+DEAL:0:3D-Test::   Number of cells = 4096
+DEAL:0:3D-Test::   Number of levels = 4
+DEAL:0:3D-Test::   Number of global levels = 4
+DEAL:0:3D-Test::   Triangulation checksum = 544659457
+DEAL:0:3D-Test::   points = 0.230000 0.230000 0.230000    0.230000 0.770000 0.230000    0.770000 0.230000 0.230000    0.770000 0.770000 0.230000    0.230000 0.230000 0.770000    0.230000 0.770000 0.770000    0.770000 0.230000 0.770000    0.770000 0.770000 0.770000    1.10000 1.10000 1.10000    
+DEAL:0:3D-Test::
+DEAL:0:3D-Test::   All point owner ranks found on all processes equal each other. TEST PASSED.
+DEAL:0:3D-Test::
+DEAL:0:3D-Test::   >>> Reached end of test <<<
+
+DEAL:1:3D-Test:: Point search on subdivided hyper cube...
+DEAL:1:3D-Test::   Number of MPI processes = 5
+DEAL:1:3D-Test::   Number of this MPI processes = 1
+DEAL:1:3D-Test::   Number of cells = 4096
+DEAL:1:3D-Test::   Number of levels = 4
+DEAL:1:3D-Test::   Number of global levels = 4
+DEAL:1:3D-Test::   Triangulation checksum = 0
+DEAL:1:3D-Test::   points = 0.230000 0.230000 0.230000    0.230000 0.770000 0.230000    0.770000 0.230000 0.230000    0.770000 0.770000 0.230000    0.230000 0.230000 0.770000    0.230000 0.770000 0.770000    0.770000 0.230000 0.770000    0.770000 0.770000 0.770000    1.10000 1.10000 1.10000    
+DEAL:1:3D-Test::   >>> Reached end of test <<<
+
+
+DEAL:2:3D-Test:: Point search on subdivided hyper cube...
+DEAL:2:3D-Test::   Number of MPI processes = 5
+DEAL:2:3D-Test::   Number of this MPI processes = 2
+DEAL:2:3D-Test::   Number of cells = 4096
+DEAL:2:3D-Test::   Number of levels = 4
+DEAL:2:3D-Test::   Number of global levels = 4
+DEAL:2:3D-Test::   Triangulation checksum = 0
+DEAL:2:3D-Test::   points = 0.230000 0.230000 0.230000    0.230000 0.770000 0.230000    0.770000 0.230000 0.230000    0.770000 0.770000 0.230000    0.230000 0.230000 0.770000    0.230000 0.770000 0.770000    0.770000 0.230000 0.770000    0.770000 0.770000 0.770000    1.10000 1.10000 1.10000    
+DEAL:2:3D-Test::   >>> Reached end of test <<<
+
+
+DEAL:3:3D-Test:: Point search on subdivided hyper cube...
+DEAL:3:3D-Test::   Number of MPI processes = 5
+DEAL:3:3D-Test::   Number of this MPI processes = 3
+DEAL:3:3D-Test::   Number of cells = 4096
+DEAL:3:3D-Test::   Number of levels = 4
+DEAL:3:3D-Test::   Number of global levels = 4
+DEAL:3:3D-Test::   Triangulation checksum = 0
+DEAL:3:3D-Test::   points = 0.230000 0.230000 0.230000    0.230000 0.770000 0.230000    0.770000 0.230000 0.230000    0.770000 0.770000 0.230000    0.230000 0.230000 0.770000    0.230000 0.770000 0.770000    0.770000 0.230000 0.770000    0.770000 0.770000 0.770000    1.10000 1.10000 1.10000    
+DEAL:3:3D-Test::   >>> Reached end of test <<<
+
+
+DEAL:4:3D-Test:: Point search on subdivided hyper cube...
+DEAL:4:3D-Test::   Number of MPI processes = 5
+DEAL:4:3D-Test::   Number of this MPI processes = 4
+DEAL:4:3D-Test::   Number of cells = 4096
+DEAL:4:3D-Test::   Number of levels = 4
+DEAL:4:3D-Test::   Number of global levels = 4
+DEAL:4:3D-Test::   Triangulation checksum = 0
+DEAL:4:3D-Test::   points = 0.230000 0.230000 0.230000    0.230000 0.770000 0.230000    0.770000 0.230000 0.230000    0.770000 0.770000 0.230000    0.230000 0.230000 0.770000    0.230000 0.770000 0.770000    0.770000 0.230000 0.770000    0.770000 0.770000 0.770000    1.10000 1.10000 1.10000    
+DEAL:4:3D-Test::   >>> Reached end of test <<<
+
diff --git a/tests/mpi/p4est_find_point_owner_rank.h b/tests/mpi/p4est_find_point_owner_rank.h
new file mode 100644 (file)
index 0000000..b414a5b
--- /dev/null
@@ -0,0 +1,153 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+
+// STL
+#include <algorithm> // std::adjacent_find
+
+/*
+ * Test if the rank of point owners returned by all MPI processes is the same.
+ */
+template <int dim>
+void
+test(const std::vector<Point<dim>> &points)
+{
+  deallog << " Point search on subdivided hyper cube..." << std::endl;
+
+  parallel::distributed::Triangulation<dim> triangulation(
+    MPI_COMM_WORLD,
+    typename Triangulation<dim>::MeshSmoothing(
+      Triangulation<dim>::smoothing_on_refinement |
+      Triangulation<dim>::smoothing_on_coarsening),
+    typename parallel::distributed::Triangulation<dim>::Settings(
+      parallel::distributed::Triangulation<
+        dim>::communicate_vertices_to_p4est));
+  GridGenerator::subdivided_hyper_cube(triangulation, 2);
+  triangulation.refine_global(3);
+
+  deallog << "   Number of MPI processes = "
+          << Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) << std::endl;
+  deallog << "   Number of this MPI processes = "
+          << Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) << std::endl;
+  deallog << "   Number of cells = " << triangulation.n_global_active_cells()
+          << std::endl;
+  deallog << "   Number of levels = " << triangulation.n_levels() << std::endl;
+  deallog << "   Number of global levels = " << triangulation.n_global_levels()
+          << std::endl;
+  const unsigned int checksum = triangulation.get_checksum();
+  deallog << "   Triangulation checksum = " << checksum << std::endl;
+
+  deallog << "   points = ";
+  for (const auto &point : points)
+    deallog << point << "    ";
+
+  deallog << std::endl;
+
+  ////////////////////////////////////////////////////////////
+  // test stuff
+  std::vector<types::subdomain_id> point_owner_ranks =
+    triangulation.find_point_owner_rank(points);
+
+  // Gather all point_owner_ranks found from all MPI ranks and compare the
+  // results on root_process. They must all be equal.
+  const unsigned int                            root_process{0};
+  std::vector<std::vector<types::subdomain_id>> all_ranks_found_on_each_process{
+    Utilities::MPI::gather(MPI_COMM_WORLD, point_owner_ranks, root_process)};
+
+  if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == root_process)
+    {
+      if (std::adjacent_find(all_ranks_found_on_each_process.begin(),
+                             all_ranks_found_on_each_process.end(),
+                             std::not_equal_to<>()) ==
+          all_ranks_found_on_each_process.end())
+        {
+          deallog << std::endl
+                  << "   All point owner ranks found on all processes equal "
+                     "each other. TEST PASSED."
+                  << std::endl
+                  << std::endl;
+        }
+      else
+        {
+          deallog
+            << std::endl
+            << "   Some point owner ranks found on all processes do not equal "
+               "each other. Check output and/or the file "
+               "tests/mpi/p4est_find_point_owner_rank.h. TEST FAILED."
+            << std::endl
+            << std::endl;
+        }
+    }
+  ////////////////////////////////////////////////////////////
+
+  deallog << "   >>> Reached end of test <<<" << std::endl;
+}
+
+
+/*
+ * Vector input version
+ */
+template <int dim>
+void
+check_error_on_invalid_mesh(const Point<dim> &point)
+{
+  deallog << " Point search on hyper shell (manifold attached)..." << std::endl;
+
+  parallel::distributed::Triangulation<dim> triangulation(
+    MPI_COMM_WORLD,
+    typename Triangulation<dim>::MeshSmoothing(
+      Triangulation<dim>::smoothing_on_refinement |
+      Triangulation<dim>::smoothing_on_coarsening),
+    typename parallel::distributed::Triangulation<dim>::Settings(
+      parallel::distributed::Triangulation<
+        dim>::communicate_vertices_to_p4est));
+
+  GridGenerator::hyper_shell(triangulation,
+                             /* center is origin */ Point<dim>(),
+                             /*        inner_radius */ 1,
+                             /*        outer_radius */ 2);
+  triangulation.refine_global(3);
+
+  deallog << "   Number of MPI processes = "
+          << Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) << std::endl;
+  deallog << "   Number of this MPI processes = "
+          << Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) << std::endl;
+  deallog << "   Number of cells = " << triangulation.n_global_active_cells()
+          << std::endl;
+  deallog << "   Number of levels = " << triangulation.n_levels() << std::endl;
+  deallog << "   Number of global levels = " << triangulation.n_global_levels()
+          << std::endl;
+  const unsigned int checksum = triangulation.get_checksum();
+  deallog << "   Triangulation checksum = " << checksum << std::endl;
+  deallog << "   point = " << point << std::endl;
+
+  deallog << std::endl;
+
+  ////////////////////////////////////////////////////////////
+  // This should result in an error.
+  types::subdomain_id point_owner_rank =
+    triangulation.find_point_owner_rank(point);
+  ////////////////////////////////////////////////////////////
+}

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.