From a431d28e2344420fa91aba7d71b4e3b34cdd97a9 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 16 Sep 2020 11:41:19 -0600 Subject: [PATCH] Make the tolerance in step-1 relative. --- examples/step-1/step-1.cc | 35 ++++++++++++++++++++++++++++++----- examples/step-2/step-2.cc | 3 ++- 2 files changed, 32 insertions(+), 6 deletions(-) diff --git a/examples/step-1/step-1.cc b/examples/step-1/step-1.cc index 38eadd26f1..76ea18a660 100644 --- a/examples/step-1/step-1.cc +++ b/examples/step-1/step-1.cc @@ -210,14 +210,39 @@ void second_grid() // If this cell is at the inner boundary, then at least one of its // vertices must sit on the inner ring and therefore have a radial // distance from the center of exactly 0.5, up to floating point - // accuracy. Compute this distance, and if we have found a vertex - // with this property flag this cell for later refinement. We can - // then also break the loop over all vertices and move on to the - // next cell. + // accuracy. So we compute this distance, and if we find a vertex + // with this property, we flag this cell for later refinement. We + // can then also break the loop over all vertices and move on to + // the next cell. + // + // Because the distance from the center is computed as a floating + // point number, we have to expect that whatever we compute is + // only accurate to within + // [round-off](https://en.wikipedia.org/wiki/Round-off_error). As + // a consequence, we can never expect to compare the distance + // with the inner radius by equality: A statement such as + // `if (distance_from_center == inner_radius)` will fail + // unless we get exceptionally lucky. Rather, we need to do this + // comparison with a certain tolerance, and the usual way to do + // this is to write it as `if (std::abs(distance_from_center - + // inner_radius) <= tolerance)` + // where `tolerance` is some small number larger + // than round-off. The question is how to choose it: We could just + // pick, say, `1e-10`, but this is only appropriate if the objects + // we compare are of size one. If we had created a mesh with cells + // of size `1e+10`, then `1e-10` would be far lower than round-off + // and, as before, the comparison will only succeed if we get + // exceptionally lucky. Rather, it is almost always useful to make + // the tolerance *relative* to a typical "scale" of the objects + // being compared. Here, the "scale" would be the inner radius, or + // maybe the diameter of cells. We choose the former and set the + // tolerance equal to $10^{-6}$ times the inner radius of the + // annulus. const double distance_from_center = center.distance(cell->vertex(v)); - if (std::fabs(distance_from_center - inner_radius) < 1e-10) + if (std::fabs(distance_from_center - inner_radius) <= + 1e-6 * inner_radius) { cell->set_refine_flag(); break; diff --git a/examples/step-2/step-2.cc b/examples/step-2/step-2.cc index e8a25a9d04..74ea05a6f5 100644 --- a/examples/step-2/step-2.cc +++ b/examples/step-2/step-2.cc @@ -80,7 +80,8 @@ void make_grid(Triangulation<2> &triangulation) const double distance_from_center = center.distance(cell->vertex(v)); - if (std::fabs(distance_from_center - inner_radius) < 1e-10) + if (std::fabs(distance_from_center - inner_radius) <= + 1e-6 * inner_radius) { cell->set_refine_flag(); break; -- 2.39.5