Make maxIterations configurable

dune/fufem/geometry/polyhedrondistance.hh

@@ -95,7 +95,8 @@ void approximate(Coordinate const &target, Dune::DynamicMatrix<double> const &A,
 template <class Coordinate>
 double distance(const Coordinate &target,
                 const ConvexPolyhedron<Coordinate> &segment,
-                double correctionTolerance) {
+                double correctionTolerance,
+                size_t maxIterations = 1000) {
   size_t nCorners = segment.vertices.size();
   if (nCorners == 0)
     return std::numeric_limits<double>::infinity();
@@ -116,9 +117,7 @@ double distance(const Coordinate &target,
   Coordinate oldCoordinates = segment.barycentricToGlobal(sol);
   Coordinate newCoordinates;
 
-  size_t maxIterations = 1000;
-  size_t iterations;
-  for (iterations = 0; iterations < maxIterations; ++iterations) {
+  for (size_t iterations = 0; iterations < maxIterations; ++iterations) {
     approximate(target, A, l, segment, sol);
 
     newCoordinates = segment.barycentricToGlobal(sol);
@@ -127,8 +126,6 @@ double distance(const Coordinate &target,
 
     oldCoordinates = newCoordinates;
   }
-  assert(iterations < maxIterations);
-
   segment.sanityCheck(sol, barycentricTolerance);
   return (target - newCoordinates).two_norm();
 }
@@ -137,7 +134,8 @@ double distance(const Coordinate &target,
 template <class Coordinate>
 double distance(const ConvexPolyhedron<Coordinate> &s1,
                 const ConvexPolyhedron<Coordinate> &s2,
-                double valueCorrectionTolerance) {
+                double valueCorrectionTolerance,
+                size_t maxIterations = 1000) {
   size_t nCorners1 = s1.vertices.size();
   size_t nCorners2 = s2.vertices.size();
   if (nCorners1 == 0 or nCorners2 == 0)
@@ -163,9 +161,7 @@ double distance(const ConvexPolyhedron<Coordinate> &s1,
   Dune::DynamicMatrix<double> A2(nCorners2, nCorners2);
   populateMatrix(s2, A2);
 
-  size_t maxIterations = 1000;
-  size_t iterations;
-  for (iterations = 0; iterations < maxIterations; ++iterations) {
+  for (size_t iterations = 0; iterations < maxIterations; ++iterations) {
     Dune::DynamicVector<double> l2(nCorners2);
     for (size_t i = 0; i < nCorners2; i++)
       l2[i] = c1 * s2.vertices[i];
@@ -180,7 +176,6 @@ double distance(const ConvexPolyhedron<Coordinate> &s1,
     c1 = s1.barycentricToGlobal(sol1);
 
     double const newDistance = (c2 - c1).two_norm();
-    assert(newDistance - oldDistance <= 1e-12); // debugging
     if (oldDistance - newDistance <= valueCorrectionTolerance) {
       s1.sanityCheck(sol1, barycentricTolerance);
       s2.sanityCheck(sol2, barycentricTolerance);
@@ -189,7 +184,7 @@ double distance(const ConvexPolyhedron<Coordinate> &s1,
 
     oldDistance = newDistance;
   }
-  assert(false);
+  return oldDistance;
 }
 
 // Point-to-Geometry convenience method
@@ -209,7 +204,8 @@ double distance(typename Geometry::GlobalCoordinate const &point,
 // Polyhedron-to-Geometry convenience method
 template <class Geometry>
 double distance(ConvexPolyhedron<typename Geometry::GlobalCoordinate> const &p1,
-                Geometry const &geo, double valueCorrectionTolerance) {
+                Geometry const &geo, double valueCorrectionTolerance,
+                size_t maxIterations = 1000) {
   using Coordinate = typename Geometry::GlobalCoordinate;
 
   ConvexPolyhedron<Coordinate> p2;
@@ -217,7 +213,7 @@ double distance(ConvexPolyhedron<typename Geometry::GlobalCoordinate> const &p1,
   for (size_t i = 0; i < p2.vertices.size(); ++i)
     p2.vertices[i] = geo.corner(i);
 
-  return distance(p2, p1, valueCorrectionTolerance);
+  return distance(p2, p1, valueCorrectionTolerance, maxIterations);
 }
 
 #endif