globalruinanonlinearity.hh

#ifndef DUNE_TECTONIC_GLOBAL_RUINA_NONLINEARITY_HH
#define DUNE_TECTONIC_GLOBAL_RUINA_NONLINEARITY_HH

#include <vector>

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

#include "globalnonlinearity.hh"
#include "localnonlinearity.hh"
#include "nicefunction.hh"

namespace Dune {
template <class MatrixType, class VectorType>
class GlobalRuinaNonlinearity
    : public GlobalNonlinearity<MatrixType, VectorType> {
private:
  typedef shared_ptr<BlockVector<FieldVector<double, 1>> const> dataptr;

public:
  using GlobalNonlinearity<MatrixType, VectorType>::dim;

  GlobalRuinaNonlinearity(dataptr nodalIntegrals, dataptr a, dataptr mu,
                          dataptr eta, dataptr normalStress, dataptr b,
                          dataptr state, dataptr L, double h)
      : nodalIntegrals(nodalIntegrals),
        a(a),
        mu(mu),
        eta(eta),
        normalStress(normalStress),
        b(b),
        state(state),
        L(L),
        h(h),
        trivialNonlinearity(
            new LocalNonlinearity<dim>(make_shared<TrivialFunction const>())),
        restrictions(nodalIntegrals->size()) // TODO: can we get the size from
                                             // another place?
  {
    for (auto &x : restrictions)
      x = shared_ptr<LocalNonlinearity<dim> const>();
  }

  /*
    Return a restriction of the outer function to the i'th node.
  */
  virtual shared_ptr<LocalNonlinearity<dim> const> restriction(int i) const {
    if ((*nodalIntegrals)[i] == 0)
      return trivialNonlinearity;

    if (restrictions[i] != nullptr)
      return restrictions[i];

    auto const func = make_shared<RuinaFunction const>(
        (*nodalIntegrals)[i], (*a)[i], (*mu)[i], (*eta)[i], (*normalStress)[i],
        (*b)[i], (*state)[i], (*L)[i], h);
    restrictions[i] = make_shared<LocalNonlinearity<dim> const>(func);
    return restrictions[i];
  }

private:
  dataptr nodalIntegrals;
  dataptr a;
  dataptr mu;
  dataptr eta;
  dataptr normalStress;
  dataptr b;
  dataptr state;
  dataptr L;
  double const h;

  shared_ptr<LocalNonlinearity<dim> const> const trivialNonlinearity;

  std::vector<shared_ptr<LocalNonlinearity<dim> const>> mutable restrictions;
};
}
#endif