--- /dev/null
+Added: Function GridTools::compute_point_locations and a test for it. The function uses the recently added Cache to improve performance.
+<br>
+(Giovanni Alzetta, 2017/10/24)
*/
/*@{*/
+ /**
+ * 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 <int dim, int spacedim>
+ std::tuple<
+ 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,
+ const double &distance_factor=0.5);
+
/**
* Return a map of index:Point<spacedim>, containing the used vertices of the
* given `container`. The key of the returned map is the global index in the
+ template <int dim, int spacedim>
+ std::tuple<
+ 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,
+ const double &distance_factor)
+ {
+ // How many points are here?
+ const unsigned int np = points.size();
+
+ std::tuple<
+ std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator >,
+ std::vector< std::vector< Point<dim> > >,
+ std::vector< std::vector<unsigned int> > >
+ 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<bool> point_flags(np, false);
+
+ // Classify the first point:
+ // find active cell returns a pair (cell,transformed point)
+ const std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
+ Point<dim> >
+ 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<spacedim> cell_center = my_pair.first->center();
+ double cell_diameter = my_pair.first->diameter()*
+ (distance_factor + std::numeric_limits<double>::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<np; ++p)
+ if (point_flags[p] == false)
+ {
+ // Assume the point is in the cell
+ flag_outside = false;
+
+ // If the point is too far from the cell center it can't be inside it
+ if ( cell_center.distance(points[p]) > cell_diameter )
+ flag_outside = true;
+ else
+ {
+ //The point is close to the cell: try transforming it
+ try
+ {
+ Point<dim> qp =cache.get_mapping().transform_real_to_unit_cell
+ (std::get<0>(cell_qpoint_map)[c], points[p]);
+ if (GeometryInfo<dim>::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<dim>::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<dim, spacedim>::active_cell_iterator, Point<dim> >
+ 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<double>::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<c; ++n)
+ {
+ Assert(std::get<1>(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<int dim, int spacedim>
std::map<unsigned int, Point<spacedim> >
extract_used_vertices(const Triangulation<dim, spacedim> &container,
const typename Triangulation<deal_II_dimension,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 > > >
+ compute_point_locations(const Cache< deal_II_dimension, deal_II_space_dimension > &,
+ const std::vector< Point< deal_II_space_dimension > > &,
+ const double &);
+ \}
#endif
}
--- /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);
+
+ 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);
+ 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