[Extend] (minor) Allow mass/density to vary in space

dune/tectonic/body.hh

@@ -15,8 +15,8 @@ template <int dimension> struct Body {
   double virtual getShearViscosity() const = 0;
   double virtual getBulkViscosity() const = 0;
 
-  double virtual getDensity() const = 0;
-  double virtual getGravity() const = 0;
+  ScalarFunction virtual const &getDensityField() const = 0;
+  VectorField virtual const &getGravityField() const = 0;
 };
 
 #endif

dune/tectonic/gravity.hh

+#ifndef GRAVITY_HH
+#define GRAVITY_HH
+
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+
+template <int dimension>
+class Gravity
+    : public Dune::VirtualFunction<Dune::FieldVector<double, dimension>,
+                                   Dune::FieldVector<double, dimension>> {
+public:
+  Gravity(
+      Dune::VirtualFunction<Dune::FieldVector<double, dimension>,
+                            Dune::FieldVector<double, 1>> const &_densityField,
+      Dune::FieldVector<double, dimension> const &_zenith, double _gravity)
+      : densityField(_densityField), zenith(_zenith), gravity(_gravity) {}
+
+  void evaluate(Dune::FieldVector<double, dimension> const &x,
+                Dune::FieldVector<double, dimension> &y) const {
+    y = zenith;
+
+    Dune::FieldVector<double, 1> density;
+    densityField.evaluate(x, density);
+
+    y *= -gravity * density;
+  }
+
+private:
+  Dune::VirtualFunction<Dune::FieldVector<double, dimension>,
+                        Dune::FieldVector<double, 1>> const &densityField;
+  Dune::FieldVector<double, dimension> const &zenith;
+  double const gravity;
+};
+#endif

src/assemblers.cc

@@ -6,11 +6,11 @@
 
 #include <dune/fufem/assemblers/localassemblers/boundarymassassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
-#include <dune/fufem/assemblers/localassemblers/massassembler.hh>
 #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/fufem/assemblers/localassemblers/weightedmassassembler.hh>
 #include <dune/fufem/boundarypatch.hh>
 #include <dune/fufem/functions/basisgridfunction.hh>
 #include <dune/fufem/functions/constantfunction.hh>
@@ -37,11 +37,18 @@ void MyAssembler<GridView, dimension>::assembleFrictionalBoundaryMass(
 }
 
 template <class GridView, int dimension>
-void MyAssembler<GridView, dimension>::assembleMass(double density, Matrix &M) {
-  MassAssembler<Grid, LocalVertexBasis, LocalVertexBasis,
-                Dune::ScaledIdentityMatrix<double, dimension>> const localMass;
-  vertexAssembler.assembleOperator(localMass, M);
-  M *= density;
+void MyAssembler<GridView, dimension>::assembleMass(
+    Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
+        densityFunction,
+    Matrix &M) {
+  // NOTE: We treat the weight as a constant function
+  QuadratureRuleKey quadKey(dimension, 0);
+
+  WeightedMassAssembler<Grid, LocalVertexBasis, LocalVertexBasis,
+                        Dune::VirtualFunction<LocalVector, LocalScalarVector>,
+                        Dune::ScaledIdentityMatrix<double, dimension>>
+  localWeightedMass(gridView.grid(), densityFunction, quadKey);
+  vertexAssembler.assembleOperator(localWeightedMass, M);
 }
 
 template <class GridView, int dimension>
@@ -62,18 +69,11 @@ void MyAssembler<GridView, dimension>::assembleViscosity(double shearViscosity,
 }
 
 template <class GridView, int dimension>
-void MyAssembler<GridView, dimension>::assembleBodyForce(double gravity,
-                                                         double density,
-                                                         LocalVector zenith,
-                                                         Vector &f) {
-  LocalVector weightedGravitationalDirection = zenith;
-  weightedGravitationalDirection *= density * gravity;
-  weightedGravitationalDirection *= -1;
-
-  ConstantFunction<LocalVector, LocalVector> const gravityFunction(
-      weightedGravitationalDirection);
+void MyAssembler<GridView, dimension>::assembleBodyForce(
+    Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
+    Vector &f) {
   L2FunctionalAssembler<Grid, LocalVertexBasis, LocalVector>
-  gravityFunctionalAssembler(gravityFunction);
+  gravityFunctionalAssembler(gravityField);
   vertexAssembler.assembleFunctional(gravityFunctionalAssembler, f);
 }
 

src/assemblers.hh

@@ -32,6 +32,7 @@ template <class GridView, int dimension> class MyAssembler {
 private:
   using Grid = typename GridView::Grid;
   using LocalVector = typename Vector::block_type;
+  using LocalScalarVector = typename ScalarVector::block_type;
 
   using LocalCellBasis = typename CellBasis::LocalFiniteElement;
   using LocalVertexBasis = typename VertexBasis::LocalFiniteElement;
@@ -47,15 +48,18 @@ template <class GridView, int dimension> class MyAssembler {
       BoundaryPatch<GridView> const &frictionalBoundary,
       ScalarMatrix &frictionalBoundaryMass);
 
-  void assembleMass(double density, Matrix &M);
+  void assembleMass(Dune::VirtualFunction<
+                        LocalVector, LocalScalarVector> const &densityFunction,
+                    Matrix &M);
 
   void assembleElasticity(double E, double nu, Matrix &A);
 
   void assembleViscosity(double shearViscosity, double bulkViscosity,
                          Matrix &C);
 
-  void assembleBodyForce(double gravity, double density, LocalVector zenith,
-                         Vector &f);
+  void assembleBodyForce(
+      Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
+      Vector &f);
 
   void assembleNeumann(BoundaryPatch<GridView> const &neumannBoundary,
                        Vector &f,

src/mybody.hh

@@ -3,7 +3,10 @@
 
 #include <dune/common/fvector.hh>
 
+#include <dune/fufem/functions/constantfunction.hh>
+
 #include <dune/tectonic/body.hh>
+#include <dune/tectonic/gravity.hh>
 
 #include "mygeometry.hh"
 
@@ -11,14 +14,18 @@ template <int dimension> class MyBody : public Body<dimension> {
   using typename Body<dimension>::ScalarFunction;
   using typename Body<dimension>::VectorField;
 
+  using MyConstantFunction = ConstantFunction<
+      Dune::FieldVector<double, dimension>, Dune::FieldVector<double, 1>>;
+
 public:
   MyBody(Dune::ParameterTree const &parset, MyGeometry const &mygeometry)
       : poissonRatio_(parset.get<double>("body.poissonRatio")),
         youngModulus_(parset.get<double>("body.youngModulus")),
         shearViscosity_(parset.get<double>("body.shearViscosity")),
         bulkViscosity_(parset.get<double>("body.bulkViscosity")),
-        density_(parset.get<double>("body.density")),
-        gravity_(parset.get<double>("gravity")) {}
+        densityField_(parset.get<double>("body.density")),
+        gravityField_(densityField_, mygeometry.zenith,
+                      parset.get<double>("gravity")) {}
 
   double getPoissonRatio() const override { return poissonRatio_; }
 
@@ -27,16 +34,18 @@ template <int dimension> class MyBody : public Body<dimension> {
   double getShearViscosity() const override { return shearViscosity_; }
   double getBulkViscosity() const override { return bulkViscosity_; }
 
-  double getDensity() const override { return density_; }
-  double getGravity() const override { return gravity_; }
+  ScalarFunction const &getDensityField() const override {
+    return densityField_;
+  }
+  VectorField const &getGravityField() const override { return gravityField_; }
 
 private:
   double const poissonRatio_;
   double const youngModulus_;
   double const shearViscosity_;
   double const bulkViscosity_;
-  double const density_;
-  double const gravity_;
+  MyConstantFunction const densityField_;
+  Gravity<dimension> const gravityField_;
 };
 
 #endif

src/one-body-sample.cc

@@ -165,7 +165,7 @@ int main(int argc, char *argv[]) {
                                    body.getPoissonRatio(), A);
     myAssembler.assembleViscosity(body.getShearViscosity(),
                                   body.getBulkViscosity(), C);
-    myAssembler.assembleMass(body.getDensity(), M);
+    myAssembler.assembleMass(body.getDensityField(), M);
     EnergyNorm<Matrix, Vector> const ANorm(A), MNorm(M);
     // Q: Does it make sense to weigh them in this manner?
     SumNorm<Vector> const AMNorm(1.0, ANorm, 1.0, MNorm);
@@ -178,8 +178,7 @@ int main(int argc, char *argv[]) {
 
     // Assemble forces
     Vector gravityFunctional;
-    myAssembler.assembleBodyForce(body.getGravity(), body.getDensity(),
-                                  myGeometry.zenith, gravityFunctional);
+    myAssembler.assembleBodyForce(body.getGravityField(), gravityFunctional);
 
     // Problem formulation: right-hand side
     auto const computeExternalForces = [&](double _relativeTime, Vector &_ell) {
@@ -204,7 +203,8 @@ int main(int argc, char *argv[]) {
       // volume * gravity * density / area    = normal stress
       // V      * g       * rho     / A       = sigma_n
       // m^d    * N/kg    * kg/m^d  / m^(d-1) = N/m^(d-1)
-      normalStress = -volume * body.getGravity() * body.getDensity() / area;
+      normalStress = -volume * parset.get<double>("gravity") *
+                     parset.get<double>("body.density") / area;
     }
     FrictionData const frictionData(parset.sub("boundary.friction"),
                                     normalStress);