]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement unification of DoFs
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 28 Jul 2006 17:17:39 +0000 (17:17 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 28 Jul 2006 17:17:39 +0000 (17:17 +0000)
git-svn-id: https://svn.dealii.org/trunk@13475 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 38f8d4aed0d186c35e4314cd62825e4a43711d45..d4b3315b6202cb6aae045c13a165405263402882 100644 (file)
@@ -17,6 +17,7 @@
 
 #include <base/config.h>
 #include <base/exceptions.h>
+#include <base/template_constraints.h>
 #include <base/smartpointer.h>
 #include <dofs/function_map.h>
 #include <dofs/dof_iterator_selector.h>
@@ -175,7 +176,7 @@ namespace hp
       virtual void clear ();
     
                                        /**
-                                        * Actually do the renumbering based on
+                                        * Renumber degrees of freedom based on
                                         * a list of new dof numbers for all the
                                         * dofs.
                                         *
@@ -185,13 +186,25 @@ namespace hp
                                         * indices after renumbering in the
                                         * order of the old indices.
                                         *
-                                        * This function is called by the
-                                        * @p renumber_dofs function after computing
-                                        * the ordering of the degrees of freedom.
-                                        * However, you can call this function
-                                        * yourself, which is necessary if a user
-                                        * wants to implement an ordering scheme
-                                        * herself, for example downwind numbering.
+                                        * This function is called by
+                                        * the functions in
+                                        * DoFRenumbering function
+                                        * after computing the ordering
+                                        * of the degrees of freedom.
+                                        * However, you can call this
+                                        * function yourself, which is
+                                        * necessary if a user wants to
+                                        * implement an ordering scheme
+                                        * herself, for example
+                                        * downwind numbering.
+                                       *
+                                       * The @p new_number array must
+                                       * have a size equal to the
+                                       * number of degrees of
+                                       * freedom. Each entry must
+                                       * state the new global DoF
+                                       * number of the degree of
+                                       * freedom referenced.
                                         */
       void renumber_dofs (const std::vector<unsigned int> &new_numbers);
 
@@ -1132,7 +1145,44 @@ namespace hp
                                        */
       void
       compute_quad_dof_identities (std::vector<unsigned int> &new_dof_indices) const;
-      
+
+                                      /**
+                                       * Renumber the objects with
+                                       * the given and all lower
+                                       * structural dimensions,
+                                       * i.e. renumber vertices by
+                                       * giving a template argument
+                                       * of zero to the int2type
+                                       * argument, lines and vertices
+                                       * with one, etc.
+                                       *
+                                       * Note that in contrast to the
+                                       * public renumber_dofs()
+                                       * function, these internal
+                                       * functions do not ensure that
+                                       * the new DoFs are
+                                       * contiguously numbered. The
+                                       * function may therefore also
+                                       * be used to assign different
+                                       * DoFs the same number, for
+                                       * example to unify hp DoFs
+                                       * corresponding to different
+                                       * finite elements but
+                                       * co-located on the same
+                                       * entity.
+                                       */
+      void renumber_dofs_internal (const std::vector<unsigned int> &new_numbers,
+                                  internal::int2type<0>);
+
+      void renumber_dofs_internal (const std::vector<unsigned int> &new_numbers,
+                                  internal::int2type<1>);
+
+      void renumber_dofs_internal (const std::vector<unsigned int> &new_numbers,
+                                  internal::int2type<2>);
+
+      void renumber_dofs_internal (const std::vector<unsigned int> &new_numbers,
+                                  internal::int2type<3>);
+
                                        /**
                                         * Space to store the DoF
                                         * numbers for the different
@@ -1351,9 +1401,6 @@ namespace hp
   template <>
   unsigned int DoFHandler<3>::distribute_dofs_on_cell (active_cell_iterator &cell,
                                                          unsigned int          next_free_dof);
-  template <> void DoFHandler<1>::renumber_dofs (const std::vector<unsigned int> &new_numbers);
-  template <> void DoFHandler<2>::renumber_dofs (const std::vector<unsigned int> &new_numbers);
-  template <> void DoFHandler<3>::renumber_dofs (const std::vector<unsigned int> &new_numbers);
   template <> void DoFHandler<1>::reserve_space ();
   template <> void DoFHandler<2>::reserve_space ();
   template <> void DoFHandler<3>::reserve_space ();
index 51005fb7e5bcae785367d7b4c7db0243145ce7aa..61c71f5b92eacd9b75f6fd45eee49974b4b0720e 100644 (file)
@@ -2444,22 +2444,54 @@ namespace hp
                                     // faces and other
                                     // lower-dimensional objects
                                     // where elements come together
-    std::vector<unsigned int> new_dof_indices (used_dofs,
-                                              deal_II_numbers::invalid_unsigned_int);
-    compute_vertex_dof_identities (new_dof_indices);
-    compute_line_dof_identities (new_dof_indices);
-    compute_quad_dof_identities (new_dof_indices);
+    std::vector<unsigned int>
+      constrained_indices (used_dofs, deal_II_numbers::invalid_unsigned_int);
+    compute_vertex_dof_identities (constrained_indices);
+    compute_line_dof_identities (constrained_indices);
+    compute_quad_dof_identities (constrained_indices);
+
+    const unsigned int used_dofs_before = used_dofs;
+
+                                    // loop over all dofs and assign
+                                    // new numbers to those which are
+                                    // not constrained
+    std::vector<unsigned int>
+      new_dof_indices (used_dofs, deal_II_numbers::invalid_unsigned_int);
+    unsigned int next_free_dof = 0;
+    for (unsigned int i=0; i<used_dofs; ++i)
+      if (constrained_indices[i] == deal_II_numbers::invalid_unsigned_int)
+       {
+         new_dof_indices[i] = next_free_dof;
+         ++next_free_dof;
+       }
+    
+                                    // then loop over all those that
+                                    // are constrained and record the
+                                    // new dof number for those:
+    for (unsigned int i=0; i<used_dofs; ++i)
+      if (constrained_indices[i] != deal_II_numbers::invalid_unsigned_int)
+       {
+         Assert (new_dof_indices[constrained_indices[i]] !=
+                 deal_II_numbers::invalid_unsigned_int,
+                 ExcInternalError());
+         
+         new_dof_indices[i] = new_dof_indices[constrained_indices[i]];
+       }
+
+    for (unsigned int i=0; i<used_dofs; ++i)
+      {
+       Assert (new_dof_indices[i] != deal_II_numbers::invalid_unsigned_int,
+               ExcInternalError());
+       Assert (new_dof_indices[i] < next_free_dof,
+               ExcInternalError());
+      }
     
-//     if (new_dof_indices.size() -
-//     std::count (new_dof_indices.begin(),
-//                 new_dof_indices.end(), deal_II_numbers::invalid_unsigned_int)
-//     != 0)
-//       std::cout << "ELIMINATING "
-//             << (new_dof_indices.size() -
-//                 std::count (new_dof_indices.begin(),
-//                             new_dof_indices.end(), deal_II_numbers::invalid_unsigned_int))
-//             << " DOFS"
-//             << std::endl;
+                                    // finally, do the renumbering
+                                    // and set the number of actually
+                                    // used dof indices
+    renumber_dofs_internal (new_dof_indices, internal::int2type<dim>());
+
+    used_dofs = next_free_dof;
     
     
                                      // finally restore the user flags
@@ -2654,10 +2686,8 @@ namespace hp
 #endif
 
 
-#if deal_II_dimension == 1
-
-  template <>
-  void DoFHandler<1>::renumber_dofs (const std::vector<unsigned int> &new_numbers)
+  template <int dim>
+  void DoFHandler<dim>::renumber_dofs (const std::vector<unsigned int> &new_numbers)
   {
     Assert (new_numbers.size() == n_dofs(), ExcRenumberingIncomplete());
 #ifdef DEBUG
@@ -2671,128 +2701,135 @@ namespace hp
         unsigned int                         i = 0;
         for (; p!=tmp.end(); ++p, ++i)
           Assert (*p == i, ExcNewNumbersNotConsecutive(i));
-      };
+      }
 #endif
 
-                                     // note that we can not use cell iterators
-                                     // in this function since then we would
-                                     // renumber the dofs on the interface of
-                                     // two cells more than once. Anyway, this
-                                     // way it's not only more correct but also
-                                     // faster; note, however, that dof numbers
-                                     // may be invalid_dof_index, namely when the appropriate
-                                     // vertex/line/etc is unused
-    for (std::vector<unsigned int>::iterator i=vertex_dofs.begin(); i!=vertex_dofs.end(); ++i)
-      if (*i != invalid_dof_index)
-        *i = new_numbers[*i];
-  
-    for (unsigned int level=0; level<levels.size(); ++level) 
-      for (std::vector<unsigned int>::iterator i=levels[level]->lines.dofs.begin();
-           i!=levels[level]->lines.dofs.end(); ++i)
-        if (*i != invalid_dof_index)
-          *i = new_numbers[*i];
+    renumber_dofs_internal (new_numbers, internal::int2type<dim>());
   }
 
-#endif
-
 
-#if deal_II_dimension == 2
 
-  template <>
-  void DoFHandler<2>::renumber_dofs (const std::vector<unsigned int> &new_numbers)
+  template <int dim>
+  void
+  DoFHandler<dim>::
+  renumber_dofs_internal (const std::vector<unsigned int> &new_numbers,
+                         internal::int2type<0>)
   {
     Assert (new_numbers.size() == n_dofs(), ExcRenumberingIncomplete());
-#ifdef DEBUG
-                                     // assert that the new indices are
-                                     // consecutively numbered
-    if (true)
-      {
-        std::vector<unsigned int> tmp(new_numbers);
-        std::sort (tmp.begin(), tmp.end());
-        std::vector<unsigned int>::const_iterator p = tmp.begin();
-        unsigned int                         i = 0;
-        for (; p!=tmp.end(); ++p, ++i)
-          Assert (*p == i, ExcNewNumbersNotConsecutive(i));
-      };
-#endif
 
-                                     // note that we can not use cell iterators
-                                     // in this function since then we would
-                                     // renumber the dofs on the interface of
-                                     // two cells more than once. Anyway, this
-                                     // way it's not only more correct but also
-                                     // faster; note, however, that dof numbers
-                                     // may be invalid_dof_index, namely when the appropriate
-                                     // vertex/line/etc is unused
-    for (std::vector<unsigned int>::iterator i=vertex_dofs.begin(); i!=vertex_dofs.end(); ++i)
-      if (*i != invalid_dof_index)
-        *i = new_numbers[*i];
-  
-    for (unsigned int level=0; level<levels.size(); ++level) 
+    for (unsigned int vertex_index=0; vertex_index<get_tria().n_vertices();
+        ++vertex_index)
       {
-        for (std::vector<unsigned int>::iterator i=levels[level]->quads.dofs.begin();
-             i!=levels[level]->quads.dofs.end(); ++i)
-          if (*i != invalid_dof_index)
-            *i = new_numbers[*i];
-      }
-    for (std::vector<unsigned int>::iterator i=faces->lines.dofs.begin();
-        i!=faces->lines.dofs.end(); ++i)
-      if (*i != invalid_dof_index)
-       *i = new_numbers[*i];
+       const unsigned int n_active_fe_indices
+         = n_active_vertex_fe_indices (vertex_index);
 
+       for (unsigned int f=0; f<n_active_fe_indices; ++f)
+         {
+           const unsigned int fe_index
+             = nth_active_vertex_fe_index (vertex_index, f);
+
+           for (unsigned int d=0; d<(*finite_elements)[fe_index].dofs_per_vertex; ++d)
+             set_vertex_dof_index (vertex_index,
+                                   fe_index,
+                                   d,
+                                   new_numbers[get_vertex_dof_index(vertex_index,
+                                                                    fe_index,
+                                                                    d)]);
+         }
+      }
   }
 
-#endif
 
 
-#if deal_II_dimension == 3
+  template <int dim>
+  void
+  DoFHandler<dim>::
+  renumber_dofs_internal (const std::vector<unsigned int> &new_numbers,
+                         internal::int2type<1>)
+  {
+    Assert (new_numbers.size() == n_dofs(), ExcRenumberingIncomplete());
 
-  template <>
-  void DoFHandler<3>::renumber_dofs (const std::vector<unsigned int> &new_numbers)
+    renumber_dofs_internal (new_numbers, internal::int2type<0>());
+    
+    for (line_iterator line=begin_line(); line!=end_line(); ++line)
+      {
+       const unsigned int n_active_fe_indices
+         = line->n_active_fe_indices ();
+
+       for (unsigned int f=0; f<n_active_fe_indices; ++f)
+         {
+           const unsigned int fe_index
+             = line->nth_active_fe_index (f);
+
+           for (unsigned int d=0; d<(*finite_elements)[fe_index].dofs_per_line; ++d)
+             line->set_dof_index (d,
+                                  new_numbers[line->dof_index(d,fe_index)],
+                                  fe_index);
+         }
+      }
+  }
+
+
+#if deal_II_dimension >= 2
+
+  template <int dim>
+  void
+  DoFHandler<dim>::
+  renumber_dofs_internal (const std::vector<unsigned int> &new_numbers,
+                         internal::int2type<2>)
   {
     Assert (new_numbers.size() == n_dofs(), ExcRenumberingIncomplete());
-#ifdef DEBUG
-                                     // assert that the new indices are
-                                     // consecutively numbered
-    if (true)
+
+    renumber_dofs_internal (new_numbers, internal::int2type<1>());
+    
+    for (quad_iterator quad=begin_quad(); quad!=end_quad(); ++quad)
       {
-        std::vector<unsigned int> tmp(new_numbers);
-        std::sort (tmp.begin(), tmp.end());
-        std::vector<unsigned int>::const_iterator p = tmp.begin();
-        unsigned int                              i = 0;
-        for (; p!=tmp.end(); ++p, ++i)
-          Assert (*p == i, ExcNewNumbersNotConsecutive(i));
-      };
+       const unsigned int n_active_fe_indices
+         = quad->n_active_fe_indices ();
+
+       for (unsigned int f=0; f<n_active_fe_indices; ++f)
+         {
+           const unsigned int fe_index
+             = quad->nth_active_fe_index (f);
+
+           for (unsigned int d=0; d<(*finite_elements)[fe_index].dofs_per_quad; ++d)
+             quad->set_dof_index (d,
+                                  new_numbers[quad->dof_index(d,fe_index)],
+                                  fe_index);
+         }
+      }
+  }
+
 #endif
 
-                                     // note that we can not use cell iterators
-                                     // in this function since then we would
-                                     // renumber the dofs on the interface of
-                                     // two cells more than once. Anyway, this
-                                     // way it's not only more correct but also
-                                     // faster; note, however, that dof numbers
-                                     // may be invalid_dof_index, namely when the appropriate
-                                     // vertex/line/etc is unused
-    for (std::vector<unsigned int>::iterator i=vertex_dofs.begin(); i!=vertex_dofs.end(); ++i)
-      if (*i != invalid_dof_index)
-        *i = new_numbers[*i];
-  
-    for (unsigned int level=0; level<levels.size(); ++level) 
+#if deal_II_dimension >= 3
+
+  template <int dim>
+  void
+  DoFHandler<dim>::
+  renumber_dofs_internal (const std::vector<unsigned int> &new_numbers,
+                         internal::int2type<3>)
+  {
+    Assert (new_numbers.size() == n_dofs(), ExcRenumberingIncomplete());
+
+    renumber_dofs_internal (new_numbers, internal::int2type<2>());
+    
+    for (hex_iterator hex=begin_hex(); hex!=end_hex(); ++hex)
       {
-        for (std::vector<unsigned int>::iterator i=levels[level]->hexes.dofs.begin();
-             i!=levels[level]->hexes.dofs.end(); ++i)
-          if (*i != invalid_dof_index)
-            *i = new_numbers[*i];
+       const unsigned int n_active_fe_indices
+         = hex->n_active_fe_indices ();
+
+       for (unsigned int f=0; f<n_active_fe_indices; ++f)
+         {
+           const unsigned int fe_index
+             = hex->nth_active_fe_index (f);
+
+           for (unsigned int d=0; d<(*finite_elements)[fe_index].dofs_per_hex; ++d)
+             hex->set_dof_index (d,
+                                 new_numbers[hex->dof_index(d,fe_index)],
+                                 fe_index);
+         }
       }
-    
-    for (std::vector<unsigned int>::iterator i=faces->lines.dofs.begin();
-        i!=faces->lines.dofs.end(); ++i)
-      if (*i != invalid_dof_index)
-       *i = new_numbers[*i];
-    for (std::vector<unsigned int>::iterator i=faces->quads.dofs.begin();
-        i!=faces->quads.dofs.end(); ++i)
-      if (*i != invalid_dof_index)
-       *i = new_numbers[*i];    
   }
 
 #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.