diff --git a/src/sand-wedge-data/mygeometry.cc b/src/sand-wedge-data/mygeometry.cc
index fcc7af22933b768c336d98aa27896067ea010a21..dedefcd8be5039bd8d6125621b38e90cea80ff13 100644
--- a/src/sand-wedge-data/mygeometry.cc
+++ b/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);
diff --git a/src/sand-wedge-data/mygeometry.hh b/src/sand-wedge-data/mygeometry.hh
index 58c97b56c65cc00b64e2992a35fc3bad199d6503..bf3baac50dd4f7c1743de28680590c69e1b53da3 100644
--- a/src/sand-wedge-data/mygeometry.hh
+++ b/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;
}
}
diff --git a/src/sand-wedge-data/patchfunction.hh b/src/sand-wedge-data/patchfunction.hh
index 871e65fdb1fc653261d6846d28a39e79cea30b70..78e90ccdd84e017445ca3bfde55500c4529bd239 100644
--- a/src/sand-wedge-data/patchfunction.hh
+++ b/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;
}
diff --git a/src/sand-wedge-data/twopiece.hh b/src/sand-wedge-data/twopiece.hh
index f243aa0fb10aaec8e68b935f12b25d6b4103310d..d969c845b1a61f691db28591820a139d26d51b95 100644
--- a/src/sand-wedge-data/twopiece.hh
+++ b/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;
}