[Cleanup] Prepare for extension to 3D

src/sand-wedge-data/mygeometry.cc

@@ -48,8 +48,8 @@ void MyGeometry::render() {
                        ((colour & 0x00FF00) >> 8) / 255.0,
                        ((colour & 0x0000FF) >> 0) / 255.0);
   };
-  auto const moveTo = [&](LocalVector const &v) { cr->move_to(v[0], -v[1]); };
-  auto const lineTo = [&](LocalVector const &v) { cr->line_to(v[0], -v[1]); };
+  auto const moveTo = [&](LocalVector2D const &v) { cr->move_to(v[0], -v[1]); };
+  auto const lineTo = [&](LocalVector2D const &v) { cr->line_to(v[0], -v[1]); };
 
   cr->scale(widthScale, heightScale);
   cr->translate(0.1, 0.1);
@@ -101,7 +101,7 @@ void MyGeometry::render() {
 
   // mark points
   {
-    auto const drawCircle = [&](LocalVector const &v) {
+    auto const drawCircle = [&](LocalVector2D const &v) {
       cr->arc(v[0], -v[1], 0.0075, -M_PI, M_PI); // x,y,radius,angle1,angle2
       cr->fill();
     };
@@ -126,7 +126,7 @@ void MyGeometry::render() {
 
   // labels
   {
-    auto const label = [&](LocalVector const &v,
+    auto const label = [&](LocalVector2D const &v,
                            std::string l) { // TODO: implicit?
       moveTo(v);
       cr->rel_move_to(0.005, -0.02);

src/sand-wedge-data/mygeometry.hh

@@ -13,18 +13,22 @@
 
 namespace MyGeometry {
 namespace {
-  using LocalVector = Dune::FieldVector<double, 2>;
-  using LocalMatrix = Dune::FieldMatrix<double, 2, 2>;
+  using LocalVector2D = Dune::FieldVector<double, 2>;
+  using LocalMatrix2D = Dune::FieldMatrix<double, 2, 2>;
 
-  LocalVector createVector(double x, double y) {
-    LocalVector ret{ 0 };
+  using LocalVector = Dune::FieldVector<double, DIM>;
+
+  // kludge because fieldvectors have no initialiser_list constructor, see
+  // https://dune-project.org/flyspray/index.php?do=details&task_id=1166
+  LocalVector2D createVector(double x, double y) {
+    LocalVector2D ret{ 0 };
     ret[0] = x;
     ret[1] = y;
     return ret;
   }
 
-  LocalMatrix createMatrix(double a11, double a12, double a21, double a22) {
-    LocalMatrix ret{ 0 };
+  LocalMatrix2D createMatrix(double a11, double a12, double a21, double a22) {
+    LocalMatrix2D ret{ 0 };
     ret[0][0] = a11;
     ret[0][1] = a12;
     ret[1][0] = a21;
@@ -32,8 +36,9 @@ namespace {
     return ret;
   }
 
-  LocalVector convexCombination(LocalVector const &x, LocalVector const &y) {
-    LocalVector ret{ 0 };
+  LocalVector2D convexCombination(LocalVector2D const &x,
+                                  LocalVector2D const &y) {
+    LocalVector2D ret{ 0 };
     Arithmetic::addProduct(ret, 0.5, x);
     Arithmetic::addProduct(ret, 0.5, y);
     return ret;
@@ -47,39 +52,42 @@ namespace reference {
   double const viscoHeight = 0.06; // Height of the viscous bottom layer
   double const weakLen = 0.20;     // Length of the weak zone
 
-  LocalVector const A = createVector(0, 0);
-  LocalVector const B = createVector(leftLeg, -rightLeg);
-  LocalVector const C = createVector(leftLeg, 0);
+  LocalVector2D const A = createVector(0, 0);
+  LocalVector2D const B = createVector(leftLeg, -rightLeg);
+  LocalVector2D const C = createVector(leftLeg, 0);
 
-  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));
+  LocalVector2D const Z = createVector(0.35, 0);
+  LocalVector2D const Y = createVector(0.35, -0.35 / leftLeg * rightLeg);
+  LocalVector2D const X = createVector(Y[0] - weakLen * std::cos(leftAngle),
+                                       Y[1] + weakLen * std::sin(leftAngle));
 
-  LocalVector const U = createVector(X[0], 0);
+  LocalVector2D const U = createVector(X[0], 0);
 
-  LocalVector const K =
+  LocalVector2D const K =
       createVector(B[0] - leftLeg * viscoHeight / rightLeg, B[1] + viscoHeight);
-  LocalVector const M = createVector(B[0], B[1] + viscoHeight);
+  LocalVector2D const M = createVector(B[0], B[1] + viscoHeight);
 
-  LocalVector const G = convexCombination(A, X);
-  LocalVector const H = convexCombination(X, Y);
-  LocalVector const J = convexCombination(Y, B);
+  LocalVector2D const G = convexCombination(A, X);
+  LocalVector2D const H = convexCombination(X, Y);
+  LocalVector2D const J = convexCombination(Y, B);
 
-  LocalVector const I = createVector(Y[0] + G[0], Y[1] + G[1]);
+  LocalVector2D const I = createVector(Y[0] + G[0], Y[1] + G[1]);
 
-  LocalVector const zenith = createVector(0, 1);
-  LocalVector const rotatedZenith = createVector(-1, 0);
+  LocalVector2D const zenith = createVector(0, 1);
+  LocalVector2D const rotatedZenith = createVector(-1, 0);
 
-  LocalMatrix const rotation =
+  LocalMatrix2D const rotation =
       createMatrix(std::cos(leftAngle), -std::sin(leftAngle),
                    std::sin(leftAngle), std::cos(leftAngle));
 }
 
 namespace {
-  LocalVector rotate(LocalVector const &x) {
+  LocalVector rotate(LocalVector2D const &x) {
+    LocalVector2D ret2D;
+    reference::rotation.mv(x, ret2D);
     LocalVector ret{ 0 };
-    reference::rotation.mv(x, ret);
+    ret[0] = ret2D[0];
+    ret[1] = ret2D[1];
     return ret;
   }
 }

src/sand-wedge-data/patchfunction.hh

@@ -6,7 +6,7 @@
 #include <dune/common/parametertree.hh>
 
 class PatchFunction
-    : public Dune::VirtualFunction<Dune::FieldVector<double, 2>,
+    : public Dune::VirtualFunction<Dune::FieldVector<double, DIM>,
                                    Dune::FieldVector<double, 1>> {
 private:
   bool static isBetween(double x, double lower, double upper) {
@@ -17,13 +17,13 @@ class PatchFunction
     return std::abs(a - b) < 1e-14;
   };
 
-  bool insideRegion2(Dune::FieldVector<double, 2> const &z) const {
+  bool insideRegion2(Dune::FieldVector<double, DIM> const &z) const {
     assert(isClose(0, z[1]));
     return isBetween(z[0], _X[0], _Y[0]);
   }
 
-  Dune::FieldVector<double, 2> const &_X;
-  Dune::FieldVector<double, 2> const &_Y;
+  Dune::FieldVector<double, DIM> const &_X;
+  Dune::FieldVector<double, DIM> const &_Y;
 
   double const _v1;
   double const _v2;
@@ -32,7 +32,7 @@ class PatchFunction
   PatchFunction(double v1, double v2)
       : _X(MyGeometry::X), _Y(MyGeometry::Y), _v1(v1), _v2(v2) {}
 
-  void evaluate(Dune::FieldVector<double, 2> const &x,
+  void evaluate(Dune::FieldVector<double, DIM> const &x,
                 Dune::FieldVector<double, 1> &y) const {
     y = insideRegion2(x) ? _v2 : _v1;
   }

src/sand-wedge-data/twopiece.hh

@@ -8,20 +8,20 @@
 #include "mygeometry.hh"
 
 class TwoPieceFunction
-    : public Dune::VirtualFunction<Dune::FieldVector<double, 2>,
+    : public Dune::VirtualFunction<Dune::FieldVector<double, DIM>,
                                    Dune::FieldVector<double, 1>> {
 private:
-  bool liesBelow(Dune::FieldVector<double, 2> const &x,
-                 Dune::FieldVector<double, 2> const &y,
-                 Dune::FieldVector<double, 2> const &z) const {
+  bool liesBelow(Dune::FieldVector<double, DIM> const &x,
+                 Dune::FieldVector<double, DIM> const &y,
+                 Dune::FieldVector<double, DIM> const &z) const {
     return (z[0] - x[0]) * (y[1] - x[1]) / (y[0] - x[0]) >= z[1] - x[1];
   };
-  bool insideRegion2(Dune::FieldVector<double, 2> const &z) const {
+  bool insideRegion2(Dune::FieldVector<double, DIM> const &z) const {
     return liesBelow(_K, _M, z);
   };
 
-  Dune::FieldVector<double, 2> const &_K;
-  Dune::FieldVector<double, 2> const &_M;
+  Dune::FieldVector<double, DIM> const &_K;
+  Dune::FieldVector<double, DIM> const &_M;
 
   double const _v1;
   double const _v2;
@@ -30,7 +30,7 @@ class TwoPieceFunction
   TwoPieceFunction(double v1, double v2)
       : _K(MyGeometry::K), _M(MyGeometry::M), _v1(v1), _v2(v2) {}
 
-  void evaluate(Dune::FieldVector<double, 2> const &x,
+  void evaluate(Dune::FieldVector<double, DIM> const &x,
                 Dune::FieldVector<double, 1> &y) const {
     y = insideRegion2(x) ? _v2 : _v1;
   }