#ifndef SRC_MULTI_BODY_PROBLEM_DATA_CUBOIDGEOMETRY_HH
#define SRC_MULTI_BODY_PROBLEM_DATA_CUBOIDGEOMETRY_HH

#include <dune/common/fvector.hh>

#include "midpoint.hh"

class CuboidGeometry {
  public:
    typedef Dune::FieldVector<double, MY_DIM> LocalVector;

    constexpr static double const lengthScale = 1.0; // scaling factor

  private:
    const double length_;
    const double width_;
    //const double weakLength_;     // length of the weak zone

  public:
    const LocalVector A;
    const LocalVector B;
    const LocalVector C;
    const LocalVector D;
    const LocalVector X;
    const LocalVector Y;

#if MY_DIM == 3
    const double depth;

    CuboidGeometry(const LocalVector& origin, const LocalVector& weakOrigin, const double weakLength = 0.20, const double length = 1.00, const double width = 0.27, const double depth = 0.60);
#else
        CuboidGeometry(const LocalVector& origin, const LocalVector& weakOrigin, const double weakLength = 0.20, const double length = 1.00, const double width = 0.27);
#endif

    void write() const;

    void render() const;
};
#endif

  /* from Elias
  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 = {0, 0};
  LocalVector2D const B = {leftLeg, -rightLeg};
  LocalVector2D const C = {leftLeg, 0};

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

  LocalVector2D const U = {X[0], 0};

  LocalVector2D const K = {B[0] - leftLeg * viscoHeight / rightLeg,
                           B[1] + viscoHeight};
  LocalVector2D const M = {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 = {Y[0] + G[0], Y[1] + G[1]};

  LocalVector2D const zenith = {0, 1};

  LocalMatrix2D const rotation = {{std::cos(leftAngle), -std::sin(leftAngle)},
                                  {std::sin(leftAngle), std::cos(leftAngle)}};
  */