Move assembly of nonlinearity to a function

src/one-body-sample.cc

@@ -135,6 +135,46 @@ void assemble_frictional(
       true); // whether to resize the output vector and zero all of its entries
 }
 
+void assemble_nonlinearity(
+    int size, Dune::ParameterTree const &parset,
+    Dune::shared_ptr<Dune::GlobalNonlinearity<dim>> &myGlobalNonlinearity,
+    std::vector<Dune::FieldVector<double, 1>> &nodalIntegrals) {
+  // {{{ Assemble terms for the nonlinearity
+  std::vector<double> mu;
+  mu.resize(size);
+  std::fill(mu.begin(), mu.end(), parset.get<double>("boundary.friction.mu"));
+
+  std::vector<double> normalStress;
+  normalStress.resize(size);
+  std::fill(normalStress.begin(), normalStress.end(),
+            parset.get<double>("boundary.friction.normalstress"));
+
+  std::string const friction_model =
+      parset.get<std::string>("boundary.friction.model");
+  if (friction_model == std::string("Ruina")) {
+    std::vector<double> a;
+    a.resize(size);
+    std::fill(a.begin(), a.end(),
+              parset.get<double>("boundary.friction.ruina.a"));
+
+    std::vector<double> eta;
+    eta.resize(size);
+    std::fill(eta.begin(), eta.end(),
+              parset.get<double>("boundary.friction.eta"));
+
+    auto const tmp = new Dune::GlobalRuinaNonlinearity<dim>(
+        nodalIntegrals, a, mu, eta, normalStress);
+    myGlobalNonlinearity = Dune::shared_ptr<Dune::GlobalNonlinearity<dim>>(tmp);
+  } else if (friction_model == std::string("Laursen")) {
+    auto const tmp =
+        new Dune::GlobalLaursenNonlinearity<dim, Dune::LinearFunction>(
+            mu, normalStress, nodalIntegrals);
+    myGlobalNonlinearity = Dune::shared_ptr<Dune::GlobalNonlinearity<dim>>(tmp);
+  } else {
+    assert(false);
+  }
+}
+
 int main(int argc, char *argv[]) {
   try {
     Dune::ParameterTree parset;
@@ -237,45 +277,9 @@ int main(int argc, char *argv[]) {
       stiffnessMatrix.mmv(u2, b2);
       stiffnessMatrix.mmv(u3, b3);
 
-      // {{{ Assemble terms for the nonlinearity
-      std::vector<double> mu;
-      mu.resize(grid.size(grid.maxLevel(), dim));
-      std::fill(mu.begin(), mu.end(),
-                parset.get<double>("boundary.friction.mu"));
-
-      std::vector<double> normalStress;
-      normalStress.resize(grid.size(grid.maxLevel(), dim));
-      std::fill(normalStress.begin(), normalStress.end(),
-                parset.get<double>("boundary.friction.normalstress"));
-
-      typedef Dune::shared_ptr<Dune::GlobalNonlinearity<dim>>
-      globalNonlinearityPtr;
-      globalNonlinearityPtr myGlobalNonlinearity;
-
-      std::string const friction_model =
-          parset.get<std::string>("boundary.friction.model");
-      if (friction_model == std::string("Ruina")) {
-        std::vector<double> a;
-        a.resize(grid.size(grid.maxLevel(), dim));
-        std::fill(a.begin(), a.end(),
-                  parset.get<double>("boundary.friction.ruina.a"));
-
-        std::vector<double> eta;
-        eta.resize(grid.size(grid.maxLevel(), dim));
-        std::fill(eta.begin(), eta.end(),
-                  parset.get<double>("boundary.friction.eta"));
-
-        auto const tmp = new Dune::GlobalRuinaNonlinearity<dim>(
-            nodalIntegrals, a, mu, eta, normalStress);
-        myGlobalNonlinearity = globalNonlinearityPtr(tmp);
-      } else if (friction_model == std::string("Laursen")) {
-        auto const tmp =
-            new Dune::GlobalLaursenNonlinearity<dim, Dune::LinearFunction>(
-                mu, normalStress, nodalIntegrals);
-        myGlobalNonlinearity = globalNonlinearityPtr(tmp);
-      } else {
-        assert(false);
-      }
+      Dune::shared_ptr<Dune::GlobalNonlinearity<dim>> myGlobalNonlinearity;
+      assemble_nonlinearity(grid.size(grid.maxLevel(), dim), parset,
+                            myGlobalNonlinearity, nodalIntegrals);
       // }}}
 
       if (parset.get<bool>("useNonlinearGS")) {