#ifndef DUNE_TECTONIC_GLOBAL_RUINA_NONLINEARITY_HH
#define DUNE_TECTONIC_GLOBAL_RUINA_NONLINEARITY_HH

#include <vector>

#include <dune/common/nullptr.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>

#include "globalnonlinearity.hh"
#include "localnonlinearity.hh"
#include "nicefunction.hh"

namespace Dune {
template <class MatrixType, class VectorType>
class GlobalRuinaNonlinearity
    : public GlobalNonlinearity<MatrixType, VectorType> {
private:
  typedef BlockVector<FieldVector<double, 1>> const &dataref;

public:
  using GlobalNonlinearity<MatrixType, VectorType>::dim;

  GlobalRuinaNonlinearity(dataref nodalIntegrals, dataref a, dataref mu,
                          dataref eta, dataref normalStress, dataref b,
                          dataref state, dataref L, double h)
      : restrictions(nodalIntegrals.size()) {
    auto trivialNonlinearity = make_shared<LocalNonlinearity<dim> const>(
        make_shared<TrivialFunction const>());
    for (size_t i = 0; i < restrictions.size(); ++i) {
      restrictions[i] =
          nodalIntegrals[i] == 0
              ? trivialNonlinearity
              : make_shared<LocalNonlinearity<dim> const>(
                    make_shared<RuinaFunction const>(
                        nodalIntegrals[i], a[i], mu[i], eta[i], normalStress[i],
                        b[i], state[i], L[i], h));
    }
  }

  /*
    Return a restriction of the outer function to the i'th node.
  */
  virtual shared_ptr<LocalNonlinearity<dim> const> restriction(int i) const {
    return restrictions[i];
  }

private:
  std::vector<shared_ptr<LocalNonlinearity<dim> const>> restrictions;
};
}
#endif