[Cleanup] Add and use intrinsic length scale

src/sand-wedge-data/mygeometry.hh

@@ -32,11 +32,13 @@ namespace {
 }
 
 namespace reference {
-  double const rightLeg = 0.27;
-  double const leftLeg = 1.00;
+  double const s = 1.0; // scaling factor
+
+  double const rightLeg = 0.27 * s;
+  double const leftLeg = 1.00 * s;
   double const leftAngle = atan(rightLeg / leftLeg);
-  double const viscoHeight = 0.06; // Height of the viscous bottom layer
-  double const weakLen = 0.20;     // Length of the weak zone
+  double const viscoHeight = 0.06 * s; // Height of the viscous bottom layer
+  double const weakLen = 0.20 * s;     // Length of the weak zone
 
   double const zDistance = 0.35;
 
@@ -44,9 +46,9 @@ namespace reference {
   LocalVector2D const B = createVector(leftLeg, -rightLeg);
   LocalVector2D const C = createVector(leftLeg, 0);
 
-  LocalVector2D const Z = createVector(zDistance, 0);
+  LocalVector2D const Z = createVector(zDistance * s, 0);
   LocalVector2D const Y =
-      createVector(zDistance, -zDistance / leftLeg * rightLeg);
+      createVector(zDistance * s, -zDistance * s / leftLeg * rightLeg);
   LocalVector2D const X = createVector(Y[0] - weakLen * std::cos(leftAngle),
                                        Y[1] + weakLen * std::sin(leftAngle));
 
@@ -81,6 +83,8 @@ namespace {
   }
 }
 
+double const lengthScale = reference::s;
+
 double const depth = 0.10;
 
 LocalVector const A = rotate(reference::A);

src/sand-wedge-data/mygrid.cc

@@ -120,7 +120,7 @@ template <class GridView>
 template <class Vector>
 bool MyFaces<GridView>::xyCollinear(Vector const &a, Vector const &b,
                                     Vector const &c) {
-  return isClose((b[0] - a[0]) * (c[1] - a[1]), (b[1] - a[1]) * (c[0] - a[0]));
+  return isClose2((b[0] - a[0]) * (c[1] - a[1]), (b[1] - a[1]) * (c[0] - a[0]));
 }
 
 template <class GridView>
@@ -130,9 +130,11 @@ bool MyFaces<GridView>::xyBoxed(Vector const &v1, Vector const &v2,
   auto const minmax0 = std::minmax(v1[0], v2[0]);
   auto const minmax1 = std::minmax(v1[1], v2[1]);
 
-  if (minmax0.first - 1e-14 > x[0] or x[0] > minmax0.second + 1e-14)
+  if (minmax0.first - 1e-14 * MyGeometry::lengthScale > x[0] or x[0] >
+      minmax0.second + 1e-14 * MyGeometry::lengthScale)
     return false;
-  if (minmax1.first - 1e-14 > x[1] or x[1] > minmax1.second + 1e-14)
+  if (minmax1.first - 1e-14 * MyGeometry::lengthScale > x[1] or x[1] >
+      minmax1.second + 1e-14 * MyGeometry::lengthScale)
     return false;
 
   return true;

src/sand-wedge-data/mygrid.hh

@@ -28,7 +28,12 @@ template <class GridView> struct MyFaces {
 
 private:
   bool isClose(double a, double b) {
-    return std::abs(a - b) < 1e-14;
+    return std::abs(a - b) < 1e-14 * MyGeometry::lengthScale;
+  };
+
+  bool isClose2(double a, double b) {
+    return std::abs(a - b) <
+           1e-14 * MyGeometry::lengthScale * MyGeometry::lengthScale;
   };
 
   template <class Vector>

src/sand-wedge-data/patchfunction.hh

@@ -14,7 +14,7 @@ class PatchFunction
   }
 
   bool static isClose(double a, double b) {
-    return std::abs(a - b) < 1e-14;
+    return std::abs(a - b) < MyGeometry::lengthScale * 1e-14;
   };
 
   bool insideRegion2(Dune::FieldVector<double, MY_DIM> const &z) const {

src/sand-wedge.cc

@@ -120,7 +120,7 @@ template <class Geometry> double distanceToWeakeningRegion(Geometry const &g) {
 }
 
 double computeAdmissibleDiameter(double distance, double smallestDiameter) {
-  return (distance / 0.0125 + 1.0) * smallestDiameter;
+  return (distance / 0.0125 / MyGeometry::lengthScale + 1.0) * smallestDiameter;
 }
 
 int main(int argc, char *argv[]) {