From 193a3baf584ecf8979891367b70e50f23d4ee7ff Mon Sep 17 00:00:00 2001
From: kanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Wed, 28 Sep 2011 05:14:05 +0000
Subject: [PATCH] make patches in DoFTools

git-svn-id: https://svn.dealii.org/trunk@24456 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/include/deal.II/dofs/dof_tools.h |  40 +++++++
 deal.II/source/dofs/dof_tools.cc         | 137 +++++++++++++++++++++++
 deal.II/source/dofs/dof_tools.inst.in    |   9 +-
 3 files changed, 185 insertions(+), 1 deletion(-)

diff --git a/deal.II/include/deal.II/dofs/dof_tools.h b/deal.II/include/deal.II/dofs/dof_tools.h
index ddb98626ac..cd27f3b21d 100644
--- a/deal.II/include/deal.II/dofs/dof_tools.h
+++ b/deal.II/include/deal.II/dofs/dof_tools.h
@@ -19,6 +19,7 @@
 #include <deal.II/base/index_set.h>
 #include <deal.II/base/point.h>
 #include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/sparsity_pattern.h>
 #include <deal.II/dofs/function_map.h>
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/fe/fe.h>
@@ -221,11 +222,19 @@ namespace DoFTools
 				   /**
 				    * Maximal number of degrees of
 				    * freedom on a cell.
+				    *
+				    * @relates DoFHandler
 				    */
   template <int dim, int spacedim>
   unsigned int
   max_dofs_per_cell (const DoFHandler<dim,spacedim> &dh);
 
+				   /**
+				    * Maximal number of degrees of
+				    * freedom on a cell.
+				    *
+				    * @relates hp::DoFHandler
+				    */
   template <int dim, int spacedim>
   unsigned int
   max_dofs_per_cell (const hp::DoFHandler<dim,spacedim> &dh);
@@ -239,6 +248,8 @@ namespace DoFTools
 				    * and hp DoFHandlers, to allow for a
 				    * uniform interface to query this
 				    * property.
+				    *
+				    * @relates DoFHandler
 				    */
   template <int dim, int spacedim>
   unsigned int
@@ -252,6 +263,8 @@ namespace DoFTools
 				    * and hp DoFHandlers, to allow for a
 				    * uniform interface to query this
 				    * property.
+				    *
+				    * @relates hp::DoFHandler
 				    */
   template <int dim, int spacedim>
   unsigned int
@@ -265,6 +278,8 @@ namespace DoFTools
 				    * and hp DoFHandlers, to allow for a
 				    * uniform interface to query this
 				    * property.
+				    *
+				    * @relates DoFHandler
 				    */
   template <int dim, int spacedim>
   unsigned int
@@ -278,6 +293,8 @@ namespace DoFTools
 				    * and hp DoFHandlers, to allow for a
 				    * uniform interface to query this
 				    * property.
+				    *
+				    * @relates hp::DoFHandler
 				    */
   template <int dim, int spacedim>
   unsigned int
@@ -292,6 +309,8 @@ namespace DoFTools
 				    * and hp DoFHandlers, to allow for a
 				    * uniform interface to query this
 				    * property.
+				    *
+				    * @relates DoFHandler
 				    */
   template <int dim, int spacedim>
   unsigned int
@@ -306,6 +325,8 @@ namespace DoFTools
 				    * and hp DoFHandlers, to allow for a
 				    * uniform interface to query this
 				    * property.
+				    *
+				    * @relates hp::DoFHandler
 				    */
   template <int dim, int spacedim>
   unsigned int
@@ -320,6 +341,8 @@ namespace DoFTools
 				    * and hp DoFHandlers, to allow for a
 				    * uniform interface to query this
 				    * property.
+				    *
+				    * @relates DoFHandler
 				    */
   template <int dim, int spacedim>
   bool
@@ -334,6 +357,8 @@ namespace DoFTools
 				    * and hp DoFHandlers, to allow for a
 				    * uniform interface to query this
 				    * property.
+				    *
+				    * @relates hp::DoFHandler
 				    */
   template <int dim, int spacedim>
   bool
@@ -1455,6 +1480,21 @@ namespace DoFTools
 			 const DH& dof_handler,
 			 const unsigned int level);
   
+				    /**
+				     * Create a sparsity pattern which
+				     * in each row lists the degrees of
+				     * freedom associated to the
+				     * corresponding cells
+				     * corresponding to a vertex.
+				     */
+   template <class DH>
+   void make_vertex_patches(SparsityPattern& block_list,
+			    const DH& dof_handler,
+			    const unsigned int level,
+			    const bool interior_dofs_only,
+			    const bool boundary_patches = true,
+			    const bool single_cell_patches = false);
+  
 				   //@}
 				   /**
 				    * Extract a vector that represents the
diff --git a/deal.II/source/dofs/dof_tools.cc b/deal.II/source/dofs/dof_tools.cc
index de0e400b78..2d4eb64933 100644
--- a/deal.II/source/dofs/dof_tools.cc
+++ b/deal.II/source/dofs/dof_tools.cc
@@ -5859,9 +5859,146 @@ namespace DoFTools
 	    }
 	}
   }
+
+  template <class DH, class Sparsity>
+  void make_cell_patches(Sparsity& block_list, const DH& dof_handler, const unsigned int level)
+  {
+    typename DH::cell_iterator cell;
+    typename DH::cell_iterator endc = dof_handler.end(level);
+    std::vector<unsigned int> indices;
+    
+    unsigned int i=0;
+    for (cell=dof_handler.begin(level); cell != endc; ++i, ++cell)
+      {
+	indices.resize(cell->get_fe().dofs_per_cell);
+	cell->get_mg_dof_indices(indices);
+	for (unsigned int j=0;j<indices.size();++j)
+	  block_list.add(i,indices[j]);
+      }
+  }
+  
+  template <class DH>
+  void make_vertex_patches(
+    SparsityPattern& block_list,
+    const DH& dof_handler,
+    const unsigned int level,
+    const bool interior_only,
+    const bool boundary_patches,
+    const bool single_cell_patches)
+  {
+    typename DH::cell_iterator cell;
+    typename DH::cell_iterator endc = dof_handler.end(level);
+
+				     // Vector mapping from vertex
+				     // index in the triangulation to
+				     // consecutive block indices on
+				     // this level
+				     // The number of cells at a vertex
+    std::vector<unsigned int> vertex_cell_count(dof_handler.get_tria().n_vertices(), 0);
+				     // Is a vertex at the boundary?
+    std::vector<bool> vertex_boundary(dof_handler.get_tria().n_vertices(), false);
+    std::vector<unsigned int> vertex_mapping(dof_handler.get_tria().n_vertices(),
+					     numbers::invalid_unsigned_int);
+				     // Estimate for the number of
+				     // dofs at this point
+    std::vector<unsigned int> vertex_dof_count(dof_handler.get_tria().n_vertices(), 0);
+    
+				     // Identify all vertices active
+				     // on this level and remember
+				     // some data about them
+    for (cell=dof_handler.begin(level); cell != endc; ++cell)
+      for (unsigned int v=0;v<GeometryInfo<DH::dimension>::vertices_per_cell;++v)
+	{
+	  const unsigned int vg = cell->vertex_index(v);
+	  vertex_dof_count[vg] += cell->get_fe().dofs_per_cell;
+	  ++vertex_cell_count[vg];
+	  for (unsigned int d=0;d<DH::dimension;++d)
+	    if (cell->at_boundary(GeometryInfo<DH::dimension>::vertex_to_face[v][d]))
+	      vertex_boundary[vg] = true;
+	}
+				     // From now on, onlly vertices
+				     // with positive dof count are
+				     // "in".
+    
+				     // Remove vertices at boundaries
+				     // or in corners
+    for (unsigned int vg=0;vg<vertex_dof_count.size();++vg)
+      if ((!single_cell_patches && vertex_cell_count[vg] < 2)
+	  ||
+	  (!boundary_patches && vertex_boundary[vg]))
+	vertex_dof_count[vg] = 0;
+
+				     // Create a mapping from all
+				     // vertices to the ones used here
+    unsigned int i=0;
+    for (unsigned int vg=0;vg<vertex_mapping.size();++vg)
+      if (vertex_dof_count[vg] != 0)
+	vertex_mapping[vg] = i++;
+
+				     // Compactify dof count
+    for (unsigned int vg=0;vg<vertex_mapping.size();++vg)
+      if (vertex_dof_count[vg] != 0)
+	vertex_dof_count[vertex_mapping[vg]] = vertex_dof_count[vg];
+    
+				     // Now that we have all the data,
+				     // we reduce it to the part we
+				     // actually want
+    vertex_dof_count.resize(i);
+    
+				     // At this point, the list of
+				     // patches is ready. Now we enter
+				     // the dofs into the sparsity
+				     // pattern.
+    block_list.reinit(vertex_dof_count.size(), dof_handler.n_dofs(level), vertex_dof_count);
+    
+    std::vector<unsigned int> indices;
+    std::vector<bool> exclude;
+    
+    for (cell=dof_handler.begin(level); cell != endc; ++cell)
+      {
+	const FiniteElement<DH::dimension>& fe = cell->get_fe();
+	indices.resize(fe.dofs_per_cell);
+	cell->get_mg_dof_indices(indices);
+	
+	for (unsigned int v=0;v<GeometryInfo<DH::dimension>::vertices_per_cell;++v)
+	  {
+	    const unsigned int vg = cell->vertex_index(v);
+	    const unsigned int block = vertex_mapping[vg];
+	    if (block == numbers::invalid_unsigned_int)
+	      continue;
+	    
+	    if (interior_only)
+	      {
+						 // Exclude degrees of
+						 // freedom on faces
+						 // opposite to the vertex
+		exclude.resize(fe.dofs_per_cell);
+		std::fill(exclude.begin(), exclude.end(), false);
+		const unsigned int dpf = fe.dofs_per_face;
+		
+		for (unsigned int d=0;d<DH::dimension;++d)
+		  {
+		    const unsigned int a_face = GeometryInfo<DH::dimension>::vertex_to_face[v][d];
+		    const unsigned int face = GeometryInfo<DH::dimension>::opposite_face[a_face];
+		    for (unsigned int i=0;i<dpf;++i)
+		      exclude[fe.face_to_cell_index(i,face)] = true;
+		  }
+		for (unsigned int j=0;j<indices.size();++j)
+		  if (!exclude[j])
+		    block_list.add(block, indices[j]);
+	      }
+	    else
+	      {
+		for (unsigned int j=0;j<indices.size();++j)
+		  block_list.add(block, indices[j]);
+	      }
+	  }
+      }
+  }
 }
 
 
+
 // explicit instantiations
 
 #include "dof_tools.inst"
diff --git a/deal.II/source/dofs/dof_tools.inst.in b/deal.II/source/dofs/dof_tools.inst.in
index 28287e1cf9..ab1df6b2e5 100644
--- a/deal.II/source/dofs/dof_tools.inst.in
+++ b/deal.II/source/dofs/dof_tools.inst.in
@@ -113,7 +113,11 @@ for (SP : SPARSITY_PATTERNS; deal_II_dimension : DIMENSIONS)
     DoFTools::make_flux_sparsity_pattern<hp::DoFHandler<deal_II_dimension>,SP>
     (const hp::DoFHandler<deal_II_dimension> &dof,
      SP    &sparsity);
-
+     
+    template void
+      DoFTools::make_cell_patches<MGDoFHandler<deal_II_dimension>,SP>
+      (SP&, const MGDoFHandler<deal_II_dimension>&, unsigned int);
+    
 #if deal_II_dimension > 1
 
     template void
@@ -236,6 +240,9 @@ for (DH : DOFHANDLERS; deal_II_dimension : DIMENSIONS)
 
 for (deal_II_dimension : DIMENSIONS)
 {
+  template
+  void DoFTools::make_vertex_patches (SparsityPattern&, const MGDoFHandler<deal_II_dimension>&,
+  unsigned int, bool, bool, bool);
 
 #if deal_II_dimension < 3
 template
-- 
2.39.5