]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Linked fe_field_function compute point locations to gridtools 5411/head
authorGiovanni Alzetta <giovannialzetta@hotmail.it>
Mon, 6 Nov 2017 11:55:06 +0000 (11:55 +0000)
committerGiovanni Alzetta <giovannialzetta@hotmail.it>
Mon, 11 Dec 2017 16:44:08 +0000 (17:44 +0100)
doc/news/changes/minor/20171106GiovanniAlzetta [new file with mode: 0644]
include/deal.II/grid/grid_tools.h
include/deal.II/numerics/fe_field_function.h
include/deal.II/numerics/fe_field_function.templates.h
source/grid/grid_tools.cc
source/grid/grid_tools.inst.in
tests/grid/compute_point_locations_01.cc
tests/grid/compute_point_locations_02.cc [new file with mode: 0644]
tests/grid/compute_point_locations_02.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20171106GiovanniAlzetta b/doc/news/changes/minor/20171106GiovanniAlzetta
new file mode 100644 (file)
index 0000000..1803b0f
--- /dev/null
@@ -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)
+<br>
+(Giovanni Alzetta, 2017/11/06)
index 8fe24838c1947f93fa663411e2ec2e6047a9b23c..40161dab24a1dc0f479b4c86b08bbf054fe3f900 100644 (file)
@@ -657,7 +657,9 @@ namespace GridTools
       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
index 119a7d43b400dbc8231283146324357669d25e85..7e821095842581d90709c7ca6f329952746c1b97 100644 (file)
@@ -24,6 +24,7 @@
 #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>
 
@@ -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<dim> &mapping;
 
+    /**
+     * The Cache object
+     */
+    GridTools::Cache<dim,DoFHandlerType::space_dimension> cache;
+
     /**
      * The latest cell hint.
      */
index 7acf5f9261e6039e50b4643c12e47afb3b8307d2..e26d997d1ef4c6928c13a973faafb4cc8f1fbc37 100644 (file)
@@ -20,6 +20,7 @@
 #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>
@@ -28,6 +29,9 @@
 #include <deal.II/numerics/fe_field_function.h>
 #include <deal.II/numerics/vector_tools.h>
 
+#include <tuple>
+
+
 
 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<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();
   }
 
 
index 4d9f8e906625cae7d7ff78698a4ef873c72a43f5..ecf5684ada805566133179dfa7831a81c78cae06 100644 (file)
@@ -5176,8 +5176,9 @@ next_cell:
   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();
@@ -5193,8 +5194,13 @@ next_cell:
 
     // 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);
index 39cb48ebfeade43588873ae7c9702644a98c19ee..51d4df9a32c8027d4cde74ef6b3efecb51c68004 100644 (file)
@@ -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
index ab0bca3fe8bd1bef9062582a51308a7d1dc598e1..981142f0593524c2a697b93b4fa63a1a912d55db 100644 (file)
@@ -61,8 +61,7 @@ void test_compute_pt_loc(unsigned int n_points)
   // 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;
diff --git a/tests/grid/compute_point_locations_02.cc b/tests/grid/compute_point_locations_02.cc
new file mode 100644 (file)
index 0000000..bbfdb40
--- /dev/null
@@ -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 <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);
+}
diff --git a/tests/grid/compute_point_locations_02.output b/tests/grid/compute_point_locations_02.output
new file mode 100644 (file)
index 0000000..0a6e0e1
--- /dev/null
@@ -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

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.