mortar_schur_precond.hpp

//==============================================================================
//==============================================================================
#ifndef SCHUR_PRECOND_HPP_
#define SCHUR_PRECOND_HPP_

#include <dune/istl/preconditioners.hh>
#include <dune/istl/matrixmatrix.hh>
#include <opm/elasticity/matrixops.hpp>
#include <opm/elasticity/mortar_utils.hpp>
#include <opm/elasticity/elasticity_preconditioners.hpp>

namespace Opm {
namespace Elasticity {

  template<class PrecondElasticityBlock>
class MortarSchurPre : public Dune::Preconditioner<Vector,Vector> {
  public:
    MortarSchurPre(const Matrix& P_, const Matrix& B_,
                   PrecondElasticityBlock& Apre_, bool symmetric_ = false) :
      Apre(Apre_), B(B_), N(B.N()), M(B.M()),
      Lpre(P_, false, false), symmetric(symmetric_)
    {
    }

    virtual ~MortarSchurPre()
    {
    }

    void pre(Vector& x, Vector& b) override
    {
      // extract displacement DOFs
      Vector tempx, tempb;
      MortarUtils::extractBlock(tempx, x, N);
      MortarUtils::extractBlock(tempb, b, N);
      Apre.pre(tempx, tempb);
      MortarUtils::injectBlock(x, tempx, N);
      MortarUtils::injectBlock(b, tempb, N);
    }

    void apply(Vector& v, const Vector& d) override
    {
      // multiplier block (second row)
      Vector temp;
      MortarUtils::extractBlock(temp, d, M, N);
      Dune::InverseOperatorResult r;
      Vector temp2;
      temp2.resize(temp.size());
      Lpre.apply(temp2, temp, r);
      MortarUtils::injectBlock(v, temp2, M, N);

      // elasticity block (first row)
      MortarUtils::extractBlock(temp, d, N);
      if (!symmetric)
        B.mmv(temp2, temp);

      temp2.resize(N);
      temp2 = 0;
      Apre.apply(temp2, temp);
      MortarUtils::injectBlock(v, temp2, N);
    }

    void post(Vector& x) override
    {
      Apre.post(x);
    }

    Dune::SolverCategory::Category category() const override
    {
      return Dune::SolverCategory::sequential;
    }

  protected:
    PrecondElasticityBlock& Apre;

    const Matrix& B;

    int N;

    int M;

    LUSolver Lpre;

    bool symmetric;
};

}
}

#endif