diff --git a/src/sand-wedge-data/mygeometry.hh b/src/sand-wedge-data/mygeometry.hh index 98bdc3d65e562313d1de07efda621f5ae4a5bfca..19e0553330e13e560cc6e906d75fffc56e0cf92a 100644 --- a/src/sand-wedge-data/mygeometry.hh +++ b/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); diff --git a/src/sand-wedge-data/mygrid.cc b/src/sand-wedge-data/mygrid.cc index d19e91d60756515b6f8c321fc62ee530d21afae3..604d2db66d02690abe24cee5cb84a9bd2745979a 100644 --- a/src/sand-wedge-data/mygrid.cc +++ b/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; diff --git a/src/sand-wedge-data/mygrid.hh b/src/sand-wedge-data/mygrid.hh index efdc989d3f238b93c04e885533496fb97579d494..c4d0768e18f7e1594eaaf38f42df48b8611eea87 100644 --- a/src/sand-wedge-data/mygrid.hh +++ b/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> diff --git a/src/sand-wedge-data/patchfunction.hh b/src/sand-wedge-data/patchfunction.hh index 197ae5cfbdfc948690d99d98ebfa9d57f44949df..864ee8c8bedae05bdd4854c6dfcea9fc386c1bb3 100644 --- a/src/sand-wedge-data/patchfunction.hh +++ b/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 { diff --git a/src/sand-wedge.cc b/src/sand-wedge.cc index 799b6509dd4fbfae9e9eed4f8c91b7ba918e3d26..f53ecbc8d683fb3586d6633b8340778bed5dec4b 100644 --- a/src/sand-wedge.cc +++ b/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[]) {