diff --git a/dune/tectonic/body.hh b/dune/tectonic/body.hh
index cadcbc52974fb274146ab9d000ef72bebc28d370..3365eac427f962180616221d4592e0bcdb0b854c 100644
--- a/dune/tectonic/body.hh
+++ b/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
diff --git a/dune/tectonic/gravity.hh b/dune/tectonic/gravity.hh
new file mode 100644
index 0000000000000000000000000000000000000000..9af62f17d86248d6deeddec1f4b0eb7f3cf65afc
--- /dev/null
+++ b/dune/tectonic/gravity.hh
@@ -0,0 +1,34 @@
+#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
diff --git a/src/assemblers.cc b/src/assemblers.cc
index 80f17a476497cdc9787d177df914a05b61e1d1df..83565cce966b7bfb67de5dcb6751bd2eef2edca8 100644
--- a/src/assemblers.cc
+++ b/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);
}
diff --git a/src/assemblers.hh b/src/assemblers.hh
index 0987f2a18625e04e8b83768084f8ba746047b443..a00681581382d72f95db4a5de6e45d6c0838600b 100644
--- a/src/assemblers.hh
+++ b/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,
diff --git a/src/mybody.hh b/src/mybody.hh
index 47d3d26c3d189b72320194f170d8041aaf126b25..ec735642edbca78cf2e0d63f02421c2f4131da77 100644
--- a/src/mybody.hh
+++ b/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
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 52b5573b1649f9c85d41027b8e4f93e23b998c7d..38659cfd41ca6f284911896470355b8a02c3182e 100644
--- a/src/one-body-sample.cc
+++ b/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);