[Cleanup] Move body info to a separate class

dune/tectonic/body.hh

+#ifndef BODY_HH
+#define BODY_HH
+
+template <int dimension> struct Body {
+  using ScalarFunction = Dune::VirtualFunction<
+      Dune::FieldVector<double, dimension>, Dune::FieldVector<double, 1>>;
+  using VectorField =
+      Dune::VirtualFunction<Dune::FieldVector<double, dimension>,
+                            Dune::FieldVector<double, dimension>>;
+
+  double virtual getPoissonRatio() const = 0;
+
+  double virtual getYoungModulus() const = 0;
+
+  double virtual getShearViscosity() const = 0;
+  double virtual getBulkViscosity() const = 0;
+
+  double virtual getDensity() const = 0;
+  double virtual getGravity() const = 0;
+};
+
+#endif

src/mybody.hh

+#ifndef MY_BODY_HH
+#define MY_BODY_HH
+
+#include <dune/common/fvector.hh>
+
+#include <dune/tectonic/body.hh>
+
+#include "mygeometry.hh"
+
+template <int dimension> class MyBody : public Body<dimension> {
+  using typename Body<dimension>::ScalarFunction;
+  using typename Body<dimension>::VectorField;
+
+public:
+  MyBody(Dune::ParameterTree const &parset, MyGeometry const &mygeometry)
+      : poissonRatio_(parset.get<double>("body.poissonRatio")),
+        youngModulus_(parset.get<double>("body.youngModulus")),
+        shearViscosity_(parset.get<double>("body.shearViscosity")),
+        bulkViscosity_(parset.get<double>("body.bulkViscosity")),
+        density_(parset.get<double>("body.density")),
+        gravity_(parset.get<double>("gravity")) {}
+
+  double getPoissonRatio() const override { return poissonRatio_; }
+
+  double getYoungModulus() const override { return youngModulus_; }
+
+  double getShearViscosity() const override { return shearViscosity_; }
+  double getBulkViscosity() const override { return bulkViscosity_; }
+
+  double getDensity() const override { return density_; }
+  double getGravity() const override { return gravity_; }
+
+private:
+  double const poissonRatio_;
+  double const youngModulus_;
+  double const shearViscosity_;
+  double const bulkViscosity_;
+  double const density_;
+  double const gravity_;
+};
+
+#endif

src/one-body-sample.cc

@@ -72,6 +72,7 @@
 #include "enum_verbosity.cc"
 #include "enums.hh"
 #include "friction_writer.hh"
+#include "mybody.hh"
 #include "mygeometry.hh"
 #include "mygrid.hh"
 #include "solverfactory.hh"
@@ -97,13 +98,6 @@ int main(int argc, char *argv[]) {
 
     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");
-    auto const bulkViscosity = parset.get<double>("body.bulkViscosity");
-    auto const density = parset.get<double>("body.density");
-    auto const gravity = parset.get<double>("gravity");
-
     // {{{ Set up grid
     using Grid = Dune::ALUGrid<dims, dims, Dune::simplex, Dune::nonconforming>;
     auto grid = constructGrid<Grid>(myGeometry); // FIXME
@@ -162,10 +156,14 @@ int main(int argc, char *argv[]) {
 
     MyAssembler myAssembler(leafView);
 
+    MyBody<dims> const body(parset, myGeometry);
+
     Matrix A, C, M;
-    myAssembler.assembleElasticity(youngModulus, poissonRatio, A);
-    myAssembler.assembleViscosity(shearViscosity, bulkViscosity, C);
-    myAssembler.assembleMass(density, M);
+    myAssembler.assembleElasticity(body.getYoungModulus(),
+                                   body.getPoissonRatio(), A);
+    myAssembler.assembleViscosity(body.getShearViscosity(),
+                                  body.getBulkViscosity(), C);
+    myAssembler.assembleMass(body.getDensity(), M);
     EnergyNorm<Matrix, Vector> const ANorm(A), MNorm(M);
     // Q: Does it make sense to weigh them in this manner?
     SumNorm<Vector> const AMNorm(1.0, ANorm, 1.0, MNorm);
@@ -178,8 +176,8 @@ int main(int argc, char *argv[]) {
 
     // Assemble forces
     Vector gravityFunctional;
-    myAssembler.assembleBodyForce(gravity, density, myGeometry.zenith,
-                                  gravityFunctional);
+    myAssembler.assembleBodyForce(body.getGravity(), body.getDensity(),
+                                  myGeometry.zenith, gravityFunctional);
 
     // Problem formulation: right-hand side
     auto const computeExternalForces = [&](double _relativeTime, Vector &_ell) {
@@ -204,7 +202,7 @@ int main(int argc, char *argv[]) {
       // volume * gravity * density / area    = normal stress
       // V      * g       * rho     / A       = sigma_n
       // m^d    * N/kg    * kg/m^d  / m^(d-1) = N/m^(d-1)
-      normalStress = -volume * gravity * density / area;
+      normalStress = -volume * body.getGravity() * body.getDensity() / area;
     }
     FrictionData const frictionData(parset.sub("boundary.friction"),
                                     normalStress);
@@ -438,8 +436,8 @@ int main(int argc, char *argv[]) {
 
       if (parset.get<bool>("io.writeVTK")) {
         ScalarVector stress;
-        myAssembler.assembleVonMisesStress(youngModulus, poissonRatio, u,
-                                           stress);
+        myAssembler.assembleVonMisesStress(body.getYoungModulus(),
+                                           body.getPoissonRatio(), u, stress);
         vtkWriter.write(u, v, alpha, stress);
       }
     }