--- /dev/null
+Improved: GridTools::find_active_cell_around_point() and
+GridTools::find_all_active_cells_around_point() now also work
+for simplices.
+<br>
+(Peter Munch, David Wells, 2023/07/25)
{
const Point<dim> p_unit =
mapping.transform_real_to_unit_cell(*cell, p);
- if (GeometryInfo<dim>::is_inside_unit_cell(p_unit,
- tolerance))
+ if ((*cell)->reference_cell().contains_point(p_unit,
+ tolerance))
{
cell_and_position.first = *cell;
cell_and_position.second = p_unit;
// outside it is and whether we want to use this cell as a
// backup if we can't find a cell within which the point
// lies.
- const double dist =
- GeometryInfo<dim>::distance_to_unit_cell(p_unit);
+ const double dist = p_unit.distance(
+ (*cell)->reference_cell().closest_point(p_unit));
if (dist < best_distance)
{
best_distance = dist;
// insert the fist cell and point into the vector
cells_and_points.push_back(first_cell);
- // check if the given point is on the surface of the unit cell. If yes,
- // need to find all neighbors
const Point<dim> 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)
+ std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
+ cells_to_add;
+
+ if (my_cell->reference_cell().is_hyper_cube())
{
- distance_to_center[d] = std::abs(unit_point[d] - 0.5);
- if (distance_to_center[d] > 0.5 - tolerance)
+ // check if the given point is on the surface of the unit cell. If yes,
+ // need to find all neighbors
+
+ 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)
{
- ++n_dirs_at_threshold;
- last_point_at_threshold = d;
+ distance_to_center[d] = std::abs(unit_point[d] - 0.5);
+ if (distance_to_center[d] > 0.5 - tolerance)
+ {
+ ++n_dirs_at_threshold;
+ last_point_at_threshold = d;
+ }
}
- }
- std::vector<typename MeshType<dim, spacedim>::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))
+ // point is within face -> only need neighbor
+ if (n_dirs_at_threshold == 1)
{
- const auto neighbor_cell = my_cell->neighbor(neighbor_index);
+ 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))
+ {
+ const auto neighbor_cell = my_cell->neighbor(neighbor_index);
- if (neighbor_cell->is_active())
- cells_to_add.push_back(neighbor_cell);
- else
- for (const auto &child_cell : neighbor_cell->child_iterators())
- {
- if (child_cell->is_active())
- cells_to_add.push_back(child_cell);
- }
+ if (neighbor_cell->is_active())
+ cells_to_add.push_back(neighbor_cell);
+ else
+ for (const auto &child_cell :
+ neighbor_cell->child_iterators())
+ {
+ if (child_cell->is_active())
+ cells_to_add.push_back(child_cell);
+ }
+ }
}
- }
- // 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;
+ // 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;
- const auto fu = [&](const auto &tentative_cells) {
- for (const auto &cell : tentative_cells)
- if (cell != my_cell)
- cells_to_add.push_back(cell);
- };
+ const auto fu = [&](const auto &tentative_cells) {
+ for (const auto &cell : tentative_cells)
+ if (cell != my_cell)
+ cells_to_add.push_back(cell);
+ };
- const auto vertex_index = my_cell->vertex_index(local_vertex_index);
+ const auto vertex_index = my_cell->vertex_index(local_vertex_index);
- if (vertex_to_cells != nullptr)
- fu((*vertex_to_cells)[vertex_index]);
- else
- fu(find_cells_adjacent_to_vertex(mesh, vertex_index));
- }
- // point on line in 3d: We cannot simply take the intersection between
- // the two vertices of cells 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<unsigned int, unsigned int> 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 (vertex_to_cells != nullptr)
+ fu((*vertex_to_cells)[vertex_index]);
+ else
+ fu(find_cells_adjacent_to_vertex(mesh, vertex_index));
+ }
+ // point on line in 3d: We cannot simply take the intersection between
+ // the two vertices of cells 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)
{
- if (distance_to_center[d] > 0.5 - tolerance)
+ std::pair<unsigned int, unsigned int> 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)
{
- vertex_indices[count_vertex_indices].first = d;
- vertex_indices[count_vertex_indices].second =
- unit_point[d] > 0.5 ? 1 : 0;
- count_vertex_indices++;
+ if (distance_to_center[d] > 0.5 - tolerance)
+ {
+ 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;
}
- else
- free_direction = d;
- }
- AssertDimension(count_vertex_indices, 2);
- Assert(free_direction != numbers::invalid_unsigned_int,
- ExcInternalError());
+ 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)
+ {
+ const auto fu = [&](const auto &tentative_cells) {
+ for (const auto &cell : tentative_cells)
+ {
+ bool cell_not_yet_present = true;
+ for (const 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);
+ }
+ };
+
+ const auto vertex_index =
+ my_cell->vertex_index(first_vertex + (d << free_direction));
+
+ if (vertex_to_cells != nullptr)
+ fu((*vertex_to_cells)[vertex_index]);
+ else
+ fu(find_cells_adjacent_to_vertex(mesh, vertex_index));
+ }
+ }
+ }
+ else
+ {
+ // Note: The non-hypercube path takes a very naive approach and
+ // checks all possible neighbors. This can be made faster by 1)
+ // checking if the point is in the inner cell and 2) identifying
+ // the right lines/vertices so that the number of potential
+ // neighbors is reduced.
- 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)
+ for (const auto v : my_cell->vertex_indices())
{
const auto fu = [&](const auto &tentative_cells) {
for (const auto &cell : tentative_cells)
}
};
- const auto vertex_index =
- my_cell->vertex_index(first_vertex + (d << free_direction));
+ const auto vertex_index = my_cell->vertex_index(v);
if (vertex_to_cells != nullptr)
fu((*vertex_to_cells)[vertex_index]);
}
}
- const double original_distance_to_unit_cell =
- my_cell->reference_cell().closest_point(unit_point).distance(unit_point);
for (const auto &cell : cells_to_add)
{
if (cell != my_cell)
{
const Point<dim> p_unit =
mapping.transform_real_to_unit_cell(cell, p);
- if (cell->reference_cell().closest_point(p_unit).distance(
- p_unit) < original_distance_to_unit_cell + tolerance)
+ if (cell->reference_cell().contains_point(p_unit, tolerance))
cells_and_points.emplace_back(cell, p_unit);
}
catch (typename Mapping<dim>::ExcTransformationFailed &)
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2023 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.
+//
+// ---------------------------------------------------------------------
+
+// Test GridTools::find_active_cell_around_point() and
+// GridTools::find_all_active_cells_around_point for simplices. These
+// functions are used in find_all_locally_owned_active_cells_around_point(),
+// which is used by distributed_compute_point_locations().
+
+#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"
+
+template <int dim>
+void
+test()
+{
+ const auto reference_cell = ReferenceCells::get_simplex<dim>();
+ const auto &mapping =
+ reference_cell.template get_default_linear_mapping<dim>();
+
+ Triangulation<dim> tria;
+ if (false)
+ GridGenerator::reference_cell(tria, reference_cell);
+ else
+ GridGenerator::subdivided_hyper_cube_with_simplices(tria, 2);
+
+ GridTools::Cache<dim> cache(tria, mapping);
+
+ const unsigned n_subdivisions = 8;
+
+ for (unsigned int i = 0; i <= n_subdivisions; ++i)
+ for (unsigned int j = 0; j <= n_subdivisions; ++j)
+ {
+ Point<dim> p(1.0 * i / n_subdivisions, 1.0 * j / n_subdivisions);
+
+ const auto first_point =
+ GridTools::find_active_cell_around_point(mapping, tria, p, {}, 1e-6);
+
+ const auto result = GridTools::find_all_active_cells_around_point(
+ mapping, tria, p, 1e-6, {first_point.first, first_point.second});
+
+ deallog << result.size() << std::endl; // should be 2
+ for (const auto &cell : result)
+ deallog << cell.first->id() << " " << cell.second << std::endl;
+ deallog << std::endl;
+ }
+}
+
+main(int argc, char **argv)
+{
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+ MPILogInitAll all;
+
+ deallog.precision(8);
+
+ test<2>();
+}
--- /dev/null
+
+DEAL:0::1
+DEAL:0::0_0: 0.0000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::0_0: 0.0000000 0.25000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::0_0: 0.0000000 0.50000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::0_0: 0.0000000 0.75000000
+DEAL:0::
+DEAL:0::3
+DEAL:0::0_0: 0.0000000 1.0000000
+DEAL:0::1_0: 1.0000000 0.0000000
+DEAL:0::4_0: 0.0000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::4_0: 0.0000000 0.25000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::4_0: 0.0000000 0.50000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::4_0: 0.0000000 0.75000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::4_0: 0.0000000 1.0000000
+DEAL:0::5_0: 1.0000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::0_0: 0.25000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::0_0: 0.25000000 0.25000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::0_0: 0.25000000 0.50000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::0_0: 0.25000000 0.75000000
+DEAL:0::1_0: 0.75000000 0.25000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::1_0: 0.75000000 0.0000000
+DEAL:0::4_0: 0.25000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::4_0: 0.25000000 0.25000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::4_0: 0.25000000 0.50000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::4_0: 0.25000000 0.75000000
+DEAL:0::5_0: 0.75000000 0.25000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::5_0: 0.75000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::0_0: 0.50000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::0_0: 0.50000000 0.25000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::0_0: 0.50000000 0.50000000
+DEAL:0::1_0: 0.50000000 0.50000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::1_0: 0.50000000 0.25000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::1_0: 0.50000000 0.0000000
+DEAL:0::4_0: 0.50000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::4_0: 0.50000000 0.25000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::4_0: 0.50000000 0.50000000
+DEAL:0::5_0: 0.50000000 0.50000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::5_0: 0.50000000 0.25000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::5_0: 0.50000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::0_0: 0.75000000 0.0000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::0_0: 0.75000000 0.25000000
+DEAL:0::1_0: 0.25000000 0.75000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::1_0: 0.25000000 0.50000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::1_0: 0.25000000 0.25000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::1_0: 0.25000000 0.0000000
+DEAL:0::4_0: 0.75000000 0.0000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::4_0: 0.75000000 0.25000000
+DEAL:0::5_0: 0.25000000 0.75000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::5_0: 0.25000000 0.50000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::5_0: 0.25000000 0.25000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::5_0: 0.25000000 0.0000000
+DEAL:0::
+DEAL:0::3
+DEAL:0::0_0: 1.0000000 0.0000000
+DEAL:0::1_0: 0.0000000 1.0000000
+DEAL:0::2_0: 0.0000000 0.0000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::1_0: 0.0000000 0.75000000
+DEAL:0::2_0: 0.0000000 0.25000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::1_0: 0.0000000 0.50000000
+DEAL:0::2_0: 0.0000000 0.50000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::1_0: 0.0000000 0.25000000
+DEAL:0::2_0: 0.0000000 0.75000000
+DEAL:0::
+DEAL:0::6
+DEAL:0::1_0: 0.0000000 0.0000000
+DEAL:0::2_0: 0.0000000 1.0000000
+DEAL:0::3_0: 1.0000000 0.0000000
+DEAL:0::4_0: 1.0000000 0.0000000
+DEAL:0::5_0: 0.0000000 1.0000000
+DEAL:0::6_0: 0.0000000 0.0000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::5_0: 0.0000000 0.75000000
+DEAL:0::6_0: 0.0000000 0.25000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::5_0: 0.0000000 0.50000000
+DEAL:0::6_0: 0.0000000 0.50000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::5_0: 0.0000000 0.25000000
+DEAL:0::6_0: 0.0000000 0.75000000
+DEAL:0::
+DEAL:0::3
+DEAL:0::5_0: 0.0000000 0.0000000
+DEAL:0::6_0: 0.0000000 1.0000000
+DEAL:0::7_0: 1.0000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::2_0: 0.25000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::2_0: 0.25000000 0.25000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::2_0: 0.25000000 0.50000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::2_0: 0.25000000 0.75000000
+DEAL:0::3_0: 0.75000000 0.25000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::3_0: 0.75000000 0.0000000
+DEAL:0::6_0: 0.25000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::6_0: 0.25000000 0.25000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::6_0: 0.25000000 0.50000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::6_0: 0.25000000 0.75000000
+DEAL:0::7_0: 0.75000000 0.25000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::7_0: 0.75000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::2_0: 0.50000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::2_0: 0.50000000 0.25000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::2_0: 0.50000000 0.50000000
+DEAL:0::3_0: 0.50000000 0.50000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::3_0: 0.50000000 0.25000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::3_0: 0.50000000 0.0000000
+DEAL:0::6_0: 0.50000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::6_0: 0.50000000 0.25000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::6_0: 0.50000000 0.50000000
+DEAL:0::7_0: 0.50000000 0.50000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::7_0: 0.50000000 0.25000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::7_0: 0.50000000 0.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::2_0: 0.75000000 0.0000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::2_0: 0.75000000 0.25000000
+DEAL:0::3_0: 0.25000000 0.75000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::3_0: 0.25000000 0.50000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::3_0: 0.25000000 0.25000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::3_0: 0.25000000 0.0000000
+DEAL:0::6_0: 0.75000000 0.0000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::6_0: 0.75000000 0.25000000
+DEAL:0::7_0: 0.25000000 0.75000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::7_0: 0.25000000 0.50000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::7_0: 0.25000000 0.25000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::7_0: 0.25000000 0.0000000
+DEAL:0::
+DEAL:0::2
+DEAL:0::2_0: 1.0000000 0.0000000
+DEAL:0::3_0: 0.0000000 1.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::3_0: 0.0000000 0.75000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::3_0: 0.0000000 0.50000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::3_0: 0.0000000 0.25000000
+DEAL:0::
+DEAL:0::3
+DEAL:0::3_0: 0.0000000 0.0000000
+DEAL:0::6_0: 1.0000000 0.0000000
+DEAL:0::7_0: 0.0000000 1.0000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::7_0: 0.0000000 0.75000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::7_0: 0.0000000 0.50000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::7_0: 0.0000000 0.25000000
+DEAL:0::
+DEAL:0::1
+DEAL:0::7_0: 0.0000000 0.0000000
+DEAL:0::
projects / dealii.git / commitdiff