MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
constraints.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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
20namespace mfem
21{
22
23class FiniteElementSpace;
24#ifdef MFEM_USE_MPI
25class ParFiniteElementSpace;
26#endif
27
28/** @brief An abstract class to solve the constrained system $Ax = f$
29 subject to the constraint $B x = r$.
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 $A$, and the Mult() interface just
43 takes $f$ as an argument. You can set $r$ 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{
55public:
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 $x$ given $f$.
79
80 If you want to set $r$, call SetConstraintRHS() before this.
81
82 If you want to get $\lambda$, 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
96protected:
99
102 mutable Vector workb;
103 mutable Vector workx;
104
105private:
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 $B_s^{-1}$ maps the lagrange space into secondary dofs, while
118 $-B_s^{-1} B_p$ maps primary dofs to secondary dofs. */
120{
121public:
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 $-B_s^{-1} B_p$.
132 void Eliminate(const Vector& in, Vector& out) const;
133
134 /// Transpose of Eliminate(), applies $-B_p^T B_s^{-T}$
135 void EliminateTranspose(const Vector& in, Vector& out) const;
136
137 /// Maps Lagrange multipliers to secondary dofs, applies $B_s^{-1}$
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 $-B_s^{-1} B_p$ explicitly assembled in mat
144 void ExplicitAssembly(DenseMatrix& mat) const;
145
146private:
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{
170public:
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 $g$, return
183 $\tilde{g}$ */
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
190private:
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 $P^T A P + Z_P$, where P is
201 EliminationProjection and Z_P is the identity on the eliminated dofs. */
203{
204public:
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
241protected:
242 /// Internal utility routine; assembles eliminated matrix explicitly
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{
262public:
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
270protected:
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
282private:
283 int dimension;
284 bool reorder;
285};
286
287/** EliminationSolver using GMRES and HypreBoomerAMG */
289{
290public:
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
298protected:
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
310private:
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{
321public:
323 real_t penalty_);
324
326 real_t 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
341protected:
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{
364public:
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
383protected:
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
395private:
396 int dimension_;
397 bool reorder_;
398};
399
400/** Uses GMRES and a HypreBoomerAMG preconditioner for the penalized system. */
402{
403public:
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
422protected:
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
434private:
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{
452public:
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
464protected:
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
477private:
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 $[ A^{-1} 0; 0 (B A^{-1} B^T)^{-1} ]$.
488
489 In the top-left block, we approximate $A^{-1}$ with HypreBoomerAMG.
490 In the bottom-right, we approximate $A^{-1}$ with the inverse of the
491 diagonal of $A$, assemble $B diag(A)^{-1} B^T$, and use
492 HypreBoomerAMG on that assembled matrix. */
494{
495public:
497 HypreParMatrix& hB_, Solver * prec = nullptr,
498 int dimension=0, bool reorder=false);
500
501private:
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
A class to handle Block diagonal preconditioners in a matrix-free implementation.
A class to handle Block systems in a matrix-free implementation.
Definition solvers.hpp:513
An abstract class to solve the constrained system subject to the constraint .
virtual void LagrangeSystemMult(const Vector &f_and_r, Vector &x_and_lambda) const
Solve for (x, lambda) given (f, r)
void GetMultiplierSolution(Vector &lambda) const
Return the Lagrange multiplier solution in lambda.
virtual void SetConstraintRHS(const Vector &r)
Set the right-hand side r for the constraint B x = r.
ConstrainedSolver(MPI_Comm comm, Operator &A_, Operator &B_)
virtual void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
virtual void Mult(const Vector &f, Vector &x) const override
Solve for given .
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
virtual Solver * BuildPreconditioner() const override
Build preconditioner for eliminated system.
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.
EliminationGMRESSolver(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 Solver * BuildPreconditioner() const override
Build preconditioner for eliminated system.
SparseMatrix * AssembleExact() const
Assemble this projector as a (processor-local) SparseMatrix.
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 ...
void RecoverMultiplier(const Vector &primalrhs, const Vector &primalvars, Vector &lm) const
void Mult(const Vector &x, Vector &y) const override
Operator application: y=A(x).
EliminationProjection(const Operator &A, Array< Eliminator * > &eliminators)
void BuildGTilde(const Vector &g, Vector &gtilde) const
Solve constrained system by eliminating the constraint; see ConstrainedSolver.
void Mult(const Vector &x, Vector &y) const override
Solve for given .
void SetPreconditioner(Solver &precond) override
This should be called before SetOperator.
IterativeSolver * krylov
Array< Eliminator * > eliminators
void BuildExplicitOperator()
Internal utility routine; assembles eliminated matrix explicitly.
EliminationSolver(HypreParMatrix &A, SparseMatrix &B, Array< int > &primary_dofs, Array< int > &secondary_dofs)
Constructor, with explicit splitting into primary/secondary dofs.
EliminationProjection * projector
virtual Solver * BuildPreconditioner() const =0
Build preconditioner for eliminated system.
virtual IterativeSolver * BuildKrylov() const =0
Select krylov solver for eliminated system.
void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
HypreParMatrix * h_explicit_operator
Perform elimination of a single constraint.
void Eliminate(const Vector &in, Vector &out) const
void ExplicitAssembly(DenseMatrix &mat) const
Return explicitly assembled in mat.
void LagrangeSecondaryTranspose(const Vector &in, Vector &out) const
Transpose of LagrangeSecondary()
void LagrangeSecondary(const Vector &in, Vector &out) const
Maps Lagrange multipliers to secondary dofs, applies .
Eliminator(const SparseMatrix &B, const Array< int > &lagrange_dofs, const Array< int > &primary_tdofs, const Array< int > &secondary_tdofs)
const Array< int > & PrimaryDofs() const
const Array< int > & SecondaryDofs() const
void EliminateTranspose(const Vector &in, Vector &out) const
Transpose of Eliminate(), applies .
const Array< int > & LagrangeDofs() const
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
GMRES method.
Definition solvers.hpp:547
The BoomerAMG solver in hypre.
Definition hypre.hpp:1691
void SetSystemsOptions(int dim, bool order_bynodes=false)
Definition hypre.cpp:5111
void SetPrintLevel(int print_level)
Definition hypre.hpp:1771
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
Abstract base class for iterative solver.
Definition solvers.hpp:67
MPI_Comm GetComm() const
Return the associated MPI communicator, or MPI_COMM_NULL if no communicator is set.
Definition solvers.hpp:302
Abstract operator.
Definition operator.hpp:25
Abstract parallel finite element space.
Definition pfespace.hpp:29
Solve constrained system with penalty method; see ConstrainedSolver.
virtual IterativeSolver * BuildKrylov() const =0
Select krylov solver for penalized system.
void Mult(const Vector &x, Vector &y) const override
Solve for given .
void SetOperator(const Operator &op) override
Also calls SetOperator for the preconditioner if there is one.
virtual Solver * BuildPreconditioner() const =0
Build preconditioner for penalized system.
void SetPreconditioner(Solver &precond) override
This should be called before SetOperator.
PenaltyConstrainedSolver(HypreParMatrix &A, SparseMatrix &B, real_t penalty_)
PenaltyGMRESSolver(HypreParMatrix &A, HypreParMatrix &B, real_t penalty_, int dimension=0, bool reorder=false)
virtual IterativeSolver * BuildKrylov() const override
Select krylov solver for penalized system.
virtual Solver * BuildPreconditioner() const override
Build preconditioner for penalized system.
PenaltyGMRESSolver(HypreParMatrix &A, HypreParMatrix &B, Vector &penalty_, int dimension=0, bool reorder=false)
PenaltyGMRESSolver(HypreParMatrix &A, SparseMatrix &B, real_t penalty_, int dimension=0, bool reorder=false)
PenaltyPCGSolver(HypreParMatrix &A, HypreParMatrix &B, Vector &penalty_, int dimension=0, bool reorder=false)
PenaltyPCGSolver(HypreParMatrix &A, HypreParMatrix &B, real_t penalty_, int dimension=0, bool reorder=false)
PenaltyPCGSolver(HypreParMatrix &A, SparseMatrix &B, real_t penalty_, int dimension=0, bool reorder=false)
virtual IterativeSolver * BuildKrylov() const override
Select krylov solver 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...
SchurConstrainedHypreSolver(MPI_Comm comm, HypreParMatrix &hA_, HypreParMatrix &hB_, Solver *prec=nullptr, int dimension=0, bool reorder=false)
Solve constrained system by solving original mixed system; see ConstrainedSolver.
BlockDiagonalPreconditioner * block_pc
SchurConstrainedSolver(MPI_Comm comm, Operator &A_, Operator &B_, Solver &primal_pc_)
virtual void LagrangeSystemMult(const Vector &x, Vector &y) const override
Solve for (x, lambda) given (f, r)
Base class for solvers.
Definition operator.hpp:683
Data type sparse matrix.
Definition sparsemat.hpp:51
The transpose of a given operator. Switches the roles of the methods Mult() and MultTranspose().
Definition operator.hpp:750
Vector data type.
Definition vector.hpp:80
constexpr int dimension
This example only works in 3D. Kernels for 2D are not implemented.
Definition hooke.cpp:45
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
SparseMatrix * BuildNormalConstraints(FiniteElementSpace &fespace, Array< int > &constrained_att, Array< int > &constraint_rowstarts, bool parallel)
Build a matrix constraining normal components to zero.
SparseMatrix * ParBuildNormalConstraints(ParFiniteElementSpace &fespace, Array< int > &constrained_att, Array< int > &constraint_rowstarts)
Parallel wrapper for BuildNormalConstraints.
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30