From b5e72f697feed5bba376968cab862f9aba81b198 Mon Sep 17 00:00:00 2001 From: Giovanni Alzetta Date: Mon, 6 Nov 2017 11:55:06 +0000 Subject: [PATCH] Linked fe_field_function compute point locations to gridtools --- .../changes/minor/20171106GiovanniAlzetta | 4 + include/deal.II/grid/grid_tools.h | 4 +- include/deal.II/numerics/fe_field_function.h | 10 ++ .../numerics/fe_field_function.templates.h | 159 ++---------------- source/grid/grid_tools.cc | 14 +- source/grid/grid_tools.inst.in | 3 +- tests/grid/compute_point_locations_01.cc | 3 +- tests/grid/compute_point_locations_02.cc | 113 +++++++++++++ tests/grid/compute_point_locations_02.output | 10 ++ 9 files changed, 168 insertions(+), 152 deletions(-) create mode 100644 doc/news/changes/minor/20171106GiovanniAlzetta create mode 100644 tests/grid/compute_point_locations_02.cc create mode 100644 tests/grid/compute_point_locations_02.output diff --git a/doc/news/changes/minor/20171106GiovanniAlzetta b/doc/news/changes/minor/20171106GiovanniAlzetta new file mode 100644 index 0000000000..1803b0f7d7 --- /dev/null +++ b/doc/news/changes/minor/20171106GiovanniAlzetta @@ -0,0 +1,4 @@ +Modified: Function Functions::FEFieldFunction::compute_point_locations now calls GridTools::compute_point_locations +Modified: Function GridTools::compute_point_locations now accepts a cell hint (a test for this hint has been added) +
+(Giovanni Alzetta, 2017/11/06) diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index 8fe24838c1..40161dab24 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -657,7 +657,9 @@ namespace GridTools std::vector< std::vector< Point > >, std::vector< std::vector > > compute_point_locations(const Cache &cache, - const std::vector > &points); + const std::vector > &points, + const typename Triangulation::active_cell_iterator &cell_hint + = typename Triangulation::active_cell_iterator()); /** * Return a map of index:Point, containing the used vertices of the diff --git a/include/deal.II/numerics/fe_field_function.h b/include/deal.II/numerics/fe_field_function.h index 119a7d43b4..7e82109584 100644 --- a/include/deal.II/numerics/fe_field_function.h +++ b/include/deal.II/numerics/fe_field_function.h @@ -24,6 +24,7 @@ #include #include #include +#include #include @@ -413,6 +414,10 @@ namespace Functions * @return This function returns the number of cells that * collectively contain the set of points give as @p * points. This also equals the lengths of the output arrays. + * + * This function simply calls GridTools::compute_point_locations : + * using the original function avoids computing a + * new Cache at every function call. */ unsigned int compute_point_locations @@ -444,6 +449,11 @@ namespace Functions */ const Mapping &mapping; + /** + * The Cache object + */ + GridTools::Cache cache; + /** * The latest cell hint. */ diff --git a/include/deal.II/numerics/fe_field_function.templates.h b/include/deal.II/numerics/fe_field_function.templates.h index 7acf5f9261..e26d997d1e 100644 --- a/include/deal.II/numerics/fe_field_function.templates.h +++ b/include/deal.II/numerics/fe_field_function.templates.h @@ -20,6 +20,7 @@ #include #include #include +#include #include #include #include @@ -28,6 +29,9 @@ #include #include +#include + + DEAL_II_NAMESPACE_OPEN @@ -44,6 +48,7 @@ namespace Functions dh(&mydh, "FEFieldFunction"), data_vector(myv), mapping(mymapping), + cache(dh->get_triangulation(),mymapping), cell_hint(dh->end()) { } @@ -438,150 +443,16 @@ namespace Functions std::vector > > &qpoints, std::vector > &maps) const { - // How many points are here? - const unsigned int np = points.size(); - - // Reset output maps. - cells.clear(); - qpoints.clear(); - maps.clear(); - - // Now the easy case. - if (np==0) return 0; - - // Keep track of the points we found - std::vector point_flags(np, false); - - // Set this to true until all points have been classified - bool left_over = true; - - // Current quadrature point - typename DoFHandlerType::active_cell_iterator cell = cell_hint.get(); - if (cell == dh->end()) - cell = dh->begin_active(); - - { - // see if the point is inside the cell. there are two ways that - // transform_real_to_unit_cell can indicate that a point is - // outside: by returning coordinates that lie outside the - // reference cell, or by throwing an exception. handle both - boost::optional > - qp = get_reference_coordinates (cell, points[0]); - if (!qp) - { - const std::pair::type, Point > - my_pair = GridTools::find_active_cell_around_point - (mapping, *dh, points[0]); - AssertThrow (!my_pair.first->is_artificial(), - VectorTools::ExcPointNotAvailableHere()); - - cell = my_pair.first; - qp = my_pair.second; - point_flags[0] = true; - } - - // check that the cell is available: - AssertThrow (!cell->is_artificial(), - VectorTools::ExcPointNotAvailableHere()); - - // Put in the first point. - cells.push_back(cell); - qpoints.emplace_back(1, qp.get()); - maps.emplace_back(1, 0); - } - - - // Check if we need to do anything else - if (points.size() > 1) - left_over = true; - else - left_over = 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 > - qp = get_reference_coordinates (cells[c], points[p]); - if (qp) - { - point_flags[p] = true; - qpoints[c].push_back(qp.get()); - maps[c].push_back(p); - } - else - { - // 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::type, Point > my_pair - = GridTools::find_active_cell_around_point (mapping, *dh, points[first_outside]); - AssertThrow (!my_pair.first->is_artificial(), - VectorTools::ExcPointNotAvailableHere()); - - cells.push_back(my_pair.first); - qpoints.emplace_back(1, my_pair.second); - maps.emplace_back(1, first_outside); - c++; - point_flags[first_outside] = true; - // And check if we can exit the loop now - if (first_outside == np-1) - left_over = false; - } - } - - // Augment of one the number of cells - ++c; - // Debug Checking - Assert(c == cells.size(), ExcInternalError()); - - Assert(c == maps.size(), - ExcDimensionMismatch(c, maps.size())); - - Assert(c == qpoints.size(), - ExcDimensionMismatch(c, qpoints.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); + cells.resize(tria_cells.size()); + unsigned int i = 0; + for (const auto &c: tria_cells) + cells[i++] = typename DoFHandlerType::cell_iterator(*c,dh); + qpoints = std::get<1>(cell_qpoint_map); + maps = std::get<2>(cell_qpoint_map); + return cells.size(); } diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 4d9f8e9066..ecf5684ada 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -5176,8 +5176,9 @@ next_cell: std::vector::active_cell_iterator >, std::vector< std::vector< Point > >, std::vector< std::vector > > - compute_point_locations(const Cache &cache, - const std::vector > &points) + compute_point_locations(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(); @@ -5193,8 +5194,13 @@ next_cell: // We begin by finding the cell/transform of the first point std::pair::active_cell_iterator, Point > - my_pair = GridTools::find_active_cell_around_point - (cache, points[0]); + my_pair; + 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]); std::get<0>(cell_qpoint_map).emplace_back(my_pair.first); std::get<1>(cell_qpoint_map).emplace_back(1, my_pair.second); diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index 39cb48ebfe..51d4df9a32 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -156,7 +156,8 @@ for (deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS 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 std::vector< Point< deal_II_space_dimension > > &, + const typename Triangulation< deal_II_dimension, deal_II_space_dimension>::active_cell_iterator &); \} #endif diff --git a/tests/grid/compute_point_locations_01.cc b/tests/grid/compute_point_locations_01.cc index ab0bca3fe8..981142f059 100644 --- a/tests/grid/compute_point_locations_01.cc +++ b/tests/grid/compute_point_locations_01.cc @@ -61,8 +61,7 @@ void test_compute_pt_loc(unsigned int n_points) // Initializing the cache GridTools::Cache 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); + auto 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; diff --git a/tests/grid/compute_point_locations_02.cc b/tests/grid/compute_point_locations_02.cc new file mode 100644 index 0000000000..bbfdb40812 --- /dev/null +++ b/tests/grid/compute_point_locations_02.cc @@ -0,0 +1,113 @@ +// --------------------------------------------------------------------- +// +// 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); + // Finding the first cell and giving it as hint + auto my_pair = GridTools::find_active_cell_around_point + (cache, points[0]); + auto cell_qpoint_map = GridTools::compute_point_locations(cache,points,my_pair.first); + 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_02.output b/tests/grid/compute_point_locations_02.output new file mode 100644 index 0000000000..0a6e0e1554 --- /dev/null +++ b/tests/grid/compute_point_locations_02.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 -- 2.39.5