#ifndef DUNE_TECTONIC_GLOBALRATESTATEFRICTION_HH
#define DUNE_TECTONIC_GLOBALRATESTATEFRICTION_HH

#include <vector>

#include <dune/common/bitsetvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>

#include <dune/tectonic/geocoordinate.hh>
#include <dune/tectonic/globalfrictiondata.hh>
#include <dune/tectonic/globalfriction.hh>
#include <dune/tectonic/index-in-sorted-range.hh>

template <class Matrix, class Vector, class ScalarFriction, class GridView>
class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
public:
  using GlobalFriction<Matrix, Vector>::block_size;
  using typename GlobalFriction<Matrix, Vector>::LocalNonlinearity;

private:
  using typename GlobalFriction<Matrix, Vector>::ScalarVector;

public:
  GlobalRateStateFriction(BoundaryPatch<GridView> const &frictionalBoundary,
                          GlobalFrictionData<block_size> const &frictionInfo,
                          ScalarVector const &weights,
                          ScalarVector const &weightedNormalStress)
      : restrictions(), localToGlobal(), zeroFriction() {
    auto const gridView = frictionalBoundary.gridView();
    Dune::MultipleCodimMultipleGeomTypeMapper<
        GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView);
    for (auto it = gridView.template begin<block_size>();
         it != gridView.template end<block_size>(); ++it) {
      auto const i = vertexMapper.index(*it);

      if (not frictionalBoundary.containsVertex(i))
        continue;

      localToGlobal.emplace_back(i);
      restrictions.emplace_back(weights[i], weightedNormalStress[i],
                                frictionInfo(geoToPoint(it->geometry())));
    }
    assert(restrictions.size() == frictionalBoundary.numVertices());
    assert(localToGlobal.size() == frictionalBoundary.numVertices());
  }

  void updateAlpha(ScalarVector const &alpha) override {
    for (size_t j = 0; j < restrictions.size(); ++j)
      restrictions[j].updateAlpha(alpha[localToGlobal[j]]);
  }

  /*
    Return a restriction of the outer function to the i'th node.
  */
  LocalNonlinearity const &restriction(size_t i) const override {
    auto const index = indexInSortedRange(localToGlobal, i);
    if (index == localToGlobal.size())
      return zeroFriction;
    return restrictions[index];
  }

private:
  std::vector<WrappedScalarFriction<block_size, ScalarFriction>> restrictions;
  std::vector<size_t> localToGlobal;
  WrappedScalarFriction<block_size, ZeroFunction> const zeroFriction;
};
#endif