From: Giovanni Alzetta Date: Tue, 24 Oct 2017 18:14:57 +0000 (+0000) Subject: Added new function compute point locations to GridTools X-Git-Tag: v9.0.0-rc1~853^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F5314%2Fhead;p=dealii.git Added new function compute point locations to GridTools --- diff --git a/doc/news/changes/minor/20171024GiovanniAlzetta b/doc/news/changes/minor/20171024GiovanniAlzetta new file mode 100644 index 0000000000..6718a083fd --- /dev/null +++ b/doc/news/changes/minor/20171024GiovanniAlzetta @@ -0,0 +1,3 @@ +Added: Function GridTools::compute_point_locations and a test for it. The function uses the recently added Cache to improve performance. +
+(Giovanni Alzetta, 2017/10/24) diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index 8b27c9f9b4..ca8f434d9d 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -605,6 +605,41 @@ namespace GridTools */ /*@{*/ + /** + * Given a @p cache and a list of @p points create the quadrature rules. This + * function returns a tuple containing the following elements: + * - The first ( get<0>), call it cells, is a vector of a vector cells of the all cells + * that contain at least one of the @p points + * - The second a vector qpoints of vector of points, containing in qpoints[i] + * the reference positions of all points that fall within the cell cells[i] + * - The third a vector indices of vector of integers, containing the mapping between + * local numbering in qpoints, and global index in points + * + * If points[a] and points[b] are the only two points that fall in cells[c], then + * qpoints[c][0] and qpoints[c][1] are the reference positions of points[a] and points[b] + * in cells[c], and indices[c][0] = a, indices[c][1] = b. The function + * Mapping::tansform_unit_to_real(qpoints[c][0]) returns points[a]. + * + * The algorithm assumes it's easier to look for a point in the cell that was used previously. + * For this reason random points are, computationally speaking, the worst case scenario while + * points grouped by the cell to which they belong are the best case. + * Pre-sorting points, trying to minimize distances, can make the algorithm much faster. + * + * Notice: given the center of a cell, points at distance greater then + * @p distance_factor * cell.diameter() are considered outside the cell. In some cases, such + * as curved meshes, this leads to cell's repetitions in the output tuple. + * To avoid this either enlarge @p distance_factor (depending on the regularity of the + * triangulation) and/or merge the repetitions contained in the output. + */ + template + std::tuple< + std::vector::active_cell_iterator >, + std::vector< std::vector< Point > >, + std::vector< std::vector > > + compute_point_locations(const Cache &cache, + const std::vector > &points, + const double &distance_factor=0.5); + /** * Return a map of index:Point, containing the used vertices of the * given `container`. The key of the returned map is the global index in the diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 196abb7329..6dd0b32123 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -4915,6 +4915,177 @@ next_cell: + template + std::tuple< + std::vector::active_cell_iterator >, + std::vector< std::vector< Point > >, + std::vector< std::vector > > + compute_point_locations(const Cache &cache, + const std::vector > &points, + const double &distance_factor) + { + // How many points are here? + const unsigned int np = points.size(); + + std::tuple< + std::vector::active_cell_iterator >, + std::vector< std::vector< Point > >, + std::vector< std::vector > > + cell_qpoint_map; + + // Now the easy case. + if (np==0) return cell_qpoint_map; + + // If distance_factor is too small we might avoid points which are inside the cell + Assert(distance_factor >= 0.5, + ExcMessage("distance_factor value must be >= 0.5")); + + // Keep track of the points we + // found + std::vector point_flags(np, false); + + // Classify the first point: + // find active cell returns a pair (cell,transformed point) + const std::pair::active_cell_iterator, + Point > + my_pair = GridTools::find_active_cell_around_point + (cache, points[0]); + + // Add data about the first point + std::get<0>(cell_qpoint_map).push_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); + + // Case only one point is now done + if (np==1) + return cell_qpoint_map; + + // Compute the information about the current cell + Point cell_center = my_pair.first->center(); + double cell_diameter = my_pair.first->diameter()* + (distance_factor + std::numeric_limits::epsilon() ); + + // Set this to true until all + // points have been classified + bool left_over = true; + + // Flag to signal that a point is not in the current cell + bool flag_outside = false; + + // This is the first index of a non processed point + unsigned int first_outside = 1; + + // And this is the index of the current cell + unsigned int c = 0; + + while (left_over == true) + { + // Assume this is the last one + left_over = false; + Assert(first_outside < np, + ExcIndexRange(first_outside, 0, np)); + + // If we found one in this cell, keep looking in the same cell + for (unsigned int p=first_outside; p cell_diameter ) + flag_outside = true; + else + { + //The point is close to the cell: try transforming it + try + { + Point qp =cache.get_mapping().transform_real_to_unit_cell + (std::get<0>(cell_qpoint_map)[c], points[p]); + if (GeometryInfo::is_inside_unit_cell(qp)) + { + point_flags[p] = true; + std::get<1>(cell_qpoint_map)[c].push_back(qp); + std::get<2>(cell_qpoint_map)[c].push_back(p); + } + else + // The point is outside the cell + flag_outside = true; + } + catch (const typename Mapping::ExcTransformationFailed &) + { + // transformation failed: assume the point is outside + flag_outside = true; + } + } + + if (flag_outside) + { + // Set things up for next round + if (left_over == false) + first_outside = p; + left_over = true; + } + } + // If we got here and there is no left over, we are + // done. Else we need to find the next cell + if (left_over == true) + { + const std::pair< + typename Triangulation::active_cell_iterator, Point > + my_pair + = GridTools::find_active_cell_around_point (cache, points[first_outside]); + + std::get<0>(cell_qpoint_map).push_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, first_outside); + // Check if we can exit the loop now because only the last point was missing + if (first_outside == np-1) + left_over = false; + + // Update the information about the current cell + cell_center = my_pair.first->center(); + cell_diameter = my_pair.first->diameter()* + (distance_factor + std::numeric_limits::epsilon() ); + c++; + point_flags[first_outside] = true; + } + } + + // Augment of one the number of cells + ++c; + // Debug Checking + Assert(c == std::get<0>(cell_qpoint_map).size(), ExcInternalError()); + + Assert(c == std::get<2>(cell_qpoint_map).size(), + ExcDimensionMismatch(c, std::get<2>(cell_qpoint_map).size())); + + Assert(c == std::get<1>(cell_qpoint_map).size(), + ExcDimensionMismatch(c, std::get<1>(cell_qpoint_map).size())); + +#ifdef DEBUG + unsigned int qps = 0; + // The number of points in all + // the cells must be the same as + // the number of points we + // started off from. + for (unsigned int n=0; n(cell_qpoint_map)[n].size() == + std::get<2>(cell_qpoint_map)[n].size(), + ExcDimensionMismatch(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)); +#endif + + return cell_qpoint_map; + } + + + template std::map > extract_used_vertices(const Triangulation &container, diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index 27020ca365..c4de24d473 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -139,8 +139,13 @@ for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS const typename Triangulation::active_cell_iterator &, const std::vector &); - - \} + 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 > > > + compute_point_locations(const Cache< deal_II_dimension, deal_II_space_dimension > &, + const std::vector< Point< deal_II_space_dimension > > &, + const double &); + \} #endif } diff --git a/tests/grid/compute_point_locations_01.cc b/tests/grid/compute_point_locations_01.cc new file mode 100644 index 0000000000..334a19a2e7 --- /dev/null +++ b/tests/grid/compute_point_locations_01.cc @@ -0,0 +1,112 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// Test Functions::FEFieldFunction::compute_point_locations + +#include "../tests.h" +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include + +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 p; + for (unsigned int d=0; d cache(tria); + + std::tuple< std::vector::active_cell_iterator >, std::vector< std::vector< Point > >, std::vector > > + cell_qpoint_map = GridTools::compute_point_locations(cache,points); + size_t n_cells = std::get<0>(cell_qpoint_map).size(); + + deallog << "Points found in " << n_cells << " cells" << std::endl; + + // testing if the transformation is correct: + // For each cell check if the quadrature points in the i-th FE + // are the same as std::get<2>(cell_qpoint_map)[i] + for (unsigned int i=0; i(cell_qpoint_map).size(); ++i) + { + auto &cell = std::get<0>(cell_qpoint_map)[i]; + auto &quad = std::get<1>(cell_qpoint_map)[i]; + auto &local_map = std::get<2>(cell_qpoint_map)[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 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; + } + } + deallog << "Test finished" << std::endl; +} + +int main() +{ + initlog(); + + deallog << "Deal.II compute_pt_loc:" << std::endl; + test_compute_pt_loc<2>(100); + test_compute_pt_loc<3>(200); +} diff --git a/tests/grid/compute_point_locations_01.output b/tests/grid/compute_point_locations_01.output new file mode 100644 index 0000000000..0a6e0e1554 --- /dev/null +++ b/tests/grid/compute_point_locations_01.output @@ -0,0 +1,10 @@ + +DEAL::Deal.II compute_pt_loc: +DEAL::Testing for dim = 2 +DEAL::Testing on: 100 points. +DEAL::Points found in 81 cells +DEAL::Test finished +DEAL::Testing for dim = 3 +DEAL::Testing on: 200 points. +DEAL::Points found in 170 cells +DEAL::Test finished