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

#include <dune/istl/scaledidmatrix.hh>

#include <dune/fufem/assemblers/localassemblers/boundarymassassembler.hh>
#include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
#include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
#include <dune/fufem/assemblers/localassemblers/viscosityassembler.hh>
#include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
#include <dune/fufem/assemblers/localassemblers/weightedmassassembler.hh>
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/functions/basisgridfunction.hh>
#include <dune/fufem/functions/constantfunction.hh>

#include "assemblers.hh"

template <class GridView, int dimension>
MyAssembler<GridView, dimension>::MyAssembler(GridView const &_gridView)
    : cellBasis(_gridView),
      vertexBasis(_gridView),
      gridView(_gridView),
      cellAssembler(cellBasis, cellBasis),
      vertexAssembler(vertexBasis, vertexBasis) {}

template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleFrictionalBoundaryMass(
    BoundaryPatch<GridView> const &frictionalBoundary,
    ScalarMatrix &frictionalBoundaryMass) {
  BoundaryMassAssembler<Grid, BoundaryPatch<GridView>, LocalVertexBasis,
                        LocalVertexBasis, Dune::FieldMatrix<double, 1, 1>> const
  frictionalBoundaryMassAssembler(frictionalBoundary);
  vertexAssembler.assembleOperator(frictionalBoundaryMassAssembler,
                                   frictionalBoundaryMass);
}

template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleMass(
    Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
        densityFunction,
    Matrix &M) {
  // NOTE: We treat the weight as a constant function
  QuadratureRuleKey quadKey(dimension, 0);

  WeightedMassAssembler<Grid, LocalVertexBasis, LocalVertexBasis,
                        Dune::VirtualFunction<LocalVector, LocalScalarVector>,
                        Dune::ScaledIdentityMatrix<double, dimension>>
  localWeightedMass(gridView.grid(), densityFunction, quadKey);
  vertexAssembler.assembleOperator(localWeightedMass, M);
}

template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleElasticity(double E, double nu,
                                                          Matrix &A) {
  StVenantKirchhoffAssembler<Grid, LocalVertexBasis, LocalVertexBasis> const
  localStiffness(E, nu);
  vertexAssembler.assembleOperator(localStiffness, A);
}

template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleViscosity(double shearViscosity,
                                                         double bulkViscosity,
                                                         Matrix &C) {
  ViscosityAssembler<Grid, LocalVertexBasis, LocalVertexBasis> const
  localViscosity(shearViscosity, bulkViscosity);
  vertexAssembler.assembleOperator(localViscosity, C);
}

template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleBodyForce(
    Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
    Vector &f) {
  L2FunctionalAssembler<Grid, LocalVertexBasis, LocalVector>
  gravityFunctionalAssembler(gravityField);
  vertexAssembler.assembleFunctional(gravityFunctionalAssembler, f);
}

template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleNeumann(
    BoundaryPatch<GridView> const &neumannBoundary, Vector &f,
    Dune::VirtualFunction<double, double> const &neumann, double relativeTime) {
  LocalVector localNeumann(0);
  neumann.evaluate(relativeTime, localNeumann[0]);
  ConstantFunction<LocalVector, LocalVector> const fNeumann(localNeumann);
  NeumannBoundaryAssembler<Grid, LocalVector> neumannBoundaryAssembler(
      fNeumann);
  vertexAssembler.assembleBoundaryFunctional(neumannBoundaryAssembler, f,
                                             neumannBoundary);
}

template <class GridView, int dimension>
auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
    BoundaryPatch<GridView> const &frictionalBoundary, FrictionData const &fd)
    -> std::shared_ptr<GlobalRuinaNonlinearity<Matrix, Vector>> {
  // Lump negative normal stress (kludge)
  ScalarVector weights;
  {
    ConstantFunction<LocalVector, typename ScalarVector::block_type> const
    constantOneFunction(1);
    NeumannBoundaryAssembler<Grid, typename ScalarVector::block_type>
    frictionalBoundaryAssembler(constantOneFunction);
    vertexAssembler.assembleBoundaryFunctional(frictionalBoundaryAssembler,
                                               weights, frictionalBoundary);
  }
  for (size_t i = 0; i < weights.size(); ++i)
    assert(weights[i] >= 0.0);

  Dune::BitSetVector<1> frictionalNodes;
  frictionalBoundary.getVertices(frictionalNodes);
  return std::make_shared<GlobalRuinaNonlinearity<Matrix, Vector>>(
      frictionalNodes, weights, fd);
}

template <class GridView, int dimension>
void MyAssembler<GridView, dimension>::assembleVonMisesStress(
    double youngModulus, double poissonRatio, Vector const &u,
    ScalarVector &stress) {
  auto const gridDisplacement =
      std::make_shared<BasisGridFunction<VertexBasis, Vector> const>(
          vertexBasis, u);
  VonMisesStressAssembler<Grid, LocalCellBasis> localStressAssembler(
      youngModulus, poissonRatio, gridDisplacement);
  cellAssembler.assembleFunctional(localStressAssembler, stress);
}

#include "assemblers_tmpl.cc"