From 0145e006fcba03676e1970e2413b3171c04a59b2 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Fri, 16 May 2014 20:39:32 +0200
Subject: [PATCH] [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 | 60 +++++++----------
src/sand-wedge-data/mygeometry.cc | 77 +++++++++++----------
src/sand-wedge-data/mygeometry.hh | 108 +++++++++++++++++++-----------
3 files changed, 131 insertions(+), 114 deletions(-)
diff --git a/src/sand-wedge-data/geometry.tex b/src/sand-wedge-data/geometry.tex
index 56254a9b..32d63b6e 100644
--- a/src/sand-wedge-data/geometry.tex
+++ b/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};
diff --git a/src/sand-wedge-data/mygeometry.cc b/src/sand-wedge-data/mygeometry.cc
index 5629772f..fcc7af22 100644
--- a/src/sand-wedge-data/mygeometry.cc
+++ b/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();
}
diff --git a/src/sand-wedge-data/mygeometry.hh b/src/sand-wedge-data/mygeometry.hh
index 0800e114..58c97b56 100644
--- a/src/sand-wedge-data/mygeometry.hh
+++ b/src/sand-wedge-data/mygeometry.hh
@@ -1,7 +1,6 @@
-#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 &);
--
GitLab