From 89cbcd997c6ad02c86657a432a8f871816eb1e8e Mon Sep 17 00:00:00 2001
From: Luca Heltai <luca.heltai@sissa.it>
Date: Thu, 7 Aug 2014 17:12:39 +0200
Subject: [PATCH] Made FlatManifold tolerance relative.

---
 include/deal.II/grid/manifold.h | 18 ++++++++----------
 source/grid/manifold.cc         | 12 +++++++-----
 2 files changed, 15 insertions(+), 15 deletions(-)

diff --git a/include/deal.II/grid/manifold.h b/include/deal.II/grid/manifold.h
index 6b1b8e3875..096702d639 100644
--- a/include/deal.II/grid/manifold.h
+++ b/include/deal.II/grid/manifold.h
@@ -259,7 +259,7 @@ public:
    * pi, but 2*pi (or zero), since, on a periodic manifold, these two
    * points are at distance 2*eps and not (2*pi-eps). Special cases
    * are taken into account, to ensure that the behavior is always as
-   * expected. The third argument is used as a tolerance when
+   * expected. The third argument is used as a relative tolerance when
    * computing distances.
    *
    * Periodicity will be intended in the following way: the domain is
@@ -270,8 +270,7 @@ public:
    * compute averages is called, an exception will be thrown if one of
    * the points which you are using for the average lies outside the
    * periodicity box. The return points are garanteed to lie in the
-   * perodicity box plus or minus the tolerance you set as the third
-   * argument.
+   * perodicity box plus or minus tolerance*periodicity.norm().
    */
   FlatManifold (const Point<spacedim> periodicity=Point<spacedim>(),
                 const double tolerance=1e-10);
@@ -312,13 +311,6 @@ public:
   virtual
   Point<spacedim> project_to_manifold (const std::vector<Point<spacedim> > &points,
                                        const Point<spacedim> &candidate) const;
-protected:
-  /**
-   * Tolerance. This tolerance is used to compute distances in double
-   * precision. Anything below this tolerance is considered zero.
-   */
-  const double tolerance;
-
 private:
   /**
    * The periodicity of this Manifold. Periodicity affects the way a
@@ -340,6 +332,12 @@ private:
                  << "The component number " << arg1 << " of the point [ " << arg2
                  << " ]  is not in the interval [ " << -arg4
                  << ", " << arg3[arg4] << "), bailing out.");
+
+  /**
+   * Relative tolerance. This tolerance is used to compute distances
+   * in double precision.
+   */
+  const double tolerance;
 };
 
 
diff --git a/source/grid/manifold.cc b/source/grid/manifold.cc
index ccd7e47c64..a2b1118a25 100644
--- a/source/grid/manifold.cc
+++ b/source/grid/manifold.cc
@@ -261,8 +261,8 @@ get_new_point_on_hex (const typename Triangulation<3, 3>::hex_iterator &hex) con
 template <int dim, int spacedim>
 FlatManifold<dim,spacedim>::FlatManifold (const Point<spacedim> periodicity,
                                           const double tolerance) :
-  tolerance(tolerance),
-  periodicity(periodicity)
+  periodicity(periodicity),
+  tolerance(tolerance)
 {}
 
 template <int dim, int spacedim>
@@ -277,6 +277,8 @@ get_new_point (const Quadrature<spacedim> &quad) const
   double sum=0;
   for (unsigned int i=0; i<weights.size(); ++i)
     sum+= weights[i];
+  // Here it is correct to use tolerance as an absolute one, since
+  // this should be relative to unity.
   Assert(std::abs(sum-1.0) < tolerance, ExcMessage("Weights should sum to 1!"));
 #endif
 
@@ -291,9 +293,9 @@ get_new_point (const Quadrature<spacedim> &quad) const
         {
           minP[d] = std::min(minP[d], surrounding_points[i][d]);
           if (periodicity[d] > 0)
-            Assert( (surrounding_points[i][d] < periodicity[d]+tolerance) ||
-                    (surrounding_points[i][d] >= -tolerance),
-                    ExcPeriodicBox(d, surrounding_points[i], periodicity, tolerance));
+            Assert( (surrounding_points[i][d] < periodicity[d]+tolerance*periodicity.norm()) ||
+                    (surrounding_points[i][d] >= -tolerance*periodicity.norm()),
+                    ExcPeriodicBox(d, surrounding_points[i], periodicity, tolerance*periodicity.norm()));
         }
 
   for (unsigned int i=0; i<surrounding_points.size(); ++i)
-- 
2.39.5