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;
+ }
}
}
}
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;
}