#ifndef DUNE_TECTONIC_SPATIAL_SOLVING_GLOBALPROBLEM_HH
#define DUNE_TECTONIC_SPATIAL_SOLVING_GLOBALPROBLEM_HH

#include <vector>
#include <memory>

#include "tnnmg/functional.hh"

template <class NBodyAssembler, class Nonlinearity, class MatrixType, class VectorType>
class FunctionalFactory {
    using Vector = VectorType;
    using Matrix = MatrixType;

    const static int dims = Vector::block_type::dimension;

public:
    using Functional = Functional<Matrix&, Vector&, Nonlinearity&, Vector&, Vector&, typename Matrix::field_type>;

    FunctionalFactory(const NBodyAssembler& nBodyAssembler, const Nonlinearity& nonlinearity, const std::vector<Matrix>& matrices, const std::vector<Vector>& rhs) :
        nBodyAssembler_(nBodyAssembler),
        nonlinearity_(nonlinearity),
        rhs_(rhs) {

        const auto nBodies = nBodyAssembler_.nGrids();

        matrices_.resize(nBodies);
        for (int i=0; i<nBodies; i++) {
            matrices_[i] = &matrices[i];
        }
    }

    FunctionalFactory(const NBodyAssembler& nBodyAssembler, const Nonlinearity& nonlinearity, const std::vector<std::shared_ptr<Matrix>>& matrices, const std::vector<Vector>& rhs) :
        nBodyAssembler_(nBodyAssembler),
        nonlinearity_(nonlinearity),
        rhs_(rhs) {

        const auto nBodies = nBodyAssembler_.nGrids();

        matrices_.resize(nBodies);
        for (size_t i=0; i<nBodies; i++) {
            matrices_[i] = matrices[i].get();
      }
    }

    void build() {
        // assemble full global contact problem
        nBodyAssembler_.assembleJacobian(matrices_, bilinearForm_);
        nBodyAssembler_.assembleRightHandSide(rhs_, totalRhs_);

        // get lower and upper obstacles
        const auto& totalObstacles = nBodyAssembler_.totalObstacles_;
        lower_.resize(totalObstacles.size());
        upper_.resize(totalObstacles.size());

        for (size_t j=0; j<totalObstacles.size(); ++j) {
            const auto& totalObstaclesj = totalObstacles[j];
            auto& lowerj = lower_[j];
            auto& upperj = upper_[j];
          for (size_t d=0; d<dims; ++d) {
              lowerj[d] = totalObstaclesj[d][0];
              upperj[d] = totalObstaclesj[d][1];
          }
        }

        // set up functional
        functional_ = std::make_shared<Functional>(bilinearForm_, totalRhs_, nonlinearity_, lower_, upper_);
    }

    std::shared_ptr<Functional> functional() const {
        return functional_;
    }

    const Matrix& bilinearForm() const {
        return bilinearForm_;
    }

    const Vector& rhs() const {
        return totalRhs_;
    }

    const Vector& lower() const {
        return lower_;
    }

    const Vector& upper() const {
        return upper_;
    }

private:
    const NBodyAssembler nBodyAssembler_;
    const Nonlinearity& nonlinearity_;
    std::vector<const Matrix*> matrices_;
    const std::vector<Vector>& rhs_;

    Matrix bilinearForm_;
    Vector totalRhs_;
    Vector lower_;
    Vector upper_;

    std::shared_ptr<Functional> functional_;
};
#endif