MFEM  v4.5.2 Finite element discretization library
constraints.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
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
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
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
Definition: solvers.hpp:493
void Mult(const Vector &x, Vector &y) const override
Solve for given .
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
void RecoverMultiplier(const Vector &primalrhs, const Vector &primalvars, Vector &lm) const
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.
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 void LagrangeSystemMult(const Vector &f_and_r, Vector &x_and_lambda) const
Solve for (x, lambda) given (f, r)
virtual Solver * BuildPreconditioner() const override
Build preconditioner for penalized system.
void ExplicitAssembly(DenseMatrix &mat) const
Return explicitly assembled in mat.
Definition: constraints.cpp:79
Basic saddle-point solver with assembled blocks (ie, the operators are assembled HypreParMatrix objec...
Abstract parallel finite element space.
Definition: pfespace.hpp:28
const Array< int > & SecondaryDofs() const
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_)
SparseMatrix * AssembleExact() const
Assemble this projector as a (processor-local) SparseMatrix.
virtual void Mult(const Vector &f, Vector &x) const override
Solve for given .
void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1590
void EliminateTranspose(const Vector &in, Vector &out) const
Transpose of Eliminate(), applies .
Definition: constraints.cpp:58
Solve constrained system by eliminating the constraint; see ConstrainedSolver.
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:1670
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)
Solve constrained system with penalty method; see ConstrainedSolver.
void BuildExplicitOperator()
Internal utility routine; assembles eliminated matrix explicitly.
void LagrangeSecondary(const Vector &in, Vector &out) const
Maps Lagrange multipliers to secondary dofs, applies .
Definition: constraints.cpp:66
PenaltyConstrainedSolver(HypreParMatrix &A, SparseMatrix &B, double penalty_)
Abstract base class for iterative solver.
Definition: solvers.hpp:66
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.
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
EliminationProjection(const Operator &A, Array< Eliminator *> &eliminators)
Definition: constraints.cpp:87
PenaltyPCGSolver(HypreParMatrix &A, HypreParMatrix &B, Vector &penalty_, int dimension=0, bool reorder=false)
GMRES method.
Definition: solvers.hpp:525
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:749
EliminationProjection * projector
const Array< int > & LagrangeDofs() const
void Mult(const Vector &x, Vector &y) const override
Solve for given .
virtual Solver * BuildPreconditioner() const override
Build preconditioner for eliminated system.
IterativeSolver * krylov
void BuildGTilde(const Vector &g, Vector &gtilde) const
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)
void Eliminate(const Vector &in, Vector &out) const
Definition: constraints.cpp:51
virtual IterativeSolver * BuildKrylov() const =0
Select krylov solver for eliminated system.
void SetPreconditioner(Solver &precond) override
This should be called before SetOperator.
MPI_Comm GetComm() const
Return the associated MPI communicator, or MPI_COMM_NULL if no communicator is set.
Definition: solvers.hpp:292
Vector data type.
Definition: vector.hpp:60
const Array< int > & PrimaryDofs() const
void SetSystemsOptions(int dim, bool order_bynodes=false)
Definition: hypre.cpp:4963
BlockDiagonalPreconditioner * block_pc
Base class for solvers.
Definition: operator.hpp:682
virtual IterativeSolver * BuildKrylov() const override
Select krylov solver for penalized system.
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.
void GetMultiplierSolution(Vector &lambda) const
Return the Lagrange multiplier solution in lambda.
Definition: constraints.hpp:76
TransposeOperator * tr_B
void LagrangeSecondaryTranspose(const Vector &in, Vector &out) const
Transpose of LagrangeSecondary()
Definition: constraints.cpp:72
virtual IterativeSolver * BuildKrylov() const override
Select krylov solver for penalized system.