diff --git a/src/assemblers.cc b/src/assemblers.cc
index 029c1767e39dea37bdf240eceff0ec5dfd53548b..160169fda8fc9fb4f676b87e64cf1137d8de3514 100644
--- a/src/assemblers.cc
+++ b/src/assemblers.cc
@@ -2,56 +2,113 @@
#include "config.h"
#endif
+#include <dune/istl/scaledidmatrix.hh>
+
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/functions/constantfunction.hh>
+#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/tectonic/globalruinanonlinearity.hh>
#include "assemblers.hh"
-// Assembles Neumann boundary term in f
-template <class GridView, class LocalVector, class Assembler>
-void assembleNeumann(GridView const &gridView, Assembler const &assembler,
- Dune::BitSetVector<1> const &neumannNodes,
- Dune::BlockVector<LocalVector> &f,
- Dune::VirtualFunction<double, double> const &neumann,
- double relativeTime) { // constant sample function on
- // neumann boundary
- BoundaryPatch<GridView> const neumannBoundary(gridView, neumannNodes);
+template <class GridView, int dimension>
+MyAssembler<GridView, dimension>::MyAssembler(GridView const &_gridView)
+ : vertexBasis(_gridView),
+ gridView(_gridView),
+ assembler(vertexBasis, vertexBasis) {}
+
+template <class GridView, int dimension>
+void MyAssembler<GridView, dimension>::assembleFrictionalBoundaryMass(
+ BoundaryPatch<GridView> const &frictionalBoundary,
+ ScalarMatrix &frictionalBoundaryMass) {
+ BoundaryMassAssembler<Grid, BoundaryPatch<GridView>, LocalVertexBasis,
+ LocalVertexBasis, Dune::FieldMatrix<double, 1, 1>> const
+ frictionalBoundaryMassAssembler(frictionalBoundary);
+ assembler.assembleOperator(frictionalBoundaryMassAssembler,
+ frictionalBoundaryMass);
+}
+
+template <class GridView, int dimension>
+void MyAssembler<GridView, dimension>::assembleMass(double density, Matrix &M) {
+ MassAssembler<Grid, LocalVertexBasis, LocalVertexBasis,
+ Dune::ScaledIdentityMatrix<double, dimension>> const localMass;
+ assembler.assembleOperator(localMass, M);
+ M *= density;
+}
+
+template <class GridView, int dimension>
+void MyAssembler<GridView, dimension>::assembleElasticity(double E, double nu,
+ Matrix &A) {
+ StVenantKirchhoffAssembler<Grid, LocalVertexBasis, LocalVertexBasis> const
+ localStiffness(E, nu);
+ assembler.assembleOperator(localStiffness, A);
+}
+
+template <class GridView, int dimension>
+void MyAssembler<GridView, dimension>::assembleViscosity(double shearViscosity,
+ double bulkViscosity,
+ Matrix &C) {
+ ViscosityAssembler<Grid, LocalVertexBasis, LocalVertexBasis> const
+ localViscosity(shearViscosity, bulkViscosity);
+ assembler.assembleOperator(localViscosity, C);
+}
+
+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);
+ L2FunctionalAssembler<Grid, LocalVertexBasis, LocalVector>
+ gravityFunctionalAssembler(gravityFunction);
+ assembler.assembleFunctional(gravityFunctionalAssembler, f);
+}
+
+template <class GridView, int dimension>
+void MyAssembler<GridView, dimension>::assembleNeumann(
+ BoundaryPatch<GridView> const &neumannBoundary, Vector &f,
+ Dune::VirtualFunction<double, double> const &neumann, double relativeTime) {
LocalVector localNeumann(0);
neumann.evaluate(relativeTime, localNeumann[0]);
ConstantFunction<LocalVector, LocalVector> const fNeumann(localNeumann);
- NeumannBoundaryAssembler<typename GridView::Grid, LocalVector>
- neumannBoundaryAssembler(fNeumann);
+ NeumannBoundaryAssembler<Grid, LocalVector> neumannBoundaryAssembler(
+ fNeumann);
assembler.assembleBoundaryFunctional(neumannBoundaryAssembler, f,
neumannBoundary);
}
-// Assembles constant 1-function on frictional boundary in nodalIntegrals
-template <class GridView, class LocalVector, class Assembler>
-Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
-assembleFrictionWeightsal(GridView const &gridView, Assembler const &assembler,
- Dune::BitSetVector<1> const &frictionalNodes) {
- using Scalar = Dune::FieldVector<double, 1>;
- BoundaryPatch<GridView> const frictionalBoundary(gridView, frictionalNodes);
- ConstantFunction<LocalVector, Scalar> const constantOneFunction(1);
- NeumannBoundaryAssembler<typename GridView::Grid, Scalar>
- frictionalBoundaryAssembler(constantOneFunction);
-
- auto const nodalIntegrals = Dune::make_shared<Dune::BlockVector<Scalar>>();
- assembler.assembleBoundaryFunctional(frictionalBoundaryAssembler,
- *nodalIntegrals, frictionalBoundary);
- return nodalIntegrals;
-}
+template <class GridView, int dimension>
+auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
+ BoundaryPatch<GridView> const &frictionalBoundary, FrictionData const &fd)
+ -> std::shared_ptr<Dune::GlobalNonlinearity<Matrix, Vector>> {
+ // Lump negative normal stress (kludge)
+ ScalarVector weights;
+ {
+ ConstantFunction<LocalVector, typename ScalarVector::block_type> const
+ constantOneFunction(1);
+ NeumannBoundaryAssembler<Grid, typename ScalarVector::block_type>
+ frictionalBoundaryAssembler(constantOneFunction);
+ assembler.assembleBoundaryFunctional(frictionalBoundaryAssembler, weights,
+ frictionalBoundary);
+ }
+ for (size_t i = 0; i < weights.size(); ++i)
+ assert(weights[i] >= 0.0);
-template <class Matrix, class Vector>
-Dune::shared_ptr<Dune::GlobalNonlinearity<Matrix, Vector>> assembleNonlinearity(
- Dune::BitSetVector<1> const &frictionalNodes,
- Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
- FrictionData const &fd) {
- return Dune::make_shared<Dune::GlobalRuinaNonlinearity<Matrix, Vector>>(
- frictionalNodes, nodalIntegrals, fd);
+ Dune::BitSetVector<1> frictionalNodes;
+ frictionalBoundary.getVertices(frictionalNodes);
+ return std::make_shared<Dune::GlobalRuinaNonlinearity<Matrix, Vector>>(
+ frictionalNodes, weights, fd);
}
#include "assemblers_tmpl.cc"
diff --git a/src/assemblers.hh b/src/assemblers.hh
index 75d3dd487b5341aae740cb4e952144df9663bb92..a7e0721fac2601e2b882f3825858b0d872704940 100644
--- a/src/assemblers.hh
+++ b/src/assemblers.hh
@@ -3,29 +3,59 @@
#include <dune/common/bitsetvector.hh>
#include <dune/common/function.hh>
-#include <dune/common/fvector.hh>
-#include <dune/common/shared_ptr.hh>
+#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/fufem/assemblers/assembler.hh>
+#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
#include <dune/tectonic/globalnonlinearity.hh>
-template <class GridView, class LocalVector, class Assembler>
-void assembleNeumann(GridView const &gridView, Assembler const &assembler,
- Dune::BitSetVector<1> const &neumannNodes,
- Dune::BlockVector<LocalVector> &f,
- Dune::VirtualFunction<double, double> const &neumann,
- double relativeTime);
-
-template <class GridView, class LocalVector, class Assembler>
-Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
-assembleFrictionWeightsal(GridView const &gridView, Assembler const &assembler,
- Dune::BitSetVector<1> const &frictionalNodes);
-
-template <class Matrix, class Vector>
-Dune::shared_ptr<Dune::GlobalNonlinearity<Matrix, Vector>> assembleNonlinearity(
- Dune::BitSetVector<1> const &frictionalNodes,
- Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
- FrictionData const &fd);
+template <class GridView, int dimension> class MyAssembler {
+public:
+ using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
+ using Matrix =
+ Dune::BCRSMatrix<Dune::FieldMatrix<double, dimension, dimension>>;
+ using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
+ using Vector = Dune::BlockVector<Dune::FieldVector<double, dimension>>;
+
+ using VertexBasis = P1NodalBasis<GridView, double>;
+ VertexBasis vertexBasis;
+
+private:
+ using Grid = typename GridView::Grid;
+ using LocalVector = typename Vector::block_type;
+ using LocalVertexBasis = typename VertexBasis::LocalFiniteElement;
+
+ GridView const &gridView;
+ Assembler<VertexBasis, VertexBasis> assembler;
+
+public:
+ MyAssembler(GridView const &gridView);
+
+ void assembleFrictionalBoundaryMass(
+ BoundaryPatch<GridView> const &frictionalBoundary,
+ ScalarMatrix &frictionalBoundaryMass);
+
+ void assembleMass(double density, 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 assembleNeumann(BoundaryPatch<GridView> const &neumannBoundary,
+ Vector &f,
+ Dune::VirtualFunction<double, double> const &neumann,
+ double relativeTime);
+
+ std::shared_ptr<Dune::GlobalNonlinearity<Matrix, Vector>>
+ assembleFrictionNonlinearity(
+ BoundaryPatch<GridView> const &frictionalBoundary,
+ FrictionData const &fd);
+};
+
#endif
diff --git a/src/assemblers_tmpl.cc b/src/assemblers_tmpl.cc
index 35cddb4a47c413fe367237e3c20c856394ffd2a2..36628ed9d1a317fe9c90dd089f081c4912d1e69e 100644
--- a/src/assemblers_tmpl.cc
+++ b/src/assemblers_tmpl.cc
@@ -2,27 +2,6 @@
#error DIM unset
#endif
-#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
-
#include "explicitgrid.hh"
-#include "explicitvectors.hh"
-
-using P1Basis = P1NodalBasis<GridView, double>;
-using MyAssembler = Assembler<P1Basis, P1Basis>;
-
-template void assembleNeumann<GridView, SmallVector, MyAssembler>(
- GridView const &gridView, MyAssembler const &assembler,
- Dune::BitSetVector<1> const &neumannNodes,
- Dune::BlockVector<SmallVector> &f,
- Dune::VirtualFunction<double, double> const &neumann, double relativeTime);
-
-template Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
-assembleFrictionWeightsal<GridView, SmallVector, MyAssembler>(
- GridView const &gridView, MyAssembler const &assembler,
- Dune::BitSetVector<1> const &frictionalNodes);
-template Dune::shared_ptr<Dune::GlobalNonlinearity<Matrix, Vector>>
-assembleNonlinearity<Matrix, Vector>(
- Dune::BitSetVector<1> const &frictionalNodes,
- Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
- FrictionData const &fd);
+template class MyAssembler<GridView, DIM>;
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 8d96bc6849043e13c7dc3dcd0840c38b5a8d2472..2d25e7d1515471662dcbc0a5ccfb0eadd2d9db73 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -49,19 +49,12 @@
#include <dune/grid/utility/structuredgridfactory.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
-#include <dune/istl/scaledidmatrix.hh>
#include <dune/fufem/assemblers/assembler.hh>
-#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/stvenantkirchhoffassembler.hh>
-#include <dune/fufem/assemblers/localassemblers/viscosityassembler.hh>
#include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/dunepython.hh>
#include <dune/fufem/functions/basisgridfunction.hh>
-#include <dune/fufem/functions/constantfunction.hh>
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wsign-compare"
@@ -163,6 +156,8 @@ int main(int argc, char *argv[]) {
auto const E = parset.get<double>("body.E");
auto const nu = parset.get<double>("body.nu");
+ auto const density = parset.get<double>("body.density");
+ double const gravity = 9.81;
// {{{ Set up grid
using Grid = Dune::ALUGrid<dims, dims, Dune::simplex, Dune::nonconforming>;
@@ -182,6 +177,8 @@ int main(int argc, char *argv[]) {
grid = Dune::StructuredGridFactory<Grid>::createSimplexGrid(
lowerLeft, upperRight, elements);
}
+ SmallVector zenith(0);
+ zenith[1] = 1;
auto const refinements = parset.get<size_t>("grid.refinements");
grid->globalRefine(refinements);
@@ -193,11 +190,8 @@ int main(int argc, char *argv[]) {
// Set up bases
using P0Basis = P0Basis<GridView, double>;
- using P1Basis = P1NodalBasis<GridView, double>;
P0Basis const p0Basis(leafView);
- P1Basis const p1Basis(leafView);
Assembler<P0Basis, P0Basis> const p0Assembler(p0Basis, p0Basis);
- Assembler<P1Basis, P1Basis> const p1Assembler(p1Basis, p1Basis);
// Set up the boundary
Dune::BitSetVector<dims> velocityDirichletNodes(fineVertexCount, false);
@@ -205,9 +199,10 @@ int main(int argc, char *argv[]) {
velocityDirichletNodes;
Dune::BitSetVector<dims> accelerationDirichletNodes(fineVertexCount, false);
Dune::BitSetVector<1> neumannNodes(fineVertexCount, false);
- Dune::BitSetVector<1> frictionalNodes(fineVertexCount, false);
+ BoundaryPatch<GridView> const neumannBoundary(leafView, neumannNodes);
Vector vertexCoordinates(fineVertexCount);
+ Dune::BitSetVector<1> frictionalNodes(fineVertexCount, false);
{
Dune::MultipleCodimMultipleGeomTypeMapper<
GridView, Dune::MCMGVertexLayout> const vertexMapper(leafView);
@@ -254,89 +249,54 @@ int main(int argc, char *argv[]) {
// Set up normal stress, mass matrix, and gravity functional
double normalStress;
- Matrix M;
- EnergyNorm<Matrix, Vector> const MNorm(M);
- Vector gravityFunctional;
{
- double const gravity = 9.81;
- double const density = parset.get<double>("body.density");
- {
- MassAssembler<Grid, P1Basis::LocalFiniteElement,
- P1Basis::LocalFiniteElement,
- Dune::ScaledIdentityMatrix<double, dims>> const localMass;
- p1Assembler.assembleOperator(localMass, M);
- M *= density;
- }
- {
- double volume = 1.0;
- for (size_t i = 0; i < dims; ++i)
- volume *= (upperRight[i] - lowerLeft[i]);
-
- double area = 1.0;
- for (size_t i = 0; i < dims; ++i)
- if (i != 1)
- area *= (upperRight[i] - lowerLeft[i]);
-
- // 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 * gravity * density / area;
- }
- {
- SmallVector weightedGravitationalDirection(0);
- weightedGravitationalDirection[1] = -density * gravity;
- ConstantFunction<SmallVector, SmallVector> const gravityFunction(
- weightedGravitationalDirection);
- L2FunctionalAssembler<Grid, P1Basis::LocalFiniteElement, SmallVector>
- gravityFunctionalAssembler(gravityFunction);
- p1Assembler.assembleFunctional(gravityFunctionalAssembler,
- gravityFunctional);
- }
+ double volume = 1.0;
+ for (size_t i = 0; i < dims; ++i)
+ volume *= (upperRight[i] - lowerLeft[i]);
+
+ double area = 1.0;
+ for (size_t i = 0; i < dims; ++i)
+ if (i != 1)
+ area *= (upperRight[i] - lowerLeft[i]);
+
+ // 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 * gravity * density / area;
}
FrictionData const frictionData(parset.sub("boundary.friction"),
normalStress);
- Matrix A;
- EnergyNorm<Matrix, Vector> const ANorm(A);
- {
- StVenantKirchhoffAssembler<Grid, P1Basis::LocalFiniteElement,
- P1Basis::LocalFiniteElement> const
- localStiffness(E, nu);
- p1Assembler.assembleOperator(localStiffness, A);
- }
- Matrix C;
- {
- ViscosityAssembler<Grid, P1Basis::LocalFiniteElement,
- P1Basis::LocalFiniteElement> const
- localViscosity(parset.get<double>("body.shearViscosity"),
- parset.get<double>("body.bulkViscosity"));
- p1Assembler.assembleOperator(localViscosity, C);
- }
+ using MyAssembler = MyAssembler<GridView, dims>;
- ScalarMatrix frictionalBoundaryMassMatrix;
- EnergyNorm<ScalarMatrix, ScalarVector> const stateEnergyNorm(
- frictionalBoundaryMassMatrix);
- {
- BoundaryMassAssembler<
- Grid, BoundaryPatch<GridView>, P1Basis::LocalFiniteElement,
- P1Basis::LocalFiniteElement, SmallScalarMatrix> const
- frictionalBoundaryMassAssembler(frictionalBoundary);
- p1Assembler.assembleOperator(frictionalBoundaryMassAssembler,
- frictionalBoundaryMassMatrix);
- }
+ MyAssembler myAssembler(leafView);
+
+ Matrix A, C, M;
+ myAssembler.assembleElasticity(E, nu, A);
+ myAssembler.assembleViscosity(parset.get<double>("body.shearViscosity"),
+ parset.get<double>("body.bulkViscosity"), C);
+ myAssembler.assembleMass(density, 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);
- auto const nodalIntegrals =
- assembleFrictionWeightsal<GridView, SmallVector>(leafView, p1Assembler,
- frictionalNodes);
- auto myGlobalNonlinearity = assembleNonlinearity<Matrix, Vector>(
- frictionalNodes, *nodalIntegrals, frictionData);
+ ScalarMatrix frictionalBoundaryMass;
+ myAssembler.assembleFrictionalBoundaryMass(frictionalBoundary,
+ frictionalBoundaryMass);
+ EnergyNorm<ScalarMatrix, ScalarVector> const stateEnergyNorm(
+ frictionalBoundaryMass);
+
+ // Assemble forces
+ Vector gravityFunctional;
+ myAssembler.assembleBodyForce(gravity, density, zenith, gravityFunctional);
+
+ auto myGlobalNonlinearity = myAssembler.assembleFrictionNonlinearity(
+ frictionalBoundary, frictionData);
// Problem formulation: right-hand side
auto const createRHS = [&](double _relativeTime, Vector &_ell) {
- assembleNeumann(leafView, p1Assembler, neumannNodes, _ell,
- neumannFunction, _relativeTime);
+ myAssembler.assembleNeumann(neumannBoundary, _ell, neumannFunction,
+ _relativeTime);
_ell += gravityFunctional;
};
Vector ell(fineVertexCount);
@@ -569,15 +529,15 @@ int main(int argc, char *argv[]) {
relaxationWriter << std::endl;
if (parset.get<bool>("io.writeVTK")) {
- auto const gridDisplacement =
- Dune::make_shared<BasisGridFunction<P1Basis, Vector> const>(p1Basis,
- u);
+ auto const gridDisplacement = Dune::make_shared<
+ BasisGridFunction<typename MyAssembler::VertexBasis, Vector> const>(
+ myAssembler.vertexBasis, u);
ScalarVector vonMisesStress;
VonMisesStressAssembler<Grid, P0Basis::LocalFiniteElement>
localStressAssembler(E, nu, gridDisplacement);
p0Assembler.assembleFunctional(localStressAssembler, vonMisesStress);
- writeVtk(p1Basis, u, alpha, p0Basis, vonMisesStress,
+ writeVtk(myAssembler.vertexBasis, u, alpha, p0Basis, vonMisesStress,
(boost::format("obs%d") % run).str());
}
}