]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Don't use the user flags any more for make_hanging_node_constraints.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 15 May 2000 12:04:17 +0000 (12:04 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 15 May 2000 12:04:17 +0000 (12:04 +0000)
git-svn-id: https://svn.dealii.org/trunk@2855 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 80b2b24b2b69f9c5ccc861452f4e6fd8acd438f8..785142f226de31a82a4f1647458957e1f86eb7a1 100644 (file)
@@ -318,6 +318,7 @@ DoFTools::make_flux_sparsity_pattern (const DoFHandler<dim> &dof,
 }
 
 
+
 #if deal_II_dimension == 1
 
 template <>
@@ -330,6 +331,7 @@ void DoFTools::make_hanging_node_constraints (const DoFHandler<1> &,
 #endif
 
 
+
 #if deal_II_dimension == 2
 
 template <>
@@ -337,91 +339,87 @@ void DoFTools::make_hanging_node_constraints (const DoFHandler<2> &dof_handler,
                                              ConstraintMatrix    &constraints)
 {
   const unsigned int dim = 2;
-  const Triangulation<dim> &tria = dof_handler.get_tria();
-  const FiniteElement<dim> &fe   = dof_handler.get_fe();
   
-                                  // first mark all faces which are subject
-                                  // to 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
-  DoFHandler<dim>::line_iterator line = dof_handler.begin_line(),
-                                endl = dof_handler.end_line();
-  for (; line!=endl; ++line)
-    line->clear_user_flag ();
+  const FiniteElement<dim> &fe   = dof_handler.get_fe();
   
-  Triangulation<dim>::active_cell_iterator cell = tria.begin_active(),
-                                          endc = tria.end();
-  for (; cell!=endc; ++cell)
-    for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-      if (cell->face(face)->has_children()) 
-       cell->face(face)->set_user_flag();
-
                                   // have space for the degrees of
                                   // freedom on mother and child
                                   // lines
-  vector<unsigned int> dofs_on_mother(2*fe.dofs_per_vertex+
-                                     fe.dofs_per_line);
-  vector<unsigned int> dofs_on_children(fe.dofs_per_vertex+
-                                       2*fe.dofs_per_line);
-
-  Assert(2*fe.dofs_per_vertex+fe.dofs_per_line ==
-        fe.constraints().n(),
-        ExcDimensionMismatch(2*fe.dofs_per_vertex+
-                             fe.dofs_per_line,
+  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;
+
+  vector<unsigned int> dofs_on_mother(n_dofs_on_mother);
+  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(fe.dofs_per_vertex+2*fe.dofs_per_line ==
-        fe.constraints().m(),
-        ExcDimensionMismatch(3*fe.dofs_per_vertex+
-                             2*fe.dofs_per_line,
+  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.
-  for (line = dof_handler.begin_line(); line != endl; ++line)
-                                    // if dofs on this line are subject
-                                    // to constraints
-    if (line->user_flag_set() == true)
-      {        
-                                        // fill the dofs indices. Use same
-                                        // enumeration scheme as in
-                                        // #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)
+                                  // 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 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
+                                          // #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_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) 
-         {
-           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));
-         };
-      };
+           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) 
+           {
+             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));
+           };
+       };
 };
 
 #endif
 
 
+
 #if deal_II_dimension == 3
 
 template <>
@@ -429,152 +427,138 @@ void DoFTools::make_hanging_node_constraints (const DoFHandler<3> &dof_handler,
                                              ConstraintMatrix    &constraints)
 {
   const unsigned int dim = 3;
-  const Triangulation<dim> &tria = dof_handler.get_tria();
+  
   const FiniteElement<dim> &fe   = dof_handler.get_fe();
   
-                                  // first mark all faces which are subject
-                                  // to 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
-  DoFHandler<dim>::face_iterator face = dof_handler.begin_face(),
-                                endf = dof_handler.end_face();
-  for (; face!=endf; ++face)
-    face->clear_user_flag ();
+                                  // 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);
+
+  vector<unsigned int> dofs_on_mother(n_dofs_on_mother);
+  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()));
 
-  Triangulation<dim>::active_cell_iterator cell = tria.begin_active(),
-                                          endc = tria.end();
+                                  // 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 face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
       if (cell->face(face)->has_children()) 
-       cell->face(face)->set_user_flag();
-
-
-                                  // allocate space for dof indices
-  vector<unsigned int> dofs_on_mother(4*fe.dofs_per_vertex+
-                                     4*fe.dofs_per_line+
-                                     fe.dofs_per_quad);
-  vector<unsigned int> dofs_on_children(5*fe.dofs_per_vertex+
-                                       12*fe.dofs_per_line+
-                                       4*fe.dofs_per_quad);
-  Assert(4*fe.dofs_per_vertex+
-        4*fe.dofs_per_line+
-        fe.dofs_per_quad
-        ==
-        fe.constraints().n(),
-        ExcDimensionMismatch(4*fe.dofs_per_vertex+
-                             4*fe.dofs_per_line+
-                             fe.dofs_per_quad,
-                             fe.constraints().n()));
-  Assert(5*fe.dofs_per_vertex+
-        12*fe.dofs_per_line+
-        4*fe.dofs_per_quad
-        ==
-        fe.constraints().m(),
-        ExcDimensionMismatch(5*fe.dofs_per_vertex+
-                             12*fe.dofs_per_line+
-                             4*fe.dofs_per_quad,
-                             fe.constraints().m()));
-       
-                                  // loop over all faces; only on faces
-                                  // there can be constraints.
-  for (face=dof_handler.begin_face(); face != endf; ++face)
-                                    // if dofs on this line are subject
-                                    // to constraints
-    if (face->user_flag_set() == true)
-      {
-       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());
-       
-                                        // dof numbers on vertex at the center
-                                        // of the face, which is vertex 2 of
-                                        // child zero, or vertex 3 of child 1
-                                        // or vertex 0 of child 2 or vertex 1
-                                        // of child 3. We're a bit cautious and
-                                        // check this (also an additional safety
-                                        // check for the internal states of the
-                                        // library)
-       Assert ((face->child(0)->vertex_dof_index(2,0) ==
-                face->child(1)->vertex_dof_index(3,0)) &&
-               (face->child(0)->vertex_dof_index(2,0) ==
-                face->child(2)->vertex_dof_index(0,0)) &&
-               (face->child(0)->vertex_dof_index(2,0) ==
-                face->child(3)->vertex_dof_index(1,0)),
-               ExcInternalError());
-       next_index = 0;
-       for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
-         dofs_on_children[next_index++] = face->child(0)->vertex_dof_index(2,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++] = face->line(line)->child(0)->vertex_dof_index(1,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(1)->line(2)->dof_index(dof);
-       for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
-         dofs_on_children[next_index++] = face->child(2)->line(3)->dof_index(dof);
-       for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
-         dofs_on_children[next_index++] = face->child(3)->line(0)->dof_index(dof);
-
-                                        // dofs on the bordering lines
-       for (unsigned int line=0; line<4; ++line)
-         for (unsigned int child=0; child<2; ++child)
+       {
+         const DoFHandler<dim>::face_iterator face = cell->face(face);
+         
+                                          // fill the dofs indices. Use same
+                                          // enumeration scheme as in
+                                          // #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_children[next_index++] = face->line(line)->child(child)->dof_index(dof);
-       
-                                        // finally, for the dofs interior
-                                        // to the four child faces
-       for (unsigned int child=0; child<4; ++child)
+             dofs_on_mother[next_index++] = face->line(line)->dof_index(dof);
          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());
-
-       Assert (dofs_on_children.size() ==
-              fe.constraints().m(),
-              ExcDimensionMismatch(dofs_on_children.size(),
-                                     fe.constraints().m()));
-       Assert (dofs_on_mother.size() ==
-              fe.constraints().n(),
-              ExcDimensionMismatch(dofs_on_mother.size(),
-                                     fe.constraints().n()));
-
-                                        // 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));
-         };
-      };
+           dofs_on_mother[next_index++] = face->dof_index(dof);
+         Assert (next_index == dofs_on_mother.size(),
+                 ExcInternalError());
+         
+         next_index = 0;
+         Assert ((face->child(0)->vertex_dof_index(2,0) ==
+                  face->child(1)->vertex_dof_index(3,0)) &&
+                 (face->child(0)->vertex_dof_index(2,0) ==
+                  face->child(2)->vertex_dof_index(0,0)) &&
+                 (face->child(0)->vertex_dof_index(2,0) ==
+                  face->child(3)->vertex_dof_index(1,0)),
+                 ExcInternalError());
+         for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
+           dofs_on_children[next_index++]
+             = face->child(0)->vertex_dof_index(2,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++]
+               = face->line(line)->child(0)->vertex_dof_index(1,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(1)->line(2)->dof_index(dof);
+         for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
+           dofs_on_children[next_index++]
+             = face->child(2)->line(3)->dof_index(dof);
+         for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
+           dofs_on_children[next_index++]
+             = face->child(3)->line(0)->dof_index(dof);
+         
+                                          // 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++]
+                 = face->line(line)->child(child)->dof_index(dof);
+         
+                                          // 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());
+         
+                                          // 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));
+           };
+       };
 };
 
 #endif
 
 
+
 template <int dim, typename Number>
 void DoFTools::distribute_cell_to_dof_vector (const DoFHandler<dim> &dof_handler,
                                              const Vector<Number>  &cell_data,
@@ -818,42 +802,26 @@ DoFTools::extract_hanging_node_dofs (const DoFHandler<2> &dof_handler,
                                   // preset all values by false
   fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false);
 
-  const Triangulation<dim> &tria = dof_handler.get_tria();
   const FiniteElement<dim> &fe   = dof_handler.get_fe();
-  
-                                  // first mark all faces which are subject
-                                  // to 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
-  DoFHandler<dim>::line_iterator line = dof_handler.begin_line(),
-                                endl = dof_handler.end_line();
-  for (; line!=endl; ++line)
-    line->clear_user_flag ();
-  
-  Triangulation<dim>::active_cell_iterator cell = tria.begin_active(),
-                                          endc = tria.end();
+
+                                  // this function is similar to the
+                                  // make_sparsity_pattern function,
+                                  // see there for more information
+  DoFHandler<dim>::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()) 
-       cell->face(face)->set_user_flag();
-
-                                  // loop over all lines; only on lines
-                                  // can there be constraints.
-  for (line = dof_handler.begin_line(); line != endl; ++line)
-                                    // if dofs on this line are subject
-                                    // to constraints
-    if (line->user_flag_set() == true)
-      {        
-       for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
-         selected_dofs[line->child(0)->vertex_dof_index(1,dof)] = true;
-       
-       for (unsigned int child=0; child<2; ++child)
-         for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
-           selected_dofs[line->child(child)->dof_index(dof)] = true;
-      };
+       {
+         const DoFHandler<dim>::line_iterator line = cell->face(face);
+
+         for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
+           selected_dofs[line->child(0)->vertex_dof_index(1,dof)] = true;
+         
+         for (unsigned int child=0; child<2; ++child)
+           for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
+             selected_dofs[line->child(child)->dof_index(dof)] = true;
+       };
 };
 
 #endif
@@ -868,69 +836,60 @@ DoFTools::extract_hanging_node_dofs (const DoFHandler<3> &dof_handler,
                                     vector<bool>        &selected_dofs)
 {
   const unsigned int dim = 3;
-  const Triangulation<dim> &tria = dof_handler.get_tria();
+
+  Assert(selected_dofs.size() == dof_handler.n_dofs(),
+        ExcDimensionMismatch(selected_dofs.size(), dof_handler.n_dofs()));
+                                  // preset all values by false
+  fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false);
+
   const FiniteElement<dim> &fe   = dof_handler.get_fe();
   
-                                  // first mark all faces which are subject
-                                  // to 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
-  DoFHandler<dim>::face_iterator face = dof_handler.begin_face(),
-                                endf = dof_handler.end_face();
-  for (; face!=endf; ++face)
-    face->clear_user_flag ();
+                                  // this function is similar to the
+                                  // make_sparsity_pattern function,
+                                  // see there for more information
 
-  Triangulation<dim>::active_cell_iterator cell = tria.begin_active(),
-                                          endc = tria.end();
+  DoFHandler<dim>::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()) 
-       cell->face(face)->set_user_flag();
-       
-                                  // loop over all faces; only on faces
-                                  // can there be constraints.
-  for (face=dof_handler.begin_face(); face != endf; ++face)
-                                    // if dofs on this line are subject
-                                    // to constraints
-    if (face->user_flag_set() == true)
-      {
-       for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
-         selected_dofs[face->child(0)->vertex_dof_index(2,dof)] = true;
-       
-                                        // dof numbers on the centers of
-                                        // the lines bounding this face
-       for (unsigned int line=0; line<4; ++line)
+       {
+         const DoFHandler<dim>::face_iterator face = cell->face(face);
+         
          for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof)
-           selected_dofs[face->line(line)->child(0)->vertex_dof_index(1,dof)] = true;
-
-                                        // 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)
-         selected_dofs[face->child(0)->line(1)->dof_index(dof)] = true;
-       for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
-         selected_dofs[face->child(1)->line(2)->dof_index(dof)] = true;
-       for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
-         selected_dofs[face->child(2)->line(3)->dof_index(dof)] = true;
-       for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
-         selected_dofs[face->child(3)->line(0)->dof_index(dof)] = true;
-
-                                        // 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)
-             selected_dofs[face->line(line)->child(child)->dof_index(dof)] = true;
+           selected_dofs[face->child(0)->vertex_dof_index(2,dof)] = true;
+         
+                                          // 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)
+             selected_dofs[face->line(line)->child(0)->vertex_dof_index(1,dof)] = true;
+         
+                                          // 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)
+           selected_dofs[face->child(0)->line(1)->dof_index(dof)] = true;
+         for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
+           selected_dofs[face->child(1)->line(2)->dof_index(dof)] = true;
+         for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
+           selected_dofs[face->child(2)->line(3)->dof_index(dof)] = true;
+         for (unsigned int dof=0; dof<fe.dofs_per_line; ++dof)
+           selected_dofs[face->child(3)->line(0)->dof_index(dof)] = true;
+         
+                                          // 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)
+               selected_dofs[face->line(line)->child(child)->dof_index(dof)] = true;
        
-                                        // 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)
-           selected_dofs[face->child(child)->dof_index(dof)] = true;
-      };
+                                          // 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)
+             selected_dofs[face->child(child)->dof_index(dof)] = true;
+       };
 };
 
 #endif

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.