]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Split out main algorithm into a function of its own.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 15 Jul 2009 14:48:17 +0000 (14:48 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 15 Jul 2009 14:48:17 +0000 (14:48 +0000)
git-svn-id: https://svn.dealii.org/trunk@19087 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/grid/grid_tools.cc

index 427ff0d35407ef40eb62a08c818ec080add24afa..3eb1644f959528ae891fa8c0dcf80be9e352def9 100644 (file)
@@ -1325,6 +1325,153 @@ namespace internal
        return objective;
       }
 
+
+
+                                      /**
+                                       * Try to fix up a single cell. Return
+                                       * whether we succeeded with this.
+                                       */
+      template <int dim, int spacedim>
+      bool
+      fix_up_cell (const typename dealii::Triangulation<dim,spacedim>::cell_iterator &cell)
+      {
+                                        // right now we can only deal
+                                        // with cells that have been
+                                        // refined isotropically
+                                        // because that is the only
+                                        // case where we have a cell
+                                        // mid-point that can be moved
+                                        // around without having to
+                                        // consider boundary
+                                        // information
+       Assert (cell->has_children(), ExcInternalError());      
+       Assert (cell->refinement_case() == RefinementCase<dim>::isotropic_refinement,
+               ExcNotImplemented());
+
+                                        // get the current location of
+                                        // the cell mid-vertex:
+       Point<spacedim> cell_mid_point
+         = cell->child(0)->vertex (GeometryInfo<dim>::max_children_per_cell-1);
+
+                                        // now do a few steepest
+                                        // descent steps to reduce the
+                                        // objective function
+       unsigned int iteration = 0;
+       do
+         {
+                                            // choose a step length
+                                            // that is initially 1/10
+                                            // of the child cells'
+                                            // diameter, and a sequence
+                                            // whose sum does not
+                                            // converge (to avoid
+                                            // premature termination of
+                                            // the iteration)
+           const double step_length = cell->diameter() / 10 / (iteration + 1);
+      
+                                            // compute the objective
+                                            // function and its derivative
+           const double val = objective_function<dim,spacedim> (cell, cell_mid_point);
+
+           Tensor<1,dim> gradient;
+           for (unsigned int d=0; d<dim; ++d)
+             {
+               Point<dim> h;
+               h[d] = step_length/2;
+               gradient[d] = (objective_function<dim,spacedim> (cell,
+                                                                cell_mid_point + h)
+                              -
+                              objective_function<dim,spacedim> (cell,
+                                                                cell_mid_point - h))
+                             /
+                             step_length;
+             }
+
+                                            // so we need to go in
+                                            // direction -gradient. the
+                                            // optimal value of the
+                                            // objective function is
+                                            // zero, so assuming that
+                                            // the model is quadratic
+                                            // we would have to go
+                                            // -2*val/||gradient|| in
+                                            // this direction, make
+                                            // sure we go at most
+                                            // step_length into this
+                                            // direction
+           cell_mid_point -= std::min(2*val / (gradient*gradient),
+                                      step_length / gradient.norm()) *
+                             gradient;
+
+           ++iteration;
+         }
+       while (iteration < 10);
+
+
+                                        // verify that the new location
+                                        // is indeed better than the
+                                        // one before in terms of the
+                                        // minimal jacobian determinant
+       double old_min_jacobian, new_min_jacobian;
+
+       for (unsigned int test=0; test<2; ++test)
+         {
+           Point<spacedim> child_vertices
+             [GeometryInfo<dim>::max_children_per_cell]
+             [GeometryInfo<dim>::vertices_per_cell];
+           double child_determinants
+             [GeometryInfo<dim>::max_children_per_cell]
+             [GeometryInfo<dim>::vertices_per_cell];
+
+           if (test == 1)
+             for (unsigned int c=0; c<cell->n_children(); ++c)
+               for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+                 child_vertices[c][i] = cell->child(c)->vertex(i);
+
+                                            // replace mid-cell
+                                            // vertex. note that for
+                                            // child i, the mid-cell
+                                            // vertex happens to have the
+                                            // number
+                                            // max_children_per_cell-i
+           for (unsigned int c=0; c<cell->n_children(); ++c)
+             child_vertices[c][GeometryInfo<dim>::max_children_per_cell-c-1]
+               = cell_mid_point;
+       
+           for (unsigned int c=0; c<cell->n_children(); ++c)
+             GeometryInfo<dim>::jacobian_determinants_at_vertices (child_vertices[c],
+                                                                   child_determinants[c]);
+       
+           double min = child_determinants[0][0];
+           for (unsigned int c=0; c<cell->n_children(); ++c)
+             for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+               min = std::min (min,
+                               child_determinants[c][i]);
+
+           if (test == 0)
+             old_min_jacobian = min;
+           else
+             new_min_jacobian = min;
+         }
+
+                                        // if new minimum jacobian is
+                                        // better than before, and if
+                                        // it is positive, then set the
+                                        // new mid point. otherwise
+                                        // return this cell as one of
+                                        // those that can't apparently
+                                        // be fixed
+       if ((new_min_jacobian > old_min_jacobian)
+           &&
+           (new_min_jacobian > 0))
+         {
+           cell->child(0)->vertex (GeometryInfo<dim>::max_children_per_cell-1)
+             = cell_mid_point;
+           return true;
+         }
+       else
+         return false;
+      }
     }
   }
 }
@@ -1344,147 +1491,8 @@ fix_up_distorted_child_cells (const typename Triangulation<dim,spacedim>::Distor
   for (typename std::list<typename Triangulation<dim,spacedim>::cell_iterator>::const_iterator
         cell_ptr = distorted_cells.distorted_cells.begin();
        cell_ptr != distorted_cells.distorted_cells.end(); ++cell_ptr)
-    {
-      const typename Triangulation<dim,spacedim>::cell_iterator
-       cell = *cell_ptr;
-
-                                      // right now we can only deal
-                                      // with cells that have been
-                                      // refined isotropically
-                                      // because that is the only
-                                      // case where we have a cell
-                                      // mid-point that can be moved
-                                      // around without having to
-                                      // consider boundary
-                                      // information
-      Assert (cell->has_children(), ExcInternalError());      
-      Assert (cell->refinement_case() == RefinementCase<dim>::isotropic_refinement,
-             ExcNotImplemented());
-
-                                      // get the current location of
-                                      // the cell mid-vertex:
-      Point<spacedim> cell_mid_point
-       = cell->child(0)->vertex (GeometryInfo<dim>::max_children_per_cell-1);
-
-                                      // now do a few steepest
-                                      // descent steps to reduce the
-                                      // objective function
-      unsigned int iteration = 0;
-      do
-       {
-                                          // choose a step length
-                                          // that is initially 1/10
-                                          // of the child cells'
-                                          // diameter, and a sequence
-                                          // whose sum does not
-                                          // converge (to avoid
-                                          // premature termination of
-                                          // the iteration)
-         const double step_length = cell->diameter() / 10 / (iteration + 1);
-      
-                                      // compute the objective
-                                      // function and its derivative
-         const double val = internal::GridTools::FixUpDistortedChildCells::
-                            objective_function<dim,spacedim> (cell, cell_mid_point);
-
-         Tensor<1,dim> gradient;
-         for (unsigned int d=0; d<dim; ++d)
-           {
-             Point<dim> h;
-             h[d] = step_length/2;
-             gradient[d] = (internal::GridTools::FixUpDistortedChildCells::
-                            objective_function<dim,spacedim> (cell,
-                                                              cell_mid_point + h)
-                            -
-                            internal::GridTools::FixUpDistortedChildCells::
-                            objective_function<dim,spacedim> (cell,
-                                                              cell_mid_point - h))
-                           /
-                           step_length;
-           }
-
-                                          // so we need to go in
-                                          // direction -gradient. the
-                                          // optimal value of the
-                                          // objective function is
-                                          // zero, so assuming that
-                                          // the model is quadratic
-                                          // we would have to go
-                                          // -2*val/||gradient|| in
-                                          // this direction, make
-                                          // sure we go at most
-                                          // step_length into this
-                                          // direction
-         cell_mid_point -= std::min(2*val / (gradient*gradient),
-                                    step_length / gradient.norm()) *
-                           gradient;
-
-         ++iteration;
-       }
-      while (iteration < 10);
-
-
-                                      // verify that the new location
-                                      // is indeed better than the
-                                      // one before in terms of the
-                                      // minimal jacobian determinant
-      double old_min_jacobian, new_min_jacobian;
-
-      for (unsigned int test=0; test<2; ++test)
-       {
-         Point<spacedim> child_vertices
-           [GeometryInfo<dim>::max_children_per_cell]
-           [GeometryInfo<dim>::vertices_per_cell];
-         double child_determinants
-           [GeometryInfo<dim>::max_children_per_cell]
-           [GeometryInfo<dim>::vertices_per_cell];
-
-         if (test == 1)
-           for (unsigned int c=0; c<cell->n_children(); ++c)
-             for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
-               child_vertices[c][i] = cell->child(c)->vertex(i);
-
-                                          // replace mid-cell
-                                          // vertex. note that for
-                                          // child i, the mid-cell
-                                          // vertex happens to have the
-                                          // number
-                                          // max_children_per_cell-i
-         for (unsigned int c=0; c<cell->n_children(); ++c)
-           child_vertices[c][GeometryInfo<dim>::max_children_per_cell-c-1]
-             = cell_mid_point;
-       
-         for (unsigned int c=0; c<cell->n_children(); ++c)
-           GeometryInfo<dim>::jacobian_determinants_at_vertices (child_vertices[c],
-                                                                 child_determinants[c]);
-       
-         double min = child_determinants[0][0];
-         for (unsigned int c=0; c<cell->n_children(); ++c)
-           for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
-             min = std::min (min,
-                             child_determinants[c][i]);
-
-         if (test == 0)
-           old_min_jacobian = min;
-         else
-           new_min_jacobian = min;
-       }
-
-                                      // if new minimum jacobian is
-                                      // better than before, and if
-                                      // it is positive, then set the
-                                      // new mid point. otherwise
-                                      // return this cell as one of
-                                      // those that can't apparently
-                                      // be fixed
-      if ((new_min_jacobian > old_min_jacobian)
-         &&
-         (new_min_jacobian > 0))
-       cell->child(0)->vertex (GeometryInfo<dim>::max_children_per_cell-1)
-         = cell_mid_point;
-      else
-       unfixable_subset.distorted_cells.push_back (cell);
-    }
+    if (! internal::GridTools::FixUpDistortedChildCells::fix_up_cell<dim,spacedim> (*cell_ptr))
+      unfixable_subset.distorted_cells.push_back (*cell_ptr);
   
   return unfixable_subset;
 }

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.