MFEM  v4.1.0 Finite element discretization library
operator.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, 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_OPERATOR
13 #define MFEM_OPERATOR
14
15 #include "vector.hpp"
16
17 namespace mfem
18 {
19
20 class ConstrainedOperator;
21 class RectangularConstrainedOperator;
22
23 /// Abstract operator
24 class Operator
25 {
26 protected:
27  int height; ///< Dimension of the output / number of rows in the matrix.
28  int width; ///< Dimension of the input / number of columns in the matrix.
29
30  /// see FormSystemOperator()
32  const Array<int> &ess_tdof_list, ConstrainedOperator* &Aout);
33
34  /// see FormRectangularSystemOperator()
36  const Array<int> &trial_tdof_list,
37  const Array<int> &test_tdof_list,
39
40  /// Returns RAP Operator of this, taking in input/output Prolongation matrices
41  Operator *SetupRAP(const Operator *Pi, const Operator *Po);
42
43 public:
44  /// Initializes memory for true vectors of linear system
45  void InitTVectors(const Operator *Po, const Operator *Ri, const Operator *Pi,
46  Vector &x, Vector &b,
47  Vector &X, Vector &B) const;
48
49  /// Construct a square Operator with given size s (default 0).
50  explicit Operator(int s = 0) { height = width = s; }
51
52  /** @brief Construct an Operator with the given height (output size) and
53  width (input size). */
54  Operator(int h, int w) { height = h; width = w; }
55
56  /// Get the height (size of output) of the Operator. Synonym with NumRows().
57  inline int Height() const { return height; }
58  /** @brief Get the number of rows (size of output) of the Operator. Synonym
59  with Height(). */
60  inline int NumRows() const { return height; }
61
62  /// Get the width (size of input) of the Operator. Synonym with NumCols().
63  inline int Width() const { return width; }
64  /** @brief Get the number of columns (size of input) of the Operator. Synonym
65  with Width(). */
66  inline int NumCols() const { return width; }
67
68  /// Return the MemoryClass preferred by the Operator.
69  /** This is the MemoryClass that will be used to access the input and output
70  vectors in the Mult() and MultTranspose() methods.
71
72  For example, classes using the MFEM_FORALL macro for implementation can
73  return the value returned by Device::GetMemoryClass().
74
75  The default implementation of this method in class Operator returns
76  MemoryClass::HOST. */
77  virtual MemoryClass GetMemoryClass() const { return MemoryClass::HOST; }
78
79  /// Operator application: y=A(x).
80  virtual void Mult(const Vector &x, Vector &y) const = 0;
81
82  /** @brief Action of the transpose operator: y=A^t(x). The default behavior
83  in class Operator is to generate an error. */
84  virtual void MultTranspose(const Vector &x, Vector &y) const
85  { mfem_error("Operator::MultTranspose() is not overloaded!"); }
86
87  /** @brief Evaluate the gradient operator at the point @a x. The default
88  behavior in class Operator is to generate an error. */
89  virtual Operator &GetGradient(const Vector &x) const
90  {
92  return const_cast<Operator &>(*this);
93  }
94
95  /** @brief Prolongation operator from linear algebra (linear system) vectors,
96  to input vectors for the operator. NULL means identity. */
97  virtual const Operator *GetProlongation() const { return NULL; }
98  /** @brief Restriction operator from input vectors for the operator to linear
99  algebra (linear system) vectors. NULL means identity. */
100  virtual const Operator *GetRestriction() const { return NULL; }
101  /** @brief Prolongation operator from linear algebra (linear system) vectors,
102  to output vectors for the operator. NULL means identity. */
103  virtual const Operator *GetOutputProlongation() const
104  {
105  return GetProlongation(); // Assume square unless specialized
106  }
107  /** @brief Restriction operator from output vectors for the operator to linear
108  algebra (linear system) vectors. NULL means identity. */
109  virtual const Operator *GetOutputRestriction() const
110  {
111  return GetRestriction(); // Assume square unless specialized
112  }
113
114  /** @brief Form a constrained linear system using a matrix-free approach.
115
116  Assuming square operator, form the operator linear system A(X)=B,
117  corresponding to it and the right-hand side @a b, by applying any
118  necessary transformations such as: parallel assembly, conforming
119  constraints for non-conforming AMR and eliminating boundary conditions.
120  @note Static condensation and hybridization are not supported for general
121  operators (cf. the analogous methods BilinearForm::FormLinearSystem() and
122  ParBilinearForm::FormLinearSystem()).
123
124  The constraints are specified through the prolongation P from
125  GetProlongation(), and restriction R from GetRestriction() methods, which
126  are e.g. available through the (parallel) finite element space of any
127  (parallel) bilinear form operator. We assume that the operator is square,
128  using the same input and output space, so we have: A(X)=[P^t (*this)
129  P](X), B=P^t(b), and X=R(x).
130
131  The vector @a x must contain the essential boundary condition values.
132  These are eliminated through the ConstrainedOperator class and the vector
133  @a X is initialized by setting its essential entries to the boundary
134  conditions and all other entries to zero (@a copy_interior == 0) or
135  copied from @a x (@a copy_interior != 0).
136
137  After solving the system A(X)=B, the (finite element) solution @a x can
138  be recovered by calling Operator::RecoverFEMSolution() with the same
139  vectors @a X, @a b, and @a x.
140
141  @note The caller is responsible for destroying the output operator @a A!
142  @note If there are no transformations, @a X simply reuses the data of @a
143  x. */
144  void FormLinearSystem(const Array<int> &ess_tdof_list,
145  Vector &x, Vector &b,
146  Operator* &A, Vector &X, Vector &B,
147  int copy_interior = 0);
148
149  /** @brief Form a column-constrained linear system using a matrix-free approach.
150
151  Form the operator linear system A(X)=B corresponding to the operator
152  and the right-hand side @a b, by applying any necessary transformations
153  such as: parallel assembly, conforming constraints for non-conforming AMR
154  and eliminating boundary conditions. @note Static condensation and
155  hybridization are not supported for general operators (cf. the method
156  MixedBilinearForm::FormRectangularLinearSystem())
157
158  The constraints are specified through the input prolongation Pi from
159  GetProlongation(), and output restriction Ro from GetOutputRestriction()
160  methods, which are e.g. available through the (parallel) finite element
161  spaces of any (parallel) mixed bilinear form operator. So we have:
162  A(X)=[Ro (*this) Pi](X), B=Ro(b), and X=Pi^T(x).
163
164  The vector @a x must contain the essential boundary condition values.
165  The "columns" in this operator corresponding to these values are
166  eliminated through the RectangularConstrainedOperator class.
167
168  After solving the system A(X)=B, the (finite element) solution @a x can
169  be recovered by calling Operator::RecoverFEMSolution() with the same
170  vectors @a X, @a b, and @a x.
171
172  @note The caller is responsible for destroying the output operator @a A!
173  @note If there are no transformations, @a X simply reuses the data of @a
174  x. */
175  void FormRectangularLinearSystem(const Array<int> &trial_tdof_list,
176  const Array<int> &test_tdof_list,
177  Vector &x, Vector &b,
178  Operator* &A, Vector &X, Vector &B);
179
180  /** @brief Reconstruct a solution vector @a x (e.g. a GridFunction) from the
181  solution @a X of a constrained linear system obtained from
182  Operator::FormLinearSystem() or Operator::FormRectangularLinearSystem().
183
184  Call this method after solving a linear system constructed using
185  Operator::FormLinearSystem() to recover the solution as an input vector,
186  @a x, for this Operator (presumably a finite element grid function). This
187  method has identical signature to the analogous method for bilinear
188  forms, though currently @a b is not used in the implementation. */
189  virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x);
190
191  /** @brief Return in @a A a parallel (on truedofs) version of this square
192  operator.
193
194  This returns the same operator as FormLinearSystem(), but does without
195  the transformations of the right-hand side and initial guess. */
196  void FormSystemOperator(const Array<int> &ess_tdof_list,
197  Operator* &A);
198
199  /** @brief Return in @a A a parallel (on truedofs) version of this
200  rectangular operator (including constraints).
201
202  This returns the same operator as FormRectangularLinearSystem(), but does
203  without the transformations of the right-hand side. */
204  void FormRectangularSystemOperator(const Array<int> &trial_tdof_list,
205  const Array<int> &test_tdof_list,
206  Operator* &A);
207
208  /** @brief Return in @a A a parallel (on truedofs) version of this
209  rectangular operator.
210
211  This is similar to FormSystemOperator(), but for dof-to-dof mappings
212  (discrete linear operators), which can also correspond to rectangular
213  matrices. The user should provide specializations of GetProlongation()
214  for the input dofs and GetOutputRestriction() for the output dofs in
215  their Operator implementation that are appropriate for the two spaces the
216  Operator maps between. These are e.g. available through the (parallel)
217  finite element space of any (parallel) bilinear form operator. We have:
218  A(X)=[Rout (*this) Pin](X). */
219  void FormDiscreteOperator(Operator* &A);
220
221  /// Prints operator with input size n and output size m in Matlab format.
222  void PrintMatlab(std::ostream & out, int n = 0, int m = 0) const;
223
224  /// Virtual destructor.
225  virtual ~Operator() { }
226
227  /// Enumeration defining IDs for some classes derived from Operator.
228  /** This enumeration is primarily used with class OperatorHandle. */
229  enum Type
230  {
231  ANY_TYPE, ///< ID for the base class Operator, i.e. any type.
232  MFEM_SPARSEMAT, ///< ID for class SparseMatrix.
233  Hypre_ParCSR, ///< ID for class HypreParMatrix.
234  PETSC_MATAIJ, ///< ID for class PetscParMatrix, MATAIJ format.
235  PETSC_MATIS, ///< ID for class PetscParMatrix, MATIS format.
236  PETSC_MATSHELL, ///< ID for class PetscParMatrix, MATSHELL format.
237  PETSC_MATNEST, ///< ID for class PetscParMatrix, MATNEST format.
238  PETSC_MATHYPRE, ///< ID for class PetscParMatrix, MATHYPRE format.
239  PETSC_MATGENERIC, ///< ID for class PetscParMatrix, unspecified format.
240  Complex_Operator, ///< ID for class ComplexOperator.
241  MFEM_ComplexSparseMat, ///< ID for class ComplexSparseMatrix.
242  Complex_Hypre_ParCSR ///< ID for class ComplexHypreParMatrix.
243  };
244
245  /// Return the type ID of the Operator class.
246  /** This method is intentionally non-virtual, so that it returns the ID of
247  the specific pointer or reference type used when calling this method. If
248  not overridden by derived classes, they will automatically use the type ID
249  of the base Operator class, ANY_TYPE. */
250  Type GetType() const { return ANY_TYPE; }
251 };
252
253
254 /// Base abstract class for first order time dependent operators.
255 /** Operator of the form: (x,t) -> f(x,t), where k = f(x,t) generally solves the
256  algebraic equation F(x,k,t) = G(x,t). The functions F and G represent the
257  _implicit_ and _explicit_ parts of the operator, respectively. For explicit
258  operators, F(x,k,t) = k, so f(x,t) = G(x,t). */
260 {
261 public:
262  enum Type
263  {
264  EXPLICIT, ///< This type assumes F(x,k,t) = k, i.e. k = f(x,t) = G(x,t).
265  IMPLICIT, ///< This is the most general type, no assumptions on F and G.
266  HOMOGENEOUS ///< This type assumes that G(x,t) = 0.
267  };
268
269  /// Evaluation mode. See SetEvalMode() for details.
270  enum EvalMode
271  {
272  /** Normal evaluation. */
274  /** Assuming additive split, f(x,t) = f1(x,t) + f2(x,t), evaluate the
275  first term, f1. */
277  /** Assuming additive split, f(x,t) = f1(x,t) + f2(x,t), evaluate the
278  second term, f2. */
280  };
281
282 protected:
283  double t; ///< Current time.
284  Type type; ///< Describes the form of the TimeDependentOperator.
285  EvalMode eval_mode; ///< Current evaluation mode.
286
287 public:
288  /** @brief Construct a "square" TimeDependentOperator y = f(x,t), where x and
289  y have the same dimension @a n. */
290  explicit TimeDependentOperator(int n = 0, double t_ = 0.0,
291  Type type_ = EXPLICIT)
292  : Operator(n) { t = t_; type = type_; eval_mode = NORMAL; }
293
294  /** @brief Construct a TimeDependentOperator y = f(x,t), where x and y have
295  dimensions @a w and @a h, respectively. */
296  TimeDependentOperator(int h, int w, double t_ = 0.0, Type type_ = EXPLICIT)
297  : Operator(h, w) { t = t_; type = type_; eval_mode = NORMAL; }
298
299  /// Read the currently set time.
300  virtual double GetTime() const { return t; }
301
302  /// Set the current time.
303  virtual void SetTime(const double _t) { t = _t; }
304
305  /// True if #type is #EXPLICIT.
306  bool isExplicit() const { return (type == EXPLICIT); }
307  /// True if #type is #IMPLICIT or #HOMOGENEOUS.
308  bool isImplicit() const { return !isExplicit(); }
309  /// True if #type is #HOMOGENEOUS.
310  bool isHomogeneous() const { return (type == HOMOGENEOUS); }
311
312  /// Return the current evaluation mode. See SetEvalMode() for details.
313  EvalMode GetEvalMode() const { return eval_mode; }
314
315  /// Set the evaluation mode of the time-dependent operator.
316  /** The evaluation mode is a switch that allows time-stepping methods to
317  request evaluation of separate components/terms of the time-dependent
318  operator. For example, IMEX methods typically assume additive split of
319  the operator: f(x,t) = f1(x,t) + f2(x,t) and they rely on the ability to
320  evaluate the two terms separately.
321
322  Generally, setting the evaluation mode should affect the behavior of all
323  evaluation-related methods in the class, such as Mult(), ImplicitSolve(),
324  etc. However, the exact list of methods that need to support a specific
325  mode will depend on the used time-stepping method. */
326  virtual void SetEvalMode(const EvalMode new_eval_mode)
327  { eval_mode = new_eval_mode; }
328
329  /** @brief Perform the action of the explicit part of the operator, G:
330  @a y = G(@a x, t) where t is the current time.
331
332  Presently, this method is used by some PETSc ODE solvers, for more
333  details, see the PETSc Manual. */
334  virtual void ExplicitMult(const Vector &x, Vector &y) const;
335
336  /** @brief Perform the action of the implicit part of the operator, F:
337  @a y = F(@a x, @a k, t) where t is the current time.
338
339  Presently, this method is used by some PETSc ODE solvers, for more
340  details, see the PETSc Manual.*/
341  virtual void ImplicitMult(const Vector &x, const Vector &k, Vector &y) const;
342
343  /** @brief Perform the action of the operator: @a y = k = f(@a x, t), where
344  k solves the algebraic equation F(@a x, k, t) = G(@a x, t) and t is the
345  current time. */
346  virtual void Mult(const Vector &x, Vector &y) const;
347
348  /** @brief Solve the equation: @a k = f(@a x + @a dt @a k, t), for the
349  unknown @a k at the current time t.
350
351  For general F and G, the equation for @a k becomes:
352  F(@a x + @a dt @a k, @a k, t) = G(@a x + @a dt @a k, t).
353
354  The input vector @a x corresponds to time index (or cycle) n, while the
355  currently set time, #t, and the result vector @a k correspond to time
356  index n+1. The time step @a dt corresponds to the time interval between
357  cycles n and n+1.
358
359  This method allows for the abstract implementation of some time
360  integration methods, including diagonal implicit Runge-Kutta (DIRK)
361  methods and the backward Euler method in particular.
362
363  If not re-implemented, this method simply generates an error. */
364  virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k);
365
366  /** @brief Return an Operator representing (dF/dk @a shift + dF/dx) at the
367  given @a x, @a k, and the currently set time.
368
369  Presently, this method is used by some PETSc ODE solvers, for more
370  details, see the PETSc Manual. */
371  virtual Operator& GetImplicitGradient(const Vector &x, const Vector &k,
372  double shift) const;
373
374  /** @brief Return an Operator representing dG/dx at the given point @a x and
375  the currently set time.
376
377  Presently, this method is used by some PETSc ODE solvers, for more
378  details, see the PETSc Manual. */
379  virtual Operator& GetExplicitGradient(const Vector &x) const;
380
381  /** @brief Setup the ODE linear system \f$A(x,t) = (I - gamma J) \f$ or
382  \f$A = (M - gamma J) \f$, where \f$J(x,t) = \frac{df}{dt(x,t)} \f$.
383
384  @param[in] x The state at which \f$A(x,t)\f$ should be evaluated.
385  @param[in] fx The current value of the ODE rhs function, \f$f(x,t)\f$.
386  @param[in] jok Flag indicating if the Jacobian should be updated.
387  @param[out] jcur Flag to signal if the Jacobian was updated.
388  @param[in] gamma The scaled time step value.
389
390  If not re-implemented, this method simply generates an error.
391
392  Presently, this method is used by SUNDIALS ODE solvers, for more
393  details, see the SUNDIALS User Guides. */
394  virtual int SUNImplicitSetup(const Vector &x, const Vector &fx,
395  int jok, int *jcur, double gamma);
396
397  /** @brief Solve the ODE linear system \f$A x = b \f$ as setup by
398  the method SUNImplicitSetup().
399
400  @param[in] b The linear system right-hand side.
401  @param[in,out] x On input, the initial guess. On output, the solution.
402  @param[in] tol Linear solve tolerance.
403
404  If not re-implemented, this method simply generates an error.
405
406  Presently, this method is used by SUNDIALS ODE solvers, for more
407  details, see the SUNDIALS User Guides. */
408  virtual int SUNImplicitSolve(const Vector &b, Vector &x, double tol);
409
410  /** @brief Setup the mass matrix in the ODE system \f$M y' = f(y,t) \f$ .
411
412  If not re-implemented, this method simply generates an error.
413
414  Presently, this method is used by SUNDIALS ARKStep integrator, for more
415  details, see the ARKode User Guide. */
416  virtual int SUNMassSetup();
417
418  /** @brief Solve the mass matrix linear system \f$M x = b \f$
419  as setup by the method SUNMassSetup().
420
421  @param[in] b The linear system right-hand side.
422  @param[in,out] x On input, the initial guess. On output, the solution.
423  @param[in] tol Linear solve tolerance.
424
425  If not re-implemented, this method simply generates an error.
426
427  Presently, this method is used by SUNDIALS ARKStep integrator, for more
428  details, see the ARKode User Guide. */
429  virtual int SUNMassSolve(const Vector &b, Vector &x, double tol);
430
431  /** @brief Compute the mass matrix-vector product \f$v = M x \f$ .
432
433  @param[in] x The vector to multiply.
434  @param[out] v The result of the matrix-vector product.
435
436  If not re-implemented, this method simply generates an error.
437
438  Presently, this method is used by SUNDIALS ARKStep integrator, for more
439  details, see the ARKode User Guide. */
440  virtual int SUNMassMult(const Vector &x, Vector &v);
441
443 };
444
445 /// Base abstract class for second order time dependent operators.
446 /** Operator of the form: (x,dxdt,t) -> f(x,dxdt,t), where k = f(x,dxdt,t)
447  generally solves the algebraic equation F(x,dxdt,k,t) = G(x,dxdt,t).
448  The functions F and G represent the_implicit_ and _explicit_ parts of
449  the operator, respectively. For explicit operators,
450  F(x,dxdt,k,t) = k, so f(x,dxdt,t) = G(x,dxdt,t). */
452 {
453 public:
454  /** @brief Construct a "square" SecondOrderTimeDependentOperator
455  y = f(x,dxdt,t), where x, dxdt and y have the same dimension @a n. */
456  explicit SecondOrderTimeDependentOperator(int n = 0, double t_ = 0.0,
457  Type type_ = EXPLICIT)
458  : TimeDependentOperator(n, t_,type_) { }
459
460  /** @brief Construct a SecondOrderTimeDependentOperator y = f(x,dxdt,t),
461  where x, dxdt and y have the same dimension @a n. */
462  SecondOrderTimeDependentOperator(int h, int w, double t_ = 0.0,
463  Type type_ = EXPLICIT)
464  : TimeDependentOperator(h, w, t_,type_) { }
465
467
468  /** @brief Perform the action of the operator: @a y = k = f(@a x,@ dxdt, t),
469  where k solves the algebraic equation
470  F(@a x,@ dxdt, k, t) = G(@a x,@ dxdt, t) and t is the current time. */
471  virtual void Mult(const Vector &x, const Vector &dxdt, Vector &y) const;
472
474  /** @brief Solve the equation:
475  @a k = f(@a x + 1/2 @a dt0^2 @a k, @a dxdt + @a dt1 @a k, t), for the
476  unknown @a k at the current time t.
477
478  For general F and G, the equation for @a k becomes:
479  F(@a x + 1/2 @a dt0^2 @a k, @a dxdt + @a dt1 @a k, t)
480  = G(@a x + 1/2 @a dt0^2 @a k, @a dxdt + @a dt1 @a k, t).
481
482  The input vector @a x corresponds to time index (or cycle) n, while the
483  currently set time, #t, and the result vector @a k correspond to time
484  index n+1. The time step @a dt corresponds to the time interval between
485  cycles n and n+1.
486
487  This method allows for the abstract implementation of some time
488  integration methods.
489
490  If not re-implemented, this method simply generates an error. */
491  virtual void ImplicitSolve(const double dt0, const double dt1,
492  const Vector &x, const Vector &dxdt, Vector &k);
493
494
496 };
497
498
499 /// Base class for solvers
500 class Solver : public Operator
501 {
502 public:
503  /// If true, use the second argument of Mult() as an initial guess.
505
506  /** @brief Initialize a square Solver with size @a s.
507
508  @warning Use a Boolean expression for the second parameter (not an int)
509  to distinguish this call from the general rectangular constructor. */
510  explicit Solver(int s = 0, bool iter_mode = false)
511  : Operator(s) { iterative_mode = iter_mode; }
512
513  /// Initialize a Solver with height @a h and width @a w.
514  Solver(int h, int w, bool iter_mode = false)
515  : Operator(h, w) { iterative_mode = iter_mode; }
516
517  /// Set/update the solver for the given operator.
518  virtual void SetOperator(const Operator &op) = 0;
519 };
520
521
522 /// Identity Operator I: x -> x.
524 {
525 public:
526  /// Create an identity operator of size @a n.
527  explicit IdentityOperator(int n) : Operator(n) { }
528
529  /// Operator application
530  virtual void Mult(const Vector &x, Vector &y) const { y = x; }
531
532  /// Application of the transpose
533  virtual void MultTranspose(const Vector &x, Vector &y) const { y = x; }
534 };
535
536 /// Returns true if P is the identity prolongation, i.e. if it is either NULL or
537 /// an IdentityOperator.
538 inline bool IsIdentityProlongation(const Operator *P)
539 {
540  return !P || dynamic_cast<const IdentityOperator*>(P);
541 }
542
543 /// Scaled Operator B: x -> a A(x).
544 class ScaledOperator : public Operator
545 {
546 private:
547  const Operator &A_;
548  double a_;
549
550 public:
551  /// Create an operator which is a scalar multiple of A.
552  explicit ScaledOperator(const Operator *A, double a)
553  : Operator(A->Width(), A->Height()), A_(*A), a_(a) { }
554
555  /// Operator application
556  virtual void Mult(const Vector &x, Vector &y) const
557  { A_.Mult(x, y); y *= a_; }
558 };
559
560
561 /** @brief The transpose of a given operator. Switches the roles of the methods
562  Mult() and MultTranspose(). */
564 {
565 private:
566  const Operator &A;
567
568 public:
569  /// Construct the transpose of a given operator @a *a.
571  : Operator(a->Width(), a->Height()), A(*a) { }
572
573  /// Construct the transpose of a given operator @a a.
575  : Operator(a.Width(), a.Height()), A(a) { }
576
577  /// Operator application. Apply the transpose of the original Operator.
578  virtual void Mult(const Vector &x, Vector &y) const
579  { A.MultTranspose(x, y); }
580
581  /// Application of the transpose. Apply the original Operator.
582  virtual void MultTranspose(const Vector &x, Vector &y) const
583  { A.Mult(x, y); }
584 };
585
586
587 /// General product operator: x -> (A*B)(x) = A(B(x)).
588 class ProductOperator : public Operator
589 {
590  const Operator *A, *B;
591  bool ownA, ownB;
592  mutable Vector z;
593
594 public:
595  ProductOperator(const Operator *A, const Operator *B, bool ownA, bool ownB);
596
597  virtual void Mult(const Vector &x, Vector &y) const
598  { B->Mult(x, z); A->Mult(z, y); }
599
600  virtual void MultTranspose(const Vector &x, Vector &y) const
601  { A->MultTranspose(x, z); B->MultTranspose(z, y); }
602
603  virtual ~ProductOperator();
604 };
605
606
607 /// The operator x -> R*A*P*x constructed through the actions of R^T, A and P
608 class RAPOperator : public Operator
609 {
610 private:
611  const Operator & Rt;
612  const Operator & A;
613  const Operator & P;
614  mutable Vector Px;
615  mutable Vector APx;
616  MemoryClass mem_class;
617
618 public:
619  /// Construct the RAP operator given R^T, A and P.
620  RAPOperator(const Operator &Rt_, const Operator &A_, const Operator &P_);
621
622  virtual MemoryClass GetMemoryClass() const { return mem_class; }
623
624  /// Operator application.
625  virtual void Mult(const Vector & x, Vector & y) const
626  { P.Mult(x, Px); A.Mult(Px, APx); Rt.MultTranspose(APx, y); }
627
628  /// Application of the transpose.
629  virtual void MultTranspose(const Vector & x, Vector & y) const
630  { Rt.Mult(x, APx); A.MultTranspose(APx, Px); P.MultTranspose(Px, y); }
631 };
632
633
634 /// General triple product operator x -> A*B*C*x, with ownership of the factors.
636 {
637  const Operator *A;
638  const Operator *B;
639  const Operator *C;
640  bool ownA, ownB, ownC;
641  mutable Vector t1, t2;
642  MemoryClass mem_class;
643
644 public:
645  TripleProductOperator(const Operator *A, const Operator *B,
646  const Operator *C, bool ownA, bool ownB, bool ownC);
647
648  virtual MemoryClass GetMemoryClass() const { return mem_class; }
649
650  virtual void Mult(const Vector &x, Vector &y) const
651  { C->Mult(x, t1); B->Mult(t1, t2); A->Mult(t2, y); }
652
653  virtual void MultTranspose(const Vector &x, Vector &y) const
654  { A->MultTranspose(x, t2); B->MultTranspose(t2, t1); C->MultTranspose(t1, y); }
655
656  virtual ~TripleProductOperator();
657 };
658
659
660 /** @brief Square Operator for imposing essential boundary conditions using only
661  the action, Mult(), of a given unconstrained Operator.
662
663  Square operator constrained by fixing certain entries in the solution to
664  given "essential boundary condition" values. This class is used by the
665  general, matrix-free system formulation of Operator::FormLinearSystem. */
667 {
668 protected:
669  Array<int> constraint_list; ///< List of constrained indices/dofs.
670  Operator *A; ///< The unconstrained Operator.
671  bool own_A; ///< Ownership flag for A.
672  mutable Vector z, w; ///< Auxiliary vectors.
674
675 public:
676  /** @brief Constructor from a general Operator and a list of essential
677  indices/dofs.
678
679  Specify the unconstrained operator @a *A and a @a list of indices to
680  constrain, i.e. each entry @a list[i] represents an essential-dof. If the
681  ownership flag @a own_A is true, the operator @a *A will be destroyed
682  when this object is destroyed. */
683  ConstrainedOperator(Operator *A, const Array<int> &list, bool own_A = false);
684
685  /// Returns the type of memory in which the solution and temporaries are stored.
686  virtual MemoryClass GetMemoryClass() const { return mem_class; }
687
688  /** @brief Eliminate "essential boundary condition" values specified in @a x
689  from the given right-hand side @a b.
690
691  Performs the following steps:
692
693  z = A((0,x_b)); b_i -= z_i; b_b = x_b;
694
695  where the "_b" subscripts denote the essential (boundary) indices/dofs of
696  the vectors, and "_i" -- the rest of the entries. */
697  void EliminateRHS(const Vector &x, Vector &b) const;
698
699  /** @brief Constrained operator action.
700
701  Performs the following steps:
702
703  z = A((x_i,0)); y_i = z_i; y_b = x_b;
704
705  where the "_b" subscripts denote the essential (boundary) indices/dofs of
706  the vectors, and "_i" -- the rest of the entries. */
707  virtual void Mult(const Vector &x, Vector &y) const;
708
709  /// Destructor: destroys the unconstrained Operator, if owned.
710  virtual ~ConstrainedOperator() { if (own_A) { delete A; } }
711 };
712
713 /** @brief Rectangular Operator for imposing essential boundary conditions on
714  the input space using only the action, Mult(), of a given unconstrained
715  Operator.
716
717  Rectangular operator constrained by fixing certain entries in the solution
718  to given "essential boundary condition" values. This class is used by the
719  general matrix-free formulation of Operator::FormRectangularLinearSystem. */
721 {
722 protected:
725  bool own_A;
726  mutable Vector z, w;
728
729 public:
730  /** @brief Constructor from a general Operator and a list of essential
731  indices/dofs.
732
733  Specify the unconstrained operator @a *A and two lists of indices to
734  constrain, i.e. each entry @a trial_list[i] represents an essential trial
735  dof. If the ownership flag @a own_A is true, the operator @a *A will be
736  destroyed when this object is destroyed. */
738  const Array<int> &test_list, bool own_A = false);
739  /// Returns the type of memory in which the solution and temporaries are stored.
740  virtual MemoryClass GetMemoryClass() const { return mem_class; }
741  /** @brief Eliminate columns corresponding to "essential boundary condition"
742  values specified in @a x from the given right-hand side @a b.
743
744  Performs the following steps:
745
746  b -= A((0,x_b));
747  b_j = 0
748
749  where the "_b" subscripts denote the essential (boundary) indices and the
750  "_j" subscript denotes the essential test indices */
751  void EliminateRHS(const Vector &x, Vector &b) const;
752  /** @brief Rectangular-constrained operator action.
753
754  Performs the following steps:
755
756  y = A((x_i,0));
757  y_j = 0
758
759  where the "_i" subscripts denote all the nonessential (boundary) trial
760  indices and the "_j" subscript denotes the essential test indices */
761  virtual void Mult(const Vector &x, Vector &y) const;
762  virtual ~RectangularConstrainedOperator() { if (own_A) { delete A; } }
763 };
764
765 }
766
767 #endif
virtual double GetTime() const
Definition: operator.hpp:300
virtual void Mult(const Vector &x, const Vector &dxdt, Vector &y) const
Perform the action of the operator: y = k = f(x,@ dxdt, t), where k solves the algebraic equation F(x...
Definition: operator.cpp:284
void FormRectangularSystemOperator(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, Operator *&A)
Return in A a parallel (on truedofs) version of this rectangular operator (including constraints)...
Definition: operator.cpp:172
bool isImplicit() const
True if type is IMPLICIT or HOMOGENEOUS.
Definition: operator.hpp:308
bool isHomogeneous() const
True if type is HOMOGENEOUS.
Definition: operator.hpp:310
virtual MemoryClass GetMemoryClass() const
Returns the type of memory in which the solution and temporaries are stored.
Definition: operator.hpp:740
Solver(int s=0, bool iter_mode=false)
Initialize a square Solver with size s.
Definition: operator.hpp:510
virtual ~ProductOperator()
Definition: operator.cpp:320
Array< int > constraint_list
List of constrained indices/dofs.
Definition: operator.hpp:669
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
Definition: operator.hpp:625
Type type
Describes the form of the TimeDependentOperator.
Definition: operator.hpp:284
Base abstract class for first order time dependent operators.
Definition: operator.hpp:259
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.cpp:230
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:89
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:63
ID for class ComplexHypreParMatrix.
Definition: operator.hpp:242
virtual void Mult(const Vector &x, Vector &y) const
Constrained operator action.
Definition: operator.cpp:446
virtual void SetTime(const double _t)
Set the current time.
Definition: operator.hpp:303
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.cpp:244
virtual MemoryClass GetMemoryClass() const
Return the MemoryClass preferred by the Operator.
Definition: operator.hpp:648
virtual const Operator * GetOutputProlongation() const
Prolongation operator from linear algebra (linear system) vectors, to output vectors for the operator...
Definition: operator.hpp:103
void EliminateRHS(const Vector &x, Vector &b) const
Eliminate columns corresponding to &quot;essential boundary condition&quot; values specified in x from the give...
Definition: operator.cpp:493
void FormDiscreteOperator(Operator *&A)
Return in A a parallel (on truedofs) version of this rectangular operator.
Definition: operator.cpp:181
virtual void MultTranspose(const Vector &x, Vector &y) const
Application of the transpose.
Definition: operator.hpp:533
EvalMode eval_mode
Current evaluation mode.
Definition: operator.hpp:285
ProductOperator(const Operator *A, const Operator *B, bool ownA, bool ownB)
Definition: operator.cpp:301
ID for class PetscParMatrix, MATHYPRE format.
Definition: operator.hpp:238
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
virtual void MultTranspose(const Vector &x, Vector &y) const
Application of the transpose. Apply the original Operator.
Definition: operator.hpp:582
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:504
SecondOrderTimeDependentOperator(int n=0, double t_=0.0, Type type_=EXPLICIT)
Construct a &quot;square&quot; SecondOrderTimeDependentOperator y = f(x,dxdt,t), where x, dxdt and y have the s...
Definition: operator.hpp:456
SecondOrderTimeDependentOperator(int h, int w, double t_=0.0, Type type_=EXPLICIT)
Construct a SecondOrderTimeDependentOperator y = f(x,dxdt,t), where x, dxdt and y have the same dimen...
Definition: operator.hpp:462
virtual MemoryClass GetMemoryClass() const
Returns the type of memory in which the solution and temporaries are stored.
Definition: operator.hpp:686
bool own_A
Ownership flag for A.
Definition: operator.hpp:671
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:84
void FormRectangularConstrainedSystemOperator(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, RectangularConstrainedOperator *&Aout)
see FormRectangularSystemOperator()
Definition: operator.cpp:147
Solver(int h, int w, bool iter_mode=false)
Initialize a Solver with height h and width w.
Definition: operator.hpp:514
double t
Current time.
Definition: operator.hpp:283
virtual int SUNImplicitSolve(const Vector &b, Vector &x, double tol)
Solve the ODE linear system as setup by the method SUNImplicitSetup().
Definition: operator.cpp:259
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.cpp:214
void FormRectangularLinearSystem(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, Vector &x, Vector &b, Operator *&A, Vector &X, Vector &B)
Form a column-constrained linear system using a matrix-free approach.
Definition: operator.cpp:68
virtual void Mult(const Vector &x, Vector &y) const
Operator application. Apply the transpose of the original Operator.
Definition: operator.hpp:578
ID for class SparseMatrix.
Definition: operator.hpp:232
virtual const Operator * GetProlongation() const
Prolongation operator from linear algebra (linear system) vectors, to input vectors for the operator...
Definition: operator.hpp:97
TransposeOperator(const Operator &a)
Construct the transpose of a given operator a.
Definition: operator.hpp:574
virtual void MultTranspose(const Vector &x, Vector &y) const
Application of the transpose.
Definition: operator.hpp:629
General product operator: x -&gt; (A*B)(x) = A(B(x)).
Definition: operator.hpp:588
ID for the base class Operator, i.e. any type.
Definition: operator.hpp:231
virtual void SetEvalMode(const EvalMode new_eval_mode)
Set the evaluation mode of the time-dependent operator.
Definition: operator.hpp:326
virtual MemoryClass GetMemoryClass() const
Return the MemoryClass preferred by the Operator.
Definition: operator.hpp:77
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:57
virtual ~ConstrainedOperator()
Destructor: destroys the unconstrained Operator, if owned.
Definition: operator.hpp:710
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:153
void FormConstrainedSystemOperator(const Array< int > &ess_tdof_list, ConstrainedOperator *&Aout)
see FormSystemOperator()
Definition: operator.cpp:134
double b
Definition: lissajous.cpp:42
This type assumes F(x,k,t) = k, i.e. k = f(x,t) = G(x,t).
Definition: operator.hpp:264
ID for class PetscParMatrix, unspecified format.
Definition: operator.hpp:239
virtual int SUNMassSetup()
Setup the mass matrix in the ODE system .
Definition: operator.cpp:265
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:296
The operator x -&gt; R*A*P*x constructed through the actions of R^T, A and P.
Definition: operator.hpp:608
virtual void Mult(const Vector &x, Vector &y) const
Rectangular-constrained operator action.
Definition: operator.cpp:518
ID for class PetscParMatrix, MATNEST format.
Definition: operator.hpp:237
bool isExplicit() const
True if type is EXPLICIT.
Definition: operator.hpp:306
RAPOperator(const Operator &Rt_, const Operator &A_, const Operator &P_)
Construct the RAP operator given R^T, A and P.
Definition: operator.cpp:327
virtual const Operator * GetRestriction() const
Restriction operator from input vectors for the operator to linear algebra (linear system) vectors...
Definition: operator.hpp:100
This is the most general type, no assumptions on F and G.
Definition: operator.hpp:265
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
Definition: operator.hpp:556
ID for class PetscParMatrix, MATSHELL format.
Definition: operator.hpp:236
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: operator.hpp:650
Operator(int s=0)
Construct a square Operator with given size s (default 0).
Definition: operator.hpp:50
IdentityOperator(int n)
Create an identity operator of size n.
Definition: operator.hpp:527
bool IsIdentityProlongation(const Operator *P)
Definition: operator.hpp:538
Type
Enumeration defining IDs for some classes derived from Operator.
Definition: operator.hpp:229
virtual int SUNMassSolve(const Vector &b, Vector &x, double tol)
Solve the mass matrix linear system as setup by the method SUNMassSetup().
Definition: operator.cpp:271
Scaled Operator B: x -&gt; a A(x).
Definition: operator.hpp:544
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:290
Vector w
Auxiliary vectors.
Definition: operator.hpp:672
int NumRows() const
Get the number of rows (size of output) of the Operator. Synonym with Height().
Definition: operator.hpp:60
void InitTVectors(const Operator *Po, const Operator *Ri, const Operator *Pi, Vector &x, Vector &b, Vector &X, Vector &B) const
Initializes memory for true vectors of linear system.
Definition: operator.cpp:22
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.cpp:236
Base abstract class for second order time dependent operators.
Definition: operator.hpp:451
int NumCols() const
Get the number of columns (size of input) of the Operator. Synonym with Width().
Definition: operator.hpp:66
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:51
ScaledOperator(const Operator *A, double a)
Create an operator which is a scalar multiple of A.
Definition: operator.hpp:552
The transpose of a given operator. Switches the roles of the methods Mult() and MultTranspose().
Definition: operator.hpp:563
Type GetType() const
Return the type ID of the Operator class.
Definition: operator.hpp:250
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:85
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:600
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
This type assumes that G(x,t) = 0.
Definition: operator.hpp:266
double a
Definition: lissajous.cpp:41
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.cpp:219
Rectangular Operator for imposing essential boundary conditions on the input space using only the act...
Definition: operator.hpp:720
Host memory; using new[] and delete[].
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:419
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:188
ID for class PetscParMatrix, MATAIJ format.
Definition: operator.hpp:234
void FormSystemOperator(const Array< int > &ess_tdof_list, Operator *&A)
Return in A a parallel (on truedofs) version of this square operator.
Definition: operator.cpp:164
General triple product operator x -&gt; A*B*C*x, with ownership of the factors.
Definition: operator.hpp:635
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:405
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
Definition: operator.hpp:530
EvalMode GetEvalMode() const
Return the current evaluation mode. See SetEvalMode() for details.
Definition: operator.hpp:313
Operator(int h, int w)
Construct an Operator with the given height (output size) and width (input size). ...
Definition: operator.hpp:54
virtual int SUNImplicitSetup(const Vector &x, const Vector &fx, int jok, int *jcur, double gamma)
Setup the ODE linear system or , where .
Definition: operator.cpp:251
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.cpp:225
TripleProductOperator(const Operator *A, const Operator *B, const Operator *C, bool ownA, bool ownB, bool ownC)
Definition: operator.cpp:361
virtual int SUNMassMult(const Vector &x, Vector &v)
Compute the mass matrix-vector product .
Definition: operator.cpp:277
ID for class ComplexSparseMatrix.
Definition: operator.hpp:241
Operator * SetupRAP(const Operator *Pi, const Operator *Po)
Returns RAP Operator of this, taking in input/output Prolongation matrices.
Definition: operator.cpp:105
virtual void SetOperator(const Operator &op)=0
Set/update the solver for the given operator.
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: operator.hpp:597
EvalMode
Evaluation mode. See SetEvalMode() for details.
Definition: operator.hpp:270
Vector data type.
Definition: vector.hpp:48
Identity Operator I: x -&gt; x.
Definition: operator.hpp:523
ID for class HypreParMatrix.
Definition: operator.hpp:233
TransposeOperator(const Operator *a)
Construct the transpose of a given operator *a.
Definition: operator.hpp:570
Base class for solvers.
Definition: operator.hpp:500
Operator * A
The unconstrained Operator.
Definition: operator.hpp:670
ID for class ComplexOperator.
Definition: operator.hpp:240
Square Operator for imposing essential boundary conditions using only the action, Mult()...
Definition: operator.hpp:666
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
Abstract operator.
Definition: operator.hpp:24
virtual void ImplicitSolve(const double dt0, const double dt1, const Vector &x, const Vector &dxdt, Vector &k)
Solve the equation: k = f(x + 1/2 dt0^2 k, dxdt + dt1 k, t), for the unknown k at the current time t...
Definition: operator.cpp:291
virtual ~Operator()
Virtual destructor.
Definition: operator.hpp:225
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:653
ID for class PetscParMatrix, MATIS format.
Definition: operator.hpp:235
MemoryClass
Memory classes identify sets of memory types.
Definition: mem_manager.hpp:57
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
RectangularConstrainedOperator(Operator *A, const Array< int > &trial_list, const Array< int > &test_list, bool own_A=false)
Constructor from a general Operator and a list of essential indices/dofs.
Definition: operator.cpp:474
virtual const Operator * GetOutputRestriction() const
Restriction operator from output vectors for the operator to linear algebra (linear system) vectors...
Definition: operator.hpp:109
virtual MemoryClass GetMemoryClass() const
Return the MemoryClass preferred by the Operator.
Definition: operator.hpp:622