From e768e320a3c9e0414849f737dabaaa17d988beeb Mon Sep 17 00:00:00 2001
From: Tobias Leicht <tobias.leicht@dlr.de>
Date: Thu, 3 Apr 2008 13:56:08 +0000
Subject: [PATCH] SolutionTransfer now has an additional template parameter for
 the dof handler to be used. Instantiations are provided for DoFHandler (the
 default) and hp::DoFHandler.

git-svn-id: https://svn.dealii.org/trunk@15980 0785d39b-7218-0410-832d-ea1e28bc413d
---
 .../include/numerics/solution_transfer.h      |  11 +-
 .../source/numerics/solution_transfer.cc      | 173 +++++++++++-------
 deal.II/doc/news/changes.h                    |   7 +
 3 files changed, 119 insertions(+), 72 deletions(-)

diff --git a/deal.II/deal.II/include/numerics/solution_transfer.h b/deal.II/deal.II/include/numerics/solution_transfer.h
index 74aa73e26d..cccef5b82b 100644
--- a/deal.II/deal.II/include/numerics/solution_transfer.h
+++ b/deal.II/deal.II/include/numerics/solution_transfer.h
@@ -201,7 +201,7 @@ DEAL_II_NAMESPACE_OPEN
  * @ingroup numerics
  * @author Ralf Hartmann, 1999
  */
-template<int dim, typename number>
+template<int dim, typename number, class DH=DoFHandler<dim> >
 class SolutionTransfer
 {
   public:
@@ -210,7 +210,7 @@ class SolutionTransfer
 				      * Constructor, takes the current DoFHandler
 				      * as argument.
 				      */
-    SolutionTransfer(const DoFHandler<dim> &dof);
+    SolutionTransfer(const DH &dof);
 
     				     /**
 				      * Destructor
@@ -374,6 +374,11 @@ class SolutionTransfer
 				      */
     DeclException0(ExcVectorsDifferFromInVectors);
 
+				     /**
+				      * Exception
+				      */
+    DeclException0(ExcNumberOfDoFsPerCellHasChanged);
+
 				     /**
 				      * Exception
 				      */
@@ -388,7 +393,7 @@ class SolutionTransfer
 				      * Pointer to the degree of freedom handler
 				      * to work with.
 				      */
-    SmartPointer<const DoFHandler<dim> > dof_handler;
+    SmartPointer<const DH> dof_handler;
     
 				     /**
 				      * Stores the number of DoFs before the
diff --git a/deal.II/deal.II/source/numerics/solution_transfer.cc b/deal.II/deal.II/source/numerics/solution_transfer.cc
index 1510619d88..8ad91606f6 100644
--- a/deal.II/deal.II/source/numerics/solution_transfer.cc
+++ b/deal.II/deal.II/source/numerics/solution_transfer.cc
@@ -26,31 +26,31 @@ DEAL_II_NAMESPACE_OPEN
 
 
 
-template<int dim, typename number>
-SolutionTransfer<dim, number>::SolutionTransfer(const DoFHandler<dim> &dof):
+template<int dim, typename number, class DH>
+SolutionTransfer<dim, number, DH>::SolutionTransfer(const DH &dof):
 		dof_handler(&dof),
 		n_dofs_old(0),
 		prepared_for(none)
 {}
 
 
-template<int dim, typename number>
-SolutionTransfer<dim, number>::~SolutionTransfer()
+template<int dim, typename number, class DH>
+SolutionTransfer<dim, number, DH>::~SolutionTransfer()
 {
   clear ();
 }
 
 
-template<int dim, typename number>
-SolutionTransfer<dim, number>::Pointerstruct::Pointerstruct():
+template<int dim, typename number, class DH>
+SolutionTransfer<dim, number, DH>::Pointerstruct::Pointerstruct():
 		indices_ptr(0),
 		dof_values_ptr(0)
 {}
 
 
 
-template<int dim, typename number>
-void SolutionTransfer<dim, number>::clear ()
+template<int dim, typename number, class DH>
+void SolutionTransfer<dim, number, DH>::clear ()
 {
   indices_on_cell.clear();
   dof_values_on_cell.clear();
@@ -60,8 +60,8 @@ void SolutionTransfer<dim, number>::clear ()
 }
 
 
-template<int dim, typename number>
-void SolutionTransfer<dim, number>::prepare_for_pure_refinement()
+template<int dim, typename number, class DH>
+void SolutionTransfer<dim, number, DH>::prepare_for_pure_refinement()
 { 
   Assert(prepared_for!=pure_refinement, ExcAlreadyPrepForRef());
   Assert(prepared_for!=coarsening_and_refinement, 
@@ -70,40 +70,34 @@ void SolutionTransfer<dim, number>::prepare_for_pure_refinement()
   clear();
 
   const unsigned int n_active_cells = dof_handler->get_tria().n_active_cells();
-  const unsigned int dofs_per_cell  = dof_handler->get_fe().dofs_per_cell;
   n_dofs_old=dof_handler->n_dofs();
 
 				   // efficient reallocation of indices_on_cell
-  std::vector<std::vector<unsigned int> > (n_active_cells,
-					   std::vector<unsigned int> (dofs_per_cell))
+  std::vector<std::vector<unsigned int> > (n_active_cells)
     .swap(indices_on_cell);
 
-  typename DoFHandler<dim>::cell_iterator cell = dof_handler->begin(),
-					  endc = dof_handler->end();
+  typename DH::active_cell_iterator cell = dof_handler->begin_active(),
+				    endc = dof_handler->end();
 
-  for (unsigned int i=0; cell!=endc; ++cell) 
+  for (unsigned int i=0; cell!=endc; ++cell, ++i) 
     {
-      if (cell->active())
-	{
-					   // on each cell store the indices
-					   // of the dofs. after refining we
-					   // get the values on the children
-					   // by taking these indices, getting
-					   // the respective values out of
-					   // the data vectors and prolonging
-					   // them to the children
-	  cell->get_dof_indices(indices_on_cell[i]);
-	  cell_map[std::make_pair(cell->level(),cell->index())].indices_ptr=&indices_on_cell[i];
-	  ++i;	  
-	}
+      indices_on_cell[i].resize(cell->get_fe().dofs_per_cell);
+				       // on each cell store the indices of the
+				       // dofs. after refining we get the values
+				       // on the children by taking these
+				       // indices, getting the respective values
+				       // out of the data vectors and prolonging
+				       // them to the children
+      cell->get_dof_indices(indices_on_cell[i]);
+      cell_map[std::make_pair(cell->level(),cell->index())].indices_ptr=&indices_on_cell[i];
     }
   prepared_for=pure_refinement;
 }
 
 
-template<int dim, typename number>  
+template<int dim, typename number, class DH>  
 void
-SolutionTransfer<dim, number>::refine_interpolate(const Vector<number> &in,
+SolutionTransfer<dim, number, DH>::refine_interpolate(const Vector<number> &in,
 						  Vector<number>       &out) const
 {
   Assert(prepared_for==pure_refinement, ExcNotPrepared());
@@ -114,11 +108,10 @@ SolutionTransfer<dim, number>::refine_interpolate(const Vector<number> &in,
          ExcMessage ("Vectors cannot be used as input and output"
                      " at the same time!"));
 
-  unsigned int dofs_per_cell=dof_handler->get_fe().dofs_per_cell;  
-  Vector<number> local_values(dofs_per_cell);
+  Vector<number> local_values(0);
 
-  typename DoFHandler<dim>::cell_iterator cell = dof_handler->begin(),
-					  endc = dof_handler->end();
+  typename DH::cell_iterator cell = dof_handler->begin(),
+			     endc = dof_handler->end();
 
   typename std::map<std::pair<unsigned int, unsigned int>, Pointerstruct>::const_iterator
     pointerstruct,
@@ -136,6 +129,18 @@ SolutionTransfer<dim, number>::refine_interpolate(const Vector<number> &in,
 					 // which is both done by one
 					 // function
 	{
+	  const unsigned int dofs_per_cell=cell->get_fe().dofs_per_cell;
+	  local_values.reinit(dofs_per_cell, true); // fast reinit, i.e. the
+						    // entries are not set to 0.
+					   // make sure that the size of the
+					   // stored indices is the same as
+					   // dofs_per_cell. this is kind of a
+					   // test if we use the same fe in the
+					   // hp case. to really do that test we
+					   // would have to store the fe_index
+					   // of all cells
+	  Assert(dofs_per_cell==(*pointerstruct->second.indices_ptr).size(),
+		 ExcNumberOfDoFsPerCellHasChanged());
 	  for (unsigned int i=0; i<dofs_per_cell; ++i)
 	    local_values(i)=in((*pointerstruct->second.indices_ptr)[i]);
 	  cell->set_dof_values_by_interpolation(local_values, out);
@@ -144,8 +149,8 @@ SolutionTransfer<dim, number>::refine_interpolate(const Vector<number> &in,
 }
 
 
-template<int dim, typename number>  
-void SolutionTransfer<dim, number>::refine_interpolate (Vector<number> &vec) const
+template<int dim, typename number, class DH>  
+void SolutionTransfer<dim, number, DH>::refine_interpolate (Vector<number> &vec) const
 {
   Assert(vec.size()==n_dofs_old, ExcWrongVectorSize(vec.size(),n_dofs_old));
 
@@ -156,9 +161,9 @@ void SolutionTransfer<dim, number>::refine_interpolate (Vector<number> &vec) con
 }
 
 
-template<int dim, typename number>
+template<int dim, typename number, class DH>
 void
-SolutionTransfer<dim, number>::
+SolutionTransfer<dim, number, DH>::
 prepare_for_coarsening_and_refinement(const std::vector<Vector<number> > &all_in)
 {
   Assert(prepared_for!=pure_refinement, ExcAlreadyPrepForRef());
@@ -171,7 +176,6 @@ prepare_for_coarsening_and_refinement(const std::vector<Vector<number> > &all_in
   clear();
 
   const unsigned int n_active_cells = dof_handler->get_tria().n_active_cells();
-  const unsigned int dofs_per_cell  = dof_handler->get_fe().dofs_per_cell;
   n_dofs_old=dof_handler->n_dofs();
 
   for (unsigned int i=0; i<in_size; ++i)
@@ -185,7 +189,7 @@ prepare_for_coarsening_and_refinement(const std::vector<Vector<number> > &all_in
 				   // and that'll stay or be refined
   unsigned int n_cells_to_coarsen=0;
   unsigned int n_cells_to_stay_or_refine=0;
-  typename DoFHandler<dim>::active_cell_iterator
+  typename DH::active_cell_iterator
     act_cell = dof_handler->begin_active(),
     endc = dof_handler->end();
   for (; act_cell!=endc; ++act_cell) 
@@ -209,25 +213,25 @@ prepare_for_coarsening_and_refinement(const std::vector<Vector<number> > &all_in
                                    // the following arrays in an efficient
                                    // way, without copying much
   std::vector<std::vector<unsigned int> >
-    (n_cells_to_stay_or_refine,
-     std::vector<unsigned int> (dofs_per_cell))
+    (n_cells_to_stay_or_refine)
     .swap(indices_on_cell);
   
   std::vector<std::vector<Vector<number> > >
     (n_coarsen_fathers,
-     std::vector<Vector<number> > (in_size,
-				   Vector<number> (dofs_per_cell)))
+     std::vector<Vector<number> > (in_size))
     .swap(dof_values_on_cell);
   
 				   // we need counters for
 				   // the 'to_stay_or_refine' cells 'n_sr' and
 				   // the 'coarsen_fathers' cells 'n_cf',
   unsigned int n_sr=0, n_cf=0;
-  typename DoFHandler<dim>::cell_iterator cell = dof_handler->begin();  
+  typename DH::cell_iterator cell = dof_handler->begin();  
   for (; cell!=endc; ++cell) 
     {
       if (cell->active() && !cell->coarsen_flag_set())
 	{
+	  const unsigned int dofs_per_cell=cell->get_fe().dofs_per_cell;
+	  indices_on_cell[n_sr].resize(dofs_per_cell);
 					   // cell will not be coarsened,
 					   // so we get away by storing the
 					   // dof indices and later
@@ -245,7 +249,10 @@ prepare_for_coarsening_and_refinement(const std::vector<Vector<number> > &all_in
 	  for (unsigned int i=1; i<cell->n_children(); ++i)
 	    Assert(cell->child(i)->coarsen_flag_set(),
 		   ExcTriaPrepCoarseningNotCalledBefore());
-	      
+	  const unsigned int dofs_per_cell=cell->get_fe().dofs_per_cell;
+	  std::vector<Vector<number> >(in_size, Vector<number>(dofs_per_cell))
+	    .swap(dof_values_on_cell[n_cf]);
+      
 	  for (unsigned int j=0; j<in_size; ++j)
 	    {
 					       // store the data of each of
@@ -264,17 +271,17 @@ prepare_for_coarsening_and_refinement(const std::vector<Vector<number> > &all_in
 }
 
 
-template<int dim, typename number>
+template<int dim, typename number, class DH>
 void
-SolutionTransfer<dim, number>::prepare_for_coarsening_and_refinement(const Vector<number> &in)
+SolutionTransfer<dim, number, DH>::prepare_for_coarsening_and_refinement(const Vector<number> &in)
 {
   std::vector<Vector<number> > all_in=std::vector<Vector<number> >(1, in);
   prepare_for_coarsening_and_refinement(all_in);
 }
 
 
-template<int dim, typename number>
-void SolutionTransfer<dim, number>::
+template<int dim, typename number, class DH>
+void SolutionTransfer<dim, number, DH>::
 interpolate (const std::vector<Vector<number> > &all_in,
 	     std::vector<Vector<number> >       &all_out) const
 {
@@ -293,16 +300,15 @@ interpolate (const std::vector<Vector<number> > &all_in,
              ExcMessage ("Vectors cannot be used as input and output"
                          " at the same time!"));
 
-  const unsigned int dofs_per_cell=dof_handler->get_fe().dofs_per_cell;  
-  Vector<number> local_values(dofs_per_cell);
-  std::vector<unsigned int> dofs(dofs_per_cell);
+  Vector<number> local_values;
+  std::vector<unsigned int> dofs;
 
   typename std::map<std::pair<unsigned int, unsigned int>, Pointerstruct>::const_iterator
     pointerstruct,
     cell_map_end=cell_map.end();
 
-  typename DoFHandler<dim>::cell_iterator cell = dof_handler->begin(),
-					  endc = dof_handler->end();
+  typename DH::cell_iterator cell = dof_handler->begin(),
+			     endc = dof_handler->end();
   for (; cell!=endc; ++cell) 
     {
       pointerstruct=cell_map.find(std::make_pair(cell->level(),cell->index()));
@@ -315,22 +321,34 @@ interpolate (const std::vector<Vector<number> > &all_in,
 	  const std::vector<Vector<number> > * const valuesptr
 	    =pointerstruct->second.dof_values_ptr;
 
+	  const unsigned int dofs_per_cell=cell->get_fe().dofs_per_cell;
+	  
 					   // cell stayed or is
 					   // refined
 	  if (indexptr)
 	    {
 	      Assert (valuesptr == 0,
 		      ExcInternalError());
-	      
+					       // make sure that the size of the
+					       // stored indices is the same as
+					       // dofs_per_cell. this is kind of
+					       // a test if we use the same fe
+					       // in the hp case. to really do
+					       // that test we would have to
+					       // store the fe_index of all
+					       // cells
+	      Assert(dofs_per_cell==(*indexptr).size(),
+		     ExcNumberOfDoFsPerCellHasChanged());
 					       // get the values of
 					       // each of the input
 					       // data vectors on this
 					       // cell and prolong it
 					       // to its children
+	      local_values.reinit(dofs_per_cell, true);
 	      for (unsigned int j=0; j<size; ++j)
 		{
 		  for (unsigned int i=0; i<dofs_per_cell; ++i)
-		    local_values(i)=all_in[j](indexptr->operator[](i));
+		    local_values(i)=all_in[j]((*indexptr)[i]);
 		  cell->set_dof_values_by_interpolation(local_values,
 							all_out[j]);
 		}
@@ -342,7 +360,7 @@ interpolate (const std::vector<Vector<number> > &all_in,
 	      Assert (!cell->has_children(), ExcInternalError());
 	      Assert (indexptr == 0,
 		      ExcInternalError());
-
+	      dofs.resize(dofs_per_cell);
 					       // get the local
 					       // indices
 	      cell->get_dof_indices(dofs);
@@ -351,8 +369,23 @@ interpolate (const std::vector<Vector<number> > &all_in,
 					       // stored data to the
 					       // new vectors
 	      for (unsigned int j=0; j<size; ++j)
-		for (unsigned int i=0; i<dofs_per_cell; ++i)
-		  all_out[j](dofs[i])=valuesptr->operator[](j)(i);
+		{
+						   // make sure that the size of
+						   // the stored indices is the
+						   // same as
+						   // dofs_per_cell. this is
+						   // kind of a test if we use
+						   // the same fe in the hp
+						   // case. to really do that
+						   // test we would have to
+						   // store the fe_index of all
+						   // cells
+		  Assert(dofs_per_cell==(*valuesptr)[j].size(),
+			 ExcNumberOfDoFsPerCellHasChanged());
+
+		  for (unsigned int i=0; i<dofs_per_cell; ++i)
+		    all_out[j](dofs[i])=(*valuesptr)[j](i);
+		}
 	    }
 					   // undefined status
 	  else
@@ -363,8 +396,8 @@ interpolate (const std::vector<Vector<number> > &all_in,
 
 
 
-template<int dim, typename number>
-void SolutionTransfer<dim, number>::interpolate(const Vector<number> &in,
+template<int dim, typename number, class DH>
+void SolutionTransfer<dim, number, DH>::interpolate(const Vector<number> &in,
 						Vector<number>       &out) const
 {
   Assert (in.size()==n_dofs_old,
@@ -383,9 +416,9 @@ void SolutionTransfer<dim, number>::interpolate(const Vector<number> &in,
 
 
 
-template<int dim, typename number>
+template<int dim, typename number, class DH>
 unsigned int
-SolutionTransfer<dim, number>::memory_consumption () const
+SolutionTransfer<dim, number, DH>::memory_consumption () const
 {
 				   // at the moment we do not include the memory
 				   // consumption of the cell_map as we have no
@@ -400,17 +433,19 @@ SolutionTransfer<dim, number>::memory_consumption () const
 
 
 
-template<int dim, typename number>
+template<int dim, typename number, class DH>
 unsigned int
-SolutionTransfer<dim, number>::Pointerstruct::memory_consumption () const
+SolutionTransfer<dim, number, DH>::Pointerstruct::memory_consumption () const
 {
   return sizeof(*this);
 }
 
 
 
-template class SolutionTransfer<deal_II_dimension, float>;
-template class SolutionTransfer<deal_II_dimension, double>;
+template class SolutionTransfer<deal_II_dimension, float, DoFHandler<deal_II_dimension> >;
+template class SolutionTransfer<deal_II_dimension, double, DoFHandler<deal_II_dimension> >;
+template class SolutionTransfer<deal_II_dimension, float, hp::DoFHandler<deal_II_dimension> >;
+template class SolutionTransfer<deal_II_dimension, double, hp::DoFHandler<deal_II_dimension> >;
 
 /*----------------------------   solution_transfer.cc     ----------------------*/
 
diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index 8566a72f00..ec1fe27f9a 100644
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -290,6 +290,13 @@ constraints individually.
 
 <ol>
 
+  <li> <p>Extendend: SolutionTransfer has an additional template parameters DH
+  for the dof handler class to use, which defaults to DoFHandler. However, an
+  instantiation is also provided for hp::DoFHandler.
+  <br>
+  (Tobias Leicht, 2008/04/03)
+  </p></li>
+
   <li> <p>Extended: DataOut has an extended interface now which enables a user
   decision, which cells shall be written as curved cells using the provided
   mapping: none, only cells at the boundary, or all cells. The last choice can
-- 
2.39.5