]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Cleanup some of the code; remove one level of indirection and re-sort functions to...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 13 Aug 2006 00:51:14 +0000 (00:51 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 13 Aug 2006 00:51:14 +0000 (00:51 +0000)
git-svn-id: https://svn.dealii.org/trunk@13690 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 8147d1a7d9dfc5dbe9fcc394325bef1da0bd859b..cb60456dff3ea7b6feda1e9a5f4fcf96d8140caa 100644 (file)
@@ -1625,20 +1625,67 @@ namespace internal
 #if deal_II_dimension == 1
     static
     void
-    do_make_hanging_node_constraints (const ::DoFHandler<1> &,
-                                     ConstraintMatrix    &,
-                                     ::internal::int2type<1>)
+    make_hp_hanging_node_constraints (const ::DoFHandler<1> &,
+                                     ConstraintMatrix    &)
     {
                                       // nothing to do for regular
                                       // dof handlers in 1d
     }
 
 
+
+    static
+    void
+    make_oldstyle_hanging_node_constraints (const ::DoFHandler<1> &,
+                                           ConstraintMatrix    &,
+                                           ::internal::int2type<1>)
+    {
+                                      // nothing to do for regular
+                                      // dof handlers in 1d
+    }
+    
+
+    static
+    void
+    make_hp_hanging_node_constraints (const ::MGDoFHandler<1> &,
+                                     ConstraintMatrix    &)
+    {
+                                      // nothing to do for regular
+                                      // dof handlers in 1d
+    }
+
+
+
+    static
+    void
+    make_oldstyle_hanging_node_constraints (const ::MGDoFHandler<1> &,
+                                           ConstraintMatrix    &,
+                                           ::internal::int2type<1>)
+    {
+                                      // nothing to do for regular
+                                      // dof handlers in 1d
+    }
+    
+
+    static
+    void
+    make_hp_hanging_node_constraints (const ::hp::DoFHandler<1> &/*dof_handler*/,
+                                     ConstraintMatrix        &/*constraints*/)
+    {
+                                      // we may have to compute
+                                      // constraints for
+                                      // vertices. gotta think about
+                                      // that a bit more
+//TODO[WB]: think about what to do here...      
+    }
+
+
+
     static
     void
-    do_make_hanging_node_constraints (const ::hp::DoFHandler<1> &/*dof_handler*/,
-                                     ConstraintMatrix        &/*constraints*/,
-                                     ::internal::int2type<1>)
+    make_oldstyle_hanging_node_constraints (const ::hp::DoFHandler<1> &/*dof_handler*/,
+                                           ConstraintMatrix        &/*constraints*/,
+                                           ::internal::int2type<1>)
     {
                                       // we may have to compute
                                       // constraints for
@@ -1649,29 +1696,15 @@ namespace internal
 #endif
 
 
+
 #if deal_II_dimension == 2
     template <class DH>
     static
     void
-    do_make_hanging_node_constraints (const DH         &dof_handler,
-                                     ConstraintMatrix &constraints,
-                                     ::internal::int2type<2>)
+    make_oldstyle_hanging_node_constraints (const DH         &dof_handler,
+                                           ConstraintMatrix &constraints,
+                                           ::internal::int2type<2>)
     {
-                                      // Decide whether to use the
-                                      // new or old make_hanging_node_constraints
-                                      // function. If all the FiniteElement
-                                      // or all elements in a FECollection support
-                                      // the new face constraint matrix, the
-                                      // new code will be used.
-                                      // Otherwise, the old implementation is used
-                                      // for the moment.
-      if (dof_handler.get_fe().hp_constraints_are_implemented ())
-       {
-         do_make_hp_hanging_node_constraints (dof_handler,
-                                              constraints);
-         return;
-       }
-         
       const unsigned int dim = 2;
   
       std::vector<unsigned int> dofs_on_mother;
@@ -1807,14 +1840,237 @@ namespace internal
                      ExcInternalError());
            } 
     }
+#endif      
+
+
+#if deal_II_dimension == 3
+    template <class DH>
+    static
+    void
+    make_oldstyle_hanging_node_constraints (const DH         &dof_handler,
+                                           ConstraintMatrix &constraints,
+                                           ::internal::int2type<3>)
+    {
+      const unsigned int dim = 3;
+  
+      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);
+         
+                                              // 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());
+         
+             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);
+         
+                                              // 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);
+         
+                                              // 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++]
+                 = 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);
+         
+                                              // 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);
+         
+                                              // 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
 
-#if deal_II_dimension > 1
+    
     template <class DH>
     static
     void
-    do_make_hp_hanging_node_constraints (const DH         &dof_handler,
-                                        ConstraintMatrix &constraints)
+    make_hp_hanging_node_constraints (const DH         &dof_handler,
+                                     ConstraintMatrix &constraints)
     {
       const unsigned int dim = DH::dimension;
       
@@ -2142,247 +2398,6 @@ namespace internal
                 }                
             }
     }
-#endif
-
-
-
-#if deal_II_dimension == 3
-    template <class DH>
-    static
-    void
-    do_make_hanging_node_constraints (const DH         &dof_handler,
-                                     ConstraintMatrix &constraints,
-                                     ::internal::int2type<3>)
-    {
-                                      // Decide whether to use the
-                                      // new or old make_hanging_node_constraints
-                                      // function. If all the FiniteElement
-                                      // or all elements in a FECollection support
-                                      // the new face constraint matrix, the
-                                      // new code will be used.
-                                      // Otherwise, the old implementation is used
-                                      // for the moment.
-      if (dof_handler.get_fe().hp_constraints_are_implemented ())
-       {
-         do_make_hp_hanging_node_constraints (dof_handler,
-                                              constraints);
-         return;
-       }      
-         
-      const unsigned int dim = 3;
-  
-      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);
-         
-                                              // 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());
-         
-             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);
-         
-                                              // 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);
-         
-                                              // 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++]
-                 = 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);
-         
-                                              // 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);
-         
-                                              // 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
-    
   }
 }
 
@@ -2393,13 +2408,23 @@ void
 DoFTools::make_hanging_node_constraints (const DH &dof_handler,
                                         ConstraintMatrix &constraints)
 {
-                                  // distribute this one template
-                                  // function to the appropriate
-                                  // functions in 1d, 2d, and 3d
-  internal::DoFTools
-    ::do_make_hanging_node_constraints (dof_handler,
-                                        constraints,
-                                        ::internal::int2type<DH::dimension>());
+                                  // Decide whether to use the
+                                  // new or old make_hanging_node_constraints
+                                  // function. If all the FiniteElement
+                                  // or all elements in a FECollection support
+                                  // the new face constraint matrix, the
+                                  // new code will be used.
+                                  // Otherwise, the old implementation is used
+                                  // for the moment.
+  if (dof_handler.get_fe().hp_constraints_are_implemented ())
+    internal::DoFTools::
+      make_hp_hanging_node_constraints (dof_handler,
+                                       constraints);
+  else
+    internal::DoFTools::
+      make_oldstyle_hanging_node_constraints (dof_handler,
+                                             constraints,
+                                             ::internal::int2type<DH::dimension>());
 }
 
 

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.