[Extend] Introduce geometry class

src/mygeometry.hh

+#ifndef MY_GEOMETRY_HH
+#define MY_GEOMETRY_HH
+
+struct MyGeometry {
+  MyGeometry() {}
+
+  Dune::FieldVector<double, 2> A = { 0, 0 };
+  Dune::FieldVector<double, 2> B = { 5, 0 };
+  Dune::FieldVector<double, 2> C = { 5, 1 };
+  Dune::FieldVector<double, 2> D = { 0, 1 };
+
+  Dune::FieldVector<double, 2> zenith = { 0, 1 };
+};
+
+#endif

src/mygrid.hh

+#include <dune/grid/common/gridfactory.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+
+#include "mygeometry.hh"
+
+template <class Grid>
+std::shared_ptr<Grid> constructGrid(MyGeometry const &myGeometry) {
+  std::array<unsigned int, 2> elements = { { 5, 1 } };
+  return Dune::StructuredGridFactory<Grid>::createSimplexGrid(
+      myGeometry.A, myGeometry.C, elements);
+}
+
+template <class GridView, class MyGeometry> class MyFaces {
+private:
+  bool isClose(double a, double b) {
+    return std::abs(a - b) < 1e-14;
+  };
+
+public:
+  MyFaces(GridView const &gridView, MyGeometry const &myGeometry)
+      : lower(gridView), upper(gridView) {
+    lower.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return isClose(myGeometry.A[1], in.geometry().center()[1]);
+    });
+
+    upper.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return isClose(myGeometry.C[1], in.geometry().center()[1]);
+    });
+  }
+
+  BoundaryPatch<GridView> lower;
+  BoundaryPatch<GridView> upper;
+};

src/one-body-sample.cc

@@ -72,6 +72,8 @@
 #include "enum_verbosity.cc"
 #include "enums.hh"
 #include "friction_writer.hh"
+#include "mygeometry.hh"
+#include "mygrid.hh"
 #include "solverfactory.hh"
 #include "state.hh"
 #include "timestepping.hh"
@@ -93,6 +95,8 @@ int main(int argc, char *argv[]) {
                                            parset);
     Dune::ParameterTreeParser::readOptions(argc, argv, parset);
 
+    MyGeometry const myGeometry;
+
     auto const youngModulus = parset.get<double>("body.youngModulus");
     auto const poissonRatio = parset.get<double>("body.poissonRatio");
     auto const shearViscosity = parset.get<double>("body.shearViscosity");
@@ -102,23 +106,7 @@ int main(int argc, char *argv[]) {
 
     // {{{ Set up grid
     using Grid = Dune::ALUGrid<dims, dims, Dune::simplex, Dune::nonconforming>;
-
-    Dune::FieldVector<size_t, dims> lowerLeft(0);
-    Dune::FieldVector<size_t, dims> upperRight(1);
-    upperRight[0] = 5;
-    upperRight[1] = 1;
-
-    std::shared_ptr<Grid> grid;
-    {
-      Dune::array<unsigned int, dims> elements;
-      for (size_t i = 0; i < dims; ++i)
-        elements[i] = upperRight[i]; // 1x1 elements
-
-      grid = Dune::StructuredGridFactory<Grid>::createSimplexGrid(
-          lowerLeft, upperRight, elements);
-    }
-    Dune::FieldVector<double, dims> zenith(0);
-    zenith[1] = 1;
+    auto grid = constructGrid<Grid>(myGeometry); // FIXME
 
     auto const refinements = parset.get<size_t>("grid.refinements");
     grid->globalRefine(refinements);
@@ -128,28 +116,14 @@ int main(int argc, char *argv[]) {
     GridView const leafView = grid->leafView();
     // }}}
 
-    // Set up faces
-    BoundaryPatch<GridView> lowerFace(leafView);
-    BoundaryPatch<GridView> upperFace(leafView);
-    {
-      auto const isClose = [](double a,
-                              double b) { return std::abs(a - b) < 1e-14; };
-      lowerFace.insertFacesByProperty([&](
-          typename GridView::Intersection const &in) {
-        return isClose(lowerLeft[1], in.geometry().center()[1]);
-      });
-
-      upperFace.insertFacesByProperty([&](
-          typename GridView::Intersection const &in) {
-        return isClose(upperRight[1], in.geometry().center()[1]);
-      });
-    }
+    // Set up myFaces
+    MyFaces<GridView, MyGeometry> myFaces(leafView, myGeometry);
 
     // Neumann boundary
     BoundaryPatch<GridView> const neumannBoundary(leafView);
 
     // Frictional Boundary
-    BoundaryPatch<GridView> const &frictionalBoundary = lowerFace;
+    BoundaryPatch<GridView> const &frictionalBoundary = myFaces.lower;
     Dune::BitSetVector<1> frictionalNodes(fineVertexCount);
     frictionalBoundary.getVertices(frictionalNodes);
 
@@ -157,10 +131,10 @@ int main(int argc, char *argv[]) {
     Dune::BitSetVector<dims> noNodes(fineVertexCount);
     Dune::BitSetVector<dims> dirichletNodes(fineVertexCount);
     for (size_t i = 0; i < fineVertexCount; ++i) {
-      if (lowerFace.containsVertex(i))
+      if (myFaces.lower.containsVertex(i))
         dirichletNodes[i][1] = true;
 
-      if (upperFace.containsVertex(i))
+      if (myFaces.upper.containsVertex(i))
         dirichletNodes[i] = true;
     }
 
@@ -204,7 +178,8 @@ int main(int argc, char *argv[]) {
 
     // Assemble forces
     Vector gravityFunctional;
-    myAssembler.assembleBodyForce(gravity, density, zenith, gravityFunctional);
+    myAssembler.assembleBodyForce(gravity, density, myGeometry.zenith,
+                                  gravityFunctional);
 
     // Problem formulation: right-hand side
     auto const computeExternalForces = [&](double _relativeTime, Vector &_ell) {
@@ -219,12 +194,12 @@ int main(int argc, char *argv[]) {
     {
       double volume = 1.0;
       for (size_t i = 0; i < dims; ++i)
-        volume *= (upperRight[i] - lowerLeft[i]);
+        volume *= (myGeometry.C[i] - myGeometry.A[i]);
 
       double area = 1.0;
       for (size_t i = 0; i < dims; ++i)
         if (i != 1)
-          area *= (upperRight[i] - lowerLeft[i]);
+          area *= (myGeometry.C[i] - myGeometry.A[i]);
 
       // volume * gravity * density / area    = normal stress
       // V      * g       * rho     / A       = sigma_n