dune-tectonic.pc.in

@@ -10,6 +10,6 @@ Name: @PACKAGE_NAME@
 Version: @VERSION@
 Description: dune-tectonic module
 URL: http://dune-project.org/
-Requires: dune-common dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid
+Requires: dune-common dune-contact dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid
 Libs: -L${libdir}
 Cflags: -I${includedir}

dune.module

@@ -5,4 +5,4 @@
 Module: dune-tectonic
 Version: 2.5-1
 Maintainer: elias.pipping@fu-berlin.de
-Depends: dune-common dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid
+Depends: dune-common dune-contact dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid

dune/tectonic/globalratestatefriction.hh

@@ -32,7 +32,7 @@ class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
       : restrictions(), localToGlobal(), zeroFriction() {
     auto const gridView = frictionalBoundary.gridView();
     Dune::MultipleCodimMultipleGeomTypeMapper<
-        GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView);
+        GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView, Dune::mcmgVertexLayout());
     for (auto it = gridView.template begin<block_size>();
          it != gridView.template end<block_size>(); ++it) {
       auto const i = vertexMapper.index(*it);
@@ -44,8 +44,8 @@ class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
       restrictions.emplace_back(weights[i], weightedNormalStress[i],
                                 frictionInfo(geoToPoint(it->geometry())));
     }
-    assert(restrictions.size() == frictionalBoundary.numVertices());
-    assert(localToGlobal.size() == frictionalBoundary.numVertices());
+    assert(restrictions.size() == (size_t) frictionalBoundary.numVertices());
+    assert(localToGlobal.size() == (size_t) frictionalBoundary.numVertices());
   }
 
   void updateAlpha(ScalarVector const &alpha) override {

dune/tectonic/myblockproblem.hh

@@ -57,7 +57,7 @@ class MyBlockProblem : /* not public */ BlockNonlinearGSProblem<ConvexProblem> {
   MyBlockProblem(Dune::ParameterTree const &parset, ConvexProblem &problem)
       : Base(parset, problem),
         maxEigenvalues_(problem.f.size()),
-        localBisection(0.0, 1.0, 0.0, true, 0.0) {
+        localBisection(0.0, 1.0, 0.0, 0.0) {
     for (size_t i = 0; i < problem.f.size(); ++i) {
       LocalVectorType eigenvalues;
       Dune::FMatrixHelp::eigenValues(problem.A[i][i], eigenvalues);

src/CMakeLists.txt

@@ -19,6 +19,29 @@ set(SW_SOURCE_FILES
   time-stepping/state.cc
   vtk.cc
 )
+
+set(MSW_SOURCE_FILES
+  assemblers.cc
+  enumparser.cc
+  hdf5/frictionalboundary-writer.cc
+  hdf5/iteration-writer.cc
+  hdf5/patchinfo-writer.cc
+  hdf5/restart-io.cc
+  hdf5/surface-writer.cc
+  hdf5/time-writer.cc
+  multi-body-problem-data/cuboidgeometry.cc
+  multi-body-problem-data/mygrids.cc
+  multi-body-problem.cc
+  spatial-solving/fixedpointiterator.cc
+  spatial-solving/solverfactory.cc
+  time-stepping/adaptivetimestepper.cc
+  time-stepping/coupledtimestepper.cc
+  time-stepping/rate.cc
+  time-stepping/rate/rateupdater.cc
+  time-stepping/state.cc
+  vtk.cc
+)
+
 set(UGW_SOURCE_FILES
   assemblers.cc # FIXME
   one-body-problem-data/mygrid.cc
@@ -28,15 +51,20 @@ set(UGW_SOURCE_FILES
 
 foreach(_dim 2 3)
   set(_sw_target one-body-problem-${_dim}D)
+  set(_msw_target multi-body-problem-${_dim}D)
   set(_ugw_target uniform-grid-writer-${_dim}D)
 
   add_executable(${_sw_target} ${SW_SOURCE_FILES})
+  add_executable(${_msw_target} ${MSW_SOURCE_FILES})
   add_executable(${_ugw_target} ${UGW_SOURCE_FILES})
   add_dune_ug_flags(${_sw_target})
+  add_dune_ug_flags(${_msw_target})
   add_dune_ug_flags(${_ugw_target})
 
   add_dune_hdf5_flags(${_sw_target})
+  add_dune_hdf5_flags(${_msw_target})
 
   set_property(TARGET ${_sw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
+  set_property(TARGET ${_msw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
   set_property(TARGET ${_ugw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
 endforeach()

src/assemblers_tmpl.cc

@@ -4,4 +4,4 @@
 
 #include "explicitgrid.hh"
 
-template class MyAssembler<GridView, MY_DIM>;
+template class MyAssembler<LeafGridView, MY_DIM>;

src/explicitgrid.hh

@@ -2,6 +2,13 @@
 #define SRC_EXPLICITGRID_HH
 
 #include "gridselector.hh"
+#include "explicitvectors.hh"
 
-using GridView = Grid::LeafGridView;
+#include <dune/contact/common/deformedcontinuacomplex.hh>
+
+using LeafGridView = Grid::LeafGridView;
+using LevelGridView = Grid::LevelGridView;
+
+using DeformedGridComplex = typename Dune::Contact::DeformedContinuaComplex<Grid, Vector>;
+using DeformedGrid = DeformedGridComplex::DeformedGridType;
 #endif

src/hdf5/frictionalboundary-writer.cc

@@ -35,23 +35,27 @@ template <class ProgramState, class GridView>
 template <class Friction>
 void FrictionalBoundaryWriter<ProgramState, GridView>::write(
     ProgramState const &programState, Friction &friction) {
+    // TODO
   auto const frictionalBoundaryDisplacements =
-      restrictToSurface(programState.u, frictionalBoundary_);
+      restrictToSurface(programState.u[0], frictionalBoundary_);
   addEntry(frictionalBoundaryDisplacementWriter_, programState.timeStep,
            frictionalBoundaryDisplacements);
 
   auto const frictionalBoundaryVelocities =
-      restrictToSurface(programState.v, frictionalBoundary_);
+          // TODO
+      restrictToSurface(programState.v[0], frictionalBoundary_);
   addEntry(frictionalBoundaryVelocityWriter_, programState.timeStep,
            frictionalBoundaryVelocities);
 
   auto const frictionalBoundaryStates =
-      restrictToSurface(programState.alpha, frictionalBoundary_);
+      // TODO
+      restrictToSurface(programState.alpha[0], frictionalBoundary_);
   addEntry(frictionalBoundaryStateWriter_, programState.timeStep,
            frictionalBoundaryStates);
 
-  friction.updateAlpha(programState.alpha);
-  auto const c = friction.coefficientOfFriction(programState.v);
+  // TODO
+  friction.updateAlpha(programState.alpha[0]);
+  auto const c = friction.coefficientOfFriction(programState.v[0]);
   auto const frictionalBoundaryCoefficient =
       restrictToSurface(c, frictionalBoundary_);
   addEntry(frictionalBoundaryCoefficientWriter_, programState.timeStep,

src/hdf5/frictionalboundary-writer_tmpl.cc

-#include "../explicitvectors.hh"
-#include "../explicitgrid.hh"
-
-#include "../program_state.hh"
-
-using MyProgramState = ProgramState<Vector, ScalarVector>;
-using MyFriction = GlobalFriction<Matrix, Vector>;
-
-template class FrictionalBoundaryWriter<MyProgramState, GridView>;
-template void FrictionalBoundaryWriter<MyProgramState, GridView>::write(
-    MyProgramState const &programState, MyFriction &friction);

src/hdf5/patchinfo-writer.cc

@@ -87,7 +87,8 @@ PatchInfoWriter<ProgramState, VertexBasis, GridView>::PatchInfoWriter(
 template <class ProgramState, class VertexBasis, class GridView>
 void PatchInfoWriter<ProgramState, VertexBasis, GridView>::write(
     ProgramState const &programState) {
-  BasisGridFunction<VertexBasis, Vector> velocity(vertexBasis_, programState.v);
+    // TODO
+  BasisGridFunction<VertexBasis, Vector> velocity(vertexBasis_, programState.v[0]);
   auto const gridVelocity = gridEvaluator_.evaluate(velocity);
   addEntry(weakPatchGridVelocityWriter_, programState.timeStep, gridVelocity);
 }

src/hdf5/patchinfo-writer_tmpl.cc

@@ -4,11 +4,11 @@
 #include "../program_state.hh"
 
 using MyProgramState = ProgramState<Vector, ScalarVector>;
-using P1Basis = P1NodalBasis<GridView, double>;
+using P1Basis = P1NodalBasis<LeafGridView, double>;
 using MyFunction = BasisGridFunction<P1Basis, Vector>;
 
-template class GridEvaluator<LocalVector, GridView>;
+template class GridEvaluator<LocalVector, LeafGridView>;
 template Dune::Matrix<typename MyFunction::RangeType>
-GridEvaluator<LocalVector, GridView>::evaluate(MyFunction const &f) const;
+GridEvaluator<LocalVector, LeafGridView>::evaluate(MyFunction const &f) const;
 
-template class PatchInfoWriter<MyProgramState, P1Basis, GridView>;
+template class PatchInfoWriter<MyProgramState, P1Basis, LeafGridView>;

src/hdf5/restart-io.cc

@@ -19,12 +19,13 @@ RestartIO<ProgramState>::RestartIO(HDF5::Grouplike &group, size_t vertexCount)
 
 template <class ProgramState>
 void RestartIO<ProgramState>::write(ProgramState const &programState) {
-  addEntry(displacementWriter_, programState.timeStep, programState.u);
-  addEntry(velocityWriter_, programState.timeStep, programState.v);
-  addEntry(accelerationWriter_, programState.timeStep, programState.a);
-  addEntry(stateWriter_, programState.timeStep, programState.alpha);
+    // TODO
+  addEntry(displacementWriter_, programState.timeStep, programState.u[0]);
+  addEntry(velocityWriter_, programState.timeStep, programState.v[0]);
+  addEntry(accelerationWriter_, programState.timeStep, programState.a[0]);
+  addEntry(stateWriter_, programState.timeStep, programState.alpha[0]);
   addEntry(weightedNormalStressWriter_, programState.timeStep,
-           programState.weightedNormalStress);
+           programState.weightedNormalStress[0]);
   addEntry(relativeTimeWriter_, programState.timeStep,
            programState.relativeTime);
   addEntry(relativeTimeIncrementWriter_, programState.timeStep,
@@ -35,12 +36,13 @@ template <class ProgramState>
 void RestartIO<ProgramState>::read(size_t timeStep,
                                    ProgramState &programState) {
   programState.timeStep = timeStep;
-  readEntry(displacementWriter_, timeStep, programState.u);
-  readEntry(velocityWriter_, timeStep, programState.v);
-  readEntry(accelerationWriter_, timeStep, programState.a);
-  readEntry(stateWriter_, timeStep, programState.alpha);
+  // TODO
+  readEntry(displacementWriter_, timeStep, programState.u[0]);
+  readEntry(velocityWriter_, timeStep, programState.v[0]);
+  readEntry(accelerationWriter_, timeStep, programState.a[0]);
+  readEntry(stateWriter_, timeStep, programState.alpha[0]);
   readEntry(weightedNormalStressWriter_, timeStep,
-            programState.weightedNormalStress);
+            programState.weightedNormalStress[0]);
   readEntry(relativeTimeWriter_, timeStep, programState.relativeTime);
   readEntry(relativeTimeIncrementWriter_, timeStep, programState.relativeTau);
 }

src/hdf5/surface-writer.cc

@@ -24,11 +24,13 @@ SurfaceWriter<ProgramState, GridView>::SurfaceWriter(
 template <class ProgramState, class GridView>
 void SurfaceWriter<ProgramState, GridView>::write(
     ProgramState const &programState) {
-  auto const surfaceDisplacements = restrictToSurface(programState.u, surface_);
+    // TODO
+  auto const surfaceDisplacements = restrictToSurface(programState.u[0], surface_);
   addEntry(surfaceDisplacementWriter_, programState.timeStep,
            surfaceDisplacements);
 
-  auto const surfaceVelocities = restrictToSurface(programState.v, surface_);
+  // TODO
+  auto const surfaceVelocities = restrictToSurface(programState.v[0], surface_);
   addEntry(surfaceVelocityWriter_, programState.timeStep, surfaceVelocities);
 }
 

src/hdf5/surface-writer_tmpl.cc

@@ -5,4 +5,4 @@
 
 using MyProgramState = ProgramState<Vector, ScalarVector>;
 
-template class SurfaceWriter<MyProgramState, GridView>;
+template class SurfaceWriter<MyProgramState, LeafGridView>;

src/matrices.hh

 #ifndef SRC_MATRICES_HH
 #define SRC_MATRICES_HH
 
-template <class Matrix> struct Matrices {
-  Matrix elasticity;
-  Matrix damping;
-  Matrix mass;
+template <class Matrix> class Matrices {
+public:
+  std::vector<std::shared_ptr<Matrix>> elasticity;
+  std::vector<std::shared_ptr<Matrix>> damping;
+  std::vector<std::shared_ptr<Matrix>> mass;
+
+  Matrices(size_t n) {
+      elasticity.resize(n);
+      damping.resize(n);
+      mass.resize(n);
+  }
 };
 #endif

src/multi-body-problem-2D.cfg

+# -*- mode:conf -*-
+[boundary.friction]
+smallestDiameter= 2e-3  # [m]
+
+[timeSteps]
+refinementTolerance = 1e-5
+
+[u0.solver]
+tolerance         = 1e-8
+
+[a0.solver]
+tolerance         = 1e-8
+
+[v.solver]
+tolerance         = 1e-8
+
+[v.fpi]
+tolerance         = 1e-5
+
+[solver.tnnmg.linear]
+tolerance          = 1e-10

src/multi-body-problem-3D.cfg

+# -*- mode:conf -*-
+[boundary.friction]
+smallestDiameter= 2e-2  # [m]
+
+[boundary.friction.weakening]
+patchType       = Trapezoidal
+
+[timeSteps]
+refinementTolerance = 1e-5
+
+[u0.solver]
+tolerance         = 1e-6
+
+[a0.solver]
+tolerance         = 1e-6
+
+[v.solver]
+tolerance         = 1e-6
+
+[v.fpi]
+tolerance         = 1e-5
+
+[solver.tnnmg.linear]
+tolerance          = 1e-10

src/multi-body-problem-data/bc.hh

+#ifndef SRC_ONE_BODY_PROBLEM_DATA_BC_HH
+#define SRC_ONE_BODY_PROBLEM_DATA_BC_HH
+
+class VelocityDirichletCondition
+    : public Dune::VirtualFunction<double, double> {
+  void evaluate(double const &relativeTime, double &y) const {
+    // Assumed to vanish at time zero
+    double const finalVelocity = -5e-5;
+    y = (relativeTime <= 0.1)
+            ? finalVelocity * (1.0 - std::cos(relativeTime * M_PI / 0.1)) / 2.0
+            : finalVelocity;
+  }
+};
+
+class NeumannCondition : public Dune::VirtualFunction<double, double> {
+  void evaluate(double const &relativeTime, double &y) const { y = 0.0; }
+};
+#endif

src/multi-body-problem-data/cuboidgeometry.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <fstream>
+
+#ifdef HAVE_CAIROMM
+#include <cairomm/context.h>
+#include <cairomm/fontface.h>
+#include <cairomm/surface.h>
+#endif
+
+#include "cuboidgeometry.hh"
+
+/*
+const CuboidGeometry::LocalVector CuboidGeometry::lift(const double v0, const double v1) {
+  vc = 0;
+  vc[0] = vc2D[0];
+  vc[1] = vc2D[1];
+}
+
+const CuboidGeometry::LocalVector& CuboidGeometry::A() const { return A_; }
+const CuboidGeometry::LocalVector& CuboidGeometry::B() const { return B_; }
+const CuboidGeometry::LocalVector& CuboidGeometry::C() const { return C_; }
+const CuboidGeometry::LocalVector& CuboidGeometry::D() const { return D_; }
+const CuboidGeometry::LocalVector& CuboidGeometry::X() const { return X_; }
+const CuboidGeometry::LocalVector& CuboidGeometry::Y() const { return Y_; }
+double CuboidGeometry::depth() const { return depth_; }
+
+void CuboidGeometry::setupWeak(const LocalVector& weakOrigin, const double weakLength) {
+    lift(weakOrigin, X_);
+    const LocalVector2D Y({X_[0]+weakLength, X_[1]});
+    lift(Y, Y_);
+}
+*/
+
+#if MY_DIM == 3
+CuboidGeometry::CuboidGeometry(const LocalVector& origin, const LocalVector& weakOrigin, const double weakLength, const double length, const double width, const double depth_):
+    length_(length*lengthScale),
+    width_(width*lengthScale),
+    A(origin),
+    B({origin[0]+length_, origin[1], 0}),
+    C({origin[0]+length_, origin[1]+width_, 0}),
+    D({origin[0], origin[1]+width_, 0}),
+    X(weakOrigin),
+    Y({X[0]+weakLength, X[1], 0}),
+    depth(depth_*lengthScale) {}
+#else
+CuboidGeometry::CuboidGeometry(const LocalVector& origin, const LocalVector& weakOrigin, const double weakLength, const double length, const double width):
+    length_(length*lengthScale),
+    width_(width*lengthScale),
+    A(origin),
+    B({origin[0]+length_, origin[1]}),
+    C({origin[0]+length_, origin[1]+width_}),
+    D({origin[0], origin[1]+width_}),
+    X(weakOrigin),
+    Y({X[0]+weakLength, X[1]}) {}
+#endif
+
+
+void CuboidGeometry::write() const {
+  std::fstream writer("geometry", std::fstream::out);
+  writer << "A = " << A << std::endl;
+  writer << "B = " << B << std::endl;
+  writer << "C = " << C << std::endl;
+  writer << "D = " << D << std::endl;
+  writer << "X = " << X << std::endl;
+  writer << "Y = " << Y << std::endl;
+}
+
+void CuboidGeometry::render() const {
+#ifdef HAVE_CAIROMM
+  std::string const filename = "geometry.png";
+  double const width = 600;
+  double const height = 400;
+  double const widthScale = 400;
+  double const heightScale = 400;
+
+  auto surface =
+      Cairo::ImageSurface::create(Cairo::FORMAT_ARGB32, width, height);
+  auto cr = Cairo::Context::create(surface);
+
+  auto const setRGBColor = [&](int colour) {
+    cr->set_source_rgb(((colour & 0xFF0000) >> 16) / 255.0,
+                       ((colour & 0x00FF00) >> 8) / 255.0,
+                       ((colour & 0x0000FF) >> 0) / 255.0);
+  };
+  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);
+  cr->set_line_width(0.0025);
+
+  // triangle
+  {
+    moveTo(reference::A);
+    lineTo(reference::B);
+    lineTo(reference::C);
+    cr->close_path();
+    cr->stroke();
+  }
+
+  // dashed lines
+  {
+    cr->save();
+    std::vector<double> dashPattern = { 0.005 };
+    cr->set_dash(dashPattern, 0);
+    moveTo(reference::Z);
+    lineTo(reference::Y);
+    moveTo(reference::U);
+    lineTo(reference::X);
+    cr->stroke();
+    cr->restore();
+  }
+
+  // fill viscoelastic region
+  {
+    cr->save();
+    setRGBColor(0x0097E0);
+    moveTo(reference::B);
+    lineTo(reference::K);
+    lineTo(reference::M);
+    cr->fill();
+    cr->restore();
+  }
+
+  // mark weakening region
+  {
+    cr->save();
+    setRGBColor(0x7AD3FF);
+    cr->set_line_width(0.005);
+    moveTo(reference::X);
+    lineTo(reference::Y);
+    cr->stroke();
+    cr->restore();
+  }
+
+  // mark points
+  {
+    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();
+    };
+
+    cr->save();
+    setRGBColor(0x002F47);
+    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();
+  }
+
+  // labels
+  {
+    auto const label = [&](LocalVector2D const &v, std::string l) {
+      moveTo(v);
+      cr->rel_move_to(0.005, -0.02);
+      cr->show_text(l);
+    };
+    auto font = Cairo::ToyFontFace::create(
+        "monospace", Cairo::FONT_SLANT_NORMAL, Cairo::FONT_WEIGHT_NORMAL);
+
+    cr->save();
+    cr->set_font_face(font);
+    cr->set_font_size(0.03);
+
+    label(A, "A");
+    label(B, "B");
+    label(C, "C");
+    label(D, "D");
+    label(X, "X");
+    label(Y, "Y");
+
+    cr->restore();
+  }
+
+  surface->write_to_png(filename);
+#endif
+}

src/multi-body-problem-data/cuboidgeometry.hh

+#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)}};
+  */