[Cleanup] Allow for variable viscosity

dune/tectonic/body.hh

@@ -12,8 +12,8 @@ template <int dimension> struct Body {
 
   double virtual getYoungModulus() const = 0;
 
-  double virtual getShearViscosity() const = 0;
-  double virtual getBulkViscosity() const = 0;
+  ScalarFunction virtual const &getShearViscosityField() const = 0;
+  ScalarFunction virtual const &getBulkViscosityField() const = 0;
 
   ScalarFunction virtual const &getDensityField() const = 0;
   VectorField virtual const &getGravityField() const = 0;

src/assemblers.cc

@@ -8,13 +8,14 @@
 #include <dune/fufem/assemblers/localassemblers/l2functionalassembler.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/variablecoefficientviscosityassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/weightedmassassembler.hh>
 #include <dune/fufem/boundarypatch.hh>
 #include <dune/fufem/computestress.hh>
 #include <dune/fufem/functions/basisgridfunction.hh>
 #include <dune/fufem/functions/constantfunction.hh>
+#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
 
 #include "assemblers.hh"
 
@@ -61,11 +62,18 @@ void MyAssembler<GridView, dimension>::assembleElasticity(double E, double nu,
 }
 
 template <class GridView, int dimension>
-void MyAssembler<GridView, dimension>::assembleViscosity(double shearViscosity,
-                                                         double bulkViscosity,
-                                                         Matrix &C) {
-  ViscosityAssembler<Grid, LocalVertexBasis, LocalVertexBasis> const
-  localViscosity(shearViscosity, bulkViscosity);
+void MyAssembler<GridView, dimension>::assembleViscosity(
+    Dune::VirtualFunction<LocalVector, LocalScalarVector> const &shearViscosity,
+    Dune::VirtualFunction<LocalVector, LocalScalarVector> const &bulkViscosity,
+    Matrix &C) {
+  // NOTE: We treat the weights as constant functions
+  QuadratureRuleKey shearViscosityKey(dimension, 0);
+  QuadratureRuleKey bulkViscosityKey(dimension, 0);
+  VariableCoefficientViscosityAssembler<
+      Grid, LocalVertexBasis, LocalVertexBasis,
+      Dune::VirtualFunction<LocalVector, LocalScalarVector>> const
+  localViscosity(gridView.grid(), shearViscosity, bulkViscosity,
+                 shearViscosityKey, bulkViscosityKey);
   vertexAssembler.assembleOperator(localViscosity, C);
 }
 

src/assemblers.hh

@@ -55,8 +55,12 @@ template <class GridView, int dimension> class MyAssembler {
 
   void assembleElasticity(double E, double nu, Matrix &A);
 
-  void assembleViscosity(double shearViscosity, double bulkViscosity,
-                         Matrix &C);
+  void assembleViscosity(
+      Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
+          shearViscosity,
+      Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
+          bulkViscosity,
+      Matrix &C);
 
   void assembleBodyForce(
       Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,

src/mybody.hh

@@ -20,8 +20,8 @@ template <int dimension> class MyBody : public Body<dimension> {
   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")),
+        shearViscosityField_(parset.get<double>("body.shearViscosity")),
+        bulkViscosityField_(parset.get<double>("body.bulkViscosity")),
         densityField_(parset.get<double>("body.density")),
         gravityField_(densityField_, mygeometry.zenith,
                       parset.get<double>("gravity")) {}
@@ -30,8 +30,12 @@ template <int dimension> class MyBody : public Body<dimension> {
 
   double getYoungModulus() const override { return youngModulus_; }
 
-  double getShearViscosity() const override { return shearViscosity_; }
-  double getBulkViscosity() const override { return bulkViscosity_; }
+  ScalarFunction const &getShearViscosityField() const override {
+    return shearViscosityField_;
+  }
+  ScalarFunction const &getBulkViscosityField() const override {
+    return bulkViscosityField_;
+  }
 
   ScalarFunction const &getDensityField() const override {
     return densityField_;
@@ -41,8 +45,8 @@ template <int dimension> class MyBody : public Body<dimension> {
 private:
   double const poissonRatio_;
   double const youngModulus_;
-  double const shearViscosity_;
-  double const bulkViscosity_;
+  MyConstantFunction const shearViscosityField_;
+  MyConstantFunction const bulkViscosityField_;
   MyConstantFunction const densityField_;
   Gravity<dimension> const gravityField_;
 };

src/one-body-sample.cc

@@ -163,8 +163,8 @@ int main(int argc, char *argv[]) {
     Matrix A, C, M;
     myAssembler.assembleElasticity(body.getYoungModulus(),
                                    body.getPoissonRatio(), A);
-    myAssembler.assembleViscosity(body.getShearViscosity(),
-                                  body.getBulkViscosity(), C);
+    myAssembler.assembleViscosity(body.getShearViscosityField(),
+                                  body.getBulkViscosityField(), C);
     myAssembler.assembleMass(body.getDensityField(), M);
     EnergyNorm<Matrix, Vector> const ANorm(A), MNorm(M);
     // Q: Does it make sense to weigh them in this manner?