diff --git a/src/quasiconvexity-test.cc b/src/quasiconvexity-test.cc
index 42df972b910993afc1e48ed5f41637c00d78244f..103950f43c8a71b9f52ef2259b85c694c49aa096 100644
--- a/src/quasiconvexity-test.cc
+++ b/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 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);
+ }
- // /////////////////////////////////////////////////
- // Create a Riemannian trust-region solver
- // /////////////////////////////////////////////////
+ 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