MFEM  v3.3.2 Finite element discretization library
operator.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
10 // Software Foundation) version 2.1 dated February 1999.
11
12 #ifndef MFEM_OPERATOR
13 #define MFEM_OPERATOR
14
15 #include "vector.hpp"
16
17 namespace mfem
18 {
19
20 /// Abstract operator
21 class Operator
22 {
23 protected:
24  int height; ///< Dimension of the output / number of rows in the matrix.
25  int width; ///< Dimension of the input / number of columns in the matrix.
26
27 public:
28  /// Construct a square Operator with given size s (default 0).
29  explicit Operator(int s = 0) { height = width = s; }
30
31  /** @brief Construct an Operator with the given height (output size) and
32  width (input size). */
33  Operator(int h, int w) { height = h; width = w; }
34
35  /// Get the height (size of output) of the Operator. Synonym with NumRows().
36  inline int Height() const { return height; }
37  /** @brief Get the number of rows (size of output) of the Operator. Synonym
38  with Height(). */
39  inline int NumRows() const { return height; }
40
41  /// Get the width (size of input) of the Operator. Synonym with NumCols().
42  inline int Width() const { return width; }
43  /** @brief Get the number of columns (size of input) of the Operator. Synonym
44  with Width(). */
45  inline int NumCols() const { return width; }
46
47  /// Operator application: `y=A(x)`.
48  virtual void Mult(const Vector &x, Vector &y) const = 0;
49
50  /** @brief Action of the transpose operator: `y=A^t(x)`. The default behavior
51  in class Operator is to generate an error. */
52  virtual void MultTranspose(const Vector &x, Vector &y) const
53  { mfem_error("Operator::MultTranspose() is not overloaded!"); }
54
55  /** @brief Evaluate the gradient operator at the point @a x. The default
56  behavior in class Operator is to generate an error. */
57  virtual Operator &GetGradient(const Vector &x) const
58  {
60  return const_cast<Operator &>(*this);
61  }
62
63  /** @brief Prolongation operator from linear algebra (linear system) vectors,
64  to input vectors for the operator. `NULL` means identity. */
65  virtual const Operator *GetProlongation() const { return NULL; }
66  /** @brief Restriction operator from input vectors for the operator to linear
67  algebra (linear system) vectors. `NULL` means identity. */
68  virtual const Operator *GetRestriction() const { return NULL; }
69
70  /** @brief Form a constrained linear system using a matrix-free approach.
71
72  Assuming square operator, form the operator linear system `A(X)=B`,
73  corresponding to it and the right-hand side @a b, by applying any
74  necessary transformations such as: parallel assembly, conforming
75  constraints for non-conforming AMR and eliminating boundary conditions.
76  @note Static condensation and hybridization are not supported for general
77  operators (cf. the analogous methods BilinearForm::FormLinearSystem() and
78  ParBilinearForm::FormLinearSystem()).
79
80  The constraints are specified through the prolongation P from
81  GetProlongation(), and restriction R from GetRestriction() methods, which
82  are e.g. available through the (parallel) finite element space of any
83  (parallel) bilinear form operator. We assume that the operator is square,
84  using the same input and output space, so we have: `A(X)=[P^t (*this)
85  P](X)`, `B=P^t(b)`, and `x=P(X)`.
86
87  The vector @a x must contain the essential boundary condition values.
88  These are eliminated through the ConstrainedOperator class and the vector
89  @a X is initialized by setting its essential entries to the boundary
90  conditions and all other entries to zero (@a copy_interior == 0) or
91  copied from @a x (@a copy_interior != 0).
92
93  After solving the system `A(X)=B`, the (finite element) solution @a x can
94  be recovered by calling Operator::RecoverFEMSolution() with the same
95  vectors @a X, @a b, and @a x.
96
97  @note The caller is responsible for destroying the output operator @a A!
98  @note If there are no transformations, @a X simply reuses the data of @a
99  x. */
100  void FormLinearSystem(const Array<int> &ess_tdof_list,
101  Vector &x, Vector &b,
102  Operator* &A, Vector &X, Vector &B,
103  int copy_interior = 0);
104
105  /** @brief Reconstruct a solution vector @a x (e.g. a GridFunction) from the
106  solution @a X of a constrained linear system obtained from
107  Operator::FormLinearSystem().
108
109  Call this method after solving a linear system constructed using
110  Operator::FormLinearSystem() to recover the solution as an input vector,
111  @a x, for this Operator (presumably a finite element grid function). This
112  method has identical signature to the analogous method for bilinear
113  forms, though currently @a b is not used in the implementation. */
114  virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x);
115
116  /// Prints operator with input size n and output size m in Matlab format.
117  void PrintMatlab(std::ostream & out, int n = 0, int m = 0) const;
118
119  /// Virtual destructor.
120  virtual ~Operator() { }
121
122  /// Enumeration defining IDs for some classes derived from Operator.
123  enum Type
124  {
133  };
134 };
135
136
137 /// Base abstract class for time dependent operators.
138 /** Operator of the form: (x,t) -> f(x,t), where k = f(x,t) generally solves the
139  algebraic equation F(x,k,t) = G(x,t). The functions F and G represent the
140  _implicit_ and _explicit_ parts of the operator, respectively. For explicit
141  operators, F(x,k,t) = k, so f(x,t) = G(x,t).*/
143 {
144 public:
145  enum Type
146  {
147  EXPLICIT, ///< This type assumes F(x,k,t) = k, i.e. k = f(x,t) = G(x,t).
148  IMPLICIT, ///< This is the most general type, no assumptions on F and G.
149  HOMOGENEOUS ///< This type assumes that G(x,t) = 0.
150  };
151
152 protected:
153  double t; ///< Current time.
154  Type type; ///< Describes the form of the TimeDependentOperator.
155
156 public:
157  /** @brief Construct a "square" TimeDependentOperator y = f(x,t), where x and
158  y have the same dimension @a n. */
159  explicit TimeDependentOperator(int n = 0, double t_ = 0.0,
160  Type type_ = EXPLICIT)
161  : Operator(n) { t = t_; type = type_; }
162
163  /** @brief Construct a TimeDependentOperator y = f(x,t), where x and y have
164  dimensions @a w and @a h, respectively. */
165  TimeDependentOperator(int h, int w, double t_ = 0.0, Type type_ = EXPLICIT)
166  : Operator(h, w) { t = t_; type = type_; }
167
168  /// Read the currently set time.
169  virtual double GetTime() const { return t; }
170
171  /// Set the current time.
172  virtual void SetTime(const double _t) { t = _t; }
173
174  /// True if #type is #EXPLICIT.
175  bool isExplicit() const { return (type == EXPLICIT); }
176  /// True if #type is #IMPLICIT or #HOMOGENEOUS.
177  bool isImplicit() const { return !isExplicit(); }
178  /// True if #type is #HOMOGENEOUS.
179  bool isHomogeneous() const { return (type == HOMOGENEOUS); }
180
181  /** @brief Perform the action of the explicit part of the operator, G:
182  @a y = G(@a x, t) where t is the current time.
183
184  Presently, this method is used by some PETSc ODE solvers, for more
185  details, see the PETSc Manual. */
186  virtual void ExplicitMult(const Vector &x, Vector &y) const
187  {
188  mfem_error("TimeDependentOperator::ExplicitMult() is not overridden!");
189  }
190
191  /** @brief Perform the action of the implicit part of the operator, F:
192  @a y = F(@a x, @a k, t) where t is the current time.
193
194  Presently, this method is used by some PETSc ODE solvers, for more
195  details, see the PETSc Manual.*/
196  virtual void ImplicitMult(const Vector &x, const Vector &k, Vector &y) const
197  {
198  mfem_error("TimeDependentOperator::ImplicitMult() is not overridden!");
199  }
200
201  /** @brief Perform the action of the operator: @a y = k = f(@a x, t), where
202  k solves the algebraic equation F(@a x, k, t) = G(@a x, t) and t is the
203  current time. */
204  virtual void Mult(const Vector &x, Vector &y) const
205  {
206  mfem_error("TimeDependentOperator::Mult() is not overridden!");
207  }
208
209  /** @brief Solve the equation: @a k = f(@a x + @a dt @a k, t), for the
210  unknown @a k at the current time t.
211
212  For general F and G, the equation for @a k becomes:
213  F(@a x + @a dt @a k, @a k, t) = G(@a x + @a dt @a k, t).
214
215  The input vector @a x corresponds to time index (or cycle) n, while the
216  currently set time, #t, and the result vector @a k correspond to time
217  index n+1. The time step @a dt corresponds to the time interval between
218  cycles n and n+1.
219
220  This method allows for the abstract implementation of some time
221  integration methods, including diagonal implicit Runge-Kutta (DIRK)
222  methods and the backward Euler method in particular.
223
224  If not re-implemented, this method simply generates an error. */
225  virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k)
226  {
227  mfem_error("TimeDependentOperator::ImplicitSolve() is not overridden!");
228  }
229
230  /** @brief Return an Operator representing (dF/dk @a shift + dF/dx) at the
231  given @a x, @a k, and the currently set time.
232
233  Presently, this method is used by some PETSc ODE solvers, for more
234  details, see the PETSc Manual. */
235  virtual Operator& GetImplicitGradient(const Vector &x, const Vector &k,
236  double shift) const
237  {
239  "not overridden!");
240  return const_cast<Operator &>(dynamic_cast<const Operator &>(*this));
241  }
242
243  /** @brief Return an Operator representing dG/dx at the given point @a x and
244  the currently set time.
245
246  Presently, this method is used by some PETSc ODE solvers, for more
247  details, see the PETSc Manual. */
248  virtual Operator& GetExplicitGradient(const Vector &x) const
249  {
251  "not overridden!");
252  return const_cast<Operator &>(dynamic_cast<const Operator &>(*this));
253  }
254
256 };
257
258 /// Base class for solvers
259 class Solver : public Operator
260 {
261 public:
262  /// If true, use the second argument of Mult() as an initial guess.
264
265  /** @brief Initialize a square Solver with size @a s.
266
267  @warning Use a Boolean expression for the second parameter (not an int)
268  to distinguish this call from the general rectangular constructor. */
269  explicit Solver(int s = 0, bool iter_mode = false)
270  : Operator(s) { iterative_mode = iter_mode; }
271
272  /// Initialize a Solver with height @a h and width @a w.
273  Solver(int h, int w, bool iter_mode = false)
274  : Operator(h, w) { iterative_mode = iter_mode; }
275
276  /// Set/update the solver for the given operator.
277  virtual void SetOperator(const Operator &op) = 0;
278 };
279
280
281 /// Identity Operator I: x -> x.
283 {
284 public:
285  /// Create an identity operator of size @a n.
286  explicit IdentityOperator(int n) : Operator(n) { }
287
288  /// Operator application
289  virtual void Mult(const Vector &x, Vector &y) const { y = x; }
290 };
291
292
293 /** @brief The transpose of a given operator. Switches the roles of the methods
294  Mult() and MultTranspose(). */
296 {
297 private:
298  const Operator &A;
299
300 public:
301  /// Construct the transpose of a given operator @a *a.
303  : Operator(a->Width(), a->Height()), A(*a) { }
304
305  /// Construct the transpose of a given operator @a a.
307  : Operator(a.Width(), a.Height()), A(a) { }
308
309  /// Operator application. Apply the transpose of the original Operator.
310  virtual void Mult(const Vector &x, Vector &y) const
311  { A.MultTranspose(x, y); }
312
313  /// Application of the transpose. Apply the original Operator.
314  virtual void MultTranspose(const Vector &x, Vector &y) const
315  { A.Mult(x, y); }
316 };
317
318
319 /// The operator x -> R*A*P*x.
320 class RAPOperator : public Operator
321 {
322 private:
323  const Operator & Rt;
324  const Operator & A;
325  const Operator & P;
326  mutable Vector Px;
327  mutable Vector APx;
328
329 public:
330  /// Construct the RAP operator given R^T, A and P.
331  RAPOperator(const Operator &Rt_, const Operator &A_, const Operator &P_)
332  : Operator(Rt_.Width(), P_.Width()), Rt(Rt_), A(A_), P(P_),
333  Px(P.Height()), APx(A.Height()) { }
334
335  /// Operator application.
336  virtual void Mult(const Vector & x, Vector & y) const
337  { P.Mult(x, Px); A.Mult(Px, APx); Rt.MultTranspose(APx, y); }
338
339  /// Application of the transpose.
340  virtual void MultTranspose(const Vector & x, Vector & y) const
341  { Rt.Mult(x, APx); A.MultTranspose(APx, Px); P.MultTranspose(Px, y); }
342 };
343
344
345 /// General triple product operator x -> A*B*C*x, with ownership of the factors.
347 {
348  Operator *A;
349  Operator *B;
350  Operator *C;
351  bool ownA, ownB, ownC;
352  mutable Vector t1, t2;
353
354 public:
356  bool ownA, bool ownB, bool ownC)
357  : Operator(A->Height(), C->Width())
358  , A(A), B(B), C(C)
359  , ownA(ownA), ownB(ownB), ownC(ownC)
360  , t1(C->Height()), t2(B->Height())
361  {}
362
363  virtual void Mult(const Vector &x, Vector &y) const
364  { C->Mult(x, t1); B->Mult(t1, t2); A->Mult(t2, y); }
365
366  virtual void MultTranspose(const Vector &x, Vector &y) const
367  { A->MultTranspose(x, t2); B->MultTranspose(t2, t1); C->MultTranspose(t1, y); }
368
370  {
371  if (ownA) { delete A; }
372  if (ownB) { delete B; }
373  if (ownC) { delete C; }
374  }
375 };
376
377
378 /** @brief Square Operator for imposing essential boundary conditions using only
379  the action, Mult(), of a given unconstrained Operator.
380
381  Square operator constrained by fixing certain entries in the solution to
382  given "essential boundary condition" values. This class is used by the
383  general, matrix-free system formulation of Operator::FormLinearSystem. */
385 {
386 protected:
387  Array<int> constraint_list; ///< List of constrained indices/dofs.
388  Operator *A; ///< The unconstrained Operator.
389  bool own_A; ///< Ownership flag for A.
390  mutable Vector z, w; ///< Auxiliary vectors.
391
392 public:
393  /** @brief Constructor from a general Operator and a list of essential
394  indices/dofs.
395
396  Specify the unconstrained operator @a *A and a @a list of indices to
397  constrain, i.e. each entry @a list[i] represents an essential-dof. If the
398  ownership flag @a own_A is true, the operator @a *A will be destroyed
399  when this object is destroyed. */
400  ConstrainedOperator(Operator *A, const Array<int> &list, bool own_A = false);
401
402  /** @brief Eliminate "essential boundary condition" values specified in @a x
403  from the given right-hand side @a b.
404
405  Performs the following steps:
406
407  z = A((0,x_b)); b_i -= z_i; b_b = x_b;
408
409  where the "_b" subscripts denote the essential (boundary) indices/dofs of
410  the vectors, and "_i" -- the rest of the entries. */
411  void EliminateRHS(const Vector &x, Vector &b) const;
412
413  /** @brief Constrained operator action.
414
415  Performs the following steps:
416
417  z = A((x_i,0)); y_i = z_i; y_b = x_b;
418
419  where the "_b" subscripts denote the essential (boundary) indices/dofs of
420  the vectors, and "_i" -- the rest of the entries. */
421  virtual void Mult(const Vector &x, Vector &y) const;
422
423  /// Destructor: destroys the unconstrained Operator @a A if @a own_A is true.
424  virtual ~ConstrainedOperator() { if (own_A) { delete A; } }
425 };
426
427 }
428
429 #endif
virtual double GetTime() const
Definition: operator.hpp:169
bool isImplicit() const
True if type is IMPLICIT or HOMOGENEOUS.
Definition: operator.hpp:177
bool isHomogeneous() const
True if type is HOMOGENEOUS.
Definition: operator.hpp:179
Solver(int s=0, bool iter_mode=false)
Initialize a square Solver with size s.
Definition: operator.hpp:269
Array< int > constraint_list
List of constrained indices/dofs.
Definition: operator.hpp:387
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
Definition: operator.hpp:336
Type type
Describes the form of the TimeDependentOperator.
Definition: operator.hpp:154
Base abstract class for time dependent operators.
Definition: operator.hpp:142
virtual Operator & GetGradient(const Vector &x) const
Evaluate the gradient operator at the point x. The default behavior in class Operator is to generate ...
Definition: operator.hpp:57
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:42
virtual void Mult(const Vector &x, Vector &y) const
Constrained operator action.
Definition: operator.cpp:126
virtual void SetTime(const double _t)
Set the current time.
Definition: operator.hpp:172
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
TripleProductOperator(Operator *A, Operator *B, Operator *C, bool ownA, bool ownB, bool ownC)
Definition: operator.hpp:355
virtual void MultTranspose(const Vector &x, Vector &y) const
Application of the transpose. Apply the original Operator.
Definition: operator.hpp:314
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:263
virtual Operator & GetExplicitGradient(const Vector &x) const
Return an Operator representing dG/dx at the given point x and the currently set time.
Definition: operator.hpp:248
bool own_A
Ownership flag for A.
Definition: operator.hpp:389
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition: operator.hpp:52
Solver(int h, int w, bool iter_mode=false)
Initialize a Solver with height h and width w.
Definition: operator.hpp:273
double t
Current time.
Definition: operator.hpp:153
virtual void Mult(const Vector &x, Vector &y) const
Operator application. Apply the transpose of the original Operator.
Definition: operator.hpp:310
virtual Operator & GetImplicitGradient(const Vector &x, const Vector &k, double shift) const
Return an Operator representing (dF/dk shift + dF/dx) at the given x, k, and the currently set time...
Definition: operator.hpp:235
virtual const Operator * GetProlongation() const
Prolongation operator from linear algebra (linear system) vectors, to input vectors for the operator...
Definition: operator.hpp:65
TransposeOperator(const Operator &a)
Construct the transpose of a given operator a.
Definition: operator.hpp:306
virtual void MultTranspose(const Vector &x, Vector &y) const
Application of the transpose.
Definition: operator.hpp:340
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:36
virtual ~ConstrainedOperator()
Destructor: destroys the unconstrained Operator A if own_A is true.
Definition: operator.hpp:424
This type assumes F(x,k,t) = k, i.e. k = f(x,t) = G(x,t).
Definition: operator.hpp:147
TimeDependentOperator(int h, int w, double t_=0.0, Type type_=EXPLICIT)
Construct a TimeDependentOperator y = f(x,t), where x and y have dimensions w and h...
Definition: operator.hpp:165
The operator x -&gt; R*A*P*x.
Definition: operator.hpp:320
bool isExplicit() const
True if type is EXPLICIT.
Definition: operator.hpp:175
virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k)
Solve the equation: k = f(x + dt k, t), for the unknown k at the current time t.
Definition: operator.hpp:225
RAPOperator(const Operator &Rt_, const Operator &A_, const Operator &P_)
Construct the RAP operator given R^T, A and P.
Definition: operator.hpp:331
virtual const Operator * GetRestriction() const
Restriction operator from input vectors for the operator to linear algebra (linear system) vectors...
Definition: operator.hpp:68
This is the most general type, no assumptions on F and G.
Definition: operator.hpp:148
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: operator.hpp:363
Operator(int s=0)
Construct a square Operator with given size s (default 0).
Definition: operator.hpp:29
IdentityOperator(int n)
Create an identity operator of size n.
Definition: operator.hpp:286
Type
Enumeration defining IDs for some classes derived from Operator.
Definition: operator.hpp:123
TimeDependentOperator(int n=0, double t_=0.0, Type type_=EXPLICIT)
Construct a &quot;square&quot; TimeDependentOperator y = f(x,t), where x and y have the same dimension n...
Definition: operator.hpp:159
Vector w
Auxiliary vectors.
Definition: operator.hpp:390
virtual void Mult(const Vector &x, Vector &y) const
Perform the action of the operator: y = k = f(x, t), where k solves the algebraic equation F(x...
Definition: operator.hpp:204
virtual void ExplicitMult(const Vector &x, Vector &y) const
Perform the action of the explicit part of the operator, G: y = G(x, t) where t is the current time...
Definition: operator.hpp:186
int NumRows() const
Get the number of rows (size of output) of the Operator. Synonym with Height().
Definition: operator.hpp:39
void mfem_error(const char *msg)
Definition: error.cpp:107
int NumCols() const
Get the number of columns (size of input) of the Operator. Synonym with Width().
Definition: operator.hpp:45
void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, Operator *&A, Vector &X, Vector &B, int copy_interior=0)
Form a constrained linear system using a matrix-free approach.
Definition: operator.cpp:21
The transpose of a given operator. Switches the roles of the methods Mult() and MultTranspose().
Definition: operator.hpp:295
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Reconstruct a solution vector x (e.g. a GridFunction) from the solution X of a constrained linear sys...
Definition: operator.cpp:57
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:24
This type assumes that G(x,t) = 0.
Definition: operator.hpp:149
void EliminateRHS(const Vector &x, Vector &b) const
Eliminate &quot;essential boundary condition&quot; values specified in x from the given right-hand side b...
Definition: operator.cpp:107
void PrintMatlab(std::ostream &out, int n=0, int m=0) const
Prints operator with input size n and output size m in Matlab format.
Definition: operator.cpp:72
General triple product operator x -&gt; A*B*C*x, with ownership of the factors.
Definition: operator.hpp:346
ConstrainedOperator(Operator *A, const Array< int > &list, bool own_A=false)
Constructor from a general Operator and a list of essential indices/dofs.
Definition: operator.cpp:98
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
Definition: operator.hpp:289
Operator(int h, int w)
Construct an Operator with the given height (output size) and width (input size). ...
Definition: operator.hpp:33
virtual void SetOperator(const Operator &op)=0
Set/update the solver for the given operator.
Vector data type.
Definition: vector.hpp:41
Identity Operator I: x -&gt; x.
Definition: operator.hpp:282
TransposeOperator(const Operator *a)
Construct the transpose of a given operator *a.
Definition: operator.hpp:302
Base class for solvers.
Definition: operator.hpp:259
Operator * A
The unconstrained Operator.
Definition: operator.hpp:388
Square Operator for imposing essential boundary conditions using only the action, Mult()...
Definition: operator.hpp:384
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:64
Abstract operator.
Definition: operator.hpp:21
virtual ~Operator()
Virtual destructor.
Definition: operator.hpp:120
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition: operator.hpp:366
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:25
virtual void ImplicitMult(const Vector &x, const Vector &k, Vector &y) const
Perform the action of the implicit part of the operator, F: y = F(x, k, t) where t is the current tim...
Definition: operator.hpp:196