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()
# P4EST_WITH_MPI
# P4EST_WITH_ZLIB
# P4EST_WITH_VTK_BINARY
+# P4EST_WITH_SEARCH_LOCAL
# P4EST_VERSION
# P4EST_VERSION_MAJOR
# P4EST_VERSION_MINOR
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:
--- /dev/null
+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 <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>
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 <>
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
};
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]);
};
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 <p8est_ghost.h>
#endif
-
DEAL_II_NAMESPACE_OPEN
#ifdef DEAL_II_WITH_P4EST
* 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
};
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.
*
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.
*
* 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;
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.
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
};
/**
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.
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;
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
#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>
}
+# 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
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)
+ 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(
{
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());
}
// 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),
// 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),
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),
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>,
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,
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,
// 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
}
}
- // 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);
+ 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()
// 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;
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,
// 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;
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;
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();
}
+ 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()
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);
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);
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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();
+}
--- /dev/null
+
+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 <<<
--- /dev/null
+
+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 <<<
+
--- /dev/null
+
+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 <<<
+
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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();
+}
--- /dev/null
+
+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 <<<
--- /dev/null
+
+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 <<<
+
--- /dev/null
+
+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 <<<
+
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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();
+}
--- /dev/null
+
+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 <<<
--- /dev/null
+
+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 <<<
+
--- /dev/null
+
+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 <<<
+
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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();
+}
--- /dev/null
+
+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 <<<
--- /dev/null
+
+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 <<<
+
--- /dev/null
+
+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 <<<
+
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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);
+ ////////////////////////////////////////////////////////////
+}
projects / dealii.git / commitdiff