mortar_schur_precond.hpp

Go to the documentation of this file.

//==============================================================================
//==============================================================================
#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