]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Re-implement the reserve_function for 2d. The existing one had the problem
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 17 Feb 2006 22:23:33 +0000 (22:23 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 17 Feb 2006 22:23:33 +0000 (22:23 +0000)
that it allocated space also for non-active cells, which is certainly not
something we want.

In the process, I have thrown out for the moment also the capability to deal
with continuous elements, but that is something that hasn't been working so
far anyway...

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

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

index 766d17561ecd8318ecfc9eb95b935861b1c003fd..c750866dc4e67c8081f51c455500d930b0fe6255 100644 (file)
@@ -1966,192 +1966,92 @@ namespace hp
   template <>
   void DoFHandler<2>::reserve_space ()
   {
-//TODO: This presently works only on the first FE of the FECollection. Fix this
     Assert (finite_elements != 0, ExcNoFESelected());
+    Assert (finite_elements->size() > 0, ExcNoFESelected());
     Assert (tria->n_levels() > 0, ExcInvalidTriangulation());
     Assert (tria->n_levels() == levels.size (), ExcInternalError ());
-  
-                                     // Backup the active_fe_indices vectors.
-                                     // The user might have put something
-                                     // important in them.
-    std::vector<std::vector<unsigned int> > active_fe_backup(levels.size ());
-    for (unsigned int i = 0; i < levels.size (); ++i)
-      std::swap (levels[i]->active_fe_indices, active_fe_backup[i]);
 
-                                     // delete all levels and set them up
-                                     // newly, since vectors are
-                                     // troublesome if you want to change
-                                     // their size
-    clear_space ();
+    Assert (finite_elements->max_dofs_per_face() == 0,
+            ExcMessage ("hp finite elements are presently only supported "
+                        "for discontinuous elements"));
+
+                                     // Release all space except the
+                                     // active_fe_indices field which
+                                     // we have to backup before
+    {
+      std::vector<std::vector<unsigned int> > active_fe_backup(levels.size ());
+      for (unsigned int level = 0; level<levels.size (); ++level)
+        std::swap (levels[level]->active_fe_indices, active_fe_backup[level]);
+
+                                       // delete all levels and set them up
+                                       // newly, since vectors are
+                                       // troublesome if you want to change
+                                       // their size
+      clear_space ();
   
-    vertex_dofs.resize(tria->n_vertices()*(*finite_elements)[0].dofs_per_vertex,
-                       invalid_dof_index);
-
-//TODO[?] Does not work for continuous FEs. Problem is at faces, which might have a different
-// number of DoFs due to the different active_fe_index of the adjacent cells.
-    for (unsigned int i=0; i<tria->n_levels(); ++i) 
+      for (unsigned int level=0; level<tria->n_levels(); ++level) 
+        {
+          levels.push_back (new internal::hp::DoFLevel<2>);
+          std::swap (active_fe_backup[level], levels[level]->active_fe_indices);
+        }
+    }
+
+                                     // count how much space we need
+                                     // on each level and set the
+                                     // dof_*_index_offset
+                                     // data. initially set the latter
+                                     // to an invalid index, and only
+                                     // later set it to something
+                                     // reasonable for active cells
+    for (unsigned int level=0; level<tria->n_levels(); ++level) 
       {
-        levels.push_back (new internal::hp::DoFLevel<2>);
-        std::swap (active_fe_backup[i], levels.back()->active_fe_indices);
-
-        levels.back()->dof_line_index_offset
-         = std::vector<unsigned int> (tria->n_raw_lines(i),
-                                      invalid_dof_index);
-        levels.back()->dof_quad_index_offset
-         = std::vector<unsigned int> (tria->n_raw_quads(i),
-                                      invalid_dof_index);
-
-                                         // Create table of active_fe_levels of cells
-                                         // adjacent to the lines.
-        std::vector<unsigned int> line_active_fe_indices
-         (tria->n_raw_lines(i) * 2, 
-          invalid_dof_index);
-
-        for (unsigned int j = 0; j < tria->n_raw_quads(i); ++j)
-          {
-            const unsigned int cell_fe_index
-             = levels.back ()->active_fe_indices[j];
-
-                                             // Referring to the following diagram from
-                                             // documentation of the triangulation class,
-                                             // it can be seen where the indices in the
-                                             // following lines of code come from.
-                                             //          2
-                                             //      3--->---2
-                                             //      |       |
-                                             //     3^       ^1
-                                             //      |       |
-                                             //      0--->---1
-                                             //          0
-                                             // Face 0 and 1 have the cell to the left,
-                                             // while Face 2 and 3 have the cell to
-                                             // the right (if the line points upward).
-            line_active_fe_indices[tria->levels[i]->
-                                   quads.quads[j].line(0) * 2] = cell_fe_index;
-            line_active_fe_indices[tria->levels[i]->
-                                   quads.quads[j].line(1) * 2] = cell_fe_index;
-            line_active_fe_indices[tria->levels[i]->
-                                   quads.quads[j].line(2) * 2 + 1] = cell_fe_index;
-            line_active_fe_indices[tria->levels[i]->
-                                   quads.quads[j].line(3) * 2 + 1] = cell_fe_index;
-          }
-
-        unsigned int dofs_for_lines = 0;
-        for (unsigned int j = 0; j < tria->n_raw_lines(i); ++j)
-          {
-            levels.back()->dof_line_index_offset[j] = dofs_for_lines;
-
-// The "line_dofs" field in the hpDoFLevel class for 2D gets a
-// slightly different meaning in 2D. Actually in 2D one line can have
-// two different active_fe_indices. Therefore it is not sufficient to
-// reserve some space for the DoFs on one active_fe_index. Instead
-// the DoF indices for both active_fe_indices will be stored in the
-// "line_dofs" field. To differentiate between these field later,
-// a simple linked list based will be maintained to allow access to
-// the correct DoFs:
-// 0: active_fe_index of line, to which the following DoFs belong.
-// 1: index of next active_fe_index for this line. Value: (x + 1)
-// 2-x: DoFs (will be filled out by distribute_dofs_on_cell
-// x+1: active_fe_index of line, to which the following DoFs belong.
-// x+2: index to next active_fe_index for this line. Value: (y + 1)
-// ...
-// In 3D, the same datastructure is used but can have up to 4 different
-// active_fe_indices on one line.
-            unsigned int active_fe1 = line_active_fe_indices[j*2],
-                         active_fe2 = line_active_fe_indices[j*2+1];
-
-                                             // Check for boundary lines, where clearly one of the line indices
-                                             // is missing.
-            if (active_fe1 == invalid_dof_index)
-             active_fe1 = active_fe2;
-            if (active_fe2 == invalid_dof_index)
-             active_fe2 = active_fe1;
-
-                                            // at this point, all
-                                            // indices should be ok
-                                            // (but weren't at one
-                                            // point, thus the
-                                            // crash_01 testcase)
-           Assert (active_fe1 < finite_elements->size(),
-                   ExcInternalError());
-           Assert (active_fe2 < finite_elements->size(),
-                   ExcInternalError());
-           
-            if (active_fe1 == active_fe2)
-             dofs_for_lines += (*finite_elements)[active_fe1].dofs_per_line + 2;
-            else
-             dofs_for_lines += (*finite_elements)[active_fe1].dofs_per_line +
-                                (*finite_elements)[active_fe2].dofs_per_line + 4;
-          }
-
-        unsigned int dofs_for_quads = 0;
-        for (unsigned int j = 0; j < tria->n_raw_quads(i); ++j)
+        levels[level]->dof_quad_index_offset
+          = std::vector<unsigned int> (tria->n_raw_quads(level),
+                                       invalid_dof_index);
+
+//TODO: do the same for line dofs        
+        unsigned int next_free_quad_dof = 0;
+        for (active_cell_iterator cell=begin_active(level);
+             cell!=end_active(level); ++cell)
+          if (!cell->has_children())
           {
-            levels.back()->dof_quad_index_offset[j] = dofs_for_quads;
-            dofs_for_quads += (*finite_elements)[
-             levels.back ()->active_fe_indices[j]].dofs_per_quad;
+            levels[level]->dof_quad_index_offset[cell->index()] = next_free_quad_dof;
+            next_free_quad_dof +=
+              (*finite_elements)[cell->active_fe_index()].dofs_per_quad;
           }
 
-        levels.back()->line_dofs = std::vector<unsigned int> (dofs_for_lines,
-                                                              invalid_dof_index);
-        levels.back()->quad_dofs = std::vector<unsigned int> (dofs_for_quads,
+        levels[level]->quad_dofs = std::vector<unsigned int> (next_free_quad_dof,
                                                               invalid_dof_index);
+      }
 
-                                         // As we already have the
-                                         // active_fe_indices for the
-                                         // lines, it is probably a
-                                         // good idea to directly use
-                                         // those to create the linked
-                                         // list structure within the
-                                         // line_dofs field.
-        for (unsigned int j = 0; j < tria->n_raw_lines(i); ++j)
-          {
-            dofs_for_lines = levels.back()->dof_line_index_offset[j];
-
-            levels.back()->dof_line_index_offset[j] = dofs_for_lines;
-            unsigned int active_fe1 = line_active_fe_indices[j*2],
-                         active_fe2 = line_active_fe_indices[j*2+1];
-
-                                             // Check for boundary
-                                             // lines, where clearly
-                                             // one of the line
-                                             // indices is missing.
-            if (active_fe1 == invalid_dof_index)
-             active_fe1 = active_fe2;
-            if (active_fe2 == invalid_dof_index)
-             active_fe2 = active_fe1;
-
-                                            // at this point, all
-                                            // indices should be ok
-                                            // (but weren't at one
-                                            // point, thus the
-                                            // crash_01 testcase)
-           Assert (active_fe1 < finite_elements->size(),
-                   ExcInternalError());
-           Assert (active_fe2 < finite_elements->size(),
-                   ExcInternalError());
-
-                                             // Now create prepare
-                                             // linked list, with
-                                             // either 1 or two
-                                             // entries.
-            levels.back()->line_dofs[dofs_for_lines] = active_fe1;
-            if (active_fe1 == active_fe2)
-             levels.back()->line_dofs[dofs_for_lines + 1] = 
-                invalid_dof_index;
-            else
-              {
-                const unsigned int start_next_fe =
-                 dofs_for_lines + (*finite_elements)[active_fe1].dofs_per_line + 2;
-
-                levels.back()->line_dofs[dofs_for_lines + 1] = start_next_fe;
 
-                levels.back()->line_dofs[start_next_fe] = active_fe2;
-                levels.back()->line_dofs[start_next_fe + 1] = 
-                 invalid_dof_index;
-              }
-          }
-      };
+                                     // safety check: make sure that
+                                     // the number of DoFs we
+                                     // allocated is actually correct
+                                     // (above we have also set the
+                                     // dof_*_index_offset field, so
+                                     // we couldn't use this simpler
+                                     // algorithm)
+#ifdef DEBUG
+    for (unsigned int level=0; level<tria->n_levels(); ++level)
+      {
+        unsigned int n_quad_dofs = 0;
+        for (cell_iterator cell=begin_active(level);
+             cell!=end_active(level); ++cell)
+          if (!cell->has_children())
+            n_quad_dofs +=
+              (*finite_elements)[cell->active_fe_index()].dofs_per_quad;
+        
+        Assert (levels[level]->quad_dofs.size() == n_quad_dofs, ExcInternalError());
+        Assert (static_cast<unsigned int>
+                (std::count (levels[level]->dof_quad_index_offset.begin(),
+                             levels[level]->dof_quad_index_offset.end(),
+                             invalid_dof_index))
+                ==
+                tria->n_raw_quads(level) - tria->n_active_quads(level),
+                ExcInternalError());
+      }
+#endif
   }
 
 #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.