My Project
mortar_schur_precond.hpp
1 //==============================================================================
11 //==============================================================================
12 #ifndef SCHUR_PRECOND_HPP_
13 #define SCHUR_PRECOND_HPP_
14 
15 #include <dune/istl/preconditioners.hh>
16 #include <dune/istl/matrixmatrix.hh>
20 
21 namespace Opm {
22 namespace Elasticity {
23 
36  template<class PrecondElasticityBlock>
37 class MortarSchurPre : public Dune::Preconditioner<Vector,Vector> {
38  public:
44  MortarSchurPre(const Matrix& P_, const Matrix& B_,
45  PrecondElasticityBlock& Apre_, bool symmetric_=false) :
46  Apre(Apre_), B(B_), N(B.N()), M(B.M()),
47  Lpre(P_, false, false), symmetric(symmetric_)
48  {
49  }
50 
52  virtual ~MortarSchurPre()
53  {
54  }
55 
57  void pre(Vector& x, Vector& b) override
58  {
59  // extract displacement DOFs
60  Vector tempx, tempb;
61  MortarUtils::extractBlock(tempx, x, N);
62  MortarUtils::extractBlock(tempb, b, N);
63  Apre.pre(tempx, tempb);
64  MortarUtils::injectBlock(x, tempx, N);
65  MortarUtils::injectBlock(b, tempb, N);
66  }
67 
71  void apply(Vector& v, const Vector& d) override
72  {
73  // multiplier block (second row)
74  Vector temp;
75  MortarUtils::extractBlock(temp, d, M, N);
76  Dune::InverseOperatorResult r;
77  Vector temp2;
78 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
79  temp2.resize(temp.size());
80 #else
81  temp2.resize(temp.size(), false);
82 #endif
83  Lpre.apply(temp2, temp, r);
84  MortarUtils::injectBlock(v, temp2, M, N);
85 
86  // elasticity block (first row)
87  MortarUtils::extractBlock(temp, d, N);
88  if (!symmetric)
89  B.mmv(temp2, temp);
90 
91 #if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 7)
92  temp2.resize(N);
93 #else
94  temp2.resize(N, false);
95 #endif
96  temp2 = 0;
97  Apre.apply(temp2, temp);
98  MortarUtils::injectBlock(v, temp2, N);
99  }
100 
102  void post(Vector& x) override
103  {
104  Apre.post(x);
105  }
106 
107  Dune::SolverCategory::Category category() const override
108  {
109  return Dune::SolverCategory::sequential;
110  }
111 
112  protected:
114  PrecondElasticityBlock& Apre;
115 
117  const Matrix& B;
118 
120  int N;
121 
123  int M;
124 
126  LUSolver Lpre;
127 
129  bool symmetric;
130 };
131 
132 }
133 }
134 
135 #endif
This implements a Schur-decomposition based preconditioner for the mortar-elasticity system [A B] [B'...
Definition: mortar_schur_precond.hpp:37
void post(Vector &x) override
Dummy post-process function.
Definition: mortar_schur_precond.hpp:102
virtual ~MortarSchurPre()
Destructor.
Definition: mortar_schur_precond.hpp:52
PrecondElasticityBlock & Apre
The preconditioner for the elasticity operator.
Definition: mortar_schur_precond.hpp:114
int M
Number of multiplier DOFs.
Definition: mortar_schur_precond.hpp:123
void pre(Vector &x, Vector &b) override
Preprocess preconditioner.
Definition: mortar_schur_precond.hpp:57
const Matrix & B
The mortar coupling matrix.
Definition: mortar_schur_precond.hpp:117
bool symmetric
Whether or not to use a symmetric preconditioner.
Definition: mortar_schur_precond.hpp:129
LUSolver Lpre
Linear solver for the multiplier block.
Definition: mortar_schur_precond.hpp:126
MortarSchurPre(const Matrix &P_, const Matrix &B_, PrecondElasticityBlock &Apre_, bool symmetric_=false)
Constructor.
Definition: mortar_schur_precond.hpp:44
int N
Number of displacement DOFs.
Definition: mortar_schur_precond.hpp:120
void apply(Vector &v, const Vector &d) override
Applies the preconditioner.
Definition: mortar_schur_precond.hpp:71
static void extractBlock(Vector &x, const Vector &y, int len, int start=0)
Extract a range of indices from a vector.
Definition: mortar_utils.hpp:25
static void injectBlock(Vector &x, const Vector &y, int len, int start=0)
Inject a range of indices into a vector.
Definition: mortar_utils.hpp:40
Preconditioners for elasticity upscaling.
Helper class with some matrix operations.
Dune::BCRSMatrix< Dune::FieldMatrix< double, 1, 1 > > Matrix
A sparse matrix holding our operator.
Definition: matrixops.hpp:27
Dune::BlockVector< Dune::FieldVector< double, 1 > > Vector
A vector holding our RHS.
Definition: matrixops.hpp:33
Mortar helper class.
Inverting small matrices.
Definition: ImplicitAssembly.hpp:43