Minimize W_magic starting from the compatible deformation

src/quasiconvexity-test-micromorphic.cc

@@ -37,6 +37,7 @@
 #include <dune/elasticity/common/trustregionsolver.hh>
 #include <dune/elasticity/assemblers/localadolcstiffness.hh>
 #include <dune/elasticity/assemblers/feassembler.hh>
+#include <dune/elasticity/materials/localintegralenergy.hh>
 #include <dune/elasticity/materials/micromorphicallyrelaxedenergy.hh>
 #include <dune/elasticity/materials/pdistanceenergy.hh>
 #include <dune/elasticity/materials/quasiconvexitytestdensities.hh>
@@ -190,7 +191,7 @@ void writeDisplacement(const Basis& basis, const Vector& x, const FieldMatrix<do
 
 
 template <typename Grid, typename IncompatibleDeformationGradientFunction>
-void constructDeformation(const std::shared_ptr<Grid>& grid,
+auto constructDeformation(const std::shared_ptr<Grid>& grid,
                           const ParameterTree& parameterSet,
                           const IncompatibleDeformationGradientFunction& incompatibleDeformationGradient)
 {
@@ -403,6 +404,207 @@ void constructDeformation(const std::shared_ptr<Grid>& grid,
   std::cout << "2-Distance to incompatible strain: " << std::sqrt(assembler->computeEnergy(x)) << std::endl;
 
   writeDisplacement(feBasis, x, F0, resultPath + "quasiconvexity-test-micromorphic-deformation");
+
+  return x;
+}
+
+template <typename Grid, typename Vector>
+void minimizeEnergyInSpaceOfDeformations(const std::shared_ptr<Grid>& grid,
+                                         const ParameterTree& parameterSet,
+                                         const Vector& initialIterate)
+{
+  // read solver settings
+  const int numLevels                   = parameterSet.get<int>("numLevels");
+  const double tolerance                = parameterSet.get<double>("tolerance");
+  const int maxTrustRegionSteps         = parameterSet.get<int>("maxTrustRegionSteps");
+  const double initialTrustRegionRadius = parameterSet.get<double>("initialTrustRegionRadius");
+  const int multigridIterations         = parameterSet.get<int>("numIt");
+  const int nu1                         = parameterSet.get<int>("nu1");
+  const int nu2                         = parameterSet.get<int>("nu2");
+  const int mu                          = parameterSet.get<int>("mu");
+  const int baseIterations              = parameterSet.get<int>("baseIt");
+  const double mgTolerance              = parameterSet.get<double>("mgTolerance");
+  const double baseTolerance            = parameterSet.get<double>("baseTolerance");
+  std::string resultPath                = parameterSet.get("resultPath", "");
+
+#if HAVE_DUNE_PARMG
+  using GridView = typename Grid::LevelGridView;
+  GridView gridView = grid->levelGridView(0);
+#else
+  using GridView = typename Grid::LeafGridView;
+  GridView gridView = grid->leafGridView();
+#endif
+
+  // FE basis spanning the FE space that we are working in
+  using namespace Dune::Functions::BasisFactory;
+  auto feBasis = makeBasis(
+      gridView,
+      power<dim>(
+        lagrange<order>(),
+        blockedInterleaved()
+    ));
+  using FEBasis = decltype(feBasis);
+
+  std::cout << "Basis has " << feBasis.size() << " degrees of freedom" << std::endl;
+
+  ///////////////////////////////////////////
+  //   Set Dirichlet values
+  ///////////////////////////////////////////
+
+  BoundaryPatch<GridView> dirichletBoundary(gridView, true);
+
+  BitSetVector<dim> dirichletDofs(feBasis.size(), false);
+  constructBoundaryDofs(dirichletBoundary,feBasis,dirichletDofs);
+
+  ////////////////////////////
+  //   Initial iterate
+  ////////////////////////////
+
+  Vector x = initialIterate;
+
+  // Set the underlying affine deformation
+  auto x0 = parameterSet.get<double>("x0");
+
+  FieldMatrix<double,2,2> F0 = { {sqrt(x0), 0},
+                                 {       0, 1.0 / sqrt(x0)} };
+
+  // Output initial iterate
+  writeDisplacement(feBasis, x, F0, resultPath + "quasiconvexity-test-micromorphic-deformation-magic-initial");
+
+  //////////////////////////////////////////////////////////////
+  //   Create an assembler for the energy functional
+  //////////////////////////////////////////////////////////////
+
+  const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
+
+  // Assembler using ADOL-C
+  std::cout << "Selected energy is: " << parameterSet.get<std::string>("energy") << std::endl;
+  std::shared_ptr<Elasticity::LocalDensity<dim,adouble> > density;
+
+  if (parameterSet.get<std::string>("energy") == "magic")
+    density = std::make_shared<Elasticity::MagicDensity<dim,adouble> >(F0,materialParameters);
+
+  if(!density)
+    DUNE_THROW(Exception, "Error: Selected energy not available!");
+
+  auto elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<typename FEBasis::LocalView,
+                                                                        adouble> >(density);
+
+  auto localADOLCStiffness = std::make_shared<Elasticity::LocalADOLCStiffness<typename FEBasis::LocalView> >(elasticEnergy);
+
+  auto assembler = std::make_shared<Elasticity::FEAssembler<FEBasis,Vector> >(feBasis, localADOLCStiffness);
+
+  ///////////////////////////////////////////////////
+  //   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,Vector> >(
+      baseTolerance,
+      baseIterations,
+      NumProc::QUIET,
+      "mumps");
+#else
+  // First create a Gauss-seidel base solver
+  TrustRegionGSStep<Matrix, Vector>* baseSolverStep = new TrustRegionGSStep<Matrix, Vector>;
+
+  // Hack: the two-norm may not scale all that well, but it is fast!
+  TwoNorm<CorrectionType>* baseNorm = new TwoNorm<Vector>;
+
+  auto baseSolver = std::make_shared<::LoopSolver<Vector>>(baseSolverStep,
+                                                           baseIterations,
+                                                           baseTolerance,
+                                                           baseNorm,
+                                                           Solver::QUIET);
+#endif
+  // Make pre and postsmoothers
+  auto presmoother  = std::make_shared< TrustRegionGSStep<Matrix, Vector> >();
+  auto postsmoother = std::make_shared< TrustRegionGSStep<Matrix, Vector> >();
+
+  auto mmgStep = std::make_shared<MonotoneMGStep<Matrix, Vector> >();
+
+  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<Vector> >());
+
+  //////////////////////////////////////
+  //   Create the transfer operators
+  //////////////////////////////////////
+
+  using TransferOperatorType = typename TruncatedCompressedMGTransfer<Vector>::TransferOperatorType;
+  std::vector<std::shared_ptr<TruncatedCompressedMGTransfer<Vector>>> transferOperators(numLevels-1);
+
+  // For the restriction operators: FE bases on all levels
+  auto dummyLevelBasis = makeBasis(
+    grid->levelGridView(0),
+    power<dim>(
+      lagrange<order>(),
+      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++)
+  {
+    levelBases[j] = std::make_shared<LevelBasis>(makeBasis(
+      grid->levelGridView(j),
+      power<dim>(
+        lagrange<order>(),
+        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<Vector> >();
+    transferOperators[j]->setMatrix(transferMatrix);
+  }
+
+  mmgStep->setTransferOperators(transferOperators);
+
+  TrustRegionSolver<FEBasis,Vector> solver;
+  solver.setup(*grid,
+               assembler,
+               x,
+               dirichletDofs,
+               tolerance,
+               maxTrustRegionSteps,
+               initialTrustRegionRadius,
+               mmgStep,
+               mgTolerance,
+               multigridIterations);
+
+
+  ///////////////////////////////////////////////////////
+  //   Solve!
+  // /////////////////////////////////////////////////////
+
+  std::cout << "Energy before minimization: " << assembler->computeEnergy(x) << std::endl;
+
+  solver.setInitialIterate(x);
+  solver.solve();
+
+  x = solver.getSol();
+
+  /////////////////////////////////
+  //   Output result
+  /////////////////////////////////
+
+  std::cout << "Energy after minimization: " << assembler->computeEnergy(x) << std::endl;
+
+  writeDisplacement(feBasis, x, F0, resultPath + "quasiconvexity-test-micromorphic-deformation-magic");
 }
 
 int main (int argc, char *argv[]) try
@@ -706,7 +908,15 @@ int main (int argc, char *argv[]) try
   // close to the one we just constructed.
   //////////////////////////////////////////////////////////////////
   auto incompatibleDeformationGradient = Functions::makeDiscreteGlobalBasisFunction<FieldMatrix<double,2,2>>(feBasis, x);
-  constructDeformation(grid, parameterSet, incompatibleDeformationGradient);
+  auto compatibleDeformation = constructDeformation(grid, parameterSet, incompatibleDeformationGradient);
+
+  //////////////////////////////////////////////////////////////////
+  //  Minimize the W_magic energy starting from the compatible deformation
+  //////////////////////////////////////////////////////////////////
+  minimizeEnergyInSpaceOfDeformations(grid,
+                                      parameterSet,
+                                      compatibleDeformation);
+
 } catch (Exception& e) {
     std::cout << e.what() << std::endl;
 }