Merge pull request #12341 from konsim83/feature_p4est_find_partition

Feature p4est find partition - Find MPI ranks of point owners in distributed meshes

cmake/configure/configure_p4est.cmake

@@ -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()
 
 

cmake/modules/FindP4EST.cmake

@@ -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:

doc/news/changes/major/20210924Simon

@@ -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)

include/deal.II/distributed/p4est_wrappers.h

@@ -26,12 +26,14 @@
 #  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]);
     };
 
 

include/deal.II/distributed/tria.h

@@ -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.

source/distributed/p4est_wrappers.cc

@@ -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