--- /dev/null
+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)
+<br>
+(Giovanni Alzetta, 2017/11/06)
std::vector< std::vector< Point<dim> > >,
std::vector< std::vector<unsigned int> > >
compute_point_locations(const Cache<dim,spacedim> &cache,
- const std::vector<Point<spacedim> > &points);
+ const std::vector<Point<spacedim> > &points,
+ const typename Triangulation<dim, spacedim>::active_cell_iterator &cell_hint
+ = typename Triangulation<dim, spacedim>::active_cell_iterator());
/**
* Return a map of index:Point<spacedim>, containing the used vertices of the
#include <deal.II/base/point.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/thread_local_storage.h>
+#include <deal.II/grid/grid_tools_cache.h>
#include <deal.II/lac/vector.h>
* @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
*/
const Mapping<dim> &mapping;
+ /**
+ * The Cache object
+ */
+ GridTools::Cache<dim,DoFHandlerType::space_dimension> cache;
+
/**
* The latest cell hint.
*/
#include <deal.II/base/utilities.h>
#include <deal.II/base/logstream.h>
#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_tools_cache.h>
#include <deal.II/hp/fe_collection.h>
#include <deal.II/hp/fe_values.h>
#include <deal.II/hp/mapping_collection.h>
#include <deal.II/numerics/fe_field_function.h>
#include <deal.II/numerics/vector_tools.h>
+#include <tuple>
+
+
DEAL_II_NAMESPACE_OPEN
dh(&mydh, "FEFieldFunction"),
data_vector(myv),
mapping(mymapping),
+ cache(dh->get_triangulation(),mymapping),
cell_hint(dh->end())
{
}
std::vector<std::vector<Point<dim> > > &qpoints,
std::vector<std::vector<unsigned int> > &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<bool> 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<Point<dim> >
- qp = get_reference_coordinates (cell, points[0]);
- if (!qp)
- {
- const std::pair<typename dealii::internal::ActiveCellIterator<dim, dim, DoFHandlerType>::type, Point<dim> >
- 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<np; ++p)
- if (point_flags[p] == false)
- {
- // same logic as above
- const boost::optional<Point<dim> >
- 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<typename dealii::internal::ActiveCellIterator<dim, dim, DoFHandlerType>::type, Point<dim> > 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<c; ++n)
- {
- Assert(qpoints[n].size() == maps[n].size(),
- ExcDimensionMismatch(qpoints[n].size(), maps[n].size()));
- qps += qpoints[n].size();
- }
- Assert(qps == np,
- ExcDimensionMismatch(qps, np));
-#endif
-
- return c;
+ // Calling the GridTools routine and preparing output
+ auto cell_qpoint_map = GridTools::compute_point_locations(cache,points,cell_hint.get());
+ const auto &tria_cells = std::get<0>(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();
}
std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator >,
std::vector< std::vector< Point<dim> > >,
std::vector< std::vector<unsigned int> > >
- compute_point_locations(const Cache<dim,spacedim> &cache,
- const std::vector<Point<spacedim> > &points)
+ compute_point_locations(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();
// We begin by finding the cell/transform of the first point
std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator, Point<dim> >
- 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);
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
// Initializing the cache
GridTools::Cache<dim,dim> cache(tria);
- std::tuple< std::vector<typename Triangulation<dim, dim>::active_cell_iterator >, std::vector< std::vector< Point<dim> > >, std::vector<std::vector<unsigned int> > >
- 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;
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <iostream>
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/dofs/dof_handler.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/dofs/dof_accessor.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/numerics/fe_field_function.h>
+
+#include <deal.II/fe/mapping_q.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)
+ {
+ Point<dim> p;
+ for (unsigned int d=0; d<dim; ++d)
+ p[d] = double(Testing::rand())/RAND_MAX; //Normalizing the value
+ points.push_back(p);
+ }
+
+ // Initializing the cache
+ GridTools::Cache<dim,dim> 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<std::get<0>(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<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;
+ }
+ }
+ 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);
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff