From b41b6cee9dfacdbd17f2f4c396738ceb1bfd8d62 Mon Sep 17 00:00:00 2001
From: heister <heister@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Sat, 4 Jan 2014 23:24:03 +0000
Subject: [PATCH] rework MGConstrainedDofs to use IndexSet as the datastructure

git-svn-id: https://svn.dealii.org/trunk@32167 0785d39b-7218-0410-832d-ea1e28bc413d
---
 .../deal.II/multigrid/mg_constrained_dofs.h   | 104 ++++++++++++++----
 deal.II/include/deal.II/multigrid/mg_tools.h  |  15 ++-
 deal.II/source/multigrid/mg_tools.cc          |  56 ++++++----
 deal.II/source/multigrid/mg_tools.inst.in     |   9 +-
 4 files changed, 133 insertions(+), 51 deletions(-)

diff --git a/deal.II/include/deal.II/multigrid/mg_constrained_dofs.h b/deal.II/include/deal.II/multigrid/mg_constrained_dofs.h
index bf736c9883..917dc33771 100644
--- a/deal.II/include/deal.II/multigrid/mg_constrained_dofs.h
+++ b/deal.II/include/deal.II/multigrid/mg_constrained_dofs.h
@@ -106,7 +106,6 @@ public:
    * Return the indices of dofs for each level that lie on the
    * boundary of the domain.
    */
-// TODO: remove
   const std::vector<std::set<types::global_dof_index> > &
   get_boundary_indices () const;
 
@@ -115,23 +114,35 @@ public:
    * refinement edge (i.e. are on faces between cells of this level
    * and cells on the level below).
    */
-  // TODO: remove
   const std::vector<std::vector<bool> > &
-  get_refinement_edge_indices () const;
+  get_refinement_edge_indices () const DEAL_II_DEPRECATED;
+
+  /**
+   * Return the indices of dofs on the given level that lie on an
+   * refinement edge (dofs on faces to neighbors that are coarser)
+   */
+  const IndexSet &
+  get_refinement_edge_indices (unsigned int level) const;
 
   /**
    * Return the indices of dofs for each level that are in the
    * intersection of the sets returned by get_boundary_indices() and
    * get_refinement_edge_indices().
    */
-  // TODO: remove
   const std::vector<std::vector<bool> > &
-  get_refinement_edge_boundary_indices () const;
+  get_refinement_edge_boundary_indices () const DEAL_II_DEPRECATED;
 
   /**
-   * Return if boundary_indices need to be set or not.
+   * Return indices of all dofs that are on boundary faces on the given level 
+   * if the cell has refinement edge indices (i.e. has a coarser neighbor).
    */
+  const IndexSet &
+  get_refinement_edge_boundary_indices (unsigned int level) const;
 
+  /**
+   * Return if boundary_indices need to be set or not.
+   */
+//TODO: remove?
   bool set_boundary_values () const;
 
 private:
@@ -145,7 +156,7 @@ private:
    * The degrees of freedom on the refinement edge between a level and
    * coarser cells.
    */
-  std::vector<std::vector<bool> > refinement_edge_indices;
+  std::vector<IndexSet> refinement_edge_indices;
 
   /**
    * The degrees of freedom on the refinement edge between a level and
@@ -153,7 +164,17 @@ private:
    *
    * This is a subset of #refinement_edge_indices.
    */
-  std::vector<std::vector<bool> > refinement_edge_boundary_indices;
+  std::vector<IndexSet> refinement_edge_boundary_indices;
+
+  /**
+   * old data structure only filled on demand
+   */
+  mutable std::vector<std::vector<bool> > refinement_edge_boundary_indices_old;
+
+  /**
+   * old data structure only filled on demand
+   */
+  mutable std::vector<std::vector<bool> > refinement_edge_indices_old;
 };
 
 
@@ -165,10 +186,12 @@ MGConstrainedDoFs::initialize(const DoFHandler<dim,spacedim> &dof)
   const unsigned int nlevels = dof.get_tria().n_global_levels();
   refinement_edge_indices.resize(nlevels);
   refinement_edge_boundary_indices.resize(nlevels);
+  refinement_edge_indices_old.clear();
+  refinement_edge_boundary_indices_old.clear();
   for (unsigned int l=0; l<nlevels; ++l)
     {
-      refinement_edge_indices[l].resize(dof.n_dofs(l));
-      refinement_edge_boundary_indices[l].resize(dof.n_dofs(l));
+      refinement_edge_indices[l] = IndexSet(dof.n_dofs(l));
+      refinement_edge_boundary_indices[l] = IndexSet(dof.n_dofs(l));
     }
   MGTools::extract_inner_interface_dofs (dof, refinement_edge_indices,
                                          refinement_edge_boundary_indices);
@@ -187,12 +210,14 @@ MGConstrainedDoFs::initialize(
   boundary_indices.resize(nlevels);
   refinement_edge_indices.resize(nlevels);
   refinement_edge_boundary_indices.resize(nlevels);
+  refinement_edge_indices_old.clear();
+  refinement_edge_boundary_indices_old.clear();
 
   for (unsigned int l=0; l<nlevels; ++l)
     {
       boundary_indices[l].clear();
-      refinement_edge_indices[l].resize(dof.n_dofs(l));
-      refinement_edge_boundary_indices[l].resize(dof.n_dofs(l));
+      refinement_edge_indices[l] = IndexSet(dof.n_dofs(l));
+      refinement_edge_boundary_indices[l] = IndexSet(dof.n_dofs(l));
     }
 
   MGTools::make_boundary_list (dof, function_map, boundary_indices, component_mask);
@@ -213,6 +238,9 @@ MGConstrainedDoFs::clear()
 
   for (unsigned int l=0; l<refinement_edge_boundary_indices.size(); ++l)
     refinement_edge_boundary_indices[l].clear();
+
+  refinement_edge_indices_old.clear();
+  refinement_edge_boundary_indices_old.clear();
 }
 
 
@@ -225,10 +253,7 @@ MGConstrainedDoFs::is_boundary_index (const unsigned int level,
     return false;
   
   AssertIndexRange(level, boundary_indices.size());
-  if (boundary_indices[level].find(index) != boundary_indices[level].end())
-    return true;
-  else
-    return false;
+  return (boundary_indices[level].find(index) != boundary_indices[level].end());
 }
 
 inline
@@ -245,9 +270,8 @@ MGConstrainedDoFs::at_refinement_edge (const unsigned int level,
                                        const types::global_dof_index index) const
 {
   AssertIndexRange(level, refinement_edge_indices.size());
-  AssertIndexRange(index, refinement_edge_indices[level].size());
 
-  return refinement_edge_indices[level][index];
+  return refinement_edge_indices[level].is_element(index);
 }
 
 
@@ -257,9 +281,8 @@ MGConstrainedDoFs::at_refinement_edge_boundary (const unsigned int level,
                                                 const types::global_dof_index index) const
 {
   AssertIndexRange(level, refinement_edge_boundary_indices.size());
-  AssertIndexRange(index, refinement_edge_boundary_indices[level].size());
 
-  return refinement_edge_boundary_indices[level][index];
+  return refinement_edge_boundary_indices[level].is_element(index);
 }
 
 inline
@@ -273,14 +296,51 @@ inline
 const std::vector<std::vector<bool> > &
 MGConstrainedDoFs::get_refinement_edge_indices () const
 {
-  return refinement_edge_indices;
+  if (refinement_edge_indices_old.size()!=refinement_edge_indices.size())
+    {
+      unsigned int n_levels = refinement_edge_indices.size();
+      refinement_edge_indices_old.resize(n_levels);
+      for (unsigned int l=0;l<n_levels;++l)
+        {
+          refinement_edge_indices_old[l].resize(refinement_edge_indices[l].size(), false);
+          refinement_edge_indices[l].fill_binary_vector(refinement_edge_indices_old[l]);
+        }
+    }
+
+  return refinement_edge_indices_old;
+}
+
+inline
+const IndexSet &
+MGConstrainedDoFs::get_refinement_edge_indices (unsigned int level) const
+{
+  AssertIndexRange(level, refinement_edge_indices.size());
+  return refinement_edge_indices[level];
 }
 
 inline
 const std::vector<std::vector<bool> > &
 MGConstrainedDoFs::get_refinement_edge_boundary_indices () const
 {
-  return refinement_edge_boundary_indices;
+  if (refinement_edge_boundary_indices_old.size()!=refinement_edge_boundary_indices.size())
+    {
+      unsigned int n_levels = refinement_edge_boundary_indices.size();
+      refinement_edge_boundary_indices_old.resize(n_levels);
+      for (unsigned int l=0;l<n_levels;++l)
+        {
+          refinement_edge_boundary_indices_old[l].resize(refinement_edge_boundary_indices[l].size(), false);
+          refinement_edge_boundary_indices[l].fill_binary_vector(refinement_edge_boundary_indices_old[l]);
+        }
+    }
+
+  return refinement_edge_boundary_indices_old;
+}
+
+inline
+const IndexSet &
+MGConstrainedDoFs::get_refinement_edge_boundary_indices (unsigned int level) const
+{
+  return refinement_edge_boundary_indices[level];
 }
 
 inline
diff --git a/deal.II/include/deal.II/multigrid/mg_tools.h b/deal.II/include/deal.II/multigrid/mg_tools.h
index e06cbc96d9..bb2b9b08d1 100644
--- a/deal.II/include/deal.II/multigrid/mg_tools.h
+++ b/deal.II/include/deal.II/multigrid/mg_tools.h
@@ -292,8 +292,8 @@ namespace MGTools
   template <int dim, int spacedim>
   void
   extract_inner_interface_dofs (const DoFHandler<dim,spacedim> &mg_dof_handler,
-                                std::vector<std::vector<bool> >  &interface_dofs,
-                                std::vector<std::vector<bool> >  &boundary_interface_dofs);
+                                std::vector<IndexSet>  &interface_dofs,
+                                std::vector<IndexSet>  &boundary_interface_dofs);
 
   /**
    * Does the same as the function above,
@@ -302,7 +302,16 @@ namespace MGTools
   template <int dim, int spacedim>
   void
   extract_inner_interface_dofs (const DoFHandler<dim,spacedim> &mg_dof_handler,
-                                std::vector<std::vector<bool> >  &interface_dofs);
+                                std::vector<IndexSet>  &interface_dofs);
+
+  /**
+   * As above but with a deprecated data structure. This makes one additional copy.
+   */
+  template <int dim, int spacedim>
+  void
+  extract_inner_interface_dofs (const DoFHandler<dim,spacedim> &mg_dof_handler,
+                                std::vector<std::vector<bool> >  &interface_dofs) DEAL_II_DEPRECATED;
+
 
   template <int dim, int spacedim>
   void
diff --git a/deal.II/source/multigrid/mg_tools.cc b/deal.II/source/multigrid/mg_tools.cc
index 0bf94ebb10..85bf569a4f 100644
--- a/deal.II/source/multigrid/mg_tools.cc
+++ b/deal.II/source/multigrid/mg_tools.cc
@@ -1416,7 +1416,7 @@ namespace MGTools
   template <int dim, int spacedim>
   void
   extract_inner_interface_dofs (const DoFHandler<dim,spacedim> &mg_dof_handler,
-                                std::vector<std::vector<bool> >  &interface_dofs)
+                                std::vector<IndexSet>  &interface_dofs)
   {
     Assert (interface_dofs.size() == mg_dof_handler.get_tria().n_global_levels(),
             ExcDimensionMismatch (interface_dofs.size(),
@@ -1427,10 +1427,7 @@ namespace MGTools
         Assert (interface_dofs[l].size() == mg_dof_handler.n_dofs(l),
                 ExcDimensionMismatch (interface_dofs[l].size(),
                                       mg_dof_handler.n_dofs(l)));
-
-        std::fill (interface_dofs[l].begin(),
-                   interface_dofs[l].end(),
-                   false);
+       interface_dofs[l].clear();
       }
 
     const FiniteElement<dim,spacedim> &fe = mg_dof_handler.get_fe();
@@ -1475,7 +1472,29 @@ namespace MGTools
 
         for (unsigned int i=0; i<dofs_per_cell; ++i)
           if (cell_dofs[i])
-            interface_dofs[level][local_dof_indices[i]] = true;
+            interface_dofs[level].add_index(local_dof_indices[i]);
+      }
+  }
+
+  template <int dim, int spacedim>
+  void
+  extract_inner_interface_dofs (const DoFHandler<dim,spacedim> &mg_dof_handler,
+                                std::vector<std::vector<bool> >  &interface_dofs)
+  {
+    std::vector<IndexSet> temp;
+    temp.resize(interface_dofs.size());
+    for (unsigned int l=0;l<interface_dofs.size();++l)
+      temp[l] = IndexSet(interface_dofs[l].size());
+
+    extract_inner_interface_dofs(mg_dof_handler, temp);
+
+    for (unsigned int l=0;l<interface_dofs.size();++l)
+      {
+        Assert (interface_dofs[l].size() == mg_dof_handler.n_dofs(l),
+                ExcDimensionMismatch (interface_dofs[l].size(),
+                                      mg_dof_handler.n_dofs(l)));
+
+        temp[l].fill_binary_vector(interface_dofs[l]);
       }
   }
 
@@ -1483,7 +1502,7 @@ namespace MGTools
   template <int dim, int spacedim>
   void
   extract_non_interface_dofs (const DoFHandler<dim,spacedim> &mg_dof_handler,
-                              std::vector<std::set<types::global_dof_index> >  &non_interface_dofs)
+                              std::vector<std::set<unsigned int> >  &non_interface_dofs)
   {
     Assert (non_interface_dofs.size() == mg_dof_handler.get_tria().n_global_levels(),
             ExcDimensionMismatch (non_interface_dofs.size(),
@@ -1552,8 +1571,8 @@ namespace MGTools
   template <int dim, int spacedim>
   void
   extract_inner_interface_dofs (const DoFHandler<dim,spacedim> &mg_dof_handler,
-                                std::vector<std::vector<bool> >  &interface_dofs,
-                                std::vector<std::vector<bool> >  &boundary_interface_dofs)
+                                std::vector<IndexSet>  &interface_dofs,
+                                std::vector<IndexSet>  &boundary_interface_dofs)
   {
     Assert (interface_dofs.size() == mg_dof_handler.get_tria().n_global_levels(),
             ExcDimensionMismatch (interface_dofs.size(),
@@ -1564,19 +1583,8 @@ namespace MGTools
 
     for (unsigned int l=0; l<mg_dof_handler.get_tria().n_global_levels(); ++l)
       {
-        Assert (interface_dofs[l].size() == mg_dof_handler.n_dofs(l),
-                ExcDimensionMismatch (interface_dofs[l].size(),
-                                      mg_dof_handler.n_dofs(l)));
-        Assert (boundary_interface_dofs[l].size() == mg_dof_handler.n_dofs(l),
-                ExcDimensionMismatch (boundary_interface_dofs[l].size(),
-                                      mg_dof_handler.n_dofs(l)));
-
-        std::fill (interface_dofs[l].begin(),
-                   interface_dofs[l].end(),
-                   false);
-        std::fill (boundary_interface_dofs[l].begin(),
-                   boundary_interface_dofs[l].end(),
-                   false);
+        interface_dofs[l].clear();
+        boundary_interface_dofs[l].clear();
       }
 
     const FiniteElement<dim,spacedim> &fe = mg_dof_handler.get_fe();
@@ -1639,10 +1647,10 @@ namespace MGTools
         for (unsigned int i=0; i<dofs_per_cell; ++i)
           {
             if (cell_dofs[i])
-              interface_dofs[level][local_dof_indices[i]] = true;
+              interface_dofs[level].add_index(local_dof_indices[i]);
 
             if (boundary_cell_dofs[i])
-              boundary_interface_dofs[level][local_dof_indices[i]] = true;
+              boundary_interface_dofs[level].add_index(local_dof_indices[i]);
           }
       }
   }
diff --git a/deal.II/source/multigrid/mg_tools.inst.in b/deal.II/source/multigrid/mg_tools.inst.in
index 8d4a68359b..7f25f462f0 100644
--- a/deal.II/source/multigrid/mg_tools.inst.in
+++ b/deal.II/source/multigrid/mg_tools.inst.in
@@ -118,12 +118,17 @@ for (deal_II_dimension : DIMENSIONS)
     template
 	void
 	extract_inner_interface_dofs (const DoFHandler<deal_II_dimension> &mg_dof_handler,
-				      std::vector<std::vector<bool> >  &interface_dofs,
-				      std::vector<std::vector<bool> >  &boundary_interface_dofs);
+				      std::vector<IndexSet>  &interface_dofs,
+				      std::vector<IndexSet>  &boundary_interface_dofs);
+      template
+	void
+	extract_inner_interface_dofs (const DoFHandler<deal_II_dimension> &mg_dof_handler,
+				      std::vector<IndexSet>  &interface_dofs);
       template
 	void
 	extract_inner_interface_dofs (const DoFHandler<deal_II_dimension> &mg_dof_handler,
 				      std::vector<std::vector<bool> >  &interface_dofs);
+
       template
         void
         extract_non_interface_dofs (const DoFHandler<deal_II_dimension> & mg_dof_handler,
-- 
2.39.5