* 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
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.
}
+
template <int dim, template <int, int> class MeshType, int spacedim>
unsigned int
find_closest_vertex (const MeshType<dim,spacedim> &mesh,
}
+
+ 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>
// 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.
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);
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
+
--- /dev/null
+
+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