]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Added new function compute point locations try all
authorGiovanni Alzetta <giovannialzetta@hotmail.it>
Wed, 9 Jan 2019 14:32:24 +0000 (15:32 +0100)
committerGiovanni Alzetta <giovannialzetta@hotmail.it>
Thu, 14 Feb 2019 08:34:15 +0000 (09:34 +0100)
include/deal.II/grid/grid_tools.h
source/grid/grid_tools.cc
source/grid/grid_tools.inst.in
tests/grid/compute_point_locations_try_all_01.cc [new file with mode: 0644]
tests/grid/compute_point_locations_try_all_01.output [new file with mode: 0644]

index 3a89f434b7519809a36afc5ce484b87d855b192d..809de1bccad4162b2397196f50788ae99857319d 100644 (file)
@@ -720,6 +720,9 @@ namespace GridTools
    * they belong are the best case. Pre-sorting points, trying to minimize
    * distances between them, might make the function extremely faster.
    *
+   * @note If a point is not found inside the mesh, or is lying inside an
+   * artificial cell of a parallel::Triangulation, an exception is thrown.
+   *
    * @note The actual return type of this function, i.e., the type referenced
    * above as @p return_type, is
    * @code
@@ -750,6 +753,44 @@ namespace GridTools
       &cell_hint =
         typename Triangulation<dim, spacedim>::active_cell_iterator());
 
+  /**
+   * This function is similar to GridTools::compute_point_locations(),
+   * but it tries to find and transform every point of @p points.
+   *
+   * @return A tuple containing four elements; the first three
+   * are documented in GridTools::compute_point_locations().
+   * The last element of the @p return_type contains the
+   * indices of points which are not not found inside the mesh
+   * or lie on artificial cells.
+   *
+   * @code
+   *   std::tuple<
+   *     std::vector<typename Triangulation<dim,
+   * spacedim>::active_cell_iterator>, std::vector<std::vector<Point<dim>>>,
+   *     std::vector<std::vector<unsigned int>>,
+   *     std::vector<unsigned int>>
+   * @endcode
+   *
+   * For a more detailed documentation see
+   * GridTools::compute_point_locations().
+   */
+  template <int dim, int spacedim>
+#  ifndef DOXYGEN
+  std::tuple<
+    std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
+    std::vector<std::vector<Point<dim>>>,
+    std::vector<std::vector<unsigned int>>,
+    std::vector<unsigned int>>
+#  else
+  return_type
+#  endif
+  compute_point_locations_try_all(
+    const Cache<dim, spacedim> &        cache,
+    const std::vector<Point<spacedim>> &points,
+    const typename Triangulation<dim, spacedim>::active_cell_iterator
+      &cell_hint =
+        typename Triangulation<dim, spacedim>::active_cell_iterator());
+
   /**
    * Given a @p cache and a list of
    * @p local_points for each process, find the points lying on the locally
index 201dfe03f9eeabff92323e8065a4437e10bea052..bd97f63dc4b95982720394bb9003030448f84414 100644 (file)
@@ -4193,14 +4193,50 @@ namespace GridTools
     const std::vector<Point<spacedim>> &points,
     const typename Triangulation<dim, spacedim>::active_cell_iterator
       &cell_hint)
+  {
+    const auto cqmp = compute_point_locations_try_all(cache, points, cell_hint);
+    // Splitting the tuple's components
+    auto &cells          = std::get<0>(cqmp);
+    auto &qpoints        = std::get<1>(cqmp);
+    auto &maps           = std::get<2>(cqmp);
+    auto &missing_points = std::get<3>(cqmp);
+    // If a point was not found, throwing an error, as the old
+    // implementation of compute_point_locations would have done
+    AssertThrow(std::get<3>(cqmp).size() == 0,
+                ExcPointNotFound<spacedim>(points[missing_points[0]]));
+
+    (void)missing_points;
+
+    return {std::move(cells), std::move(qpoints), std::move(maps)};
+  }
+
+
+
+  template <int dim, int spacedim>
+#ifndef DOXYGEN
+  std::tuple<
+    std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
+    std::vector<std::vector<Point<dim>>>,
+    std::vector<std::vector<unsigned int>>,
+    std::vector<unsigned int>>
+#else
+  return_type
+#endif
+  compute_point_locations_try_all(
+    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();
+    const unsigned int        np = points.size();
+    std::vector<unsigned int> points_outside;
 
     std::tuple<
       std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>,
       std::vector<std::vector<Point<dim>>>,
-      std::vector<std::vector<unsigned int>>>
+      std::vector<std::vector<unsigned int>>,
+      std::vector<unsigned int>>
       cell_qpoint_map;
 
     // Now the easy case.
@@ -4211,18 +4247,53 @@ namespace GridTools
     std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
               Point<dim>>
       my_pair;
+
+    bool         found          = false;
+    unsigned int points_checked = 0;
     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]);
+      {
+        try
+          {
+            my_pair = GridTools::find_active_cell_around_point(cache,
+                                                               points[0],
+                                                               cell_hint);
+            found   = true;
+          }
+        catch (const GridTools::ExcPointNotFound<dim> &)
+          {
+            points_outside.emplace_back(0);
+          }
+        ++points_checked;
+      }
+
 
-    std::get<0>(cell_qpoint_map).emplace_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);
+    while (!found && points_checked < np)
+      {
+        try
+          {
+            my_pair =
+              GridTools::find_active_cell_around_point(cache,
+                                                       points[points_checked]);
+            found = true;
+          }
+        catch (const GridTools::ExcPointNotFound<dim> &)
+          {
+            points_outside.emplace_back(points_checked);
+          }
+        // Updating the position of the analyzed points
+        ++points_checked;
+      }
+
+    // If the point has been found in a cell, adding it
+    if (found)
+      {
+        std::get<0>(cell_qpoint_map).emplace_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, points_checked - 1);
+      }
 
     // Now the second easy case.
-    if (np == 1)
+    if (np == points_outside.size())
       return cell_qpoint_map;
     // Computing the cell center and diameter
     Point<spacedim> cell_center = std::get<0>(cell_qpoint_map)[0]->center();
@@ -4230,15 +4301,24 @@ namespace GridTools
                            (0.5 + std::numeric_limits<double>::epsilon());
 
     // Cycle over all points left
-    for (unsigned int p = 1; p < np; ++p)
+    for (unsigned int p = points_checked; p < np; ++p)
       {
         // Checking if the point is close to the cell center, in which
         // case calling find active cell with a cell hint
-        if (cell_center.distance(points[p]) < cell_diameter)
-          my_pair = GridTools::find_active_cell_around_point(
-            cache, points[p], std::get<0>(cell_qpoint_map).back());
-        else
-          my_pair = GridTools::find_active_cell_around_point(cache, points[p]);
+        try
+          {
+            if (cell_center.distance(points[p]) < cell_diameter)
+              my_pair = GridTools::find_active_cell_around_point(
+                cache, points[p], std::get<0>(cell_qpoint_map).back());
+            else
+              my_pair =
+                GridTools::find_active_cell_around_point(cache, points[p]);
+          }
+        catch (const GridTools::ExcPointNotFound<dim> &)
+          {
+            points_outside.push_back(p);
+            continue;
+          }
 
         // Assuming the cell is probably the last cell added
         if (my_pair.first == std::get<0>(cell_qpoint_map).back())
@@ -4298,7 +4378,8 @@ namespace GridTools
     // The number of points in all
     // the cells must be the same as
     // the number of points we
-    // started off from.
+    // started off from,
+    // plus the points which were ignored
     for (unsigned int n = 0; n < c; ++n)
       {
         Assert(std::get<1>(cell_qpoint_map)[n].size() ==
@@ -4307,9 +4388,12 @@ namespace GridTools
                                     std::get<2>(cell_qpoint_map)[n].size()));
         qps += std::get<1>(cell_qpoint_map)[n].size();
       }
-    Assert(qps == np, ExcDimensionMismatch(qps, np));
+
+    Assert(qps + points_outside.size() == np,
+           ExcDimensionMismatch(qps + points_outside.size(), np));
 #endif
 
+    std::get<3>(cell_qpoint_map) = points_outside;
     return cell_qpoint_map;
   }
 
index c8b0ce9120e1095c4bcd418326caae361e2ab886..23dd04b907b6debaae230217d13a17a10d5c1b03 100644 (file)
@@ -93,6 +93,19 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
           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>>,
+                          std::vector<unsigned int>>
+      compute_point_locations_try_all(
+        const Cache<deal_II_dimension, 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 &);
+
       template std::tuple<std::vector<typename Triangulation<
                             deal_II_dimension,
                             deal_II_space_dimension>::active_cell_iterator>,
diff --git a/tests/grid/compute_point_locations_try_all_01.cc b/tests/grid/compute_point_locations_try_all_01.cc
new file mode 100644 (file)
index 0000000..ddb72f2
--- /dev/null
@@ -0,0 +1,143 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// Test GridTools::compute_point_locations_try_all
+// Test for corner case: vector of points, some of which
+// lie outside the mesh
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q.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/grid/tria.h>
+
+#include <deal.II/numerics/fe_field_function.h>
+
+#include <iostream>
+
+#include "../tests.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)
+    {
+      points.push_back(random_point<dim>());
+      // One every three points is
+      // shifted outside the mesh
+      if (points.size() % 3 == dim - 1)
+        {
+          points.back()[0] += 5.0;
+        }
+    }
+
+  // Initializing the cache
+  GridTools::Cache<dim, dim> cache(tria);
+
+  auto cell_qpoint_map =
+    GridTools::compute_point_locations_try_all(cache, points);
+  const auto &cells   = std::get<0>(cell_qpoint_map);
+  const auto &qpoints = std::get<1>(cell_qpoint_map);
+  const auto &maps    = std::get<2>(cell_qpoint_map);
+  const auto &other_p = std::get<3>(cell_qpoint_map);
+  size_t      n_cells = cells.size();
+
+  deallog << "Points found in " << n_cells << " cells" << std::endl;
+
+  unsigned int checked_points = 0;
+  // testing if the transformation is correct:
+  for (unsigned int i = 0; i < n_cells; ++i)
+    {
+      const auto &cell      = cells[i];
+      const auto &quad      = qpoints[i];
+      const auto &local_map = maps[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;
+          else
+            ++checked_points;
+        }
+    }
+  for (auto i : other_p)
+    if (i % 3 != dim - 2)
+      {
+        deallog << "Error: point index " << i << " was not found!" << std::endl;
+      }
+
+  if (checked_points + other_p.size() == points.size())
+    deallog << "All points where processed correctly" << std::endl;
+  else
+    {
+      deallog << "Error: some points where not processed correctly"
+              << std::endl;
+      deallog << "Points found inside the mesh: " << checked_points
+              << std::endl;
+      deallog << "Discarded points: " << points.size() - checked_points
+              << std::endl;
+      deallog << "Points outside the mesh: " << other_p.size() << std::endl;
+      deallog << "Total points processed: " << points.size() << std::endl;
+    }
+  deallog << "Test finished" << std::endl;
+}
+
+int
+main()
+{
+  initlog();
+
+  deallog << "Deal.II compute_pt_loc_try_all:" << std::endl;
+  test_compute_pt_loc<2>(101);
+  test_compute_pt_loc<3>(202);
+}
diff --git a/tests/grid/compute_point_locations_try_all_01.output b/tests/grid/compute_point_locations_try_all_01.output
new file mode 100644 (file)
index 0000000..1272686
--- /dev/null
@@ -0,0 +1,12 @@
+
+DEAL::Deal.II compute_pt_loc_try_all:
+DEAL::Testing for dim = 2
+DEAL::Testing on: 101 points.
+DEAL::Points found in 58 cells
+DEAL::All points where processed correctly
+DEAL::Test finished
+DEAL::Testing for dim = 3
+DEAL::Testing on: 202 points.
+DEAL::Points found in 119 cells
+DEAL::All points where processed correctly
+DEAL::Test finished

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.