mygeometry.hh

#ifndef SRC_MYGEOMETRY_HH
#define SRC_MYGEOMETRY_HH

#include <dune/common/fvector.hh>

#include "midpoint.hh"

namespace MyGeometry {
namespace {
  using LocalVector2D = Dune::FieldVector<double, 2>;
  using LocalMatrix2D = Dune::FieldMatrix<double, 2, 2>;

  using LocalVector = Dune::FieldVector<double, MY_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;
  }

  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;
    ret[1][1] = a22;
    return ret;
  }
}

namespace reference {
  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 * s; // Height of the viscous bottom layer
  double const weakLen = 0.20 * s;     // Length of the weak zone

  double const zDistance = 0.35;

  LocalVector2D const A = createVector(0, 0);
  LocalVector2D const B = createVector(leftLeg, -rightLeg);
  LocalVector2D const C = createVector(leftLeg, 0);

  LocalVector2D const Z = createVector(zDistance * s, 0);
  LocalVector2D const Y =
      createVector(zDistance * s, -zDistance * s / leftLeg * rightLeg);
  LocalVector2D const X = createVector(Y[0] - weakLen * std::cos(leftAngle),
                                       Y[1] + weakLen * std::sin(leftAngle));

  LocalVector2D const U = createVector(X[0], 0);

  LocalVector2D const K =
      createVector(B[0] - leftLeg * viscoHeight / rightLeg, B[1] + viscoHeight);
  LocalVector2D const M = createVector(B[0], B[1] + viscoHeight);

  LocalVector2D const G = midPoint(A, X);
  LocalVector2D const H = midPoint(X, Y);
  LocalVector2D const J = midPoint(Y, B);

  LocalVector2D const I = createVector(Y[0] + G[0], Y[1] + G[1]);

  LocalVector2D const zenith = createVector(0, 1);
  LocalVector2D const orthogonalZenith = createVector(-1, 0);

  LocalMatrix2D const rotation =
      createMatrix(std::cos(leftAngle), -std::sin(leftAngle),
                   std::sin(leftAngle), std::cos(leftAngle));
}

namespace {
  LocalVector rotate(LocalVector2D const &x) {
    LocalVector2D ret2D;
    reference::rotation.mv(x, ret2D);
    LocalVector ret{ 0 };
    ret[0] = ret2D[0];
    ret[1] = ret2D[1];
    return ret;
  }
}

double const lengthScale = reference::s;

double const depth = 0.60 * lengthScale;

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 orthogonalZenith = rotate(reference::orthogonalZenith);

double verticalProjection(LocalVector const &);
double horizontalProjection(LocalVector const &);

void write();

void render();
}
#endif