diff --git a/src/quasiconvexity-test-micromorphic.cc b/src/quasiconvexity-test-micromorphic.cc
index f4df37024a429100b164102f514fb868a8c90b39..c2abca48cfe40a1e34aa242b9c8a37e58434c4b1 100644
--- a/src/quasiconvexity-test-micromorphic.cc
+++ b/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;
}