Add gravity functional

src/one-body-sample.cc

@@ -47,12 +47,14 @@
 #include <dune/istl/bvector.hh>
 
 #include <dune/fufem/assemblers/functionalassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/massassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
 #include <dune/fufem/assemblers/operatorassembler.hh>
 #include <dune/fufem/dunepython.hh>
 #include <dune/fufem/functions/basisgridfunction.hh>
+#include <dune/fufem/functions/constantfunction.hh>
 #include <dune/fufem/functionspacebases/p0basis.hh>
 #include <dune/fufem/functionspacebases/p1nodalbasis.hh>
 #include <dune/fufem/sharedpointermap.hh>
@@ -167,6 +169,7 @@ int main(int argc, char *argv[]) {
     P1Basis const p1Basis(leafView);
 
     MatrixType massMatrix;
+    VectorType gravityFunctional;
     {
       timer.reset();
 
@@ -199,6 +202,16 @@ int main(int argc, char *argv[]) {
       if (parset.get<bool>("enable_timer"))
         std::cout << "Assembled mass matrix in " << timer.elapsed() << "s"
                   << std::endl;
+
+      // Compute gravitational body force
+      SmallVector weightedGravitationalDirection(0);
+      weightedGravitationalDirection[1] = -density * gravity;
+      ConstantFunction<SmallVector, SmallVector> const gravityFunction(
+          weightedGravitationalDirection);
+      L2FunctionalAssembler<GridType, SmallVector> gravityFunctionalAssembler(
+          gravityFunction);
+      FunctionalAssembler<P1Basis>(p1Basis)
+          .assemble(gravityFunctionalAssembler, gravityFunctional, true);
     }
 
     // Assemble elastic force on the body
@@ -296,6 +309,7 @@ int main(int argc, char *argv[]) {
         VectorType ell(finestSize);
         assemble_neumann<GridType, GridView, SmallVector, P1Basis>(
             leafView, p1Basis, neumannNodes, ell, neumannFunction, time);
+        ell += gravityFunctional;
 
         MatrixType problem_A;
         VectorType problem_rhs(finestSize);