[Cleanup] Outsource vonMisesStress assembly

src/assemblers.cc

@@ -5,6 +5,7 @@
 #include <dune/istl/scaledidmatrix.hh>
 
 #include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/functions/basisgridfunction.hh>
 #include <dune/fufem/functions/constantfunction.hh>
 #include <dune/fufem/assemblers/localassemblers/boundarymassassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
@@ -12,6 +13,7 @@
 #include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/viscosityassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
 
 #include <dune/tectonic/globalruinanonlinearity.hh>
 
@@ -19,8 +21,10 @@
 
 template <class GridView, int dimension>
 MyAssembler<GridView, dimension>::MyAssembler(GridView const &_gridView)
-    : vertexBasis(_gridView),
+    : cellBasis(_gridView),
+      vertexBasis(_gridView),
       gridView(_gridView),
+      cellAssembler(cellBasis, cellBasis),
       vertexAssembler(vertexBasis, vertexBasis) {}
 
 template <class GridView, int dimension>
@@ -111,4 +115,16 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
       frictionalNodes, weights, fd);
 }
 
+template <class GridView, int dimension>
+void MyAssembler<GridView, dimension>::assembleVonMisesStress(
+    double youngModulus, double poissonRatio, Vector const &u,
+    ScalarVector &stress) {
+  auto const gridDisplacement =
+      std::make_shared<BasisGridFunction<VertexBasis, Vector> const>(
+          vertexBasis, u);
+  VonMisesStressAssembler<Grid, LocalCellBasis> localStressAssembler(
+      youngModulus, poissonRatio, gridDisplacement);
+  cellAssembler.assembleFunctional(localStressAssembler, stress);
+}
+
 #include "assemblers_tmpl.cc"

src/assemblers.hh

@@ -7,6 +7,10 @@
 #include <dune/istl/bvector.hh>
 
 #include <dune/fufem/assemblers/assembler.hh>
+#pragma clang diagnostic push
+#pragma clang diagnostic ignored "-Wsign-compare"
+#include <dune/fufem/functionspacebases/p0basis.hh>
+#pragma clang diagnostic pop
 #include <dune/fufem/functionspacebases/p1nodalbasis.hh>
 
 #include <dune/tectonic/globalnonlinearity.hh>
@@ -19,15 +23,21 @@ template <class GridView, int dimension> class MyAssembler {
   using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
   using Vector = Dune::BlockVector<Dune::FieldVector<double, dimension>>;
 
+  using CellBasis = P0Basis<GridView, double>;
   using VertexBasis = P1NodalBasis<GridView, double>;
+
+  CellBasis cellBasis;
   VertexBasis vertexBasis;
 
 private:
   using Grid = typename GridView::Grid;
   using LocalVector = typename Vector::block_type;
+
+  using LocalCellBasis = typename CellBasis::LocalFiniteElement;
   using LocalVertexBasis = typename VertexBasis::LocalFiniteElement;
 
   GridView const &gridView;
+  Assembler<CellBasis, CellBasis> cellAssembler;
   Assembler<VertexBasis, VertexBasis> vertexAssembler;
 
 public:
@@ -56,6 +66,9 @@ template <class GridView, int dimension> class MyAssembler {
   assembleFrictionNonlinearity(
       BoundaryPatch<GridView> const &frictionalBoundary,
       FrictionData const &fd);
+
+  void assembleVonMisesStress(double youngModulus, double poissonRatio,
+                              Vector const &u, ScalarVector &stress);
 };
 
 #endif

src/one-body-sample.cc

@@ -49,17 +49,10 @@
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
 
-#include <dune/fufem/assemblers/assembler.hh>
-#include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
 #include <dune/fufem/boundarypatch.hh>
 #include <dune/fufem/dunepython.hh>
 #include <dune/fufem/functions/basisgridfunction.hh>
 
-#pragma clang diagnostic push
-#pragma clang diagnostic ignored "-Wsign-compare"
-#include <dune/fufem/functionspacebases/p0basis.hh>
-#pragma clang diagnostic pop
-
 #include <dune/fufem/functionspacebases/p1nodalbasis.hh>
 #include <dune/fufem/sharedpointermap.hh>
 #include <dune/solvers/norms/energynorm.hh>
@@ -145,11 +138,6 @@ int main(int argc, char *argv[]) {
     GridView const leafView = grid->leafView();
     // }}}
 
-    // Set up bases
-    using P0Basis = P0Basis<GridView, double>;
-    P0Basis const p0Basis(leafView);
-    Assembler<P0Basis, P0Basis> const p0Assembler(p0Basis, p0Basis);
-
     // Set up the boundary
     Dune::BitSetVector<dims> velocityDirichletNodes(fineVertexCount, false);
     Dune::BitSetVector<dims> const &displacementDirichletNodes =
@@ -470,16 +458,11 @@ int main(int argc, char *argv[]) {
       relaxationWriter << std::endl;
 
       if (parset.get<bool>("io.writeVTK")) {
-        auto const gridDisplacement = std::make_shared<
-            BasisGridFunction<typename MyAssembler::VertexBasis, Vector> const>(
-            myAssembler.vertexBasis, u);
         ScalarVector vonMisesStress;
-        VonMisesStressAssembler<Grid, P0Basis::LocalFiniteElement>
-        localStressAssembler(youngModulus, poissonRatio, gridDisplacement);
-        p0Assembler.assembleFunctional(localStressAssembler, vonMisesStress);
-
-        writeVtk(myAssembler.vertexBasis, u, alpha, p0Basis, vonMisesStress,
-                 (boost::format("obs%d") % run).str());
+        myAssembler.assembleVonMisesStress(youngModulus, poissonRatio, u,
+                                           vonMisesStress);
+        writeVtk(myAssembler.vertexBasis, u, alpha, myAssembler.cellBasis,
+                 vonMisesStress, (boost::format("obs%d") % run).str());
       }
     }
     iterationWriter.close();