.

src/data-structures/contactnetwork.cc

@@ -26,7 +26,7 @@ ContactNetwork<HostGridType, VectorType>::ContactNetwork(
     couplings_(nCouplings),
     frictionBoundaries_(nBodies),
     stateEnergyNorms_(nBodies),
-    nBodyAssembler_(nBodies, nCouplings) //, false, 0.0)
+    nBodyAssembler_(nBodies, nCouplings, false, 0.0)
 {}
 
 template <class HostGridType, class VectorType>
@@ -324,6 +324,7 @@ void ContactNetwork<HostGridType, VectorType>::totalNodes(
         std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> boundaryConditions;
         body->boundaryConditions(tag, boundaryConditions);
 
+        if (boundaryConditions.size()>0) {
             const int idxBackup = idx;
             for (size_t bc = 0; bc<boundaryConditions.size(); bc++) {
                 const auto& nodes = boundaryConditions[bc]->boundaryNodes();
@@ -332,6 +333,9 @@ void ContactNetwork<HostGridType, VectorType>::totalNodes(
                         totalNodes[idx][k] = (*nodes)[j][k];
                 idx = (bc==boundaryConditions.size()-1 ? idx : idxBackup);
             }
+        } else {
+            idx += body->nVertices();
+        }
     }
 }
 

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

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

src/multi-body-problem-data/bc.hh

@@ -9,12 +9,12 @@ class VelocityDirichletCondition
     // Assumed to vanish at time zero
       double const finalVelocity = -5e-5;
     
-    std::cout << "VelocityDirichletCondition::evaluate()" << std::endl;
+    /*std::cout << "VelocityDirichletCondition::evaluate()" << std::endl;
     
     if (relativeTime <= 0.1)
         std::cout << "- loading phase" << std::endl;
     else
-        std::cout << "- final velocity reached" << std::endl;
+        std::cout << "- final velocity reached" << std::endl;*/
     
     y = (relativeTime <= 0.1)
             ? finalVelocity * (1.0 - std::cos(relativeTime * M_PI / 0.1)) / 2.0

src/multi-body-problem.cfg

@@ -11,7 +11,7 @@ vtk.write       = true
 
 [problem]
 finalTime       = 100  # [s] #1000
-bodyCount       = 2
+bodyCount       = 4
 
 [body]
 bulkModulus     = 0.5e5 # [Pa]

src/solverfactorytest.cc

@@ -57,7 +57,7 @@
 
 const int dim = MY_DIM;
 
-const int timeSteps = 1;
+const int timeSteps = 5;
 
 Dune::ParameterTree parset;
 size_t bodyCount;
@@ -86,6 +86,11 @@ void solveProblem(const ContactNetwork& contactNetwork,
     //using LocalSolver = LocalBisectionSolver;
     //using Linearization = Linearization<Functional, BitVector>;
 
+    /*print(ignore, "ignore:");
+    for (size_t i=0; i<x.size(); i++) {
+        std::cout << x[i] << std::endl;
+    }*/
+
     // set up reference solver
     Vector refX = x;
     using ContactFunctional = Dune::TNNMG::BoxConstrainedQuadraticFunctional<Matrix&, Vector&, Vector&, Vector&, double>;
@@ -438,8 +443,8 @@ void run(Updaters& updaters, ContactNetwork& contactNetwork,
   nBodyAssembler.concatenateVectors(u0, totalu0);
   nBodyAssembler.concatenateVectors(v0, totalv0);
 
-  updaters.rate_->velocityObstacles(totalu0, lower, totalv0, vLower);
-  updaters.rate_->velocityObstacles(totalu0, upper, totalv0, vUpper);
+  //updaters.rate_->velocityObstacles(totalu0, lower, totalv0, vLower);
+  //updaters.rate_->velocityObstacles(totalu0, upper, totalv0, vUpper);
 
   const auto& errorNorms = contactNetwork.stateEnergyNorms();
 
@@ -465,7 +470,7 @@ void run(Updaters& updaters, ContactNetwork& contactNetwork,
     Vector totalVelocityIterate;
     nBodyAssembler.nodalToTransformed(velocityIterates, totalVelocityIterate);
 
-    solveProblem(contactNetwork, bilinearForm, totalRhs, totalVelocityIterate, totalDirichletNodes, vLower, vUpper);
+    solveProblem(contactNetwork, bilinearForm, totalRhs, totalVelocityIterate, totalDirichletNodes, lower, upper);
 
     nBodyAssembler.postprocess(totalVelocityIterate, velocityIterates);
 
@@ -496,8 +501,8 @@ void run(Updaters& updaters, ContactNetwork& contactNetwork,
         if (alpha[i].size()==0 || newAlpha[i].size()==0)
             continue;
 
-        print(alpha[i], "alpha i:");
-        print(newAlpha[i], "new alpha i:");
+        //print(alpha[i], "alpha i:");
+        //print(newAlpha[i], "new alpha i:");
         if (errorNorms[i]->diff(alpha[i], newAlpha[i]) >= fixedPointTolerance) {
             breakCriterion = false;
             std::cout << "fixedPoint error: " << errorNorms[i]->diff(alpha[i], newAlpha[i]) << std::endl;
@@ -674,10 +679,10 @@ int main(int argc, char *argv[]) { try {
     // Set up TNNMG solver
     // -------------------
 
-    BitVector totalDirichletNodes;
+    /*BitVector totalDirichletNodes;
     contactNetwork.totalNodes("dirichlet", totalDirichletNodes);
 
-    //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>;

src/time-stepping/rate/newmark.cc

@@ -95,20 +95,20 @@ void Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::setup(
             double dirichletValue;
             bodyDirichletFunctions[bc]->evaluate(relativeTime, dirichletValue);
 
-            std::cout << "dirichletValue: " << dirichletValue << std::endl;
+            //std::cout << "dirichletValue: " << dirichletValue << std::endl;
 
                   
             for (size_t k=0; k<bcDirichletNodes.size() ; ++k) {
                 for (size_t j=0; j<dim; ++j) {
                     if (bcDirichletNodes[k][j]) {
-                              std::cout << k << " " << j << std::endl;
+                              //std::cout << k << " " << j << std::endl;
                         bodyIterate[k][j] = dirichletValue;
                     }
                 }
             }
         }
     }
-    print(iterate, "iterate:");
+    //sprint(iterate, "iterate:");
 }
 
 template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>

src/time-stepping/state/ageinglawstateupdater.cc

@@ -66,7 +66,7 @@ template <class ScalarVector, class Vector>
 void AgeingLawStateUpdater<ScalarVector, Vector>::extractAlpha(
         ScalarVector& alpha) {
 
-    std::cout << "alpha size: " << alpha.size() << " nodes_.size() " << nodes_.size() << std::endl;
+    //std::cout << "alpha size: " << alpha.size() << " nodes_.size() " << nodes_.size() << std::endl;
     if (alpha.size() != nodes_.size()) {
         alpha.resize(nodes_.size());
     }