Move assembler invocations to separate file

src/assemblers.cc

+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/functions/constantfunction.hh>
+#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
+#include <dune/fufem/assemblers/boundaryfunctionalassembler.hh>
+
+#include <dune/tectonic/globallaursennonlinearity.hh>
+#include <dune/tectonic/globalruinanonlinearity.hh>
+
+// Assembles Neumann boundary term in f
+template <class GridType, class GridView, class LocalVectorType, class FEBasis>
+void assemble_neumann(GridView const &gridView, FEBasis const &feBasis,
+                      Dune::BitSetVector<1> const &neumannNodes,
+                      Dune::BlockVector<LocalVectorType> &f,
+                      Dune::VirtualFunction<double, double> const &neumann,
+                      double time) { // constant sample function on neumann
+                                     // boundary
+  BoundaryPatch<GridView> const neumannBoundary(gridView, neumannNodes);
+  LocalVectorType SampleVector(0);
+  neumann.evaluate(time, SampleVector[0]);
+  SampleVector[1] = 0;
+  ConstantFunction<LocalVectorType, LocalVectorType> const fNeumann(
+      SampleVector);
+  NeumannBoundaryAssembler<GridType, LocalVectorType> neumannBoundaryAssembler(
+      fNeumann);
+
+  BoundaryFunctionalAssembler<FEBasis>(feBasis, neumannBoundary).assemble(
+      neumannBoundaryAssembler, f, true); // resize and zero output vector
+}
+
+// Assembles constant 1-function on frictional boundary in nodalIntegrals
+template <class GridType, class GridView, class LocalVectorType, class FEBasis>
+Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
+assemble_frictional(GridView const &gridView, FEBasis const &feBasis,
+                    Dune::BitSetVector<1> const &frictionalNodes) {
+  typedef Dune::FieldVector<double, 1> Singleton;
+  BoundaryPatch<GridView> const frictionalBoundary(gridView, frictionalNodes);
+  ConstantFunction<LocalVectorType, Singleton> const constantOneFunction(1);
+  NeumannBoundaryAssembler<GridType, Singleton> frictionalBoundaryAssembler(
+      constantOneFunction);
+
+  auto const nodalIntegrals = Dune::make_shared<Dune::BlockVector<Singleton>>();
+  BoundaryFunctionalAssembler<FEBasis>(feBasis, frictionalBoundary)
+      .assemble(frictionalBoundaryAssembler, *nodalIntegrals,
+                true); // resize and zero output vector
+  return nodalIntegrals;
+}
+
+template <class VectorType, class MatrixType>
+Dune::shared_ptr<Dune::GlobalNonlinearity<VectorType, MatrixType> const>
+assemble_nonlinearity(
+    int size, Dune::ParameterTree const &parset,
+    Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
+        nodalIntegrals,
+    Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>> state,
+    double h) {
+  typedef Dune::BlockVector<Dune::FieldVector<double, 1>> SingletonVectorType;
+  // {{{ Assemble terms for the nonlinearity
+  auto mu = Dune::make_shared<SingletonVectorType>(size);
+  *mu = parset.get<double>("boundary.friction.mu");
+
+  auto normalStress = Dune::make_shared<SingletonVectorType>(size);
+  *normalStress = parset.get<double>("boundary.friction.normalstress");
+
+  std::string const friction_model =
+      parset.get<std::string>("boundary.friction.model");
+  if (friction_model == std::string("Ruina")) {
+    auto a = Dune::make_shared<SingletonVectorType>(size);
+    *a = parset.get<double>("boundary.friction.ruina.a");
+
+    auto eta = Dune::make_shared<SingletonVectorType>(size);
+    *eta = parset.get<double>("boundary.friction.eta");
+
+    auto b = Dune::make_shared<SingletonVectorType>(size);
+    *b = parset.get<double>("boundary.friction.ruina.b");
+
+    auto L = Dune::make_shared<SingletonVectorType>(size);
+    *L = parset.get<double>("boundary.friction.ruina.L");
+
+    return Dune::make_shared<
+        Dune::GlobalRuinaNonlinearity<VectorType, MatrixType> const>(
+        nodalIntegrals, a, mu, eta, normalStress, b, state, L, h);
+  } else if (friction_model == std::string("Laursen")) {
+    return
+        // TODO: take state and h into account
+        Dune::make_shared<
+            Dune::GlobalLaursenNonlinearity<Dune::LinearFunction, VectorType,
+                                            MatrixType> const>(mu, normalStress,
+                                                               nodalIntegrals);
+  } else {
+    assert(false);
+  }
+}

src/assemblers.hh

+#ifndef ASSEMBLERS_HH
+#define ASSEMBLERS_HH
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/shared_ptr.hh>
+#include <dune/istl/bvector.hh>
+
+#include <dune/tectonic/globalnonlinearity.hh>
+
+template <class GridType, class GridView, class LocalVectorType, class FEBasis>
+void assemble_neumann(GridView const &gridView, FEBasis const &feBasis,
+                      Dune::BitSetVector<1> const &neumannNodes,
+                      Dune::BlockVector<LocalVectorType> &f,
+                      Dune::VirtualFunction<double, double> const &neumann,
+                      double time);
+
+template <class GridType, class GridView, class LocalVectorType, class FEBasis>
+Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
+assemble_frictional(GridView const &gridView, FEBasis const &feBasis,
+                    Dune::BitSetVector<1> const &frictionalNodes);
+
+template <class VectorType, class MatrixType>
+Dune::shared_ptr<Dune::GlobalNonlinearity<VectorType, MatrixType> const>
+assemble_nonlinearity(
+    int size, Dune::ParameterTree const &parset,
+    Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
+        nodalIntegrals,
+    Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>> state,
+    double h);
+
+#include "assemblers.cc"
+#endif

src/one-body-sample.cc

@@ -41,17 +41,13 @@
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
 
-#include <dune/fufem/assemblers/boundaryfunctionalassembler.hh>
 #include <dune/fufem/assemblers/functionalassembler.hh>
-#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
 #include <dune/fufem/assemblers/operatorassembler.hh>
 #include <dune/fufem/assemblers/transferoperatorassembler.hh>
-#include <dune/fufem/boundarypatch.hh>
 #include <dune/fufem/dunepython.hh>
 #include <dune/fufem/functions/basisgridfunction.hh>
-#include <dune/fufem/functions/constantfunction.hh>
 #include <dune/fufem/functions/vtkbasisgridfunction.hh>
 #include <dune/fufem/functionspacebases/p0basis.hh>
 #include <dune/fufem/functionspacebases/p1nodalbasis.hh>
@@ -66,12 +62,12 @@
 #include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
 #include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
 
-#include <dune/tectonic/globallaursennonlinearity.hh>
-#include <dune/tectonic/globalruinanonlinearity.hh>
 #include <dune/tectonic/myblockproblem.hh>
 #include <dune/tectonic/myconvexproblem.hh>
 #include <dune/tectonic/mydirectionalconvexfunction.hh>
 
+#include "assemblers.hh"
+
 int const dim = 2;
 
 template <class GridView>
@@ -112,91 +108,6 @@ void setup_boundary(GridView const &gridView,
   std::cout << "Number of frictional nodes: " << dirichlet_nodes << std::endl;
 }
 
-// Assembles Neumann boundary term in f
-template <class GridType, class GridView, class LocalVectorType, class FEBasis>
-void assemble_neumann(GridView const &gridView, FEBasis const &feBasis,
-                      Dune::BitSetVector<1> const &neumannNodes,
-                      Dune::BlockVector<LocalVectorType> &f,
-                      Dune::VirtualFunction<double, double> const &neumann,
-                      double time) { // constant sample function on neumann
-                                     // boundary
-  BoundaryPatch<GridView> const neumannBoundary(gridView, neumannNodes);
-  LocalVectorType SampleVector(0);
-  neumann.evaluate(time, SampleVector[0]);
-  SampleVector[1] = 0;
-  ConstantFunction<LocalVectorType, LocalVectorType> const fNeumann(
-      SampleVector);
-  NeumannBoundaryAssembler<GridType, LocalVectorType> neumannBoundaryAssembler(
-      fNeumann);
-
-  BoundaryFunctionalAssembler<FEBasis>(feBasis, neumannBoundary).assemble(
-      neumannBoundaryAssembler, f, true); // resize and zero output vector
-}
-
-// Assembles constant 1-function on frictional boundary in nodalIntegrals
-template <class GridType, class GridView, class LocalVectorType, class FEBasis>
-Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
-assemble_frictional(GridView const &gridView, FEBasis const &feBasis,
-                    Dune::BitSetVector<1> const &frictionalNodes) {
-  typedef Dune::FieldVector<double, 1> Singleton;
-  BoundaryPatch<GridView> const frictionalBoundary(gridView, frictionalNodes);
-  ConstantFunction<LocalVectorType, Singleton> const constantOneFunction(1);
-  NeumannBoundaryAssembler<GridType, Singleton> frictionalBoundaryAssembler(
-      constantOneFunction);
-
-  auto const nodalIntegrals = Dune::make_shared<Dune::BlockVector<Singleton>>();
-  BoundaryFunctionalAssembler<FEBasis>(feBasis, frictionalBoundary)
-      .assemble(frictionalBoundaryAssembler, *nodalIntegrals,
-                true); // resize and zero output vector
-  return nodalIntegrals;
-}
-
-template <class VectorType, class MatrixType>
-Dune::shared_ptr<Dune::GlobalNonlinearity<VectorType, MatrixType> const>
-assemble_nonlinearity(
-    int size, Dune::ParameterTree const &parset,
-    Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
-        nodalIntegrals,
-    Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>> state,
-    double h) {
-  typedef Dune::BlockVector<Dune::FieldVector<double, 1>> SingletonVectorType;
-  // {{{ Assemble terms for the nonlinearity
-  auto mu = Dune::make_shared<SingletonVectorType>(size);
-  *mu = parset.get<double>("boundary.friction.mu");
-
-  auto normalStress = Dune::make_shared<SingletonVectorType>(size);
-  *normalStress = parset.get<double>("boundary.friction.normalstress");
-
-  std::string const friction_model =
-      parset.get<std::string>("boundary.friction.model");
-  if (friction_model == std::string("Ruina")) {
-    auto a = Dune::make_shared<SingletonVectorType>(size);
-    *a = parset.get<double>("boundary.friction.ruina.a");
-
-    auto eta = Dune::make_shared<SingletonVectorType>(size);
-    *eta = parset.get<double>("boundary.friction.eta");
-
-    auto b = Dune::make_shared<SingletonVectorType>(size);
-    *b = parset.get<double>("boundary.friction.ruina.b");
-
-    auto L = Dune::make_shared<SingletonVectorType>(size);
-    *L = parset.get<double>("boundary.friction.ruina.L");
-
-    return Dune::make_shared<
-        Dune::GlobalRuinaNonlinearity<VectorType, MatrixType> const>(
-        nodalIntegrals, a, mu, eta, normalStress, b, state, L, h);
-  } else if (friction_model == std::string("Laursen")) {
-    return
-        // TODO: take state and h into account
-        Dune::make_shared<
-            Dune::GlobalLaursenNonlinearity<Dune::LinearFunction, VectorType,
-                                            MatrixType> const>(mu, normalStress,
-                                                               nodalIntegrals);
-  } else {
-    assert(false);
-  }
-}
-
 // The nonlinearity exp(-x)
 class DecayingExponential {
 public: