globalfriction.hh

#ifndef DUNE_TECTONIC_GLOBALFRICTION_HH
#define DUNE_TECTONIC_GLOBALFRICTION_HH

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

#include <dune/solvers/common/interval.hh>

#include "localfriction.hh"

template <class Matrix, class Vector>
class GlobalFriction {
protected:
  using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;

public:
  using IndexSet = Dune::MatrixIndexSet;
  using MatrixType = Matrix;
  using VectorType = Vector;
  using LocalMatrix = typename Matrix::block_type;
  using LocalVectorType = typename Vector::block_type;
  size_t static const block_size = LocalVectorType::dimension;
  using LocalNonlinearity = LocalFriction<block_size>;

  double operator()(Vector const &x) const {
    double tmp = 0;
    for (size_t i = 0; i < x.size(); ++i) {
      tmp += restriction(i)(x[i]);
    }
    return tmp;
  }

  /*
    Return a restriction of the outer function to the i'th node.
  */
  LocalNonlinearity const virtual &restriction(size_t i) const = 0;

  void addHessian(Vector const &v, Matrix &hessian) const {
    for (size_t i = 0; i < v.size(); ++i) {
      restriction(i).addHessian(v[i], hessian[i][i]);
    }
  }

  void directionalDomain(Vector const &, Vector const &,
                         Dune::Solvers::Interval<double> &dom) const {
    dom[0] = -std::numeric_limits<double>::max();
    dom[1] = std::numeric_limits<double>::max();
  }

  void directionalSubDiff(
      Vector const &u, Vector const &v,
      Dune::Solvers::Interval<double> &subdifferential) const {
    subdifferential[0] = subdifferential[1] = 0;
    for (size_t i = 0; i < u.size(); ++i) {
      Dune::Solvers::Interval<double> D;
      restriction(i).directionalSubDiff(u[i], v[i], D);
      subdifferential[0] += D[0];
      subdifferential[1] += D[1];
    }
  }

  void addHessianIndices(Dune::MatrixIndexSet &indices) const {
    for (size_t i = 0; i < indices.rows(); ++i)
      indices.add(i, i);
  }

  void addGradient(Vector const &v, Vector &gradient) const {
    for (size_t i = 0; i < v.size(); ++i)
      restriction(i).addGradient(v[i], gradient[i]);
  }

  double regularity(size_t i, typename Vector::block_type const &x) const {
    return restriction(i).regularity(x);
  }

  ScalarVector coefficientOfFriction(Vector const &x) const {
    ScalarVector ret(x.size());
    for (size_t i = 0; i < x.size(); ++i) {
      ret[i] = restriction(i).coefficientOfFriction(x[i]);
    }
    return ret;
  }

  void virtual updateAlpha(const std::vector<ScalarVector>& alpha) = 0;

};
#endif