Construct our own MMGStep object

Because we need to build specialized multigrid transfer operators
for the periodic bases.

src/quasiconvexity-test.cc

@@ -31,6 +31,7 @@
 #include <dune/fufem/functiontools/boundarydofs.hh>
 #include <dune/fufem/dunepython.hh>
 
+#include <dune/solvers/iterationsteps/mmgstep.hh>
 #include <dune/solvers/solvers/iterativesolver.hh>
 #include <dune/solvers/norms/energynorm.hh>
 
@@ -47,6 +48,38 @@ const int order = 1;
 
 using namespace Dune;
 
+template <class GridView>
+Functions::BasisFactory::PeriodicIndexSet getPeriodicIndices(const GridView& gridView)
+{
+  // Check whether two points are equal on R/Z x R/Z
+  auto equivalent = [](const FieldVector<double,2>& x, const FieldVector<double,2>& y)
+  {
+    return ( (FloatCmp::eq(x[0],y[0]) or FloatCmp::eq(x[0]+1,y[0]) or FloatCmp::eq(x[0]-1,y[0]))
+           and (FloatCmp::eq(x[1],y[1]) or FloatCmp::eq(x[1]+1,y[1]) or FloatCmp::eq(x[1]-1,y[1])) );
+  };
+
+  // Extract all boundary vertices
+  std::vector<std::pair<std::size_t, FieldVector<double,dim> > > boundaryVertices;
+  BoundaryPatch<GridView> domainBoundary(gridView, true);
+
+  for (const auto& v1 : vertices(gridView))
+    if (domainBoundary.containsVertex(gridView.indexSet().index(v1)))
+      boundaryVertices.push_back(std::make_pair(gridView.indexSet().index(v1), v1.geometry().corner(0)));
+
+  std::cout << "Grid has " << boundaryVertices.size() << " boundary vertices" << std::endl;
+
+  Functions::BasisFactory::PeriodicIndexSet periodicIndices;
+  if (order!=1)
+    DUNE_THROW(NotImplemented, "Periodicity detection only implemented for order==1");
+
+  for (const auto& v1 : boundaryVertices)
+    for (const auto& v2 : boundaryVertices)
+      if (equivalent(v1.second, v2.second))
+        periodicIndices.unifyIndexPair(v1.first, v2.first);
+
+  return periodicIndices;
+}
+
 template <class Basis, class Vector>
 void writeDisplacement(const Basis& basis, const Vector& periodicX, const FieldMatrix<double,2,2>& F0, std::string filename)
 {
@@ -230,6 +263,8 @@ int main (int argc, char *argv[]) try
   GridView gridView = grid->leafGridView();
 #endif
 
+  bool periodic = parameterSet.get<bool>("periodic", false);
+
   // Check whether two points are equal on R/Z x R/Z
   auto equivalent = [](const FieldVector<double,2>& x, const FieldVector<double,2>& y)
   {
@@ -243,31 +278,10 @@ int main (int argc, char *argv[]) try
     return FloatCmp::eq(x[0],y[0]) && FloatCmp::eq(x[1],y[1]);
   };
 
-  // Extract all boundary vertices
-  std::vector<std::pair<std::size_t, FieldVector<double,dim> > > boundaryVertices;
-  BoundaryPatch<GridView> domainBoundary(gridView, true);
-
-  for (const auto& v1 : vertices(gridView))
-    if (domainBoundary.containsVertex(gridView.indexSet().index(v1)))
-      boundaryVertices.push_back(std::make_pair(gridView.indexSet().index(v1), v1.geometry().corner(0)));
-
-  std::cout << "Grid has " << boundaryVertices.size() << " boundary vertices" << std::endl;
-
-  bool periodic = parameterSet.get<bool>("periodic", false);
-
-  //using namespace Dune::Fufem::BasisFactory;
   using namespace Dune::Functions::BasisFactory;
   Functions::BasisFactory::PeriodicIndexSet periodicIndices;
   if (periodic)
-  {
-    if (order!=1)
-      DUNE_THROW(NotImplemented, "Periodicity detection only implemented for order==1");
-
-    for (const auto& v1 : boundaryVertices)
-      for (const auto& v2 : boundaryVertices)
-        if (equivalent(v1.second, v2.second))
-          periodicIndices.unifyIndexPair(v1.first, v2.first);
-  }
+    periodicIndices = getPeriodicIndices(gridView);
 
   // FE basis spanning the FE space that we are working in
   auto feBasis = makeBasis(
@@ -339,7 +353,6 @@ int main (int argc, char *argv[]) try
 
   SolutionType x(feBasis.size());
 
-  using namespace Dune::Functions::BasisFactory;
   auto powerBasis = makeBasis(
     gridView,
     power<dim>(
@@ -508,27 +521,114 @@ int main (int argc, char *argv[]) try
 
     auto localADOLCStiffness = std::make_shared<Elasticity::LocalADOLCStiffness<FEBasis::LocalView> >(elasticEnergy);
 
-    Elasticity::FEAssembler<FEBasis,SolutionType> assembler(feBasis, localADOLCStiffness);
+    auto assembler = std::make_shared<Elasticity::FEAssembler<FEBasis,SolutionType> >(feBasis, localADOLCStiffness);
 
     ///////////////////////////////////////////////////
-    //   Create a Riemannian trust-region solver
+    //   Create a trust-region solver
     ///////////////////////////////////////////////////
 
+    using Matrix = BCRSMatrix<FieldMatrix<double, dim, dim> >;
+
+#ifdef HAVE_IPOPT
+    // First create an IPOpt base solver
+    auto baseSolver = std::make_shared<QuadraticIPOptSolver<Matrix,SolutionType> >(
+        baseTolerance,
+        baseIterations,
+        NumProc::REDUCED,
+        "mumps");
+#else
+    // First create a Gauss-seidel base solver
+    TrustRegionGSStep<Matrix, CorrectionType>* baseSolverStep = new TrustRegionGSStep<Matrix, SolutionType>;
+
+    // Hack: the two-norm may not scale all that well, but it is fast!
+    TwoNorm<CorrectionType>* baseNorm = new TwoNorm<CorrectionType>;
+
+    auto baseSolver = std::make_shared<::LoopSolver<CorrectionType>>(baseSolverStep,
+                                                                     baseIterations,
+                                                                     baseTolerance,
+                                                                     baseNorm,
+                                                                     Solver::QUIET);
+#endif
+    // Make pre and postsmoothers
+    auto presmoother  = std::make_shared< TrustRegionGSStep<Matrix, SolutionType> >();
+    auto postsmoother = std::make_shared< TrustRegionGSStep<Matrix, SolutionType> >();
+
+    auto mmgStep = std::make_shared<MonotoneMGStep<Matrix, SolutionType> >();
+
+    mmgStep->setVerbosity(NumProc::FULL);
+    mmgStep->setMGType(mu, nu1, nu2);
+    mmgStep->ignoreNodes_ = &dirichletDofs;
+    mmgStep->setBaseSolver(baseSolver);
+    mmgStep->setSmoother(presmoother, postsmoother);
+    mmgStep->setObstacleRestrictor(std::make_shared<MandelObstacleRestrictor<SolutionType> >());
+
+    //////////////////////////////////////
+    //   Create the transfer operators
+    //////////////////////////////////////
+
+    using TransferOperatorType = typename TruncatedCompressedMGTransfer<SolutionType>::TransferOperatorType;
+    std::vector<std::shared_ptr<TruncatedCompressedMGTransfer<SolutionType>>> transferOperators(numLevels-1);
+
+    // For the restriction operators: FE bases on all levels
+
+    // This particular object is only needed to get the type of the level basis
+    PeriodicIndexSet dummyPeriodicIndices;
+    auto dummyLevelBasis = makeBasis(
+      grid->levelGridView(0),
+      power<dim>(
+        Functions::BasisFactory::periodic(lagrange<order>(), dummyPeriodicIndices),
+        blockedInterleaved()
+    ));
+
+    using LevelBasis = decltype(dummyLevelBasis);
+
+    // Construct the actual level bases
+    std::vector<std::shared_ptr<LevelBasis> > levelBases(numLevels);
+    for (int j=0; j<numLevels; j++)
+    {
+      PeriodicIndexSet periodicIndices;
+      if (periodic)
+        periodicIndices = getPeriodicIndices(grid->levelGridView(j));
+      levelBases[j] = std::make_shared<LevelBasis>(makeBasis(
+        grid->levelGridView(j),
+        power<dim>(
+          Functions::BasisFactory::periodic(lagrange<order>(), periodicIndices),
+          blockedInterleaved()
+      )));
+    }
+
+    for (int j=0; j<numLevels-1; j++)
+    {
+      auto transferMatrix = std::make_shared<TransferOperatorType>();
+      assembleGlobalBasisTransferMatrix(*transferMatrix, *levelBases[j], *levelBases[j+1]);
+
+      // Construct the local multigrid transfer matrix
+      transferOperators[j] = std::make_shared<TruncatedCompressedMGTransfer<SolutionType> >();
+      transferOperators[j]->setMatrix(transferMatrix);
+    }
+
+    {
+      auto transferMatrix = std::make_shared<TransferOperatorType>();
+      assembleGlobalBasisTransferMatrix(*transferMatrix, *levelBases.back(), *levelBases.back());
+
+      // Construct the local multigrid transfer matrix
+      transferOperators.push_back(std::make_shared<TruncatedCompressedMGTransfer<SolutionType> >());
+      transferOperators.back()->setMatrix(transferMatrix);
+    }
+
+    mmgStep->setTransferOperators(transferOperators);
+
     TrustRegionSolver<FEBasis,SolutionType> solver;
     solver.setup(*grid,
-                 &assembler,
+                 assembler,
                  x,
                  dirichletDofs,
                  tolerance,
                  maxTrustRegionSteps,
                  initialTrustRegionRadius,
-                 multigridIterations,
+                 mmgStep,
                  mgTolerance,
-                 mu, nu1, nu2,
-                 baseIterations,
-                 baseTolerance,
-                 1.0
-                );
+                 multigridIterations);
 
     ////////////////////////////////////////////////////////
     //   Set Dirichlet values