From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Thu, 4 May 2006 16:47:26 +0000 (+0000)
Subject: Implement make_hanging_node_constraints for hp in the case of all the same elements... 
X-Git-Tag: v8.0.0~11765
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=494435ca8a9ceb5ba30d206f0a6afaf2f38d1e19;p=dealii.git

Implement make_hanging_node_constraints for hp in the case of all the same elements on all cells.


git-svn-id: https://svn.dealii.org/trunk@13053 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/deal.II/deal.II/include/dofs/dof_tools.h b/deal.II/deal.II/include/dofs/dof_tools.h
index 1ed08d628f..c264a55a95 100644
--- a/deal.II/deal.II/include/dofs/dof_tools.h
+++ b/deal.II/deal.II/include/dofs/dof_tools.h
@@ -590,77 +590,63 @@ class DoFTools
 				      */
     
 				     /**
-				      * Make up the constraints which
-				      * is result from the use of hanging
-				      * nodes. The object into which these
-				      * are inserted is later
-				      * used to condensate the global
-				      * system matrices and to prolong
-				      * the solution vectors from the true
+				      * Compute the constraints
+				      * resulting from the presence of
+				      * hanging nodes. The object into
+				      * which these are inserted is
+				      * later used to condense the
+				      * global system matrix and right
+				      * hand side, and to extend the
+				      * solution vectors from the true
 				      * degrees of freedom also to the
-				      * constraint nodes.
+				      * constraint nodes. This
+				      * function is explained in
+				      * detail in the @ref step_6
+				      * "step-6" tutorial program and
+				      * is used in almost all
+				      * following programs as well.
 				      *
-				      * Since this method does not make sense in
-				      * one dimension, the function returns
-				      * immediately. The object is not cleared
-				      * before use, so you should make sure
-				      * it containts only constraints you still
-				      * want; otherwise call the @p clear
-				      * function.
+				      * This function does not clear
+				      * the constraint matrix object
+				      * before use, in order to allow
+				      * adding constraints from
+				      * different sources to the same
+				      * object. You therefore need to
+				      * make sure it contains only
+				      * constraints you still want;
+				      * otherwise call the
+				      * ConstraintMatrix::clear()
+				      * function.  Likewise, this
+				      * function does not close the
+				      * object since you may want to
+				      * enter other constraints later
+				      * on yourself.
 				      *
-				      * To condense a given sparsity pattern,
-				      * use ConstraintMatrix@p ::condense.
-				      * Before doing so, you need to close
-				      * the constraint object, which must be
-				      * done after all constraints are entered.
-				      * This function does not close the object
-				      * since you may want to enter other
-				      * constraints later on yourself.
-				      */
-    static void
-    make_hanging_node_constraints (const DoFHandler<1> &dof_handler,
-				   ConstraintMatrix      &constraints);
-
-				     /**
-				      * Declaration of same function
-				      * for different space dimension.
-				      */
-    static void
-    make_hanging_node_constraints (const DoFHandler<2> &dof_handler,
-				   ConstraintMatrix    &constraints);
-
-				     /**
-				      * Declaration of same function
-				      * for different space dimension.
-				      */
-    static void
-    make_hanging_node_constraints (const DoFHandler<3> &dof_handler,
-				   ConstraintMatrix    &constraints);
-
-				     /**
-				      * Declaration of same function
-				      * for hp::DoFHandler
-				      */
-    static void
-    make_hanging_node_constraints (const hp::DoFHandler<1> &dof_handler,
-				   ConstraintMatrix      &constraints);
-
-				     /**
-				      * Declaration of same function
-				      * for different space dimension.
-				      */
-    static void
-    make_hanging_node_constraints (const hp::DoFHandler<2> &dof_handler,
-				   ConstraintMatrix      &constraints);
-
-				     /**
-				      * Declaration of same function
-				      * for different space dimension.
+				      * In the hp-case, i.e. when the
+				      * argument is of type
+				      * hp::DoFHandler, we consider
+				      * constraints due to different
+				      * finite elements used on two
+				      * sides of a face between cells
+				      * as hanging nodes as well. In
+				      * other words, for hp finite
+				      * elements, this function
+				      * computes all constraints due
+				      * to differing mesh sizes (h) or
+				      * polynomial degrees (p) between
+				      * adjacent cells.
+				      *
+				      * The template argument (and by
+				      * consequence the type of the
+				      * first argument to this
+				      * function) can be either a
+				      * ::DoFHandler, hp::DoFHandler,
+				      * or MGDoFHandler.
 				      */
+    template <class DH>
     static void
-    make_hanging_node_constraints (const hp::DoFHandler<3> &dof_handler,
-				   ConstraintMatrix      &constraints);
-
+    make_hanging_node_constraints (const DH         &dof_handler,
+				   ConstraintMatrix &constraints);
 				     //@}
     
 				     /**
diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc
index 6e2683de6a..3f84e789e5 100644
--- a/deal.II/deal.II/source/dofs/dof_tools.cc
+++ b/deal.II/deal.II/source/dofs/dof_tools.cc
@@ -1648,275 +1648,424 @@ DoFTools::make_flux_sparsity_pattern (
 
 
 
-#if deal_II_dimension == 1
-
-void DoFTools::make_hanging_node_constraints (
-  const DoFHandler<1> &,
-  ConstraintMatrix &)
+namespace internal
 {
-				   // nothing to be done here
-}
+  namespace DoFTools
+  {
+    template <int> class int2type 
+    {};
 
 
-void DoFTools::make_hanging_node_constraints (
-  const hp::DoFHandler<1> &,
-  ConstraintMatrix &)
-{
-                                   // we may have to compute
-                                   // constraints for vertices
-  Assert (false, ExcNotImplemented());
-}
+#if deal_II_dimension == 1
+    void
+    do_make_hanging_node_constraints (const ::DoFHandler<1> &,
+				      ConstraintMatrix    &,
+				      int2type<1>)
+    {
+				       // nothing to do for regular
+				       // dof handlers in 1d
+    }
 
-#endif
 
+    void
+    do_make_hanging_node_constraints (const ::hp::DoFHandler<1> &/*dof_handler*/,
+				      ConstraintMatrix        &/*constraints*/,
+				      int2type<1>)
+    {
+				       // we may have to compute
+				       // constraints for
+				       // vertices. gotta think about
+				       // that a bit more
+//TODO[WB]: think about what to do here...      
+    }
+#endif
 
 
 #if deal_II_dimension == 2
-
-void DoFTools::make_hanging_node_constraints (
-  const DoFHandler<2> &dof_handler,
-  ConstraintMatrix    &constraints)
-{
-  const unsigned int dim = 2;
-  
-  const FiniteElement<dim> &fe   = dof_handler.get_fe();
+    template <class DH>
+    void
+    do_make_hanging_node_constraints (const DH         &dof_handler,
+				      ConstraintMatrix &constraints,
+				      int2type<2>)
+    {
+      const unsigned int dim = 2;
   
-				   // have space for the degrees of
-				   // freedom on mother and child
-				   // lines
-  const unsigned int n_dofs_on_mother   = 2*fe.dofs_per_vertex + fe.dofs_per_line,
-		     n_dofs_on_children = fe.dofs_per_vertex + 2*fe.dofs_per_line;
-
-  std::vector<unsigned int> dofs_on_mother(n_dofs_on_mother);
-  std::vector<unsigned int> dofs_on_children(n_dofs_on_children);
-
-  Assert(n_dofs_on_mother == fe.constraints().n(),
-	 ExcDimensionMismatch(n_dofs_on_mother,
-			      fe.constraints().n()));
-  Assert(n_dofs_on_children == fe.constraints().m(),
-	 ExcDimensionMismatch(n_dofs_on_children,
-			      fe.constraints().m()));
-
-				   // loop over all lines; only on
-				   // lines there can be constraints.
-				   // We do so by looping over all
-				   // active cells and checking
-				   // whether any of the faces are
-				   // refined which can only be from
-				   // the neighboring cell because
-				   // this one is active. In that
-				   // case, the face is subject to
-				   // constraints
-				   //
-				   // note that even though we may
-				   // visit a face twice if the
-				   // neighboring cells are equally
-				   // refined, we can only visit each
-				   // face with hanging nodes once
-  DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+      std::vector<unsigned int> dofs_on_mother;
+      std::vector<unsigned int> dofs_on_children;
+
+				       // loop over all lines; only on
+				       // lines there can be constraints.
+				       // We do so by looping over all
+				       // active cells and checking
+				       // whether any of the faces are
+				       // refined which can only be from
+				       // the neighboring cell because
+				       // this one is active. In that
+				       // case, the face is subject to
+				       // constraints
+				       //
+				       // note that even though we may
+				       // visit a face twice if the
+				       // neighboring cells are equally
+				       // refined, we can only visit each
+				       // face with hanging nodes once
+      typename DH::active_cell_iterator cell = dof_handler.begin_active(),
 					endc = dof_handler.end();
-  for (; cell!=endc; ++cell)
-    for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-      if (cell->face(face)->has_children()) 
-	{
-	  const DoFHandler<dim>::line_iterator line = cell->face(face);
-	  
-					   // fill the dofs indices. Use same
-					   // enumeration scheme as in
-					   // @p{FiniteElement::constraints()}
-	  unsigned int next_index = 0;
-	  for (unsigned int vertex=0; vertex<2; ++vertex)
-	    for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
-	      dofs_on_mother[next_index++] = line->vertex_dof_index(vertex,dof);
-	  for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
-	    dofs_on_mother[next_index++] = line->dof_index(dof);
-	  Assert (next_index == dofs_on_mother.size(),
-		  ExcInternalError());
-	  
-	  next_index = 0;
-	  for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
-	    dofs_on_children[next_index++] = line->child(0)->vertex_dof_index(1,dof);
-	  for (unsigned int child=0; child<2; ++child)
-	    for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
-	      dofs_on_children[next_index++] = line->child(child)->dof_index(dof);
-	  Assert (next_index == dofs_on_children.size(),
-		  ExcInternalError());
-	  
-					   // for each row in the constraint
-					   // matrix for this line:
-	  for (unsigned int row=0; row!=dofs_on_children.size(); ++row) 
+      for (; cell!=endc; ++cell)
+	for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+	  if (cell->face(face)->has_children()) 
 	    {
-	      constraints.add_line (dofs_on_children[row]);
-	      for (unsigned int i=0; i!=dofs_on_mother.size(); ++i)
-		constraints.add_entry (dofs_on_children[row],
-				       dofs_on_mother[i],
-				       fe.constraints()(row,i));
-	    }
-	}
-}
-
+					       // so now we've found a
+					       // face of an active
+					       // cell that has
+					       // children. that means
+					       // that there are
+					       // hanging nodes here.
+
+					       // in any case, faces
+					       // can have at most two
+					       // active fe indices,
+					       // but here the face
+					       // can have only one
+					       // (namely the same as
+					       // that from the cell
+					       // we're sitting on),
+					       // and each of the
+					       // children can have
+					       // only one as
+					       // well. check this
+	      Assert (cell->face(face)->n_active_fe_indices() == 1,
+		      ExcInternalError());
+	      Assert (cell->face(face)->fe_index_is_active(cell->active_fe_index())
+		      == true,
+		      ExcInternalError());
+	      for (unsigned int c=0; c<GeometryInfo<dim>::subfaces_per_face; ++c)
+		Assert (cell->face(face)->child(c)->n_active_fe_indices() == 1,
+			ExcInternalError());
 
-void DoFTools::make_hanging_node_constraints (
-    const hp::DoFHandler<2> &/*dof_handler*/,
-    ConstraintMatrix    &/*constraints*/)
-{
-//TODO[?]: Implement (required for continuous elements)
-    Assert (false, ExcNotImplemented());
-}
+					       // right now, all that
+					       // is implemented is
+					       // the case that both
+					       // sides use the same
+					       // fe
+	      for (unsigned int c=0; c<GeometryInfo<dim>::subfaces_per_face; ++c)
+		Assert (cell->face(face)->child(c)
+			->fe_index_is_active(cell->active_fe_index()) == true,
+			ExcNotImplemented());
+
+					       // ok, start up the work
+	      const FiniteElement<dim> &fe       = cell->get_fe();
+	      const unsigned int        fe_index = cell->active_fe_index();
+		
+	      const unsigned int
+		n_dofs_on_mother   = 2*fe.dofs_per_vertex + fe.dofs_per_line,
+		n_dofs_on_children = fe.dofs_per_vertex + 2*fe.dofs_per_line;
 
+	      dofs_on_mother.resize (n_dofs_on_mother);
+	      dofs_on_children.resize (n_dofs_on_children);
+	      
+	      Assert(n_dofs_on_mother == fe.constraints().n(),
+		     ExcDimensionMismatch(n_dofs_on_mother,
+					  fe.constraints().n()));
+	      Assert(n_dofs_on_children == fe.constraints().m(),
+		     ExcDimensionMismatch(n_dofs_on_children,
+					  fe.constraints().m()));	      
+
+	      const typename DH::line_iterator this_face = cell->face(face);
+	      
+					       // fill the dofs indices. Use same
+					       // enumeration scheme as in
+					       // @p{FiniteElement::constraints()}
+	      unsigned int next_index = 0;
+	      for (unsigned int vertex=0; vertex<2; ++vertex)
+		for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
+		  dofs_on_mother[next_index++] = this_face->vertex_dof_index(vertex,dof,
+									     fe_index);
+	      for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
+		dofs_on_mother[next_index++] = this_face->dof_index(dof, fe_index);
+	      Assert (next_index == dofs_on_mother.size(),
+		      ExcInternalError());
+	  
+	      next_index = 0;
+	      for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
+		dofs_on_children[next_index++]
+		  = this_face->child(0)->vertex_dof_index(1,dof,fe_index);
+	      for (unsigned int child=0; child<2; ++child)
+		for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
+		  dofs_on_children[next_index++]
+		    = this_face->child(child)->dof_index(dof, fe_index);
+	      Assert (next_index == dofs_on_children.size(),
+		      ExcInternalError());
+	  
+					       // for each row in the constraint
+					       // matrix for this line:
+	      for (unsigned int row=0; row!=dofs_on_children.size(); ++row) 
+		{
+		  constraints.add_line (dofs_on_children[row]);
+		  for (unsigned int i=0; i!=dofs_on_mother.size(); ++i)
+		    constraints.add_entry (dofs_on_children[row],
+					   dofs_on_mother[i],
+					   fe.constraints()(row,i));
+		}
+	    }
+	  else
+	    {
+					       // this face has no
+					       // children, but it
+					       // could still be that
+					       // it is shared by two
+					       // cells that use a
+					       // different fe index
+	      Assert (cell->face(face)->n_active_fe_indices() == 1,
+		      ExcNotImplemented());
+	      Assert (cell->face(face)
+		      ->fe_index_is_active(cell->active_fe_index()) == true,
+		      ExcInternalError());
+	    } 
+    }
 #endif
 
 
 
 #if deal_II_dimension == 3
-
-void DoFTools::make_hanging_node_constraints (
-  const DoFHandler<3> &dof_handler,
-  ConstraintMatrix    &constraints)
-{
-  const unsigned int dim = 3;
-  
-  const FiniteElement<dim> &fe   = dof_handler.get_fe();
+    template <class DH>
+    void
+    do_make_hanging_node_constraints (const DH         &dof_handler,
+				      ConstraintMatrix &constraints,
+				      int2type<3>)
+    {
+      const unsigned int dim = 3;
   
-				   // have space for the degrees of
-				   // freedom on mother and child
-				   // lines
-  const unsigned int
-    n_dofs_on_mother   = (4*fe.dofs_per_vertex+
-			  4*fe.dofs_per_line+
-			  fe.dofs_per_quad),
-    n_dofs_on_children = (5*fe.dofs_per_vertex+
-			  12*fe.dofs_per_line+
-			  4*fe.dofs_per_quad);
-
-  std::vector<unsigned int> dofs_on_mother(n_dofs_on_mother);
-  std::vector<unsigned int> dofs_on_children(n_dofs_on_children);
-
-  Assert(n_dofs_on_mother == fe.constraints().n(),
-	 ExcDimensionMismatch(n_dofs_on_mother,
-			      fe.constraints().n()));
-  Assert(n_dofs_on_children == fe.constraints().m(),
-	 ExcDimensionMismatch(n_dofs_on_children,
-			      fe.constraints().m()));
-
-				   // loop over all lines; only on
-				   // lines there can be constraints.
-				   // We do so by looping over all
-				   // active cells and checking
-				   // whether any of the faces are
-				   // refined which can only be from
-				   // the neighboring cell because
-				   // this one is active. In that
-				   // case, the face is subject to
-				   // constraints
-				   //
-				   // note that even though we may
-				   // visit a face twice if the
-				   // neighboring cells are equally
-				   // refined, we can only visit each
-				   // face with hanging nodes once
-  DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
-					endc = dof_handler.end();
-  for (; cell!=endc; ++cell)
-    for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
-      if (cell->face(f)->has_children()) 
-	{
-	  const DoFHandler<dim>::face_iterator face = cell->face(f);
-	  
-					   // fill the dofs indices. Use same
-					   // enumeration scheme as in
-					   // @p{FiniteElement::constraints()}
-	  unsigned int next_index = 0;
-	  for (unsigned int vertex=0; vertex<4; ++vertex)
-	    for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
-	      dofs_on_mother[next_index++] = face->vertex_dof_index(vertex,dof);
-	  for (unsigned int line=0; line<4; ++line)
-	    for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
-	      dofs_on_mother[next_index++] = face->line(line)->dof_index(dof);
-	  for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof)
-	    dofs_on_mother[next_index++] = face->dof_index(dof);
-	  Assert (next_index == dofs_on_mother.size(),
-		  ExcInternalError());
+      std::vector<unsigned int> dofs_on_mother;
+      std::vector<unsigned int> dofs_on_children;
+
+				       // loop over all quads; only on
+				       // quads there can be constraints.
+				       // We do so by looping over all
+				       // active cells and checking
+				       // whether any of the faces are
+				       // refined which can only be from
+				       // the neighboring cell because
+				       // this one is active. In that
+				       // case, the face is subject to
+				       // constraints
+				       //
+				       // note that even though we may
+				       // visit a face twice if the
+				       // neighboring cells are equally
+				       // refined, we can only visit each
+				       // face with hanging nodes once
+      typename DH::active_cell_iterator cell = dof_handler.begin_active(),
+					    endc = dof_handler.end();
+      for (; cell!=endc; ++cell)
+	for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+	  if (cell->face(face)->has_children()) 
+	    {
+					       // so now we've found a
+					       // face of an active
+					       // cell that has
+					       // children. that means
+					       // that there are
+					       // hanging nodes here.
+
+					       // in any case, faces
+					       // can have at most two
+					       // active fe indices,
+					       // but here the face
+					       // can have only one
+					       // (namely the same as
+					       // that from the cell
+					       // we're sitting on),
+					       // and each of the
+					       // children can have
+					       // only one as
+					       // well. check this
+	      Assert (cell->face(face)->n_active_fe_indices() == 1,
+		      ExcInternalError());
+	      Assert (cell->face(face)->fe_index_is_active(cell->active_fe_index())
+		      == true,
+		      ExcInternalError());
+	      for (unsigned int c=0; c<GeometryInfo<dim>::subfaces_per_face; ++c)
+		Assert (cell->face(face)->child(c)->n_active_fe_indices() == 1,
+			ExcInternalError());
+
+					       // right now, all that
+					       // is implemented is
+					       // the case that both
+					       // sides use the same
+					       // fe, and not only
+					       // that but also that
+					       // all lines bounding
+					       // this face and the
+					       // children have the
+					       // same fe
+	      for (unsigned int c=0; c<GeometryInfo<dim>::subfaces_per_face; ++c)
+		{
+		  Assert (cell->face(face)->child(c)
+			  ->fe_index_is_active(cell->active_fe_index()) == true,
+			  ExcNotImplemented());
+		  for (unsigned int e=0; e<4; ++e)
+		    {
+		      Assert (cell->face(face)->child(c)->line(e)
+			      ->n_active_fe_indices() == 1,
+			      ExcNotImplemented());
+		      Assert (cell->face(face)->child(c)->line(e)
+			      ->fe_index_is_active(cell->active_fe_index()) == true,
+			      ExcNotImplemented());
+		    }
+		}
+	      for (unsigned int e=0; e<4; ++e)
+		{
+		  Assert (cell->face(face)->line(e)
+			  ->n_active_fe_indices() == 1,
+			  ExcNotImplemented());
+		  Assert (cell->face(face)->line(e)
+			  ->fe_index_is_active(cell->active_fe_index()) == true,
+			  ExcNotImplemented());
+		}
+	      
+					       // ok, start up the work
+	      const FiniteElement<dim> &fe       = cell->get_fe();
+	      const unsigned int        fe_index = cell->active_fe_index();
+
+	      const unsigned int
+		n_dofs_on_mother   = (4*fe.dofs_per_vertex+
+				      4*fe.dofs_per_line+
+				      fe.dofs_per_quad),
+		n_dofs_on_children = (5*fe.dofs_per_vertex+
+				      12*fe.dofs_per_line+
+				      4*fe.dofs_per_quad);	      
+	      
+	      dofs_on_mother.resize (n_dofs_on_mother);
+	      dofs_on_children.resize (n_dofs_on_children);
+	      
+	      Assert(n_dofs_on_mother == fe.constraints().n(),
+		     ExcDimensionMismatch(n_dofs_on_mother,
+					  fe.constraints().n()));
+	      Assert(n_dofs_on_children == fe.constraints().m(),
+		     ExcDimensionMismatch(n_dofs_on_children,
+					  fe.constraints().m()));	      
+
+	      const typename DH::face_iterator this_face = cell->face(face);
 	  
-	  next_index = 0;
-
-					   // assert some consistency
-					   // assumptions
-	  Assert ((face->child(0)->vertex_index(3) ==
-		   face->child(1)->vertex_index(2)) &&
-		  (face->child(0)->vertex_index(3) ==
-		   face->child(2)->vertex_index(1)) &&
-		  (face->child(0)->vertex_index(3) ==
-		   face->child(3)->vertex_index(0)),
-		  ExcInternalError());
-	  for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
-	    dofs_on_children[next_index++]
-	      = face->child(0)->vertex_dof_index(3,dof);
+					       // fill the dofs indices. Use same
+					       // enumeration scheme as in
+					       // @p{FiniteElement::constraints()}
+	      unsigned int next_index = 0;
+	      for (unsigned int vertex=0; vertex<4; ++vertex)
+		for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
+		  dofs_on_mother[next_index++] = this_face->vertex_dof_index(vertex,dof,
+									     fe_index);
+	      for (unsigned int line=0; line<4; ++line)
+		for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
+		  dofs_on_mother[next_index++]
+		    = this_face->line(line)->dof_index(dof, fe_index);
+	      for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof)
+		dofs_on_mother[next_index++] = this_face->dof_index(dof, fe_index);
+	      Assert (next_index == dofs_on_mother.size(),
+		      ExcInternalError());
 	  
-					   // dof numbers on the centers of
-					   // the lines bounding this face
-	  for (unsigned int line=0; line<4; ++line)
-	    for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
-	      dofs_on_children[next_index++]
-		= face->line(line)->child(0)->vertex_dof_index(1,dof);
+	      next_index = 0;
+
+					       // assert some consistency
+					       // assumptions
+	      Assert ((this_face->child(0)->vertex_index(3) ==
+		       this_face->child(1)->vertex_index(2)) &&
+		      (this_face->child(0)->vertex_index(3) ==
+		       this_face->child(2)->vertex_index(1)) &&
+		      (this_face->child(0)->vertex_index(3) ==
+		       this_face->child(3)->vertex_index(0)),
+		      ExcInternalError());
+	      for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
+		dofs_on_children[next_index++]
+		  = this_face->child(0)->vertex_dof_index(3,dof);
 	  
-					   // next the dofs on the lines interior
-					   // to the face; the order of these
-					   // lines is laid down in the
-					   // FiniteElement class documentation
-	  for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
-	    dofs_on_children[next_index++]
-	      = face->child(0)->line(1)->dof_index(dof);
-	  for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
-	    dofs_on_children[next_index++]
-	      = face->child(2)->line(1)->dof_index(dof);
-	  for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
-	    dofs_on_children[next_index++]
-	      = face->child(0)->line(3)->dof_index(dof);
-	  for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
-	    dofs_on_children[next_index++]
-	      = face->child(1)->line(3)->dof_index(dof);
+					       // dof numbers on the centers of
+					       // the lines bounding this face
+	      for (unsigned int line=0; line<4; ++line)
+		for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
+		  dofs_on_children[next_index++]
+		    = this_face->line(line)->child(0)->vertex_dof_index(1,dof, fe_index);
 	  
-					   // dofs on the bordering lines
-	  for (unsigned int line=0; line<4; ++line)
-	    for (unsigned int child=0; child<2; ++child)
-	      for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
+					       // next the dofs on the lines interior
+					       // to the face; the order of these
+					       // lines is laid down in the
+					       // FiniteElement class documentation
+	      for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
+		dofs_on_children[next_index++]
+		  = this_face->child(0)->line(1)->dof_index(dof, fe_index);
+	      for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
+		dofs_on_children[next_index++]
+		  = this_face->child(2)->line(1)->dof_index(dof, fe_index);
+	      for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
 		dofs_on_children[next_index++]
-		  = face->line(line)->child(child)->dof_index(dof);
+		  = this_face->child(0)->line(3)->dof_index(dof, fe_index);
+	      for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
+		dofs_on_children[next_index++]
+		  = this_face->child(1)->line(3)->dof_index(dof, fe_index);
 	  
-					   // finally, for the dofs interior
-					   // to the four child faces
-	  for (unsigned int child=0; child<4; ++child)
-	    for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof)
-	      dofs_on_children[next_index++]
-		= face->child(child)->dof_index(dof);
-	  Assert (next_index == dofs_on_children.size(),
-		  ExcInternalError());
+					       // dofs on the bordering lines
+	      for (unsigned int line=0; line<4; ++line)
+		for (unsigned int child=0; child<2; ++child)
+		  for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
+		    dofs_on_children[next_index++]
+		      = this_face->line(line)->child(child)->dof_index(dof, fe_index);
 	  
-					   // for each row in the constraint
-					   // matrix for this line:
-	  for (unsigned int row=0; row!=dofs_on_children.size(); ++row) 
-	    {
-	      constraints.add_line (dofs_on_children[row]);
-	      for (unsigned int i=0; i!=dofs_on_mother.size(); ++i)
-		constraints.add_entry (dofs_on_children[row],
-				       dofs_on_mother[i],
-				       fe.constraints()(row,i));
+					       // finally, for the dofs interior
+					       // to the four child faces
+	      for (unsigned int child=0; child<4; ++child)
+		for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof)
+		  dofs_on_children[next_index++]
+		    = this_face->child(child)->dof_index(dof, fe_index);
+	      Assert (next_index == dofs_on_children.size(),
+		      ExcInternalError());
+	  
+					       // for each row in the constraint
+					       // matrix for this line:
+	      for (unsigned int row=0; row!=dofs_on_children.size(); ++row) 
+		{
+		  constraints.add_line (dofs_on_children[row]);
+		  for (unsigned int i=0; i!=dofs_on_mother.size(); ++i)
+		    constraints.add_entry (dofs_on_children[row],
+					   dofs_on_mother[i],
+					   fe.constraints()(row,i));
+		}
 	    }
-	}
+	  else
+	    {
+					       // this face has no
+					       // children, but it
+					       // could still be that
+					       // it is shared by two
+					       // cells that use a
+					       // different fe index
+	      Assert (cell->face(face)->n_active_fe_indices() == 1,
+		      ExcNotImplemented());
+	      Assert (cell->face(face)
+		      ->fe_index_is_active(cell->active_fe_index()) == true,
+		      ExcInternalError());
+	    } 
+    }
+#endif
+    
+  }
 }
 
 
-void DoFTools::make_hanging_node_constraints (
-    const hp::DoFHandler<3> &/*dof_handler*/,
-    ConstraintMatrix    &/*constraints*/)
+
+template <class DH>
+void
+DoFTools::make_hanging_node_constraints (const DH &dof_handler,
+					 ConstraintMatrix &constraints)
 {
-//TODO:[?] Implement (required for continuous elements)
-    Assert (false, ExcNotImplemented());
+				   // distribute this one template
+				   // function to the appropriate
+				   // functions in 1d, 2d, and 3d
+  do_make_hanging_node_constraints (dof_handler,
+				    constraints,
+				    internal::DoFTools::int2type<DH::dimension>());
 }
 
-#endif
-
 
 
 template <class DH, typename Number>
@@ -4237,6 +4386,20 @@ DoFTools::make_flux_sparsity_pattern<DoFHandler<deal_II_dimension>,CompressedBlo
 #endif
 
 
+template
+void
+DoFTools::make_hanging_node_constraints (const DoFHandler<deal_II_dimension> &dof_handler,
+					 ConstraintMatrix &constraints);
+template
+void
+DoFTools::make_hanging_node_constraints (const MGDoFHandler<deal_II_dimension> &dof_handler,
+					 ConstraintMatrix &constraints);
+template
+void
+DoFTools::make_hanging_node_constraints (const hp::DoFHandler<deal_II_dimension> &dof_handler,
+					 ConstraintMatrix &constraints);
+
+
 
 template
 void