From: Giovanni Alzetta Date: Wed, 9 Jan 2019 14:32:24 +0000 (+0100) Subject: Added new function compute point locations try all X-Git-Tag: v9.1.0-rc1~331^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b2bac494a46b5322fa7d2507c937b1b200bf5491;p=dealii.git Added new function compute point locations try all --- diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index 3a89f434b7..809de1bcca 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -720,6 +720,9 @@ namespace GridTools * 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 @@ -750,6 +753,44 @@ namespace GridTools &cell_hint = typename Triangulation::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::active_cell_iterator>, std::vector>>, + * std::vector>, + * std::vector> + * @endcode + * + * For a more detailed documentation see + * GridTools::compute_point_locations(). + */ + template +# ifndef DOXYGEN + std::tuple< + std::vector::active_cell_iterator>, + std::vector>>, + std::vector>, + std::vector> +# else + return_type +# endif + compute_point_locations_try_all( + const Cache & cache, + const std::vector> &points, + 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 diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 201dfe03f9..bd97f63dc4 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -4193,14 +4193,50 @@ namespace GridTools const std::vector> &points, const typename Triangulation::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(points[missing_points[0]])); + + (void)missing_points; + + return {std::move(cells), std::move(qpoints), std::move(maps)}; + } + + + + template +#ifndef DOXYGEN + std::tuple< + std::vector::active_cell_iterator>, + std::vector>>, + std::vector>, + std::vector> +#else + return_type +#endif + compute_point_locations_try_all( + const Cache & cache, + const std::vector> &points, + const typename Triangulation::active_cell_iterator + &cell_hint) { // How many points are here? - const unsigned int np = points.size(); + const unsigned int np = points.size(); + std::vector points_outside; std::tuple< std::vector::active_cell_iterator>, std::vector>>, - std::vector>> + std::vector>, + std::vector> cell_qpoint_map; // Now the easy case. @@ -4211,18 +4247,53 @@ namespace GridTools std::pair::active_cell_iterator, Point> 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 &) + { + 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 &) + { + 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 cell_center = std::get<0>(cell_qpoint_map)[0]->center(); @@ -4230,15 +4301,24 @@ namespace GridTools (0.5 + std::numeric_limits::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 &) + { + 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()) @@ -4298,7 +4378,8 @@ namespace GridTools // 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() == @@ -4307,9 +4388,12 @@ namespace GridTools 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; } diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index c8b0ce9120..23dd04b907 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -93,6 +93,19 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) deal_II_space_dimension>::active_cell_iterator &, const std::vector &); + template std::tuple::active_cell_iterator>, + std::vector>>, + std::vector>, + std::vector> + compute_point_locations_try_all( + const Cache &, + const std::vector> &, + const typename Triangulation< + deal_II_dimension, + deal_II_space_dimension>::active_cell_iterator &); + template std::tuple::active_cell_iterator>, diff --git a/tests/grid/compute_point_locations_try_all_01.cc b/tests/grid/compute_point_locations_try_all_01.cc new file mode 100644 index 0000000000..ddb72f2e3f --- /dev/null +++ b/tests/grid/compute_point_locations_try_all_01.cc @@ -0,0 +1,143 @@ +// --------------------------------------------------------------------- +// +// 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 + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +template +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 tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(std::max(6 - dim, 2)); + + // Creating the finite elements needed: + FE_Q fe(1); + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + // Creating the random points + std::vector> points; + + for (size_t i = 0; i < n_points; ++i) + { + points.push_back(random_point()); + // 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 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 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); +} diff --git a/tests/grid/compute_point_locations_try_all_01.output b/tests/grid/compute_point_locations_try_all_01.output new file mode 100644 index 0000000000..1272686f85 --- /dev/null +++ b/tests/grid/compute_point_locations_try_all_01.output @@ -0,0 +1,12 @@ + +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