From fcbc08d4de547cb647b0d67b280a6701d2c5a7ce Mon Sep 17 00:00:00 2001 From: Giovanni Alzetta Date: Fri, 24 Nov 2017 00:30:49 +0100 Subject: [PATCH] Added distributed compute point locations and two tests --- .../changes/minor/20180119GiovanniAlzetta | 5 + include/deal.II/grid/grid_tools.h | 56 + source/grid/grid_tools.cc | 1055 ++++++++++++++++- source/grid/grid_tools.inst.in | 13 + .../distributed_compute_point_locations_01.cc | 132 +++ ...compute_point_locations_01.mpirun=1.output | 10 + .../distributed_compute_point_locations_02.cc | 229 ++++ ...h_mpi=true.with_p4est=true.mpirun=2.output | 37 + ...h_mpi=true.with_p4est=true.mpirun=3.output | 56 + 9 files changed, 1592 insertions(+), 1 deletion(-) create mode 100644 doc/news/changes/minor/20180119GiovanniAlzetta create mode 100644 tests/grid/distributed_compute_point_locations_01.cc create mode 100644 tests/grid/distributed_compute_point_locations_01.mpirun=1.output create mode 100644 tests/grid/distributed_compute_point_locations_02.cc create mode 100644 tests/grid/distributed_compute_point_locations_02.with_mpi=true.with_p4est=true.mpirun=2.output create mode 100644 tests/grid/distributed_compute_point_locations_02.with_mpi=true.with_p4est=true.mpirun=3.output diff --git a/doc/news/changes/minor/20180119GiovanniAlzetta b/doc/news/changes/minor/20180119GiovanniAlzetta new file mode 100644 index 0000000000..eab078b6c9 --- /dev/null +++ b/doc/news/changes/minor/20180119GiovanniAlzetta @@ -0,0 +1,5 @@ +New: New function GridTools::distributed_compute_point_locations ; similarly to GridTools::compute_point_locations , given +a vector of points, it returns vectors containing them, their reference position and the process owning them as it works +with shared and distributed meshes. +
+(Giovanni Alzetta, 2018/01/19) diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index b22ed339d1..670782e358 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -661,6 +661,62 @@ namespace GridTools const typename Triangulation::active_cell_iterator &cell_hint = typename Triangulation::active_cell_iterator()); + /** + * Given a @p cache and a list of + * @p local_points for each process, find the points lying on the locally owned + * part of the mesh and compute the quadrature rules for them. + * Distributed compute point locations is a function similar to + * GridTools::compute_point_locations but working for parallel::Triangulation + * objects and, unlike its serial version, also for a distributed triangulation + * (see parallel::distributed::Triangulation). + * + * @param[in] cache a GridTools::Cache object + * @param[in] local_points the array of points owned by the current process. Every + * process can have a different array of points which can be empty and not + * contained within the locally owned part of the triangulation + * @param[in] local_bbox the description of the locally owned part of the mesh made + * with bounding boxes. It can be obtained from + * GridTools::compute_mesh_predicate_bounding_box + * @param[out] tuple containing the quadrature information + * + * The elements of the output tuple are: + * - cells : a vector of cells of the all cells containing at + * least a point. + * - qpoints : a vector of vector of points; containing in @p qpoints[i] + * the reference positions of all points that fall within the cell @P cells[i] . + * - maps : a vector of vector of integers, containing the mapping between + * the numbering in qpoints (previous element of the tuple), and the vector + * of local points of the process owning the points. + * - points : a vector of vector of points. @p points[i][j] is the point in the + * real space corresponding. + * to @p qpoints[i][j] . Notice @p points are the points lying on the locally + * owned part of the mesh; thus these can be either copies of @p local_points + * or points received from other processes i.e. local_points for other processes + * - owners : a vector of vectors; @p owners[i][j] contains the rank of + * the process owning the point[i][j] (previous element of the tuple). + * + * The function uses the triangulation's mpi communicator: for this reason it + * throws an assert error if the Triangulation is not derived from + * parallel::Triangulation . + * + * In a serial execution the first three elements of the tuple are the same + * as in GridTools::compute_point_locations . + * + * @author Giovanni Alzetta, 2017-2018 + */ + template + std::tuple< + std::vector< typename Triangulation::active_cell_iterator >, + std::vector< std::vector< Point > >, + std::vector< std::vector< unsigned int > >, + std::vector< std::vector< Point > >, + std::vector< std::vector< unsigned int > > + > + distributed_compute_point_locations + (const GridTools::Cache &cache, + const std::vector > &local_points, + const std::vector< BoundingBox > &local_bbox); + /** * Return a map of index:Point, containing the used vertices of the * given `container`. The key of the returned map is the global index in the diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 6eef34a4b5..a5df281fec 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -15,6 +15,8 @@ #include #include +#include +#include #include #include @@ -59,7 +61,8 @@ #include #include #include - +#include +#include DEAL_II_NAMESPACE_OPEN @@ -1071,6 +1074,479 @@ namespace GridTools + template class MeshType, int spacedim> + unsigned int + find_closest_vertex (const MeshType &mesh, + const Point &p, + const std::vector &marked_vertices) + { + // first get the underlying + // triangulation from the + // mesh and determine vertices + // and used vertices + const Triangulation &tria = mesh.get_triangulation(); + + const std::vector< Point > &vertices = tria.get_vertices(); + + Assert ( tria.get_vertices().size() == marked_vertices.size() || marked_vertices.size() ==0, + ExcDimensionMismatch(tria.get_vertices().size(), marked_vertices.size())); + + // If p is an element of marked_vertices, + // and q is that of used_Vertices, + // the vector marked_vertices does NOT + // contain unused vertices if p implies q. + // I.e., if p is true q must be true + // (if p is false, q could be false or true). + // p implies q logic is encapsulated in ~p|q. + Assert( marked_vertices.size()==0 + || + std::equal( marked_vertices.begin(), + marked_vertices.end(), + tria.get_used_vertices().begin(), + [](bool p, bool q) + { + return !p || q; + }), + ExcMessage("marked_vertices should be a subset of used vertices in the triangulation " + "but marked_vertices contains one or more vertices that are not used vertices!") ); + + // In addition, if a vector bools + // is specified (marked_vertices) + // marking all the vertices which + // could be the potentially closest + // vertex to the point, use it instead + // of used vertices + const std::vector &used = + (marked_vertices.size()==0) ? tria.get_used_vertices() : marked_vertices; + + // At the beginning, the first + // used vertex is the closest one + std::vector::const_iterator first = + std::find(used.begin(), used.end(), true); + + // Assert that at least one vertex + // is actually used + Assert(first != used.end(), ExcInternalError()); + + unsigned int best_vertex = std::distance(used.begin(), first); + double best_dist = (p - vertices[best_vertex]).norm_square(); + + // For all remaining vertices, test + // whether they are any closer + for (unsigned int j = best_vertex+1; j < vertices.size(); j++) + if (used[j]) + { + double dist = (p - vertices[j]).norm_square(); + if (dist < best_dist) + { + best_vertex = j; + best_dist = dist; + } + } + + return best_vertex; + } + + + + template class MeshType, int spacedim> + unsigned int + find_closest_vertex (const Mapping &mapping, + const MeshType &mesh, + const Point &p, + const std::vector &marked_vertices) + { + // Take a shortcut in the simple case. + if (mapping.preserves_vertex_locations() == true) + return find_closest_vertex(mesh, p, marked_vertices); + + // first get the underlying + // triangulation from the + // mesh and determine vertices + // and used vertices + const Triangulation &tria = mesh.get_triangulation(); + + auto vertices = extract_used_vertices(tria, mapping); + + Assert ( tria.get_vertices().size() == marked_vertices.size() || marked_vertices.size() ==0, + ExcDimensionMismatch(tria.get_vertices().size(), marked_vertices.size())); + + // If p is an element of marked_vertices, + // and q is that of used_Vertices, + // the vector marked_vertices does NOT + // contain unused vertices if p implies q. + // I.e., if p is true q must be true + // (if p is false, q could be false or true). + // p implies q logic is encapsulated in ~p|q. + Assert( marked_vertices.size()==0 + || + std::equal( marked_vertices.begin(), + marked_vertices.end(), + tria.get_used_vertices().begin(), + [](bool p, bool q) + { + return !p || q; + }), + ExcMessage("marked_vertices should be a subset of used vertices in the triangulation " + "but marked_vertices contains one or more vertices that are not used vertices!") ); + + // Remove from the map unwanted elements. + if (marked_vertices.size()) + for (auto it = vertices.begin(); it != vertices.end(); ) + { + if (marked_vertices[it->first] == false) + { + vertices.erase(it++); + } + else + { + ++it; + } + } + + return find_closest_vertex(vertices, p); + } + + + + template class MeshType, int spacedim> +#ifndef _MSC_VER + std::vector::active_cell_iterator> +#else + std::vector >::type> +#endif + find_cells_adjacent_to_vertex(const MeshType &mesh, + const unsigned int vertex) + { + // make sure that the given vertex is + // an active vertex of the underlying + // triangulation + Assert(vertex < mesh.get_triangulation().n_vertices(), + ExcIndexRange(0,mesh.get_triangulation().n_vertices(),vertex)); + Assert(mesh.get_triangulation().get_used_vertices()[vertex], + ExcVertexNotUsed(vertex)); + + // use a set instead of a vector + // to ensure that cells are inserted only + // once + std::set >::type> adjacent_cells; + + typename dealii::internal::ActiveCellIterator >::type + cell = mesh.begin_active(), + endc = mesh.end(); + + // go through all active cells and look if the vertex is part of that cell + // + // in 1d, this is all we need to care about. in 2d/3d we also need to worry + // that the vertex might be a hanging node on a face or edge of a cell; in + // this case, we would want to add those cells as well on whose faces the + // vertex is located but for which it is not a vertex itself. + // + // getting this right is a lot simpler in 2d than in 3d. in 2d, a hanging + // node can only be in the middle of a face and we can query the neighboring + // cell from the current cell. on the other hand, in 3d a hanging node + // vertex can also be on an edge but there can be many other cells on + // this edge and we can not access them from the cell we are currently + // on. + // + // so, in the 3d case, if we run the algorithm as in 2d, we catch all + // those cells for which the vertex we seek is on a *subface*, but we + // miss the case of cells for which the vertex we seek is on a + // sub-edge for which there is no corresponding sub-face (because the + // immediate neighbor behind this face is not refined), see for example + // the bits/find_cells_adjacent_to_vertex_6 testcase. thus, if we + // haven't yet found the vertex for the current cell we also need to + // look at the mid-points of edges + // + // as a final note, deciding whether a neighbor is actually coarser is + // simple in the case of isotropic refinement (we just need to look at + // the level of the current and the neighboring cell). however, this + // isn't so simple if we have used anisotropic refinement since then + // the level of a cell is not indicative of whether it is coarser or + // not than the current cell. ultimately, we want to add all cells on + // which the vertex is, independent of whether they are coarser or + // finer and so in the 2d case below we simply add *any* *active* neighbor. + // in the worst case, we add cells multiple times to the adjacent_cells + // list, but std::set throws out those cells already entered + for (; cell != endc; ++cell) + { + for (unsigned int v = 0; v < GeometryInfo::vertices_per_cell; v++) + if (cell->vertex_index(v) == vertex) + { + // OK, we found a cell that contains + // the given vertex. We add it + // to the list. + adjacent_cells.insert(cell); + + // as explained above, in 2+d we need to check whether + // this vertex is on a face behind which there is a + // (possibly) coarser neighbor. if this is the case, + // then we need to also add this neighbor + if (dim >= 2) + for (unsigned int vface = 0; vface < dim; vface++) + { + const unsigned int face = + GeometryInfo::vertex_to_face[v][vface]; + + if (!cell->at_boundary(face) + && + cell->neighbor(face)->active()) + { + // there is a (possibly) coarser cell behind a + // face to which the vertex belongs. the + // vertex we are looking at is then either a + // vertex of that coarser neighbor, or it is a + // hanging node on one of the faces of that + // cell. in either case, it is adjacent to the + // vertex, so add it to the list as well (if + // the cell was already in the list then the + // std::set makes sure that we get it only + // once) + adjacent_cells.insert (cell->neighbor(face)); + } + } + + // in any case, we have found a cell, so go to the next cell + goto next_cell; + } + + // in 3d also loop over the edges + if (dim >= 3) + { + for (unsigned int e=0; e::lines_per_cell; ++e) + if (cell->line(e)->has_children()) + // the only place where this vertex could have been + // hiding is on the mid-edge point of the edge we + // are looking at + if (cell->line(e)->child(0)->vertex_index(1) == vertex) + { + adjacent_cells.insert(cell); + + // jump out of this tangle of nested loops + goto next_cell; + } + } + + // in more than 3d we would probably have to do the same as + // above also for even lower-dimensional objects + Assert (dim <= 3, ExcNotImplemented()); + + // move on to the next cell if we have found the + // vertex on the current one +next_cell: + ; + } + + // if this was an active vertex then there needs to have been + // at least one cell to which it is adjacent! + Assert (adjacent_cells.size() > 0, ExcInternalError()); + + // return the result as a vector, rather than the set we built above + return + std::vector >::type> + (adjacent_cells.begin(), adjacent_cells.end()); + } + + + + namespace + { + template class MeshType, int spacedim> + void find_active_cell_around_point_internal + (const MeshType &mesh, +#ifndef _MSC_VER + std::set::active_cell_iterator> &searched_cells, + std::set::active_cell_iterator> &adjacent_cells) +#else + std::set >::type> &searched_cells, + std::set >::type> &adjacent_cells) +#endif + { +#ifndef _MSC_VER + typedef typename MeshType::active_cell_iterator cell_iterator; +#else + typedef typename dealii::internal::ActiveCellIterator >::type cell_iterator; +#endif + + // update the searched cells + searched_cells.insert(adjacent_cells.begin(), adjacent_cells.end()); + // now we to collect all neighbors + // of the cells in adjacent_cells we + // have not yet searched. + std::set adjacent_cells_new; + + typename std::set::const_iterator + cell = adjacent_cells.begin(), + endc = adjacent_cells.end(); + for (; cell != endc; ++cell) + { + std::vector active_neighbors; + get_active_neighbors >(*cell, active_neighbors); + for (unsigned int i=0; i class MeshType, int spacedim> +#ifndef _MSC_VER + typename MeshType::active_cell_iterator +#else + typename dealii::internal::ActiveCellIterator >::type +#endif + find_active_cell_around_point (const MeshType &mesh, + const Point &p, + const std::vector &marked_vertices) + { + return + find_active_cell_around_point + (StaticMappingQ1::mapping, + mesh, p, marked_vertices).first; + } + + + template class MeshType, int spacedim> +#ifndef _MSC_VER + std::pair::active_cell_iterator, Point > +#else + std::pair >::type, Point > +#endif + find_active_cell_around_point (const Mapping &mapping, + const MeshType &mesh, + const Point &p, + const std::vector &marked_vertices) + { + typedef typename dealii::internal::ActiveCellIterator >::type active_cell_iterator; + + // The best distance is set to the + // maximum allowable distance from + // the unit cell; we assume a + // max. deviation of 1e-10 + double best_distance = 1e-10; + int best_level = -1; + std::pair > best_cell; + + // Find closest vertex and determine + // all adjacent cells + std::vector adjacent_cells_tmp + = find_cells_adjacent_to_vertex(mesh, + find_closest_vertex(mapping, mesh, p, marked_vertices)); + + // Make sure that we have found + // at least one cell adjacent to vertex. + Assert(adjacent_cells_tmp.size()>0, ExcInternalError()); + + // Copy all the cells into a std::set + std::set adjacent_cells (adjacent_cells_tmp.begin(), + adjacent_cells_tmp.end()); + std::set searched_cells; + + // Determine the maximal number of cells + // in the grid. + // As long as we have not found + // the cell and have not searched + // every cell in the triangulation, + // we keep on looking. + const unsigned int n_active_cells = mesh.get_triangulation().n_active_cells(); + bool found = false; + unsigned int cells_searched = 0; + while (!found && cells_searched < n_active_cells) + { + typename std::set::const_iterator + cell = adjacent_cells.begin(), + endc = adjacent_cells.end(); + for (; cell != endc; ++cell) + { + try + { + const Point p_cell = mapping.transform_real_to_unit_cell(*cell, p); + + // calculate the infinity norm of + // the distance vector to the unit cell. + const double dist = GeometryInfo::distance_to_unit_cell(p_cell); + + // We compare if the point is inside the + // unit cell (or at least not too far + // outside). If it is, it is also checked + // that the cell has a more refined state + if ((dist < best_distance) + || + ((dist == best_distance) + && + ((*cell)->level() > best_level))) + { + found = true; + best_distance = dist; + best_level = (*cell)->level(); + best_cell = std::make_pair(*cell, p_cell); + } + } + catch (typename MappingQGeneric::ExcTransformationFailed &) + { + // ok, the transformation + // failed presumably + // because the point we + // are looking for lies + // outside the current + // cell. this means that + // the current cell can't + // be the cell around the + // point, so just ignore + // this cell and move on + // to the next + } + } + + // update the number of cells searched + cells_searched += adjacent_cells.size(); + + // if the user provided a custom mask for vertices, + // terminate the search without trying to expand the search + // to all cells of the triangulation, as done below. + if (marked_vertices.size() > 0) + cells_searched = n_active_cells; + + // if we have not found the cell in + // question and have not yet searched every + // cell, we expand our search to + // all the not already searched neighbors of + // the cells in adjacent_cells. This is + // what find_active_cell_around_point_internal + // is for. + if (!found && cells_searched < n_active_cells) + { + find_active_cell_around_point_internal + (mesh, searched_cells, adjacent_cells); + } + } + + AssertThrow (best_cell.first.state() == IteratorState::valid, + ExcPointNotFound(p)); + + return best_cell; + } + + + template std::vector > > vertex_to_cell_centers_directions(const Triangulation &mesh, @@ -3315,6 +3791,580 @@ namespace GridTools + namespace internal + { + // Functions are needed for distributed compute point locations + namespace distributed_cptloc + { + // Hash function for cells; needed for unordered maps/multimaps + template < int dim, int spacedim> + struct cell_hash + { + std::size_t operator()(const typename Triangulation::active_cell_iterator &k) const + { + // Return active cell index, which is faster than CellId to compute + return k->active_cell_index(); + } + }; + + + + // Compute point locations; internal version which returns an unordered map + // The algorithm is the same as GridTools::compute_point_locations + template + std::unordered_map< typename Triangulation::active_cell_iterator, + std::pair >,std::vector >, cell_hash > + compute_point_locations_unmap(const GridTools::Cache &cache, + const std::vector > &points) + { + // How many points are here? + const unsigned int np = points.size(); + // Creating the output tuple + std::unordered_map< typename Triangulation::active_cell_iterator, + std::pair >,std::vector >, cell_hash > + cell_qpoint_map; + + // Now the easy case. + if (np==0) return cell_qpoint_map; + // We begin by finding the cell/transform of the first point + auto my_pair = GridTools::find_active_cell_around_point + (cache, points[0]); + + auto last_cell = cell_qpoint_map.emplace( + std::make_pair(my_pair.first, std::make_pair( + std::vector > {my_pair.second}, + std::vector {0}))); + // Now the second easy case. + if (np==1) return cell_qpoint_map; + // Computing the cell center and diameter + Point cell_center = my_pair.first->center(); + double cell_diameter = my_pair.first->diameter()* + (0.5 + std::numeric_limits::epsilon() ); + + // Cycle over all points left + for (unsigned int p=1; p< np; ++p) + { + // Checking if the point is close to the cell center, in which + // case calling find active cell with a cell hint + if ( cell_center.distance(points[p]) < cell_diameter ) + my_pair = GridTools::find_active_cell_around_point + (cache, points[p],last_cell.first->first); + else + my_pair = GridTools::find_active_cell_around_point + (cache, points[p]); + + if ( last_cell.first->first == my_pair.first) + { + last_cell.first->second.first.emplace_back(my_pair.second); + last_cell.first->second.second.emplace_back(p); + } + else + { + // Check if it is in another cell already found + last_cell = cell_qpoint_map.emplace(std::make_pair(my_pair.first, std::make_pair( + std::vector > {my_pair.second}, + std::vector {p}))); + + if ( last_cell.second == false ) + { + // Cell already present: adding the new point + last_cell.first->second.first.emplace_back(my_pair.second); + last_cell.first->second.second.emplace_back(p); + } + else + { + // New cell was added, updating center and diameter + cell_center = my_pair.first->center(); + cell_diameter = my_pair.first->diameter()* + (0.5 + std::numeric_limits::epsilon() ); + } + } + } + +#ifdef DEBUG + unsigned int qps = 0; + // The number of points in all + // the cells must be the same as + // the number of points we + // started off from. + for (const auto &m: cell_qpoint_map) + { + Assert(m.second.second.size() == + m.second.first.size(), + ExcDimensionMismatch(m.second.second.size(), + m.second.first.size())); + qps += m.second.second.size(); + } + Assert(qps == np, + ExcDimensionMismatch(qps, np)); +#endif + return cell_qpoint_map; + } + + + + // Merging the output means to add data to a previous output, here contained + // in output unmap: + // if the cell is already present: add information about new points + // if the cell is not present: add the cell with all information + // + // Notice we call "information" the data associated with a point of the sort: + // cell containing it, transformed point on reference cell, index, + // rank of the owner etc. + template + void + merge_cptloc_outputs( + std::unordered_map< typename Triangulation::active_cell_iterator, + std::tuple< + std::vector< Point >, + std::vector< unsigned int >, + std::vector< Point >, + std::vector< unsigned int > + >, + cell_hash> &output_unmap, + const std::vector::active_cell_iterator > &in_cells, + const std::vector< std::vector< Point > > &in_qpoints, + const std::vector< std::vector > &in_maps, + const std::vector< std::vector< Point > > &in_points, + const unsigned int in_rank + ) + { + // Adding cells, one by one + for (unsigned int c=0; c< in_cells.size(); ++c) + { + // Attempt to add a new cell with its relative data + auto current_c = output_unmap.emplace( + std::make_pair(in_cells[c], + std::make_tuple(in_qpoints[c], + in_maps[c], + in_points[c], + std::vector + (in_points[c].size(),in_rank)))); + // If the flag is false no new cell was added: + if ( current_c.second == false ) + { + // Cell in output map at current_c.first: + // Adding the information to it + auto &cell_qpts = std::get<0>(current_c.first->second); + auto &cell_maps = std::get<1>(current_c.first->second); + auto &cell_pts = std::get<2>(current_c.first->second); + auto &cell_ranks = std::get<3>(current_c.first->second); + cell_qpts.insert(cell_qpts.end(), + in_qpoints[c].begin(), + in_qpoints[c].end()); + cell_maps.insert(cell_maps.end(), + in_maps[c].begin(), + in_maps[c].end()); + cell_pts.insert(cell_pts.end(), + in_points[c].begin(), + in_points[c].end()); + std::vector< unsigned int > ranks_tmp(in_points[c].size(),in_rank); + cell_ranks.insert(cell_ranks.end(), + ranks_tmp.begin(), + ranks_tmp.end()); + } + } + } + + + + // This function initializes the output by calling compute point locations + // on local points; vector containing points which are probably local. + // Its output is then sorted in the following manner: + // - output unmap: points, with relative information, inside locally onwed cells, + // - ghost loc pts: points, with relative information, inside ghost cells, + // - classified pts: vector of all points returned in output map and ghost loc pts + // (these are stored as indices) + template + void + compute_and_classify_points( + const GridTools::Cache &cache, + const std::vector > &local_points, + const std::vector< unsigned int > &local_points_idx, + std::unordered_map< + typename Triangulation::active_cell_iterator, + std::tuple< + std::vector< Point >, + std::vector< unsigned int >, + std::vector< Point >, + std::vector< unsigned int > + >, + cell_hash> &output_unmap, + std::map< unsigned int, + std::tuple< + std::vector< CellId >, + std::vector< std::vector< Point > >, + std::vector< std::vector< unsigned int > >, + std::vector< std::vector< Point > > + > > &ghost_loc_pts, + std::vector< unsigned int > &classified_pts + ) + { + auto cpt_loc_pts = compute_point_locations_unmap(cache,local_points); + + // Alayzing the output discarding artificial cell + // and storing in the proper container locally owned and ghost cells + for (auto const &cell_tuples : cpt_loc_pts) + { + auto &cell_loc = cell_tuples.first; + auto &q_loc = std::get<0>(cell_tuples.second); + auto &indices_loc = std::get<1>(cell_tuples.second); + if (cell_loc->is_locally_owned() ) + { + // Point inside locally owned cell: storing all its data + std::vector < Point > cell_points(indices_loc.size()); + for (unsigned int i=0; i< indices_loc.size(); ++i) + { + // Adding the point to the cell points + cell_points[i] = local_points[indices_loc[i]]; + // Storing the index: notice indices loc refer to the local points + // vector, but we need to return the index with respect of + // the points owned by the current process + classified_pts.emplace_back(local_points_idx[indices_loc[i]]); + } + output_unmap.emplace(std::make_pair(cell_loc, + std::make_tuple(q_loc, + indices_loc, + cell_points, + std::vector + (indices_loc.size(),cell_loc->subdomain_id())))); + } + else if ( cell_loc->is_ghost() ) + { + // Point inside ghost cell: storing all its information and preparing + // it to be sent + std::vector < Point > cell_points(indices_loc.size()); + for (unsigned int i=0; i< indices_loc.size(); ++i) + { + cell_points[i] = local_points[indices_loc[i]]; + classified_pts.emplace_back(local_points_idx[indices_loc[i]]); + } + // Each key of the following map represent a process, + // each mapped value is a tuple containing the information to be sent: + // preparing the output for the owner, which has rank subdomain id + auto &map_tuple_owner = ghost_loc_pts[cell_loc->subdomain_id()]; + // To identify the cell on the other process we use the cell id + std::get<0>(map_tuple_owner).emplace_back(cell_loc->id()); + std::get<1>(map_tuple_owner).emplace_back(q_loc); + std::get<2>(map_tuple_owner).emplace_back(indices_loc); + std::get<3>(map_tuple_owner).emplace_back(cell_points); + } + // else: the cell is artificial, nothing to do + } + } + + + + // Given the map obtained from a communication, where the key is rank and the mapped + // value is a pair of (points,indices), calls compute point locations; its output + // is then merged with output tuple + // if check_owned is set to true only points + // lying inside locally onwed cells shall be merged, otherwise all points shall be merged. + template + void + compute_and_merge_from_map( + const GridTools::Cache &cache, + const std::map< unsigned int, + std::pair< + std::vector < Point >, + std::vector < unsigned int > > + > &map_pts, + std::unordered_map< typename Triangulation::active_cell_iterator, + std::tuple< + std::vector< Point >, + std::vector< unsigned int >, + std::vector< Point >, + std::vector< unsigned int > + >, + cell_hash> &output_unmap, + const bool &check_owned + ) + { + bool no_check = !check_owned; + + // rank and points is a pair: first rank, then a pair of vectors (points, indices) + for (auto const &rank_and_points : map_pts) + { + // Rewriting the contents of the map in human readable format + const auto &received_process = rank_and_points.first; + const auto &received_points = rank_and_points.second.first; + const auto &received_ranks = rank_and_points.second.second; + + // Initializing the vectors needed to store the result of compute point location + std::vector< typename Triangulation::active_cell_iterator > in_cell; + std::vector< std::vector< Point > > in_qpoints; + std::vector< std::vector< unsigned int > > in_maps; + std::vector< std::vector< Point > > in_points; + + auto cpt_loc_pts = compute_point_locations_unmap(cache,rank_and_points.second.first); + for (const auto &map_c_pt_idx: cpt_loc_pts) + { + // Human-readable variables: + const auto &proc_cell = map_c_pt_idx.first; + const auto &proc_qpoints = map_c_pt_idx.second.first; + const auto &proc_maps = map_c_pt_idx.second.second; + + // This is stored either if we're not checking if the cell is owned or + // if the cell is locally owned + if ( no_check || proc_cell->is_locally_owned() ) + { + in_cell.emplace_back(proc_cell); + in_qpoints.emplace_back(proc_qpoints); + // The other two vectors need to be built + unsigned int loc_size = proc_qpoints.size(); + std::vector< unsigned int > cell_maps(loc_size); + std::vector< Point > cell_points(loc_size); + for (unsigned int pt=0; pt + std::tuple< + std::vector< typename Triangulation::active_cell_iterator >, + std::vector< std::vector< Point > >, + std::vector< std::vector< unsigned int > >, + std::vector< std::vector< Point > >, + std::vector< std::vector< unsigned int > > + > + distributed_compute_point_locations + (const GridTools::Cache &cache, + const std::vector > &local_points, + const std::vector< BoundingBox > &local_bbox) + { +#ifndef DEAL_II_WITH_MPI + (void)cache; + (void)local_points; + (void)local_bbox; + Assert(false, ExcMessage("GridTools::distributed_compute_point_locations() requires MPI.")); + std::tuple< + std::vector< typename Triangulation::active_cell_iterator >, + std::vector< std::vector< Point > >, + std::vector< std::vector< unsigned int > >, + std::vector< std::vector< Point > >, + std::vector< std::vector< unsigned int > > + > tup; + return tup; +#else + // Recovering the mpi communicator used to create the triangulation + const auto &tria_mpi = + dynamic_cast< const parallel::Triangulation< dim, spacedim >*>(&cache.get_triangulation()); + // If the dynamic cast failed we can't recover the mpi communicator: throwing an assertion error + Assert(tria_mpi, ExcMessage("GridTools::distributed_compute_point_locations() requires a parallel triangulation.")); + auto mpi_communicator = tria_mpi->get_communicator(); + // Preparing the output tuple + std::tuple< + std::vector< typename Triangulation::active_cell_iterator >, + std::vector< std::vector< Point > >, + std::vector< std::vector< unsigned int > >, + std::vector< std::vector< Point > >, + std::vector< std::vector< unsigned int > > + > output_tuple; + + // Preparing the temporary unordered map + std::unordered_map< typename Triangulation::active_cell_iterator, + std::tuple< + std::vector< Point >, + std::vector< unsigned int >, + std::vector< Point >, + std::vector< unsigned int > + >, + internal::distributed_cptloc::cell_hash > + temporary_unmap; + + // Obtaining the global mesh description through an all to all communication + std::vector< std::vector< BoundingBox > > global_bounding_boxes; + global_bounding_boxes = Utilities::MPI::all_gather(mpi_communicator,local_bbox); + + // Step 1 (part 1): Using the bounding boxes to guess the owner of each points + // in local_points + unsigned int my_rank = Utilities::MPI::this_mpi_process(mpi_communicator); + + // Using global bounding boxes to guess/find owner/s of each point + std::tuple< std::vector< std::vector< unsigned int > >, std::map< unsigned int, unsigned int >, + std::map< unsigned int, std::vector< unsigned int > > > guessed_points; + guessed_points = + GridTools::guess_point_owner(global_bounding_boxes, local_points); + + // Preparing to call compute point locations on points which are/might be + // local + const auto &guess_loc_idx = std::get<0>(guessed_points)[my_rank]; + const unsigned int n_local_guess = guess_loc_idx.size(); + // Vector containing points which are probably local + std::vector< Point > guess_local_pts(n_local_guess); + for (unsigned int i=0; i, + std::vector< std::vector< Point > >, + std::vector< std::vector< unsigned int > >, + std::vector< std::vector< Point > > > > ghost_loc_pts; + // Vector containing indices of points lying either on locally owned + // cells or ghost cells, to avoid computing them more than once + std::vector< unsigned int > classified_pts; + + // Thread used to call compute point locations on guess local pts + Threads::Task + cpt_loc_tsk + = Threads::new_task ( + &internal::distributed_cptloc::compute_and_classify_points, + cache, + guess_local_pts, + guess_loc_idx, + temporary_unmap, + ghost_loc_pts, + classified_pts); + + // Step 1 (part 2): communicate point which are owned by a certain process + // Preparing the map with points whose owner is known with certainty: + const auto &other_owned_idx = std::get<1>(guessed_points); + std::map< + unsigned int, + std::pair< std::vector> , std::vector > > + other_owned_pts; + + for (const auto &indices: other_owned_idx) + if (indices.second != my_rank) + { + // Finding/adding in the map the current process + auto ¤t_pts = other_owned_pts[indices.second]; + current_pts.first.emplace_back(local_points[indices.first]); + current_pts.second.emplace_back(indices.first); + } + + // Communicating the points whose owner is sure + auto owned_rank_pts = Utilities::MPI::some_to_some(mpi_communicator,other_owned_pts); + // Waiting for part 1 to finish to avoid concurrency problems + cpt_loc_tsk.join(); + + // Step 2 (part 1): compute received points which are owned + Threads::Task + owned_pts_tsk + = Threads::new_task (&internal::distributed_cptloc::compute_and_merge_from_map, + cache, + owned_rank_pts, + temporary_unmap, + false); + + // Step 2 (part 2): communicate info on points lying on ghost cells + auto cpt_ghost = Utilities::MPI::some_to_some(mpi_communicator,ghost_loc_pts); + + // Step 3: construct vectors containing uncertain points i.e. those whose owner + // is known among few guesses + std::map< + unsigned int, + std::pair< std::vector < Point >, + std::vector > > + other_check_pts; + + const auto &other_check_idx = std::get<2>(guessed_points); + + // Points in classified pts need not to be communicated; + // sorting the array classified pts in order to use + // binary search when checking if the points needs to be + // communicated + // Notice classified pts is a vector of integer indexes + std::sort (classified_pts.begin(), classified_pts.end()); + + for (const auto &pt_to_guesses: other_check_idx) + { + if ( !std::binary_search( + classified_pts.begin(), classified_pts.end(),pt_to_guesses.first) ) + // The point wasn't found in ghost or locally owned cells: adding it to the map + for (unsigned int rank=0; rank(rank_vals.second); + unsigned int n_cells = cell_ids.size(); + std::vector< typename Triangulation::active_cell_iterator > + cell_iter(n_cells); + for (unsigned int c=0; c(rank_vals.second), + std::get<2>(rank_vals.second), + std::get<3>(rank_vals.second), + rank_vals.first); + } + + // Step 6: use compute point locations on the uncertain points and + // merge output + internal::distributed_cptloc::compute_and_merge_from_map( + cache, + check_pts, + temporary_unmap, + true); + + // Copying data from the unordered map to the tuple + // and returning output + unsigned int size_output = temporary_unmap.size(); + auto &out_cells = std::get<0>(output_tuple); + auto &out_qpoints = std::get<1>(output_tuple); + auto &out_maps = std::get<2>(output_tuple); + auto &out_points = std::get<3>(output_tuple); + auto &out_ranks = std::get<4>(output_tuple); + + out_cells.resize(size_output); + out_qpoints.resize(size_output); + out_maps.resize(size_output); + out_points.resize(size_output); + out_ranks.resize(size_output); + + unsigned int c = 0; + for (const auto &rank_and_tuple: temporary_unmap) + { + out_cells[c] = rank_and_tuple.first; + out_qpoints[c] = std::get<0>(rank_and_tuple.second); + out_maps[c] = std::get<1>(rank_and_tuple.second); + out_points[c] = std::get<2>(rank_and_tuple.second); + out_ranks[c] = std::get<3>(rank_and_tuple.second); + ++c; + } + + return output_tuple; +#endif + } + + template std::map > extract_used_vertices(const Triangulation &container, @@ -3332,6 +4382,7 @@ namespace GridTools return result; } + template unsigned int find_closest_vertex(const std::map > &vertices, @@ -3347,6 +4398,7 @@ namespace GridTools return id_and_v->first; } + template std::pair::active_cell_iterator, Point > find_active_cell_around_point(const Cache &cache, @@ -3443,6 +4495,7 @@ namespace GridTools #endif // DEAL_II_WITH_MPI } + } /* namespace GridTools */ diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index 22d8a2afdb..f8e55bc93b 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -75,6 +75,19 @@ for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS compute_point_locations(const Cache< deal_II_dimension, deal_II_space_dimension > &, const std::vector< Point< deal_II_space_dimension > > &, const typename Triangulation< deal_II_dimension, deal_II_space_dimension>::active_cell_iterator &); + + template + std::tuple< + std::vector< typename Triangulation< deal_II_dimension, deal_II_space_dimension>::active_cell_iterator >, + std::vector< std::vector< Point > >, + std::vector< std::vector< unsigned int > >, + std::vector< std::vector< Point > >, + std::vector< std::vector< unsigned int > > + > + distributed_compute_point_locations + (const Cache< deal_II_dimension, deal_II_space_dimension > &, + const std::vector< Point< deal_II_space_dimension > > &, + const std::vector< BoundingBox< deal_II_space_dimension > > &); \} #endif diff --git a/tests/grid/distributed_compute_point_locations_01.cc b/tests/grid/distributed_compute_point_locations_01.cc new file mode 100644 index 0000000000..8633f21be2 --- /dev/null +++ b/tests/grid/distributed_compute_point_locations_01.cc @@ -0,0 +1,132 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017-2018 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// Test GridTools::distributed_compute_point_locations for the serial case + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace dealii; + +template +void test_compute_pt_loc(unsigned int n_points) +{ + MPI_Comm mpi_communicator = MPI_COMM_WORLD; + deallog << "Testing for dim = " << dim << std::endl; + deallog << "Testing on: " << n_points << " points." << std::endl; + + // Creating a grid in the square [0,1]x[0,1] + parallel::distributed::Triangulation tria(mpi_communicator); + GridGenerator::hyper_cube(tria); + tria.refine_global(std::max(6-dim,2)); + + // Creating the random points + std::vector> points; + + for (size_t i=0; i()); + + // Initializing the cache + GridTools::Cache cache(tria); + + // Computing the description of the locally owned part of the mesh + IteratorFilters::LocallyOwnedCell locally_owned_cell_predicate; + std::vector< BoundingBox > local_bbox = GridTools::compute_mesh_predicate_bounding_box + (cache.get_triangulation(), locally_owned_cell_predicate, + 1, true, 4); // These options should be passed + // Using the distributed version of compute point location + + // Using the distributed version + auto output_tuple = distributed_compute_point_locations(cache,points,local_bbox); + // Testing in serial against the serial version + auto cell_qpoint_map = GridTools::compute_point_locations(cache,points); + + auto &serial_cells = std::get<0>(cell_qpoint_map); + auto &serial_qpoints = std::get<1>(cell_qpoint_map); + size_t n_cells = std::get<0>(output_tuple).size(); + + deallog << "Points found in " << n_cells << " cells" << std::endl; + + // testing if the result coincides with + // the serial one + for (unsigned int c=0; c(output_tuple)[c]; + auto &quad = std::get<1>(output_tuple)[c]; + auto &local_map = std::get<2>(output_tuple)[c]; + auto &original_points = std::get<3>(output_tuple)[c]; + auto &ranks = std::get<4>(output_tuple)[c]; + + auto pos_cell = std::find(serial_cells.begin(),serial_cells.end(),cell); + for (auto r: ranks) + if (r!=0) + deallog << "ERROR: rank is not 0 but " << std::to_string(r) << std::endl; + + if (pos_cell == serial_cells.end()) + deallog << "ERROR: cell not found" << std::endl; + else + { + auto serial_cell_idx = pos_cell - serial_cells.begin(); + if ( original_points.size() != serial_qpoints[serial_cell_idx].size()) + deallog << "ERROR: in the number of points for cell" << std::to_string(serial_cell_idx) << std::endl; + if ( quad.size() != serial_qpoints[serial_cell_idx].size()) + deallog << "ERROR: in the number of points for cell" << std::to_string(serial_cell_idx) << std::endl; + + unsigned int pt_num = 0; + for (const auto &p_idx: local_map) + { + auto serial_pt_pos = std::find(local_map.begin(),local_map.end(),p_idx); + auto serial_pt_idx = serial_pt_pos-local_map.begin(); + if ( serial_pt_pos == local_map.end()) + deallog << "ERROR: point index not found for " << std::to_string(serial_pt_idx) << std::endl; + else + { + if ( (original_points[pt_num] - points[p_idx]).norm() > 1e-12 ) + { + deallog << "ERROR: Point in serial : " << points[p_idx] << " Point in distributed: " << original_points[pt_num] << std::endl; + } + + if ( (quad[pt_num] - serial_qpoints[serial_cell_idx][serial_pt_idx]).norm() > 1e-10 ) + { + deallog << " ERROR: Transformation of qpoint to point is not correct" << std::endl; + deallog << "qpoint in serial : " << quad[pt_num] << " Point in distributed: " << serial_qpoints[serial_cell_idx][serial_pt_idx] << std::endl; + } + } + ++pt_num; + } + } + } + + deallog << "Test finished" << std::endl; +} + +int main (int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); + MPILogInitAll log; + + deallog << "Deal.II distributed_compute_point_locations:" << std::endl; + test_compute_pt_loc<2>(100); + test_compute_pt_loc<3>(200); +} diff --git a/tests/grid/distributed_compute_point_locations_01.mpirun=1.output b/tests/grid/distributed_compute_point_locations_01.mpirun=1.output new file mode 100644 index 0000000000..c8ceeada23 --- /dev/null +++ b/tests/grid/distributed_compute_point_locations_01.mpirun=1.output @@ -0,0 +1,10 @@ + +DEAL:0::Deal.II distributed_compute_point_locations: +DEAL:0::Testing for dim = 2 +DEAL:0::Testing on: 100 points. +DEAL:0::Points found in 81 cells +DEAL:0::Test finished +DEAL:0::Testing for dim = 3 +DEAL:0::Testing on: 200 points. +DEAL:0::Points found in 170 cells +DEAL:0::Test finished diff --git a/tests/grid/distributed_compute_point_locations_02.cc b/tests/grid/distributed_compute_point_locations_02.cc new file mode 100644 index 0000000000..08208c6eb3 --- /dev/null +++ b/tests/grid/distributed_compute_point_locations_02.cc @@ -0,0 +1,229 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017-2018 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// Test GridTools::distributed_compute_point_locations for the parallel case: +// Inside a distributed hypercube there's a shared sphere: +// call distributed point locations on the sphere's cells centers and check +// the result. + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace dealii; + +template +void test_compute_pt_loc(unsigned int ref_cube, unsigned int ref_sphere) +{ + MPI_Comm mpi_communicator = MPI_COMM_WORLD; + unsigned int n_procs = Utilities::MPI::n_mpi_processes(mpi_communicator); + unsigned int my_rank = Utilities::MPI::this_mpi_process(mpi_communicator); + + deallog << "Testing for dim = " << dim << " on " << n_procs << " processes" << std::endl; + deallog << "Cube refinements: " << ref_cube << std::endl; + deallog << "Sphere refinements:" << ref_sphere << std::endl; + + // Initializing and refining meshes + parallel::distributed::Triangulation cube(mpi_communicator); + GridGenerator::hyper_cube(cube); + cube.refine_global(ref_cube); + + parallel::shared::Triangulation sphere(mpi_communicator); + Point sphere_center; + // Defining center and radius + for (unsigned int i=0; i < dim; ++i) + sphere_center[i] = 0.47 - i*0.05; + double radius = 0.4 - dim*0.05; + GridGenerator::hyper_ball( sphere, sphere_center,radius); + static SphericalManifold surface_description(sphere_center); + sphere.set_manifold(0, surface_description); + sphere.refine_global(ref_sphere); + + deallog << "Sphere center:" << sphere_center << std::endl; + deallog << "Sphere radius:" << radius << std::endl; + + // Initializing the cache + GridTools::Cache cache(cube); + + // Centers of locally owned cells + std::vector> loc_owned_points; + // Building by hand the output of distributed points (see the function's + // description for more details, this is in fact compute point location + // code with the addition of rank storing) + std::vector< typename Triangulation::active_cell_iterator > computed_cells; + std::vector< std::vector > > computed_qpoints; + std::vector< std::vector > > computed_points; + std::vector< std::vector< unsigned int > > computed_ranks; + + unsigned int computed_pts = 0; + for (auto cell: sphere.active_cell_iterators()) + { + // The points we consider are the cell centers + auto center_pt = cell->center(); + // Store the point only if it is inside a locally owned sphere cell + if (cell->subdomain_id()==my_rank) + loc_owned_points.emplace_back(center_pt); + // Find the cube cell where center pt lies + auto my_pair = GridTools::find_active_cell_around_point + (cache, center_pt); + // If it is inside a locally owned cell it shall be returned + // from distributed compute point locations + if ( my_pair.first->is_locally_owned() ) + { + computed_pts++; + auto cells_it = + std::find(computed_cells.begin(),computed_cells.end(),my_pair.first); + + if ( cells_it == computed_cells.end() ) + { + // Cell not found: adding a new cell + computed_cells.emplace_back(my_pair.first); + computed_qpoints.emplace_back(1, my_pair.second); + computed_points.emplace_back(1, center_pt); + computed_ranks.emplace_back(1, cell->subdomain_id()); + } + else + { + // Cell found: just adding the point index and qpoint to the list + unsigned int current_cell = cells_it - computed_cells.begin(); + computed_qpoints[current_cell].emplace_back(my_pair.second); + computed_points[current_cell].emplace_back(center_pt); + computed_ranks[current_cell].emplace_back(cell->subdomain_id()); + } + } + } + + // Computing bounding boxes describing the locally owned part of the mesh + IteratorFilters::LocallyOwnedCell locally_owned_cell_predicate; + std::vector< BoundingBox > local_bbox = GridTools::compute_mesh_predicate_bounding_box + (cache.get_triangulation(), locally_owned_cell_predicate, + 1, true, 4); + + // Using the distributed version of compute point location + auto output_tuple = distributed_compute_point_locations + (cache,loc_owned_points,local_bbox); + deallog << "Comparing results" << std::endl; + const auto &output_cells = std::get<0>(output_tuple); + const auto &output_qpoints = std::get<1>(output_tuple); + const auto &output_points = std::get<3>(output_tuple); + const auto &output_ranks = std::get<4>(output_tuple); + + // Comparing the output with the previously computed computed result + bool test_passed = true; + if (output_cells.size() != computed_cells.size() ) + { + test_passed = false; + deallog << "ERROR: non-matching number of cell found" << std::endl; + } + + unsigned int output_computed_pts = 0; + for (unsigned int c=0; c< output_cells.size(); c++) + { + output_computed_pts += output_points[c].size(); + const auto &cell = output_cells[c]; + auto cell_it = + std::find(computed_cells.begin(),computed_cells.end(),cell); + if ( cell_it == computed_cells.end() ) + { + deallog << "ERROR: active cell " << cell->active_cell_index() << " not found" << std::endl; + test_passed = false; + } + else + { + unsigned int c_cell = cell_it - computed_cells.begin(); + if (output_points[c].size() != computed_points[c_cell].size() ) + { + test_passed = false; + deallog << "ERROR: non-matching number of points for cell " << cell->active_cell_index() << std::endl; + deallog << "Distributed compute point location output:" << std::endl; + for (unsigned int pt_idx=0; pt_idx< output_points[c].size(); pt_idx++) + deallog << output_points[c][pt_idx] << " from process " << output_ranks[c][pt_idx] << " to " << my_rank << std::endl; + deallog << "Expected points:" << std::endl; + for (unsigned int pt_idx=0; pt_idx< computed_points[c_cell].size(); pt_idx++) + deallog << computed_points[c_cell][pt_idx] << std::endl; + + } + else + { + // Checking if the points inside are the same + for (unsigned int pt_idx=0; pt_idx< output_points[c].size(); pt_idx++) + { + const auto &pt = output_points[c][pt_idx]; + auto pt_it = + std::find(computed_points[c_cell].begin(),computed_points[c_cell].end(),pt); + if ( pt_it == computed_points[c_cell].end() ) + { + deallog << "ERROR: point " << pt << " not found" << std::endl; + test_passed = false; + } + else + { + unsigned int c_pt = pt_it - computed_points[c_cell].begin(); + // Checking the value of the transformed point + if ( (output_qpoints[c][pt_idx] - computed_qpoints[c_cell][c_pt]).norm() > 1e-12 ) + { + // Cell not found: adding a new cell + deallog << "ERROR: qpoint " << c_pt << " not matching" << std::endl; + test_passed = false; + } + // Checking the rank of the owner + if ( output_ranks[c][pt_idx] != computed_ranks[c_cell][c_pt]) + { + // Cell not found: adding a new cell + deallog << "ERROR: rank of point " << c_pt << " not matching" << std::endl; + test_passed = false; + } + } + } + } + } + } + + + + if (output_computed_pts != computed_pts) + { + deallog << "ERROR: the number of points is different from expected: " << std::endl; + deallog << "Number of locally computed points: " << computed_pts << std::endl; + deallog << "Number of points from distributed: " << output_computed_pts << std::endl; + } + + if (test_passed) + deallog << "Test passed" << std::endl; + else + deallog << "Test FAILED" << std::endl; +} + +int main (int argc, char *argv[]) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1); + MPILogInitAll log; + + deallog << "Deal.II distributed_compute_point_locations:" << std::endl; + deallog << "Test on parallel setting 2D:" << std::endl; + test_compute_pt_loc<2>(3,3); + deallog << "Test on parallel setting 3D:" << std::endl; + test_compute_pt_loc<3>(3,2); +} diff --git a/tests/grid/distributed_compute_point_locations_02.with_mpi=true.with_p4est=true.mpirun=2.output b/tests/grid/distributed_compute_point_locations_02.with_mpi=true.with_p4est=true.mpirun=2.output new file mode 100644 index 0000000000..f42d916e44 --- /dev/null +++ b/tests/grid/distributed_compute_point_locations_02.with_mpi=true.with_p4est=true.mpirun=2.output @@ -0,0 +1,37 @@ + +DEAL:0::Deal.II distributed_compute_point_locations: +DEAL:0::Test on parallel setting 2D: +DEAL:0::Testing for dim = 2 on 2 processes +DEAL:0::Cube refinements: 3 +DEAL:0::Sphere refinements:3 +DEAL:0::Sphere center:0.470000 0.420000 +DEAL:0::Sphere radius:0.300000 +DEAL:0::Comparing results +DEAL:0::Test passed +DEAL:0::Test on parallel setting 3D: +DEAL:0::Testing for dim = 3 on 2 processes +DEAL:0::Cube refinements: 3 +DEAL:0::Sphere refinements:2 +DEAL:0::Sphere center:0.470000 0.420000 0.370000 +DEAL:0::Sphere radius:0.250000 +DEAL:0::Comparing results +DEAL:0::Test passed + +DEAL:1::Deal.II distributed_compute_point_locations: +DEAL:1::Test on parallel setting 2D: +DEAL:1::Testing for dim = 2 on 2 processes +DEAL:1::Cube refinements: 3 +DEAL:1::Sphere refinements:3 +DEAL:1::Sphere center:0.470000 0.420000 +DEAL:1::Sphere radius:0.300000 +DEAL:1::Comparing results +DEAL:1::Test passed +DEAL:1::Test on parallel setting 3D: +DEAL:1::Testing for dim = 3 on 2 processes +DEAL:1::Cube refinements: 3 +DEAL:1::Sphere refinements:2 +DEAL:1::Sphere center:0.470000 0.420000 0.370000 +DEAL:1::Sphere radius:0.250000 +DEAL:1::Comparing results +DEAL:1::Test passed + diff --git a/tests/grid/distributed_compute_point_locations_02.with_mpi=true.with_p4est=true.mpirun=3.output b/tests/grid/distributed_compute_point_locations_02.with_mpi=true.with_p4est=true.mpirun=3.output new file mode 100644 index 0000000000..371fd7d9b1 --- /dev/null +++ b/tests/grid/distributed_compute_point_locations_02.with_mpi=true.with_p4est=true.mpirun=3.output @@ -0,0 +1,56 @@ + +DEAL:0::Deal.II distributed_compute_point_locations: +DEAL:0::Test on parallel setting 2D: +DEAL:0::Testing for dim = 2 on 3 processes +DEAL:0::Cube refinements: 3 +DEAL:0::Sphere refinements:3 +DEAL:0::Sphere center:0.470000 0.420000 +DEAL:0::Sphere radius:0.300000 +DEAL:0::Comparing results +DEAL:0::Test passed +DEAL:0::Test on parallel setting 3D: +DEAL:0::Testing for dim = 3 on 3 processes +DEAL:0::Cube refinements: 3 +DEAL:0::Sphere refinements:2 +DEAL:0::Sphere center:0.470000 0.420000 0.370000 +DEAL:0::Sphere radius:0.250000 +DEAL:0::Comparing results +DEAL:0::Test passed + +DEAL:1::Deal.II distributed_compute_point_locations: +DEAL:1::Test on parallel setting 2D: +DEAL:1::Testing for dim = 2 on 3 processes +DEAL:1::Cube refinements: 3 +DEAL:1::Sphere refinements:3 +DEAL:1::Sphere center:0.470000 0.420000 +DEAL:1::Sphere radius:0.300000 +DEAL:1::Comparing results +DEAL:1::Test passed +DEAL:1::Test on parallel setting 3D: +DEAL:1::Testing for dim = 3 on 3 processes +DEAL:1::Cube refinements: 3 +DEAL:1::Sphere refinements:2 +DEAL:1::Sphere center:0.470000 0.420000 0.370000 +DEAL:1::Sphere radius:0.250000 +DEAL:1::Comparing results +DEAL:1::Test passed + + +DEAL:2::Deal.II distributed_compute_point_locations: +DEAL:2::Test on parallel setting 2D: +DEAL:2::Testing for dim = 2 on 3 processes +DEAL:2::Cube refinements: 3 +DEAL:2::Sphere refinements:3 +DEAL:2::Sphere center:0.470000 0.420000 +DEAL:2::Sphere radius:0.300000 +DEAL:2::Comparing results +DEAL:2::Test passed +DEAL:2::Test on parallel setting 3D: +DEAL:2::Testing for dim = 3 on 3 processes +DEAL:2::Cube refinements: 3 +DEAL:2::Sphere refinements:2 +DEAL:2::Sphere center:0.470000 0.420000 0.370000 +DEAL:2::Sphere radius:0.250000 +DEAL:2::Comparing results +DEAL:2::Test passed + -- 2.39.5