src/spatial-solving/fixedpointiterator.cc

@@ -4,10 +4,24 @@
 
 #include <dune/common/exceptions.hh>
 
-#include <dune/solvers/common/arithmetic.hh>
+#include <dune/matrix-vector/axpy.hh>
+
 #include <dune/solvers/norms/energynorm.hh>
 #include <dune/solvers/solvers/loopsolver.hh>
 
+#include <dune/contact/assemblers/nbodyassembler.hh>
+#include <dune/contact/common/dualbasisadapter.hh>
+
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+
+#include <dune/functions/gridfunctions/gridfunction.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+#include <dune/geometry/type.hh>
+#include <dune/geometry/referenceelements.hh>
+
+#include <dune/fufem/functions/basisgridfunction.hh>
+
 #include "../enums.hh"
 #include "../enumparser.hh"
 
@@ -22,8 +36,9 @@ void FixedPointIterationCounter::operator+=(
 template <class Factory, class Updaters, class ErrorNorm>
 FixedPointIterator<Factory, Updaters, ErrorNorm>::FixedPointIterator(
     Factory &factory, Dune::ParameterTree const &parset,
-    std::shared_ptr<Nonlinearity> globalFriction, ErrorNorm const &errorNorm)
-    : step_(factory.getStep()),
+    std::vector<std::shared_ptr<Nonlinearity>>& globalFriction, const std::vector<const ErrorNorm* >& errorNorms)
+    : factory_(factory),
+      step_(factory_.getStep()),
       parset_(parset),
       globalFriction_(globalFriction),
       fixedPointMaxIterations_(parset.get<size_t>("v.fpi.maximumIterations")),
@@ -32,47 +47,100 @@ FixedPointIterator<Factory, Updaters, ErrorNorm>::FixedPointIterator(
       velocityMaxIterations_(parset.get<size_t>("v.solver.maximumIterations")),
       velocityTolerance_(parset.get<double>("v.solver.tolerance")),
       verbosity_(parset.get<Solver::VerbosityMode>("v.solver.verbosity")),
-      errorNorm_(errorNorm) {}
+      errorNorms_(errorNorms) {}
 
 template <class Factory, class Updaters, class ErrorNorm>
 FixedPointIterationCounter
 FixedPointIterator<Factory, Updaters, ErrorNorm>::run(
-    Updaters updaters, Matrix const &velocityMatrix, Vector const &velocityRHS,
-    Vector &velocityIterate) {
-  EnergyNorm<Matrix, Vector> energyNorm(velocityMatrix);
+    Updaters updaters, const std::vector<Matrix>& velocityMatrices, const std::vector<Vector>& velocityRHSs,
+    std::vector<Vector>& velocityIterates) {
+
+  auto contactAssembler = factory_.getNBodyAssembler();
+
+  Matrix mat;
+  std::vector<const Matrix*> matrices_ptr(velocityMatrices.size());
+  for (size_t i=0; i<matrices_ptr.size(); i++) {
+      matrices_ptr[i] = &velocityMatrices[i];
+  }
+  contactAssembler.assembleJacobian(matrices_ptr, mat);
+
+  EnergyNorm<Matrix, Vector> energyNorm(mat);
   LoopSolver<Vector> velocityProblemSolver(step_.get(), velocityMaxIterations_,
                                            velocityTolerance_, &energyNorm,
                                            verbosity_, false); // absolute error
 
   size_t fixedPointIteration;
   size_t multigridIterations = 0;
-  ScalarVector alpha;
+  std::vector<ScalarVector> alpha;
   updaters.state_->extractAlpha(alpha);
   for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations_;
        ++fixedPointIteration) {
+
+    // compute relative velocities
+    std::vector<Vector> v_rel;
+    relativeVelocities(velocityIterates, v_rel);
+
+    // contribution from nonlinearity
+    for (size_t i=0; i<alpha.size(); i++) {
+      globalFriction_[i]->updateAlpha(alpha[i]);
+    }
+
+    std::vector<Matrix> matrices(velocityMatrices.size());
+    std::vector<Vector> rhs(velocityRHSs.size());
+    for (size_t i=0; i<globalFriction_.size(); i++) {
+      matrices[i] = velocityMatrices[i];
+      rhs[i] = velocityRHSs[i];
+
+      globalFriction_[i]->addHessian(v_rel[i], matrices[i]);
+      globalFriction_[i]->addGradient(v_rel[i], rhs[i]);
+
+      matrices_ptr[i] = &matrices[i];
+    }
+
+    // assemble full global contact problem
+    Matrix bilinearForm;
+    contactAssembler.assembleJacobian(matrices_ptr, bilinearForm);
+
+    Vector totalRhs;
+    contactAssembler.assembleRightHandSide(rhs, totalRhs);
+
+    Vector totalVelocityIterate;
+    contactAssembler.nodalToTransformed(velocityIterates, totalVelocityIterate);
+
     // solve a velocity problem
-    globalFriction_->updateAlpha(alpha);
-    ConvexProblem convexProblem(1.0, velocityMatrix, *globalFriction_,
-                                velocityRHS, velocityIterate);
-    BlockProblem velocityProblem(parset_, convexProblem);
-    step_->setProblem(velocityIterate, velocityProblem);
+    step_->setProblem(bilinearForm, totalVelocityIterate, totalRhs);
+
     velocityProblemSolver.preprocess();
     velocityProblemSolver.solve();
+    contactAssembler.postprocess(totalVelocityIterate, velocityIterates);
 
     multigridIterations += velocityProblemSolver.getResult().iterations;
 
-    Vector v_m;
+    std::vector<Vector> v_m;
     updaters.rate_->extractOldVelocity(v_m);
-    v_m *= 1.0 - lambda_;
-    Arithmetic::addProduct(v_m, lambda_, velocityIterate);
+
+    for (size_t i=0; i<v_m.size(); i++) {
+      v_m[i] *= 1.0 - lambda_;
+      Dune::MatrixVector::addProduct(v_m[i], lambda_, velocityIterates[i]);
+    }
+
+    // compute relative velocities on contact boundaries
+    relativeVelocities(v_m, v_rel);
 
     // solve a state problem
-    updaters.state_->solve(v_m);
-    ScalarVector newAlpha;
+    updaters.state_->solve(v_rel);
+    std::vector<ScalarVector> newAlpha;
     updaters.state_->extractAlpha(newAlpha);
 
-    if (lambda_ < 1e-12 or
-        errorNorm_.diff(alpha, newAlpha) < fixedPointTolerance_) {
+    bool breakCriterion = true;
+    for (size_t i=0; i<alpha.size(); i++) {
+        if (errorNorms_[i]->diff(alpha[i], newAlpha[i]) >= fixedPointTolerance_) {
+            breakCriterion = false;
+            break;
+        }
+    }
+
+    if (lambda_ < 1e-12 or breakCriterion) {
       fixedPointIteration++;
       break;
     }
@@ -81,7 +149,7 @@ FixedPointIterator<Factory, Updaters, ErrorNorm>::run(
   if (fixedPointIteration == fixedPointMaxIterations_)
     DUNE_THROW(Dune::Exception, "FPI failed to converge");
 
-  updaters.rate_->postProcess(velocityIterate);
+  updaters.rate_->postProcess(velocityIterates);
 
   // Cannot use return { fixedPointIteration, multigridIterations };
   // with gcc 4.9.2, see also http://stackoverflow.com/a/37777814/179927
@@ -97,4 +165,190 @@ std::ostream &operator<<(std::ostream &stream,
                 << ")";
 }
 
+template <class Factory, class Updaters, class ErrorNorm>
+void FixedPointIterator<Factory, Updaters, ErrorNorm>::relativeVelocities(const std::vector<Vector>& v, std::vector<Vector>& v_rel) const {
+  // init result
+  v_rel.resize(v.size());
+  for (size_t i=0; i<v_rel.size(); i++) {
+      v_rel[i].resize(v[i].size());
+      v_rel[i] = 0;
+  }
+
+  // needs assemblers to obtain basis
+    /*
+  std::vector<std::shared_ptr<MyAssembler>> assemblers(bodyCount);
+
+  using field_type = typename Factory::Matrix::field_type;
+
+
+  // adaptation of DualMortarCoupling::setup()
+
+  const size_t dim = DeformedGrid::dimension;
+  typedef typename DeformedGrid::LeafGridView GridView;
+
+  //cache of local bases
+  typedef Dune::PQkLocalFiniteElementCache<typename DeformedGrid::ctype, field_type, dim,1> FiniteElementCache1;
+  FiniteElementCache1 cache1;
+
+  // cache for the dual functions on the boundary
+  using DualCache = Dune::Contact::DualBasisAdapter<GridView, field_type>;
+  std::unique_ptr<DualCache> dualCache;
+  dualCache = std::make_unique< Dune::Contact::DualBasisAdapterGlobal<GridView, field_type> >();
+
+  // define FE grid functions
+  std::vector<BasisGridFunction<VertexBasis, Vector>* > gridFunctions(v_m.size());
+  for (size_t i=0; i<gridFunctions.size(); i++) {
+    gridFunctions[i] = new BasisGridFunction<MyAssembler::VertexBasis, Vector>(assemblers[i]->vertexBasis, v_m[i]);
+  }
+ */
+
+  /*
+  for (size_t i=0; i<nBodyAssembler_.nCouplings(); i++) {
+    const auto& coupling = nBodyAssembler_.getCoupling(i);
+    auto glue = coupling.backend();
+
+    const std::array<int, 2> gridIdx = coupling.gridIdx_;
+    const int nonmortarGridIdx = ;
+    const int mortarGridIdx = ;
+
+    // loop over all intersections
+    for (const auto& rIs : intersections(glue)) {
+        const auto& inside = rIs.inside();
+
+        if (!nonmortarBoundary_.contains(rIs.inside(),rIs.indexInInside()))
+            continue;
+
+        const auto& outside = rIs.outside();
+
+        // types of the elements supporting the boundary segments in question
+        Dune::GeometryType nonmortarEType = inside.type();
+        Dune::GeometryType mortarEType    = outside.type();
+
+        const auto& domainRefElement = Dune::ReferenceElements<ctype, dim>::general(nonmortarEType);
+        const auto& targetRefElement = Dune::ReferenceElements<ctype, dim>::general(mortarEType);
+
+        int noOfMortarVec = targetRefElement.size(dim);
+
+        Dune::GeometryType nmFaceType = domainRefElement.type(rIs.indexInInside(),1);
+        Dune::GeometryType mFaceType  = targetRefElement.type(rIs.indexInOutside(),1);
+
+        // Select a quadrature rule
+        // 2 in 2d and for integration over triangles in 3d.  If one (or both) of the two faces involved
+        // are quadrilaterals, then the quad order has to be risen to 3 (4).
+        int quadOrder = 2 + (!nmFaceType.isSimplex()) + (!mFaceType.isSimplex());
+        const auto& quadRule = Dune::QuadratureRules<ctype, dim-1>::rule(rIs.type(), quadOrder);
+
+        const auto& mortarFiniteElement = cache1.get(mortarEType);
+        dualCache->bind(inside, rIs.indexInInside());
+
+        std::vector<Dune::FieldVector<field_type,1> > mortarQuadValues, dualQuadValues;
+
+        const auto& rGeom = rIs.geometry();
+        const auto& rGeomOutside = rIs.geometryOutside();
+        const auto& rGeomInInside = rIs.geometryInInside();
+        const auto& rGeomInOutside = rIs.geometryInOutside();
+
+        int nNonmortarFaceNodes = domainRefElement.size(rIs.indexInInside(),1,dim);
+        std::vector<int> nonmortarFaceNodes;
+        for (int i=0; i<nNonmortarFaceNodes; i++) {
+          int faceIdxi = domainRefElement.subEntity(rIs.indexInInside(), 1, i, dim);
+          nonmortarFaceNodes.push_back(faceIdxi);
+        }
+
+        for (const auto& quadPt : quadRule) {
+
+            // compute integration element of overlap
+            ctype integrationElement = rGeom.integrationElement(quadPt.position());
+
+            // quadrature point positions on the reference element
+            Dune::FieldVector<ctype,dim> nonmortarQuadPos = rGeomInInside.global(quadPt.position());
+            Dune::FieldVector<ctype,dim> mortarQuadPos    = rGeomInOutside.global(quadPt.position());
+
+            // The current quadrature point in world coordinates
+            Dune::FieldVector<field_type,dim> nonmortarQpWorld = rGeom.global(quadPt.position());
+            Dune::FieldVector<field_type,dim> mortarQpWorld    = rGeomOutside.global(quadPt.position());;
+
+            // the gap direction (normal * gapValue)
+            Dune::FieldVector<field_type,dim> gapVector = mortarQpWorld  - nonmortarQpWorld;
+
+            //evaluate all shapefunctions at the quadrature point
+            //nonmortarFiniteElement.localBasis().evaluateFunction(nonmortarQuadPos,nonmortarQuadValues);
+            mortarFiniteElement.localBasis().evaluateFunction(mortarQuadPos,mortarQuadValues);
+            dualCache->evaluateFunction(nonmortarQuadPos,dualQuadValues);
+
+            // loop over all Lagrange multiplier shape functions
+            for (int j=0; j<nNonmortarFaceNodes; j++) {
+
+                int globalDomainIdx = indexSet0.subIndex(inside,nonmortarFaceNodes[j],dim);
+                int rowIdx = globalToLocal[globalDomainIdx];
+
+                weakObstacle_[rowIdx][0] += integrationElement * quadPt.weight()
+                    * dualQuadValues[nonmortarFaceNodes[j]] * (gapVector*avNormals[globalDomainIdx]);
+
+                // loop over all mortar shape functions
+                for (int k=0; k<noOfMortarVec; k++) {
+
+                    int colIdx  = indexSet1.subIndex(outside, k, dim);
+                    if (!mortarBoundary_.containsVertex(colIdx))
+                        continue;
+
+                    // Integrate over the product of two shape functions
+                    field_type mortarEntry =  integrationElement* quadPt.weight()* dualQuadValues[nonmortarFaceNodes[j]]* mortarQuadValues[k];
+
+                    Dune::MatrixVector::addToDiagonal(mortarLagrangeMatrix_[rowIdx][colIdx], mortarEntry);
+
+                }
+
+            }
+
+        }
+
+
+
+
+
+
+  }
+
+
+      // Create mapping from the global set of block dofs to the ones on the contact boundary
+      std::vector<int> globalToLocal;
+      nonmortarBoundary_.makeGlobalToLocal(globalToLocal);
+
+      // loop over all intersections
+      for (const auto& rIs : intersections(glue)) {
+
+          if (!nonmortarBoundary_.contains(rIs.inside(),rIs.indexInInside()))
+              continue;
+
+          const auto& inside = rIs.inside();
+          const auto& outside = rIs.outside();
+
+          const auto& domainRefElement = Dune::ReferenceElements<ctype, dim>::general(inside.type());
+          const auto& targetRefElement = Dune::ReferenceElements<ctype, dim>::general(outside.type());
+
+          int nDomainVertices = domainRefElement.size(dim);
+          int nTargetVertices = targetRefElement.size(dim);
+
+          for (int j=0; j<nDomainVertices; j++) {
+
+              int localDomainIdx = globalToLocal[indexSet0.subIndex(inside,j,dim)];
+
+              // if the vertex is not contained in the restricted contact boundary then dismiss it
+              if (localDomainIdx == -1)
+                  continue;
+
+              for (int k=0; k<nTargetVertices; k++) {
+                  int globalTargetIdx = indexSet1.subIndex(outside,k,dim);
+                  if (!mortarBoundary_.containsVertex(globalTargetIdx))
+                      continue;
+
+                  mortarIndices.add(localDomainIdx, globalTargetIdx);
+              }
+          }
+      }
+    */
+
+}
+
 #include "fixedpointiterator_tmpl.cc"

src/spatial-solving/fixedpointiterator.hh

@@ -8,6 +8,8 @@
 #include <dune/solvers/norms/norm.hh>
 #include <dune/solvers/solvers/solver.hh>
 
+#include <dune/contact/assemblers/nbodyassembler.hh>
+
 struct FixedPointIterationCounter {
   size_t iterations = 0;
   size_t multigridIterations = 0;
@@ -23,24 +25,30 @@ class FixedPointIterator {
   using ScalarVector = typename Updaters::StateUpdater::ScalarVector;
   using Vector = typename Factory::Vector;
   using Matrix = typename Factory::Matrix;
-  using ConvexProblem = typename Factory::ConvexProblem;
-  using BlockProblem = typename Factory::BlockProblem;
-  using Nonlinearity = typename ConvexProblem::NonlinearityType;
+  using Nonlinearity = typename Factory::Nonlinearity;
+
+ // using Nonlinearity = typename ConvexProblem::NonlinearityType;
+
+  using DeformedGrid = typename Factory::DeformedGrid;
+
+private:
+  void relativeVelocities(const std::vector<Vector>& v, std::vector<Vector>& v_rel) const;
 
 public:
-  FixedPointIterator(Factory &factory, Dune::ParameterTree const &parset,
-                     std::shared_ptr<Nonlinearity> globalFriction,
-                     ErrorNorm const &errorNorm_);
+  FixedPointIterator(Factory& factory, const Dune::ParameterTree& parset,
+                     std::vector<std::shared_ptr<Nonlinearity>>& globalFriction,
+                     const std::vector<const ErrorNorm* >& errorNorms);
 
   FixedPointIterationCounter run(Updaters updaters,
-                                 Matrix const &velocityMatrix,
-                                 Vector const &velocityRHS,
-                                 Vector &velocityIterate);
+                                 const std::vector<Matrix>& velocityMatrices,
+                                 const std::vector<Vector>& velocityRHSs,
+                                 std::vector<Vector>& velocityIterates);
 
 private:
+  Factory& factory_;
   std::shared_ptr<typename Factory::Step> step_;
   Dune::ParameterTree const &parset_;
-  std::shared_ptr<Nonlinearity> globalFriction_;
+  std::vector<std::shared_ptr<Nonlinearity>>& globalFriction_;
 
   size_t fixedPointMaxIterations_;
   double fixedPointTolerance_;
@@ -48,6 +56,6 @@ class FixedPointIterator {
   size_t velocityMaxIterations_;
   double velocityTolerance_;
   Solver::VerbosityMode verbosity_;
-  ErrorNorm const &errorNorm_;
+  const std::vector<const ErrorNorm* >& errorNorms_;
 };
 #endif

src/spatial-solving/fixedpointiterator_tmpl.cc

@@ -10,6 +10,8 @@
 #include <dune/solvers/norms/energynorm.hh>
 #include <dune/tnnmg/problem-classes/convexproblem.hh>
 
+#include <dune/contact/common/deformedcontinuacomplex.hh>
+
 #include <dune/tectonic/globalfriction.hh>
 #include <dune/tectonic/myblockproblem.hh>
 
@@ -19,10 +21,8 @@
 #include "solverfactory.hh"
 
 using Function = Dune::VirtualFunction<double, double>;
-using Factory = SolverFactory<
-    MY_DIM,
-    MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
-    Grid>;
+using Factory = SolverFactory<DeformedGrid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
+
 using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
 using MyRateUpdater = RateUpdater<Vector, Matrix, Function, MY_DIM>;
 using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;

src/spatial-solving/solverfactory.cc

@@ -6,54 +6,80 @@
 #undef HAVE_IPOPT
 #endif
 
+#include <dune/common/bitsetvector.hh>
+
 #include <dune/fufem/assemblers/transferoperatorassembler.hh>
 #include <dune/solvers/solvers/solver.hh>
 
 #include "solverfactory.hh"
 
-template <size_t dim, class BlockProblem, class Grid>
-SolverFactory<dim, BlockProblem, Grid>::SolverFactory(
-    Dune::ParameterTree const &parset, Grid const &grid,
-    Dune::BitSetVector<dim> const &ignoreNodes)
-    : baseEnergyNorm(linearBaseSolverStep),
-      linearBaseSolver(&linearBaseSolverStep,
+template <class DeformedGrid, class Nonlinearity, class Matrix, class Vector>
+SolverFactory<DeformedGrid, Nonlinearity, Matrix, Vector>::SolverFactory(
+    const Dune::ParameterTree& parset, const Dune::Contact::NBodyAssembler<DeformedGrid, Vector>& nBodyAssembler,
+    const std::vector<Dune::BitSetVector<DeformedGrid::dimension>>& ignoreNodes)
+    : nBodyAssembler_(nBodyAssembler),
+      baseEnergyNorm_(baseSolverStep_),
+      baseSolver_(&baseSolverStep_,
                        parset.get<size_t>("linear.maxiumumIterations"),
-                       parset.get<double>("linear.tolerance"), &baseEnergyNorm,
+                       parset.get<double>("linear.tolerance"), &baseEnergyNorm_,
                        Solver::QUIET),
-      transferOperators(grid.maxLevel()),
-      multigridStep(
-          std::make_shared<Step>(linearIterationStep, nonlinearSmoother)) {
-  // linear iteration step
-  linearIterationStep.setMGType(parset.get<int>("linear.cycle"),
-                                parset.get<int>("linear.pre"),
-                                parset.get<int>("linear.post"));
-  linearIterationStep.basesolver_ = &linearBaseSolver;
-  linearIterationStep.setSmoother(&linearPresmoother, &linearPostsmoother);
-
-  // transfer operators
-  for (auto &&x : transferOperators)
-    x = new CompressedMultigridTransfer<Vector>;
-  TransferOperatorAssembler<Grid>(grid)
-      .assembleOperatorPointerHierarchy(transferOperators);
-  linearIterationStep.setTransferOperators(transferOperators);
+      transferOperators_(nBodyAssembler.getGrids().at(0)->maxLevel()) {
 
+  multigridStep_ = std::make_shared<Step>();
   // tnnmg iteration step
-  multigridStep->setSmoothingSteps(parset.get<int>("main.pre"),
-                                   parset.get<int>("main.multi"),
-                                   parset.get<int>("main.post"));
-  multigridStep->ignoreNodes_ = &ignoreNodes;
+  multigridStep_->setMGType(parset.get<int>("main.multi"), parset.get<int>("main.pre"), parset.get<int>("main.post"));
+  //multigridStep_->setIgnore(ignoreNodes);
+  multigridStep_->setBaseSolver(baseSolver_);
+  multigridStep_->setSmoother(&presmoother_, &postsmoother_);
+ // multigridStep_->setHasObstacle(nBodyAssembler_.totalHasObstacle_);
+ // multigridStep_->setObstacles(nBodyAssembler_.totalObstacles_);
+
+  // create the transfer operators
+  const int nCouplings = nBodyAssembler_.nCouplings();
+  const auto grids = nBodyAssembler_.getGrids();
+  for (size_t i=0; i<transferOperators_.size(); i++) {
+    std::vector<const Dune::BitSetVector<1>*> coarseHasObstacle(nCouplings);
+    std::vector<const Dune::BitSetVector<1>*> fineHasObstacle(nCouplings);
+
+    std::vector<const Matrix*> mortarTransfer(nCouplings);
+    std::vector<std::array<int,2> > gridIdx(nCouplings);
+
+    /*
+    for (int j=0; j<nCouplings; j++) {
+      coarseHasObstacle[j]  = nBodyAssembler_.contactCoupling_[j]->nonmortarBoundary_[i].getVertices();
+      fineHasObstacle[j]    = nBodyAssembler_.nonmortarBoundary_[j][i+1].getVertices();
+
+      mortarTransfer[j] = &nBodyAssembler_.contactCoupling_[j]->mortarLagrangeMatrix(i);
+      gridIdx[j]        = nBodyAssembler_.coupling_[j].gridIdx_;
+    }
+
+    transferOperators_[i] = new Dune::Contact::ContactMGTransfer<Vector>;
+    transferOperators_[i]->setup(grids, i, i+1,
+                                    nBodyAssembler_.localCoordSystems_[i],
+                                    nBodyAssembler_.localCoordSystems_[i+1],
+                                    coarseHasObstacle, fineHasObstacle,
+                                    mortarTransfer,
+                                    gridIdx); */
+  }
+/*
+  grids[0], *grids[1], i,
+                          contactAssembler.contactCoupling_[0]->mortarLagrangeMatrix(),
+                          contactAssembler.localCoordSystems_,
+                          *contactAssembler.contactCoupling_[0]->nonmortarBoundary().getVertices());
+*/
+  multigridStep_->setTransferOperators(transferOperators_);
 }
 
-template <size_t dim, class BlockProblem, class Grid>
-SolverFactory<dim, BlockProblem, Grid>::~SolverFactory() {
-  for (auto &&x : transferOperators)
+template <class DeformedGrid, class Nonlinearity, class Matrix, class Vector>
+SolverFactory<DeformedGrid, Nonlinearity, Matrix, Vector>::~SolverFactory() {
+  for (auto &&x : transferOperators_)
     delete x;
 }
 
-template <size_t dim, class BlockProblem, class Grid>
-auto SolverFactory<dim, BlockProblem, Grid>::getStep()
+template <class DeformedGrid, class Nonlinearity, class Matrix, class Vector>
+auto SolverFactory<DeformedGrid, Nonlinearity, Matrix, Vector>::getStep()
     -> std::shared_ptr<Step> {
-  return multigridStep;
+  return multigridStep_;
 }
 
 #include "solverfactory_tmpl.cc"

src/spatial-solving/solverfactory.hh

@@ -5,44 +5,61 @@
 #include <dune/common/parametertree.hh>
 
 #include <dune/solvers/iterationsteps/multigridstep.hh>
-#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
+#include <dune/solvers/iterationsteps/projectedblockgsstep.hh>
 #include <dune/solvers/norms/energynorm.hh>
 #include <dune/solvers/solvers/loopsolver.hh>
-#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
-#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
-#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
+//#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
+//#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
+//#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
 
-template <size_t dim, class BlockProblemTEMPLATE, class Grid>
+#include <dune/contact/assemblers/nbodyassembler.hh>
+
+//#include <dune/contact/estimators/hierarchiccontactestimator.hh>
+#include <dune/contact/solvers/nsnewtonmgstep.hh>
+#include <dune/contact/solvers/contacttransfer.hh>
+#include <dune/contact/solvers/contactobsrestrict.hh>
+#include <dune/contact/solvers/contacttransferoperatorassembler.hh>
+#include <dune/contact/solvers/nsnewtoncontacttransfer.hh>
+
+#define USE_OLD_TNNMG //needed for old bisection.hh in tnnmg
+
+template <class DeformedGridTEMPLATE, class NonlinearityTEMPLATE, class MatrixType, class VectorType>
 class SolverFactory {
+//protected:
+//  const size_t dim = DeformedGrid::dimension;
+
 public:
-  using BlockProblem = BlockProblemTEMPLATE;
-  using ConvexProblem = typename BlockProblem::ConvexProblemType;
-  using Vector = typename BlockProblem::VectorType;
-  using Matrix = typename BlockProblem::MatrixType;
+  using Vector = VectorType;
+  using Matrix = MatrixType;
 
-private:
-  using NonlinearSmoother = GenericNonlinearGS<BlockProblem>;
+  using DeformedGrid = DeformedGridTEMPLATE;
+  using Nonlinearity = NonlinearityTEMPLATE;
 
 public:
-  using Step =
-      TruncatedNonsmoothNewtonMultigrid<BlockProblem, NonlinearSmoother>;
+  using Step = Dune::Contact::NonSmoothNewtonMGStep<Matrix, Vector>;
 
-  SolverFactory(Dune::ParameterTree const &parset, Grid const &grid,
-                Dune::BitSetVector<dim> const &ignoreNodes);
+  SolverFactory(Dune::ParameterTree const &parset, const Dune::Contact::NBodyAssembler<DeformedGrid, Vector>& nBodyAssembler,
+                const std::vector<Dune::BitSetVector<DeformedGrid::dimension>>& ignoreNodes);
 
   ~SolverFactory();
 
   std::shared_ptr<Step> getStep();
 
+  const Dune::Contact::NBodyAssembler<DeformedGrid, Vector>& getNBodyAssembler() const {
+    return nBodyAssembler_;
+  }
+
 private:
-  TruncatedBlockGSStep<Matrix, Vector> linearBaseSolverStep;
-  EnergyNorm<Matrix, Vector> baseEnergyNorm;
-  LoopSolver<Vector> linearBaseSolver;
-  TruncatedBlockGSStep<Matrix, Vector> linearPresmoother;
-  TruncatedBlockGSStep<Matrix, Vector> linearPostsmoother;
-  MultigridStep<Matrix, Vector> linearIterationStep;
-  std::vector<CompressedMultigridTransfer<Vector> *> transferOperators;
-  NonlinearSmoother nonlinearSmoother;
-  std::shared_ptr<Step> multigridStep;
+  const Dune::Contact::NBodyAssembler<DeformedGrid, Vector>& nBodyAssembler_;
+
+  ProjectedBlockGSStep<Matrix, Vector> baseSolverStep_;
+  EnergyNorm<Matrix, Vector> baseEnergyNorm_;
+  LoopSolver<Vector> baseSolver_;
+
+  ProjectedBlockGSStep<Matrix, Vector> presmoother_;
+  ProjectedBlockGSStep<Matrix, Vector> postsmoother_;
+  std::shared_ptr<Step> multigridStep_;
+
+  std::vector<Dune::Contact::ContactMGTransfer<Vector>* > transferOperators_;
 };
 #endif

src/spatial-solving/solverfactory_tmpl.cc

@@ -12,11 +12,9 @@
 #include <dune/tectonic/globalfriction.hh>
 #include <dune/tectonic/myblockproblem.hh>
 
-template class SolverFactory<
-    MY_DIM,
-    MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
-    Grid>;
-template class SolverFactory<
+template class SolverFactory<DeformedGrid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
+
+/*template class SolverFactory<
     MY_DIM, BlockNonlinearTNNMGProblem<ConvexProblem<
                 ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
-    Grid>;
+    Grid>;*/

src/time-stepping/adaptivetimestepper.cc

@@ -20,10 +20,10 @@ void IterationRegister::reset() {
 template <class Factory, class Updaters, class ErrorNorm>
 AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::AdaptiveTimeStepper(
     Factory &factory, Dune::ParameterTree const &parset,
-    std::shared_ptr<Nonlinearity> globalFriction, Updaters &current,
+    std::vector<std::shared_ptr<Nonlinearity>>& globalFriction, Updaters &current,
     double relativeTime, double relativeTau,
-    std::function<void(double, Vector &)> externalForces,
-    ErrorNorm const &errorNorm,
+    std::vector<std::function<void(double, Vector &)>>& externalForces,
+    const std::vector<const ErrorNorm* >& errorNorms,
     std::function<bool(Updaters &, Updaters &)> mustRefine)
     : relativeTime_(relativeTime),
       relativeTau_(relativeTau),
@@ -35,7 +35,7 @@ AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::AdaptiveTimeStepper(
       R1_(),
       externalForces_(externalForces),
       mustRefine_(mustRefine),
-      errorNorm_(errorNorm) {}
+      errorNorms_(errorNorms) {}
 
 template <class Factory, class Updaters, class ErrorNorm>
 bool AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::reachedEnd() {
@@ -100,7 +100,7 @@ AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::step(
   UpdatersWithCount newUpdatersAndCount = {oldUpdaters.clone(), {}};
   newUpdatersAndCount.count =
       MyCoupledTimeStepper(finalTime_, factory_, parset_, globalFriction_,
-                           newUpdatersAndCount.updaters, errorNorm_,
+                           newUpdatersAndCount.updaters, errorNorms_,
                            externalForces_)
           .step(rTime, rTau);
   iterationRegister_.registerCount(newUpdatersAndCount.count);

src/time-stepping/adaptivetimestepper.hh

@@ -23,18 +23,18 @@ class AdaptiveTimeStepper {
   };
 
   using Vector = typename Factory::Vector;
-  using ConvexProblem = typename Factory::ConvexProblem;
-  using Nonlinearity = typename ConvexProblem::NonlinearityType;
+  //using ConvexProblem = typename Factory::ConvexProblem;
+  using Nonlinearity = typename Factory::Nonlinearity;
 
   using MyCoupledTimeStepper = CoupledTimeStepper<Factory, Updaters, ErrorNorm>;
 
 public:
   AdaptiveTimeStepper(Factory &factory, Dune::ParameterTree const &parset,
-                      std::shared_ptr<Nonlinearity> globalFriction,
+                      std::vector<std::shared_ptr<Nonlinearity>>& globalFriction,
                       Updaters &current, double relativeTime,
                       double relativeTau,
-                      std::function<void(double, Vector &)> externalForces,
-                      ErrorNorm const &errorNorm,
+                      std::vector<std::function<void(double, Vector &)>>& externalForces,
+                      const std::vector<const ErrorNorm* >& errorNorms,
                       std::function<bool(Updaters &, Updaters &)> mustRefine);
 
   bool reachedEnd();
@@ -50,12 +50,12 @@ class AdaptiveTimeStepper {
   double finalTime_;
   Factory &factory_;
   Dune::ParameterTree const &parset_;
-  std::shared_ptr<Nonlinearity> globalFriction_;
+  std::vector<std::shared_ptr<Nonlinearity>>& globalFriction_;
   Updaters &current_;
   UpdatersWithCount R1_;
-  std::function<void(double, Vector &)> externalForces_;
+  std::vector<std::function<void(double, Vector &)>>& externalForces_;
   std::function<bool(Updaters &, Updaters &)> mustRefine_;
-  ErrorNorm const &errorNorm_;
+  const std::vector<const ErrorNorm* >& errorNorms_;
 
   IterationRegister iterationRegister_;
 };

src/time-stepping/adaptivetimestepper_tmpl.cc

@@ -19,10 +19,7 @@
 #include "updaters.hh"
 
 using Function = Dune::VirtualFunction<double, double>;
-using Factory = SolverFactory<
-    MY_DIM,
-    MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
-    Grid>;
+using Factory = SolverFactory<DeformedGrid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
 using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
 using MyRateUpdater = RateUpdater<Vector, Matrix, Function, MY_DIM>;
 using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;

src/time-stepping/coupledtimestepper.cc

@@ -7,16 +7,16 @@
 template <class Factory, class Updaters, class ErrorNorm>
 CoupledTimeStepper<Factory, Updaters, ErrorNorm>::CoupledTimeStepper(
     double finalTime, Factory &factory, Dune::ParameterTree const &parset,
-    std::shared_ptr<Nonlinearity> globalFriction, Updaters updaters,
-    ErrorNorm const &errorNorm,
-    std::function<void(double, Vector &)> externalForces)
+    std::vector<std::shared_ptr<Nonlinearity>>& globalFriction, Updaters updaters,
+    const std::vector<const ErrorNorm* >& errorNorms,
+    std::vector<std::function<void(double, Vector &)>>& externalForces)
     : finalTime_(finalTime),
       factory_(factory),
       parset_(parset),
       globalFriction_(globalFriction),
       updaters_(updaters),
       externalForces_(externalForces),
-      errorNorm_(errorNorm) {}
+      errorNorms_(errorNorms) {}
 
 template <class Factory, class Updaters, class ErrorNorm>
 FixedPointIterationCounter
@@ -26,19 +26,21 @@ CoupledTimeStepper<Factory, Updaters, ErrorNorm>::step(double relativeTime,
   updaters_.rate_->nextTimeStep();
 
   auto const newRelativeTime = relativeTime + relativeTau;
-  Vector ell;
-  externalForces_(newRelativeTime, ell);
+  std::vector<Vector> ell(externalForces_.size());
+  for (size_t i=0; i<externalForces_.size(); i++) {
+    externalForces_[i](newRelativeTime, ell[i]);
+  }
 
-  Matrix velocityMatrix;
-  Vector velocityRHS;
-  Vector velocityIterate;
+  std::vector<Matrix> velocityMatrix;
+  std::vector<Vector> velocityRHS;
+  std::vector<Vector> velocityIterate;
 
   auto const tau = relativeTau * finalTime_;
-  updaters_.state_->setup(tau);
-  updaters_.rate_->setup(ell, tau, newRelativeTime, velocityRHS,
-                         velocityIterate, velocityMatrix);
+  updaters_.state_->setup(tau); 
+  updaters_.rate_->setup(ell, tau, newRelativeTime, velocityRHS, velocityIterate, velocityMatrix);
+
   FixedPointIterator<Factory, Updaters, ErrorNorm> fixedPointIterator(
-      factory_, parset_, globalFriction_, errorNorm_);
+      factory_, parset_, globalFriction_, errorNorms_);
   auto const iterations = fixedPointIterator.run(updaters_, velocityMatrix,
                                                  velocityRHS, velocityIterate);
   return iterations;

src/time-stepping/coupledtimestepper.hh

@@ -12,15 +12,14 @@ template <class Factory, class Updaters, class ErrorNorm>
 class CoupledTimeStepper {
   using Vector = typename Factory::Vector;
   using Matrix = typename Factory::Matrix;
-  using ConvexProblem = typename Factory::ConvexProblem;
-  using Nonlinearity = typename ConvexProblem::NonlinearityType;
+  using Nonlinearity = typename Factory::Nonlinearity;
 
 public:
   CoupledTimeStepper(double finalTime, Factory &factory,
                      Dune::ParameterTree const &parset,
-                     std::shared_ptr<Nonlinearity> globalFriction,
-                     Updaters updaters, ErrorNorm const &errorNorm,
-                     std::function<void(double, Vector &)> externalForces);
+                     std::vector<std::shared_ptr<Nonlinearity>>& globalFriction,
+                     Updaters updaters, const std::vector<const ErrorNorm* >& errorNorms,
+                     std::vector<std::function<void(double, Vector &)>>& externalForces);
 
   FixedPointIterationCounter step(double relativeTime, double relativeTau);
 
@@ -28,9 +27,9 @@ class CoupledTimeStepper {
   double finalTime_;
   Factory &factory_;
   Dune::ParameterTree const &parset_;
-  std::shared_ptr<Nonlinearity> globalFriction_;
+  std::vector<std::shared_ptr<Nonlinearity>>& globalFriction_;
   Updaters updaters_;
-  std::function<void(double, Vector &)> externalForces_;
-  ErrorNorm const &errorNorm_;
+  std::vector<std::function<void(double, Vector &)>>& externalForces_;
+  const std::vector<const ErrorNorm* >& errorNorms_;
 };
 #endif

src/time-stepping/coupledtimestepper_tmpl.cc

@@ -19,10 +19,8 @@
 #include "updaters.hh"
 
 using Function = Dune::VirtualFunction<double, double>;
-using Factory = SolverFactory<
-    MY_DIM,
-    MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
-    Grid>;
+using Factory = SolverFactory<DeformedGrid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
+
 using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
 using MyRateUpdater = RateUpdater<Vector, Matrix, Function, MY_DIM>;
 using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;

src/time-stepping/rate.cc

@@ -9,20 +9,22 @@
 template <class Vector, class Matrix, class Function, int dimension>
 std::shared_ptr<RateUpdater<Vector, Matrix, Function, dimension>>
 initRateUpdater(Config::scheme scheme,
-                Function const &velocityDirichletFunction,
-                Dune::BitSetVector<dimension> const &velocityDirichletNodes,
-                Matrices<Matrix> const &matrices, Vector const &u_initial,
-                Vector const &v_initial, Vector const &a_initial) {
+                const std::vector<const Function*>& velocityDirichletFunctions,
+                const std::vector<Dune::BitSetVector<dimension>>& velocityDirichletNodes,
+                const Matrices<Matrix>& matrices,
+                const std::vector<Vector>& u_initial,
+                const std::vector<Vector>& v_initial,
+                const std::vector<Vector>& a_initial) {
   switch (scheme) {
     case Config::Newmark:
       return std::make_shared<Newmark<Vector, Matrix, Function, dimension>>(
           matrices, u_initial, v_initial, a_initial, velocityDirichletNodes,
-          velocityDirichletFunction);
+          velocityDirichletFunctions);
     case Config::BackwardEuler:
       return std::make_shared<
           BackwardEuler<Vector, Matrix, Function, dimension>>(
           matrices, u_initial, v_initial, a_initial, velocityDirichletNodes,
-          velocityDirichletFunction);
+          velocityDirichletFunctions);
     default:
       assert(false);
   }

src/time-stepping/rate.hh

@@ -9,8 +9,10 @@
 template <class Vector, class Matrix, class Function, int dimension>
 std::shared_ptr<RateUpdater<Vector, Matrix, Function, dimension>>
 initRateUpdater(Config::scheme scheme,
-                Function const &velocityDirichletFunction,
-                Dune::BitSetVector<dimension> const &velocityDirichletNodes,
-                Matrices<Matrix> const &matrices, Vector const &u_initial,
-                Vector const &v_initial, Vector const &a_initial);
+                const std::vector<const Function* >& velocityDirichletFunctions,
+                const std::vector<Dune::BitSetVector<dimension>>& velocityDirichletNodes,
+                const Matrices<Matrix>& matrices,
+                const std::vector<Vector>& u_initial,
+                const std::vector<Vector>& v_initial,
+                const std::vector<Vector>& a_initial);
 #endif

src/time-stepping/rate/backward_euler.cc

-#include <dune/solvers/common/arithmetic.hh>
+#include <dune/matrix-vector/axpy.hh>
+#include <dune/istl/matrixindexset.hh>
 
 #include "backward_euler.hh"
 
 template <class Vector, class Matrix, class Function, size_t dim>
 BackwardEuler<Vector, Matrix, Function, dim>::BackwardEuler(
-    Matrices<Matrix> const &_matrices, Vector const &_u_initial,
-    Vector const &_v_initial, Vector const &_a_initial,
-    Dune::BitSetVector<dim> const &_dirichletNodes,
-    Function const &_dirichletFunction)
+    const Matrices<Matrix>& _matrices, const std::vector<Vector>& _u_initial,
+    const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
+    const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
+    const std::vector<const Function*>& _dirichletFunctions)
     : RateUpdater<Vector, Matrix, Function, dim>(
           _matrices, _u_initial, _v_initial, _a_initial, _dirichletNodes,
-          _dirichletFunction) {}
+          _dirichletFunctions) {}
 
 template <class Vector, class Matrix, class Function, size_t dim>
-void BackwardEuler<Vector, Matrix, Function, dim>::setup(
-    Vector const &ell, double _tau, double relativeTime, Vector &rhs,
-    Vector &iterate, Matrix &AM) {
-  this->dirichletFunction.evaluate(relativeTime, this->dirichletValue);
-  this->tau = _tau;
-
-  /* We start out with the formulation
-
-       M a + C v + A u = ell
-
-     Backward Euler means
-
-       a1 = 1.0/tau ( v1 - v0 )
-       u1 = tau v1 + u0
-
-     in summary, we get at time t=1
-
-       M [1.0/tau ( v1 - v0 )] + C v1
-       + A [tau v1 + u0] = ell
-
-     or
-
-       1.0/tau M v1 + C v1 + tau A v1
-       = [1.0/tau M + C + tau A] v1
-       = ell + 1.0/tau M v0 - A u0
-  */
-
-  // set up LHS (for fixed tau, we'd only really have to do this once)
-  {
-    Dune::MatrixIndexSet indices(this->matrices.elasticity.N(),
-                                 this->matrices.elasticity.M());
-    indices.import(this->matrices.elasticity);
-    indices.import(this->matrices.mass);
-    indices.import(this->matrices.damping);
-    indices.exportIdx(AM);
+void BackwardEuler<Vector, Matrix, Function, dim>::setup(const std::vector<Vector>& ell,
+                                                         double _tau,
+                                                         double relativeTime,
+                                                         std::vector<Vector>& rhs, std::vector<Vector>& iterate,
+                                                         std::vector<Matrix>& AM) {
+  for (size_t i=0; i<this->u.size(); i++) {
+    this->dirichletFunctions[i]->evaluate(relativeTime, this->dirichletValue);
+    this->tau = _tau;
+
+    /* We start out with the formulation
+
+           M a + C v + A u = ell
+
+         Backward Euler means
+
+           a1 = 1.0/tau ( v1 - v0 )
+           u1 = tau v1 + u0
+
+         in summary, we get at time t=1
+
+           M [1.0/tau ( v1 - v0 )] + C v1
+           + A [tau v1 + u0] = ell
+
+         or
+
+           1.0/tau M v1 + C v1 + tau A v1
+           = [1.0/tau M + C + tau A] v1
+           = ell + 1.0/tau M v0 - A u0
+    */
+
+    // set up LHS (for fixed tau, we'd only really have to do this once)
+    Matrix& LHS = AM[i];
+    {
+      Dune::MatrixIndexSet indices(this->matrices.elasticity[i]->N(),
+                                   this->matrices.elasticity[i]->M());
+      indices.import(*this->matrices.elasticity[i]);
+      indices.import(*this->matrices.mass[i]);
+      indices.import(*this->matrices.damping[i]);
+      indices.exportIdx(LHS);
+    }
+    LHS = 0.0;
+    Dune::MatrixVector::addProduct(LHS, 1.0 / this->tau, *this->matrices.mass[i]);
+    Dune::MatrixVector::addProduct(LHS, 1.0, *this->matrices.damping[i]);
+    Dune::MatrixVector::addProduct(LHS, this->tau, *this->matrices.elasticity[i]);
+
+    // set up RHS
+    {
+      Vector& rhss = rhs[i];
+      rhss = ell[i];
+      Dune::MatrixVector::addProduct(rhss, 1.0 / this->tau, *this->matrices.mass[i],
+                           this->v_o[i]);
+      Dune::MatrixVector::subtractProduct(rhss, *this->matrices.elasticity[i], this->u_o[i]);
+    }
+
+    iterate = this->v_o;
+
+    const Dune::BitSetVector<dim>& dirichletNodess = this->dirichletNodes[i];
+    for (size_t k = 0; k < dirichletNodess.size(); ++k)
+      for (size_t j = 0; j < dim; ++j)
+        if (dirichletNodess[k][j])
+          iterate[k][j] = (j == 0) ? this->dirichletValue : 0;
   }
-  AM = 0.0;
-  Arithmetic::addProduct(AM, 1.0 / this->tau, this->matrices.mass);
-  Arithmetic::addProduct(AM, 1.0, this->matrices.damping);
-  Arithmetic::addProduct(AM, this->tau, this->matrices.elasticity);
-
-  // set up RHS
-  {
-    rhs = ell;
-    Arithmetic::addProduct(rhs, 1.0 / this->tau, this->matrices.mass,
-                           this->v_o);
-    Arithmetic::subtractProduct(rhs, this->matrices.elasticity, this->u_o);
-  }
-
-  iterate = this->v_o;
-
-  for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
-    for (size_t j = 0; j < dim; ++j)
-      if (this->dirichletNodes[i][j])
-        iterate[i][j] = (j == 0) ? this->dirichletValue : 0;
 }
 
 template <class Vector, class Matrix, class Function, size_t dim>
 void BackwardEuler<Vector, Matrix, Function, dim>::postProcess(
-    Vector const &iterate) {
+    const std::vector<Vector>& iterate) {
   this->postProcessCalled = true;
 
   this->v = iterate;
 
   this->u = this->u_o;
-  Arithmetic::addProduct(this->u, this->tau, this->v);
 
-  this->a = this->v;
-  this->a -= this->v_o;
-  this->a /= this->tau;
+  for (size_t i=0; i<this->u.size(); i++) {
+    Dune::MatrixVector::addProduct(this->u[i], this->tau, this->v[i]);
+
+    Vector& ai = this->a[i];
+    ai = this->v[i];
+    ai -= this->v_o[i];
+    ai /= this->tau;
+  }
 }
 
 template <class Vector, class Matrix, class Function, size_t dim>

src/time-stepping/rate/backward_euler.hh

@@ -4,14 +4,15 @@
 template <class Vector, class Matrix, class Function, size_t dim>
 class BackwardEuler : public RateUpdater<Vector, Matrix, Function, dim> {
 public:
-  BackwardEuler(Matrices<Matrix> const &_matrices, Vector const &_u_initial,
-                Vector const &_v_initial, Vector const &_a_initial,
-                Dune::BitSetVector<dim> const &_dirichletNodes,
-                Function const &_dirichletFunction);
+  BackwardEuler(const Matrices<Matrix> &_matrices, const std::vector<Vector> &_u_initial,
+                const std::vector<Vector> &_v_initial, const std::vector<Vector> &_a_initial,
+                const std::vector<Dune::BitSetVector<dim> > &_dirichletNodes,
+                const std::vector<const Function*> &_dirichletFunctions);
 
-  void setup(Vector const &, double, double, Vector &, Vector &,
-             Matrix &) override;
-  void postProcess(Vector const &) override;
+  void setup(const std::vector<Vector>&, double, double, std::vector<Vector>&, std::vector<Vector>&,
+             std::vector<Matrix>&) override;
+
+  void postProcess(const std::vector<Vector>&) override;
 
   std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>> clone()
       const override;

src/time-stepping/rate/newmark.cc

-#include <dune/solvers/common/arithmetic.hh>
+#include <dune/matrix-vector/axpy.hh>
+#include <dune/istl/matrixindexset.hh>
 
 #include "newmark.hh"
 
 template <class Vector, class Matrix, class Function, size_t dim>
 Newmark<Vector, Matrix, Function, dim>::Newmark(
-    Matrices<Matrix> const &_matrices, Vector const &_u_initial,
-    Vector const &_v_initial, Vector const &_a_initial,
-    Dune::BitSetVector<dim> const &_dirichletNodes,
-    Function const &_dirichletFunction)
+    const Matrices<Matrix>& _matrices, const std::vector<Vector>& _u_initial,
+    const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
+    const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
+    const std::vector<const Function*>& _dirichletFunctions)
     : RateUpdater<Vector, Matrix, Function, dim>(
           _matrices, _u_initial, _v_initial, _a_initial, _dirichletNodes,
-          _dirichletFunction) {}
+          _dirichletFunctions) {}
 
 template <class Vector, class Matrix, class Function, size_t dim>
-void Newmark<Vector, Matrix, Function, dim>::setup(Vector const &ell,
+void Newmark<Vector, Matrix, Function, dim>::setup(const std::vector<Vector>& ell,
                                                    double _tau,
                                                    double relativeTime,
-                                                   Vector &rhs, Vector &iterate,
-                                                   Matrix &AM) {
-  this->dirichletFunction.evaluate(relativeTime, this->dirichletValue);
-  this->tau = _tau;
+                                                   std::vector<Vector>& rhs, std::vector<Vector>& iterate,
+                                                   std::vector<Matrix>& AM) {
+  for (size_t i=0; i<this->u.size(); i++) {
+    this->dirichletFunctions[i]->evaluate(relativeTime, this->dirichletValue);
+    this->tau = _tau;
 
-  /* We start out with the formulation
+    /* We start out with the formulation
 
        M a + C v + A u = ell
 
@@ -41,58 +43,64 @@ void Newmark<Vector, Matrix, Function, dim>::setup(Vector const &ell,
        = [2/tau M + C + tau/2 A] v1
        = ell + 2/tau M v0 + M a0
        - tau/2 A v0 - A u0
-  */
-
-  // set up LHS (for fixed tau, we'd only really have to do this once)
-  {
-    Dune::MatrixIndexSet indices(this->matrices.elasticity.N(),
-                                 this->matrices.elasticity.M());
-    indices.import(this->matrices.elasticity);
-    indices.import(this->matrices.mass);
-    indices.import(this->matrices.damping);
-    indices.exportIdx(AM);
+    */
+
+    // set up LHS (for fixed tau, we'd only really have to do this once)
+    Matrix& LHS = AM[i];
+    {
+      Dune::MatrixIndexSet indices(this->matrices.elasticity[i]->N(),
+                                 this->matrices.elasticity[i]->M());
+      indices.import(*this->matrices.elasticity[i]);
+      indices.import(*this->matrices.mass[i]);
+      indices.import(*this->matrices.damping[i]);
+      indices.exportIdx(LHS);
+    }
+    LHS = 0.0;
+    Dune::MatrixVector::addProduct(LHS, 2.0 / this->tau, *this->matrices.mass[i]);
+    Dune::MatrixVector::addProduct(LHS, 1.0, *this->matrices.damping[i]);
+    Dune::MatrixVector::addProduct(LHS, this->tau / 2.0, *this->matrices.elasticity[i]);
+
+    // set up RHS
+    {
+      Vector& rhss = rhs[i];
+      rhss = ell[i];
+      Dune::MatrixVector::addProduct(rhss, 2.0 / this->tau, *this->matrices.mass[i],
+                           this->v_o[i]);
+      Dune::MatrixVector::addProduct(rhss, *this->matrices.mass[i], this->a_o[i]);
+      Dune::MatrixVector::subtractProduct(rhss, this->tau / 2.0, *this->matrices.elasticity[i],
+                                this->v_o[i]);
+      Dune::MatrixVector::subtractProduct(rhss, *this->matrices.elasticity[i], this->u_o[i]);
+    }
+
+    iterate = this->v_o;
+
+    const Dune::BitSetVector<dim>& dirichletNodess = this->dirichletNodes[i];
+    for (size_t k = 0; k < dirichletNodess.size(); ++k)
+      for (size_t j = 0; j < dim; ++j)
+        if (dirichletNodess[k][j])
+          iterate[k][j] = (j == 0) ? this->dirichletValue : 0;
   }
-  AM = 0.0;
-  Arithmetic::addProduct(AM, 2.0 / this->tau, this->matrices.mass);
-  Arithmetic::addProduct(AM, 1.0, this->matrices.damping);
-  Arithmetic::addProduct(AM, this->tau / 2.0, this->matrices.elasticity);
-
-  // set up RHS
-  {
-    rhs = ell;
-    Arithmetic::addProduct(rhs, 2.0 / this->tau, this->matrices.mass,
-                           this->v_o);
-    Arithmetic::addProduct(rhs, this->matrices.mass, this->a_o);
-    Arithmetic::subtractProduct(rhs, this->tau / 2.0, this->matrices.elasticity,
-                                this->v_o);
-    Arithmetic::subtractProduct(rhs, this->matrices.elasticity, this->u_o);
-  }
-
-  iterate = this->v_o;
-
-  for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
-    for (size_t j = 0; j < dim; ++j)
-      if (this->dirichletNodes[i][j])
-        iterate[i][j] = (j == 0) ? this->dirichletValue : 0;
 }
 
 template <class Vector, class Matrix, class Function, size_t dim>
-void Newmark<Vector, Matrix, Function, dim>::postProcess(
-    Vector const &iterate) {
+void Newmark<Vector, Matrix, Function, dim>::postProcess(const std::vector<Vector>& iterate) {
   this->postProcessCalled = true;
 
   this->v = iterate;
 
   // u1 = tau/2 ( v1 + v0 ) + u0
   this->u = this->u_o;
-  Arithmetic::addProduct(this->u, this->tau / 2.0, this->v);
-  Arithmetic::addProduct(this->u, this->tau / 2.0, this->v_o);
-
-  // a1 = 2/tau ( v1 - v0 ) - a0
-  this->a = 0.0;
-  Arithmetic::addProduct(this->a, 2.0 / this->tau, this->v);
-  Arithmetic::subtractProduct(this->a, 2.0 / this->tau, this->v_o);
-  Arithmetic::subtractProduct(this->a, 1.0, this->a_o);
+
+  for (size_t i=0; i<this->u.size(); i++) {
+    Dune::MatrixVector::addProduct(this->u[i], this->tau / 2.0, this->v[i]);
+    Dune::MatrixVector::addProduct(this->u[i], this->tau / 2.0, this->v_o[i]);
+
+    // a1 = 2/tau ( v1 - v0 ) - a0
+    this->a[i] = 0.0;
+    Dune::MatrixVector::addProduct(this->a[i], 2.0 / this->tau, this->v[i]);
+    Dune::MatrixVector::subtractProduct(this->a[i], 2.0 / this->tau, this->v_o[i]);
+    Dune::MatrixVector::subtractProduct(this->a[i], 1.0, this->a_o[i]);
+  }
 }
 
 template <class Vector, class Matrix, class Function, size_t dim>

src/time-stepping/rate/newmark.hh

@@ -4,14 +4,15 @@
 template <class Vector, class Matrix, class Function, size_t dim>
 class Newmark : public RateUpdater<Vector, Matrix, Function, dim> {
 public:
-  Newmark(Matrices<Matrix> const &_matrices, Vector const &_u_initial,
-          Vector const &_v_initial, Vector const &_a_initial,
-          Dune::BitSetVector<dim> const &_dirichletNodes,
-          Function const &_dirichletFunction);
+  Newmark(const Matrices<Matrix>& _matrices, const std::vector<Vector>& _u_initial,
+          const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
+          const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
+          const std::vector<const Function*>& _dirichletFunctions);
 
-  void setup(Vector const &, double, double, Vector &, Vector &,
-             Matrix &) override;
-  void postProcess(Vector const &) override;
+  void setup(const std::vector<Vector>&, double, double, std::vector<Vector>&, std::vector<Vector>&,
+             std::vector<Matrix>&) override;
+
+  void postProcess(const std::vector<Vector>&) override;
 
   std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>> clone()
       const override;

src/time-stepping/rate/rateupdater.cc

@@ -6,16 +6,16 @@
 
 template <class Vector, class Matrix, class Function, size_t dim>
 RateUpdater<Vector, Matrix, Function, dim>::RateUpdater(
-    Matrices<Matrix> const &_matrices, Vector const &_u_initial,
-    Vector const &_v_initial, Vector const &_a_initial,
-    Dune::BitSetVector<dim> const &_dirichletNodes,
-    Function const &_dirichletFunction)
+    const Matrices<Matrix>& _matrices, const std::vector<Vector>& _u_initial,
+    const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
+    const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
+    const std::vector<const Function*>& _dirichletFunctions)
     : matrices(_matrices),
       u(_u_initial),
       v(_v_initial),
       a(_a_initial),
       dirichletNodes(_dirichletNodes),
-      dirichletFunction(_dirichletFunction) {}
+      dirichletFunctions(_dirichletFunctions) {}
 
 template <class Vector, class Matrix, class Function, size_t dim>
 void RateUpdater<Vector, Matrix, Function, dim>::nextTimeStep() {
@@ -26,17 +26,15 @@ void RateUpdater<Vector, Matrix, Function, dim>::nextTimeStep() {
 }
 
 template <class Vector, class Matrix, class Function, size_t dim>
-void RateUpdater<Vector, Matrix, Function, dim>::extractDisplacement(
-    Vector &displacement) const {
+void RateUpdater<Vector, Matrix, Function, dim>::extractDisplacement(std::vector<Vector>& displacements) const {
   if (!postProcessCalled)
     DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
 
-  displacement = u;
+  displacements = u;
 }
 
 template <class Vector, class Matrix, class Function, size_t dim>
-void RateUpdater<Vector, Matrix, Function, dim>::extractVelocity(
-    Vector &velocity) const {
+void RateUpdater<Vector, Matrix, Function, dim>::extractVelocity(std::vector<Vector>& velocity) const {
   if (!postProcessCalled)
     DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
 
@@ -44,14 +42,12 @@ void RateUpdater<Vector, Matrix, Function, dim>::extractVelocity(
 }
 
 template <class Vector, class Matrix, class Function, size_t dim>
-void RateUpdater<Vector, Matrix, Function, dim>::extractOldVelocity(
-    Vector &oldVelocity) const {
+void RateUpdater<Vector, Matrix, Function, dim>::extractOldVelocity(std::vector<Vector>& oldVelocity) const {
   oldVelocity = v_o;
 }
 
 template <class Vector, class Matrix, class Function, size_t dim>
-void RateUpdater<Vector, Matrix, Function, dim>::extractAcceleration(
-    Vector &acceleration) const {
+void RateUpdater<Vector, Matrix, Function, dim>::extractAcceleration(std::vector<Vector>& acceleration) const {
   if (!postProcessCalled)
     DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
 

src/time-stepping/rate/rateupdater.hh

@@ -10,33 +10,33 @@
 template <class Vector, class Matrix, class Function, size_t dim>
 class RateUpdater {
 protected:
-  RateUpdater(Matrices<Matrix> const &_matrices, Vector const &_u_initial,
-              Vector const &_v_initial, Vector const &_a_initial,
-              Dune::BitSetVector<dim> const &_dirichletNodes,
-              Function const &_dirichletFunction);
+  RateUpdater(const Matrices<Matrix>& _matrices, const std::vector<Vector>& _u_initial,
+              const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
+              const std::vector<Dune::BitSetVector<dim>>& _dirichletNodes,
+              const std::vector<const Function*>& _dirichletFunctions);
 
 public:
   void nextTimeStep();
-  void virtual setup(Vector const &ell, double _tau, double relativeTime,
-                     Vector &rhs, Vector &iterate, Matrix &AB) = 0;
+  void virtual setup(const std::vector<Vector>& ell, double _tau, double relativeTime,
+                     std::vector<Vector>& rhs, std::vector<Vector>& iterate, std::vector<Matrix>& AB) = 0;
 
-  void virtual postProcess(Vector const &iterate) = 0;
-  void extractDisplacement(Vector &displacement) const;
-  void extractVelocity(Vector &velocity) const;
-  void extractOldVelocity(Vector &velocity) const;
-  void extractAcceleration(Vector &acceleration) const;
+  void virtual postProcess(const std::vector<Vector>& iterate) = 0;
+  void extractDisplacement(std::vector<Vector>& displacements) const;
+  void extractVelocity(std::vector<Vector>& velocity) const;
+  void extractOldVelocity(std::vector<Vector>& velocity) const;
+  void extractAcceleration(std::vector<Vector>& acceleration) const;
 
   std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>> virtual clone()
       const = 0;
 
 protected:
-  Matrices<Matrix> const &matrices;
-  Vector u, v, a;
-  Dune::BitSetVector<dim> const &dirichletNodes;
-  Function const &dirichletFunction;
+  const Matrices<Matrix>& matrices;
+  std::vector<Vector> u, v, a;
+  const std::vector<Dune::BitSetVector<dim>>& dirichletNodes;
+  const std::vector<const Function*>& dirichletFunctions;
   double dirichletValue;
 
-  Vector u_o, v_o, a_o;
+  std::vector<Vector> u_o, v_o, a_o;
   double tau;
 
   bool postProcessCalled = true;