]> https://gitweb.dealii.org/ - dealii.git/commitdiff
GridTools::find_closest_vertex with mapping argument.
authorLuca Heltai <luca.heltai@sissa.it>
Tue, 19 Sep 2017 14:37:50 +0000 (16:37 +0200)
committerLuca Heltai <luca.heltai@sissa.it>
Wed, 27 Sep 2017 13:07:33 +0000 (15:07 +0200)
include/deal.II/grid/grid_tools.h
source/grid/grid_tools.cc
source/grid/grid_tools.inst.in
tests/grid/grid_tools_find_closest_vertex_01.cc [new file with mode: 0644]
tests/grid/grid_tools_find_closest_vertex_01.output [new file with mode: 0644]

index f6fdc81b454fbc412aa8b4d29cd175f9b889fda1..9586e715c861ee23e722e1c36a0faccea403025f 100644 (file)
@@ -601,6 +601,8 @@ namespace GridTools
    *        extract_used_vertices()
    * @param p The target point
    * @return The index of the vertex that is closest to the target point `p`
+   *
+   * @author Luca Heltai, 2017.
    */
   template<int spacedim>
   unsigned int
@@ -633,6 +635,35 @@ namespace GridTools
                        const Point<spacedim>         &p,
                        const std::vector<bool>       &marked_vertices = std::vector<bool>());
 
+  /**
+   * Find and return the number of the used vertex (or marked vertex) in a
+   * given mesh that is located closest to a given point. Use the given
+   * mapping to compute the actual location of the vertices.
+   *
+   * @param mapping A mapping used to compute the vertex locations
+   * @param mesh A variable of a type that satisfies the requirements of the
+   * @ref ConceptMeshType "MeshType concept".
+   * @param p The point for which we want to find the closest vertex.
+   * @param marked_vertices An array of bools indicating which
+   * vertices of @p mesh will be considered within the search
+   * as the potentially closest vertex. On receiving a non-empty
+   * @p marked_vertices, the function will
+   * only search among @p marked_vertices for the closest vertex.
+   * The size of this array should be equal to the value returned by
+   * Triangulation::n_vertices() for the triangulation underlying the given mesh
+   * (as opposed to the value returned by Triangulation::n_used_vertices()).
+   * @return The index of the closest vertex found.
+   *
+   * @author Luca Heltai, 2017
+   */
+  template <int dim, template <int, int> class MeshType, int spacedim>
+  unsigned int
+  find_closest_vertex (const Mapping<dim,spacedim>   &mapping,
+                       const MeshType<dim, spacedim> &mesh,
+                       const Point<spacedim>         &p,
+                       const std::vector<bool>       &marked_vertices = std::vector<bool>());
+
+
   /**
    * Find and return a vector of iterators to active cells that surround a
    * given vertex with index @p vertex_index.
index cac93c628a8d5ab37e2497bf21a4ffc7eee260c5..612d0391261bd5ca6468fb17535717bd5165664c 100644 (file)
@@ -1048,6 +1048,7 @@ namespace GridTools
   }
 
 
+
   template <int dim, template <int, int> class MeshType, int spacedim>
   unsigned int
   find_closest_vertex (const MeshType<dim,spacedim> &mesh,
@@ -1122,6 +1123,67 @@ namespace GridTools
   }
 
 
+
+  template <int dim, template <int, int> class MeshType, int spacedim>
+  unsigned int
+  find_closest_vertex (const Mapping<dim,spacedim>  &mapping,
+                       const MeshType<dim,spacedim> &mesh,
+                       const Point<spacedim>        &p,
+                       const std::vector<bool>      &marked_vertices)
+  {
+    // Take a shortcut in the simple case.
+    if (mapping.preserves_vertex_locations() == true)
+      return find_closest_vertex(mesh, p, marked_vertices);
+
+    // first get the underlying
+    // triangulation from the
+    // mesh and determine vertices
+    // and used vertices
+    const Triangulation<dim, spacedim> &tria = mesh.get_triangulation();
+
+    auto vertices = extract_used_vertices(tria, mapping);
+
+    Assert ( tria.get_vertices().size() == marked_vertices.size() || marked_vertices.size() ==0,
+             ExcDimensionMismatch(tria.get_vertices().size(), marked_vertices.size()));
+
+    // If p is an element of marked_vertices,
+    // and q is that of used_Vertices,
+    // the vector marked_vertices does NOT
+    // contain unused vertices if p implies q.
+    // I.e., if p is true q must be true
+    // (if p is false, q could be false or true).
+    // p implies q logic is encapsulated in ~p|q.
+    Assert( marked_vertices.size()==0
+            ||
+            std::equal( marked_vertices.begin(),
+                        marked_vertices.end(),
+                        tria.get_used_vertices().begin(),
+                        [](bool p, bool q)
+    {
+      return !p || q;
+    }),
+    ExcMessage("marked_vertices should be a subset of used vertices in the triangulation "
+               "but marked_vertices contains one or more vertices that are not used vertices!") );
+
+    // Remove from the map unwanted elements.
+    if (marked_vertices.size())
+      for (auto it = vertices.begin(); it != vertices.end(); )
+        {
+          if (marked_vertices[it->first] == false)
+            {
+              vertices.erase(it++);
+            }
+          else
+            {
+              ++it;
+            }
+        }
+
+    return find_closest_vertex(vertices, p);
+  }
+
+
+
   template <int dim, template <int, int> class MeshType, int spacedim>
 #ifndef _MSC_VER
   std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
@@ -1362,7 +1424,7 @@ next_cell:
     // all adjacent cells
     std::vector<active_cell_iterator> adjacent_cells_tmp
       = find_cells_adjacent_to_vertex(mesh,
-                                      find_closest_vertex(mesh, p, marked_vertices));
+                                      find_closest_vertex(mapping, mesh, p, marked_vertices));
 
     // Make sure that we have found
     // at least one cell adjacent to vertex.
index 3938a71118085483d6b6117653c3ade982490314..7c47ceaf4e740e1799a732050592e9b8ab66629e 100644 (file)
@@ -45,6 +45,13 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS ; deal_II
                          const Point<deal_II_space_dimension> &,
                          const std::vector<bool> &);
 
+    template
+    unsigned int
+    find_closest_vertex (const Mapping<deal_II_dimension, deal_II_space_dimension>&,
+                         const X &,
+                         const Point<deal_II_space_dimension> &,
+                         const std::vector<bool> &);
+
     template
     std::vector<dealii::internal::ActiveCellIterator<deal_II_dimension, deal_II_space_dimension, X>::type>
     find_cells_adjacent_to_vertex(const X &, const unsigned int);
diff --git a/tests/grid/grid_tools_find_closest_vertex_01.cc b/tests/grid/grid_tools_find_closest_vertex_01.cc
new file mode 100644 (file)
index 0000000..01fcac6
--- /dev/null
@@ -0,0 +1,76 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2001 - 2015 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.
+//
+// ---------------------------------------------------------------------
+
+// Check extract used_vertices and find_closest_vertex when a mapping
+// is present
+
+#include "../tests.h"
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q_eulerian.h>
+#include <deal.II/dofs/dof_handler.h>
+
+
+template <int dim, int spacedim>
+void test (const Point<spacedim> &p)
+{
+  deallog << "dim: " << dim << ", spacedim: "
+          << spacedim << std::endl;
+
+  Triangulation<dim,spacedim> tria;
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(1);
+
+  // Vector fe
+  FESystem<dim,spacedim> fe(FE_Q<dim, spacedim>(1)^spacedim);
+  DoFHandler<dim,spacedim> dh(tria);
+  dh.distribute_dofs(fe);
+
+  // Displacement vector.
+  Vector<double> disp(dh.n_dofs());
+  for (unsigned int i=0; i<dh.n_dofs(); ++i)
+    disp[i] = .5;
+
+  // mapping
+  MappingQEulerian<dim, Vector<double>, spacedim> mapping(1, dh, disp);
+
+  auto m = GridTools::extract_used_vertices(tria, mapping);
+
+  std::vector<bool> selected_vertices(tria.n_vertices(), false);
+  selected_vertices[0] = true;
+  selected_vertices[5] = true;
+
+  for (auto &e: m)
+    deallog << "Vertex: " << e.first << ": " << e.second << std::endl;
+
+  auto i = GridTools::find_closest_vertex(mapping,tria,p,selected_vertices);
+  deallog << "Closest vertex to " << p
+          << ", v["<< i << "] :"
+          << m[i] << std::endl;
+};
+
+
+int main ()
+{
+  initlog();
+  test<2,2> (Point<2>(.2,.2));
+  test<2,2> (Point<2>(.6,.9));
+  return 0;
+}
+
diff --git a/tests/grid/grid_tools_find_closest_vertex_01.output b/tests/grid/grid_tools_find_closest_vertex_01.output
new file mode 100644 (file)
index 0000000..0331bd7
--- /dev/null
@@ -0,0 +1,23 @@
+
+DEAL::dim: 2, spacedim: 2
+DEAL::Vertex: 0: 0.500000 0.500000
+DEAL::Vertex: 1: 1.50000 0.500000
+DEAL::Vertex: 2: 0.500000 1.50000
+DEAL::Vertex: 3: 1.50000 1.50000
+DEAL::Vertex: 4: 1.00000 0.500000
+DEAL::Vertex: 5: 0.500000 1.00000
+DEAL::Vertex: 6: 1.50000 1.00000
+DEAL::Vertex: 7: 1.00000 1.50000
+DEAL::Vertex: 8: 1.00000 1.00000
+DEAL::Closest vertex to 0.200000 0.200000, v[0] :0.500000 0.500000
+DEAL::dim: 2, spacedim: 2
+DEAL::Vertex: 0: 0.500000 0.500000
+DEAL::Vertex: 1: 1.50000 0.500000
+DEAL::Vertex: 2: 0.500000 1.50000
+DEAL::Vertex: 3: 1.50000 1.50000
+DEAL::Vertex: 4: 1.00000 0.500000
+DEAL::Vertex: 5: 0.500000 1.00000
+DEAL::Vertex: 6: 1.50000 1.00000
+DEAL::Vertex: 7: 1.00000 1.50000
+DEAL::Vertex: 8: 1.00000 1.00000
+DEAL::Closest vertex to 0.600000 0.900000, v[5] :0.500000 1.00000

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.