From: Martin Kronbichler Date: Mon, 28 May 2018 16:52:13 +0000 (+0200) Subject: Implement find_all_active_cells_around_point. X-Git-Tag: v9.1.0-rc1~1066^2~5 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7a44823745e820b32485c369c4e6df04d6827742;p=dealii.git Implement find_all_active_cells_around_point. --- diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index bf167c889d..2ac920356f 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -1164,6 +1164,40 @@ namespace GridTools cell_hint = typename Triangulation::active_cell_iterator(), const std::vector &marked_vertices = std::vector()); + /** + * A variant of the previous find_active_cell_around_point() function that, + * instead of returning only the first matching cell, identifies all cells + * around a point for a given tolerance level `tolerance` in terms of unit + * coordinates. More precisely, whenever the point returned by + * find_active_cell_around_point() is within the given tolerance from the + * surface of the unit cell, all corresponding neighbors are also + * identified, including the location of the point in unit coordinates on + * any of these cells. + * + * This function is useful e.g. for discontinuous function spaces where, for + * the case the given point `p` coincides with a vertex or an edge, several + * cells might hold independent values of the solution that get combined in + * some way in a user code. + * + * @author Niklas Fehn, Martin Kronbichler, 2018 + */ + template class MeshType, int spacedim> +# ifndef _MSC_VER + std::vector::active_cell_iterator, + Point>> +# else + std::vector>::type, + Point>> +# endif + find_all_active_cells_around_point( + const Mapping & mapping, + const MeshType &mesh, + const Point & p, + const double tolerance = 1e-12, + const std::vector & marked_vertices = std::vector()); + /** * Return a list of all descendants of the given cell that are active. For * example, if the current cell is once refined but none of its children are diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 76b38257d4..7a961e9f96 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -1493,8 +1493,144 @@ namespace GridTools } } } + + + + template class MeshType, int spacedim> +#ifndef _MSC_VER + std::pair::active_cell_iterator, + Point> +#else + std::pair< + typename dealii::internal:: + ActiveCellIterator>::type, + Point> +#endif + find_active_cell_around_point_tolerance( + const Mapping & mapping, + const MeshType &mesh, + const Point & p, + const std::vector & marked_vertices, + const double tolerance) + { + 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 the given tolerance + double best_distance = tolerance; + 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; + } } // namespace + + template class MeshType, int spacedim> #ifndef _MSC_VER typename MeshType::active_cell_iterator @@ -1512,6 +1648,7 @@ namespace GridTools } + template class MeshType, int spacedim> #ifndef _MSC_VER std::pair::active_cell_iterator, Point> @@ -1525,118 +1662,159 @@ namespace GridTools 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) + return find_active_cell_around_point_tolerance( + mapping, mesh, p, marked_vertices, 1e-10); + } + + + + template class MeshType, int spacedim> +#ifndef _MSC_VER + std::vector::active_cell_iterator, + Point>> +#else + std::vector>::type, + Point>> +#endif + find_all_active_cells_around_point(const Mapping & mapping, + const MeshType &mesh, + const Point & p, + const double tolerance, + const std::vector &marked_vertices) + { + // first use the result of the single point function as a guess. In order + // not to make the other find_all_active_cells_around_point more expensive + // and avoid some additional logic there, we first start with one cell as + // given by that other function (that possibly goes through a larger set + // of cells) and later add a list of more cells as appropriate. + std::vector< + std::pair::active_cell_iterator, + Point>> + cells_and_points; + cells_and_points.push_back(find_active_cell_around_point_tolerance( + mapping, mesh, p, marked_vertices, tolerance)); + + if (!cells_and_points.empty()) { - typename std::set::const_iterator - cell = adjacent_cells.begin(), - endc = adjacent_cells.end(); - for (; cell != endc; ++cell) + // check if the given point is on the surface of the unit cell. if yes, + // need to find all neighbors + const Point unit_point = cells_and_points.front().second; + const auto my_cell = cells_and_points.front().first; + Tensor<1, dim> distance_to_center; + unsigned int n_dirs_at_threshold = 0; + unsigned int last_point_at_threshold = numbers::invalid_unsigned_int; + for (unsigned int d = 0; d < dim; ++d) { - try + distance_to_center[d] = std::abs(unit_point[d] - 0.5); + if (distance_to_center[d] > 0.5 - tolerance) { - const Point p_cell = - mapping.transform_real_to_unit_cell(*cell, p); + ++n_dirs_at_threshold; + last_point_at_threshold = d; + } + } - // 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))) + std::vector::active_cell_iterator> + cells_to_add; + // point is within face -> only need neighbor + if (n_dirs_at_threshold == 1) + { + unsigned int neighbor_index = + 2 * last_point_at_threshold + + (unit_point[last_point_at_threshold] > 0.5 ? 1 : 0); + if (!my_cell->at_boundary(neighbor_index)) + cells_to_add.push_back(my_cell->neighbor(neighbor_index)); + } + // corner point -> use all neighbors + else if (n_dirs_at_threshold == dim) + { + unsigned int local_vertex_index = 0; + for (unsigned int d = 0; d < dim; ++d) + local_vertex_index += (unit_point[d] > 0.5 ? 1 : 0) << d; + std::vector::active_cell_iterator> + cells = find_cells_adjacent_to_vertex( + mesh, my_cell->vertex_index(local_vertex_index)); + for (auto cell : cells) + if (cell != my_cell) + cells_to_add.push_back(cell); + } + // point on line in 3D: We cannot simply take the intersection between + // the two vertices of cels because of hanging nodes. So instead we + // list the vertices around both points and then select the + // appropriate cells according to the result of read_to_unit_cell + // below. + else if (n_dirs_at_threshold == 2) + { + std::pair vertex_indices[3]; + unsigned int count_vertex_indices = 0; + unsigned int free_direction = numbers::invalid_unsigned_int; + for (unsigned int d = 0; d < dim; ++d) + { + if (distance_to_center[d] > 0.5 - tolerance) { - found = true; - best_distance = dist; - best_level = (*cell)->level(); - best_cell = std::make_pair(*cell, p_cell); + vertex_indices[count_vertex_indices].first = d; + vertex_indices[count_vertex_indices].second = + unit_point[d] > 0.5 ? 1 : 0; + count_vertex_indices++; } + else + free_direction = d; } - catch ( - typename MappingQGeneric::ExcTransformationFailed - &) + + AssertDimension(count_vertex_indices, 2); + Assert(free_direction != numbers::invalid_unsigned_int, + ExcInternalError()); + + const unsigned int first_vertex = + (vertex_indices[0].second << vertex_indices[0].first) + + (vertex_indices[1].second << vertex_indices[1].first); + for (unsigned int d = 0; d < 2; ++d) { - // 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 + auto tentative_cells = find_cells_adjacent_to_vertex( + mesh, + my_cell->vertex_index(first_vertex + (d << free_direction))); + for (auto cell : tentative_cells) + { + bool cell_not_yet_present = true; + for (auto other_cell : cells_to_add) + if (cell == other_cell) + { + cell_not_yet_present = false; + break; + } + if (cell_not_yet_present) + cells_to_add.push_back(cell); + } } } - // 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) + const double original_distance_to_unit_cell = + GeometryInfo::distance_to_unit_cell(unit_point); + for (auto cell : cells_to_add) { - find_active_cell_around_point_internal( - mesh, searched_cells, adjacent_cells); + if (cell != my_cell) + try + { + const Point p_unit = + mapping.transform_real_to_unit_cell(cell, p); + if (GeometryInfo::distance_to_unit_cell(p_unit) < + original_distance_to_unit_cell + tolerance) + cells_and_points.emplace_back(cell, p_unit); + } + catch (typename Mapping::ExcTransformationFailed &) + {} } } - - AssertThrow(best_cell.first.state() == IteratorState::valid, - ExcPointNotFound(p)); - - return best_cell; + std::sort( + cells_and_points.begin(), + cells_and_points.end(), + [](const std::pair::active_cell_iterator, + Point> &a, + const std::pair::active_cell_iterator, + Point> &b) { return a.first < b.first; }); + + return cells_and_points; } diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index 27c26eb78e..3158b88c8d 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -92,6 +92,17 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) deal_II_space_dimension>::active_cell_iterator &, const std::vector &); + template std::vector::active_cell_iterator, + Point>> + find_all_active_cells_around_point( + const Mapping &, + const Triangulation &, + const Point &, + const double, + const std::vector &); + template std::tuple::active_cell_iterator>, @@ -145,7 +156,22 @@ for (deal_II_space_dimension : SPACE_DIMENSIONS) find_active_cell_around_point( const Mapping &, const parallel::distributed::Triangulation &, - const Point &); + const Point &, + const std::vector &); + + std::vector>::type, + Point>> + find_all_active_cells_around_point( + const Mapping &, + const parallel::distributed::Triangulation &, + const Point &, + const double, + const std::vector &); template unsigned int GridTools::find_closest_vertex( const std::map> &vertices,