use makelinearsolver now instead of constructing locally

dune/tectonic/time-stepping/step.hh

@@ -12,7 +12,7 @@
 #include <dune/solvers/iterationsteps/cgstep.hh>
 #include <dune/solvers/solvers/loopsolver.hh>
 
-#include "../spatial-solving/preconditioners/multilevelpatchpreconditioner.hh"
+#include "../spatial-solving/makelinearsolver.hh"
 
 #include <dune/tectonic/utils/reductionfactors.hh>
 
@@ -130,71 +130,7 @@ class Step : protected StepBase<Factory, ContactNetwork, Updaters, ErrorNorms> {
         using Vector = typename Factory::Vector;
         using Matrix = typename Factory::Matrix;
 
-        /* old, pre multi threading, was unused
-        using Norm =  EnergyNorm<Matrix, Vector>;
-        using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, Matrix, Vector>;
-        using LinearSolver = typename Dune::Solvers::LoopSolver<Vector>;
-
-        const auto& preconditionerParset = parset_.sub("solver.tnnmg.linear.preconditioner");
-
-        Dune::BitSetVector<1> activeLevels(contactNetwork_.nLevels(), true);
-        Preconditioner preconditioner(preconditionerParset, contactNetwork_, activeLevels);
-        preconditioner.setPatchDepth(preconditionerParset.template get<size_t>("patchDepth"));
-        preconditioner.build();
-
-        auto cgStep = std::make_shared<Dune::Solvers::CGStep<Matrix, Vector>>();
-        cgStep->setPreconditioner(preconditioner);
-
-        Norm norm(*cgStep);
-
-        return std::make_shared<LinearSolver>(cgStep, parset_.template get<int>("solver.tnnmg.main.multi"), parset_.template get<double>("solver.tnnmg.linear.tolerance"), norm, Solver::QUIET);
-        */
-
-        // patch preconditioner only needs to be computed once per advance()
-        // make linear solver for linear correction in TNNMGStep
-        using Norm =  EnergyNorm<Matrix, Vector>;
-        using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, Matrix, Vector>;
-        using LinearSolver = typename Dune::Solvers::LoopSolver<Vector>;
-
-        /*const auto& preconditionerParset = parset_.sub("solver.tnnmg.preconditioner");
-
-        Dune::BitSetVector<1> activeLevels(contactNetwork_.nLevels(), true);
-        Preconditioner preconditioner(preconditionerParset, contactNetwork_, activeLevels);
-        preconditioner.setPatchDepth(preconditionerParset.template get<size_t>("patchDepth"));
-        preconditioner.build();
-
-        auto cgStep = std::make_shared<Dune::Solvers::CGStep<Matrix, Vector>>();
-        cgStep->setPreconditioner(preconditioner);
-
-        Norm norm(*cgStep);
-
-        auto linearSolver = std::make_shared<LinearSolver>(cgStep, parset_.template get<int>("solver.tnnmg.main.multi"), parset_.template get<double>("solver.tnnmg.preconditioner.basesolver.tolerance"), norm, Solver::QUIET);
-        */
-
-        // set multigrid solver
-        auto smoother = TruncatedBlockGSStep<Matrix, Vector>();
-
-        // transfer operators need to be recomputed on change due to setDeformation()
-        using TransferOperator = NBodyContactTransfer<ContactNetwork, Vector>;
-        using TransferOperators = std::vector<std::shared_ptr<TransferOperator>>;
-
-        TransferOperators transfer(this->contactNetwork_.nLevels()-1);
-        for (size_t i=0; i<transfer.size(); i++) {
-            transfer[i] = std::make_shared<TransferOperator>();
-            transfer[i]->setup(this->contactNetwork_, i, i+1);
-        }
-
-        // Remove any recompute filed so that initially the full transferoperator is assembled
-        for (size_t i=0; i<transfer.size(); i++)
-            std::dynamic_pointer_cast<TruncatedMGTransfer<Vector> >(transfer[i])->setRecomputeBitField(nullptr);
-
-        auto linearMultigridStep = std::make_shared<Dune::Solvers::MultigridStep<Matrix, Vector> >();
-        linearMultigridStep->setMGType(1, 3, 3);
-        linearMultigridStep->setSmoother(smoother);
-        linearMultigridStep->setTransferOperators(transfer);
-
-        Norm norm(*linearMultigridStep);
-        auto linearSolver = std::make_shared<LinearSolver>(linearMultigridStep, this->parset_.template get<int>("solver.tnnmg.main.multi"), this->parset_.template get<double>("solver.tnnmg.preconditioner.basesolver.tolerance"), norm, Solver::QUIET);
+        auto linearSolver = makeLinearSolver<ContactNetwork, Vector>(this->parset_, this->contactNetwork_);
 
         Vector x;
         x.resize(nBodyAssembler_.totalHasObstacle_.size());