towards recreating 2013 pipping results with rigid foundation

src/foam/foam.cc

@@ -32,6 +32,7 @@
 #include <dune/fufem/geometry/convexpolyhedron.hh>
 #include <dune/fufem/formatstring.hh>
 #include <dune/fufem/hdf5/file.hh>
+#include <dune/fufem/assemblers/transferoperatorassembler.hh>
 
 #include <dune/solvers/norms/energynorm.hh>
 #include <dune/solvers/solvers/loopsolver.hh>
@@ -41,6 +42,8 @@
 #include <dune/contact/common/couplingpair.hh>
 #include <dune/contact/projections/normalprojection.hh>
 
+#include <dune/tnnmg/functionals/quadraticfunctional.hh>
+
 
 #include <dune/tectonic/assemblers.hh>
 #include <dune/tectonic/gridselector.hh>
@@ -196,173 +199,173 @@ int main(int argc, char *argv[]) {
       programState.setupInitialConditions(parset, contactNetwork);
     }
 
-    /*
-    // setup initial conditions in program state
-    programState.alpha[0] = 0;
-    programState.alpha[1] = parset.get<double>("boundary.friction.initialAlpha"); */
-
-    /*auto myGlobalNonlinearity = myAssembler.assembleFrictionNonlinearity(
-          frictionalBoundary, frictionData);
-      myGlobalNonlinearity->updateLogState(alpha_initial);
-*/
 
-    /*
-    // Solving a linear problem with a multigrid solver
-    auto linearSolver = makeLinearSolver<ContactNetwork, Vector>(parset, contactNetwork);
     const auto& nBodyAssembler = contactNetwork.nBodyAssembler();
+   /* auto linearSolver = makeLinearSolver<ContactNetwork, Vector>(parset, contactNetwork);
+    auto nonlinearity = ZeroNonlinearity();
 
+    // Solving a linear problem with a multigrid solver
     auto const solveLinearProblem = [&](
         const BitVector& _dirichletNodes, const std::vector<std::shared_ptr<Matrix>>& _matrices,
         const std::vector<Vector>& _rhs, std::vector<Vector>& _x,
         Dune::ParameterTree const &_localParset) {
 
-      std::vector<const Matrix*> matrices_ptr(_matrices.size());
-      for (size_t i=0; i<matrices_ptr.size(); i++) {
-            matrices_ptr[i] = _matrices[i].get();
-      }
+        Vector totalX;
+        nBodyAssembler.nodalToTransformed(_x, totalX);
 
-      // assemble full global contact problem
-      Matrix bilinearForm;
-      nBodyAssembler.assembleJacobian(matrices_ptr, bilinearForm);
-
-      Vector totalRhs, oldTotalRhs;
-      nBodyAssembler.assembleRightHandSide(_rhs, totalRhs);
-      oldTotalRhs = totalRhs;
-
-      Vector totalX, oldTotalX;
-      nBodyAssembler.nodalToTransformed(_x, totalX);
-      oldTotalX = totalX;
-
-      // get lower and upper obstacles
-      const auto& totalObstacles = nBodyAssembler.totalObstacles_;
-      Vector lower(totalObstacles.size());
-      Vector upper(totalObstacles.size());
-
-      for (size_t j=0; j<totalObstacles.size(); ++j) {
-          const auto& totalObstaclesj = totalObstacles[j];
-          auto& lowerj = lower[j];
-          auto& upperj = upper[j];
-        for (size_t d=0; d<dims; ++d) {
-            lowerj[d] = totalObstaclesj[d][0];
-            upperj[d] = totalObstaclesj[d][1];
-        }
-      }
+        FunctionalFactory<std::decay_t<decltype(nBodyAssembler)>, decltype(nonlinearity), Matrix, Vector> functionalFactory(nBodyAssembler, nonlinearity, _matrices, _rhs);
+        functionalFactory.build();
+        auto functional = functionalFactory.functional();
 
-      // set up functional
-      using Functional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, typename Matrix::field_type>;
-      Functional J(bilinearForm, totalRhs, ZeroNonlinearity(), lower, upper);
+        NonlinearSolver<std::remove_reference_t<decltype(*functional)>, BitVector> solver(parset.sub("solver.tnnmg"), linearSolver, functional, _dirichletNodes);
+        solver.solve(_localParset, totalX);
 
-      // set up TNMMG solver
-      using Factory = SolverFactory<Functional, BitVector>;
-      //Factory factory(parset.sub("solver.tnnmg"), J, linearSolver, _dirichletNodes);
-      Factory factory(parset.sub("solver.tnnmg"), J, _dirichletNodes);
-      factory.build(linearSolver); */
+        nBodyAssembler.postprocess(totalX, _x);
+    };
 
-     /* std::vector<BitVector> bodyDirichletNodes;
-      nBodyAssembler.postprocess(_dirichletNodes, bodyDirichletNodes);
-      for (size_t i=0; i<bodyDirichletNodes.size(); i++) {
-        print(bodyDirichletNodes[i], "bodyDirichletNodes_" + std::to_string(i) + ": ");
-      }*/
+    programState.timeStep = parset.get<size_t>("initialTime.timeStep");
+    programState.relativeTime = parset.get<double>("initialTime.relativeTime");
+    programState.relativeTau = parset.get<double>("initialTime.relativeTau");
 
-    /*
-      auto tnnmgStep = factory.step();
-      factory.setProblem(totalX);
+    std::vector<Vector> ell0(bodyCount);
+    for (size_t i=0; i<bodyCount; i++) {
+      programState.u[i] = 0.0;
+      programState.v[i] = 0.0;
+      programState.a[i] = 0.0;
 
-      const EnergyNorm<Matrix, Vector> norm(bilinearForm);
+      ell0[i].resize(programState.u[i].size());
+      ell0[i] = 0.0;
 
-      LoopSolver<Vector> solver(
-          *tnnmgStep.get(), _localParset.get<size_t>("maximumIterations"),
-          _localParset.get<double>("tolerance"), norm,
-          _localParset.get<Solver::VerbosityMode>("verbosity")); // absolute error
+      contactNetwork.body(i)->externalForce()(programState.relativeTime, ell0[i]);
+    }
+    */
+    // Initial displacement: Start from a situation of minimal stress,
+    // which is automatically attained in the case [v = 0 = a].
+    // Assuming dPhi(v = 0) = 0, we thus only have to solve Au = ell0
+    BitVector totalDirichletNodes;
+    contactNetwork.totalNodes("dirichlet", totalDirichletNodes);
+    size_t dof = 0;
+    for (size_t i=0; i<bodyCount; i++) {
+        const auto& body = *contactNetwork.body(i);
+
+        if (body.data()->getYoungModulus() == 0.0) {
+            for (; dof<body.nVertices(); dof++) {
+                totalDirichletNodes[dof] = true;
+            }
+        } else {
+            dof += body.nVertices();
+        }
+    }
 
-      solver.preprocess();
-      solver.solve();
+    using BoundaryNodes = typename ContactNetwork::BoundaryNodes;
+    BoundaryNodes dirichletNodes;
+    contactNetwork.boundaryNodes("dirichlet", dirichletNodes);
 
-      nBodyAssembler.postprocess(tnnmgStep->getSol(), _x);
-    }; */
+    std::vector<const Dune::BitSetVector<1>*> frictionNodes;
+    contactNetwork.frictionNodes(frictionNodes);
 
-     /*using LinearFactory = SolverFactory<
-          dims, BlockNonlinearTNNMGProblem<ConvexProblem<
-                    ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
-          Grid>;
-      ZeroNonlinearity<LocalVector, LocalMatrix> zeroNonlinearity;*/
+    /*
+    std::cout << "solving linear problem for u..." << std::endl;
+
+    {
+        int nVertices = contactNetwork.body(1)->nVertices();
+        BitVector uDirichletNodes(nVertices, false);
+
+        int idx=0;
+        const auto& body = contactNetwork.body(1);
+        using LeafBody = typename ContactNetwork::LeafBody;
+        std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> boundaryConditions;
+        body->boundaryConditions("dirichlet", boundaryConditions);
+
+        if (boundaryConditions.size()>0) {
+            const int idxBackup = idx;
+            for (size_t bc = 0; bc<boundaryConditions.size(); bc++) {
+                const auto& nodes = boundaryConditions[bc]->boundaryNodes();
+                for (size_t j=0; j<nodes->size(); j++, idx++)
+                    for (int k=0; k<dims; k++)
+                        uDirichletNodes[idx][k] = uDirichletNodes[idx][k] or (*nodes)[j][k];
+                idx = (bc==boundaryConditions.size()-1 ? idx : idxBackup);
+            }
+        }
 
-  /*
-    // TODO: clean up once generic lambdas arrive
-      auto const solveLinearProblem = [&](
-          Dune::BitSetVector<dims> const &_dirichletNodes, Matrix const &_matrix,
-          Vector const &_rhs, Vector &_x, EnergyNorm<Matrix, Vector> _norm,
-          Dune::ParameterTree const &_localParset) {
+        for (auto i=0; i<nVertices; i++) {
+            if ((*frictionNodes[1])[i][0]) {
+                uDirichletNodes[i][1] = true;
+            }
+        }
 
+        using Norm = EnergyNorm<Matrix, Vector>;
+        using LinearSolver = typename Dune::Solvers::LoopSolver<Vector>;
+
+        // transfer operators need to be recomputed on change due to setDeformation()
+        using TransferOperator = CompressedMultigridTransfer<Vector>;
+        using TransferOperators = std::vector<std::shared_ptr<TransferOperator>>;
+
+        const auto& grid = contactNetwork.body(1)->grid();
+        TransferOperators transfer(grid->maxLevel());
+        for (size_t i = 0; i < transfer.size(); ++i)
+        {
+            // create transfer operators from level i to i+1 (note that this will only work for either uniformly refined grids or adaptive grids with RefinementType=COPY)
+            auto t = std::make_shared<TransferOperator>();
+            t->setup(*grid, i, i+1);
+            transfer[i] = t;
+        }
 
+        auto linearMultigridStep = std::make_shared<Dune::Solvers::MultigridStep<Matrix, Vector> >(
+                    *contactNetwork.matrices().elasticity[1],
+                    programState.u[1],
+                    ell0[1]);
+        linearMultigridStep->setIgnore(uDirichletNodes);
+        linearMultigridStep->setMGType(1, 3, 3);
+        linearMultigridStep->setSmoother(TruncatedBlockGSStep<Matrix, Vector>());
+        linearMultigridStep->setTransferOperators(transfer);
 
-        typename LinearFactory::ConvexProblem convexProblem(
-            1.0, _matrix, zeroNonlinearity, _rhs, _x);
-        typename LinearFactory::BlockProblem problem(parset, convexProblem);
+        const auto& solverParset = parset.sub("u0.solver");
+        Dune::Solvers::LoopSolver<Vector> solver(linearMultigridStep, solverParset.get<size_t>("maximumIterations"),
+            solverParset.get<double>("tolerance"), Norm(*linearMultigridStep),
+            solverParset.get<Solver::VerbosityMode>("verbosity")); // absolute error
 
-        auto multigridStep = factory.getSolver();
-        multigridStep->setProblem(_x, problem);
-        LoopSolver<Vector> solver(
-            multigridStep, _localParset.get<size_t>("maximumIterations"),
-            _localParset.get<double>("tolerance"), &_norm,
-            _localParset.get<Solver::VerbosityMode>("verbosity"),
-            false); // absolute error
 
         solver.preprocess();
         solver.solve();
-      };
-
-      // Solve the stationary problem
-      programState.u[0] = 0.0;
-      programState.u[1] = 0.0;
-
-      solveLinearProblem(dirichletNodes, contactNetwork.matrices().elasticity, ell0, programState.u,
-                         parset.sub("u0.solver"));
-
-      programState.v[0] = 0.0;
-      programState.v[1] = 0.0;
-
-
-      programState.a[0] = 0.0;
-      programState.a[1] = 0.0;
-      {
-          // Initial acceleration: Computed in agreement with Ma = ell0 - Au
-          // (without Dirichlet constraints), again assuming dPhi(v = 0) = 0
-          std::vector<Vector> accelerationRHS = ell0;
-
-          for (size_t i=0; i<bodyCount; i++) {
-            const auto& body = contactNetwork.body(i);
-            Dune::MatrixVector::subtractProduct(accelerationRHS[i], *body->matrices().elasticity, u[i]);
-          }
-
-          BitVector noNodes(dirichletNodes.size(), false);
-          solveLinearProblem(noNodes, contactNetwork.matrices().mass, accelerationRHS, programState.a,
-                             parset.sub("a0.solver"));
+    }
 
-      }
+    //print(u, "initial u:");
 
-      for (size_t i=0; i<contactNetwork.nCouplings(); i++) {
-        const auto& coupling = contactNetwork.coupling(i);
-        const auto& contactCoupling = contactNetwork.nBodyAssembler().getContactCouplings()[i];
+    // Initial acceleration: Computed in agreement with Ma = ell0 - Au
+    // (without Dirichlet constraints), again assuming dPhi(v = 0) = 0
+    std::vector<Vector> accelerationRHS = ell0;
+    for (size_t i=0; i<bodyCount; i++) {
+      const auto& body = contactNetwork.body(i);
+      Dune::MatrixVector::subtractProduct(accelerationRHS[i], *body->matrices().elasticity, programState.u[i]);
+    }
 
-        const auto nonmortarIdx = coupling->gridIdx_[0];
-        const auto& body = contactNetwork.body(nonmortarIdx);
+    std::cout << "solving linear problem for a..." << std::endl;
 
-        ScalarVector frictionBoundaryStress(weightedNormalStress[nonmortarIdx].size());
+    BitVector noNodes(totalDirichletNodes.size(), false);
+    solveLinearProblem(noNodes, contactNetwork.matrices().mass, accelerationRHS, programState.a,
+                       parset.sub("a0.solver"));
 
-        body->assembler()->assembleWeightedNormalStress(
-          contactCoupling->nonmortarBoundary(), frictionBoundaryStress, body->data()->getYoungModulus(),
-          body->data()->getPoissonRatio(), u[nonmortarIdx]);
+    // setup initial conditions in program state
+    programState.alpha[0] = 0;
+    programState.alpha[1] = parset.get<double>("boundary.friction.initialAlpha");
 
-        weightedNormalStress[nonmortarIdx] += frictionBoundaryStress;
-      }
+    for (size_t i=0; i<contactNetwork.nCouplings(); i++) {
+      const auto& coupling = contactNetwork.coupling(i);
+      const auto& contactCoupling = contactNetwork.nBodyAssembler().getContactCouplings()[i];
 
+      const auto nonmortarIdx = coupling->gridIdx_[0];
+      const auto& body = contactNetwork.body(nonmortarIdx);
 
-*/
+      ScalarVector frictionBoundaryStress(programState.weightedNormalStress[nonmortarIdx].size());
 
+      body->assembler()->assembleWeightedNormalStress(
+        contactCoupling->nonmortarBoundary(), frictionBoundaryStress, body->data()->getYoungModulus(),
+        body->data()->getPoissonRatio(), programState.u[nonmortarIdx]);
 
+      programState.weightedNormalStress[nonmortarIdx] += frictionBoundaryStress;
+    } */
 
-    auto& nBodyAssembler = contactNetwork.nBodyAssembler();
     for (size_t i=0; i<bodyCount; i++) {
       contactNetwork.body(i)->setDeformation(programState.u[i]);
     }
@@ -387,46 +390,23 @@ int main(int argc, char *argv[]) {
     // Set up TNNMG solver
     // -------------------
 
-    BitVector totalDirichletNodes;
-    contactNetwork.totalNodes("dirichlet", totalDirichletNodes);
-
-    /*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];
-        }
-
-        totalDirichletNodes[i] = val;
-        for (size_t d=0; d<dims; d++) {
-            totalDirichletNodes[i][d] = val;
-        }
-    }*/
-
-    //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>;
     using NonlinearFactory = SolverFactory<Functional, BitVector>;
 
     using BoundaryFunctions = typename ContactNetwork::BoundaryFunctions;
-    using BoundaryNodes = typename ContactNetwork::BoundaryNodes;
     using Updaters = Updaters<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>,
                                StateUpdater<ScalarVector, Vector>>;
 
     BoundaryFunctions velocityDirichletFunctions;
     contactNetwork.boundaryFunctions("dirichlet", velocityDirichletFunctions);
 
-    BoundaryNodes dirichletNodes;
-    contactNetwork.boundaryNodes("dirichlet", dirichletNodes);
-
     /*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++) {
         print(*frictionNodes[i], "frictionNodes_body_" + std::to_string(i));
@@ -491,7 +471,7 @@ int main(int argc, char *argv[]) {
         timeStepper(stepBase, contactNetwork, current,
                             programState.relativeTime, programState.relativeTau);
 
-    size_t timeSteps = parset.get<size_t>("timeSteps.timeSteps");
+    size_t timeSteps = std::round(parset.get<double>("timeSteps.timeSteps"));
 
     while (!timeStepper.reachedEnd()) {
       programState.timeStep++;