.

src/multi-body-problem.cc

@@ -81,20 +81,6 @@
 
 size_t const dims = MY_DIM;
 
-template <class GridView>
-void writeToVTK(const GridView& gridView, const std::string path, const std::string name) {
-    Dune::VTKWriter<GridView> vtkwriter(gridView);
-
-    //std::ofstream lStream( "garbage.txt" );
-    std::streambuf* lBufferOld = std::cout.rdbuf();
-    //std::cout.rdbuf(lStream.rdbuf());
-
-    vtkwriter.pwrite(name, path, path);
-
-    std::cout.rdbuf( lBufferOld );
-}
-
-
 Dune::ParameterTree getParameters(int argc, char *argv[]) {
   Dune::ParameterTree parset;
   Dune::ParameterTreeParser::readINITree("/home/mi/podlesny/software/dune/dune-tectonic/src/multi-body-problem.cfg", parset);
@@ -160,7 +146,9 @@ int main(int argc, char *argv[]) {
     }
 #endif
 
+    // ------------
     // set up grids
+    // ------------
     GridsConstructor<Grid> gridConstructor(cuboidGeometries);
     auto& grids = gridConstructor.getGrids();
 
@@ -174,8 +162,6 @@ int main(int argc, char *argv[]) {
         weakPatch = getWeakPatch<LocalVector>(parset.sub("boundary.friction.weakening"), cuboidGeometry);
         refine(*grids[i], weakPatch, parset.get<double>("boundary.friction.smallestDiameter"), cuboidGeometry.lengthScale);
 
-        writeToVTK(grids[i]->leafGridView(), "", "grid" + std::to_string(i));
-
         // determine minDiameter and maxDiameter
         double minDiameter = std::numeric_limits<double>::infinity();
         double maxDiameter = 0.0;
@@ -197,18 +183,9 @@ int main(int argc, char *argv[]) {
 
     const auto deformedGrids = deformedGridComplex.gridVector();
 
-    /*
-    // test deformed grids
-    for (size_t i=0; i<deformedGrids.size(); i++) {
-      Vector def(deformedGrids[i]->size(dims));
-      def = 1;
-      deformedGridComplex.setDeformation(def, i);
-
-      writeToVTK(deformedGrids[i]->leafGridView(), "", "deformedGrid" + std::to_string(i));
-    }
-    // test deformed grids
-    */
-
+    // ------------------------
+    // set up contact couplings
+    // ------------------------
     std::vector<std::shared_ptr<DefLeafGridView>> leafViews(deformedGrids.size());
     std::vector<std::shared_ptr<DefLevelGridView>> levelViews(deformedGrids.size());
     std::vector<size_t> leafVertexCounts(deformedGrids.size());
@@ -247,12 +224,9 @@ int main(int argc, char *argv[]) {
       couplings[i]->set(i, i+1, nonmortarPatch, mortarPatch, 0.1, Dune::Contact::CouplingPairBase::STICK_SLIP, contactProjection, backend);
     }
 
+    // ----------------------------------
     // set up dune-contact nBodyAssembler
-    /*std::vector<const DeformedGrid*> const_grids(bodyCount);
-    for (size_t i=0; i<const_grids.size(); i++) {
-        const_grids[i] = &deformedGrids.grid("body" + std::to_string(i));
-    }*/
-
+    // ----------------------------------
     Dune::Contact::NBodyAssembler<DeformedGrid, Vector> nBodyAssembler(bodyCount, bodyCount-1);
     nBodyAssembler.setGrids(deformedGrids);
 
@@ -263,7 +237,9 @@ int main(int argc, char *argv[]) {
     nBodyAssembler.assembleTransferOperator();
     nBodyAssembler.assembleObstacle();
 
+    // --------------------------
     // set up boundary conditions
+    // --------------------------
     std::vector<BoundaryPatch<DefLeafGridView>> neumannBoundaries(bodyCount);
     std::vector<BoundaryPatch<DefLeafGridView>> surfaces(bodyCount);
 
@@ -273,7 +249,7 @@ int main(int argc, char *argv[]) {
 
     // Dirichlet boundary
     std::vector<Dune::BitSetVector<dims>> noNodes(bodyCount);
-    std::vector<Dune::BitSetVector<dims>> dirichletNodes(bodyCount);
+    std::vector<std::vector<Dune::BitSetVector<dims>>> dirichletNodes(bodyCount);
 
     // set up functions for time-dependent boundary conditions
     using Function = Dune::VirtualFunction<double, double>;
@@ -306,8 +282,10 @@ int main(int argc, char *argv[]) {
       velocityDirichletFunctions[i] = new VelocityDirichletCondition();
 
       noNodes[i] = Dune::BitSetVector<dims>(leafVertexCount);
-      dirichletNodes[i] = Dune::BitSetVector<dims>(leafVertexCount);
+
       auto& gridDirichletNodes = dirichletNodes[i];
+      gridDirichletNodes = Dune::BitSetVector<dims>(leafVertexCount);
+
       for (int j=0; j<leafVertexCount; j++) {
         if (leafFaces[i]->right.containsVertex(j))
           gridDirichletNodes[j][0] = true;
@@ -322,7 +300,42 @@ int main(int argc, char *argv[]) {
       }
     }
 
+    // extend Dirichlet information to all grid levels
+    /*            BoundaryPatchProlongator<GridType>::prolong(dirichletBoundary[i]);
+                dirichletNodes[i].resize(toplevel+1);
+
+                for (int j=0; j<=toplevel; j++) {
+
+                    int fSSize = grids[i]->size(j,dim);
+                    dirichletNodes[i][j].resize(fSSize);
+                    for (int k=0; k<fSSize; k++)
+                        for (int l=0; l<dim; l++)
+                            dirichletNodes[i][j][k][l] = dirichletBoundary[i][j].containsVertex(k);
+
+                }
+
+                sampleOnBitField(*grids[i], dirichletValues[i], dirichletNodes[i]);*/
+
+    const int toplevel = deformedGrids[0]->maxLevel();
+    std::vector<Dune::BitSetVector<dims> > totalDirichletNodes(toplevel+1);
+
+    for (int i=0; i<=toplevel; i++) {
+      int totalSize = 0;
+      for (int j=0; j<deformedGrids.size(); j++)
+        totalSize += deformedGrids[j]->size(i, dims);
+
+        totalDirichletNodes[i].resize(totalSize);
+
+        int idx=0;
+        for (int j=0; j<deformedGrids.size(); j++)
+          for (size_t k=0; k<dirichletNodes[j][i].size(); k++)
+            for (int l=0; l<dims; l++)
+              totalDirichletNodes[i][idx++][l] = dirichletNodes[j][i][k][l];
+     }
+
+    // ------------------------------------------------------------------------------------
     // set up individual assemblers for each body, assemble problem (matrices, forces, rhs)
+    // ------------------------------------------------------------------------------------
     std::vector<std::shared_ptr<Assembler>> assemblers(bodyCount);
     Matrices<Matrix> matrices(bodyCount);
 
@@ -355,40 +368,9 @@ int main(int argc, char *argv[]) {
       };
     }
 
-    /* Jonny 2bcontact
-    // make dirichlet bitfields containing dirichlet information for both grids
-           int size = rhs[0].size() + rhs[1].size();
-
-           Dune::BitSetVector<dims> totalDirichletNodes(size);
-
-           for (size_t i=0; i<dirichletNodes[0].size(); i++)
-               for (int j=0; j<dims; j++)
-                   totalDirichletNodes[i][j] = dirichletNodes[0][i][j];
-
-           int offset = rhs[0].size();
-           for (size_t i=0; i<dirichletNodes[1].size(); i++)
-               for (int j=0; j<dims; j++)
-                   totalDirichletNodes[offset + i][j] = dirichletNodes[1][i][j];
-
-           // assemble separate linear elasticity problems
-           std::array<MatrixType,2> stiffnessMatrix;
-           std::array<const MatrixType*, 2> submat;
-
-           for (size_t i=0; i<2; i++) {
-               OperatorAssembler<P1Basis,P1Basis> globalAssembler(*p1Basis[i],*p1Basis[i]);
-
-               double s = (i==0) ? E0 : E1;
-               StVenantKirchhoffAssembler<GridType, P1Basis::LocalFiniteElement, P1Basis::LocalFiniteElement> localAssembler(s, nu);
-
-               globalAssembler.assemble(localAssembler, stiffnessMatrix[i]);
-
-               submat[i] = &stiffnessMatrix[i];
-           }
-
-           MatrixType bilinearForm;
-           contactAssembler.assembleJacobian(submat, bilinearForm);
- */
-
+    // -----------------
+    // init input/output
+    // -----------------
     using MyProgramState = ProgramState<Vector, ScalarVector>;
     MyProgramState programState(leafVertexCounts);
     auto const firstRestart = parset.get<size_t>("io.restarts.first");
@@ -418,11 +400,10 @@ int main(int argc, char *argv[]) {
       restartIO->read(firstRestart, programState);
     else
      programState.setupInitialConditions(parset, nBodyAssembler, externalForces,
-                                          matrices, assemblers, dirichletNodes,
+                                          matrices, assemblers, totalDirichletNodes,
                                           noNodes, frictionalBoundaries, body);
 
 
-
     // assemble friction
     std::vector<std::shared_ptr<MyGlobalFrictionData<LocalVector>>> frictionInfo(weakPatches.size());
     std::vector<std::shared_ptr<GlobalFriction<Matrix, Vector>>> globalFriction(weakPatches.size());
@@ -432,7 +413,6 @@ int main(int argc, char *argv[]) {
         globalFriction[i]->updateAlpha(programState.alpha[i]);
     }
 
-
     using MyVertexBasis = typename Assembler::VertexBasis;
     using MyCellBasis = typename Assembler::CellBasis;
     std::vector<Vector> vertexCoordinates(leafVertexCounts.size());
@@ -499,9 +479,11 @@ int main(int argc, char *argv[]) {
     };
     report(true);
 
+    // -------------------
     // Set up TNNMG solver
+    // -------------------
     using NonlinearFactory = SolverFactory<DeformedGrid, GlobalFriction<Matrix, Vector>, Matrix, Vector>;
-    NonlinearFactory factory(parset.sub("solver.tnnmg"), nBodyAssembler, dirichletNodes);
+    NonlinearFactory factory(parset.sub("solver.tnnmg"), nBodyAssembler, totalDirichletNodes);
 
     using MyUpdater = Updaters<RateUpdater<Vector, Matrix, Function, dims>,
                                StateUpdater<ScalarVector, Vector>>;
@@ -518,8 +500,7 @@ int main(int argc, char *argv[]) {
             L, V0));
 
 
-    auto const refinementTolerance =
-        parset.get<double>("timeSteps.refinementTolerance");
+    auto const refinementTolerance = parset.get<double>("timeSteps.refinementTolerance");
     auto const mustRefine = [&](MyUpdater &coarseUpdater,
                                 MyUpdater &fineUpdater) {
       std::vector<ScalarVector> coarseAlpha;
@@ -540,11 +521,11 @@ int main(int argc, char *argv[]) {
     std::signal(SIGINT, handleSignal);
     std::signal(SIGTERM, handleSignal);
 
-    AdaptiveTimeStepper<NonlinearFactory, MyUpdater,
-                        EnergyNorm<ScalarMatrix, ScalarVector>>
+    AdaptiveTimeStepper<NonlinearFactory, MyUpdater, EnergyNorm<ScalarMatrix, ScalarVector>>
         adaptiveTimeStepper(factory, parset, globalFriction, current,
                             programState.relativeTime, programState.relativeTau,
                             externalForces, stateEnergyNorms, mustRefine);
+
     while (!adaptiveTimeStepper.reachedEnd()) {
       programState.timeStep++;
 

src/spatial-solving/solverfactory.cc

@@ -16,7 +16,7 @@
 template <class DeformedGrid, class Nonlinearity, class Matrix, class Vector>
 SolverFactory<DeformedGrid, Nonlinearity, Matrix, Vector>::SolverFactory(
     const Dune::ParameterTree& parset, const Dune::Contact::NBodyAssembler<DeformedGrid, Vector>& nBodyAssembler,
-    const std::vector<Dune::BitSetVector<DeformedGrid::dimension>>& ignoreNodes)
+    const std::vector<std::vector<Dune::BitSetVector<DeformedGrid::dimension>>>& ignoreNodes)
     : nBodyAssembler_(nBodyAssembler),
       baseEnergyNorm_(baseSolverStep_),
       baseSolver_(&baseSolverStep_,
@@ -28,11 +28,11 @@ SolverFactory<DeformedGrid, Nonlinearity, Matrix, Vector>::SolverFactory(
   multigridStep_ = std::make_shared<Step>();
   // tnnmg iteration step
   multigridStep_->setMGType(parset.get<int>("main.multi"), parset.get<int>("main.pre"), parset.get<int>("main.post"));
-  //multigridStep_->setIgnore(ignoreNodes);
+  multigridStep_->setIgnore(ignoreNodes);
   multigridStep_->setBaseSolver(baseSolver_);
   multigridStep_->setSmoother(&presmoother_, &postsmoother_);
- // multigridStep_->setHasObstacle(nBodyAssembler_.totalHasObstacle_);
- // multigridStep_->setObstacles(nBodyAssembler_.totalObstacles_);
+  multigridStep_->setHasObstacle(nBodyAssembler_.totalHasObstacle_);
+  multigridStep_->setObstacles(nBodyAssembler_.totalObstacles_);
 
   // create the transfer operators
   const int nCouplings = nBodyAssembler_.nCouplings();