.

dune/tectonic/spatial-solving/contact/dualmortarcoupling.cc

@@ -149,7 +149,10 @@ void DualMortarCoupling<field_type, GridType0, GridType1>::assembleNonmortarLagr
             // we can use subEntity here because we add all indices anyway
             int v = globalToLocal_[indexSet.subIndex(inside,refElement.subEntity(nIt.indexInInside(),1,
                                                                           i,dim),dim)];
-                Dune::MatrixVector::addToDiagonal(nonmortarLagrangeMatrix_[v][v],sfI);
+            if (v == -1)
+               continue;
+
+            Dune::MatrixVector::addToDiagonal(nonmortarLagrangeMatrix_[v][v],sfI);
         }
 
     }
@@ -178,9 +181,13 @@ void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
     const auto& vertices = *nonmortarBoundary_.getVertices();
 
     globalToLocal_.resize(vertices.size());
+    hasObstacle_.resize(vertices.size());
     int idx = 0;
-    for (size_t i=0; i<globalToLocal_.size(); i++)
-        globalToLocal_[i] = (vertices[i][0] and crosspoints_.count(i)==0) ? idx++ : -1;
+    for (size_t i=0; i<globalToLocal_.size(); i++) {
+        bool nonmortarVertex = vertices[i][0] and crosspoints_.count(i)==0;
+        globalToLocal_[i] = nonmortarVertex ? idx++ : -1;
+        hasObstacle_[i][0] = nonmortarVertex ? true : false;
+    }
 
     // Assemble the diagonal matrix coupling the nonmortar side with the lagrange multiplyers there
     assembleNonmortarLagrangeMatrix();
@@ -249,8 +256,6 @@ void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
 
     // loop over all intersections and compute the matrix entries
     for (const auto& rIs : intersections(*glue_)) {
-
-
         const auto& inside = rIs.inside();
 
         if (!nonmortarBoundary_.contains(rIs.inside(),rIs.indexInInside()))
@@ -291,9 +296,8 @@ void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
         const auto& rGeomInOutside = rIs.geometryInOutside();
 
         bool isCrosspointFace = false;
-        int nNonmortarFaceNodes = domainRefElement.size(rIs.indexInInside(),1,dim);
         std::vector<int> nonmortarFaceNodes;
-        for (int i=0; i<nNonmortarFaceNodes; i++) {
+        for (int i=0; i<domainRefElement.size(rIs.indexInInside(),1,dim); i++) {
           int faceIdxi = domainRefElement.subEntity(rIs.indexInInside(), 1, i, dim);
           int globalIdx = indexSet0.subIndex(inside,faceIdxi,dim);
 
@@ -326,35 +330,59 @@ void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
 
             if (isCrosspointFace) {
               constantDualFE.localBasis().evaluateFunction(nonmortarQuadPos,dualQuadValues);
+
+              // loop over all Lagrange multiplier shape functions
+              for (size_t j=0; j<nonmortarFaceNodes.size(); j++) {
+
+                  int globalDomainIdx = indexSet0.subIndex(inside,nonmortarFaceNodes[j],dim);
+                  int rowIdx = globalToLocal_[globalDomainIdx];
+
+                  weakObstacle_[rowIdx][0] += integrationElement * quadPt.weight()
+                      * dualQuadValues[0] * (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[0]* mortarQuadValues[k];
+
+                      Dune::MatrixVector::addToDiagonal(mortarLagrangeMatrix_[rowIdx][colIdx], mortarEntry);
+
+                  }
+
+               }
             } else {
               dualCache->evaluateFunction(nonmortarQuadPos,dualQuadValues);
-            }
 
-            // loop over all Lagrange multiplier shape functions
-            for (int j=0; j<nNonmortarFaceNodes; j++) {
+              // loop over all Lagrange multiplier shape functions
+              for (size_t j=0; j<nonmortarFaceNodes.size(); j++) {
 
-                int globalDomainIdx = indexSet0.subIndex(inside,nonmortarFaceNodes[j],dim);
-                int rowIdx = globalToLocal_[globalDomainIdx];
+                  int globalDomainIdx = indexSet0.subIndex(inside,nonmortarFaceNodes[j],dim);
+                  int rowIdx = globalToLocal_[globalDomainIdx];
 
-                weakObstacle_[rowIdx][0] += integrationElement * quadPt.weight()
-                    * dualQuadValues[nonmortarFaceNodes[j]] * (gapVector*avNormals[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++) {
+                  // 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;
+                      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];
+                      // 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);
+                      Dune::MatrixVector::addToDiagonal(mortarLagrangeMatrix_[rowIdx][colIdx], mortarEntry);
 
-                }
+                  }
 
+               }
             }
-
         }
 
     }

dune/tectonic/spatial-solving/contact/dualmortarcoupling.hh

@@ -144,6 +144,10 @@ class DualMortarCoupling {
         return weakObstacle_;
     }
 
+    const Dune::BitSetVector<1>* hasObstacle() const {
+        return &hasObstacle_;
+    }
+
     const std::vector<int>& globalToLocal() const {
         return globalToLocal_;
     }
@@ -181,6 +185,8 @@ class DualMortarCoupling {
     /** \brief The weak obstacles. */
     ObstacleVector weakObstacle_;
 
+    Dune::BitSetVector<1> hasObstacle_;
+
     /** \brief Allow some small penetration of the bodies. */
     field_type overlap_;
 

dune/tectonic/spatial-solving/contact/nbodyassembler.cc

@@ -15,7 +15,7 @@
 
 template <class GridType, class VectorType>
 void NBodyAssembler<GridType, VectorType>::setCrosspoints() {
-    std::vector<Dune::BitSetVector<1>> bodywiseBoundaries(nGrids());
+   /* std::vector<Dune::BitSetVector<1>> bodywiseBoundaries(nGrids());
     for (size_t i=0; i<nGrids(); i++) {
         bodywiseBoundaries[i] = *dirichletVertices_[i];
     }
@@ -36,6 +36,29 @@ void NBodyAssembler<GridType, VectorType>::setCrosspoints() {
             bodywiseBoundary[j][0] = bodywiseBoundary[j][0] or nmBoundaryVertices[j][0];
         }
 
+    }*/
+
+    for (int i=0; i<nCouplings(); i++) {
+        auto& coupling = coupling_[i];
+        const auto& contactCoupling = contactCoupling_[i]; // dual mortar object holding boundary patches
+        const auto nonmortarIdx = coupling.gridIdx_[0];
+
+        auto& crosspoints = crosspoints_[nonmortarIdx];
+
+        Dune::BitSetVector<1> nonmortarPatchBoundary;
+        contactCoupling->nonmortarBoundary().getPatchBoundaryVertices(nonmortarPatchBoundary);
+
+        for (size_t j=0; j<nonmortarPatchBoundary.size(); j++) {
+            if (nonmortarPatchBoundary[j][0]) {
+                crosspoints.emplace(j);
+            }
+        }
+
+    }
+
+
+    for (size_t i=0; i<crosspoints_.size(); i++) {
+        print(crosspoints_[i], "crosspoints " + std::to_string(i));
     }
 }
 
@@ -73,13 +96,16 @@ void NBodyAssembler<GridType, VectorType>::assembleTransferOperator()
 
     setCrosspoints();
 
-    for (size_t i=0; i<nGrids(); i++) {
+    for (int i=0; i<nGrids(); i++) {
         print(crosspoints_[i], "crosspoints " + std::to_string(i));
     }
 
     for (size_t i=0; i<contactCoupling_.size(); i++) {
         contactCoupling_[i]->setCrosspoints(crosspoints_[coupling_[i].gridIdx_[0]]);
         contactCoupling_[i]->setup();
+
+        print(contactCoupling_[i]->nonmortarLagrangeMatrix() , "nonmortarLagrangeMatrix " + std::to_string(i) + ": ");
+                print(contactCoupling_[i]->mortarLagrangeMatrix() , "mortarLagrangeMatrix " + std::to_string(i) + ": ");
     }
 
     // setup Householder reflections
@@ -135,14 +161,15 @@ void NBodyAssembler<GridType, VectorType>::assembleHouseholderReflections()
         for (int k=0;k<grid0Idx;k++)
             idx += grids_[k]->size(dim);
 
-        const auto& nonmortarBoundary = contactCoupling_[j]->nonmortarBoundary();
         const auto& globalToLocal = contactCoupling_[j]->globalToLocal();
 
         // if a body is non-mortar more than once, one has to be careful with not overwriting entries
         for (std::size_t k=0; k<coordinateSystems[j].size(); k++)
-          if(globalToLocal[k] > -1)
+          if(globalToLocal[k] > -1 ) //or crosspoints.count(k) > 0)
             this->localCoordSystems_[idx+k] = coordinateSystems[j][k];
     }
+
+    //print(this->localCoordSystems_, "localCoordSystems:");
 }
 
 template <class GridType, class VectorType>
@@ -162,7 +189,6 @@ void NBodyAssembler<GridType, VectorType>::assembleObstacle()
     for (int j=0; j<nCouplings(); j++) {
 
         int grid0Idx = coupling_[j].gridIdx_[0];
-        const auto& nonmortarBoundary = contactCoupling_[j]->nonmortarBoundary();
         const auto& globalToLocal = globalToLocals[j];
 
         // compute offset
@@ -243,6 +269,11 @@ void NBodyAssembler<GridType, VectorType>::assembleObstacle()
         }
     }
 
+    //totalHasObstacle_[35] = true;
+    //totalObstacles_[35].lower(0) = 0;
+   // totalObstacles_[35].lower(1) = 0;
+    //totalObstacles_[35].upper(0) = 0;
+   // totalObstacles_[35].upper(1) = 0;
 }
 
 template <class GridType, class VectorType>

dune/tectonic/spatial-solving/contact/nbodycontacttransfer.cc

@@ -81,7 +81,7 @@ void NBodyContactTransfer<ContactNetwork, VectorType>::setup(const ContactNetwor
 
         for (size_t i=0; i<nCouplings; i++) {
             mortarTransferOperators[i] = &contactCouplings[i]->mortarLagrangeMatrix();
-            fineHasObstacle[i] = contactCouplings[i]->nonmortarBoundary().getVertices();
+            fineHasObstacle[i] = contactCouplings[i]->hasObstacle();
             gridIdx[i] = nBodyAssembler.getCoupling(i).gridIdx_;
         }
 

src/foam/foam-2D.cfg

 # -*- mode:conf -*-
 [body0]
-smallestDiameter = 0.01  # 2e-3 [m]
+smallestDiameter = 0.02  # 2e-3 [m]
 
 [body1]
-smallestDiameter = 0.01  # 2e-3 [m]
+smallestDiameter = 0.02  # 2e-3 [m]
 
 [timeSteps]
 refinementTolerance = 1e-5 # 1e-5

src/foam/foam.cc

@@ -198,9 +198,6 @@ int main(int argc, char *argv[]) {
     else
      programState.setupInitialConditions(parset, contactNetwork);
 
-
-    //DUNE_THROW(Dune::Exception, "Just need to stop here!");
-
     auto& nBodyAssembler = contactNetwork.nBodyAssembler();
     for (size_t i=0; i<bodyCount; i++) {
       contactNetwork.body(i)->setDeformation(programState.u[i]);
@@ -343,6 +340,8 @@ int main(int argc, char *argv[]) {
         print(*frictionNodes[i], "frictionNodes_body_" + std::to_string(i));
     }
 
+    DUNE_THROW(Dune::Exception, "Just need to stop here!");
+
     Updaters current(
         initRateUpdater(
             parset.get<Config::scheme>("timeSteps.scheme"),

src/foam/foam.cfg

@@ -2,19 +2,19 @@
 gravity         = 9.81     # [m/s^2]
 
 [body0]
-length          = 0.04      # [m]
+length          = 0.08      # [m]
 height          = 0.04     # [m]
 depth           = 0.04     # [m]
 bulkModulus     = 0.5e5    # [Pa] #2190
 poissonRatio    = 0.30     # [1]  #0.11
 [body0.elastic]
 density         = 900      # [kg/m^3] #750
-shearViscosity  = 1e3      # [Pas]
-bulkViscosity   = 1e3      # [Pas]
+shearViscosity  = 1e20     # [Pas]
+bulkViscosity   = 1e20     # [Pas]
 [body0.viscoelastic]
 density         = 1000     # [kg/m^3]
-shearViscosity  = 1e4      # [Pas]
-bulkViscosity   = 1e4      # [Pas]
+shearViscosity  = 1e20     # [Pas]
+bulkViscosity   = 1e20     # [Pas]
 
 [body1]
 length          = 0.04     # [m]
@@ -24,12 +24,12 @@ bulkModulus     = 0.5e5    # [Pa]
 poissonRatio    = 0.30     # [1]
 [body1.elastic]
 density         = 900      # [kg/m^3]
-shearViscosity  = 1e3      # [Pas]
-bulkViscosity   = 1e3      # [Pas]
+shearViscosity  = 1e20     # [Pas]
+bulkViscosity   = 1e20     # [Pas]
 [body1.viscoelastic]
 density         = 1000     # [kg/m^3]
-shearViscosity  = 1e4      # [Pas]
-bulkViscosity   = 1e4      # [Pas]
+shearViscosity  = 1e20     # [Pas]
+bulkViscosity   = 1e20     # [Pas]
 
 [boundary.friction]
 C               = 10       # [Pa]

src/multi-body-problem/multi-body-problem-2D.cfg

 # -*- mode:conf -*-
 [boundary.friction]
-smallestDiameter = 0.05  # 2e-3 [m]
+smallestDiameter = 0.1  # 2e-3 [m]
 
 [timeSteps]
 refinementTolerance = 1e-5 # 1e-5

src/multi-body-problem/multi-body-problem.cc

@@ -137,20 +137,21 @@ int main(int argc, char *argv[]) {
     const size_t bodyCount = contactNetwork.nBodies();
 
     for (size_t i=0; i<contactNetwork.nLevels(); i++) {
-        // printDofLocation(contactNetwork.body(i)->gridView());
+        //printDofLocation(contactNetwork.body(i)->gridView());
 
         //Vector def(contactNetwork.deformedGrids()[i]->size(dims));
         //def = 1;
         //deformedGridComplex.setDeformation(def, i);
 
-        const auto& level = *contactNetwork.level(i);
+        /*const auto& level = *contactNetwork.level(i);
 
         for (size_t j=0; j<level.nBodies(); j++) {
             writeToVTK(level.body(j)->gridView(), "../debug_print/bodies/", "body_" + std::to_string(j) + "_level_" + std::to_string(i));
-        }
+        }*/
     }
 
     for (size_t i=0; i<bodyCount; i++) {
+        printDofLocation(contactNetwork.body(i)->gridView());
         writeToVTK(contactNetwork.body(i)->gridView(), "../debug_print/bodies/", "body_" + std::to_string(i) + "_leaf");
     }
 
@@ -199,8 +200,7 @@ int main(int argc, char *argv[]) {
     else
      programState.setupInitialConditions(parset, contactNetwork);
 
-
-    //DUNE_THROW(Dune::Exception, "Just need to stop here!");
+    //print(programState.u, "u:");
 
     auto& nBodyAssembler = contactNetwork.nBodyAssembler();
     for (size_t i=0; i<bodyCount; i++) {
@@ -296,6 +296,8 @@ int main(int argc, char *argv[]) {
     };
     report(true);
 
+    //DUNE_THROW(Dune::Exception, "Just need to stop here!");
+
     // -------------------
     // Set up TNNMG solver
     // -------------------
@@ -303,7 +305,7 @@ int main(int argc, char *argv[]) {
     BitVector totalDirichletNodes;
     contactNetwork.totalNodes("dirichlet", totalDirichletNodes);
 
-    /*for (size_t i=0; i<totalDirichletNodes.size(); i++) {
+    for (size_t i=0; i<totalDirichletNodes.size(); i++) {
         bool val = false;
         for (size_t d=0; d<dims; d++) {
             val = val || totalDirichletNodes[i][d];
@@ -313,9 +315,9 @@ int main(int argc, char *argv[]) {
         for (size_t d=0; d<dims; d++) {
             totalDirichletNodes[i][d] = val;
         }
-    }*/
+    }
 
-    //print(totalDirichletNodes, "totalDirichletNodes:");
+    print(totalDirichletNodes, "totalDirichletNodes:");
 
     //using Functional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, field_type>;
     using Functional = Functional<Matrix&, Vector&, GlobalFriction<Matrix, Vector>&, Vector&, Vector&, field_type>;
@@ -332,18 +334,18 @@ int main(int argc, char *argv[]) {
     BoundaryNodes dirichletNodes;
     contactNetwork.boundaryNodes("dirichlet", dirichletNodes);
 
-    /*for (size_t i=0; i<dirichletNodes.size(); i++) {
+    for (size_t i=0; i<dirichletNodes.size(); i++) {
         for (size_t j=0; j<dirichletNodes[i].size(); j++) {
         print(*dirichletNodes[i][j], "dirichletNodes_body_" + std::to_string(i) + "_boundary_" + std::to_string(j));
         }
-    }*/
+    }
 
     std::vector<const Dune::BitSetVector<1>*> frictionNodes;
     contactNetwork.frictionNodes(frictionNodes);
 
-    /*for (size_t i=0; i<frictionNodes.size(); i++) {
+    for (size_t i=0; i<frictionNodes.size(); i++) {
         print(*frictionNodes[i], "frictionNodes_body_" + std::to_string(i));
-    }*/
+    }
 
     Updaters current(
         initRateUpdater(