save weights in programstate member

dune/tectonic/assemblers.cc

@@ -144,6 +144,7 @@ template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleWeightedNormalStress(
         BoundaryPatch<GridView> const &frictionalBoundary,
         ScalarVector &weightedNormalStress,
+        ScalarVector &weights,
         double youngModulus,
         double poissonRatio,
         Vector const &displacement) const {
@@ -159,7 +160,6 @@ void MyAssembler<GridView, dimension>::assembleWeightedNormalStress(
   auto const nodalTractionAverage =
       interpolateP0ToP1(frictionalBoundary, traction);
 
-  ScalarVector weights;
   {
     NeumannBoundaryAssembler<Grid, typename ScalarVector::block_type>
         frictionalBoundaryAssembler(
@@ -168,6 +168,7 @@ void MyAssembler<GridView, dimension>::assembleWeightedNormalStress(
     vertexAssembler.assembleBoundaryFunctional(frictionalBoundaryAssembler,
                                                weights, frictionalBoundary);
   }
+
   auto const normals = frictionalBoundary.getNormals();
   for (size_t i = 0; i < vertexBasis.size(); ++i)
     weightedNormalStress[i] =

dune/tectonic/assemblers.hh

@@ -89,6 +89,7 @@ class MyAssembler {
     void assembleWeightedNormalStress(
             BoundaryPatch<GridView> const &frictionalBoundary,
             ScalarVector &weightedNormalStress,
+            ScalarVector &weights,
             double youngModulus,
             double poissonRatio,
             Vector const &displacement) const;

dune/tectonic/data-structures/program_state.hh

@@ -74,6 +74,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
       v(bodyCount_),
       a(bodyCount_),
       alpha(bodyCount_),
+      weights(bodyCount_),
       weightedNormalStress(bodyCount_) {
     for (size_t i=0; i<bodyCount_; i++) {
       size_t leafVertexCount = leafVertexCounts[i];
@@ -82,6 +83,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
       v[i].resize(leafVertexCount),
       a[i].resize(leafVertexCount),
       alpha[i].resize(leafVertexCount),
+      weights[i].resize(leafVertexCount),
       weightedNormalStress[i].resize(leafVertexCount),
 
       bodyStates[i] = new BodyState(&u[i], &v[i], &a[i], &alpha[i], &weightedNormalStress[i]);
@@ -205,12 +207,14 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
       const auto& body = contactNetwork.body(nonmortarIdx);
 
       ScalarVector frictionBoundaryStress(weightedNormalStress[nonmortarIdx].size());
+      ScalarVector nmWeights(frictionBoundaryStress);
 
       body->assembler()->assembleWeightedNormalStress(
-        contactCoupling->nonmortarBoundary(), frictionBoundaryStress, body->data()->getYoungModulus(),
+        contactCoupling->nonmortarBoundary(), frictionBoundaryStress, nmWeights, body->data()->getYoungModulus(),
         body->data()->getPoissonRatio(), u[nonmortarIdx]);
 
       weightedNormalStress[nonmortarIdx] += frictionBoundaryStress;
+      weights[nonmortarIdx] += nmWeights;
     }
 
     std::cout << "solving linear problem for a..." << std::endl;
@@ -231,6 +235,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
   std::vector<Vector> v;
   std::vector<Vector> a;
   std::vector<ScalarVector> alpha;
+  std::vector<ScalarVector> weights;
   std::vector<ScalarVector> weightedNormalStress;
   double relativeTime;
   double relativeTau;