Make assemble_frictional a factory method

src/one-body-sample.cc

@@ -122,21 +122,21 @@ void assemble_neumann(GridView const &gridView, FEBasis const &feBasis,
 
 // Assembles constant 1-function on frictional boundary in nodalIntegrals
 template <class GridType, class GridView, class LocalVectorType, class FEBasis>
-void assemble_frictional(
-    GridView const &gridView, FEBasis const &feBasis,
-    Dune::BitSetVector<1> const &frictionalNodes,
-    Dune::BlockVector<Dune::FieldVector<double, 1>> &nodalIntegrals) {
+Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
+assemble_frictional(GridView const &gridView, FEBasis const &feBasis,
+                    Dune::BitSetVector<1> const &frictionalNodes) {
   BoundaryPatch<GridView> frictionalBoundary(gridView, frictionalNodes);
   ConstantFunction<LocalVectorType, Dune::FieldVector<double, 1>>
   constantOneFunction(1);
   NeumannBoundaryAssembler<GridType, Dune::FieldVector<double, 1>>
   frictionalBoundaryAssembler(constantOneFunction);
 
-  BoundaryFunctionalAssembler<FEBasis> boundaryFunctionalAssembler(
-      feBasis, frictionalBoundary);
-  boundaryFunctionalAssembler.assemble(
-      frictionalBoundaryAssembler, nodalIntegrals,
-      true); // whether to resize the output vector and zero all of its entries
+  auto nodalIntegrals =
+      Dune::make_shared<Dune::BlockVector<Dune::FieldVector<double, 1>>>();
+  BoundaryFunctionalAssembler<FEBasis>(feBasis, frictionalBoundary)
+      .assemble(frictionalBoundaryAssembler, *nodalIntegrals,
+                true); // resize and zero output vector
+  return nodalIntegrals;
 }
 
 template <class VectorType, class MatrixType>
@@ -267,9 +267,8 @@ int main(int argc, char *argv[]) {
     VectorType b4;
 
     auto nodalIntegrals =
-        Dune::make_shared<Dune::BlockVector<Dune::FieldVector<double, 1>>>();
-    assemble_frictional<GridType, GridView, SmallVector, P1Basis>(
-        leafView, p1Basis, frictionalNodes, *nodalIntegrals);
+        assemble_frictional<GridType, GridView, SmallVector, P1Basis>(
+            leafView, p1Basis, frictionalNodes);
 
     auto myGlobalNonlinearity = assemble_nonlinearity<VectorType, OperatorType>(
         grid->size(grid->maxLevel(), dim), parset, nodalIntegrals);