.

src/explicitgrid.hh

@@ -2,7 +2,13 @@
 #define SRC_EXPLICITGRID_HH
 
 #include "gridselector.hh"
+#include "explicitvectors.hh"
+
+#include <dune/contact/common/deformedcontinuacomplex.hh>
 
 using LeafGridView = Grid::LeafGridView;
 using LevelGridView = Grid::LevelGridView;
+
+using DeformedGridComplex = typename Dune::Contact::DeformedContinuaComplex<Grid, Vector>;
+using DeformedGrid = DeformedGridComplex::DeformedGridType;
 #endif

src/multi-body-problem.cc

@@ -65,11 +65,11 @@
 #include "multi-body-problem-data/myglobalfrictiondata.hh"
 #include "multi-body-problem-data/mygrids.hh"
 #include "multi-body-problem-data/weakpatch.hh"
-#include "spatial-solving/solverfactory.hh"
-#include "time-stepping/adaptivetimestepper.hh"
+//#include "spatial-solving/solverfactory.hh"
+//#include "time-stepping/adaptivetimestepper.hh"
 #include "time-stepping/rate.hh"
 #include "time-stepping/state.hh"
-#include "time-stepping/updaters.hh"
+//#include "time-stepping/updaters.hh"
 #include "vtk.hh"
 
 // for getcwd
@@ -194,8 +194,21 @@ int main(int argc, char *argv[]) {
     for (size_t i=0; i<grids.size(); i++) {
         deformedGridComplex.addGrid(*grids[i], "body" + std::to_string(i));
     }
+
     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
+    */
+
     std::vector<std::shared_ptr<DefLeafGridView>> leafViews(deformedGrids.size());
     std::vector<std::shared_ptr<DefLevelGridView>> levelViews(deformedGrids.size());
     std::vector<int> leafVertexCounts(deformedGrids.size());

src/spatial-solving/fixedpointiterator.cc

@@ -32,8 +32,8 @@ void FixedPointIterationCounter::operator+=(
   multigridIterations += other.multigridIterations;
 }
 
-template <class DeformedGrid, class Factory, class Updaters, class ErrorNorm>
-FixedPointIterator<DeformedGrid, Factory, Updaters, ErrorNorm>::FixedPointIterator(
+template <class Factory, class Updaters, class ErrorNorm>
+FixedPointIterator<Factory, Updaters, ErrorNorm>::FixedPointIterator(
     const Dune::Contact::NBodyAssembler<DeformedGrid, Vector>& nBodyAssembler,
     Factory &factory, Dune::ParameterTree const &parset,
     std::shared_ptr<Nonlinearity> globalFriction, ErrorNorm const &errorNorm)
@@ -49,9 +49,9 @@ FixedPointIterator<DeformedGrid, Factory, Updaters, ErrorNorm>::FixedPointIterat
       verbosity_(parset.get<Solver::VerbosityMode>("v.solver.verbosity")),
       errorNorm_(errorNorm) {}
 
-template <class DeformedGrid,class Factory, class Updaters, class ErrorNorm>
+template <class Factory, class Updaters, class ErrorNorm>
 FixedPointIterationCounter
-FixedPointIterator<DeformedGrid, Factory, Updaters, ErrorNorm>::run(
+FixedPointIterator<Factory, Updaters, ErrorNorm>::run(
     Updaters updaters, Matrix const &velocityMatrix, Vector const &velocityRHS,
     Vector &velocityIterate) {
   EnergyNorm<Matrix, Vector> energyNorm(velocityMatrix);
@@ -66,7 +66,7 @@ FixedPointIterator<DeformedGrid, Factory, Updaters, ErrorNorm>::run(
   for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations_;
        ++fixedPointIteration) {
     // solve a velocity problem
-    globalFriction_->updateAlpha(alpha);
+   // globalFriction_->updateAlpha(alpha);
     ConvexProblem convexProblem(1.0, velocityMatrix, *globalFriction_,
                                 velocityRHS, velocityIterate);
     BlockProblem velocityProblem(parset_, convexProblem);
@@ -82,7 +82,7 @@ FixedPointIterator<DeformedGrid, Factory, Updaters, ErrorNorm>::run(
 
     for (size_t i=0; i<v_m.size(); i++) {
       v_m[i] *= 1.0 - lambda_;
-      Arithmetic::addProduct(v_m[i], lambda_, velocityIterate[i]);
+      //Arithmetic::addProduct(v_m[i], lambda_, velocityIterate[i]);
     }
 
     // compute relative velocities on contact boundaries
@@ -91,14 +91,14 @@ FixedPointIterator<DeformedGrid, Factory, Updaters, ErrorNorm>::run(
     // solve a state problem
     updaters.state_->solve(v_m);
     ScalarVector newAlpha;
-    updaters.state_->extractAlpha(newAlpha);
+   /* updaters.state_->extractAlpha(newAlpha);
 
     if (lambda_ < 1e-12 or
         errorNorm_.diff(alpha, newAlpha) < fixedPointTolerance_) {
       fixedPointIteration++;
       break;
     }
-    alpha = newAlpha;
+    alpha = newAlpha;*/
   }
   if (fixedPointIteration == fixedPointMaxIterations_)
     DUNE_THROW(Dune::Exception, "FPI failed to converge");
@@ -119,8 +119,10 @@ std::ostream &operator<<(std::ostream &stream,
                 << ")";
 }
 
-template <class DeformedGrid, class Factory, class Updaters, class ErrorNorm>
-void FixedPointIterator<DeformedGrid, Factory, Updaters, ErrorNorm>::relativeVelocities(std::vector<Vector>& v_m) const {
+template <class Factory, class Updaters, class ErrorNorm>
+void FixedPointIterator<Factory, Updaters, ErrorNorm>::relativeVelocities(std::vector<Vector>& v_m) const {
+  using field_type = typename Factory::Matrix::field_type;
+
   // adaptation of DualMortarCoupling::setup()
 
   const size_t dim = DeformedGrid::dimension;
@@ -134,7 +136,7 @@ void FixedPointIterator<DeformedGrid, Factory, Updaters, ErrorNorm>::relativeVel
   using DualCache = Dune::Contact::DualBasisAdapter<GridView, field_type>;
   std::unique_ptr<DualCache> dualCache;
   dualCache = std::make_unique< Dune::Contact::DualBasisAdapterGlobal<GridView, field_type> >();
-
+/*
   // define FE grid functions
   std::vector<BasisGridFunction<> > gridFunctions(nBodyAssembler_.nGrids());
   for (size_t i=0; i<gridFunctions.size(); i++) {
@@ -245,9 +247,6 @@ void FixedPointIterator<DeformedGrid, Factory, Updaters, ErrorNorm>::relativeVel
 
   }
 
-
-
-
       ///////////////////////////////////
       //  reducing nonmortar boundary
       /////////////////////////////////
@@ -278,7 +277,7 @@ void FixedPointIterator<DeformedGrid, Factory, Updaters, ErrorNorm>::relativeVel
       //writeBoundary(boundaryWithMapping,debugPath_ + "relevantNonmortar");
 
 
-      /** \todo replace by all fine grid segments which are totally covered by the RemoteIntersections. */
+      // \todo replace by all fine grid segments which are totally covered by the RemoteIntersections.
       //for (const auto& rIs : intersections(glue))
       //    boundaryWithMapping.addFace(rIs.inside(),rIs.indexInInside());
 
@@ -303,7 +302,7 @@ void FixedPointIterator<DeformedGrid, Factory, Updaters, ErrorNorm>::relativeVel
       //   Get the occupation structure for the mortar matrix
       // ///////////////////////////////////////////////////////////
 
-      /** \todo Also restrict mortar indices and don't use the whole grid level. */
+      // todo Also restrict mortar indices and don't use the whole grid level.
       Dune::MatrixIndexSet mortarIndices(nonmortarBoundary_.numVertices(), grid1_->size(dim));
 
       // Create mapping from the global set of block dofs to the ones on the contact boundary
@@ -382,7 +381,7 @@ void FixedPointIterator<DeformedGrid, Factory, Updaters, ErrorNorm>::relativeVel
       for(size_t rowIdx=0; rowIdx<weakObstacle_.size(); rowIdx++)
           weakObstacle_[rowIdx] *= D[rowIdx][rowIdx][0][0];
 
-      gridGlueBackend_->clear();
+      gridGlueBackend_->clear(); */
 }
 
 #include "fixedpointiterator_tmpl.cc"

src/spatial-solving/fixedpointiterator.hh

@@ -8,6 +8,8 @@
 #include <dune/solvers/norms/norm.hh>
 #include <dune/solvers/solvers/solver.hh>
 
+#include <dune/contact/assemblers/nbodyassembler.hh>
+
 struct FixedPointIterationCounter {
   size_t iterations = 0;
   size_t multigridIterations = 0;
@@ -27,10 +29,16 @@ class FixedPointIterator {
   using BlockProblem = typename Factory::BlockProblem;
   using Nonlinearity = typename ConvexProblem::NonlinearityType;
 
+  using DeformedGrid = typename Factory::DeformedGrid;
+
+private:
+  void relativeVelocities(std::vector<Vector>& v_m) const;
+
 public:
-  FixedPointIterator(Factory &factory, Dune::ParameterTree const &parset,
+  FixedPointIterator(const Dune::Contact::NBodyAssembler<DeformedGrid, Vector>& nBodyAssembler,
+                     Factory &factory, const Dune::ParameterTree& parset,
                      std::shared_ptr<Nonlinearity> globalFriction,
-                     ErrorNorm const &errorNorm_);
+                     const ErrorNorm& errorNorm);
 
   FixedPointIterationCounter run(Updaters updaters,
                                  Matrix const &velocityMatrix,
@@ -38,6 +46,7 @@ class FixedPointIterator {
                                  Vector &velocityIterate);
 
 private:
+  const Dune::Contact::NBodyAssembler<DeformedGrid, Vector>& nBodyAssembler_;
   std::shared_ptr<typename Factory::Step> step_;
   Dune::ParameterTree const &parset_;
   std::shared_ptr<Nonlinearity> globalFriction_;

src/spatial-solving/fixedpointiterator_tmpl.cc

@@ -10,6 +10,8 @@
 #include <dune/solvers/norms/energynorm.hh>
 #include <dune/tnnmg/problem-classes/convexproblem.hh>
 
+#include <dune/contact/common/deformedcontinuacomplex.hh>
+
 #include <dune/tectonic/globalfriction.hh>
 #include <dune/tectonic/myblockproblem.hh>
 
@@ -22,7 +24,7 @@ using Function = Dune::VirtualFunction<double, double>;
 using Factory = SolverFactory<
     MY_DIM,
     MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
-    Grid>;
+    DeformedGrid>;
 using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
 using MyRateUpdater = RateUpdater<Vector, Matrix, Function, MY_DIM>;
 using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;

src/spatial-solving/solverfactory.cc

@@ -11,9 +11,9 @@
 
 #include "solverfactory.hh"
 
-template <size_t dim, class BlockProblem, class Grid>
-SolverFactory<dim, BlockProblem, Grid>::SolverFactory(
-    Dune::ParameterTree const &parset, Grid const &grid,
+template <size_t dim, class BlockProblem, class DeformedGrid>
+SolverFactory<dim, BlockProblem, DeformedGrid>::SolverFactory(
+    Dune::ParameterTree const &parset, const DeformedGrid& grid,
     Dune::BitSetVector<dim> const &ignoreNodes)
     : baseEnergyNorm(linearBaseSolverStep),
       linearBaseSolver(&linearBaseSolverStep,
@@ -33,7 +33,7 @@ SolverFactory<dim, BlockProblem, Grid>::SolverFactory(
   // transfer operators
   for (auto &&x : transferOperators)
     x = new CompressedMultigridTransfer<Vector>;
-  TransferOperatorAssembler<Grid>(grid)
+  TransferOperatorAssembler<DeformedGrid>(grid)
       .assembleOperatorPointerHierarchy(transferOperators);
   linearIterationStep.setTransferOperators(transferOperators);
 
@@ -44,14 +44,14 @@ SolverFactory<dim, BlockProblem, Grid>::SolverFactory(
   multigridStep->ignoreNodes_ = &ignoreNodes;
 }
 
-template <size_t dim, class BlockProblem, class Grid>
-SolverFactory<dim, BlockProblem, Grid>::~SolverFactory() {
+template <size_t dim, class BlockProblem, class DeformedGrid>
+SolverFactory<dim, BlockProblem, DeformedGrid>::~SolverFactory() {
   for (auto &&x : transferOperators)
     delete x;
 }
 
-template <size_t dim, class BlockProblem, class Grid>
-auto SolverFactory<dim, BlockProblem, Grid>::getStep()
+template <size_t dim, class BlockProblem, class DeformedGrid>
+auto SolverFactory<dim, BlockProblem, DeformedGrid>::getStep()
     -> std::shared_ptr<Step> {
   return multigridStep;
 }

src/spatial-solving/solverfactory.hh

@@ -12,7 +12,7 @@
 #include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
 #include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
 
-template <size_t dim, class BlockProblemTEMPLATE, class Grid>
+template <size_t dim, class BlockProblemTEMPLATE, class DeformedGridTEMPLATE>
 class SolverFactory {
 public:
   using BlockProblem = BlockProblemTEMPLATE;
@@ -20,6 +20,8 @@ class SolverFactory {
   using Vector = typename BlockProblem::VectorType;
   using Matrix = typename BlockProblem::MatrixType;
 
+  using DeformedGrid = DeformedGridTEMPLATE;
+
 private:
   using NonlinearSmoother = GenericNonlinearGS<BlockProblem>;
 
@@ -27,7 +29,7 @@ class SolverFactory {
   using Step =
       TruncatedNonsmoothNewtonMultigrid<BlockProblem, NonlinearSmoother>;
 
-  SolverFactory(Dune::ParameterTree const &parset, Grid const &grid,
+  SolverFactory(Dune::ParameterTree const &parset, const DeformedGrid& grid,
                 Dune::BitSetVector<dim> const &ignoreNodes);
 
   ~SolverFactory();