From: Wolfgang Bangerth Date: Wed, 15 Jul 2009 14:48:17 +0000 (+0000) Subject: Split out main algorithm into a function of its own. X-Git-Tag: v8.0.0~7473 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=50132224cc1d79d87359362262edafceb9a83aed;p=dealii.git Split out main algorithm into a function of its own. git-svn-id: https://svn.dealii.org/trunk@19087 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/source/grid/grid_tools.cc b/deal.II/deal.II/source/grid/grid_tools.cc index 427ff0d354..3eb1644f95 100644 --- a/deal.II/deal.II/source/grid/grid_tools.cc +++ b/deal.II/deal.II/source/grid/grid_tools.cc @@ -1325,6 +1325,153 @@ namespace internal return objective; } + + + /** + * Try to fix up a single cell. Return + * whether we succeeded with this. + */ + template + bool + fix_up_cell (const typename dealii::Triangulation::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::isotropic_refinement, + ExcNotImplemented()); + + // get the current location of + // the cell mid-vertex: + Point cell_mid_point + = cell->child(0)->vertex (GeometryInfo::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 (cell, cell_mid_point); + + Tensor<1,dim> gradient; + for (unsigned int d=0; d h; + h[d] = step_length/2; + gradient[d] = (objective_function (cell, + cell_mid_point + h) + - + objective_function (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 child_vertices + [GeometryInfo::max_children_per_cell] + [GeometryInfo::vertices_per_cell]; + double child_determinants + [GeometryInfo::max_children_per_cell] + [GeometryInfo::vertices_per_cell]; + + if (test == 1) + for (unsigned int c=0; cn_children(); ++c) + for (unsigned int i=0; i::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; cn_children(); ++c) + child_vertices[c][GeometryInfo::max_children_per_cell-c-1] + = cell_mid_point; + + for (unsigned int c=0; cn_children(); ++c) + GeometryInfo::jacobian_determinants_at_vertices (child_vertices[c], + child_determinants[c]); + + double min = child_determinants[0][0]; + for (unsigned int c=0; cn_children(); ++c) + for (unsigned int i=0; i::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::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::Distor for (typename std::list::cell_iterator>::const_iterator cell_ptr = distorted_cells.distorted_cells.begin(); cell_ptr != distorted_cells.distorted_cells.end(); ++cell_ptr) - { - const typename Triangulation::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::isotropic_refinement, - ExcNotImplemented()); - - // get the current location of - // the cell mid-vertex: - Point cell_mid_point - = cell->child(0)->vertex (GeometryInfo::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 (cell, cell_mid_point); - - Tensor<1,dim> gradient; - for (unsigned int d=0; d h; - h[d] = step_length/2; - gradient[d] = (internal::GridTools::FixUpDistortedChildCells:: - objective_function (cell, - cell_mid_point + h) - - - internal::GridTools::FixUpDistortedChildCells:: - objective_function (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 child_vertices - [GeometryInfo::max_children_per_cell] - [GeometryInfo::vertices_per_cell]; - double child_determinants - [GeometryInfo::max_children_per_cell] - [GeometryInfo::vertices_per_cell]; - - if (test == 1) - for (unsigned int c=0; cn_children(); ++c) - for (unsigned int i=0; i::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; cn_children(); ++c) - child_vertices[c][GeometryInfo::max_children_per_cell-c-1] - = cell_mid_point; - - for (unsigned int c=0; cn_children(); ++c) - GeometryInfo::jacobian_determinants_at_vertices (child_vertices[c], - child_determinants[c]); - - double min = child_determinants[0][0]; - for (unsigned int c=0; cn_children(); ++c) - for (unsigned int i=0; i::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::max_children_per_cell-1) - = cell_mid_point; - else - unfixable_subset.distorted_cells.push_back (cell); - } + if (! internal::GridTools::FixUpDistortedChildCells::fix_up_cell (*cell_ptr)) + unfixable_subset.distorted_cells.push_back (*cell_ptr); return unfixable_subset; }