]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement make_hanging_node_constraints for hp in the case of all the same elements...
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 4 May 2006 16:47:26 +0000 (16:47 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 4 May 2006 16:47:26 +0000 (16:47 +0000)
git-svn-id: https://svn.dealii.org/trunk@13053 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/dofs/dof_tools.h
deal.II/deal.II/source/dofs/dof_tools.cc

index 1ed08d628f57cfeb6eae88447ba8f0258afad3d8..c264a55a952771abacc7fb69338f0443dbbc85e4 100644 (file)
@@ -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);
                                     //@}
     
                                     /**
index 6e2683de6a236ade517950cc7df3f8eb6609770c..3f84e789e51613aec077c3fd6ba8ebaa6a31c56e 100644 (file)
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.