Add and Externalise MySolver class (mostly)

src/mysolver.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#ifdef HAVE_IPOPT
+#undef HAVE_IPOPT
+#endif
+
+#include <dune/solvers/common/numproc.hh> // Solver::FULL
+
+#include "mysolver.hh"
+
+template <int dim, class VectorType, class OperatorType, class GridType,
+          class BlockProblemType>
+MySolver<dim, VectorType, OperatorType, GridType, BlockProblemType>::MySolver(
+    Dune::ParameterTree parset, int refinements, double solver_tolerance,
+    GridType const &grid, Dune::BitSetVector<dim> const &ignoreNodes)
+    : baseEnergyNorm(linearBaseSolverStep),
+      linearBaseSolver(&linearBaseSolverStep,
+                       parset.get<int>("linear.maxiterations"),
+                       solver_tolerance, &baseEnergyNorm, Solver::QUIET),
+      transferOperators(refinements),
+      multigridStep(new SolverType(linearIterationStep, nonlinearSmoother)) {
+  // linear iteration step
+  linearIterationStep.setNumberOfLevels(refinements + 1);
+  linearIterationStep.setMGType(parset.get<int>("linear.mu"),
+                                parset.get<int>("linear.nu1"),
+                                parset.get<int>("linear.nu2"));
+  linearIterationStep.basesolver_ = &linearBaseSolver;
+  linearIterationStep.setSmoother(&linearPresmoother, &linearPostsmoother);
+
+  // transfer operators
+  for (auto &x : transferOperators)
+    x = new CompressedMultigridTransfer<VectorType>;
+  TransferOperatorAssembler<GridType>(grid)
+      .assembleOperatorPointerHierarchy(transferOperators);
+  linearIterationStep.setTransferOperators(transferOperators);
+
+  // tnnmg iteration step
+  multigridStep->setSmoothingSteps(parset.get<int>("main.nu1"),
+                                   parset.get<int>("main.mu"),
+                                   parset.get<int>("main.nu2"));
+  multigridStep->ignoreNodes_ = &ignoreNodes;
+}
+
+template <int dim, class VectorType, class OperatorType, class GridType,
+          class BlockProblemType>
+MySolver<dim, VectorType, OperatorType, GridType,
+         BlockProblemType>::~MySolver() {
+  for (auto &x : transferOperators)
+    delete x;
+
+  delete multigridStep;
+}
+
+template <int dim, class VectorType, class OperatorType, class GridType,
+          class BlockProblemType>
+typename MySolver<dim, VectorType, OperatorType, GridType,
+                  BlockProblemType>::SolverType *
+MySolver<dim, VectorType, OperatorType, GridType,
+         BlockProblemType>::getSolver() {
+  return multigridStep;
+}
+
+//#include "mysolver_tmpl.cc"

src/mysolver.hh

+#ifndef MYSOLVER_HH
+#define MYSOLVER_HH
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/parametertree.hh>
+
+#include <dune/solvers/iterationsteps/multigridstep.hh>
+#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
+#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
+#include <dune/fufem/assemblers/transferoperatorassembler.hh>
+#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
+
+template <int dim, class VectorType, class OperatorType, class GridType,
+          class BlockProblemType>
+class MySolver {
+private:
+  typedef GenericNonlinearGS<BlockProblemType> NonlinearSmootherType;
+  typedef TruncatedNonsmoothNewtonMultigrid<BlockProblemType,
+                                            NonlinearSmootherType> SolverType;
+
+public:
+  MySolver(Dune::ParameterTree parset, int refinements, double solver_tolerance,
+           GridType const &grid, Dune::BitSetVector<dim> const &ignoreNodes);
+
+  ~MySolver();
+
+  SolverType *getSolver();
+
+private:
+  TruncatedBlockGSStep<OperatorType, VectorType> linearBaseSolverStep;
+  EnergyNorm<OperatorType, VectorType> baseEnergyNorm;
+  LoopSolver<VectorType> linearBaseSolver;
+  TruncatedBlockGSStep<OperatorType, VectorType> linearPresmoother;
+  TruncatedBlockGSStep<OperatorType, VectorType> linearPostsmoother;
+  MultigridStep<OperatorType, VectorType> linearIterationStep;
+  std::vector<CompressedMultigridTransfer<VectorType> *> transferOperators;
+  NonlinearSmootherType nonlinearSmoother;
+  SolverType *multigridStep;
+};
+
+#endif

src/mysolver_tmpl.cc

+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+#include <dune/grid/yaspgrid.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+#include <dune/tectonic/myblockproblem.hh>
+#include <dune/tectonic/myconvexproblem.hh>
+
+// {{{ 2D
+typedef Dune::FieldVector<double, 2> SmallVector2;
+typedef Dune::FieldMatrix<double, 2, 2> SmallMatrix2;
+typedef Dune::BlockVector<SmallVector2> VectorType2;
+typedef Dune::BCRSMatrix<SmallMatrix2> OperatorType2;
+
+typedef MyConvexProblem<OperatorType2, VectorType2> MyConvexProblemType2;
+typedef MyBlockProblem<MyConvexProblemType2> MyBlockProblemType2;
+
+typedef Dune::YaspGrid<2> GridType2;
+
+template class MySolver<2, VectorType2, OperatorType2, GridType2,
+                        MyBlockProblemType2>;
+// }}}
+
+// {{{ 3D
+typedef Dune::FieldVector<double, 3> SmallVector3;
+typedef Dune::FieldMatrix<double, 3, 3> SmallMatrix3;
+typedef Dune::BlockVector<SmallVector3> VectorType3;
+typedef Dune::BCRSMatrix<SmallMatrix3> OperatorType3;
+
+typedef MyConvexProblem<OperatorType3, VectorType3> MyConvexProblemType3;
+typedef MyBlockProblem<MyConvexProblemType3> MyBlockProblemType3;
+
+typedef Dune::YaspGrid<3> GridType3;
+
+template class MySolver<3, VectorType3, OperatorType3, GridType3,
+                        MyBlockProblemType3>;
+// }}}

src/one-body-sample.cc

@@ -6,10 +6,6 @@
 #error srcdir unset
 #endif
 
-#ifdef HAVE_IPOPT
-#undef HAVE_IPOPT
-#endif
-
 #ifndef HAVE_PYTHON
 #error Python is required
 #endif
@@ -42,20 +38,14 @@
 #include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
 #include <dune/fufem/assemblers/operatorassembler.hh>
-#include <dune/fufem/assemblers/transferoperatorassembler.hh>
 #include <dune/fufem/dunepython.hh>
 #include <dune/fufem/functions/basisgridfunction.hh>
 #include <dune/fufem/functionspacebases/p0basis.hh>
 #include <dune/fufem/functionspacebases/p1nodalbasis.hh>
 #include <dune/fufem/sharedpointermap.hh>
 #include <dune/solvers/common/numproc.hh> // Solver::FULL
-#include <dune/solvers/iterationsteps/multigridstep.hh>
-#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
 #include <dune/solvers/norms/energynorm.hh>
 #include <dune/solvers/solvers/loopsolver.hh>
-#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
-#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
-#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
 
 #include <dune/tectonic/myblockproblem.hh>
 #include <dune/tectonic/myconvexproblem.hh>
@@ -65,6 +55,8 @@
 #include "print.hh"
 #include "vtk.hh"
 
+#include "mysolver.cc" // TODO
+
 int const dim = 2;
 
 template <class GridView>
@@ -139,7 +131,6 @@ int main(int argc, char *argv[]) {
     auto const solver_tolerance = parset.get<double>("solver.tolerance");
 
     auto const refinements = parset.get<size_t>("grid.refinements");
-    size_t const levels = refinements + 1;
 
     auto const verbose = parset.get<bool>("verbose");
     Solver::VerbosityMode const verbosity =
@@ -218,45 +209,11 @@ int main(int argc, char *argv[]) {
     auto const nodalIntegrals =
         assemble_frictional<GridType, GridView, SmallVector, P1Basis>(
             leafView, p1Basis, frictionalNodes);
+
     // {{{ Set up TNNMG solver
-    // linear iteration step components
-    TruncatedBlockGSStep<OperatorType, VectorType> linearBaseSolverStep;
-    EnergyNorm<OperatorType, VectorType> baseEnergyNorm(linearBaseSolverStep);
-    LoopSolver<VectorType> linearBaseSolver(
-        &linearBaseSolverStep,
-        parset.get<int>("solver.tnnmg.linear.maxiterations"), solver_tolerance,
-        &baseEnergyNorm, Solver::QUIET);
-    TruncatedBlockGSStep<OperatorType, VectorType> linearPresmoother;
-    TruncatedBlockGSStep<OperatorType, VectorType> linearPostsmoother;
-
-    // linear iteration step
-    MultigridStep<OperatorType, VectorType> linearIterationStep;
-    linearIterationStep.setNumberOfLevels(levels);
-    linearIterationStep.setMGType(parset.get<int>("solver.tnnmg.linear.mu"),
-                                  parset.get<int>("solver.tnnmg.linear.nu1"),
-                                  parset.get<int>("solver.tnnmg.linear.nu2"));
-    linearIterationStep.basesolver_ = &linearBaseSolver;
-    linearIterationStep.setSmoother(&linearPresmoother, &linearPostsmoother);
-
-    // transfer operators
-    std::vector<CompressedMultigridTransfer<VectorType> *> transferOperators(
-        refinements);
-    for (auto &x : transferOperators)
-      x = new CompressedMultigridTransfer<VectorType>;
-    TransferOperatorAssembler<GridType>(*grid)
-        .assembleOperatorPointerHierarchy(transferOperators);
-    linearIterationStep.setTransferOperators(transferOperators);
-
-    // tnnmg iteration step
-    typedef GenericNonlinearGS<MyBlockProblemType> NonlinearSmootherType;
-    NonlinearSmootherType nonlinearSmoother;
-    TruncatedNonsmoothNewtonMultigrid<MyBlockProblemType, NonlinearSmootherType>
-    multigridStep(linearIterationStep, nonlinearSmoother);
-    multigridStep.setSmoothingSteps(parset.get<int>("solver.tnnmg.main.nu1"),
-                                    parset.get<int>("solver.tnnmg.main.mu"),
-                                    parset.get<int>("solver.tnnmg.main.nu2"));
-    multigridStep.ignoreNodes_ = &ignoreNodes;
-    // }}}
+    MySolver<dim, VectorType, OperatorType, GridType, MyBlockProblemType>
+    mySolver(parset.sub("solver.tnnmg"), refinements, solver_tolerance, *grid,
+             ignoreNodes);
 
     std::fstream octave_writer("data", std::fstream::out);
     timer.reset();
@@ -286,10 +243,11 @@ int main(int argc, char *argv[]) {
                                                     *myGlobalNonlinearity, b4);
 
           MyBlockProblemType myBlockProblem(parset, myConvexProblem);
-          multigridStep.setProblem(u4_diff, myBlockProblem);
+          auto multigridStep = mySolver.getSolver();
+          multigridStep->setProblem(u4_diff, myBlockProblem);
 
           LoopSolver<VectorType> overallSolver(
-              &multigridStep, parset.get<size_t>("solver.tnnmg.maxiterations"),
+              multigridStep, parset.get<size_t>("solver.tnnmg.maxiterations"),
               solver_tolerance, &energyNorm, verbosity);
           overallSolver.solve();
 
@@ -330,10 +288,11 @@ int main(int argc, char *argv[]) {
                                                     *myGlobalNonlinearity, b5);
 
           MyBlockProblemType myBlockProblem(parset, myConvexProblem);
-          multigridStep.setProblem(u5_diff, myBlockProblem);
+          auto multigridStep = mySolver.getSolver();
+          multigridStep->setProblem(u5_diff, myBlockProblem);
 
           LoopSolver<VectorType> overallSolver(
-              &multigridStep, parset.get<size_t>("solver.tnnmg.maxiterations"),
+              multigridStep, parset.get<size_t>("solver.tnnmg.maxiterations"),
               solver_tolerance, &energyNorm, verbosity);
           overallSolver.solve();
 
@@ -376,10 +335,6 @@ int main(int argc, char *argv[]) {
     std::cout << "Making " << timesteps << " time steps took "
               << timer.elapsed() << "s" << std::endl;
 
-    // Clean up after the TNNMG solver
-    for (auto &x : transferOperators)
-      delete x;
-
     octave_writer.close();
 
     if (parset.get<bool>("printFrictionalBoundary")) {