Compute von Mises stress

src/one-body-sample.cc

@@ -22,12 +22,16 @@
 #include <dune/grid/common/mcmgmapper.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/boundarypatch.hh>
+#include <dune/fufem/functions/basisgridfunction.hh>
 #include <dune/fufem/functions/constantfunction.hh>
 #include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+#include <dune/fufem/functionspacebases/p0basis.hh>
 #include <dune/grid/io/file/vtk/vtkwriter.hh>
 
 #include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
@@ -158,8 +162,10 @@ int main() {
     GridView const leafView = grid.leafView();
     // }}}
 
-    // Set up nodal basis
+    // Set up bases
+    typedef P0Basis<GridType::LeafGridView, double> P0Basis;
     typedef P1NodalBasis<GridType::LeafGridView, double> P1Basis;
+    P0Basis const p0Basis(leafView);
     P1Basis const p1Basis(leafView);
 
     // Assemble elastic force on the body
@@ -231,13 +237,27 @@ int main() {
         solver.solve();
       }
 
+      typedef Dune::BlockVector<Dune::FieldVector<double, 1>> CellVectorType;
+
+      CellVectorType vonMisesStress(grid.size(grid.maxLevel(), dim - 1));
+      auto *displacement =
+          new BasisGridFunction<P1Basis, VectorType>(p1Basis, u1);
+      VonMisesStressAssembler<GridType> localStressAssembler(E, nu,
+                                                             displacement);
+      FunctionalAssembler<P0Basis>(p0Basis)
+          .assemble(localStressAssembler, vonMisesStress, true);
+
       {
         Dune::VTKWriter<GridType::LeafGridView> writer(leafView);
         std::string filename((boost::format("obs%d") % run).str());
-        Dune::shared_ptr<Dune::VTKBasisGridFunction<P1Basis, VectorType>> ptr(
-            new Dune::VTKBasisGridFunction<P1Basis, VectorType>(
-                p1Basis, u1, "displacement"));
-        writer.addVertexData(ptr);
+        Dune::shared_ptr<Dune::VTKBasisGridFunction<P1Basis, VectorType>>
+        displacement_ptr(new Dune::VTKBasisGridFunction<P1Basis, VectorType>(
+            p1Basis, u1, "displacement"));
+        Dune::shared_ptr<Dune::VTKBasisGridFunction<P0Basis, CellVectorType>>
+        vonmises_ptr(new Dune::VTKBasisGridFunction<P0Basis, CellVectorType>(
+            p0Basis, vonMisesStress, "stress"));
+        writer.addVertexData(displacement_ptr);
+        writer.addCellData(vonmises_ptr);
         writer.write(filename.c_str());
       }