* they belong are the best case. Pre-sorting points, trying to minimize
* distances between them, might make the function extremely faster.
*
+ * @note If a point is not found inside the mesh, or is lying inside an
+ * artificial cell of a parallel::Triangulation, an exception is thrown.
+ *
* @note The actual return type of this function, i.e., the type referenced
* above as @p return_type, is
* @code
&cell_hint =
typename Triangulation<dim, spacedim>::active_cell_iterator());
+ /**
+ * This function is similar to GridTools::compute_point_locations(),
+ * but it tries to find and transform every point of @p points.
+ *
+ * @return A tuple containing four elements; the first three
+ * are documented in GridTools::compute_point_locations().
+ * The last element of the @p return_type contains the
+ * indices of points which are not not found inside the mesh
+ * or lie on artificial cells.
+ *
+ * @code
+ * std::tuple<
+ * std::vector<typename Triangulation<dim,
+ * spacedim>::active_cell_iterator>, std::vector<std::vector<Point<dim>>>,
+ * std::vector<std::vector<unsigned int>>,
+ * std::vector<unsigned int>>
+ * @endcode
+ *
+ * For a more detailed documentation see
+ * GridTools::compute_point_locations().
+ */
+ template <int dim, int spacedim>
+# ifndef DOXYGEN
+ std::tuple<
+ std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
+ std::vector<std::vector<Point<dim>>>,
+ std::vector<std::vector<unsigned int>>,
+ std::vector<unsigned int>>
+# else
+ return_type
+# endif
+ compute_point_locations_try_all(
+ const Cache<dim, spacedim> & cache,
+ const std::vector<Point<spacedim>> &points,
+ const typename Triangulation<dim, spacedim>::active_cell_iterator
+ &cell_hint =
+ typename Triangulation<dim, spacedim>::active_cell_iterator());
+
/**
* Given a @p cache and a list of
* @p local_points for each process, find the points lying on the locally
const std::vector<Point<spacedim>> &points,
const typename Triangulation<dim, spacedim>::active_cell_iterator
&cell_hint)
+ {
+ const auto cqmp = compute_point_locations_try_all(cache, points, cell_hint);
+ // Splitting the tuple's components
+ auto &cells = std::get<0>(cqmp);
+ auto &qpoints = std::get<1>(cqmp);
+ auto &maps = std::get<2>(cqmp);
+ auto &missing_points = std::get<3>(cqmp);
+ // If a point was not found, throwing an error, as the old
+ // implementation of compute_point_locations would have done
+ AssertThrow(std::get<3>(cqmp).size() == 0,
+ ExcPointNotFound<spacedim>(points[missing_points[0]]));
+
+ (void)missing_points;
+
+ return {std::move(cells), std::move(qpoints), std::move(maps)};
+ }
+
+
+
+ template <int dim, int spacedim>
+#ifndef DOXYGEN
+ std::tuple<
+ std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
+ std::vector<std::vector<Point<dim>>>,
+ std::vector<std::vector<unsigned int>>,
+ std::vector<unsigned int>>
+#else
+ return_type
+#endif
+ compute_point_locations_try_all(
+ const Cache<dim, spacedim> & cache,
+ const std::vector<Point<spacedim>> &points,
+ const typename Triangulation<dim, spacedim>::active_cell_iterator
+ &cell_hint)
{
// How many points are here?
- const unsigned int np = points.size();
+ const unsigned int np = points.size();
+ std::vector<unsigned int> points_outside;
std::tuple<
std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
std::vector<std::vector<Point<dim>>>,
- std::vector<std::vector<unsigned int>>>
+ std::vector<std::vector<unsigned int>>,
+ std::vector<unsigned int>>
cell_qpoint_map;
// Now the easy case.
std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
Point<dim>>
my_pair;
+
+ bool found = false;
+ unsigned int points_checked = 0;
if (cell_hint.state() == IteratorState::valid)
- my_pair =
- GridTools::find_active_cell_around_point(cache, points[0], cell_hint);
- else
- my_pair = GridTools::find_active_cell_around_point(cache, points[0]);
+ {
+ try
+ {
+ my_pair = GridTools::find_active_cell_around_point(cache,
+ points[0],
+ cell_hint);
+ found = true;
+ }
+ catch (const GridTools::ExcPointNotFound<dim> &)
+ {
+ points_outside.emplace_back(0);
+ }
+ ++points_checked;
+ }
+
- std::get<0>(cell_qpoint_map).emplace_back(my_pair.first);
- std::get<1>(cell_qpoint_map).emplace_back(1, my_pair.second);
- std::get<2>(cell_qpoint_map).emplace_back(1, 0);
+ while (!found && points_checked < np)
+ {
+ try
+ {
+ my_pair =
+ GridTools::find_active_cell_around_point(cache,
+ points[points_checked]);
+ found = true;
+ }
+ catch (const GridTools::ExcPointNotFound<dim> &)
+ {
+ points_outside.emplace_back(points_checked);
+ }
+ // Updating the position of the analyzed points
+ ++points_checked;
+ }
+
+ // If the point has been found in a cell, adding it
+ if (found)
+ {
+ std::get<0>(cell_qpoint_map).emplace_back(my_pair.first);
+ std::get<1>(cell_qpoint_map).emplace_back(1, my_pair.second);
+ std::get<2>(cell_qpoint_map).emplace_back(1, points_checked - 1);
+ }
// Now the second easy case.
- if (np == 1)
+ if (np == points_outside.size())
return cell_qpoint_map;
// Computing the cell center and diameter
Point<spacedim> cell_center = std::get<0>(cell_qpoint_map)[0]->center();
(0.5 + std::numeric_limits<double>::epsilon());
// Cycle over all points left
- for (unsigned int p = 1; p < np; ++p)
+ for (unsigned int p = points_checked; 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], std::get<0>(cell_qpoint_map).back());
- else
- my_pair = GridTools::find_active_cell_around_point(cache, points[p]);
+ try
+ {
+ if (cell_center.distance(points[p]) < cell_diameter)
+ my_pair = GridTools::find_active_cell_around_point(
+ cache, points[p], std::get<0>(cell_qpoint_map).back());
+ else
+ my_pair =
+ GridTools::find_active_cell_around_point(cache, points[p]);
+ }
+ catch (const GridTools::ExcPointNotFound<dim> &)
+ {
+ points_outside.push_back(p);
+ continue;
+ }
// Assuming the cell is probably the last cell added
if (my_pair.first == std::get<0>(cell_qpoint_map).back())
// The number of points in all
// the cells must be the same as
// the number of points we
- // started off from.
+ // started off from,
+ // plus the points which were ignored
for (unsigned int n = 0; n < c; ++n)
{
Assert(std::get<1>(cell_qpoint_map)[n].size() ==
std::get<2>(cell_qpoint_map)[n].size()));
qps += std::get<1>(cell_qpoint_map)[n].size();
}
- Assert(qps == np, ExcDimensionMismatch(qps, np));
+
+ Assert(qps + points_outside.size() == np,
+ ExcDimensionMismatch(qps + points_outside.size(), np));
#endif
+ std::get<3>(cell_qpoint_map) = points_outside;
return cell_qpoint_map;
}
deal_II_space_dimension>::active_cell_iterator &,
const std::vector<bool> &);
+ template std::tuple<std::vector<typename Triangulation<
+ deal_II_dimension,
+ deal_II_space_dimension>::active_cell_iterator>,
+ std::vector<std::vector<Point<deal_II_dimension>>>,
+ std::vector<std::vector<unsigned int>>,
+ std::vector<unsigned int>>
+ compute_point_locations_try_all(
+ 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>,
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 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::compute_point_locations_try_all
+// Test for corner case: vector of points, some of which
+// lie outside the mesh
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q.h>
+
+#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 <deal.II/numerics/fe_field_function.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+void
+test_compute_pt_loc(unsigned int n_points)
+{
+ 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]
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(std::max(6 - dim, 2));
+
+ // Creating the finite elements needed:
+ FE_Q<dim> fe(1);
+ DoFHandler<dim> dof_handler(tria);
+ dof_handler.distribute_dofs(fe);
+
+ // Creating the random points
+ std::vector<Point<dim>> points;
+
+ for (size_t i = 0; i < n_points; ++i)
+ {
+ points.push_back(random_point<dim>());
+ // One every three points is
+ // shifted outside the mesh
+ if (points.size() % 3 == dim - 1)
+ {
+ points.back()[0] += 5.0;
+ }
+ }
+
+ // Initializing the cache
+ GridTools::Cache<dim, dim> cache(tria);
+
+ auto cell_qpoint_map =
+ GridTools::compute_point_locations_try_all(cache, points);
+ const auto &cells = std::get<0>(cell_qpoint_map);
+ const auto &qpoints = std::get<1>(cell_qpoint_map);
+ const auto &maps = std::get<2>(cell_qpoint_map);
+ const auto &other_p = std::get<3>(cell_qpoint_map);
+ size_t n_cells = cells.size();
+
+ deallog << "Points found in " << n_cells << " cells" << std::endl;
+
+ unsigned int checked_points = 0;
+ // testing if the transformation is correct:
+ for (unsigned int i = 0; i < n_cells; ++i)
+ {
+ const auto &cell = cells[i];
+ const auto &quad = qpoints[i];
+ const auto &local_map = maps[i];
+
+ // Given the std::get<1>(cell_qpoint_map) of the current cell, compute the
+ // real points
+ FEValues<dim> fev(fe, quad, update_quadrature_points);
+ fev.reinit(cell);
+ const auto &real_quad = fev.get_quadrature_points();
+
+ for (unsigned int q = 0; q < real_quad.size(); ++q)
+ {
+ // Check if points are the same as real points
+ if (real_quad[q].distance(points[local_map[q]]) > 1e-10)
+ deallog << "Error on cell : " << cell << " at local point " << i
+ << ", corresponding to real point " << points[local_map[q]]
+ << ", that got transformed to " << real_quad[q]
+ << " instead." << std::endl;
+ else
+ ++checked_points;
+ }
+ }
+ for (auto i : other_p)
+ if (i % 3 != dim - 2)
+ {
+ deallog << "Error: point index " << i << " was not found!" << std::endl;
+ }
+
+ if (checked_points + other_p.size() == points.size())
+ deallog << "All points where processed correctly" << std::endl;
+ else
+ {
+ deallog << "Error: some points where not processed correctly"
+ << std::endl;
+ deallog << "Points found inside the mesh: " << checked_points
+ << std::endl;
+ deallog << "Discarded points: " << points.size() - checked_points
+ << std::endl;
+ deallog << "Points outside the mesh: " << other_p.size() << std::endl;
+ deallog << "Total points processed: " << points.size() << std::endl;
+ }
+ deallog << "Test finished" << std::endl;
+}
+
+int
+main()
+{
+ initlog();
+
+ deallog << "Deal.II compute_pt_loc_try_all:" << std::endl;
+ test_compute_pt_loc<2>(101);
+ test_compute_pt_loc<3>(202);
+}
--- /dev/null
+
+DEAL::Deal.II compute_pt_loc_try_all:
+DEAL::Testing for dim = 2
+DEAL::Testing on: 101 points.
+DEAL::Points found in 58 cells
+DEAL::All points where processed correctly
+DEAL::Test finished
+DEAL::Testing for dim = 3
+DEAL::Testing on: 202 points.
+DEAL::Points found in 119 cells
+DEAL::All points where processed correctly
+DEAL::Test finished
projects / dealii.git / commitdiff