#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/functions/constantfunction.hh>
#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
#include <dune/fufem/assemblers/boundaryfunctionalassembler.hh>

#include <dune/tectonic/globalruinanonlinearity.hh>

#include "assemblers.hh"

// Assembles Neumann boundary term in f
template <class GridType, class GridView, class LocalVectorType, class FEBasis>
void assemble_neumann(GridView const &gridView, FEBasis const &feBasis,
                      Dune::BitSetVector<1> const &neumannNodes,
                      Dune::BlockVector<LocalVectorType> &f,
                      Dune::VirtualFunction<double, double> const &neumann,
                      double time) { // constant sample function on neumann
                                     // boundary
  BoundaryPatch<GridView> const neumannBoundary(gridView, neumannNodes);
  LocalVectorType SampleVector(0);
  neumann.evaluate(time, SampleVector[0]);
  SampleVector[1] = 0;
  ConstantFunction<LocalVectorType, LocalVectorType> const fNeumann(
      SampleVector);
  NeumannBoundaryAssembler<GridType, LocalVectorType> neumannBoundaryAssembler(
      fNeumann);
  BoundaryFunctionalAssembler<FEBasis>(feBasis, neumannBoundary)
      .assemble(neumannBoundaryAssembler, f, true);
}

// Assembles constant 1-function on frictional boundary in nodalIntegrals
template <class GridType, class GridView, class LocalVectorType, class FEBasis>
Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
assemble_frictional(GridView const &gridView, FEBasis const &feBasis,
                    Dune::BitSetVector<1> const &frictionalNodes) {
  typedef Dune::FieldVector<double, 1> Singleton;
  BoundaryPatch<GridView> const frictionalBoundary(gridView, frictionalNodes);
  ConstantFunction<LocalVectorType, Singleton> const constantOneFunction(1);
  NeumannBoundaryAssembler<GridType, Singleton> frictionalBoundaryAssembler(
      constantOneFunction);

  auto const nodalIntegrals = Dune::make_shared<Dune::BlockVector<Singleton>>();
  BoundaryFunctionalAssembler<FEBasis>(feBasis, frictionalBoundary)
      .assemble(frictionalBoundaryAssembler, *nodalIntegrals, true);
  return nodalIntegrals;
}

template <class MatrixType, class VectorType>
Dune::shared_ptr<Dune::GlobalNonlinearity<MatrixType, VectorType> const>
assemble_nonlinearity(
    Dune::ParameterTree const &parset,
    Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
    Dune::BlockVector<Dune::FieldVector<double, 1>> const &state,
    Dune::BlockVector<Dune::FieldVector<double, 1>> const &normalStress) {
  auto const size = nodalIntegrals.size();

  typedef Dune::BlockVector<Dune::FieldVector<double, 1>> SingletonVectorType;
  // {{{ Assemble terms for the nonlinearity
  SingletonVectorType mu(size);
  mu = parset.get<double>("mu");

  SingletonVectorType a(size);
  a = parset.get<double>("a");

  SingletonVectorType V0(size);
  V0 = parset.get<double>("V0");

  SingletonVectorType b(size);
  b = parset.get<double>("b");

  SingletonVectorType L(size);
  L = parset.get<double>("L");

  return Dune::make_shared<
      Dune::GlobalRuinaNonlinearity<MatrixType, VectorType> const>(
      nodalIntegrals, a, mu, V0, normalStress, b, state, L);
}

#include "assemblers_tmpl.cc"