Make Neumann boundary term assembly into a func.

src/one-body-sample.cc

@@ -114,6 +114,29 @@ void setup_boundary(GridType const &grid,
   std::cout << "Number of nodes: " << overall_nodes << std::endl;
 }
 
+// Assembles Neumann boundary term in f
+template <class GridType, class LocalVectorType, class FEBasis>
+void assemble_neumann(GridType const &grid, FEBasis const &feBasis,
+                      Dune::BitSetVector<1> const &neumannNodes,
+                      Dune::BlockVector<LocalVectorType> &
+                          f) { // constant sample function on neumann boundary
+  BoundaryPatch<typename GridType::LeafGridView> neumannBoundary(
+      grid.leafView(), neumannNodes);
+  LocalVectorType SampleVector;
+  // FIXME: random values
+  SampleVector[0] = 1;
+  SampleVector[1] = 2;
+  ConstantFunction<LocalVectorType, LocalVectorType> fNeumann(SampleVector);
+  NeumannBoundaryAssembler<GridType, LocalVectorType> neumannBoundaryAssembler(
+      fNeumann);
+
+  BoundaryFunctionalAssembler<FEBasis> boundaryFunctionalAssembler(
+      feBasis, neumannBoundary);
+  boundaryFunctionalAssembler.assemble(
+      neumannBoundaryAssembler, f,
+      true); // resize the output vector and zero all of its entries
+}
+
 int main() {
   try {
     typedef Dune::FieldVector<double, dim> SmallVector;
@@ -176,22 +199,8 @@ int main() {
     normalStress.resize(grid.size(grid.maxLevel(), dim));
     std::fill(normalStress.begin(), normalStress.end(), 0.0);
 
-    { // constant sample function on neumann boundary
-      BoundaryPatch<GridType::LeafGridView> neumannBoundary(grid.leafView(),
-                                                            neumannNodes);
-      SmallVector SampleVector;
-      SampleVector[0] = 1;
-      SampleVector[1] = 2;
-      ConstantFunction<SmallVector, SmallVector> fNeumann(SampleVector);
-      NeumannBoundaryAssembler<GridType, SmallVector> neumannBoundaryAssembler(
-          fNeumann);
-
-      BoundaryFunctionalAssembler<P1Basis> boundaryFunctionalAssembler(
-          p1Basis, neumannBoundary);
-      boundaryFunctionalAssembler.assemble(
-          neumannBoundaryAssembler, f,
-          true); // resize the output vector and zero all of its entries
-    }
+    assemble_neumann<GridType, SmallVector, P1Basis>(grid, p1Basis,
+                                                     neumannNodes, f);
 
     {
       // constant 1-function on frictional boundary