MFEM  v3.3 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  {
131  };
132 };
133
134
135 /// Base abstract class for time dependent operators.
136 /** Operator of the form: (x,t) -> f(x,t), where k = f(x,t) generally solves the
137  algebraic equation F(x,k,t) = G(x,t). The functions F and G represent the
138  _implicit_ and _explicit_ parts of the operator, respectively. For explicit
139  operators, F(x,k,t) = k, so f(x,t) = G(x,t).*/
141 {
142 public:
143  enum Type
144  {
145  EXPLICIT, ///< This type assumes F(x,k,t) = k, i.e. k = f(x,t) = G(x,t).
146  IMPLICIT, ///< This is the most general type, no assumptions on F and G.
147  HOMOGENEOUS ///< This type assumes that G(x,t) = 0.
148  };
149
150 protected:
151  double t; ///< Current time.
152  Type type; ///< Describes the form of the TimeDependentOperator.
153
154 public:
155  /** @brief Construct a "square" TimeDependentOperator y = f(x,t), where x and
156  y have the same dimension @a n. */
157  explicit TimeDependentOperator(int n = 0, double t_ = 0.0,
158  Type type_ = EXPLICIT)
159  : Operator(n) { t = t_; type = type_; }
160
161  /** @brief Construct a TimeDependentOperator y = f(x,t), where x and y have
162  dimensions @a w and @a h, respectively. */
163  TimeDependentOperator(int h, int w, double t_ = 0.0, Type type_ = EXPLICIT)
164  : Operator(h, w) { t = t_; type = type_; }
165
166  /// Read the currently set time.
167  virtual double GetTime() const { return t; }
168
169  /// Set the current time.
170  virtual void SetTime(const double _t) { t = _t; }
171
172  /// True if #type is #EXPLICIT.
173  bool isExplicit() const { return (type == EXPLICIT); }
174  /// True if #type is #IMPLICIT or #HOMOGENEOUS.
175  bool isImplicit() const { return !isExplicit(); }
176  /// True if #type is #HOMOGENEOUS.
177  bool isHomogeneous() const { return (type == HOMOGENEOUS); }
178
179  /** @brief Perform the action of the explicit part of the operator, G:
180  @a y = G(@a x, t) where t is the current time.
181
182  Presently, this method is used by some PETSc ODE solvers, for more
183  details, see the PETSc Manual. */
184  virtual void ExplicitMult(const Vector &x, Vector &y) const
185  {
186  mfem_error("TimeDependentOperator::ExplicitMult() is not overridden!");
187  }
188
189  /** @brief Perform the action of the implicit part of the operator, F:
190  @a y = F(@a x, @a k, t) where t is the current time.
191
192  Presently, this method is used by some PETSc ODE solvers, for more
193  details, see the PETSc Manual.*/
194  virtual void ImplicitMult(const Vector &x, const Vector &k, Vector &y) const
195  {
196  mfem_error("TimeDependentOperator::ImplicitMult() is not overridden!");
197  }
198
199  /** @brief Perform the action of the operator: @a y = k = f(@a x, t), where
200  k solves the algebraic equation F(@a x, k, t) = G(@a x, t) and t is the
201  current time. */
202  virtual void Mult(const Vector &x, Vector &y) const
203  {
204  mfem_error("TimeDependentOperator::Mult() is not overridden!");
205  }
206
207  /** @brief Solve the equation: @a k = f(@a x + @a dt @a k, t), for the
208  unknown @a k at the current time t.
209
210  For general F and G, the equation for @a k becomes:
211  F(@a x + @a dt @a k, @a k, t) = G(@a x + @a dt @a k, t).
212
213  The input vector @a x corresponds to time index (or cycle) n, while the
214  currently set time, #t, and the result vector @a k correspond to time
215  index n+1. The time step @a dt corresponds to the time interval between
216  cycles n and n+1.
217
218  This method allows for the abstract implementation of some time
219  integration methods, including diagonal implicit Runge-Kutta (DIRK)
220  methods and the backward Euler method in particular.
221
222  If not re-implemented, this method simply generates an error. */
223  virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k)
224  {
225  mfem_error("TimeDependentOperator::ImplicitSolve() is not overridden!");
226  }
227
228  /** @brief Return an Operator representing (dF/dk @a shift + dF/dx) at the
229  given @a x, @a k, and the currently set time.
230
231  Presently, this method is used by some PETSc ODE solvers, for more
232  details, see the PETSc Manual. */
233  virtual Operator& GetImplicitGradient(const Vector &x, const Vector &k,
234  double shift) const
235  {
237  "not overridden!");
238  return const_cast<Operator &>(dynamic_cast<const Operator &>(*this));
239  }
240
241  /** @brief Return an Operator representing dG/dx at the given point @a x and
242  the currently set time.
243
244  Presently, this method is used by some PETSc ODE solvers, for more
245  details, see the PETSc Manual. */
246  virtual Operator& GetExplicitGradient(const Vector &x) const
247  {
249  "not overridden!");
250  return const_cast<Operator &>(dynamic_cast<const Operator &>(*this));
251  }
252
254 };
255
256 /// Base class for solvers
257 class Solver : public Operator
258 {
259 public:
260  /// If true, use the second argument of Mult() as an initial guess.
262
263  /** @brief Initialize a square Solver with size @a s.
264
265  @warning Use a Boolean expression for the second parameter (not an int)
266  to distinguish this call from the general rectangular constructor. */
267  explicit Solver(int s = 0, bool iter_mode = false)
268  : Operator(s) { iterative_mode = iter_mode; }
269
270  /// Initialize a Solver with height @a h and width @a w.
271  Solver(int h, int w, bool iter_mode = false)
272  : Operator(h, w) { iterative_mode = iter_mode; }
273
274  /// Set/update the solver for the given operator.
275  virtual void SetOperator(const Operator &op) = 0;
276 };
277
278
279 /// Identity Operator I: x -> x.
281 {
282 public:
283  /// Create an identity operator of size @a n.
284  explicit IdentityOperator(int n) : Operator(n) { }
285
286  /// Operator application
287  virtual void Mult(const Vector &x, Vector &y) const { y = x; }
288 };
289
290
291 /** @brief The transpose of a given operator. Switches the roles of the methods
292  Mult() and MultTranspose(). */
294 {
295 private:
296  const Operator &A;
297
298 public:
299  /// Construct the transpose of a given operator @a *a.
301  : Operator(a->Width(), a->Height()), A(*a) { }
302
303  /// Construct the transpose of a given operator @a a.
305  : Operator(a.Width(), a.Height()), A(a) { }
306
307  /// Operator application. Apply the transpose of the original Operator.
308  virtual void Mult(const Vector &x, Vector &y) const
309  { A.MultTranspose(x, y); }
310
311  /// Application of the transpose. Apply the original Operator.
312  virtual void MultTranspose(const Vector &x, Vector &y) const
313  { A.Mult(x, y); }
314 };
315
316
317 /// The operator x -> R*A*P*x.
318 class RAPOperator : public Operator
319 {
320 private:
321  const Operator & Rt;
322  const Operator & A;
323  const Operator & P;
324  mutable Vector Px;
325  mutable Vector APx;
326
327 public:
328  /// Construct the RAP operator given R^T, A and P.
329  RAPOperator(const Operator &Rt_, const Operator &A_, const Operator &P_)
330  : Operator(Rt_.Width(), P_.Width()), Rt(Rt_), A(A_), P(P_),
331  Px(P.Height()), APx(A.Height()) { }
332
333  /// Operator application.
334  virtual void Mult(const Vector & x, Vector & y) const
335  { P.Mult(x, Px); A.Mult(Px, APx); Rt.MultTranspose(APx, y); }
336
337  /// Application of the transpose.
338  virtual void MultTranspose(const Vector & x, Vector & y) const
339  { Rt.Mult(x, APx); A.MultTranspose(APx, Px); P.MultTranspose(Px, y); }
340 };
341
342
343 /// General triple product operator x -> A*B*C*x, with ownership of the factors.
345 {
346  Operator *A;
347  Operator *B;
348  Operator *C;
349  bool ownA, ownB, ownC;
350  mutable Vector t1, t2;
351
352 public:
354  bool ownA, bool ownB, bool ownC)
355  : Operator(A->Height(), C->Width())
356  , A(A), B(B), C(C)
357  , ownA(ownA), ownB(ownB), ownC(ownC)
358  , t1(C->Height()), t2(B->Height())
359  {}
360
361  virtual void Mult(const Vector &x, Vector &y) const
362  { C->Mult(x, t1); B->Mult(t1, t2); A->Mult(t2, y); }
363
364  virtual void MultTranspose(const Vector &x, Vector &y) const
365  { A->MultTranspose(x, t2); B->MultTranspose(t2, t1); C->MultTranspose(t1, y); }
366
368  {
369  if (ownA) { delete A; }
370  if (ownB) { delete B; }
371  if (ownC) { delete C; }
372  }
373 };
374
375
376 /** @brief Square Operator for imposing essential boundary conditions using only
377  the action, Mult(), of a given unconstrained Operator.
378
379  Square operator constrained by fixing certain entries in the solution to
380  given "essential boundary condition" values. This class is used by the
381  general, matrix-free system formulation of Operator::FormLinearSystem. */
383 {
384 protected:
385  Array<int> constraint_list; ///< List of constrained indices/dofs.
386  Operator *A; ///< The unconstrained Operator.
387  bool own_A; ///< Ownership flag for A.
388  mutable Vector z, w; ///< Auxiliary vectors.
389
390 public:
391  /** @brief Constructor from a general Operator and a list of essential
392  indices/dofs.
393
394  Specify the unconstrained operator @a *A and a @a list of indices to
395  constrain, i.e. each entry @a list[i] represents an essential-dof. If the
396  ownership flag @a own_A is true, the operator @a *A will be destroyed
397  when this object is destroyed. */
398  ConstrainedOperator(Operator *A, const Array<int> &list, bool own_A = false);
399
400  /** @brief Eliminate "essential boundary condition" values specified in @a x
401  from the given right-hand side @a b.
402
403  Performs the following steps:
404
405  z = A((0,x_b)); b_i -= z_i; b_b = x_b;
406
407  where the "_b" subscripts denote the essential (boundary) indices/dofs of
408  the vectors, and "_i" -- the rest of the entries. */
409  void EliminateRHS(const Vector &x, Vector &b) const;
410
411  /** @brief Constrained operator action.
412
413  Performs the following steps:
414
415  z = A((x_i,0)); y_i = z_i; y_b = x_b;
416
417  where the "_b" subscripts denote the essential (boundary) indices/dofs of
418  the vectors, and "_i" -- the rest of the entries. */
419  virtual void Mult(const Vector &x, Vector &y) const;
420
421  /// Destructor: destroys the unconstrained Operator @a A if @a own_A is true.
422  virtual ~ConstrainedOperator() { if (own_A) { delete A; } }
423 };
424
425 }
426
427 #endif
virtual double GetTime() const
Definition: operator.hpp:167
bool isImplicit() const
True if type is IMPLICIT or HOMOGENEOUS.
Definition: operator.hpp:175
bool isHomogeneous() const
True if type is HOMOGENEOUS.
Definition: operator.hpp:177
Solver(int s=0, bool iter_mode=false)
Initialize a square Solver with size s.
Definition: operator.hpp:267
Array< int > constraint_list
List of constrained indices/dofs.
Definition: operator.hpp:385
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
Definition: operator.hpp:334
Type type
Describes the form of the TimeDependentOperator.
Definition: operator.hpp:152
Base abstract class for time dependent operators.
Definition: operator.hpp:140
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:170
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:353
virtual void MultTranspose(const Vector &x, Vector &y) const
Application of the transpose. Apply the original Operator.
Definition: operator.hpp:312
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:261
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:246
bool own_A
Ownership flag for A.
Definition: operator.hpp:387
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:271
double t
Current time.
Definition: operator.hpp:151
virtual void Mult(const Vector &x, Vector &y) const
Operator application. Apply the transpose of the original Operator.
Definition: operator.hpp:308
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:233
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:304
virtual void MultTranspose(const Vector &x, Vector &y) const
Application of the transpose.
Definition: operator.hpp:338
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:422
This type assumes F(x,k,t) = k, i.e. k = f(x,t) = G(x,t).
Definition: operator.hpp:145
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:163
The operator x -&gt; R*A*P*x.
Definition: operator.hpp:318
bool isExplicit() const
True if type is EXPLICIT.
Definition: operator.hpp:173
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:223
RAPOperator(const Operator &Rt_, const Operator &A_, const Operator &P_)
Construct the RAP operator given R^T, A and P.
Definition: operator.hpp:329
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:146
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: operator.hpp:361
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:284
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:157
Vector w
Auxiliary vectors.
Definition: operator.hpp:388
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:202
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:184
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:106
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:293
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:147
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:344
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:287
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:36
Identity Operator I: x -&gt; x.
Definition: operator.hpp:280
TransposeOperator(const Operator *a)
Construct the transpose of a given operator *a.
Definition: operator.hpp:300
Base class for solvers.
Definition: operator.hpp:257
Operator * A
The unconstrained Operator.
Definition: operator.hpp:386
Square Operator for imposing essential boundary conditions using only the action, Mult()...
Definition: operator.hpp:382
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:364
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:194