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:
    // define the category
    enum {
      category=Dune::SolverCategory::sequential
    };

    MortarEvaluator(const Matrix& A_,
                    const Matrix& B_) :
      A(A_), B(B_)
    {
    }

    void apply(const Vector& x, Vector& y) const
    {
      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
    {
      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()];
    }
  protected:
    const Matrix& A; 
    const Matrix& B; 
};

}
}

#endif