MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
constraints.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #ifndef MFEM_CONSTRAINED
13 #define MFEM_CONSTRAINED
14 
15 #include "solvers.hpp"
16 #include "blockoperator.hpp"
17 #include "sparsemat.hpp"
18 #include "hypre.hpp"
19 
20 namespace mfem
21 {
22 
23 class FiniteElementSpace;
24 #ifdef MFEM_USE_MPI
25 class ParFiniteElementSpace;
26 #endif
27 
28 /** @brief An abstract class to solve the constrained system \f$ Ax = f \f$
29  subject to the constraint \f$ B x = r \f$.
30 
31  Although implementations may not use the below formulation, for
32  understanding some of its methods and notation you can think of
33  it as solving the saddle-point system
34 
35  ( A B^T ) ( x ) ( f )
36  ( B ) ( lambda ) = ( r )
37 
38  Do not confuse with ConstrainedOperator, which handles only simple
39  pointwise constraints and is not a Solver.
40 
41  The height and width of this object as an IterativeSolver are the same as
42  just the unconstrained operator \f$ A \f$, and the Mult() interface just
43  takes \f$ f \f$ as an argument. You can set \f$ r \f$ with
44  SetConstraintRHS() (it defaults to zero) and get the Lagrange multiplier
45  solution with GetMultiplierSolution().
46 
47  Alternatively, you can use LagrangeSystemMult() to solve the block system
48  shown above.
49 
50  This abstract object unifies this interface so that derived classes can
51  solve whatever linear system makes sense and the interface will provide
52  uniform access to solutions, Lagrange multipliers, etc. */
54 {
55 public:
56 #ifdef MFEM_USE_MPI
57  ConstrainedSolver(MPI_Comm comm, Operator& A_, Operator& B_);
58 #endif
60 
61  virtual ~ConstrainedSolver() { }
62 
63  virtual void SetOperator(const Operator& op) override { }
64 
65  /** @brief Set the right-hand side r for the constraint B x = r
66 
67  (r defaults to zero if you don't call this) */
68  virtual void SetConstraintRHS(const Vector& r);
69 
70  /** @brief Return the Lagrange multiplier solution in lambda
71 
72  Mult() only gives you x, this provides access to lambda
73 
74  Does not make sense unless you've already solved the constrained
75  system with Mult() or LagrangeSystemMult() */
76  void GetMultiplierSolution(Vector& lambda) const { lambda = multiplier_sol; }
77 
78  /** @brief Solve for \f$ x \f$ given \f$ f \f$.
79 
80  If you want to set \f$ r \f$, call SetConstraintRHS() before this.
81 
82  If you want to get \f$ \lambda \f$, call GetMultiplierSolution() after
83  this.
84 
85  The base class implementation calls LagrangeSystemMult(), so derived
86  classes must implement either this or LagrangeSystemMult() */
87  virtual void Mult(const Vector& f, Vector& x) const override;
88 
89  /** @brief Solve for (x, lambda) given (f, r)
90 
91  The base class implementation calls Mult(), so derived classes
92  must implement either this or Mult() */
93  virtual void LagrangeSystemMult(const Vector& f_and_r,
94  Vector& x_and_lambda) const;
95 
96 protected:
99 
102  mutable Vector workb;
103  mutable Vector workx;
104 
105 private:
106  void Initialize();
107 };
108 
109 
110 /** @brief Perform elimination of a single constraint.
111 
112  See EliminationProjection, EliminationCGSolver
113 
114  This keeps track of primary / secondary tdofs and does small dense block
115  solves to eliminate constraints from a global system.
116 
117  \f$ B_s^{-1} \f$ maps the lagrange space into secondary dofs, while
118  \f$ -B_s^{-1} B_p \f$ maps primary dofs to secondary dofs. */
120 {
121 public:
122  Eliminator(const SparseMatrix& B, const Array<int>& lagrange_dofs,
123  const Array<int>& primary_tdofs,
124  const Array<int>& secondary_tdofs);
125 
126  const Array<int>& LagrangeDofs() const { return lagrange_tdofs; }
127  const Array<int>& PrimaryDofs() const { return primary_tdofs; }
128  const Array<int>& SecondaryDofs() const { return secondary_tdofs; }
129 
130  /// Given primary dofs in in, return secondary dofs in out
131  /// This applies \f$ -B_s^{-1} B_p \f$.
132  void Eliminate(const Vector& in, Vector& out) const;
133 
134  /// Transpose of Eliminate(), applies \f$ -B_p^T B_s^{-T} \f$
135  void EliminateTranspose(const Vector& in, Vector& out) const;
136 
137  /// Maps Lagrange multipliers to secondary dofs, applies \f$ B_s^{-1} \f$
138  void LagrangeSecondary(const Vector& in, Vector& out) const;
139 
140  /// Transpose of LagrangeSecondary()
141  void LagrangeSecondaryTranspose(const Vector& in, Vector& out) const;
142 
143  /// Return \f$ -B_s^{-1} B_p \f$ explicitly assembled in mat
144  void ExplicitAssembly(DenseMatrix& mat) const;
145 
146 private:
147  Array<int> lagrange_tdofs;
148  Array<int> primary_tdofs;
149  Array<int> secondary_tdofs;
150 
151  DenseMatrix Bp;
152  DenseMatrix Bs; // gets inverted in place
153  LUFactors Bsinverse;
154  DenseMatrix BsT; // gets inverted in place
155  LUFactors BsTinverse;
156  Array<int> ipiv;
157  Array<int> ipivT;
158 };
159 
160 
161 /** Collects action of several Eliminator objects to perform elimination of
162  constraints.
163 
164  Works in parallel, but each Eliminator must be processor local, and must
165  operate on disjoint degrees of freedom (ie, the primary and secondary dofs
166  for one Eliminator must have no overlap with any dofs from a different
167  Eliminator). */
169 {
170 public:
171  EliminationProjection(const Operator& A, Array<Eliminator*>& eliminators);
172 
173  void Mult(const Vector& x, Vector& y) const override;
174 
175  void MultTranspose(const Vector& x, Vector& y) const override;
176 
177  /** @brief Assemble this projector as a (processor-local) SparseMatrix.
178 
179  Some day we may also want to try approximate variants. */
180  SparseMatrix * AssembleExact() const;
181 
182  /** Given Lagrange multiplier right-hand-side \f$ g \f$, return
183  \f$ \tilde{g} \f$ */
184  void BuildGTilde(const Vector& g, Vector& gtilde) const;
185 
186  /** After a solve, recover the Lagrange multiplier. */
187  void RecoverMultiplier(const Vector& primalrhs,
188  const Vector& primalvars, Vector& lm) const;
189 
190 private:
191  const Operator& Aop;
192  Array<Eliminator*> eliminators;
193 };
194 
195 
196 #ifdef MFEM_USE_MPI
197 /** @brief Solve constrained system by eliminating the constraint; see
198  ConstrainedSolver
199 
200  Solves the system with the operator \f$ P^T A P + Z_P \f$, where P is
201  EliminationProjection and Z_P is the identity on the eliminated dofs. */
203 {
204 public:
205  /** @brief Constructor, with explicit splitting into primary/secondary dofs.
206 
207  This constructor uses a single elimination block (per processor), which
208  provides the most general algorithm but is also not scalable
209 
210  The secondary_dofs are eliminated from the system in this algorithm,
211  as they can be written in terms of the primary_dofs.
212 
213  Both primary_dofs and secondary_dofs are in the local truedof numbering;
214  All elimination has to be done locally on processor, though the global
215  system can be parallel. */
217  Array<int>& primary_dofs,
218  Array<int>& secondary_dofs);
219 
220  /** @brief Constructor, elimination is by blocks.
221 
222  Each block is eliminated independently; if the blocks are reasonably
223  small this can be reasonably efficient.
224 
225  The nonzeros in B are assumed to be in disjoint rows and columns; the
226  rows are identified with the constraint_rowstarts array, the secondary
227  dofs are assumed to be the first nonzeros in the rows. */
229  Array<int>& constraint_rowstarts);
230 
232 
233  void Mult(const Vector& x, Vector& y) const override;
234 
235  void SetOperator(const Operator& op) override
236  { MFEM_ABORT("Operator cannot be reset!"); }
237 
238  void SetPreconditioner(Solver& precond) override
239  { prec = &precond; }
240 
241 protected:
242  /// Internal utility routine; assembles eliminated matrix explicitly
243  void BuildExplicitOperator();
244 
245  /// Build preconditioner for eliminated system
246  virtual Solver* BuildPreconditioner() const = 0;
247  /// Select krylov solver for eliminated system
248  virtual IterativeSolver* BuildKrylov() const = 0;
249 
255  mutable Solver* prec;
256 };
257 
258 
259 /** EliminationSolver using CG and HypreBoomerAMG */
261 {
262 public:
264  Array<int>& constraint_rowstarts,
265  int dimension_=0, bool reorder_=false) :
266  EliminationSolver(A, B, constraint_rowstarts),
267  dimension(dimension_), reorder(reorder_)
268  { }
269 
270 protected:
271  virtual Solver* BuildPreconditioner() const override
272  {
274  h_prec->SetPrintLevel(0);
275  if (dimension > 0) { h_prec->SetSystemsOptions(dimension, reorder); }
276  return h_prec;
277  }
278 
279  virtual IterativeSolver* BuildKrylov() const override
280  { return new CGSolver(GetComm()); }
281 
282 private:
283  int dimension;
284  bool reorder;
285 };
286 
287 /** EliminationSolver using GMRES and HypreBoomerAMG */
289 {
290 public:
292  Array<int>& constraint_rowstarts,
293  int dimension_=0, bool reorder_=false) :
294  EliminationSolver(A, B, constraint_rowstarts),
295  dimension(dimension_), reorder(reorder_)
296  { }
297 
298 protected:
299  virtual Solver* BuildPreconditioner() const override
300  {
302  h_prec->SetPrintLevel(0);
303  if (dimension > 0) { h_prec->SetSystemsOptions(dimension, reorder); }
304  return h_prec;
305  }
306 
307  virtual IterativeSolver* BuildKrylov() const override
308  { return new GMRESSolver(GetComm()); }
309 
310 private:
311  int dimension;
312  bool reorder;
313 };
314 
315 /** @brief Solve constrained system with penalty method; see ConstrainedSolver.
316 
317  Only approximates the solution, better approximation with higher penalty,
318  but with higher penalty the preconditioner is less effective. */
320 {
321 public:
323  double penalty_);
324 
326  double penalty_);
327 
329  Vector& penalty_);
330 
332 
333  void Mult(const Vector& x, Vector& y) const override;
334 
335  void SetOperator(const Operator& op) override
336  { MFEM_ABORT("Operator cannot be reset!"); }
337 
338  void SetPreconditioner(Solver& precond) override
339  { prec = &precond; }
340 
341 protected:
342  /// Initialize the matrix A + B*D*B^T for the constrained linear system
343  /// A - original matrix (N x N square matrix)
344  /// B - constraint matrix (M x N rectangular matrix)
345  /// D - diagonal matrix of penalty values (M x M square matrix)
346  void Initialize(HypreParMatrix& A, HypreParMatrix& B, HypreParMatrix& D);
347 
348  /// Build preconditioner for penalized system
349  virtual Solver* BuildPreconditioner() const = 0;
350  /// Select krylov solver for penalized system
351  virtual IterativeSolver* BuildKrylov() const = 0;
352 
357  mutable Solver * prec;
358 };
359 
360 
361 /** Uses CG and a HypreBoomerAMG preconditioner for the penalized system. */
363 {
364 public:
366  int dimension=0, bool reorder=false) :
367  PenaltyConstrainedSolver(A, B, penalty_),
368  dimension_(dimension), reorder_(reorder)
369  { }
370 
372  int dimension=0, bool reorder=false) :
373  PenaltyConstrainedSolver(A, B, penalty_),
374  dimension_(dimension), reorder_(reorder)
375  { }
376 
378  int dimension=0, bool reorder=false) :
379  PenaltyConstrainedSolver(A, B, penalty_),
380  dimension_(dimension), reorder_(reorder)
381  { }
382 
383 protected:
384  virtual Solver* BuildPreconditioner() const override
385  {
387  h_prec->SetPrintLevel(0);
388  if (dimension_ > 0) { h_prec->SetSystemsOptions(dimension_, reorder_); }
389  return h_prec;
390  }
391 
392  virtual IterativeSolver* BuildKrylov() const override
393  { return new CGSolver(GetComm()); }
394 
395 private:
396  int dimension_;
397  bool reorder_;
398 };
399 
400 /** Uses GMRES and a HypreBoomerAMG preconditioner for the penalized system. */
402 {
403 public:
405  int dimension=0, bool reorder=false) :
406  PenaltyConstrainedSolver(A, B, penalty_),
407  dimension_(dimension), reorder_(reorder)
408  { }
409 
411  int dimension=0, bool reorder=false) :
412  PenaltyConstrainedSolver(A, B, penalty_),
413  dimension_(dimension), reorder_(reorder)
414  { }
415 
417  int dimension=0, bool reorder=false) :
418  PenaltyConstrainedSolver(A, B, penalty_),
419  dimension_(dimension), reorder_(reorder)
420  { }
421 
422 protected:
423  virtual Solver* BuildPreconditioner() const override
424  {
426  h_prec->SetPrintLevel(0);
427  if (dimension_ > 0) { h_prec->SetSystemsOptions(dimension_, reorder_); }
428  return h_prec;
429  }
430 
431  virtual IterativeSolver* BuildKrylov() const override
432  { return new GMRESSolver(GetComm()); }
433 
434 private:
435  int dimension_;
436  bool reorder_;
437 };
438 
439 #endif
440 
441 /** @brief Solve constrained system by solving original mixed system;
442  see ConstrainedSolver.
443 
444  Solves the saddle-point problem with a block-diagonal preconditioner, with
445  user-provided preconditioner in the top-left block and (by default) an
446  identity matrix in the bottom-right.
447 
448  This is the most general ConstrainedSolver, needing only Operator objects
449  to function. But in general it is not very efficient or scalable. */
451 {
452 public:
453  /// Setup constrained system, with primal_pc a user-provided preconditioner
454  /// for the top-left block.
455 #ifdef MFEM_USE_MPI
456  SchurConstrainedSolver(MPI_Comm comm, Operator& A_, Operator& B_,
457  Solver& primal_pc_);
458 #endif
459  SchurConstrainedSolver(Operator& A_, Operator& B_, Solver& primal_pc_);
460  virtual ~SchurConstrainedSolver();
461 
462  virtual void LagrangeSystemMult(const Vector& x, Vector& y) const override;
463 
464 protected:
465 #ifdef MFEM_USE_MPI
466  SchurConstrainedSolver(MPI_Comm comm, Operator& A_, Operator& B_);
467 #endif
469 
473  Solver * primal_pc; // NOT owned
475  Solver * dual_pc; // owned
476 
477 private:
478  void Initialize();
479 };
480 
481 
482 #ifdef MFEM_USE_MPI
483 /** @brief Basic saddle-point solver with assembled blocks (ie, the
484  operators are assembled HypreParMatrix objects.)
485 
486  This uses a block-diagonal preconditioner that approximates
487  \f$ [ A^{-1} 0; 0 (B A^{-1} B^T)^{-1} ] \f$.
488 
489  In the top-left block, we approximate \f$ A^{-1} \f$ with HypreBoomerAMG.
490  In the bottom-right, we approximate \f$ A^{-1} \f$ with the inverse of the
491  diagonal of \f$ A \f$, assemble \f$ B diag(A)^{-1} B^T \f$, and use
492  HypreBoomerAMG on that assembled matrix. */
494 {
495 public:
496  SchurConstrainedHypreSolver(MPI_Comm comm, HypreParMatrix& hA_,
497  HypreParMatrix& hB_, Solver * prec = nullptr,
498  int dimension=0, bool reorder=false);
500 
501 private:
502  HypreParMatrix& hA;
503  HypreParMatrix& hB;
504  HypreParMatrix * schur_mat;
505 };
506 #endif
507 
508 /** @brief Build a matrix constraining normal components to zero.
509 
510  Given a vector space fespace, and the array constrained_att that
511  includes the boundary *attributes* that are constrained to have normal
512  component zero, this returns a SparseMatrix representing the
513  constraints that need to be imposed.
514 
515  Each row of the returned matrix corresponds to a node that is
516  constrained. The rows are arranged in (contiguous) blocks corresponding
517  to a physical constraint; in 3D, a one-row constraint means the node
518  is free to move along a plane, a two-row constraint means it is free
519  to move along a line (e.g. the intersection of two normal-constrained
520  planes), and a three-row constraint is fully constrained (equivalent
521  to MFEM's usual essential boundary conditions).
522 
523  The constraint_rowstarts array is filled in to describe the structure of
524  these constraints, so that (block) constraint k is encoded in rows
525  constraint_rowstarts[k] to constraint_rowstarts[k + 1] - 1, inclusive,
526  of the returned matrix.
527 
528  Constraints are imposed on "true" degrees of freedom, which are different
529  in serial and parallel, so we need different numbering systems for the
530  serial and parallel versions of this function.
531 
532  When two attributes intersect, this version will combine constraints,
533  so in 2D the point at the intersection is fully constrained (ie,
534  fixed in both directions). This is the wrong thing to do if the
535  two boundaries are (close to) parallel at that point.
536 
537  @param[in] fespace A vector finite element space
538  @param[in] constrained_att Boundary attributes to constrain
539  @param[out] constraint_rowstarts The rowstarts for separately
540  eliminated constraints, possible
541  input to EliminationCGSolver
542  @param[in] parallel Indicate that fespace is actually a
543  ParFiniteElementSpace and the numbering
544  in the returned matrix should be based
545  on truedofs.
546 
547  @return a constraint matrix
548 */
550  Array<int>& constrained_att,
551  Array<int>& constraint_rowstarts,
552  bool parallel=false);
553 
554 #ifdef MFEM_USE_MPI
555 /// Parallel wrapper for BuildNormalConstraints
557  Array<int>& constrained_att,
558  Array<int>& constraint_rowstarts);
559 #endif
560 
561 }
562 
563 #endif
HypreParMatrix & hA
void Eliminate(const Vector &in, Vector &out) const
Definition: constraints.cpp:51
Conjugate gradient method.
Definition: solvers.hpp:465
void Mult(const Vector &x, Vector &y) const override
Solve for given .
MPI_Comm GetComm() const
Return the associated MPI communicator, or MPI_COMM_NULL if no communicator is set.
Definition: solvers.hpp:264
virtual void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
Definition: constraints.hpp:63
SparseMatrix * ParBuildNormalConstraints(ParFiniteElementSpace &fespace, Array< int > &constrained_att, Array< int > &constraint_rowstarts)
Parallel wrapper for BuildNormalConstraints.
Solve constrained system by solving original mixed system; see ConstrainedSolver. ...
EliminationSolver(HypreParMatrix &A, SparseMatrix &B, Array< int > &primary_dofs, Array< int > &secondary_dofs)
Constructor, with explicit splitting into primary/secondary dofs.
void Mult(const Vector &x, Vector &y) const override
Operator application: y=A(x).
Definition: constraints.cpp:96
SparseMatrix * BuildNormalConstraints(FiniteElementSpace &fespace, Array< int > &constrained_att, Array< int > &constraint_rowstarts, bool parallel)
Build a matrix constraining normal components to zero.
SchurConstrainedHypreSolver(MPI_Comm comm, HypreParMatrix &hA_, HypreParMatrix &hB_, Solver *prec=nullptr, int dimension=0, bool reorder=false)
void SetPreconditioner(Solver &precond) override
This should be called before SetOperator.
EliminationProjection(const Operator &A, Array< Eliminator * > &eliminators)
Definition: constraints.cpp:87
An abstract class to solve the constrained system subject to the constraint .
Definition: constraints.hpp:53
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
virtual Solver * BuildPreconditioner() const override
Build preconditioner for penalized system.
virtual Solver * BuildPreconditioner() const override
Build preconditioner for penalized system.
Basic saddle-point solver with assembled blocks (ie, the operators are assembled HypreParMatrix objec...
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual IterativeSolver * BuildKrylov() const override
Select krylov solver for eliminated system.
HypreParMatrix * h_explicit_operator
Eliminator(const SparseMatrix &B, const Array< int > &lagrange_dofs, const Array< int > &primary_tdofs, const Array< int > &secondary_tdofs)
Definition: constraints.cpp:22
void MultTranspose(const Vector &x, Vector &y) const override
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
virtual Solver * BuildPreconditioner() const override
Build preconditioner for eliminated system.
PenaltyGMRESSolver(HypreParMatrix &A, SparseMatrix &B, double penalty_, int dimension=0, bool reorder=false)
virtual Solver * BuildPreconditioner() const =0
Build preconditioner for penalized system.
Array< Eliminator * > eliminators
ConstrainedSolver(MPI_Comm comm, Operator &A_, Operator &B_)
virtual void Mult(const Vector &f, Vector &x) const override
Solve for given .
void LagrangeSecondary(const Vector &in, Vector &out) const
Maps Lagrange multipliers to secondary dofs, applies .
Definition: constraints.cpp:66
void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1579
Solve constrained system by eliminating the constraint; see ConstrainedSolver.
void BuildGTilde(const Vector &g, Vector &gtilde) const
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
A class to handle Block diagonal preconditioners in a matrix-free implementation. ...
Data type sparse matrix.
Definition: sparsemat.hpp:50
PenaltyGMRESSolver(HypreParMatrix &A, HypreParMatrix &B, Vector &penalty_, int dimension=0, bool reorder=false)
Perform elimination of a single constraint.
virtual IterativeSolver * BuildKrylov() const =0
Select krylov solver for penalized system.
PenaltyPCGSolver(HypreParMatrix &A, SparseMatrix &B, double penalty_, int dimension=0, bool reorder=false)
void SetPrintLevel(int print_level)
Definition: hypre.hpp:1659
const Array< int > & PrimaryDofs() const
void RecoverMultiplier(const Vector &primalrhs, const Vector &primalvars, Vector &lm) const
EliminationCGSolver(HypreParMatrix &A, SparseMatrix &B, Array< int > &constraint_rowstarts, int dimension_=0, bool reorder_=false)
virtual IterativeSolver * BuildKrylov() const override
Select krylov solver for eliminated system.
virtual void SetConstraintRHS(const Vector &r)
Set the right-hand side r for the constraint B x = r.
EliminationGMRESSolver(HypreParMatrix &A, SparseMatrix &B, Array< int > &constraint_rowstarts, int dimension_=0, bool reorder_=false)
void ExplicitAssembly(DenseMatrix &mat) const
Return explicitly assembled in mat.
Definition: constraints.cpp:79
Solve constrained system with penalty method; see ConstrainedSolver.
void EliminateTranspose(const Vector &in, Vector &out) const
Transpose of Eliminate(), applies .
Definition: constraints.cpp:58
void BuildExplicitOperator()
Internal utility routine; assembles eliminated matrix explicitly.
PenaltyConstrainedSolver(HypreParMatrix &A, SparseMatrix &B, double penalty_)
Abstract base class for iterative solver.
Definition: solvers.hpp:66
void GetMultiplierSolution(Vector &lambda) const
Return the Lagrange multiplier solution in lambda.
Definition: constraints.hpp:76
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
PenaltyPCGSolver(HypreParMatrix &A, HypreParMatrix &B, Vector &penalty_, int dimension=0, bool reorder=false)
virtual void LagrangeSystemMult(const Vector &f_and_r, Vector &x_and_lambda) const
Solve for (x, lambda) given (f, r)
GMRES method.
Definition: solvers.hpp:497
SchurConstrainedSolver(MPI_Comm comm, Operator &A_, Operator &B_, Solver &primal_pc_)
virtual Solver * BuildPreconditioner() const =0
Build preconditioner for eliminated system.
The transpose of a given operator. Switches the roles of the methods Mult() and MultTranspose().
Definition: operator.hpp:722
EliminationProjection * projector
void Mult(const Vector &x, Vector &y) const override
Solve for given .
virtual Solver * BuildPreconditioner() const override
Build preconditioner for eliminated system.
IterativeSolver * krylov
PenaltyGMRESSolver(HypreParMatrix &A, HypreParMatrix &B, double penalty_, int dimension=0, bool reorder=false)
PenaltyPCGSolver(HypreParMatrix &A, HypreParMatrix &B, double penalty_, int dimension=0, bool reorder=false)
constexpr int dimension
This example only works in 3D. Kernels for 2D are not implemented.
Definition: hooke.cpp:45
virtual void LagrangeSystemMult(const Vector &x, Vector &y) const override
Solve for (x, lambda) given (f, r)
virtual IterativeSolver * BuildKrylov() const =0
Select krylov solver for eliminated system.
void SetPreconditioner(Solver &precond) override
This should be called before SetOperator.
Vector data type.
Definition: vector.hpp:60
const Array< int > & LagrangeDofs() const
void SetSystemsOptions(int dim, bool order_bynodes=false)
Definition: hypre.cpp:4962
BlockDiagonalPreconditioner * block_pc
Base class for solvers.
Definition: operator.hpp:655
SparseMatrix * AssembleExact() const
Assemble this projector as a (processor-local) SparseMatrix.
virtual IterativeSolver * BuildKrylov() const override
Select krylov solver for penalized system.
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
A class to handle Block systems in a matrix-free implementation.
TransposeOperator * tr_B
void LagrangeSecondaryTranspose(const Vector &in, Vector &out) const
Transpose of LagrangeSecondary()
Definition: constraints.cpp:72
const Array< int > & SecondaryDofs() const
virtual IterativeSolver * BuildKrylov() const override
Select krylov solver for penalized system.