#ifndef ASSEMBLERS_HH
#define ASSEMBLERS_HH

#include <dune/common/bitsetvector.hh>
#include <dune/common/function.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>

#include <dune/fufem/assemblers/assembler.hh>
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>

#include <dune/tectonic/globalnonlinearity.hh>

template <class GridView, int dimension> class MyAssembler {
public:
  using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
  using Matrix =
      Dune::BCRSMatrix<Dune::FieldMatrix<double, dimension, dimension>>;
  using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
  using Vector = Dune::BlockVector<Dune::FieldVector<double, dimension>>;

  using VertexBasis = P1NodalBasis<GridView, double>;
  VertexBasis vertexBasis;

private:
  using Grid = typename GridView::Grid;
  using LocalVector = typename Vector::block_type;
  using LocalVertexBasis = typename VertexBasis::LocalFiniteElement;

  GridView const &gridView;
  Assembler<VertexBasis, VertexBasis> vertexAssembler;

public:
  MyAssembler(GridView const &gridView);

  void assembleFrictionalBoundaryMass(
      BoundaryPatch<GridView> const &frictionalBoundary,
      ScalarMatrix &frictionalBoundaryMass);

  void assembleMass(double density, Matrix &M);

  void assembleElasticity(double E, double nu, Matrix &A);

  void assembleViscosity(double shearViscosity, double bulkViscosity,
                         Matrix &C);

  void assembleBodyForce(double gravity, double density, LocalVector zenith,
                         Vector &f);

  void assembleNeumann(BoundaryPatch<GridView> const &neumannBoundary,
                       Vector &f,
                       Dune::VirtualFunction<double, double> const &neumann,
                       double relativeTime);

  std::shared_ptr<GlobalNonlinearity<Matrix, Vector>>
  assembleFrictionNonlinearity(
      BoundaryPatch<GridView> const &frictionalBoundary,
      FrictionData const &fd);
};

#endif