From dafb624c8d169d02fe4942cbdb2b79f93ddcdc21 Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Fri, 14 Apr 2000 17:09:32 +0000
Subject: [PATCH] Move the generation of intergrid constraints to the library.

git-svn-id: https://svn.dealii.org/trunk@2728 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/deal.II/include/dofs/dof_tools.h |  40 ++-
 deal.II/deal.II/source/dofs/dof_tools.cc | 437 +++++++++++++++++++++++
 2 files changed, 466 insertions(+), 11 deletions(-)

diff --git a/deal.II/deal.II/include/dofs/dof_tools.h b/deal.II/deal.II/include/dofs/dof_tools.h
index 5bcb95a0a4..985cf0b36d 100644
--- a/deal.II/deal.II/include/dofs/dof_tools.h
+++ b/deal.II/deal.II/include/dofs/dof_tools.h
@@ -248,8 +248,9 @@ class DoFTools
 				      * considered.
 				      */
     template<int dim>
-    static void make_flux_sparsity_pattern (const DoFHandler<dim> &dof_handler,
-					    SparsityPattern       &sparsity_pattern);
+    static void
+    make_flux_sparsity_pattern (const DoFHandler<dim> &dof_handler,
+				SparsityPattern       &sparsity_pattern);
     
 				     /**
 				      * Make up the constraints which
@@ -353,22 +354,35 @@ class DoFTools
 				      * the number of components in the
 				      * finite element used by #dof#.
 				      */
-    template<int dim>
-    static void extract_dofs(const DoFHandler<dim> &dof_handler,
-			     const vector<bool>    &select,
-			     vector<bool>          &selected_dofs);
+    template <int dim>
+    static void
+    extract_dofs(const DoFHandler<dim> &dof_handler,
+		 const vector<bool>    &select,
+		 vector<bool>          &selected_dofs);
 
 				     /**
 				      * Do the same thing as #extract_dofs#
 				      * for one level of a multi-grid DoF
 				      * numbering.
 				      */
-    template<int dim>
-    static void extract_level_dofs(const unsigned int       level,
-				   const MGDoFHandler<dim> &dof,
-				   const vector<bool>      &select,
-				   vector<bool>            &selected_dofs);
+    template <int dim>
+    static void
+    extract_level_dofs(const unsigned int       level,
+		       const MGDoFHandler<dim> &dof,
+		       const vector<bool>      &select,
+		       vector<bool>            &selected_dofs);
 
+
+    template <int dim>
+    static void
+    compute_intergrid_constraints (const DoFHandler<dim>              &coarse_grid,
+				   const unsigned int                  coarse_component,
+				   const DoFHandler<dim>              &fine_grid,
+				   const unsigned int                  fine_component,
+				   const InterGridMap<DoFHandler,dim> &coarse_to_fine_grid_map,
+				   ConstraintMatrix                   &constraints);
+
+				   
 				     /**
 				      * Exception
 				      */
@@ -385,6 +399,10 @@ class DoFTools
 		    << "is invalid with respect to the number "
 		    << "of components in the finite element "
 		    << "(" << arg2 << ")");
+				     /**
+				      * Exception
+				      */
+    DeclException0 (ExcFiniteElementsDontMatch);
 };
 
 
diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc
index e7a841a06d..1a7cceb263 100644
--- a/deal.II/deal.II/source/dofs/dof_tools.cc
+++ b/deal.II/deal.II/source/dofs/dof_tools.cc
@@ -14,6 +14,7 @@
 
 #include <grid/tria.h>
 #include <grid/tria_iterator.h>
+#include <grid/intergrid_map.h>
 #include <dofs/dof_handler.h>
 #include <dofs/dof_accessor.h>
 #include <dofs/dof_constraints.h>
@@ -709,6 +710,433 @@ DoFTools::extract_level_dofs(const unsigned int       level,
 }
 
 
+
+
+template <int dim>
+void
+DoFTools::compute_intergrid_constraints (const DoFHandler<dim>              &coarse_grid,
+					 const unsigned int                  coarse_component,
+					 const DoFHandler<dim>              &fine_grid,
+					 const unsigned int                  fine_component,
+					 const InterGridMap<DoFHandler,dim> &coarse_to_fine_grid_map,
+					 ConstraintMatrix                   &constraints)
+{
+				   // aliases to the finite elements
+				   // used by the dof handlers:
+  const FiniteElement<dim> &coarse_fe = coarse_grid.get_fe(),
+			   &fine_fe   = fine_grid.get_fe();
+
+				   // global numbers of dofs
+  const unsigned int n_coarse_dofs = coarse_grid.n_dofs(),
+		     n_fine_dofs   = fine_grid.n_dofs();
+
+				   // local numbers of dofs
+  const unsigned int coarse_dofs_per_cell = coarse_fe.dofs_per_cell,
+		     fine_dofs_per_cell   = fine_fe.dofs_per_cell;
+
+				   // alias the number of dofs per
+				   // cell belonging to the
+				   // coarse_component which is to be
+				   // the restriction of the fine
+				   // grid:
+  const unsigned int coarse_dofs_per_cell_component
+    = coarse_fe.base_element(coarse_fe.component_to_base(coarse_component)).dofs_per_cell;
+  
+  
+				   // check whether component numbers
+				   // are valid
+  Assert (coarse_component < coarse_fe.n_components(),
+	  ExcInvalidComponent (coarse_component, coarse_fe.n_components()));
+  Assert (fine_component < fine_fe.n_components(),
+	  ExcInvalidComponent (fine_component, fine_fe.n_components()));
+				   // check whether respective finite
+				   // elements are equal
+  Assert (coarse_fe.base_element (coarse_fe.component_to_base(coarse_component))
+	  ==
+	  fine_fe.base_element (fine_fe.component_to_base(fine_component)),
+	  ExcFiniteElementsDontMatch());
+  
+
+/*
+ * From here on: the term `parameter' refers to the selected component
+ * on the coarse grid and its analogon on the fine grid. The naming of
+ * variables containing this term is due to the fact that
+ * `selected_component' is longer, but also due to the fact that the
+ * code of this function was initially written for a program where the
+ * component which we wanted to match between grids was actually the
+ * `parameter' variable.
+ *
+ * Likewise, the terms `parameter grid' and `state grid' refer to the
+ * coarse and fine grids, respectively.
+ *
+ * Changing the names of variables would in principle be a good idea,
+ * but would not make things simpler and would be another source of
+ * errors. If anyone feels like doing so: patches would be welcome!
+ */
+
+
+  
+				   // set up vectors of cell-local
+				   // data; each vector represents one
+				   // degree of freedom of the
+				   // coarse-grid variable in the
+				   // fine-grid element
+  vector<Vector<double> > parameter_dofs (coarse_dofs_per_cell_component,
+					  Vector<double>(fine_dofs_per_cell));
+				   // for each coarse dof: find its
+				   // position within the fine element
+				   // and set this value to one in the
+				   // respective vector (all other values
+				   // are zero by construction)
+  for (unsigned int local_coarse_dof=0;
+       local_coarse_dof<coarse_dofs_per_cell_component;
+       ++local_coarse_dof)
+    parameter_dofs[local_coarse_dof]
+      (fine_fe.component_to_system_index (fine_component,local_coarse_dof)) = 1.;
+
+
+				   // find out how many DoFs there are
+				   // on the grids belonging to the
+				   // components we want to match
+  unsigned int n_parameters_on_fine_grid=0;
+  if (true)
+    {
+				       // have a flag for each dof on
+				       // the fine grid and set it
+				       // to true if this is an
+				       // interesting dof. finally count
+				       // how many true's there
+      vector<bool>          dof_is_interesting (fine_grid.n_dofs(), false);
+      vector<unsigned int>  local_dof_indices (fine_fe.dofs_per_cell);
+      
+      for (DoFHandler<dim>::active_cell_iterator cell=fine_grid.begin_active();
+	   cell!=fine_grid.end(); ++cell)
+	{
+	  cell->get_dof_indices (local_dof_indices);
+	  for (unsigned int i=0; i<fine_fe.dofs_per_cell; ++i)
+	    if (fine_fe.system_to_component_index(i).first == fine_component)
+	      dof_is_interesting[local_dof_indices[i]] = true;
+	};
+
+      n_parameters_on_fine_grid = count (dof_is_interesting.begin(),
+					    dof_is_interesting.end(),
+					    true);
+    };
+
+				   // vector to hold the representation of
+				   // a single degree of freedom on the
+				   // coarse grid (for the selected fe)
+				   // on the fine grid
+  Vector<double> global_parameter_representation (n_fine_dofs);
+
+				   // store the weights with which a dof
+				   // on the parameter grid contributes
+				   // to a dof on the fine grid. see the
+				   // long doc below for more info
+				   //
+				   // allocate as many rows as there are
+				   // parameter dofs on the coarse grid
+				   // and as many columns as there are
+				   // parameter dofs on the fine grid.
+				   //
+				   // weight_mapping is used to map the
+				   // global (fine grid) parameter dof
+				   // indices to the columns
+  FullMatrix<double> weights (n_coarse_dofs,
+			      n_parameters_on_fine_grid);
+				   // this is this mapping. there is one
+				   // entry for each dof on the fine grid;
+				   // if it is a parameter dof, then its
+				   // value is the column in weights for
+				   // that parameter dof, if it is any
+				   // other dof, then its value is -1,
+				   // indicating an error
+  vector<int> weight_mapping (n_fine_dofs, -1);
+
+				   // set up this mapping
+  if (true)
+    {
+				       // list of local dof numbers
+				       // belonging to the parameter
+				       // part of the FE on the fine grid
+      vector<unsigned int> parameter_dofs;
+      for (unsigned int i=0; i<coarse_dofs_per_cell; ++i)
+	parameter_dofs.push_back(fine_fe.component_to_system_index(coarse_component,i));
+
+      vector<unsigned int> local_dof_indices(fine_fe.dofs_per_cell);
+      unsigned int next_free_index=0;
+      for (typename DoFHandler<dim>::active_cell_iterator cell=fine_grid.begin_active();
+	   cell != fine_grid.end(); ++cell)
+	{
+	  cell->get_dof_indices (local_dof_indices);
+	  for (unsigned int i=0; i<fine_fe.dofs_per_cell; ++i)
+					     // if this DoF is a
+					     // parameter dof and has
+					     // not yet been numbered,
+					     // then do so
+	    if ((fine_fe.system_to_component_index(i).first == fine_component) &&
+		(weight_mapping[local_dof_indices[i]] == -1))
+	      {
+		weight_mapping[local_dof_indices[i]] = next_free_index;
+		++next_free_index;
+	      };
+	};
+
+      Assert (next_free_index == n_parameters_on_fine_grid,
+	      ExcInternalError());
+    };
+  
+  
+				   // for each cell on the parameter grid:
+				   // find out which degrees of freedom on the
+				   // fine grid correspond in which way to
+				   // the degrees of freedom on the parameter
+				   // grid
+				   //
+				   // since for continuous FEs some
+				   // dofs exist on more than one
+				   // cell, we have to track which
+				   // ones were already visited. the
+				   // problem is that if we visit a
+				   // dof first on one cell and
+				   // compute its weight with respect
+				   // to some global dofs to be
+				   // non-zero, and later visit the
+				   // dof again on another cell and
+				   // (since we are on another cell)
+				   // recompute the weights with
+				   // respect to the same dofs as
+				   // above to be zero now, we have to
+				   // preserve them. we therefore
+				   // overwrite all weights if they
+				   // are nonzero and do not enforce
+				   // zero weights since that might be
+				   // only due to the fact that we are
+				   // on another cell.
+				   //
+				   // example:
+				   // coarse grid
+				   //  |     |     |
+				   //  *-----*-----*
+				   //  | cell|cell |
+				   //  |  1  |  2  |
+				   //  |     |     |
+				   //  0-----1-----*
+				   //
+				   // fine grid
+				   //  |  |  |  |  |
+				   //  *--*--*--*--*
+				   //  |  |  |  |  |
+				   //  *--*--*--*--*
+				   //  |  |  |  |  |
+				   //  *--x--y--*--*
+				   //
+				   // when on cell 1, we compute the
+				   // weights of dof 'x' to be 1/2
+				   // from parameter dofs 0 and 1,
+				   // respectively. however, when
+				   // later we are on cell 2, we again
+				   // compute the prolongation of
+				   // shape function 1 restricted to
+				   // cell 2 to the globla grid and
+				   // find that the weight of global
+				   // dof 'x' now is zero. however, we
+				   // should not overwrite the old
+				   // value.
+				   //
+				   // we therefore always only set
+				   // nonzero values. why adding up is
+				   // not useful: dof 'y' would get
+				   // weight 1 from parameter dof 1 on
+				   // both cells 1 and 2, but the
+				   // correct weight is nevertheless
+				   // only 1.
+  typename DoFHandler<dim>::active_cell_iterator cell;
+  vector<unsigned int> parameter_dof_indices (coarse_fe.dofs_per_cell);
+  
+  for (cell=coarse_grid.begin_active(); cell!=coarse_grid.end(); ++cell)
+    {
+				       // get the global indices of the
+				       // parameter dofs on this parameter
+				       // grid cell
+      cell->get_dof_indices (parameter_dof_indices);
+
+				       // loop over all dofs on this
+				       // cell and check whether they
+				       // are interesting for us
+      for (unsigned int local_parameter_dof=0;
+	   local_parameter_dof<coarse_fe.dofs_per_cell;
+	   ++local_parameter_dof)
+	if (coarse_fe.system_to_component_index(local_parameter_dof).first
+	    ==
+	    coarse_component)
+	  {
+	    global_parameter_representation.clear ();
+	    
+					     // distribute the representation of
+					     // #local_parameter_dof# on the
+					     // parameter grid cell #cell# to
+					     // the global data space
+	    coarse_to_fine_grid_map[cell]->
+	      set_dof_values_by_interpolation (parameter_dofs[local_parameter_dof],
+					       global_parameter_representation);
+					     // now that we've got the global
+					     // representation of each parameter
+					     // dof, we've only got to clobber the
+					     // non-zero entries in that vector and
+					     // store the result
+					     //
+					     // what we have learned: if entry #i#
+					     // of the global vector holds the value
+					     // #v[i]#, then this is the weight with
+					     // which the present dof contributes
+					     // to #i#. there may be several such
+					     // #i#s and their weights' sum should
+					     // be one. Then, #v[i]# should be
+					     // equal to #\sum_j w_{ij} p[j]# with
+					     // #p[j]# be the values of the degrees
+					     // of freedom on the coarse grid. we
+					     // can thus compute constraints which
+					     // link the degrees of freedom #v[i]#
+					     // on the fine grid to those on the
+					     // coarse grid, #p[j]#. Now to use
+					     // these as real constraints, rather
+					     // than as additional equations, we
+					     // have to identify representants
+					     // among the #i# for each #j#. this will
+					     // be done by simply taking the first
+					     // #i# for which #w_{ij}==1#.
+	    for (unsigned int i=0; i<global_parameter_representation.size(); ++i)
+					       // set this weight if it belongs
+					       // to a parameter dof.
+	      if (weight_mapping[i] != -1)
+		{
+						   // only overwrite old
+						   // value if not by
+						   // zero
+		  if (global_parameter_representation(i) != 0)
+		    {
+		      const unsigned int wi = parameter_dof_indices[local_parameter_dof],
+					 wj = weight_mapping[i];
+		      weights(wi,wj) = global_parameter_representation(i);
+		    };
+		}
+	      else
+		Assert (global_parameter_representation(i) == 0,
+			ExcInternalError());
+	  };
+    };
+  
+				   // ok, now we have all weights for each
+				   // dof on the fine grid. if in debug
+				   // mode lets see if everything went smooth,
+				   // i.e. each dof has sum of weights one
+				   //
+				   // in other words this means that
+				   // if the sum of all shape
+				   // functions on the parameter grid
+				   // is one (which is always the
+				   // case), then the representation
+				   // on the state grid should be as
+				   // well (division of unity)
+				   //
+				   // if the parameter grid has more
+				   // than one component, then the
+				   // respective dofs of the other
+				   // components have sum of weights
+				   // zero, of course. we do not
+				   // explicitely ask which component
+				   // a dof belongs to, but this at
+				   // least tests some errors
+#ifdef DEBUG
+  for (unsigned int col=0; col<weights.n(); ++col)
+    {
+      double sum=0;
+      for (unsigned int row=0; row<weights.m(); ++row)
+	sum += weights(row,col);
+      Assert ((sum==1) ||
+	      ((coarse_fe.n_components()>1) && (sum==0)), ExcInternalError());
+    };
+#endif
+
+
+				   // now we know that the weights in each
+				   // row constitute a constraint. enter this
+				   // into the constraints object
+				   //
+				   // first task: for each parameter dof on
+				   // the parameter grid, find a representant
+				   // on the fine, global grid. this is
+				   // possible since we use conforming
+				   // finite element. we take this representant
+				   // to be the first element in this row with
+				   // weight identical to one
+  vector<int> representants(weights.m(), -1);
+  for (unsigned int parameter_dof=0; parameter_dof<weights.m(); ++parameter_dof)
+    {
+      unsigned int column=0;
+      for (; column<weights.n(); ++column)
+	if (weights(parameter_dof,column) == 1)
+	  break;
+      Assert (column < weights.n(), ExcInternalError());
+
+				       // now we know in which column of
+				       // weights the representant is, but
+				       // we don't know its global index. get
+				       // it using the inverse operation of
+				       // the weight_mapping
+      unsigned int global_dof=0;
+      for (; global_dof<weight_mapping.size(); ++global_dof)
+	if (weight_mapping[global_dof] == static_cast<int>(column))
+	  break;
+      Assert (global_dof < weight_mapping.size(), ExcInternalError());
+
+				       // now enter the representants global
+				       // index into our list
+      representants[parameter_dof] = global_dof;
+    };
+  
+  for (unsigned int global_dof=0; global_dof<n_fine_dofs; ++global_dof)
+    if (weight_mapping[global_dof] != -1)
+				       // this global dof is a parameter
+				       // dof, so it may carry a constraint
+				       // note that for each global dof,
+				       // the sum of weights shall be one,
+				       // so we can find out whether this
+				       // dof is constrained in the following
+				       // way: if the only weight in this row
+				       // is a one, and the representant for
+				       // the parameter dof of the line in
+				       // which this one is is the present
+				       // dof, then we consider this dof
+				       // to be unconstrained. otherwise,
+				       // all other dofs are constrained
+      {
+	unsigned int first_used_row=0;
+	for (; first_used_row<weights.m(); ++first_used_row)
+	  if (weights(first_used_row,weight_mapping[global_dof]) != 0)
+	    break;
+
+	if ((weights(first_used_row,weight_mapping[global_dof]) == 1) &&
+	    (representants[first_used_row] == static_cast<int>(global_dof)))
+					   // dof unconstrained
+	  continue;
+
+					 // otherwise enter all constraints
+	constraints.add_line (global_dof);
+	for (unsigned int row=first_used_row; row<weights.m(); ++row)
+	  if (weights(row,weight_mapping[global_dof]) != 0)
+	    constraints.add_entry (global_dof,
+				   representants[row],
+				   weights(row,weight_mapping[global_dof]));
+      };
+};
+
+
+
+
+
+
 // explicit instantiations
 #if deal_II_dimension > 1
 template void
@@ -755,3 +1183,12 @@ template void DoFTools::extract_level_dofs(unsigned int level,
 					   const vector<bool>& local_select,
 					   vector<bool>& selected_dofs);
 
+
+template
+void
+DoFTools::compute_intergrid_constraints (const DoFHandler<deal_II_dimension> &,
+					 const unsigned int                   ,
+					 const DoFHandler<deal_II_dimension> &,
+					 const unsigned int                   ,
+					 const InterGridMap<DoFHandler,deal_II_dimension> &,
+					 ConstraintMatrix                    &);
-- 
2.39.5