bool found_cell = false;
bool approx_cell = false;
- unsigned int closest_vertex_index = 0;
+ unsigned int closest_vertex_index = 0;
+ // ensure closest vertex index is a marked one, otherwise cell (with vertex
+ // 0) might be found even though it is not marked. This is only relevant if
+ // searching with rtree, using find_closest_vertex already can manage not
+ // finding points
+ if (marked_vertices.size() && !used_vertices_rtree.empty())
+ {
+ const auto itr =
+ std::find(marked_vertices.begin(), marked_vertices.end(), true);
+ Assert(itr != marked_vertices.end(),
+ dealii::ExcMessage("No vertex has been marked!"));
+ closest_vertex_index = std::distance(marked_vertices.begin(), itr);
+ }
+
Tensor<1, spacedim> vertex_to_point;
auto current_cell = cell_hint;
+ // check whether cell has at least one marked vertex
+ const auto cell_marked = [&mesh, &marked_vertices](const auto &cell) {
+ if (marked_vertices.size() == 0)
+ return true;
+
+ if (cell != mesh.active_cell_iterators().end())
+ for (unsigned int i = 0; i < cell->n_vertices(); ++i)
+ if (marked_vertices[cell->vertex_index(i)])
+ return true;
+
+ return false;
+ };
+
+ // check whether any cell in collection is marked
+ const auto any_cell_marked = [&cell_marked](const auto &cells) {
+ return std::any_of(cells.begin(),
+ cells.end(),
+ [&cell_marked](const auto &cell) {
+ return cell_marked(cell);
+ });
+ };
+
while (found_cell == false)
{
// First look at the vertices of the cell cell_hint. If it's an
// invalid cell, then query for the closest global vertex
- if (current_cell.state() == IteratorState::valid)
+ if (current_cell.state() == IteratorState::valid &&
+ cell_marked(cell_hint))
{
const auto cell_vertices = mapping.get_vertices(current_cell);
const unsigned int closest_vertex =
boost::geometry::index::satisfies(marked),
std::back_inserter(res));
- // We should have one and only one result
- AssertDimension(res.size(), 1);
- closest_vertex_index = res[0].second;
+ // Searching for a point which is located outside the
+ // triangulation results in res.size() = 0
+ Assert(res.size() < 2,
+ dealii::ExcMessage("There can not be multiple results"));
+
+ if (res.size() > 0)
+ if (any_cell_marked(vertex_to_cells[res[0].second]))
+ closest_vertex_index = res[0].second;
}
else
{
vertex_to_point = p - mesh.get_vertices()[closest_vertex_index];
}
+#ifdef DEBUG
+ {
+ // Double-check if found index is at marked cell
+ Assert(any_cell_marked(vertex_to_cells[closest_vertex_index]),
+ dealii::ExcMessage("Found non-marked vertex"));
+ }
+#endif
+
const double vertex_point_norm = vertex_to_point.norm();
if (vertex_point_norm > 0)
vertex_to_point /= vertex_point_norm;
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2001 - 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.
+//
+// ---------------------------------------------------------------------
+
+//
+// check find_active_cell_around_point finds only marked cells.
+// If there are no marked cell containing the point it should not find
+// anything. Latter is especially important if run with MPI.
+//
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_tools_cache.h>
+#include <deal.II/grid/tria.h>
+
+#include "../tests.h"
+
+/*
+ * Generate a non-matching grid consisting of two elements which
+ * do not share any face
+ */
+
+void
+generate_grid(Triangulation<2> &triangulation)
+{
+ // +-------------------+ +-------------------+
+ // | | | | |
+ // | 2 | 3 | | |
+ // | | | | |
+ // 10|_________|_________|8 9| 4 |
+ // | | | | |
+ // | | | | |
+ // | 0 | 1 | p | |
+ // | | | | |
+ // +-------------------+ +-------------------+
+
+ Triangulation<2> triangulation0;
+ Triangulation<2> triangulation1;
+ GridGenerator::hyper_rectangle(triangulation0, {-.5, -.5}, {.5, .5}, true);
+ for (const auto &face : triangulation0.active_face_iterators())
+ if (face->at_boundary() && face->boundary_id() == 1)
+ face->set_boundary_id(8);
+ for (const auto &face : triangulation0.active_face_iterators())
+ if (face->at_boundary() && face->boundary_id() == 0)
+ face->set_boundary_id(10);
+ triangulation0.refine_global(1);
+ GridGenerator::flatten_triangulation(triangulation0, triangulation1);
+
+ Triangulation<2> triangulation2;
+ GridGenerator::hyper_rectangle(triangulation2, {.5, -.5}, {1.5, .5}, true);
+ for (const auto &face : triangulation2.active_face_iterators())
+ if (face->at_boundary() && face->boundary_id() == 0)
+ face->set_boundary_id(9);
+
+ // tolerance 0. to ensure vertices are not merged
+ GridGenerator::merge_triangulations(
+ triangulation1, triangulation2, triangulation, 0., false);
+
+ // make sure boundary ids are kept
+ std::vector<types::boundary_id> boundary_ids;
+ for (const auto &face : triangulation1.active_face_iterators())
+ if (face->at_boundary())
+ boundary_ids.emplace_back(face->boundary_id());
+ for (const auto &face : triangulation2.active_face_iterators())
+ if (face->at_boundary())
+ boundary_ids.emplace_back(face->boundary_id());
+
+ unsigned int i = 0;
+ for (const auto &face : triangulation.active_face_iterators())
+ if (face->at_boundary())
+ face->set_boundary_id(boundary_ids[i++]);
+}
+
+std::vector<unsigned int>
+mark_vertices_at_boundary(const types::boundary_id boundary_id,
+ const Triangulation<2> & triangulation,
+ std::vector<bool> & marked_vertices)
+{
+ std::vector<unsigned int> marked_cell_idxs;
+
+ // mark_vertices_of_cells_at_boundary
+ for (const auto &cell : triangulation.active_cell_iterators())
+ {
+ for (unsigned int face = 0; face < cell->n_faces(); ++face)
+ {
+ if (cell->face(face)->at_boundary() &&
+ cell->face(face)->boundary_id() == boundary_id)
+ {
+ marked_cell_idxs.push_back(cell->index());
+ for (unsigned int v = 0; v < cell->face(face)->n_vertices(); ++v)
+ {
+ marked_vertices[cell->face(face)->vertex_index(v)] = true;
+ }
+ continue;
+ }
+ }
+ }
+ return marked_cell_idxs;
+}
+
+int
+find_cell_at_point(const Point<2> & p,
+ const Triangulation<2> & triangulation,
+ const std::vector<bool> &marked_vertices,
+ const Triangulation<2>::active_cell_iterator &cell_hint,
+ const double tolerance = 1e-6)
+{
+ const MappingQ<2> mapping{1};
+ const GridTools::Cache<2, 2> cache{triangulation, mapping};
+ const auto cell_and_pnt = GridTools::find_active_cell_around_point(
+ cache.get_mapping(),
+ cache.get_triangulation(),
+ p,
+ cache.get_vertex_to_cell_map(),
+ cache.get_vertex_to_cell_centers_directions(),
+ cell_hint,
+ marked_vertices,
+ cache.get_used_vertices_rtree(),
+ tolerance,
+ &cache.get_locally_owned_cell_bounding_boxes_rtree());
+
+ if (cell_and_pnt.first != triangulation.active_cell_iterators().end())
+ return cell_and_pnt.first->index();
+ return -1;
+}
+
+int
+main()
+{
+ // Setup
+ Triangulation<2> triangulation;
+ generate_grid(triangulation);
+
+ const Point<2> p{.5, -.25}; // point located between cell 4 and 1
+ const auto cell_hint = ++triangulation.begin_active(); // hint cell 1
+
+ initlog();
+ deallog << std::setprecision(4);
+
+ // Test1
+ deallog << "Test if cell not found:\n\n";
+ deallog << "Searching for point with marked cells that do not hold the "
+ "point:\n";
+ {
+ std::vector<bool> marked_vertices(triangulation.n_vertices(), false);
+ const auto cell_idxs =
+ mark_vertices_at_boundary(10, triangulation, marked_vertices);
+ const int found_cell =
+ find_cell_at_point(p, triangulation, marked_vertices, cell_hint);
+
+ deallog << "Marked cells ";
+ for (const auto &cell_idx : cell_idxs)
+ deallog << cell_idx << " ";
+ deallog << "\nFound cell: " << found_cell << "\n";
+ deallog << "\n---------------\n\n";
+ }
+
+ // Test2
+ deallog << "Test if marked cell is found with cell hint that has the "
+ "point but is not marked:\n\n";
+ {
+ std::vector<bool> marked_vertices(triangulation.n_vertices(), false);
+
+ // mark cell 4
+ const auto cell_idxs =
+ mark_vertices_at_boundary(9, triangulation, marked_vertices);
+ deallog << "Searching for cells ";
+ for (const auto &cell_idx : cell_idxs)
+ deallog << cell_idx << " ";
+ deallog << "with cell_hint " << cell_hint->index() << ":\n";
+
+ const int found_cell =
+ find_cell_at_point(p, triangulation, marked_vertices, cell_hint);
+ deallog << "Found cell: " << found_cell << "\n";
+ }
+ deallog << std::endl << std::flush;
+
+ return 0;
+}
projects / dealii.git / commitdiff