DoFHandler<deal_II_dimension, deal_II_space_dimension>;
hp::DoFHandler<deal_II_dimension, deal_II_space_dimension> }
+// concrete Triangulation types (all types you can iterate over cells from)
+// with hard-coded <deal_II_dimension, deal_II_space_dimension>
+TRIANGULATIONS := { Triangulation<deal_II_dimension, deal_II_space_dimension>;
+ parallel::shared::Triangulation<deal_II_dimension, deal_II_space_dimension>;
+ parallel::distributed::Triangulation<deal_II_dimension, deal_II_space_dimension>;}
+
// serial and parallel sparsity pattern types
SPARSITY_PATTERNS := { SparsityPattern;
DynamicSparsityPattern;
--- /dev/null
+New: A new version of GridTools::find_active_cell_around_point()
+has been added that exploits local maps constructed
+using standard GridTools utilities.
+<br>
+(Rene Gassmoeller, Luca Heltai, 2017/08/13)
const Point<spacedim> &p,
const std::vector<bool> &marked_vertices = std::vector<bool>());
+ /**
+ * A version of the previous function that exploits an already existing
+ * map between vertices and cells, constructed using the function
+ * GridTools::vertex_to_cell_map, a map of vertex_to_cell_centers, obtained
+ * through GridTools::vertex_to_cells_center_directions, and a guess `cell_hint`.
+ *
+ * @author Luca Heltai, Rene Gassmoeller, 2017
+ */
+ template <int dim, template <int, int> class MeshType, int spacedim>
+#ifndef _MSC_VER
+ std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim> >
+#else
+ std::pair<typename dealii::internal::ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim> >::type, Point<dim> >
+#endif
+ find_active_cell_around_point (const Mapping<dim,spacedim> &mapping,
+ const MeshType<dim,spacedim> &mesh,
+ const Point<spacedim> &p,
+ const std::vector<std::set<typename MeshType<dim,spacedim>::active_cell_iterator > > &vertex_to_cell_map,
+ const std::vector<std::vector<Tensor<1,spacedim> > > &vertex_to_cell_centers,
+ const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint=typename MeshType<dim, spacedim>::active_cell_iterator(),
+ const std::vector<bool> &marked_vertices = std::vector<bool>());
+
/**
* A version of the previous function where we use that mapping on a given
* cell that corresponds to the active finite element index of that cell.
vertex_to_cell_map(const Triangulation<dim,spacedim> &triangulation);
/**
- * Returns a vector that contains a tensor for every vertex-cell combination
- * of the output of GridTools::vertex_to_cell_map() (which is expected as
- * input parameter for this function). Each tensor represents a geometric
- * vector from the vertex to the respective cell center.
+ * Returns a vector of normalized tensors for each vertex-cell combination of
+ * the output of GridTools::vertex_to_cell_map() (which is expected as input
+ * parameter for this function). Each tensor represents a geometric vector
+ * from the vertex to the respective cell center.
+ *
+ * An assertion will be thrown if the size of the input vector is not equal to
+ * the number of vertices of the triangulation.
+ *
+ * result[v][c] is a unit Tensor for vertex index v, indicating the direction of
+ * the center of the c-th cell with respect to the vertex v.
*
* @author Rene Gassmoeller, Luca Heltai, 2017.
*/
*/
template <int dim, int spacedim>
unsigned int
- get_closest_vertex_of_cell(const typename Triangulation<dim,spacedim>::active_cell_iterator &cell,
- const Point<spacedim> &position);
+ find_closest_vertex_of_cell(const typename Triangulation<dim,spacedim>::active_cell_iterator &cell,
+ const Point<spacedim> &position);
/**
* Compute a globally unique index for each vertex and hanging node
const std::vector<Point<spacedim> > &vertices = mesh.get_vertices();
const unsigned int n_vertices = vertex_to_cells.size();
+ AssertDimension(vertices.size(), n_vertices);
+
+
std::vector<std::vector<Tensor<1,spacedim> > > vertex_to_cell_centers(n_vertices);
for (unsigned int vertex=0; vertex<n_vertices; ++vertex)
if (mesh.vertex_used(vertex))
}
+ namespace
+ {
+ template <int spacedim>
+ bool
+ compare_point_association(const unsigned int a,
+ const unsigned int b,
+ const Tensor<1,spacedim> &point_direction,
+ const std::vector<Tensor<1,spacedim> > ¢er_directions)
+ {
+ const double scalar_product_a = center_directions[a] * point_direction;
+ const double scalar_product_b = center_directions[b] * point_direction;
+
+ // The function is supposed to return if a is before b. We are looking
+ // for the alignment of point direction and center direction, therefore
+ // return if the scalar product of a is larger.
+ return (scalar_product_a > scalar_product_b);
+ }
+ }
+
+ template <int dim, template <int, int> class MeshType, int spacedim>
+#ifndef _MSC_VER
+ std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim> >
+#else
+ std::pair<typename dealii::internal::ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim> >::type, Point<dim> >
+#endif
+ find_active_cell_around_point (const Mapping<dim,spacedim> &mapping,
+ const MeshType<dim,spacedim> &mesh,
+ const Point<spacedim> &p,
+ const std::vector<std::set<typename MeshType<dim,spacedim>::active_cell_iterator > > &vertex_to_cells,
+ const std::vector<std::vector<Tensor<1,spacedim> > > &vertex_to_cell_centers,
+ const typename MeshType<dim, spacedim>::active_cell_iterator &cell_hint ,
+ const std::vector<bool> &marked_vertices)
+ {
+ std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim> > cell_and_position;
+
+ bool found_cell = false;
+
+ unsigned int closest_vertex_index = 0;
+ Tensor<1,spacedim> vertex_to_point;
+ auto current_cell = cell_hint;
+
+ while (found_cell == false)
+ {
+ // First look at the vertices of the cell cell_hint. If it's an
+ // invalid cell, then query for the closest global vertex
+ if (current_cell.state() == IteratorState::valid)
+ {
+ const unsigned int closest_vertex = find_closest_vertex_of_cell<dim,spacedim>(current_cell , p);
+ vertex_to_point = p - current_cell ->vertex(closest_vertex);
+ closest_vertex_index = current_cell ->vertex_index(closest_vertex);
+ }
+ else
+ {
+ closest_vertex_index = GridTools::find_closest_vertex(mesh,p,marked_vertices);
+ vertex_to_point = p - mesh.get_vertices()[closest_vertex_index];
+ }
+
+ vertex_to_point /= vertex_to_point.norm();
+ const unsigned int n_neighbor_cells = vertex_to_cells[closest_vertex_index].size();
+
+ // Create a corresponding map of vectors from vertex to cell center
+ std::vector<unsigned int> neighbor_permutation(n_neighbor_cells);
+
+ for (unsigned int i=0; i<n_neighbor_cells; ++i)
+ neighbor_permutation[i] = i;
+
+ auto comp = [&](const unsigned int a, const unsigned int b) -> bool
+ {
+ return compare_point_association<spacedim>(a,b,vertex_to_point,vertex_to_cell_centers[closest_vertex_index]);
+ };
+
+ std::sort(neighbor_permutation.begin(),
+ neighbor_permutation.end(),
+ comp);
+
+ // Search all of the cells adjacent to the closest vertex of the cell hint
+ // Most likely we will find the point in them.
+ for (unsigned int i=0; i<n_neighbor_cells; ++i)
+ {
+ try
+ {
+ auto cell = vertex_to_cells[closest_vertex_index].begin();
+ std::advance(cell,neighbor_permutation[i]);
+ const Point<dim> p_unit = mapping.transform_real_to_unit_cell(*cell, p);
+ if (GeometryInfo<dim>::is_inside_unit_cell(p_unit))
+ {
+ cell_and_position.first = *cell;
+ cell_and_position.second = p_unit;
+ found_cell = true;
+ break;
+ }
+ }
+ catch (typename Mapping<dim>::ExcTransformationFailed &)
+ {}
+ }
+
+ if (found_cell == true)
+ return cell_and_position;
+
+ // The first time around, we check for vertices in the hint_cell. If that
+ // does not work, we set the cell iterator to an invalid one, and look
+ // for a global vertex close to the point. If that does not work, we are in
+ // trouble, and just throw an exception.
+ //
+ // If we got here, then we did not find the point. If the
+ // current_cell.state() here is not IteratorState::valid, it means that
+ // the user did not provide a hint_cell, and at the beginning of the
+ // while loop we performed an actual global search on the mesh
+ // vertices. Not finding the point then means the point is outside the
+ // domain.
+ AssertThrow(current_cell.state() == IteratorState::valid,
+ ExcPointNotFound<spacedim>(p));
+
+ current_cell = typename MeshType<dim,spacedim>::active_cell_iterator();
+ }
+ return cell_and_position;
+ }
+
+
template <int dim, int spacedim>
unsigned int
- get_closest_vertex_of_cell(const typename Triangulation<dim,spacedim>::active_cell_iterator &cell,
- const Point<spacedim> &position)
+ find_closest_vertex_of_cell(const typename Triangulation<dim,spacedim>::active_cell_iterator &cell,
+ const Point<spacedim> &position)
{
- double minimum_distance = std::numeric_limits<double>::max();
- unsigned int closest_vertex = numbers::invalid_unsigned_int;
+ double minimum_distance = position.distance_square(cell->vertex(0));
+ unsigned int closest_vertex = 0;
- for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+ for (unsigned int v=1; v<GeometryInfo<dim>::vertices_per_cell; ++v)
{
- const double vertex_distance = position.distance(cell->vertex(v));
+ const double vertex_distance = position.distance_square(cell->vertex(v));
if (vertex_distance < minimum_distance)
{
closest_vertex = v;
//
// ---------------------------------------------------------------------
+for (X : TRIANGULATIONS; deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS)
+{
+
+#if deal_II_dimension <= deal_II_space_dimension
+ namespace GridTools \{
+ template
+ std::pair<dealii::internal::ActiveCellIterator<deal_II_dimension, deal_II_space_dimension, X>::type, Point<deal_II_dimension> >
+ find_active_cell_around_point (const Mapping<deal_II_dimension, deal_II_space_dimension> &,
+ const X &,
+ const Point<deal_II_space_dimension> &,
+ const std::vector<std::set<typename dealii::internal::ActiveCellIterator<deal_II_dimension, deal_II_space_dimension, X>::type> > &,
+ const std::vector<std::vector<Tensor<1,deal_II_space_dimension> > > &,
+ const dealii::internal::ActiveCellIterator<deal_II_dimension, deal_II_space_dimension, X>::type &,
+ const std::vector<bool> &);
+ \}
+
+#endif
+}
for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II_space_dimension : SPACE_DIMENSIONS)
template
std::map<unsigned int,types::global_vertex_index>
compute_local_to_global_vertex_index_map(const parallel::distributed::Triangulation<deal_II_dimension,deal_II_space_dimension> &triangulation);
+
+
\}
#endif
std::vector<std::set<Triangulation<deal_II_dimension,deal_II_space_dimension>::active_cell_iterator> >
vertex_to_cell_map(const Triangulation<deal_II_dimension,deal_II_space_dimension> &triangulation);
+ template
+ std::vector<std::vector<Tensor<1,deal_II_space_dimension> > >
+ vertex_to_cell_centers_directions(const Triangulation<deal_II_dimension,deal_II_space_dimension> &mesh,
+ const std::vector<std::set<typename Triangulation<deal_II_dimension,
+ deal_II_space_dimension>::active_cell_iterator> > &vertex_to_cells);
+
#if deal_II_dimension == deal_II_space_dimension
# if deal_II_dimension > 1
template
--- /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.
+//
+// ---------------------------------------------------------------------
+
+
+/*
+ * Given the number of refinements and the number of random points
+ * it benchmarks the time needed to run the function FCT
+ * which can be point_locator_D2 (or point_locator when it shall be written)
+ */
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_generator.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/hp/q_collection.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/fe_field_function.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/numerics/vector_tools.h>
+
+
+template <int dim, int spacedim>
+double test(unsigned int n_ref, unsigned int n_points)
+{
+ deallog << "Testing " << dim << ", " << spacedim << std::endl;
+
+ Triangulation<dim,spacedim> tria;
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(n_ref);
+ DoFHandler<dim,spacedim> dof_handler(tria);
+
+ std::vector<Point<spacedim>> points;
+
+ deallog << "Points in study: " << n_points << std::endl;
+ for (size_t i=0; i<n_points; ++i)
+ {
+ //We need points in the square [0,1]x[0,1]: this is achieved normalizing with RAND_MAX
+ //RAND_MAX for the used algorithm is 2147483647
+ Point<spacedim> p;
+ for (unsigned int d=0; d<spacedim; ++d)
+ p[d] = double(Testing::rand())/RAND_MAX; //Normalizing the value
+ points.push_back(p);
+ }
+
+ auto v_to_c = GridTools::vertex_to_cell_map (tria);
+ auto v_to_c_d = GridTools::vertex_to_cell_centers_directions (tria, v_to_c);
+
+ auto &mapping = StaticMappingQ1<dim,spacedim>::mapping;
+ auto cell = tria.begin_active();
+ for (auto &p : points)
+ {
+ auto c_and_p = GridTools::find_active_cell_around_point(mapping, tria, p,
+ v_to_c, v_to_c_d);
+ auto p2 = mapping.transform_unit_to_real_cell(c_and_p.first, c_and_p.second);
+ if (p2.distance(p)>1e-10)
+ deallog << "NOT OK!" << p << ", " << p2
+ << ", " << c_and_p.first << std::endl;
+ cell = c_and_p.first;
+ }
+ deallog << "OK" << std::endl;
+}
+
+int main ()
+{
+ initlog();
+
+ test<1,1> (4,10);
+ test<2,2> (3,20);
+ test<3,3> (2,30);
+}
--- /dev/null
+
+DEAL::Testing 1, 1
+DEAL::Points in study: 10
+DEAL::OK
+DEAL::Testing 2, 2
+DEAL::Points in study: 20
+DEAL::OK
+DEAL::Testing 3, 3
+DEAL::Points in study: 30
+DEAL::OK
projects / dealii.git / commitdiff