From: Bruno Turcksin <bruno.turcksin@gmail.com>
Date: Wed, 30 Sep 2015 17:16:21 +0000 (-0500)
Subject: Add vertex to cells map.
X-Git-Tag: v8.4.0-rc2~354^2~3
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8d357ba88e5fb96c267a1162d0fce3e7542e5775;p=dealii.git

Add vertex to cells map.
---

diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h
index a3155cf56d..5daf06dd66 100644
--- a/include/deal.II/grid/grid_tools.h
+++ b/include/deal.II/grid/grid_tools.h
@@ -26,6 +26,7 @@
 
 #include <bitset>
 #include <list>
+#include <set>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -630,6 +631,15 @@ namespace GridTools
   compute_ghost_cell_halo_layer (const Container &container);
 
 
+  /**
+   * Return the adjacent cells of all the vertices. The vertices are ordered by
+   * their ID.
+   */
+  template <int dim, int spacedim>
+  std::vector<std::set<typename Triangulation<dim,spacedim>::active_cell_iterator> >
+  vertex_to_cell_map(const Triangulation<dim,spacedim> &triangulation);
+
+
   /*@}*/
   /**
    * @name Partitions and subdomains of triangulations
diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc
index dce6a7f6c7..828d0ccd92 100644
--- a/source/grid/grid_tools.cc
+++ b/source/grid/grid_tools.cc
@@ -1594,6 +1594,37 @@ next_cell:
 
 
 
+  template <int dim, int spacedim>
+  std::vector<std::set<typename Triangulation<dim,spacedim>::active_cell_iterator> >
+  vertex_to_cell_map(const Triangulation<dim,spacedim> &triangulation)
+  {
+    std::vector<std::set<typename Triangulation<dim,spacedim>::active_cell_iterator> >
+    vertex_to_cell_map(triangulation.n_vertices());
+    typename Triangulation<dim,spacedim>::active_cell_iterator cell = triangulation.begin_active(),
+                                                               endc = triangulation.end();
+    for (; cell!=endc; ++cell)
+      for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+        vertex_to_cell_map[cell->vertex_index(i)].insert(cell);
+
+    // Take care of hanging nodes
+    cell = triangulation.begin_active();
+    for (; cell!=endc; ++cell)
+      for (unsigned int i=0; i<GeometryInfo<dim>::faces_per_cell; ++i)
+        {
+          if ((cell->neighbor_level(i)>=0) && (cell->neighbor_level(i)<cell->level()))
+            {
+              typename Triangulation<dim,spacedim>::active_cell_iterator adjacent_cell =
+                cell->neighbor(i);
+              for (unsigned int j=0; j<GeometryInfo<dim>::vertices_per_face; ++j)
+                vertex_to_cell_map[cell->face(i)->vertex_index(j)].insert(adjacent_cell);
+            }
+        }
+
+    return vertex_to_cell_map;
+  }
+
+
+
   template <int dim, int spacedim>
   void
   get_face_connectivity_of_cells (const Triangulation<dim,spacedim> &triangulation,
diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in
index e6a85a1827..5e63a31814 100644
--- a/source/grid/grid_tools.inst.in
+++ b/source/grid/grid_tools.inst.in
@@ -234,6 +234,10 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
     std::map<unsigned int,Point<deal_II_space_dimension> >
     get_all_vertices_at_boundary (const Triangulation<deal_II_dimension,deal_II_space_dimension> &tria);
 
+    template
+    std::vector<std::set<typename 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);
+
 #if deal_II_dimension == deal_II_space_dimension
 #  if deal_II_dimension > 1
     template
diff --git a/tests/grid/grid_tools_vertex_to_cell_map.cc b/tests/grid/grid_tools_vertex_to_cell_map.cc
new file mode 100644
index 0000000000..3e66bf9902
--- /dev/null
+++ b/tests/grid/grid_tools_vertex_to_cell_map.cc
@@ -0,0 +1,93 @@
+//
+// Copyright (C) 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/distributed/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/grid/cell_id.h>
+
+#include <vector>
+
+void test()
+{
+  const MPI_Comm mpi_comm = MPI_COMM_WORLD;
+
+  parallel::distributed::Triangulation<2> tria(mpi_comm);
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(2);
+  unsigned int rank = Utilities::MPI::this_mpi_process(mpi_comm);
+  if (rank==0)
+    {
+      typename Triangulation<2>::active_cell_iterator cell = tria.begin_active(),
+                                                      end_cell = tria.end();
+      for (; cell!=end_cell; ++cell)
+        if (cell->is_locally_owned())
+          cell->set_refine_flag();
+    }
+
+  tria.execute_coarsening_and_refinement();
+
+  std::vector<std::set<typename Triangulation<2>::active_cell_iterator> >
+  vertex_to_cell = GridTools::vertex_to_cell_map(tria);
+
+  AssertThrow(tria.n_vertices()==vertex_to_cell.size(),
+              ExcMessage("Wrong number of vertices"));
+  std::vector<unsigned int> n_cells;
+  for (unsigned int i=0; i<vertex_to_cell.size(); ++i)
+    n_cells.push_back(vertex_to_cell[i].size());
+
+  if (rank==0)
+    {
+      std::vector<unsigned int> histogram(5,0);
+      for (unsigned int i=0; i<n_cells.size(); ++i)
+        histogram[n_cells[i]] += 1;
+
+      AssertThrow(histogram[0]==0, ExcMessage("Wrong cell distribution"));
+      AssertThrow(histogram[1]==4, ExcMessage("Wrong cell distribution"));
+      AssertThrow(histogram[2]==20, ExcMessage("Wrong cell distribution"));
+      AssertThrow(histogram[3]==4, ExcMessage("Wrong cell distribution"));
+      AssertThrow(histogram[4]==27, ExcMessage("Wrong cell distribution"));
+    }
+  if (rank==1)
+    {
+      std::vector<unsigned int> histogram(5,0);
+      for (unsigned int i=0; i<n_cells.size(); ++i)
+        histogram[n_cells[i]] += 1;
+
+      AssertThrow(histogram[0]==0, ExcMessage("Wrong cell distribution"));
+      AssertThrow(histogram[1]==4, ExcMessage("Wrong cell distribution"));
+      AssertThrow(histogram[2]==18, ExcMessage("Wrong cell distribution"));
+      AssertThrow(histogram[3]==6, ExcMessage("Wrong cell distribution"));
+      AssertThrow(histogram[4]==24, ExcMessage("Wrong cell distribution"));
+    }
+}
+
+
+int main (int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1);
+  MPILogInitAll log;
+  deallog.depth_console(0);
+
+  test();
+
+  deallog << "OK" << std::endl;
+
+  return 0;
+}
diff --git a/tests/grid/grid_tools_vertex_to_cell_map.mpirun=2.output b/tests/grid/grid_tools_vertex_to_cell_map.mpirun=2.output
new file mode 100644
index 0000000000..52103938b0
--- /dev/null
+++ b/tests/grid/grid_tools_vertex_to_cell_map.mpirun=2.output
@@ -0,0 +1,5 @@
+
+DEAL:0::OK
+
+DEAL:1::OK
+