Agglomerate using connected components (#130)

* Add utility to compute errors

* Agglomerate connected components

include/agglomeration_handler.h

@@ -49,6 +49,7 @@
 
 #include <agglomeration_iterator.h>
 #include <agglomerator.h>
+#include <utils.h>
 
 #include <fstream>
 
@@ -283,6 +284,14 @@ class AgglomerationHandler : public Subscriptor
   agglomeration_iterator
   define_agglomerate(const AgglomerationContainer &cells);
 
+  /**
+   * Same as above, but checking that every vector of cells is connected. If
+   * not, each connected component is agglomerated by calling the
+   * define_agglomerate() function defined above.
+   */
+  void
+  define_agglomerate_with_check(const AgglomerationContainer &cells);
+
   inline const Triangulation<dim, spacedim> &
   get_triangulation() const;
 

include/multigrid_amg.h

@@ -19,6 +19,8 @@
 
 #include <deal.II/base/mg_level_object.h>
 
+#include <deal.II/fe/fe_dgq.h>
+
 #include <deal.II/lac/la_parallel_vector.h>
 
 #include <deal.II/matrix_free/operators.h>

include/utils.h

@@ -56,6 +56,30 @@ namespace Utils
 
 
 
+  /**
+   * We describe a undirect graph through a vector of nodes and its adjacency
+  information
+   */
+  struct Graph
+  {
+    std::vector<types::global_cell_index>              nodes;
+    std::vector<std::vector<types::global_cell_index>> adjacency;
+
+    void
+    print_graph() const
+    {
+      std::cout << "Graph information" << std::endl;
+      for (const auto &node : nodes)
+        {
+          std::cout << "Neighbors for node " << node << std::endl;
+          for (const auto &neigh_node : adjacency[node])
+            std::cout << neigh_node << std::endl;
+        }
+    }
+  };
+
+
+
   /**
    * Given a coarse AgglomerationHandler @p coarse_ah and a fine
    * AgglomerationHandler @p fine_ah, this function fills the injection matrix
@@ -1054,6 +1078,96 @@ namespace Utils
 
 
 
+  /**
+   * Helper function to compute the position of index @p index in vector @p v.
+   */
+  inline types::global_cell_index
+  get_index(const std::vector<types::global_cell_index> &v,
+            const types::global_cell_index               index)
+  {
+    return std::distance(v.begin(), std::find(v.begin(), v.end(), index));
+  }
+
+
+
+  template <int dim>
+  void
+  create_graph_from_agglomerate(
+    const std::vector<typename Triangulation<dim>::active_cell_iterator> &cells,
+    Graph                                                                &g)
+  {
+    Assert(cells.size() > 0, ExcMessage("No cells to be agglomerated."));
+    const unsigned int n_faces = cells[0]->n_faces();
+
+    std::vector<types::global_cell_index> vec_cells(cells.size());
+    for (size_t i = 0; i < cells.size(); i++)
+      vec_cells[i] = cells[i]->active_cell_index();
+
+    g.adjacency.resize(cells.size());
+    for (const auto &cell : cells)
+      {
+        // std::cout << "Cell idx: " << cell->active_cell_index() << std::endl;
+        // std::cout << "new idx: "
+        //           << get_index(vec_cells, cell->active_cell_index())
+        //           << std::endl;
+        g.nodes.push_back(get_index(vec_cells, cell->active_cell_index()));
+        for (unsigned int f = 0; f < n_faces; ++f)
+          {
+            const auto &neigh = cell->neighbor(f);
+            if (neigh.state() == IteratorState::IteratorStates::valid &&
+                std::find(cells.begin(), cells.end(), neigh) != std::end(cells))
+              g.adjacency[get_index(vec_cells, cell->active_cell_index())]
+                .push_back(get_index(vec_cells, neigh->active_cell_index()));
+          }
+      }
+  }
+
+
+
+  inline void
+  dfs(std::vector<types::global_cell_index> &comp,
+      std::vector<bool>                     &visited,
+      const Graph                           &g,
+      const types::global_cell_index         v)
+  {
+    visited[v] = true;
+    comp.push_back(v);
+    for (const types::global_cell_index u : g.adjacency[v])
+      {
+        if (!visited[u])
+          dfs(comp, visited, g, u);
+      }
+  }
+
+
+
+  void
+  compute_connected_components(
+    Graph                                              &g,
+    std::vector<std::vector<types::global_cell_index>> &connected_components)
+  {
+    Assert(g.nodes.size() > 0, ExcMessage("No nodes in this graph."));
+    Assert(
+      connected_components.size() == 0,
+      ExcMessage(
+        "Connected components have to be computed by the present function."));
+
+    std::vector<bool> visited(g.nodes.size()); // register visited node
+    std::fill(visited.begin(), visited.end(), 0);
+
+
+    for (types::global_cell_index v : g.nodes)
+      {
+        if (!visited[v])
+          {
+            connected_components.emplace_back();
+            dfs(connected_components.back(), visited, g, v);
+          }
+      }
+  }
+
+
+
 } // namespace Utils
 
 

source/agglomeration_handler.cc

@@ -104,6 +104,45 @@ AgglomerationHandler<dim, spacedim>::define_agglomerate(
 
 
 
+template <int dim, int spacedim>
+void
+AgglomerationHandler<dim, spacedim>::define_agglomerate_with_check(
+  const AgglomerationContainer &cells)
+{
+  Assert(cells.size() > 0, ExcMessage("No cells to be agglomerated."));
+
+  // Compute the graph associated to the agglomerate
+  Utils::Graph g;
+  Utils::create_graph_from_agglomerate<dim>(cells, g);
+
+  // The following vector will be filled with connected components
+  std::vector<std::vector<types::global_cell_index>> connected_components;
+  Utils::compute_connected_components(g, connected_components);
+
+  // Be verbose in debug mode
+#ifdef AGGLO_DEBUG
+  if (connected_components.size() > 1)
+    std::cout << "Disconnected elements. Connected components will be "
+                 "computed and agglomerated together."
+              << std::endl;
+#endif
+
+  // Get connected components and define one agglomerate for each one of them.
+  std::vector<typename Triangulation<dim>::active_cell_iterator> agglomerate;
+  for (const auto &comp : connected_components)
+    {
+      agglomerate.reserve(comp.size());
+      for (const types::global_cell_index local_idx : comp)
+        agglomerate.push_back(cells[local_idx]);
+
+      // Perform agglomeration
+      define_agglomerate(agglomerate);
+      agglomerate.clear();
+    }
+}
+
+
+
 template <int dim, int spacedim>
 void
 AgglomerationHandler<dim, spacedim>::initialize_fe_values(
@@ -547,8 +586,8 @@ AgglomerationHandler<dim, spacedim>::reinit(
 
   // First check if the polytope is made just by a single cell. If so, use
   // classical FEValues
-  if (polytope->n_background_cells() == 1)
-    return standard_scratch->reinit(deal_cell);
+  // if (polytope->n_background_cells() == 1)
+  //   return standard_scratch->reinit(deal_cell);
 
   const auto &agglo_cells = polytope->get_agglomerate();
 
@@ -972,9 +1011,13 @@ namespace dealii
                                                   &master_cell,
         const AgglomerationHandler<dim, spacedim> &handler)
       {
-        Assert(handler.master_slave_relationships
-                   [master_cell->global_active_cell_index()] == -1,
-               ExcMessage("The present cell is not a master one."));
+        Assert(
+          handler.master_slave_relationships[master_cell
+                                               ->global_active_cell_index()] ==
+            -1,
+          ExcMessage("The present cell with index " +
+                     std::to_string(master_cell->global_active_cell_index()) +
+                     "is not a master one."));
 
         const auto &agglomeration = handler.get_agglomerate(master_cell);
         const types::global_cell_index current_polytope_index =

source/multigrid_amg.cc

@@ -10,8 +10,8 @@
 //
 // -----------------------------------------------------------------------------
 
+#include <deal.II/lac/trilinos_sparse_matrix.h>
 
-#include <agglomeration_handler.h>
 #include <multigrid_amg.h>
 
 namespace dealii