[Problem] Use the original coordinate system

The original coordinates are used and then rotated

This also reveals that the left angle is closer to 15.1 degrees than to
15; new calculations thus yield slightly different results

src/sand-wedge-data/geometry.tex

@@ -4,58 +4,50 @@
 \usetikzlibrary{calc}
 \usetikzlibrary{decorations.pathreplacing}
 
-
 \usepackage{siunitx}
 
 \begin{document}
-\begin{tikzpicture}[scale=12]
-  \pgfmathsetmacro{\rightleg}{0.27};
-  \pgfmathsetmacro{\leftleg}{1.00};
-  \pgfmathsetmacro{\hypot}{sqrt(\leftleg^2+\rightleg^2)};
-  \pgfmathsetmacro{\myangle}{atan(\rightleg/\leftleg)};
-  \pgfmathsetmacro{\ssin}{sin(\myangle)};
-  \pgfmathsetmacro{\scos}{cos(\myangle)};
-  \pgfmathsetmacro{\lsin}{sin(90-\myangle)};
+\pgfmathsetmacro{\rightleg}{0.27}
+\pgfmathsetmacro{\leftleg}{1.00}
+\pgfmathsetmacro{\leftangle}{atan(\rightleg/\leftleg)}
+\begin{tikzpicture}[scale=12, rotate=\leftangle]
+  \pgfmathsetmacro{\mysin}{sin(\leftangle)}
+  \pgfmathsetmacro{\mycos}{cos(\leftangle)}
+  \pgfmathsetmacro{\viscoheight}{0.06}
+  \pgfmathsetmacro{\Zx}{0.35}
+  \pgfmathsetmacro{\weaklen}{0.20}
 
   \coordinate (A) at (0,0);
   \node at (A) [left] {A};
-  \coordinate (B) at (\hypot,0);
+  \coordinate (B) at (\leftleg,-\rightleg);
   \node at (B) [right] {B};
-  \coordinate (C) at (\lsin, \rightleg*\lsin);
-  \node at (C) [above] {C};
+  \coordinate (C) at (\leftleg,0);
+  \node at (C) [right] {C};
 
-  \draw (A)
-  -- (B)
-  -- (C)
-  -- node [above=.5cm, sloped] {$\overline{AC}=\SI{100}{cm}$} (A);
+  \draw (A) -- (B) -- (C) -- node [above=.5cm, sloped] {$\overline{AC}=\SI{100}{cm}$} (A);
 
-  \coordinate (Y) at ($(0.35/\lsin, 0)$);
-  \node at (Y) [below right] {Y};
-  \coordinate (Z) at ($(0.35*\lsin,0.35*\ssin)$);
+  \coordinate (Z) at (\Zx,0);
   \node at (Z) [above] {Z};
-  \coordinate (X) at ($(Y) - (0.20,0)$);
-  \node at (X) [below left] {X};
-  \path let \p1 = (X), \p2 = (Y), \p3 = (Z) in
-    coordinate (U) at ($(\x1 * \x3 / \x2, \x1 * \y3 / \x2)$);
+  \coordinate (Y) at ($(\Zx,-\Zx/\leftleg * \rightleg)$);
+  \node at (Y) [below] {Y};
+  \coordinate (X) at ($(Y) + (-\weaklen*\mycos,\weaklen*\mysin)$);
+  \node at (X) [below] {X};
+  \path let \p1 = (X) in coordinate (U) at ($(\x1, 0)$);
   \node at (U) [above] {U};
 
   \path (A) -- node [above=.25cm, sloped] {$\overline{AZ} = \SI{35}{cm}$} (Z);
 
-  \draw[color=red, thick] (X)
-  -- node [below=.25cm, sloped] {$\overline{XY}=\SI{20}{cm}$} (Y);
+  \draw[color=red, thick] (X) -- node [below=.25cm] {$\overline{XY}=\SI{20}{cm}$} (Y);
   \draw[dashed] (Y) -- (Z);
   \draw[dashed] (U) -- (X);
 
-  \coordinate (K) at ($(B) - (0.06 / \ssin,0)$);
+  \coordinate (K) at ($(B) + (-\leftleg * \viscoheight / \rightleg,\viscoheight)$);
   \node at (K) [below] {K};
-  \coordinate (M) at ($(B) + (-0.06 * \ssin, 0.06 * \lsin)$);
+  \coordinate (M) at ($(B) + (0, \viscoheight)$);
   \node at (M) [right] {M};
-
   \path (C) -- node [right=.5cm] {$\overline{CM} = \SI{21}{cm}$} (M);
 
-  \draw[thick] (K)
-  -- (B)
-  -- node [right=.75cm] {$\overline{BM}=\SI{6}{cm}$} (M) -- cycle;
+  \path[fill=blue] (K) -- (B) -- node [right=.75cm] {$\overline{BM}=\SI{6}{cm}$} (M) -- cycle;
 
   \coordinate (G) at ($(A) ! 0.5 ! (X)$);
   \node at (G) [below] {G};
@@ -67,11 +59,7 @@
   \coordinate (I) at ($(Y) + (G)$);
   \node at (I) [below] {I};
 
-  % Note: Lengths and position are arbitrarily chosen here
-  \draw[->] (0.4,0.45) -- node [above, sloped] {vertical} +(-0.15*\ssin,0.15*\scos);
-  \draw[->] (0.4,0.45) -- node [above, sloped] {horizontal} +(-0.15*\scos,-0.15*\ssin);
-
-  \node[align=left] at (0.5,-0.125) {
+  \node[align=left] at (0.5,-0.225) {
     $Z$: coast line\\
     $\overline{XY}$: velocity weakening zone\\
     $BKM$: visco-elastic domain};

src/sand-wedge-data/mygeometry.cc

@@ -53,14 +53,13 @@ void MyGeometry::render() {
 
   cr->scale(widthScale, heightScale);
   cr->translate(0.1, 0.1);
-  cr->rotate(leftAngle);
   cr->set_line_width(0.0025);
 
   // triangle
   {
-    moveTo(A);
-    lineTo(B);
-    lineTo(C);
+    moveTo(reference::A);
+    lineTo(reference::B);
+    lineTo(reference::C);
     cr->close_path();
     cr->stroke();
   }
@@ -70,10 +69,10 @@ void MyGeometry::render() {
     cr->save();
     std::vector<double> dashPattern = { 0.005 };
     cr->set_dash(dashPattern, 0);
-    moveTo(Z);
-    lineTo(Y);
-    moveTo(U);
-    lineTo(X);
+    moveTo(reference::Z);
+    lineTo(reference::Y);
+    moveTo(reference::U);
+    lineTo(reference::X);
     cr->stroke();
     cr->restore();
   }
@@ -82,9 +81,9 @@ void MyGeometry::render() {
   {
     cr->save();
     setRGBColor(0x0097E0);
-    moveTo(B);
-    lineTo(K);
-    lineTo(M);
+    moveTo(reference::B);
+    lineTo(reference::K);
+    lineTo(reference::M);
     cr->fill();
     cr->restore();
   }
@@ -94,8 +93,8 @@ void MyGeometry::render() {
     cr->save();
     setRGBColor(0x7AD3FF);
     cr->set_line_width(0.005);
-    moveTo(X);
-    lineTo(Y);
+    moveTo(reference::X);
+    lineTo(reference::Y);
     cr->stroke();
     cr->restore();
   }
@@ -109,19 +108,19 @@ void MyGeometry::render() {
 
     cr->save();
     setRGBColor(0x002F47);
-    drawCircle(A);
-    drawCircle(B);
-    drawCircle(C);
-    drawCircle(Y);
-    drawCircle(X);
-    drawCircle(Z);
-    drawCircle(U);
-    drawCircle(K);
-    drawCircle(M);
-    drawCircle(G);
-    drawCircle(H);
-    drawCircle(J);
-    drawCircle(I);
+    drawCircle(reference::A);
+    drawCircle(reference::B);
+    drawCircle(reference::C);
+    drawCircle(reference::Y);
+    drawCircle(reference::X);
+    drawCircle(reference::Z);
+    drawCircle(reference::U);
+    drawCircle(reference::K);
+    drawCircle(reference::M);
+    drawCircle(reference::G);
+    drawCircle(reference::H);
+    drawCircle(reference::J);
+    drawCircle(reference::I);
     cr->restore();
   }
 
@@ -140,19 +139,19 @@ void MyGeometry::render() {
     cr->set_font_face(font);
     cr->set_font_size(0.03);
 
-    label(A, "A");
-    label(B, "B");
-    label(C, "C");
-    label(K, "K");
-    label(M, "M");
-    label(U, "U");
-    label(X, "X");
-    label(Y, "Y");
-    label(Z, "Z");
-    label(G, "G");
-    label(H, "H");
-    label(J, "J");
-    label(I, "I");
+    label(reference::A, "A");
+    label(reference::B, "B");
+    label(reference::C, "C");
+    label(reference::K, "K");
+    label(reference::M, "M");
+    label(reference::U, "U");
+    label(reference::X, "X");
+    label(reference::Y, "Y");
+    label(reference::Z, "Z");
+    label(reference::G, "G");
+    label(reference::H, "H");
+    label(reference::J, "J");
+    label(reference::I, "I");
     cr->restore();
   }
 

src/sand-wedge-data/mygeometry.hh

-#ifndef SRC_SAND_WEDGE_DATA_MYGEOMETRY_HH
-#define SRC_SAND_WEDGE_DATA_MYGEOMETRY_HH
+#ifndef SRC_MYGEOMETRY_HH
+#define SRC_MYGEOMETRY_HH
 
-#include <dune/common/fassign.hh>
 #include <dune/common/fvector.hh>
 
 #include <dune/fufem/arithmetic.hh>
@@ -15,62 +14,93 @@
 namespace MyGeometry {
 namespace {
   using LocalVector = Dune::FieldVector<double, 2>;
+  using LocalMatrix = Dune::FieldMatrix<double, 2, 2>;
 
-  // kludge because fieldvectors have no initialiser_list constructor, see
-  // https://dune-project.org/flyspray/index.php?do=details&task_id=1166
-  LocalVector generateVector(double x, double y) {
-    LocalVector tmp;
-    tmp <<= x, y;
-    return tmp;
+  LocalVector createVector(double x, double y) {
+    LocalVector ret{ 0 };
+    ret[0] = x;
+    ret[1] = y;
+    return ret;
   }
 
-  double degreesToRadians(double degrees) {
-    return M_PI * degrees / 180;
-  };
+  LocalMatrix createMatrix(double a11, double a12, double a21, double a22) {
+    LocalMatrix ret{ 0 };
+    ret[0][0] = a11;
+    ret[0][1] = a12;
+    ret[1][0] = a21;
+    ret[1][1] = a22;
+    return ret;
+  }
 
   LocalVector convexCombination(LocalVector const &x, LocalVector const &y) {
-    LocalVector ret;
-    ret = 0.0;
+    LocalVector ret{ 0 };
     Arithmetic::addProduct(ret, 0.5, x);
     Arithmetic::addProduct(ret, 0.5, y);
     return ret;
   }
+}
 
-  double leftLeg = 1.00;
-  double rightLeg = 0.27;
-  double leftAngle = degreesToRadians(15);
-  double rightAngle = degreesToRadians(75);
+namespace reference {
+  double const rightLeg = 0.27;
+  double const leftLeg = 1.00;
+  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 depth = 0.10;
-}
+  LocalVector const A = createVector(0, 0);
+  LocalVector const B = createVector(leftLeg, -rightLeg);
+  LocalVector const C = createVector(leftLeg, 0);
 
-LocalVector const A = generateVector(0, 0);
-LocalVector const B = generateVector(std::hypot(leftLeg, rightLeg), 0);
-LocalVector const C = generateVector(leftLeg * std::sin(rightAngle),
-                                     leftLeg * std::sin(leftAngle));
+  LocalVector const Z = createVector(0.35, 0);
+  LocalVector const Y = createVector(0.35, -0.35 / leftLeg * rightLeg);
+  LocalVector const X = createVector(Y[0] - weakLen * std::cos(leftAngle),
+                                     Y[1] + weakLen * std::sin(leftAngle));
 
-LocalVector const Y = generateVector(0.35 / std::sin(rightAngle), 0);
-LocalVector const X = generateVector(Y[0] - 0.20, 0);
-LocalVector const Z =
-    generateVector(0.35 * std::sin(rightAngle), 0.35 * std::sin(leftAngle));
+  LocalVector const U = createVector(X[0], 0);
 
-LocalVector const U = generateVector(Z[0] * X[0] / Y[0], Z[1] * X[0] / Y[0]);
+  LocalVector const K =
+      createVector(B[0] - leftLeg * viscoHeight / rightLeg, B[1] + viscoHeight);
+  LocalVector const M = createVector(B[0], B[1] + viscoHeight);
 
-LocalVector const K = generateVector(B[0] - 0.06 / std::sin(leftAngle), 0);
-LocalVector const M = generateVector(B[0] - 0.06 * std::sin(leftAngle),
-                                     0.06 * std::sin(rightAngle));
+  LocalVector const G = convexCombination(A, X);
+  LocalVector const H = convexCombination(X, Y);
+  LocalVector const J = convexCombination(Y, B);
 
-LocalVector const G = convexCombination(A, X);
-LocalVector const H = convexCombination(X, Y);
-LocalVector const J = convexCombination(Y, B);
+  LocalVector const I = createVector(Y[0] + G[0], Y[1] + G[1]);
 
-LocalVector const I = generateVector(Y[0] + G[0], 0);
+  LocalVector const zenith = createVector(0, 1);
+  LocalVector const rotatedZenith = createVector(-1, 0);
 
-LocalVector const zenith =
-    generateVector(-std::sin(leftAngle), std::cos(leftAngle));
+  LocalMatrix const rotation =
+      createMatrix(std::cos(leftAngle), -std::sin(leftAngle),
+                   std::sin(leftAngle), std::cos(leftAngle));
+}
+
+namespace {
+  LocalVector rotate(LocalVector const &x) {
+    LocalVector ret{ 0 };
+    reference::rotation.mv(x, ret);
+    return ret;
+  }
+}
 
+LocalVector const A = rotate(reference::A);
+LocalVector const B = rotate(reference::B);
+LocalVector const C = rotate(reference::C);
+LocalVector const G = rotate(reference::G);
+LocalVector const H = rotate(reference::H);
+LocalVector const I = rotate(reference::I);
+LocalVector const J = rotate(reference::J);
+LocalVector const K = rotate(reference::K);
+LocalVector const M = rotate(reference::M);
+LocalVector const U = rotate(reference::U);
+LocalVector const X = rotate(reference::X);
+LocalVector const Y = rotate(reference::Y);
+LocalVector const Z = rotate(reference::Z);
+
+LocalVector const zenith = rotate(reference::zenith);
 LocalVector const rotatedZenith =
-    generateVector(-std::cos(leftAngle), -std::sin(leftAngle));
+    rotate(reference::rotatedZenith); // FIXME: misnomer
 
 double verticalProjection(LocalVector const &);
 double horizontalProjection(LocalVector const &);