From: Simon, Dr. Konrad Date: Sun, 23 May 2021 08:40:37 +0000 (+0200) Subject: Implement communication-free owner rank search (p4est interface). X-Git-Tag: v9.4.0-rc1~939^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9791aa1564a26dabf74920ef683b82f7a2ea2cfd;p=dealii.git Implement communication-free owner rank search (p4est interface). --- diff --git a/cmake/configure/configure_p4est.cmake b/cmake/configure/configure_p4est.cmake index cc9c3c0016..1f0465aca1 100644 --- a/cmake/configure/configure_p4est.cmake +++ b/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() diff --git a/cmake/modules/FindP4EST.cmake b/cmake/modules/FindP4EST.cmake index eb87d2f47e..a91ef26876 100644 --- a/cmake/modules/FindP4EST.cmake +++ b/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: diff --git a/doc/news/changes/major/20210924Simon b/doc/news/changes/major/20210924Simon new file mode 100644 index 0000000000..04418f7814 --- /dev/null +++ b/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. +
+(Konrad Simon, 2021/09/24) diff --git a/include/deal.II/distributed/p4est_wrappers.h b/include/deal.II/distributed/p4est_wrappers.h index 21410aa792..bf63f43697 100644 --- a/include/deal.II/distributed/p4est_wrappers.h +++ b/include/deal.II/distributed/p4est_wrappers.h @@ -26,12 +26,14 @@ # include # include # include +# include # include # include # include # include # include # include +# include # include # include @@ -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]); }; diff --git a/include/deal.II/distributed/tria.h b/include/deal.II/distributed/tria.h index 51f22b9712..dd4f34e903 100644 --- a/include/deal.II/distributed/tria.h +++ b/include/deal.II/distributed/tria.h @@ -50,7 +50,6 @@ # include #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 &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> &points)` + * (requires p4est v2.2 and higher). Please see the documentation of + * `find_point_owner_rank(const std::vector> &points)`. + */ + types::subdomain_id + find_point_owner_rank(const Point &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 + find_point_owner_rank(const std::vector> &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. diff --git a/source/distributed/p4est_wrappers.cc b/source/distributed/p4est_wrappers.cc index 996a0fcc32..85540da99b 100644 --- a/source/distributed/p4est_wrappers.cc +++ b/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 void diff --git a/source/distributed/tria.cc b/source/distributed/tria.cc index 648aab7efb..e53164f288 100644 --- a/source/distributed/tria.cc +++ b/source/distributed/tria.cc @@ -16,6 +16,7 @@ #include #include +#include #include #include @@ -567,6 +568,583 @@ namespace } +# ifdef P4EST_SEARCH_LOCAL + template + class PartitionSearch + { + public: + PartitionSearch() + { + Assert(dim > 1, ExcNotImplemented()); + } + + PartitionSearch(const PartitionSearch &other) = delete; + + PartitionSearch & + operator=(const PartitionSearch &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::forest *forest, + typename internal::p4est::types::topidx which_tree, + typename internal::p4est::types::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::forest * forest, + typename internal::p4est::types::topidx which_tree, + typename internal::p4est::types::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::forest * forest, + typename internal::p4est::types::topidx which_tree, + typename internal::p4est::types::quadrant *quadrant, + const typename internal::p4est::types::quadrant_coord + quad_length_on_level); + + void + initialize_mapping(); + + Point + map_real_to_unit_cell(const Point &p) const; + + bool + is_in_this_quadrant(const Point &p) const; + + private: + std::vector> cell_vertices; + + /** + * Matrix holds coefficients mapping from this physical cell to unit + * cell. + */ + FullMatrix 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 + PartitionSearch::local_quadrant_fn( + typename internal::p4est::types::forest * forest, + typename internal::p4est::types::topidx which_tree, + typename internal::p4est::types::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 *this_object = + reinterpret_cast *>(forest->user_pointer); + + // Avoid p4est macros, instead do bitshifts manually with fixed size types + const typename internal::p4est::types::quadrant_coord + quad_length_on_level = + 1 << (static_cast::quadrant_coord>( + (dim == 2 ? P4EST_MAXLEVEL : P8EST_MAXLEVEL)) - + static_cast::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 + PartitionSearch::local_point_fn( + typename internal::p4est::types::forest *forest, + typename internal::p4est::types::topidx /* which_tree */, + typename internal::p4est::types::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 *this_object = + reinterpret_cast *>(forest->user_pointer); + + // point with rank as double pointer + double *this_point_dptr = static_cast(point); + + Point this_point = + (dim == 2 ? Point(this_point_dptr[0], this_point_dptr[1]) : + Point(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(rank_begin); + + // stop recursion. + return /* false */ 0; + } + + + + template + bool + PartitionSearch::QuadrantData::is_in_this_quadrant( + const Point &p) const + { + const Point p_ref = map_real_to_unit_cell(p); + + return GeometryInfo::is_inside_unit_cell(p_ref); + } + + + + template + Point + PartitionSearch::QuadrantData::map_real_to_unit_cell( + const Point &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 p_out; + + if (dim == 2) + { + for (unsigned int alpha = 0; + alpha < GeometryInfo::vertices_per_cell; + ++alpha) + { + const Point &p_ref = + GeometryInfo::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::vertices_per_cell; + ++alpha) + { + const Point &p_ref = + GeometryInfo::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 + PartitionSearch::QuadrantData::QuadrantData() + : cell_vertices(GeometryInfo::vertices_per_cell) + , quadrant_mapping_matrix(GeometryInfo::vertices_per_cell, + GeometryInfo::vertices_per_cell) + , are_vertices_initialized(false) + , is_reference_mapping_initialized(false) + {} + + + + template + void + PartitionSearch::QuadrantData::initialize_mapping() + { + Assert( + are_vertices_initialized, + dealii::ExcMessage( + "Cell vertices must be initialized before the cell mapping can be filled.")); + + FullMatrix point_matrix(GeometryInfo::vertices_per_cell, + GeometryInfo::vertices_per_cell); + + if (dim == 2) + { + for (unsigned int alpha = 0; + alpha < GeometryInfo::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::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::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::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::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::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::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::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::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::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::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::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::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::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::limit_level_difference_at_vertices) : smooth_grid, false) +# ifdef DEBUG + // Always communicate vertices to p4est in debug mode + , settings( + static_cast(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 + bool + Triangulation::are_vertices_communicated_to_p4est() const + { + return settings & + Triangulation::communicate_vertices_to_p4est; + } + + + template void Triangulation::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::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::cell_iterator; typename std::map, @@ -2009,8 +2594,8 @@ namespace parallel v < GeometryInfo::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::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> 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 + types::subdomain_id + Triangulation::find_point_owner_rank(const Point &p) + { + // Call the other function + std::vector> point{p}; + std::vector owner = find_point_owner_rank(point); + + return owner[0]; + } + + + + template + std::vector + Triangulation::find_point_owner_rank( + const std::vector> &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(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 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 &p = points[i]; + // get a non-const view of the array + double *this_sc_point = + static_cast(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::search_partition( + parallel_forest, + /* execute quadrant function when leaving quadrant */ + static_cast(false), + &PartitionSearch::local_quadrant_fn, + &PartitionSearch::local_point_fn, + point_sc_array); + + // copy the points found to an std::array + std::vector 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(sc_array_index_ssize_t(point_sc_array, i)); + owner_rank[i] = static_cast(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 void Triangulation::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::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::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::gloidx> previous_global_first_quadrant; @@ -2856,8 +3554,8 @@ namespace parallel PartitionWeights 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 + bool + Triangulation<1, spacedim>::are_vertices_communicated_to_p4est() const + { + Assert(false, ExcNotImplemented()); + return false; + } + + + template void Triangulation<1, spacedim>::update_cell_relations() diff --git a/tests/distributed_grids/coarse_grid_common.h b/tests/distributed_grids/coarse_grid_common.h index a72d7ca901..3deab18dd8 100644 --- a/tests/distributed_grids/coarse_grid_common.h +++ b/tests/distributed_grids/coarse_grid_common.h @@ -35,6 +35,11 @@ void write_vtk(const parallel::distributed::Triangulation &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/codim_01.cc b/tests/mpi/codim_01.cc index 65b8e95962..ac15dffe55 100644 --- a/tests/mpi/codim_01.cc +++ b/tests/mpi/codim_01.cc @@ -38,6 +38,11 @@ void write_vtk(const parallel::distributed::Triangulation &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 index 0000000000..4da881e38d --- /dev/null +++ b/tests/mpi/p4est_2d_find_point_owner_rank_01.cc @@ -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 +#include +#include + +#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(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 index 0000000000..13b37b6605 --- /dev/null +++ b/tests/mpi/p4est_2d_find_point_owner_rank_01.mpirun=1.with_p4est=true.with_p4est_with_search_local=true.output @@ -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 index 0000000000..44f9155678 --- /dev/null +++ b/tests/mpi/p4est_2d_find_point_owner_rank_01.mpirun=3.with_p4est=true.with_p4est_with_search_local=true.output @@ -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 index 0000000000..3373add304 --- /dev/null +++ b/tests/mpi/p4est_2d_find_point_owner_rank_01.mpirun=5.with_p4est=true.with_p4est_with_search_local=true.output @@ -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 index 0000000000..19b82ee2af --- /dev/null +++ b/tests/mpi/p4est_2d_find_point_owner_rank_02.cc @@ -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 +#include +#include + +#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> 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 index 0000000000..6c9cd9e899 --- /dev/null +++ b/tests/mpi/p4est_2d_find_point_owner_rank_02.mpirun=1.with_p4est=true.with_p4est_with_search_local=true.output @@ -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 index 0000000000..718bd08b4c --- /dev/null +++ b/tests/mpi/p4est_2d_find_point_owner_rank_02.mpirun=3.with_p4est=true.with_p4est_with_search_local=true.output @@ -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 index 0000000000..ed2fd7c8d0 --- /dev/null +++ b/tests/mpi/p4est_2d_find_point_owner_rank_02.mpirun=5.with_p4est=true.with_p4est_with_search_local=true.output @@ -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 index 0000000000..1387b1836c --- /dev/null +++ b/tests/mpi/p4est_3d_find_point_owner_rank_01.cc @@ -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 +#include +#include + +#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(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 index 0000000000..82161b9a1a --- /dev/null +++ b/tests/mpi/p4est_3d_find_point_owner_rank_01.mpirun=1.with_p4est=true.with_p4est_with_search_local=true.output @@ -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 index 0000000000..298f2e2219 --- /dev/null +++ b/tests/mpi/p4est_3d_find_point_owner_rank_01.mpirun=3.with_p4est=true.with_p4est_with_search_local=true.output @@ -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 index 0000000000..3aa30967c8 --- /dev/null +++ b/tests/mpi/p4est_3d_find_point_owner_rank_01.mpirun=5.with_p4est=true.with_p4est_with_search_local=true.output @@ -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 index 0000000000..b5cb75d0b3 --- /dev/null +++ b/tests/mpi/p4est_3d_find_point_owner_rank_02.cc @@ -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 +#include +#include + +#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> 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 index 0000000000..32d2da35b8 --- /dev/null +++ b/tests/mpi/p4est_3d_find_point_owner_rank_02.mpirun=1.with_p4est=true.with_p4est_with_search_local=true.output @@ -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 index 0000000000..19f8b54e13 --- /dev/null +++ b/tests/mpi/p4est_3d_find_point_owner_rank_02.mpirun=3.with_p4est=true.with_p4est_with_search_local=true.output @@ -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 index 0000000000..b9341893f9 --- /dev/null +++ b/tests/mpi/p4est_3d_find_point_owner_rank_02.mpirun=5.with_p4est=true.with_p4est_with_search_local=true.output @@ -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 index 0000000000..b414a5be8b --- /dev/null +++ b/tests/mpi/p4est_find_point_owner_rank.h @@ -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 +#include +#include + +#include + +#include +#include +#include +#include + +// STL +#include // std::adjacent_find + +/* + * Test if the rank of point owners returned by all MPI processes is the same. + */ +template +void +test(const std::vector> &points) +{ + deallog << " Point search on subdivided hyper cube..." << std::endl; + + parallel::distributed::Triangulation triangulation( + MPI_COMM_WORLD, + typename Triangulation::MeshSmoothing( + Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening), + typename parallel::distributed::Triangulation::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 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> 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 +void +check_error_on_invalid_mesh(const Point &point) +{ + deallog << " Point search on hyper shell (manifold attached)..." << std::endl; + + parallel::distributed::Triangulation triangulation( + MPI_COMM_WORLD, + typename Triangulation::MeshSmoothing( + Triangulation::smoothing_on_refinement | + Triangulation::smoothing_on_coarsening), + typename parallel::distributed::Triangulation::Settings( + parallel::distributed::Triangulation< + dim>::communicate_vertices_to_p4est)); + + GridGenerator::hyper_shell(triangulation, + /* center is origin */ Point(), + /* 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); + //////////////////////////////////////////////////////////// +}