]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix some really nasty bugs which came to the surface with the implementation of coars...
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 17 Jun 1998 09:47:50 +0000 (09:47 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 17 Jun 1998 09:47:50 +0000 (09:47 +0000)
git-svn-id: https://svn.dealii.org/trunk@406 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/dofs/dof_handler.cc
deal.II/deal.II/source/grid/tria.cc

index bcd59f9096cfa0f4d57ce9978c15719031fb27f7..bc9b73ee0ea44f21107abf876d3012bce7540c00 100644 (file)
@@ -779,6 +779,7 @@ void DoFHandler<dim>::distribute_dofs (const FiniteElementBase<dim> &fe) {
   unsigned int next_free_dof = 0;   
   active_cell_iterator cell = begin_active(),
                       endc = end();
+
   for (; cell != endc; ++cell) 
     next_free_dof = distribute_dofs_on_cell (cell, next_free_dof);
   
@@ -857,23 +858,19 @@ unsigned int DoFHandler<2>::distribute_dofs_on_cell (active_cell_iterator &cell,
          cell->set_vertex_dof_index (vertex, d, next_free_dof++);
     
                                   // for the four sides
-  for (unsigned int side=0; side<GeometryInfo<2>::faces_per_cell; ++side)
-    {
-      line_iterator line = cell->line(side);
-
-                                      // distribute dofs if necessary
-      if ((selected_fe->dofs_per_line > 0) &&
+  if (selected_fe->dofs_per_line > 0)
+    for (unsigned int side=0; side<GeometryInfo<2>::faces_per_cell; ++side)
+      {
+       line_iterator line = cell->line(side);
+       
+                                        // distribute dofs if necessary:
                                         // check whether line dof is already
                                         // numbered (check only first dof)
-         (line->dof_index(0) == -1))
-                                        // if not: distribute dofs
-       for (unsigned int d=0; d<selected_fe->dofs_per_line; ++d)
-         line->set_dof_index (d, next_free_dof++);         
-
-                                      // note if line is subject to constraints
-      if (line->child_index(0) != -1) 
-       line->set_user_flag ();
-    };
+       if (line->dof_index(0) == -1)
+                                          // if not: distribute dofs
+         for (unsigned int d=0; d<selected_fe->dofs_per_line; ++d)
+           line->set_dof_index (d, next_free_dof++);       
+      };
   
 
                                       // dofs of quad
@@ -1156,6 +1153,7 @@ void DoFHandler<2>::make_constraint_matrix (ConstraintMatrix &constraints) const
     for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
       if (cell->face(face)->has_children()) 
        cell->face(face)->set_user_flag();
+         
   
 
   
@@ -1834,16 +1832,17 @@ void DoFHandler<1>::reserve_space () {
     delete levels[i];
   levels.resize (0);
 
-  vertex_dofs.erase (vertex_dofs.begin(), vertex_dofs.end());
-  vertex_dofs.resize (tria->vertices.size(), -1);
+  vertex_dofs = vector<int>(tria->vertices.size()*
+                           selected_fe->dofs_per_vertex,
+                           -1);
     
   for (unsigned int i=0; i<tria->n_levels(); ++i) 
     {
       levels.push_back (new DoFLevel<1>);
 
-      levels.back()->line_dofs.resize (tria->levels[i]->lines.lines.size() *
-                                      selected_fe->dofs_per_line,
-                                      -1);
+      levels.back()->line_dofs = vector<int>(tria->levels[i]->lines.lines.size() *
+                                            selected_fe->dofs_per_line,
+                                            -1);
     };
 };
 
@@ -1865,20 +1864,19 @@ void DoFHandler<2>::reserve_space () {
     delete levels[i];
   levels.resize (0);
 
-  vertex_dofs.erase (vertex_dofs.begin(), vertex_dofs.end());
-  vertex_dofs.resize (tria->vertices.size(),
-                     -1);
-    
+  vertex_dofs = vector<int>(tria->vertices.size()*
+                           selected_fe->dofs_per_vertex,
+                           -1);
   for (unsigned int i=0; i<tria->n_levels(); ++i) 
     {
       levels.push_back (new DoFLevel<2>);
 
-      levels.back()->line_dofs.resize (tria->levels[i]->lines.lines.size() *
-                                      selected_fe->dofs_per_line,
-                                      -1);
-      levels.back()->quad_dofs.resize (tria->levels[i]->quads.quads.size() *
-                                      selected_fe->dofs_per_quad,
-                                      -1);
+      levels.back()->line_dofs = vector<int> (tria->levels[i]->lines.lines.size() *
+                                             selected_fe->dofs_per_line,
+                                             -1);
+      levels.back()->quad_dofs = vector<int> (tria->levels[i]->quads.quads.size() *
+                                             selected_fe->dofs_per_quad,
+                                             -1);
     };
 };
 
index e3a6d5f2038d6214e1153df15861e1ac33700ac7..7ab1b8987fbfe57bdceee146044d3e42f78a6bf6 100644 (file)
@@ -605,7 +605,7 @@ void Triangulation<2>::create_hyper_L (const double a, const double b) {
 template <>
 void Triangulation<2>::create_hyper_ball (const Point<2> &p, const double radius) {
   const unsigned int dim=2;
-  const double a = 0.453;         // equilibrate cell sizes
+  const double a = 1./(1+sqrt(2)/2);         // equilibrate cell sizes
   const Point<dim> vertices[8] = { p+Point<dim>(-1,-1)*(radius/sqrt(2)),
                                   p+Point<dim>(+1,-1)*(radius/sqrt(2)),
                                   p+Point<dim>(-1,-1)*(radius/sqrt(2)*a),
@@ -1815,7 +1815,26 @@ unsigned int Triangulation<dim>::n_active_quads (const unsigned int level) const
 
 template <int dim>
 unsigned int Triangulation<dim>::n_levels () const {
-  return levels.size();
+  if (levels.size() == 0)
+    return 0;
+                                  // check whether there are
+                                  // cells on the highest levels
+                                  // (there need not be, since they
+                                  // might all have been coarsened
+                                  // away)
+  raw_cell_iterator cell = last_raw (levels.size()-1),
+                   endc = end();
+  for (; cell!=endc; --cell)
+    if (cell->used())
+                                      // return level of most refined
+                                      // existing cell (+1 because of
+                                      // counting conventions)
+      return cell->level()+1;
+
+                                  // no cells at all?
+  Assert (false, ExcInternalError());
+                                  // just to make the compiler happy:
+  return 0;
 };
 
 
@@ -1922,6 +1941,14 @@ void Triangulation<dim>::refine (const dVector &criteria,
   Assert (criteria.size() == n_active_cells(),
          ExcInvalidVectorSize(criteria.size(), n_active_cells()));
 
+                                  // nothing to do; especially we
+                                  // do not want to flag with zero
+                                  // error since then we may get
+                                  // into conflict with coarsening
+                                  // in some cases
+  if (threshold==0)
+    return;
+  
   active_cell_iterator cell = begin_active();
   const unsigned int n_cells = criteria.size();
   
@@ -2026,8 +2053,6 @@ Triangulation<dim>::refine_and_coarsen_fixed_fraction (const dVector &criteria,
   else
     bottom_threshold = sqrt(*p - *(p-1));
 
-  Assert (bottom_threshold<=top_threshold, ExcInternalError());
-
                                   // in some rare cases it may happen that
                                   // both thresholds are the same (e.g. if
                                   // there are many cells with the same
@@ -2038,10 +2063,16 @@ Triangulation<dim>::refine_and_coarsen_fixed_fraction (const dVector &criteria,
                                   // which the indicator equals the
                                   // thresholds. This is forbidden, however.
                                   //
-                                  // In that case we arbitrarily reduce the
-                                  // bottom threshold by one permille.
-  if (bottom_threshold==top_threshold)
-    bottom_threshold *= 0.999;
+                                  // In some rare cases with very few cells
+                                  // we also could get integer round off
+                                  // errors and get problems with
+                                  // the top and bottom fractions.
+                                  //
+                                  // In these case we arbitrarily reduce the
+                                  // bottom threshold by one permille below
+                                  // the top threshold
+  if (bottom_threshold>=top_threshold)
+    bottom_threshold = 0.999*top_threshold;
   
                                   // actually flag cells
   refine (criteria, top_threshold);
@@ -2111,11 +2142,10 @@ void Triangulation<1>::execute_refinement () {
 
                                       // count number of used cells on
                                       // the next higher level
-      unsigned int used_cells = 0;
-      count_if (levels[level+1]->lines.used.begin(),
-               levels[level+1]->lines.used.end(),
-               bind2nd (equal_to<bool>(), true),
-               used_cells);
+      const unsigned int used_cells
+       =  count_if (levels[level+1]->lines.used.begin(),
+                    levels[level+1]->lines.used.end(),
+                    bind2nd (equal_to<bool>(), true));
 
                                       // reserve space for the used_cells
                                       // cells already existing on the next
@@ -2137,12 +2167,18 @@ void Triangulation<1>::execute_refinement () {
 
                                   // add to needed vertices how
                                   // many vertices are already in use
-  count_if (vertices_used.begin(), vertices_used.end(),
-           bind2nd (equal_to<bool>(), true),
-           needed_vertices);
-
-  vertices.resize (needed_vertices, Point<1>());
-  vertices_used.resize (needed_vertices, false);
+  needed_vertices += count_if (vertices_used.begin(), vertices_used.end(),
+                              bind2nd (equal_to<bool>(), true));
+                                  // if we need more vertices: create them,
+                                  // if not: leave the array as is, since
+                                  // shrinking is not really possible because
+                                  // some of the vertices at the end may be
+                                  // in use
+  if (needed_vertices > vertices.size())
+    {
+      vertices.resize (needed_vertices, Point<1>());
+      vertices_used.resize (needed_vertices, false);
+    };
 
 
                                   // Do REFINEMENT
@@ -2401,21 +2437,12 @@ void Triangulation<2>::execute_refinement () {
              };
          };
       
-                                      // count number of used lines on
-                                      // the next higher level
-      unsigned int used_lines = 0;
-      count_if (levels[level+1]->lines.used.begin(),
-               levels[level+1]->lines.used.end(),
-               bind2nd (equal_to<bool>(), true),
-               used_lines);
-
                                       // count number of used cells on
                                       // the next higher level
-      unsigned int used_cells = 0;
-      count_if (levels[level+1]->quads.used.begin(),
+      const unsigned int used_cells
+       = count_if (levels[level+1]->quads.used.begin(),
                levels[level+1]->quads.used.end(),
-               bind2nd (equal_to<bool>(), true),
-               used_cells);
+               bind2nd (equal_to<bool>(), true));
       
       
                                       // reserve space for the used_cells
@@ -2438,14 +2465,18 @@ void Triangulation<2>::execute_refinement () {
 
                                   // add to needed vertices how
                                   // many vertices are already in use
-  count_if (vertices_used.begin(), vertices_used.end(),
-           bind2nd (equal_to<bool>(), true),
-           needed_vertices);
-
-                                  // reserve enough space for all vertices
-  vertices.resize (needed_vertices, Point<2>());
-  vertices_used.resize (needed_vertices, false);
-
+  needed_vertices += count_if (vertices_used.begin(), vertices_used.end(),
+                              bind2nd (equal_to<bool>(), true));
+                                  // if we need more vertices: create them,
+                                  // if not: leave the array as is, since
+                                  // shrinking is not really possible because
+                                  // some of the vertices at the end may be
+                                  // in use
+  if (needed_vertices > vertices.size())
+    {
+      vertices.resize (needed_vertices, Point<dim>());
+      vertices_used.resize (needed_vertices, false);
+    };
 
 
                                   // Do REFINEMENT
@@ -3537,14 +3568,10 @@ void TriangulationLevel<0>::monitor_memory (const unsigned int true_dimension) c
 
 
 void TriangulationLevel<1>::reserve_space (const unsigned int new_lines) {
-  vector<bool>::iterator u = lines.used.begin(),
-                        e = lines.used.end();
-
-  unsigned int used_lines = 0;
-  for (; u!=e; ++u)
-    ++used_lines;
-
-  const unsigned int new_size = used_lines + new_lines;
+  const unsigned int new_size = new_lines +
+                               count_if (lines.used.begin(),
+                                         lines.used.end(),
+                                         bind2nd (equal_to<bool>(), true));
 
                                   // same as in #reserve_space<0>#: only
                                   // allocate space if necessary
@@ -3616,14 +3643,10 @@ void TriangulationLevel<1>::monitor_memory (const unsigned int true_dimension) c
 #if deal_II_dimension >= 2
 
 void TriangulationLevel<2>::reserve_space (const unsigned int new_quads) {
-  vector<bool>::iterator u = quads.used.begin(),
-                        e = quads.used.end();
-
-  unsigned int used_quads = 0;
-  for (; u!=e; ++u)
-    ++used_quads;
-
-  const unsigned int new_size = used_quads + new_quads;
+  const unsigned int new_size = new_quads +
+                               count_if (quads.used.begin(),
+                                         quads.used.end(),
+                                         bind2nd (equal_to<bool>(), true));
 
                                   // see above...
   if (new_size>quads.quads.size())

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.