[Cleanup] Outsource more assembly; Turn header into class

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"

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

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>;

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());
       }
     }