My Project
uzawa_solver.hpp
Go to the documentation of this file.
1 #ifndef UZAWA_SOLVER_HPP_
2 #define UZAWA_SOLVER_HPP_
3 //==============================================================================
13 //==============================================================================
14 
15 #include <dune/istl/solvers.hh>
16 
17 namespace Opm {
18  namespace Elasticity {
19 
23  template<class X, class Y>
24 class UzawaSolver : public Dune::InverseOperator<X,Y>
25 {
26  public:
27  typedef std::shared_ptr<Dune::InverseOperator<X,Y> > OperatorPtr;
32  UzawaSolver(OperatorPtr& innersolver_,
33  OperatorPtr& outersolver_,
34  const Matrix& B_) :
35  innersolver(innersolver_), outersolver(outersolver_), B(B_)
36  {
37  }
38 
44  void apply(X& x, Y& b, double /* reduction */,
45  Dune::InverseOperatorResult& res) override
46  {
47  apply(x, b, res);
48  }
49 
54  void apply(X& x, Y& b, Dune::InverseOperatorResult& res) override
55  {
56  Vector lambda, lambda2, u, u2;
57  MortarUtils::extractBlock(u, b, B.N());
58  Dune::InverseOperatorResult res2;
59  u2.resize(u.size());
60  u2 = 0;
61  innersolver->apply(u2, u, res2);
62  lambda.resize(B.M());
63  B.mtv(u2, lambda);
64  lambda2.resize(lambda.size());
65  lambda2 = 0;
66  outersolver->apply(lambda2, lambda, res2);
67  MortarUtils::extractBlock(u, b, B.N());
68  B.usmv(-1.0, lambda2, u);
69  u2 = 0;
70  innersolver->apply(u2, u, res);
71  MortarUtils::injectBlock(x, u2, B.N());
72  MortarUtils::injectBlock(x, lambda2, B.M(), B.N());
73  }
74 
75  Dune::SolverCategory::Category category() const override
76  {
77  return Dune::SolverCategory::sequential;
78  }
79 
80  protected:
81  OperatorPtr innersolver;
82  OperatorPtr outersolver;
83  const Matrix& B;
84 };
85 
86 }
87 }
88 
89 #endif
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
Template implementing an Uzawa scheme (block Gaussian-elimination) for a (symmetric indefinite) saddl...
Definition: uzawa_solver.hpp:25
OperatorPtr outersolver
The outer solver.
Definition: uzawa_solver.hpp:82
const Matrix & B
The coupling matrix.
Definition: uzawa_solver.hpp:83
void apply(X &x, Y &b, double, Dune::InverseOperatorResult &res) override
Apply the scheme to a vector.
Definition: uzawa_solver.hpp:44
void apply(X &x, Y &b, Dune::InverseOperatorResult &res) override
Apply the scheme to a vector.
Definition: uzawa_solver.hpp:54
UzawaSolver(OperatorPtr &innersolver_, OperatorPtr &outersolver_, const Matrix &B_)
Default constructor.
Definition: uzawa_solver.hpp:32
OperatorPtr innersolver
The inner solver.
Definition: uzawa_solver.hpp:81
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
Inverting small matrices.
Definition: ImplicitAssembly.hpp:43