#ifndef DUNE_TECTONIC_SPATIAL_SOLVING_MAKELINEARSOLVER_HH
#define DUNE_TECTONIC_SPATIAL_SOLVING_MAKELINEARSOLVER_HH

#include <dune/common/parametertree.hh>

#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/iterationsteps/multigridstep.hh>
#include <dune/solvers/iterationsteps/cgstep.hh>
#include <dune/solvers/solvers/loopsolver.hh>

#include "preconditioners/multilevelpatchpreconditioner.hh"

template <class ContactNetwork, class VectorType>
auto makeLinearSolver(const Dune::ParameterTree& parset, const ContactNetwork& contactNetwork) {
    // make linear solver for linear correction in TNNMGStep
    using Vector = VectorType;
    using Matrix = typename ContactNetwork::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);
        */




        // 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(contactNetwork.nLevels()-1);
        for (size_t i=0; i<transfer.size(); i++) {
            transfer[i] = std::make_shared<TransferOperator>();
            transfer[i]->setup(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(TruncatedBlockGSStep<Matrix, Vector>());
        linearMultigridStep->setTransferOperators(transfer);

        return std::make_shared<LinearSolver>(linearMultigridStep, parset.template get<int>("solver.tnnmg.main.multi"), parset.template get<double>("solver.tnnmg.preconditioner.basesolver.tolerance"), Norm(*linearMultigridStep), Solver::QUIET);
    }
#endif