mortar_evaluator.hpp

Go to the documentation of this file.

//==============================================================================
//==============================================================================
#ifndef MORTAR_EVALUATOR_HPP_
#define MORTAR_EVALUATOR_HPP_
 
#include <dune/istl/matrixmatrix.hh>
#include <dune/istl/solvers.hh>
#include <opm/elasticity/matrixops.hpp>
#include <opm/elasticity/mortar_utils.hpp>
 
namespace Opm {
namespace Elasticity {
 
class MortarEvaluator : public Dune::LinearOperator<Vector, Vector> {
  public:
    MortarEvaluator(const Matrix& A_,
                    const Matrix& B_) :
      A(A_), B(B_)
    {
    }
 
    void apply(const Vector& x, Vector& y) const override
    {
      Vector lambda, l2;
 
      MortarUtils::extractBlock(lambda, x, B.M(), A.N());
      A.mv(x, y);
      B.umv(lambda, y);
      l2.resize(lambda.size());
      B.mtv(x, l2);
      MortarUtils::injectBlock(y, l2, B.M(), A.N());
    }
 
    void applyscaleadd(field_type alpha, const Vector& x, Vector& y) const override
    {
      Vector lambda, l2;
 
      A.usmv(alpha, x, y);
      MortarUtils::extractBlock(lambda, x, B.M(), A.N());
      B.umv(lambda, y);
      l2.resize(lambda.size());
      B.mtv(x, l2);
      for (size_t i=A.N(); i < y.size(); ++i)
        y[i] += alpha*l2[i-A.N()];
    }
 
    Dune::SolverCategory::Category category() const override
    {
      return Dune::SolverCategory::sequential;
    }
 
  protected:
    const Matrix& A; 
    const Matrix& B; 
};
 
}
}
 
#endif