From 7fcefaf505c22c5394ca739ec93236972a44b214 Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@mi.fu-berlin.de>
Date: Thu, 9 May 2019 15:48:45 +0200
Subject: [PATCH] .

---
 dune/tectonic/globalratestatefriction.hh      |  2 +-
 src/data-structures/levelcontactnetwork.cc    | 35 ++++++++++++-------
 src/data-structures/levelcontactnetwork.hh    |  1 +
 src/data-structures/program_state.hh          |  9 +++--
 src/multi-body-problem.cc                     |  8 ++---
 src/spatial-solving/fixedpointiterator.cc     | 24 +++++++++++--
 src/spatial-solving/fixedpointiterator.hh     |  4 +++
 src/spatial-solving/solverfactory.cc          |  4 ++-
 src/spatial-solving/solverfactory.hh          |  3 +-
 src/time-stepping/adaptivetimestepper.cc      |  4 ++-
 src/time-stepping/adaptivetimestepper.hh      |  4 +++
 src/time-stepping/coupledtimestepper.cc       |  4 ++-
 src/time-stepping/coupledtimestepper.hh       |  4 +++
 src/time-stepping/state.cc                    |  7 +++-
 .../state/ageinglawstateupdater.cc            |  4 ++-
 src/time-stepping/state/stateupdater.hh       |  2 +-
 16 files changed, 87 insertions(+), 32 deletions(-)

diff --git a/dune/tectonic/globalratestatefriction.hh b/dune/tectonic/globalratestatefriction.hh
index edfb05e0..078bd85d 100644
--- a/dune/tectonic/globalratestatefriction.hh
+++ b/dune/tectonic/globalratestatefriction.hh
@@ -88,7 +88,7 @@ class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
 
                 localToGlobal_.emplace_back(i);
                 restrictions_.emplace_back(weights[bodyIdx][i], weightedNormalStress[bodyIdx][i],
-                                          couplings[i]->frictionData()(geoToPoint(it->geometry())));
+                                          couplings[j]->frictionData()(geoToPoint(it->geometry())));
                 break;
             }
 
diff --git a/src/data-structures/levelcontactnetwork.cc b/src/data-structures/levelcontactnetwork.cc
index 9709a867..aa336513 100644
--- a/src/data-structures/levelcontactnetwork.cc
+++ b/src/data-structures/levelcontactnetwork.cc
@@ -27,6 +27,8 @@ LevelContactNetwork<GridType, dimension>::LevelContactNetwork(
     leafViews_(nBodies),
     levelViews_(nBodies),
     leafVertexCounts_(nBodies),
+    frictionBoundaries_(nCouplings),
+    stateEnergyNorms_(nBodies),
     nBodyAssembler_(nBodies, nCouplings)
 {}
 
@@ -34,6 +36,7 @@ template <class GridType, int dimension>
 LevelContactNetwork<GridType, dimension>::~LevelContactNetwork() {
     for (size_t i=0; i<nBodies(); i++) {
         delete frictionBoundaries_[i];
+        delete frictionBoundaryMass_[i];
         delete stateEnergyNorms_[i];
     }
 }
@@ -46,33 +49,37 @@ void LevelContactNetwork<GridType, dimension>::assemble() {
         bodies_[i]->assemble();
         //printBasisDofLocation(bodies_[i]->assembler()->vertexBasis); //TODO remove after debugging
 
-        frictionBoundaries_[i] = new BoundaryPatch(bodies_[i]->leafView(),false);
+        frictionBoundaries_[i] = new BoundaryPatch(bodies_[i]->leafView(), false);
     }
 
     // set up dune-contact nBodyAssembler
     nBodyAssembler_.setGrids(deformedGrids_);
 
     for (size_t i=0; i<nCouplings(); i++) {
-      auto& coupling = couplings_[i];
-      nBodyAssembler_.setCoupling(*coupling, i);
+      nBodyAssembler_.setCoupling(*couplings_[i], i);
+    }
+
+    nBodyAssembler_.assembleTransferOperator();
+    nBodyAssembler_.assembleObstacle();
 
+    for (size_t i=0; i<nCouplings(); i++) {
+      auto& coupling = couplings_[i];
       const auto& contactCoupling = nBodyAssembler_.getContactCouplings()[i]; // dual mortar object holding boundary patches
       const auto nonmortarIdx = coupling->gridIdx_[0];
-      const auto mortarIdx = coupling->gridIdx_[1];
+      //const auto mortarIdx = coupling->gridIdx_[1];
 
       frictionBoundaries_[nonmortarIdx]->addPatch(contactCoupling->nonmortarBoundary());
-      frictionBoundaries_[mortarIdx]->addPatch(contactCoupling->mortarBoundary());
+      //frictionBoundaries_[mortarIdx]->addPatch(contactCoupling->mortarBoundary());
     }
 
-    nBodyAssembler_.assembleTransferOperator();
-    nBodyAssembler_.assembleObstacle();
-
     // assemble state energy norm
+    frictionBoundaryMass_.resize(nBodies());
     for (size_t i=0; i<nBodies(); i++) {
-        ScalarMatrix relativeFrictionBoundaryMass;
-        bodies_[i]->assembler()->assembleFrictionalBoundaryMass(*frictionBoundaries_[i], relativeFrictionBoundaryMass);
-        relativeFrictionBoundaryMass /= frictionBoundaries_[i]->area(); // TODO: weight by individual friction patches?
-        stateEnergyNorms_[i] = new EnergyNorm<ScalarMatrix, ScalarVector>(relativeFrictionBoundaryMass);
+        frictionBoundaryMass_[i] = new ScalarMatrix();
+        bodies_[i]->assembler()->assembleFrictionalBoundaryMass(*frictionBoundaries_[i], *frictionBoundaryMass_[i]);
+        *frictionBoundaryMass_[i] /= frictionBoundaries_[i]->area(); // TODO: weight by individual friction patches?
+
+        stateEnergyNorms_[i] = new EnergyNorm<ScalarMatrix, ScalarVector>(*frictionBoundaryMass_[i]);
     }
 }
 
@@ -103,12 +110,15 @@ void LevelContactNetwork<GridType, dimension>::assembleFriction(
           globalFriction_ = std::make_shared<GlobalRateStateFriction<
               Matrix, Vector, TruncatedRateState, DeformedGridType>>(
               nBodyAssembler_.getContactCouplings(), couplings_, weights, weightedNormalStress);
+          break;
         case Config::Regularised:
           globalFriction_ = std::make_shared<GlobalRateStateFriction<
               Matrix, Vector, RegularisedRateState, DeformedGridType>>(
               nBodyAssembler_.getContactCouplings(), couplings_, weights, weightedNormalStress);
+          break;
         default:
           assert(false);
+          break;
       }
 }
 
@@ -266,6 +276,7 @@ void LevelContactNetwork<GridType, dimension>::frictionNodes(std::vector<const D
 template <class GridType, int dimension>
 auto LevelContactNetwork<GridType, dimension>::stateEnergyNorms() const
 -> const StateEnergyNorms& {
+
     return stateEnergyNorms_;
 }
 
diff --git a/src/data-structures/levelcontactnetwork.hh b/src/data-structures/levelcontactnetwork.hh
index 4721208b..97792f2e 100644
--- a/src/data-structures/levelcontactnetwork.hh
+++ b/src/data-structures/levelcontactnetwork.hh
@@ -143,6 +143,7 @@ template <class GridType, int dimension> class LevelContactNetwork {
         std::vector<size_t> leafVertexCounts_;
 
         std::vector<DeformedLeafBoundaryPatch*> frictionBoundaries_; // union of all boundary patches per body
+        std::vector<ScalarMatrix*> frictionBoundaryMass_;
         StateEnergyNorms stateEnergyNorms_;
 
         NBodyAssembler nBodyAssembler_;
diff --git a/src/data-structures/program_state.hh b/src/data-structures/program_state.hh
index ec0311fe..298f40f7 100644
--- a/src/data-structures/program_state.hh
+++ b/src/data-structures/program_state.hh
@@ -143,19 +143,19 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
       }
 
       // print problem
-      print(bilinearForm, "bilinearForm");
+     /* print(bilinearForm, "bilinearForm");
       print(totalRhs, "totalRhs");
       print(_dirichletNodes, "ignore");
       print(totalObstacles, "totalObstacles");
       print(lower, "lower");
-      print(upper, "upper");
+      print(upper, "upper");*/
 
       // set up functional and TNMMG solver
       using Functional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, typename Matrix::field_type>;
       Functional J(bilinearForm, totalRhs, ZeroNonlinearity(), lower, upper);
 
       using Factory = SolverFactory<Functional, Dune::BitSetVector<dims>>;
-      Factory factory(parset.sub("solver.tnnmg"), J);
+      Factory factory(parset.sub("solver.tnnmg"), J, _dirichletNodes);
 
    /*   std::vector<Dune::BitSetVector<dims>> bodyDirichletNodes;
       nBodyAssembler.postprocess(_dirichletNodes, bodyDirichletNodes);
@@ -168,7 +168,6 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
       print(totalRhs, "totalRhs: ");*/
 
       auto tnnmgStep = factory.step();
-      tnnmgStep->setIgnore(_dirichletNodes);
       factory.setProblem(totalX);
 
       const EnergyNorm<Matrix, Vector> norm(bilinearForm);
@@ -209,7 +208,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
       bilinearForm.mmv(oldStep->getSol(), oldRes);
       std::cout << "Old Res - energy norm: " << norm.operator ()(oldRes) << std::endl;*/
 
-      print(tnnmgStep->getSol(), "TNNMG Solution: ");
+   //   print(tnnmgStep->getSol(), "TNNMG Solution: ");
    /*   print(oldStep->getSol(), "Old Solution: ");
       auto diff = tnnmgStep->getSol();
       diff -= oldStep->getSol();
diff --git a/src/multi-body-problem.cc b/src/multi-body-problem.cc
index 709adff5..fda8fc8f 100644
--- a/src/multi-body-problem.cc
+++ b/src/multi-body-problem.cc
@@ -215,8 +215,8 @@ int main(int argc, char *argv[]) {
         vertexCoords[vertexMapper.index(v)] = geoToPoint(v.geometry());
     }
 
-    typename LevelContactNetwork::BoundaryPatches frictionBoundaries;
-    levelContactNetwork.boundaryPatches("friction", frictionBoundaries);
+    //typename LevelContactNetwork::BoundaryPatches frictionBoundaries;
+    //levelContactNetwork.boundaryPatches("friction", frictionBoundaries);
 
     /*
     auto dataWriter =
@@ -273,7 +273,7 @@ int main(int argc, char *argv[]) {
     // -------------------
     auto& nBodyAssembler = levelContactNetwork.nBodyAssembler();
 
-    Dune::BitSetVector<dims> totalDirichletNodes;
+    BitVector totalDirichletNodes;
     levelContactNetwork.totalNodes("dirichlet", totalDirichletNodes);
 
     using Functional = Functional<Matrix&, Vector&, GlobalFriction<Matrix, Vector>&, Vector&, Vector&, field_type>;
@@ -338,7 +338,7 @@ int main(int argc, char *argv[]) {
     levelContactNetwork.externalForces(externalForces);
 
     AdaptiveTimeStepper<NonlinearFactory, std::decay_t<decltype(nBodyAssembler)>, Updaters, EnergyNorm<ScalarMatrix, ScalarVector>>
-        adaptiveTimeStepper(parset, nBodyAssembler, globalFriction, frictionNodes, current,
+        adaptiveTimeStepper(parset, nBodyAssembler, totalDirichletNodes, globalFriction, frictionNodes, current,
                             programState.relativeTime, programState.relativeTau,
                             externalForces, stateEnergyNorms, mustRefine);
 
diff --git a/src/spatial-solving/fixedpointiterator.cc b/src/spatial-solving/fixedpointiterator.cc
index 5ca84926..ccf52350 100644
--- a/src/spatial-solving/fixedpointiterator.cc
+++ b/src/spatial-solving/fixedpointiterator.cc
@@ -27,8 +27,10 @@
 
 #include "fixedpointiterator.hh"
 
+#include "../utils/tobool.hh"
 #include "../utils/debugutils.hh"
 
+
 void FixedPointIterationCounter::operator+=(
     FixedPointIterationCounter const &other) {
   iterations += other.iterations;
@@ -39,11 +41,13 @@ template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
 FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::FixedPointIterator(
     Dune::ParameterTree const &parset,
     NBodyAssembler& nBodyAssembler,
+    const IgnoreVector& ignoreNodes,
     GlobalFriction& globalFriction,
     const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
     const std::vector<const ErrorNorm* >& errorNorms)
     : parset_(parset),
       nBodyAssembler_(nBodyAssembler),
+      ignoreNodes_(ignoreNodes),
       globalFriction_(globalFriction),
       bodywiseNonmortarBoundaries_(bodywiseNonmortarBoundaries),
       fixedPointMaxIterations_(parset.get<size_t>("v.fpi.maximumIterations")),
@@ -94,7 +98,7 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::run(
 
   // set up functional and TNMMG solver
   Functional J(bilinearForm, totalRhs, globalFriction_, lower, upper);
-  Factory solverFactory(parset_.sub("solver.tnnmg"), J);
+  Factory solverFactory(parset_.sub("solver.tnnmg"), J, ignoreNodes_);
   auto step = solverFactory.step();
 
   EnergyNorm<Matrix, Vector> energyNorm(bilinearForm);
@@ -183,10 +187,24 @@ void FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::relativeV
     const size_t nBodies = nBodyAssembler_.nGrids();
    // const auto& contactCouplings = nBodyAssembler_.getContactCouplings();
 
-
+    size_t globalIdx = 0;
     v_rel.resize(nBodies);
+    for (size_t bodyIdx=0; bodyIdx<nBodies; bodyIdx++) {
+        const auto& nonmortarBoundary = *bodywiseNonmortarBoundaries_[bodyIdx];
+        auto& v_rel_ = v_rel[bodyIdx];
+
+        v_rel_.resize(nonmortarBoundary.size());
+
+        for (size_t i=0; i<v_rel_.size(); i++) {
+            if (toBool(nonmortarBoundary[i])) {
+                v_rel_[i] = v[globalIdx];
+            }
+            globalIdx++;
+        }
+    }
+
 
-   /* for (size_t bodyIdx=0; bodyIdx<nBodies; bodyIdx++) {
+   /*
         boundaryNodes =
 
         const auto gridView = contactCouplings[couplingIndices[0]]->nonmortarBoundary().gridView();
diff --git a/src/spatial-solving/fixedpointiterator.hh b/src/spatial-solving/fixedpointiterator.hh
index 078b1d20..fa93c087 100644
--- a/src/spatial-solving/fixedpointiterator.hh
+++ b/src/spatial-solving/fixedpointiterator.hh
@@ -33,6 +33,7 @@ class FixedPointIterator {
 
   const static int dims = Vector::block_type::dimension;
 
+  using IgnoreVector = typename Factory::BitVector;
  // using Nonlinearity = typename ConvexProblem::NonlinearityType;
  //  using DeformedGrid = typename Factory::DeformedGrid;
 
@@ -46,6 +47,7 @@ class FixedPointIterator {
 public:
   FixedPointIterator(const Dune::ParameterTree& parset,
                      NBodyAssembler& nBodyAssembler,
+                     const IgnoreVector& ignoreNodes,
                      GlobalFriction& globalFriction,
                      const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
                      const std::vector<const ErrorNorm* >& errorNorms);
@@ -58,6 +60,8 @@ class FixedPointIterator {
 private:
   const Dune::ParameterTree& parset_;
   NBodyAssembler& nBodyAssembler_;
+  const IgnoreVector& ignoreNodes_;
+
   GlobalFriction& globalFriction_;
   const std::vector<const BitVector*>& bodywiseNonmortarBoundaries_;
 
diff --git a/src/spatial-solving/solverfactory.cc b/src/spatial-solving/solverfactory.cc
index 25bc7575..4804b788 100644
--- a/src/spatial-solving/solverfactory.cc
+++ b/src/spatial-solving/solverfactory.cc
@@ -10,7 +10,8 @@
 template <class Functional, class BitVector>
 SolverFactory<Functional, BitVector>::SolverFactory(
     const Dune::ParameterTree& parset,
-    Functional& J) :
+    Functional& J,
+    const BitVector& ignoreNodes) :
       J_(Dune::Solvers::wrap_own_share<const Functional>(std::forward<Functional>(J))),
       //linear solver
       //preconditioner_(),
@@ -24,6 +25,7 @@ SolverFactory<Functional, BitVector>::SolverFactory(
 
     tnnmgStep_ = std::make_shared<Step>(*J_, dummyIterate_, nonlinearSmoother_, linearSolver_, DefectProjection(), LineSearchSolver());
     tnnmgStep_->setPreSmoothingSteps(parset.get<int>("main.pre"));
+    tnnmgStep_->setIgnore(ignoreNodes);
 }
 
 template <class Functional, class BitVector>
diff --git a/src/spatial-solving/solverfactory.hh b/src/spatial-solving/solverfactory.hh
index c9a3d3b1..ebf30349 100644
--- a/src/spatial-solving/solverfactory.hh
+++ b/src/spatial-solving/solverfactory.hh
@@ -41,7 +41,8 @@ class SolverFactory {
 
 
   SolverFactory(Dune::ParameterTree const &,
-                Functional&);
+                Functional&,
+                const BitVector&);
 
   void setProblem(Vector& x);
 
diff --git a/src/time-stepping/adaptivetimestepper.cc b/src/time-stepping/adaptivetimestepper.cc
index bfde1aeb..c4a7371e 100644
--- a/src/time-stepping/adaptivetimestepper.cc
+++ b/src/time-stepping/adaptivetimestepper.cc
@@ -21,6 +21,7 @@ template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
 AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::AdaptiveTimeStepper(
         Dune::ParameterTree const &parset,
         NBodyAssembler& nBodyAssembler,
+        const IgnoreVector& ignoreNodes,
         GlobalFriction& globalFriction,
         const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
         Updaters &current,
@@ -34,6 +35,7 @@ AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::AdaptiveTimeS
       finalTime_(parset.get<double>("problem.finalTime")),
       parset_(parset),
       nBodyAssembler_(nBodyAssembler),
+      ignoreNodes_(ignoreNodes),
       globalFriction_(globalFriction),
       bodywiseNonmortarBoundaries_(bodywiseNonmortarBoundaries),
       current_(current),
@@ -104,7 +106,7 @@ AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::step(
     Updaters const &oldUpdaters, double rTime, double rTau) {
   UpdatersWithCount newUpdatersAndCount = {oldUpdaters.clone(), {}};
   newUpdatersAndCount.count =
-      MyCoupledTimeStepper(finalTime_, parset_, nBodyAssembler_, globalFriction_, bodywiseNonmortarBoundaries_,
+      MyCoupledTimeStepper(finalTime_, parset_, nBodyAssembler_, ignoreNodes_, globalFriction_, bodywiseNonmortarBoundaries_,
                            newUpdatersAndCount.updaters, errorNorms_,
                            externalForces_)
           .step(rTime, rTau);
diff --git a/src/time-stepping/adaptivetimestepper.hh b/src/time-stepping/adaptivetimestepper.hh
index d0a014d0..f1f9b07c 100644
--- a/src/time-stepping/adaptivetimestepper.hh
+++ b/src/time-stepping/adaptivetimestepper.hh
@@ -23,6 +23,7 @@ class AdaptiveTimeStepper {
   };
 
   using Vector = typename Factory::Vector;
+  using IgnoreVector = typename Factory::BitVector;
   //using ConvexProblem = typename Factory::ConvexProblem;
   //using Nonlinearity = typename Factory::Nonlinearity;
 
@@ -36,6 +37,7 @@ class AdaptiveTimeStepper {
   AdaptiveTimeStepper(
                       Dune::ParameterTree const &parset,
                       NBodyAssembler& nBodyAssembler,
+                      const IgnoreVector& ignoreNodes,
                       GlobalFriction& globalFriction,
                       const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
                       Updaters &current,
@@ -58,6 +60,8 @@ class AdaptiveTimeStepper {
   double finalTime_;
   Dune::ParameterTree const &parset_;
   NBodyAssembler& nBodyAssembler_;
+  const IgnoreVector& ignoreNodes_;
+
   GlobalFriction& globalFriction_;
   const std::vector<const BitVector*>& bodywiseNonmortarBoundaries_;
 
diff --git a/src/time-stepping/coupledtimestepper.cc b/src/time-stepping/coupledtimestepper.cc
index 0ee30c60..a102bb04 100644
--- a/src/time-stepping/coupledtimestepper.cc
+++ b/src/time-stepping/coupledtimestepper.cc
@@ -8,6 +8,7 @@ template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
 CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::CoupledTimeStepper(
     double finalTime, Dune::ParameterTree const &parset,
     NBodyAssembler& nBodyAssembler,
+    const IgnoreVector& ignoreNodes,
     GlobalFriction& globalFriction,
     const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
     Updaters updaters,
@@ -16,6 +17,7 @@ CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::CoupledTimeSte
     : finalTime_(finalTime),
       parset_(parset),
       nBodyAssembler_(nBodyAssembler),
+      ignoreNodes_(ignoreNodes),
       globalFriction_(globalFriction),
       bodywiseNonmortarBoundaries_(bodywiseNonmortarBoundaries),
       updaters_(updaters),
@@ -44,7 +46,7 @@ CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorm>::step(double re
   updaters_.rate_->setup(ell, tau, newRelativeTime, velocityRHS, velocityIterate, velocityMatrix);
 
   FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm> fixedPointIterator(
-      parset_, nBodyAssembler_, globalFriction_, bodywiseNonmortarBoundaries_, errorNorms_);
+      parset_, nBodyAssembler_, ignoreNodes_, globalFriction_, bodywiseNonmortarBoundaries_, errorNorms_);
   auto const iterations = fixedPointIterator.run(updaters_, velocityMatrix,
                                                  velocityRHS, velocityIterate);
   return iterations;
diff --git a/src/time-stepping/coupledtimestepper.hh b/src/time-stepping/coupledtimestepper.hh
index bc2b1f67..4610c7a7 100644
--- a/src/time-stepping/coupledtimestepper.hh
+++ b/src/time-stepping/coupledtimestepper.hh
@@ -12,6 +12,7 @@ template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
 class CoupledTimeStepper {
   using Vector = typename Factory::Vector;
   using Matrix = typename Factory::Matrix;
+  using IgnoreVector = typename Factory::BitVector;
   //using Nonlinearity = typename Factory::Nonlinearity;
 public:
   using GlobalFriction = typename FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::GlobalFriction;
@@ -22,6 +23,7 @@ class CoupledTimeStepper {
   CoupledTimeStepper(double finalTime,
                      Dune::ParameterTree const &parset,
                      NBodyAssembler& nBodyAssembler,
+                     const IgnoreVector& ignoreNodes,
                      GlobalFriction& globalFriction,
                      const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
                      Updaters updaters, const std::vector<const ErrorNorm* >& errorNorms,
@@ -33,6 +35,8 @@ class CoupledTimeStepper {
   double finalTime_;
   Dune::ParameterTree const &parset_;
   NBodyAssembler& nBodyAssembler_;
+  const IgnoreVector& ignoreNodes_;
+
   GlobalFriction& globalFriction_;
   const std::vector<const BitVector*>& bodywiseNonmortarBoundaries_;
 
diff --git a/src/time-stepping/state.cc b/src/time-stepping/state.cc
index 1e01e8a3..4184f3d7 100644
--- a/src/time-stepping/state.cc
+++ b/src/time-stepping/state.cc
@@ -26,9 +26,11 @@ auto initStateUpdater(
 
         auto localUpdater =  std::make_shared<AgeingLawStateUpdater<ScalarVector, Vector>>(
                     alpha_initial[nonmortarIdx], *contactCouplings[i]->nonmortarBoundary().getVertices(), coupling->frictionData().L(), coupling->frictionData().V0());
+        localUpdater->setBodyIndex(nonmortarIdx);
 
         stateUpdater->addLocalUpdater(localUpdater);
       }
+      break;
     case Config::SlipLaw:
       for (size_t i=0; i<couplings.size(); i++) {
         const auto& coupling = couplings[i];
@@ -36,14 +38,17 @@ auto initStateUpdater(
 
         auto localUpdater =  std::make_shared<SlipLawStateUpdater<ScalarVector, Vector>>(
                   alpha_initial[nonmortarIdx], *contactCouplings[i]->nonmortarBoundary().getVertices(), coupling->frictionData().L(), coupling->frictionData().V0());
+        localUpdater->setBodyIndex(nonmortarIdx);
 
         stateUpdater->addLocalUpdater(localUpdater);
       }
+      break;
     default:
       assert(false);
+      break;
+    }
 
     return stateUpdater;
-  }
 }
 
 #include "state_tmpl.cc"
diff --git a/src/time-stepping/state/ageinglawstateupdater.cc b/src/time-stepping/state/ageinglawstateupdater.cc
index 6862029e..a4c47377 100644
--- a/src/time-stepping/state/ageinglawstateupdater.cc
+++ b/src/time-stepping/state/ageinglawstateupdater.cc
@@ -66,7 +66,9 @@ template <class ScalarVector, class Vector>
 void AgeingLawStateUpdater<ScalarVector, Vector>::extractAlpha(
         ScalarVector& alpha) {
 
-    assert(alpha.size() == nodes_.size());
+    if (alpha.size() != nodes_.size()) {
+        alpha.resize(nodes_.size());
+    }
 
     for (size_t i=0; i<localToGlobal_.size(); i++) {
         alpha[localToGlobal_[i]] = alpha_[i];
diff --git a/src/time-stepping/state/stateupdater.hh b/src/time-stepping/state/stateupdater.hh
index ed7a5241..458fc81a 100644
--- a/src/time-stepping/state/stateupdater.hh
+++ b/src/time-stepping/state/stateupdater.hh
@@ -57,7 +57,7 @@ template <class ScalarVectorTEMPLATE, class Vector> class StateUpdater {
     }
   }
 
-  void extractAlpha(std::vector<ScalarVector>& alpha) {
+  void extractAlpha(std::vector<ScalarVector>& alpha) {      
     for (size_t i=0; i<localStateUpdaters_.size(); i++) {
         auto& localStateUpdater = localStateUpdaters_[i];
         localStateUpdater->extractAlpha(alpha[localStateUpdater->bodyIndex()]);
-- 
GitLab