 MFEM  v4.4.0 Finite element discretization library
operator.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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  /** @brief Returns RAP Operator of this, using input/output Prolongation matrices
41  @a Pi corresponds to "P", @a Po corresponds to "Rt" */
42  Operator *SetupRAP(const Operator *Pi, const Operator *Po);
43
44 public:
45  /// Defines operator diagonal policy upon elimination of rows and/or columns.
47  {
48  DIAG_ZERO, ///< Set the diagonal value to zero
49  DIAG_ONE, ///< Set the diagonal value to one
50  DIAG_KEEP ///< Keep the diagonal value
51  };
52
53  /// Initializes memory for true vectors of linear system
54  void InitTVectors(const Operator *Po, const Operator *Ri, const Operator *Pi,
55  Vector &x, Vector &b,
56  Vector &X, Vector &B) const;
57
58  /// Construct a square Operator with given size s (default 0).
59  explicit Operator(int s = 0) { height = width = s; }
60
61  /** @brief Construct an Operator with the given height (output size) and
62  width (input size). */
63  Operator(int h, int w) { height = h; width = w; }
64
65  /// Get the height (size of output) of the Operator. Synonym with NumRows().
66  inline int Height() const { return height; }
67  /** @brief Get the number of rows (size of output) of the Operator. Synonym
68  with Height(). */
69  inline int NumRows() const { return height; }
70
71  /// Get the width (size of input) of the Operator. Synonym with NumCols().
72  inline int Width() const { return width; }
73  /** @brief Get the number of columns (size of input) of the Operator. Synonym
74  with Width(). */
75  inline int NumCols() const { return width; }
76
77  /// Return the MemoryClass preferred by the Operator.
78  /** This is the MemoryClass that will be used to access the input and output
79  vectors in the Mult() and MultTranspose() methods.
80
81  For example, classes using the MFEM_FORALL macro for implementation can
82  return the value returned by Device::GetMemoryClass().
83
84  The default implementation of this method in class Operator returns
85  MemoryClass::HOST. */
86  virtual MemoryClass GetMemoryClass() const { return MemoryClass::HOST; }
87
88  /// Operator application: y=A(x).
89  virtual void Mult(const Vector &x, Vector &y) const = 0;
90
91  /** @brief Action of the transpose operator: y=A^t(x). The default behavior
92  in class Operator is to generate an error. */
93  virtual void MultTranspose(const Vector &x, Vector &y) const
94  { mfem_error("Operator::MultTranspose() is not overloaded!"); }
95
96  /** @brief Evaluate the gradient operator at the point @a x. The default
97  behavior in class Operator is to generate an error. */
98  virtual Operator &GetGradient(const Vector &x) const
99  {
101  return const_cast<Operator &>(*this);
102  }
103
104  /** @brief Computes the diagonal entries into @a diag. Typically, this
105  operation only makes sense for linear Operator%s. In some cases, only an
106  approximation of the diagonal is computed. */
107  virtual void AssembleDiagonal(Vector &diag) const
108  {
109  MFEM_CONTRACT_VAR(diag);
110  MFEM_ABORT("Not relevant or not implemented for this Operator.");
111  }
112
113  /** @brief Prolongation operator from linear algebra (linear system) vectors,
114  to input vectors for the operator. NULL means identity. */
115  virtual const Operator *GetProlongation() const { return NULL; }
116  /** @brief Restriction operator from input vectors for the operator to linear
117  algebra (linear system) vectors. NULL means identity. */
118  virtual const Operator *GetRestriction() const { return NULL; }
119  /** @brief Prolongation operator from linear algebra (linear system) vectors,
120  to output vectors for the operator. NULL means identity. */
121  virtual const Operator *GetOutputProlongation() const
122  {
123  return GetProlongation(); // Assume square unless specialized
124  }
125  /** @brief Transpose of GetOutputRestriction, directly available in this
126  form to facilitate matrix-free RAP-type operators.
127
128  NULL means identity. */
129  virtual const Operator *GetOutputRestrictionTranspose() const { return NULL; }
130  /** @brief Restriction operator from output vectors for the operator to linear
131  algebra (linear system) vectors. NULL means identity. */
132  virtual const Operator *GetOutputRestriction() const
133  {
134  return GetRestriction(); // Assume square unless specialized
135  }
136
137  /** @brief Form a constrained linear system using a matrix-free approach.
138
139  Assuming square operator, form the operator linear system A(X)=B,
140  corresponding to it and the right-hand side @a b, by applying any
141  necessary transformations such as: parallel assembly, conforming
142  constraints for non-conforming AMR and eliminating boundary conditions.
143  @note Static condensation and hybridization are not supported for general
144  operators (cf. the analogous methods BilinearForm::FormLinearSystem() and
145  ParBilinearForm::FormLinearSystem()).
146
147  The constraints are specified through the prolongation P from
148  GetProlongation(), and restriction R from GetRestriction() methods, which
149  are e.g. available through the (parallel) finite element space of any
150  (parallel) bilinear form operator. We assume that the operator is square,
151  using the same input and output space, so we have: A(X)=[P^t (*this)
152  P](X), B=P^t(b), and X=R(x).
153
154  The vector @a x must contain the essential boundary condition values.
155  These are eliminated through the ConstrainedOperator class and the vector
156  @a X is initialized by setting its essential entries to the boundary
157  conditions and all other entries to zero (@a copy_interior == 0) or
158  copied from @a x (@a copy_interior != 0).
159
160  After solving the system A(X)=B, the (finite element) solution @a x can
161  be recovered by calling Operator::RecoverFEMSolution() with the same
162  vectors @a X, @a b, and @a x.
163
164  @note The caller is responsible for destroying the output operator @a A!
165  @note If there are no transformations, @a X simply reuses the data of @a
166  x. */
167  void FormLinearSystem(const Array<int> &ess_tdof_list,
168  Vector &x, Vector &b,
169  Operator* &A, Vector &X, Vector &B,
170  int copy_interior = 0);
171
172  /** @brief Form a column-constrained linear system using a matrix-free approach.
173
174  Form the operator linear system A(X)=B corresponding to the operator
175  and the right-hand side @a b, by applying any necessary transformations
176  such as: parallel assembly, conforming constraints for non-conforming AMR
177  and eliminating boundary conditions. @note Static condensation and
178  hybridization are not supported for general operators (cf. the method
179  MixedBilinearForm::FormRectangularLinearSystem())
180
181  The constraints are specified through the input prolongation Pi from
182  GetProlongation(), and output restriction Ro from GetOutputRestriction()
183  methods, which are e.g. available through the (parallel) finite element
184  spaces of any (parallel) mixed bilinear form operator. So we have:
185  A(X)=[Ro (*this) Pi](X), B=Ro(b), and X=Pi^T(x).
186
187  The vector @a x must contain the essential boundary condition values.
188  The "columns" in this operator corresponding to these values are
189  eliminated through the RectangularConstrainedOperator class.
190
191  After solving the system A(X)=B, the (finite element) solution @a x can
192  be recovered by calling Operator::RecoverFEMSolution() with the same
193  vectors @a X, @a b, and @a x.
194
195  @note The caller is responsible for destroying the output operator @a A!
196  @note If there are no transformations, @a X simply reuses the data of @a
197  x. */
198  void FormRectangularLinearSystem(const Array<int> &trial_tdof_list,
199  const Array<int> &test_tdof_list,
200  Vector &x, Vector &b,
201  Operator* &A, Vector &X, Vector &B);
202
203  /** @brief Reconstruct a solution vector @a x (e.g. a GridFunction) from the
204  solution @a X of a constrained linear system obtained from
205  Operator::FormLinearSystem() or Operator::FormRectangularLinearSystem().
206
207  Call this method after solving a linear system constructed using
208  Operator::FormLinearSystem() to recover the solution as an input vector,
209  @a x, for this Operator (presumably a finite element grid function). This
210  method has identical signature to the analogous method for bilinear
211  forms, though currently @a b is not used in the implementation. */
212  virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x);
213
214  /** @brief Return in @a A a parallel (on truedofs) version of this square
215  operator.
216
217  This returns the same operator as FormLinearSystem(), but does without
218  the transformations of the right-hand side and initial guess. */
219  void FormSystemOperator(const Array<int> &ess_tdof_list,
220  Operator* &A);
221
222  /** @brief Return in @a A a parallel (on truedofs) version of this
223  rectangular operator (including constraints).
224
225  This returns the same operator as FormRectangularLinearSystem(), but does
226  without the transformations of the right-hand side. */
227  void FormRectangularSystemOperator(const Array<int> &trial_tdof_list,
228  const Array<int> &test_tdof_list,
229  Operator* &A);
230
231  /** @brief Return in @a A a parallel (on truedofs) version of this
232  rectangular operator.
233
234  This is similar to FormSystemOperator(), but for dof-to-dof mappings
235  (discrete linear operators), which can also correspond to rectangular
236  matrices. The user should provide specializations of GetProlongation()
237  for the input dofs and GetOutputRestriction() for the output dofs in
238  their Operator implementation that are appropriate for the two spaces the
239  Operator maps between. These are e.g. available through the (parallel)
240  finite element space of any (parallel) bilinear form operator. We have:
241  A(X)=[Rout (*this) Pin](X). */
242  void FormDiscreteOperator(Operator* &A);
243
244  /// Prints operator with input size n and output size m in Matlab format.
245  void PrintMatlab(std::ostream & out, int n, int m = 0) const;
246
247  /// Prints operator in Matlab format.
248  virtual void PrintMatlab(std::ostream & out) const;
249
250  /// Virtual destructor.
251  virtual ~Operator() { }
252
253  /// Enumeration defining IDs for some classes derived from Operator.
254  /** This enumeration is primarily used with class OperatorHandle. */
255  enum Type
256  {
257  ANY_TYPE, ///< ID for the base class Operator, i.e. any type.
258  MFEM_SPARSEMAT, ///< ID for class SparseMatrix.
259  Hypre_ParCSR, ///< ID for class HypreParMatrix.
260  PETSC_MATAIJ, ///< ID for class PetscParMatrix, MATAIJ format.
261  PETSC_MATIS, ///< ID for class PetscParMatrix, MATIS format.
262  PETSC_MATSHELL, ///< ID for class PetscParMatrix, MATSHELL format.
263  PETSC_MATNEST, ///< ID for class PetscParMatrix, MATNEST format.
264  PETSC_MATHYPRE, ///< ID for class PetscParMatrix, MATHYPRE format.
265  PETSC_MATGENERIC, ///< ID for class PetscParMatrix, unspecified format.
266  Complex_Operator, ///< ID for class ComplexOperator.
267  MFEM_ComplexSparseMat, ///< ID for class ComplexSparseMatrix.
268  Complex_Hypre_ParCSR ///< ID for class ComplexHypreParMatrix.
269  };
270
271  /// Return the type ID of the Operator class.
272  /** This method is intentionally non-virtual, so that it returns the ID of
273  the specific pointer or reference type used when calling this method. If
274  not overridden by derived classes, they will automatically use the type ID
275  of the base Operator class, ANY_TYPE. */
276  Type GetType() const { return ANY_TYPE; }
277 };
278
279
280 /// Base abstract class for first order time dependent operators.
281 /** Operator of the form: (x,t) -> f(x,t), where k = f(x,t) generally solves the
282  algebraic equation F(x,k,t) = G(x,t). The functions F and G represent the
283  _implicit_ and _explicit_ parts of the operator, respectively. For explicit
284  operators, F(x,k,t) = k, so f(x,t) = G(x,t). */
286 {
287 public:
288  enum Type
289  {
290  EXPLICIT, ///< This type assumes F(x,k,t) = k, i.e. k = f(x,t) = G(x,t).
291  IMPLICIT, ///< This is the most general type, no assumptions on F and G.
292  HOMOGENEOUS ///< This type assumes that G(x,t) = 0.
293  };
294
295  /// Evaluation mode. See SetEvalMode() for details.
296  enum EvalMode
297  {
298  /** Normal evaluation. */
300  /** Assuming additive split, f(x,t) = f1(x,t) + f2(x,t), evaluate the
301  first term, f1. */
303  /** Assuming additive split, f(x,t) = f1(x,t) + f2(x,t), evaluate the
304  second term, f2. */
306  };
307
308 protected:
309  double t; ///< Current time.
310  Type type; ///< Describes the form of the TimeDependentOperator.
311  EvalMode eval_mode; ///< Current evaluation mode.
312
313 public:
314  /** @brief Construct a "square" TimeDependentOperator y = f(x,t), where x and
315  y have the same dimension @a n. */
316  explicit TimeDependentOperator(int n = 0, double t_ = 0.0,
317  Type type_ = EXPLICIT)
318  : Operator(n) { t = t_; type = type_; eval_mode = NORMAL; }
319
320  /** @brief Construct a TimeDependentOperator y = f(x,t), where x and y have
321  dimensions @a w and @a h, respectively. */
322  TimeDependentOperator(int h, int w, double t_ = 0.0, Type type_ = EXPLICIT)
323  : Operator(h, w) { t = t_; type = type_; eval_mode = NORMAL; }
324
325  /// Read the currently set time.
326  virtual double GetTime() const { return t; }
327
328  /// Set the current time.
329  virtual void SetTime(const double t_) { t = t_; }
330
331  /// True if #type is #EXPLICIT.
332  bool isExplicit() const { return (type == EXPLICIT); }
333  /// True if #type is #IMPLICIT or #HOMOGENEOUS.
334  bool isImplicit() const { return !isExplicit(); }
335  /// True if #type is #HOMOGENEOUS.
336  bool isHomogeneous() const { return (type == HOMOGENEOUS); }
337
338  /// Return the current evaluation mode. See SetEvalMode() for details.
339  EvalMode GetEvalMode() const { return eval_mode; }
340
341  /// Set the evaluation mode of the time-dependent operator.
342  /** The evaluation mode is a switch that allows time-stepping methods to
343  request evaluation of separate components/terms of the time-dependent
344  operator. For example, IMEX methods typically assume additive split of
345  the operator: f(x,t) = f1(x,t) + f2(x,t) and they rely on the ability to
346  evaluate the two terms separately.
347
348  Generally, setting the evaluation mode should affect the behavior of all
349  evaluation-related methods in the class, such as Mult(), ImplicitSolve(),
350  etc. However, the exact list of methods that need to support a specific
351  mode will depend on the used time-stepping method. */
352  virtual void SetEvalMode(const EvalMode new_eval_mode)
353  { eval_mode = new_eval_mode; }
354
355  /** @brief Perform the action of the explicit part of the operator, G:
356  @a y = G(@a x, t) where t is the current time.
357
358  Presently, this method is used by some PETSc ODE solvers, for more
359  details, see the PETSc Manual. */
360  virtual void ExplicitMult(const Vector &x, Vector &y) const;
361
362  /** @brief Perform the action of the implicit part of the operator, F:
363  @a y = F(@a x, @a k, t) where t is the current time.
364
365  Presently, this method is used by some PETSc ODE solvers, for more
366  details, see the PETSc Manual.*/
367  virtual void ImplicitMult(const Vector &x, const Vector &k, Vector &y) const;
368
369  /** @brief Perform the action of the operator: @a y = k = f(@a x, t), where
370  k solves the algebraic equation F(@a x, k, t) = G(@a x, t) and t is the
371  current time. */
372  virtual void Mult(const Vector &x, Vector &y) const;
373
374  /** @brief Solve the equation: @a k = f(@a x + @a dt @a k, t), for the
375  unknown @a k at the current time t.
376
377  For general F and G, the equation for @a k becomes:
378  F(@a x + @a dt @a k, @a k, t) = G(@a x + @a dt @a k, t).
379
380  The input vector @a x corresponds to time index (or cycle) n, while the
381  currently set time, #t, and the result vector @a k correspond to time
382  index n+1. The time step @a dt corresponds to the time interval between
383  cycles n and n+1.
384
385  This method allows for the abstract implementation of some time
386  integration methods, including diagonal implicit Runge-Kutta (DIRK)
387  methods and the backward Euler method in particular.
388
389  If not re-implemented, this method simply generates an error. */
390  virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k);
391
392  /** @brief Return an Operator representing (dF/dk @a shift + dF/dx) at the
393  given @a x, @a k, and the currently set time.
394
395  Presently, this method is used by some PETSc ODE solvers, for more
396  details, see the PETSc Manual. */
397  virtual Operator& GetImplicitGradient(const Vector &x, const Vector &k,
398  double shift) const;
399
400  /** @brief Return an Operator representing dG/dx at the given point @a x and
401  the currently set time.
402
403  Presently, this method is used by some PETSc ODE solvers, for more
404  details, see the PETSc Manual. */
405  virtual Operator& GetExplicitGradient(const Vector &x) const;
406
407  /** @brief Setup the ODE linear system \f$A(x,t) = (I - gamma J) \f$ or
408  \f$A = (M - gamma J) \f$, where \f$J(x,t) = \frac{df}{dt(x,t)} \f$.
409
410  @param[in] x The state at which \f$A(x,t)\f$ should be evaluated.
411  @param[in] fx The current value of the ODE rhs function, \f$f(x,t)\f$.
412  @param[in] jok Flag indicating if the Jacobian should be updated.
413  @param[out] jcur Flag to signal if the Jacobian was updated.
414  @param[in] gamma The scaled time step value.
415
416  If not re-implemented, this method simply generates an error.
417
418  Presently, this method is used by SUNDIALS ODE solvers, for more
419  details, see the SUNDIALS User Guides. */
420  virtual int SUNImplicitSetup(const Vector &x, const Vector &fx,
421  int jok, int *jcur, double gamma);
422
423  /** @brief Solve the ODE linear system \f$A x = b \f$ as setup by
424  the method SUNImplicitSetup().
425
426  @param[in] b The linear system right-hand side.
427  @param[in,out] x On input, the initial guess. On output, the solution.
428  @param[in] tol Linear solve tolerance.
429
430  If not re-implemented, this method simply generates an error.
431
432  Presently, this method is used by SUNDIALS ODE solvers, for more
433  details, see the SUNDIALS User Guides. */
434  virtual int SUNImplicitSolve(const Vector &b, Vector &x, double tol);
435
436  /** @brief Setup the mass matrix in the ODE system \f$M y' = f(y,t) \f$ .
437
438  If not re-implemented, this method simply generates an error.
439
440  Presently, this method is used by SUNDIALS ARKStep integrator, for more
441  details, see the ARKode User Guide. */
442  virtual int SUNMassSetup();
443
444  /** @brief Solve the mass matrix linear system \f$M x = b \f$
445  as setup by the method SUNMassSetup().
446
447  @param[in] b The linear system right-hand side.
448  @param[in,out] x On input, the initial guess. On output, the solution.
449  @param[in] tol Linear solve tolerance.
450
451  If not re-implemented, this method simply generates an error.
452
453  Presently, this method is used by SUNDIALS ARKStep integrator, for more
454  details, see the ARKode User Guide. */
455  virtual int SUNMassSolve(const Vector &b, Vector &x, double tol);
456
457  /** @brief Compute the mass matrix-vector product \f$v = M x \f$ .
458
459  @param[in] x The vector to multiply.
460  @param[out] v The result of the matrix-vector product.
461
462  If not re-implemented, this method simply generates an error.
463
464  Presently, this method is used by SUNDIALS ARKStep integrator, for more
465  details, see the ARKode User Guide. */
466  virtual int SUNMassMult(const Vector &x, Vector &v);
467
469 };
470
471
473  equations to be used with CVODESSolver. */
475 {
476 public:
477
478  /**
479  \brief The TimedependentAdjointOperator extends the TimeDependentOperator
480  class to use features in SUNDIALS CVODESSolver for computing quadratures
482
484  method to tell CVODES what the adjoint rate equation is.
485
486  QuadratureIntegration (optional) can be used to compute values over the
487  forward problem
488
489  QuadratureSensitivityMult (optional) can be used to find the sensitivity
491
492  SUNImplicitSetupB (optional) can be used to setup custom solvers for the
493  newton solve for the adjoint problem.
494
495  SUNImplicitSolveB (optional) actually uses the solvers from
496  SUNImplicitSetupB to solve the adjoint problem.
497
498  See SUNDIALS user manuals for specifics.
499
500  \param[in] dim Dimension of the forward operator
501  \param[in] adjdim Dimension of the adjoint operator. Typically it is the
502  same size as dim. However, SUNDIALS allows users to specify the size if
503  one wants to perform custom operations.
504  \param[in] t Starting time to set
505  \param[in] type The TimeDependentOperator type
506  */
508  Type type = EXPLICIT) :
511  {}
512
513  /// Destructor
515
516  /**
517  \brief Provide the operator integration of a quadrature equation
518
519  \param[in] y The current value at time t
520  \param[out] qdot The current quadrature rate value at t
521  */
522  virtual void QuadratureIntegration(const Vector &y, Vector &qdot) const {};
523
524  /** @brief Perform the action of the operator:
525  @a yBdot = k = f(@a y,@2 yB, t), where
526
527  @param[in] y The primal solution at time t
528  @param[in] yB The adjoint solution at time t
529  @param[out] yBdot the rate at time t
530  */
531  virtual void AdjointRateMult(const Vector &y, Vector & yB,
532  Vector &yBdot) const = 0;
533
534  /**
535  \brief Provides the sensitivity of the quadrature w.r.t to primal and
537
538  \param[in] y the value of the primal solution at time t
539  \param[in] yB the value of the adjoint solution at time t
540  \param[out] qBdot the value of the sensitivity of the quadrature rate at
541  time t
542  */
543  virtual void QuadratureSensitivityMult(const Vector &y, const Vector &yB,
544  Vector &qBdot) const {}
545
546  /** @brief Setup the ODE linear system \f$A(x,t) = (I - gamma J) \f$ or
547  \f$A = (M - gamma J) \f$, where \f$J(x,t) = \frac{df}{dt(x,t)} \f$.
548
549  @param[in] t The current time
550  @param[in] x The state at which \f$A(x,xB,t)\f$ should be evaluated.
551  @param[in] xB The state at which \f$A(x,xB,t)\f$ should be evaluated.
552  @param[in] fxB The current value of the ODE rhs function, \f$f(x,t)\f$.
553  @param[in] jokB Flag indicating if the Jacobian should be updated.
554  @param[out] jcurB Flag to signal if the Jacobian was updated.
555  @param[in] gammaB The scaled time step value.
556
557  If not re-implemented, this method simply generates an error.
558
559  Presently, this method is used by SUNDIALS ODE solvers, for more details,
560  see the SUNDIALS User Guides.
561  */
562  virtual int SUNImplicitSetupB(const double t, const Vector &x,
563  const Vector &xB, const Vector &fxB,
564  int jokB, int *jcurB, double gammaB)
565  {
567  "overridden!");
568  return (-1);
569  }
570
571  /** @brief Solve the ODE linear system \f$A(x,xB,t) xB = b \f$ as setup by
572  the method SUNImplicitSetup().
573
574  @param[in] b The linear system right-hand side.
575  @param[in,out] x On input, the initial guess. On output, the solution.
576  @param[in] tol Linear solve tolerance.
577
578  If not re-implemented, this method simply generates an error.
579
580  Presently, this method is used by SUNDIALS ODE solvers, for more details,
581  see the SUNDIALS User Guides. */
582  virtual int SUNImplicitSolveB(Vector &x, const Vector &b, double tol)
583  {
585  "overridden!");
586  return (-1);
587  }
588
589  /// Returns the size of the adjoint problem state space
591
592 protected:
594 };
595
596
597 /// Base abstract class for second order time dependent operators.
598 /** Operator of the form: (x,dxdt,t) -> f(x,dxdt,t), where k = f(x,dxdt,t)
599  generally solves the algebraic equation F(x,dxdt,k,t) = G(x,dxdt,t).
600  The functions F and G represent the_implicit_ and _explicit_ parts of
601  the operator, respectively. For explicit operators,
602  F(x,dxdt,k,t) = k, so f(x,dxdt,t) = G(x,dxdt,t). */
604 {
605 public:
606  /** @brief Construct a "square" SecondOrderTimeDependentOperator
607  y = f(x,dxdt,t), where x, dxdt and y have the same dimension @a n. */
608  explicit SecondOrderTimeDependentOperator(int n = 0, double t_ = 0.0,
609  Type type_ = EXPLICIT)
610  : TimeDependentOperator(n, t_,type_) { }
611
612  /** @brief Construct a SecondOrderTimeDependentOperator y = f(x,dxdt,t),
613  where x, dxdt and y have the same dimension @a n. */
614  SecondOrderTimeDependentOperator(int h, int w, double t_ = 0.0,
615  Type type_ = EXPLICIT)
616  : TimeDependentOperator(h, w, t_,type_) { }
617
619
620  /** @brief Perform the action of the operator: @a y = k = f(@a x,@ dxdt, t),
621  where k solves the algebraic equation
622  F(@a x,@ dxdt, k, t) = G(@a x,@ dxdt, t) and t is the current time. */
623  virtual void Mult(const Vector &x, const Vector &dxdt, Vector &y) const;
624
626  /** @brief Solve the equation:
627  @a k = f(@a x + @a fac0 @a k, @a dxdt + @a fac1 @a k, t), for the
628  unknown @a k at the current time t.
629
630  For general F and G, the equation for @a k becomes:
631  F(@a x + @a fac0 @a k, @a dxdt + @a fac1 @a k, t)
632  = G(@a x + @a fac0 @a k, @a dxdt + @a fac1 @a k, t).
633
634  The input vectors @a x and @a dxdt corresponds to time index (or cycle) n, while the
635  currently set time, #t, and the result vector @a k correspond to time
636  index n+1.
637
638  This method allows for the abstract implementation of some time
639  integration methods.
640
641  If not re-implemented, this method simply generates an error. */
642  virtual void ImplicitSolve(const double fac0, const double fac1,
643  const Vector &x, const Vector &dxdt, Vector &k);
644
645
647 };
648
649
650 /// Base class for solvers
651 class Solver : public Operator
652 {
653 public:
654  /// If true, use the second argument of Mult() as an initial guess.
656
657  /** @brief Initialize a square Solver with size @a s.
658
659  @warning Use a Boolean expression for the second parameter (not an int)
660  to distinguish this call from the general rectangular constructor. */
661  explicit Solver(int s = 0, bool iter_mode = false)
662  : Operator(s) { iterative_mode = iter_mode; }
663
664  /// Initialize a Solver with height @a h and width @a w.
665  Solver(int h, int w, bool iter_mode = false)
666  : Operator(h, w) { iterative_mode = iter_mode; }
667
668  /// Set/update the solver for the given operator.
669  virtual void SetOperator(const Operator &op) = 0;
670 };
671
672
673 /// Identity Operator I: x -> x.
675 {
676 public:
677  /// Create an identity operator of size @a n.
678  explicit IdentityOperator(int n) : Operator(n) { }
679
680  /// Operator application
681  virtual void Mult(const Vector &x, Vector &y) const { y = x; }
682
683  /// Application of the transpose
684  virtual void MultTranspose(const Vector &x, Vector &y) const { y = x; }
685 };
686
687 /// Returns true if P is the identity prolongation, i.e. if it is either NULL or
688 /// an IdentityOperator.
689 inline bool IsIdentityProlongation(const Operator *P)
690 {
691  return !P || dynamic_cast<const IdentityOperator*>(P);
692 }
693
694 /// Scaled Operator B: x -> a A(x).
695 class ScaledOperator : public Operator
696 {
697 private:
698  const Operator &A_;
699  double a_;
700
701 public:
702  /// Create an operator which is a scalar multiple of A.
703  explicit ScaledOperator(const Operator *A, double a)
704  : Operator(A->Height(), A->Width()), A_(*A), a_(a) { }
705
706  /// Operator application
707  virtual void Mult(const Vector &x, Vector &y) const
708  { A_.Mult(x, y); y *= a_; }
709
710  /// Application of the transpose.
711  virtual void MultTranspose(const Vector &x, Vector &y) const
712  { A_.MultTranspose(x, y); y *= a_; }
713 };
714
715
716 /** @brief The transpose of a given operator. Switches the roles of the methods
717  Mult() and MultTranspose(). */
719 {
720 private:
721  const Operator &A;
722
723 public:
724  /// Construct the transpose of a given operator @a *a.
726  : Operator(a->Width(), a->Height()), A(*a) { }
727
728  /// Construct the transpose of a given operator @a a.
730  : Operator(a.Width(), a.Height()), A(a) { }
731
732  /// Operator application. Apply the transpose of the original Operator.
733  virtual void Mult(const Vector &x, Vector &y) const
734  { A.MultTranspose(x, y); }
735
736  /// Application of the transpose. Apply the original Operator.
737  virtual void MultTranspose(const Vector &x, Vector &y) const
738  { A.Mult(x, y); }
739 };
740
741
742 /// General product operator: x -> (A*B)(x) = A(B(x)).
743 class ProductOperator : public Operator
744 {
745  const Operator *A, *B;
746  bool ownA, ownB;
747  mutable Vector z;
748
749 public:
750  ProductOperator(const Operator *A, const Operator *B, bool ownA, bool ownB);
751
752  virtual void Mult(const Vector &x, Vector &y) const
753  { B->Mult(x, z); A->Mult(z, y); }
754
755  virtual void MultTranspose(const Vector &x, Vector &y) const
756  { A->MultTranspose(x, z); B->MultTranspose(z, y); }
757
758  virtual ~ProductOperator();
759 };
760
761
762 /// The operator x -> R*A*P*x constructed through the actions of R^T, A and P
763 class RAPOperator : public Operator
764 {
765 private:
766  const Operator & Rt;
767  const Operator & A;
768  const Operator & P;
769  mutable Vector Px;
770  mutable Vector APx;
771  MemoryClass mem_class;
772
773 public:
774  /// Construct the RAP operator given R^T, A and P.
775  RAPOperator(const Operator &Rt_, const Operator &A_, const Operator &P_);
776
777  virtual MemoryClass GetMemoryClass() const { return mem_class; }
778
779  /// Operator application.
780  virtual void Mult(const Vector & x, Vector & y) const
781  { P.Mult(x, Px); A.Mult(Px, APx); Rt.MultTranspose(APx, y); }
782
783  /// Approximate diagonal of the RAP Operator.
784  /** Returns the diagonal of A, as returned by its AssembleDiagonal method,
785  multiplied be P^T.
786
787  When P is the FE space prolongation operator on a mesh without hanging
788  nodes and Rt = P, the returned diagonal is exact, as long as the diagonal
789  of A is also exact. */
790  virtual void AssembleDiagonal(Vector &diag) const
791  {
792  A.AssembleDiagonal(APx);
793  P.MultTranspose(APx, diag);
794
795  // TODO: For an AMR mesh, a convergent diagonal can be assembled with
796  // |P^T| APx, where |P^T| has entry-wise absolute values of the conforming
797  // prolongation transpose operator. See BilinearForm::AssembleDiagonal.
798  }
799
800  /// Application of the transpose.
801  virtual void MultTranspose(const Vector & x, Vector & y) const
802  { Rt.Mult(x, APx); A.MultTranspose(APx, Px); P.MultTranspose(Px, y); }
803 };
804
805
806 /// General triple product operator x -> A*B*C*x, with ownership of the factors.
808 {
809  const Operator *A;
810  const Operator *B;
811  const Operator *C;
812  bool ownA, ownB, ownC;
813  mutable Vector t1, t2;
814  MemoryClass mem_class;
815
816 public:
817  TripleProductOperator(const Operator *A, const Operator *B,
818  const Operator *C, bool ownA, bool ownB, bool ownC);
819
820  virtual MemoryClass GetMemoryClass() const { return mem_class; }
821
822  virtual void Mult(const Vector &x, Vector &y) const
823  { C->Mult(x, t1); B->Mult(t1, t2); A->Mult(t2, y); }
824
825  virtual void MultTranspose(const Vector &x, Vector &y) const
826  { A->MultTranspose(x, t2); B->MultTranspose(t2, t1); C->MultTranspose(t1, y); }
827
828  virtual ~TripleProductOperator();
829 };
830
831
832 /** @brief Square Operator for imposing essential boundary conditions using only
833  the action, Mult(), of a given unconstrained Operator.
834
835  Square operator constrained by fixing certain entries in the solution to
836  given "essential boundary condition" values. This class is used by the
837  general, matrix-free system formulation of Operator::FormLinearSystem.
838
839  Do not confuse with ConstrainedSolver, which despite the name has very
840  different functionality. */
842 {
843 protected:
844  Array<int> constraint_list; ///< List of constrained indices/dofs.
845  Operator *A; ///< The unconstrained Operator.
846  bool own_A; ///< Ownership flag for A.
847  mutable Vector z, w; ///< Auxiliary vectors.
849  DiagonalPolicy diag_policy; ///< Diagonal policy for constrained dofs
850
851 public:
852  /** @brief Constructor from a general Operator and a list of essential
853  indices/dofs.
854
855  Specify the unconstrained operator @a *A and a @a list of indices to
856  constrain, i.e. each entry @a list[i] represents an essential dof. If the
857  ownership flag @a own_A is true, the operator @a *A will be destroyed
858  when this object is destroyed. The @a diag_policy determines how the
859  operator sets entries corresponding to essential dofs. */
860  ConstrainedOperator(Operator *A, const Array<int> &list, bool own_A = false,
862
863  /// Returns the type of memory in which the solution and temporaries are stored.
864  virtual MemoryClass GetMemoryClass() const { return mem_class; }
865
866  /// Set the diagonal policy for the constrained operator.
867  void SetDiagonalPolicy(const DiagonalPolicy diag_policy_)
868  { diag_policy = diag_policy_; }
869
870  /// Diagonal of A, modified according to the used DiagonalPolicy.
871  virtual void AssembleDiagonal(Vector &diag) const;
872
873  /** @brief Eliminate "essential boundary condition" values specified in @a x
874  from the given right-hand side @a b.
875
876  Performs the following steps:
877
878  z = A((0,x_b)); b_i -= z_i; b_b = x_b;
879
880  where the "_b" subscripts denote the essential (boundary) indices/dofs of
881  the vectors, and "_i" -- the rest of the entries. */
882  void EliminateRHS(const Vector &x, Vector &b) const;
883
884  /** @brief Constrained operator action.
885
886  Performs the following steps:
887
888  z = A((x_i,0)); y_i = z_i; y_b = x_b;
889
890  where the "_b" subscripts denote the essential (boundary) indices/dofs of
891  the vectors, and "_i" -- the rest of the entries. */
892  virtual void Mult(const Vector &x, Vector &y) const;
893
894  /// Destructor: destroys the unconstrained Operator, if owned.
895  virtual ~ConstrainedOperator() { if (own_A) { delete A; } }
896 };
897
898 /** @brief Rectangular Operator for imposing essential boundary conditions on
899  the input space using only the action, Mult(), of a given unconstrained
900  Operator.
901
902  Rectangular operator constrained by fixing certain entries in the solution
903  to given "essential boundary condition" values. This class is used by the
904  general matrix-free formulation of Operator::FormRectangularLinearSystem. */
906 {
907 protected:
910  bool own_A;
911  mutable Vector z, w;
913
914 public:
915  /** @brief Constructor from a general Operator and a list of essential
916  indices/dofs.
917
918  Specify the unconstrained operator @a *A and two lists of indices to
919  constrain, i.e. each entry @a trial_list[i] represents an essential trial
920  dof. If the ownership flag @a own_A is true, the operator @a *A will be
921  destroyed when this object is destroyed. */
923  const Array<int> &test_list, bool own_A = false);
924  /// Returns the type of memory in which the solution and temporaries are stored.
925  virtual MemoryClass GetMemoryClass() const { return mem_class; }
926  /** @brief Eliminate columns corresponding to "essential boundary condition"
927  values specified in @a x from the given right-hand side @a b.
928
929  Performs the following steps:
930
931  b -= A((0,x_b));
932  b_j = 0
933
934  where the "_b" subscripts denote the essential (boundary) indices and the
935  "_j" subscript denotes the essential test indices */
936  void EliminateRHS(const Vector &x, Vector &b) const;
937  /** @brief Rectangular-constrained operator action.
938
939  Performs the following steps:
940
941  y = A((x_i,0));
942  y_j = 0
943
944  where the "_i" subscripts denote all the nonessential (boundary) trial
945  indices and the "_j" subscript denotes the essential test indices */
946  virtual void Mult(const Vector &x, Vector &y) const;
947  virtual void MultTranspose(const Vector &x, Vector &y) const;
948  virtual ~RectangularConstrainedOperator() { if (own_A) { delete A; } }
949 };
950
951 /** @brief PowerMethod helper class to estimate the largest eigenvalue of an
952  operator using the iterative power method. */
954 {
955  Vector v1;
956 #ifdef MFEM_USE_MPI
957  MPI_Comm comm;
958 #endif
959
960 public:
961
962 #ifdef MFEM_USE_MPI
963  PowerMethod() : comm(MPI_COMM_NULL) {}
964 #else
966 #endif
967
968 #ifdef MFEM_USE_MPI
969  PowerMethod(MPI_Comm comm_) : comm(comm_) {}
970 #endif
971
972  /// @brief Returns an estimate of the largest eigenvalue of the operator \p opr
973  /// using the iterative power method.
974  /** \p v0 is being used as the vector for the iterative process and will contain
975  the eigenvector corresponding to the largest eigenvalue after convergence.
976  The maximum number of iterations may set with \p numSteps, the relative
977  tolerance with \p tolerance and the seed of the random initialization of
978  \p v0 with \p seed. If \p seed is 0 \p v0 will not be random-initialized. */
979  double EstimateLargestEigenvalue(Operator& opr, Vector& v0,
980  int numSteps = 10, double tolerance = 1e-8,
981  int seed = 12345);
982 };
983
984 }
985
986 #endif
PowerMethod helper class to estimate the largest eigenvalue of an operator using the iterative power ...
Definition: operator.hpp:953
virtual double GetTime() const
Definition: operator.hpp:326
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:289
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
ConstrainedOperator(Operator *A, const Array< int > &list, bool own_A=false, DiagonalPolicy diag_policy=DIAG_ONE)
Constructor from a general Operator and a list of essential indices/dofs.
Definition: operator.cpp:410
bool isImplicit() const
True if type is IMPLICIT or HOMOGENEOUS.
Definition: operator.hpp:334
bool isHomogeneous() const
True if type is HOMOGENEOUS.
Definition: operator.hpp:336
virtual MemoryClass GetMemoryClass() const
Returns the type of memory in which the solution and temporaries are stored.
Definition: operator.hpp:925
virtual int SUNImplicitSolveB(Vector &x, const Vector &b, double tol)
Solve the ODE linear system as setup by the method SUNImplicitSetup().
Definition: operator.hpp:582
Solver(int s=0, bool iter_mode=false)
Initialize a square Solver with size s.
Definition: operator.hpp:661
virtual ~ProductOperator()
Definition: operator.cpp:325
virtual void AssembleDiagonal(Vector &diag) const
Diagonal of A, modified according to the used DiagonalPolicy.
Definition: operator.cpp:426
virtual void QuadratureSensitivityMult(const Vector &y, const Vector &yB, Vector &qBdot) const
Provides the sensitivity of the quadrature w.r.t to primal and adjoint solutions. ...
Definition: operator.hpp:543
Array< int > constraint_list
List of constrained indices/dofs.
Definition: operator.hpp:844
virtual void ImplicitSolve(const double fac0, const double fac1, const Vector &x, const Vector &dxdt, Vector &k)
Solve the equation: k = f(x + fac0 k, dxdt + fac1 k, t), for the unknown k at the current time t...
Definition: operator.cpp:296
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
Definition: operator.hpp:780
Type type
Describes the form of the TimeDependentOperator.
Definition: operator.hpp:310
Base abstract class for first order time dependent operators.
Definition: operator.hpp:285
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:235
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:98
virtual void AdjointRateMult(const Vector &y, Vector &yB, Vector &yBdot) const =0
Perform the action of the operator: yBdot = k = f(y,@2 yB, t), where.
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
ID for class ComplexHypreParMatrix.
Definition: operator.hpp:268
virtual void Mult(const Vector &x, Vector &y) const
Constrained operator action.
Definition: operator.cpp:484
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:249
virtual MemoryClass GetMemoryClass() const
Return the MemoryClass preferred by the Operator.
Definition: operator.hpp:820
virtual const Operator * GetOutputProlongation() const
Prolongation operator from linear algebra (linear system) vectors, to output vectors for the operator...
Definition: operator.hpp:121
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:550
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:684
EvalMode eval_mode
Current evaluation mode.
Definition: operator.hpp:311
ProductOperator(const Operator *A, const Operator *B, bool ownA, bool ownB)
Definition: operator.cpp:306
ID for class PetscParMatrix, MATHYPRE format.
Definition: operator.hpp:264
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
virtual void SetTime(const double t_)
Set the current time.
Definition: operator.hpp:329
virtual void MultTranspose(const Vector &x, Vector &y) const
Application of the transpose. Apply the original Operator.
Definition: operator.hpp:737
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:655
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:608
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:614
virtual MemoryClass GetMemoryClass() const
Returns the type of memory in which the solution and temporaries are stored.
Definition: operator.hpp:864
Destructor.
Definition: operator.hpp:514
bool own_A
Ownership flag for A.
Definition: operator.hpp:846
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:93
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:665
double t
Current time.
Definition: operator.hpp:309
virtual void AssembleDiagonal(Vector &diag) const
Approximate diagonal of the RAP Operator.
Definition: operator.hpp:790
virtual int SUNImplicitSolve(const Vector &b, Vector &x, double tol)
Solve the ODE linear system as setup by the method SUNImplicitSetup().
Definition: operator.cpp:264
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:219
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:733
ID for class SparseMatrix.
Definition: operator.hpp:258
virtual const Operator * GetProlongation() const
Prolongation operator from linear algebra (linear system) vectors, to input vectors for the operator...
Definition: operator.hpp:115
TransposeOperator(const Operator &a)
Construct the transpose of a given operator a.
Definition: operator.hpp:729
virtual void MultTranspose(const Vector &x, Vector &y) const
Application of the transpose.
Definition: operator.hpp:801
General product operator: x -&gt; (A*B)(x) = A(B(x)).
Definition: operator.hpp:743
ID for the base class Operator, i.e. any type.
Definition: operator.hpp:257
void SetDiagonalPolicy(const DiagonalPolicy diag_policy_)
Set the diagonal policy for the constrained operator.
Definition: operator.hpp:867
virtual void SetEvalMode(const EvalMode new_eval_mode)
Set the evaluation mode of the time-dependent operator.
Definition: operator.hpp:352
virtual MemoryClass GetMemoryClass() const
Return the MemoryClass preferred by the Operator.
Definition: operator.hpp:86
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
virtual ~ConstrainedOperator()
Destructor: destroys the unconstrained Operator, if owned.
Definition: operator.hpp:895
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:154
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:290
ID for class PetscParMatrix, unspecified format.
Definition: operator.hpp:265
virtual int SUNMassSetup()
Setup the mass matrix in the ODE system .
Definition: operator.cpp:270
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:322
Returns the size of the adjoint problem state space.
Definition: operator.hpp:590
The operator x -&gt; R*A*P*x constructed through the actions of R^T, A and P.
Definition: operator.hpp:763
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.cpp:603
virtual void Mult(const Vector &x, Vector &y) const
Rectangular-constrained operator action.
Definition: operator.cpp:575
PowerMethod(MPI_Comm comm_)
Definition: operator.hpp:969
ID for class PetscParMatrix, MATNEST format.
Definition: operator.hpp:263
bool isExplicit() const
True if type is EXPLICIT.
Definition: operator.hpp:332
Set the diagonal value to one.
Definition: operator.hpp:49
RAPOperator(const Operator &Rt_, const Operator &A_, const Operator &P_)
Construct the RAP operator given R^T, A and P.
Definition: operator.cpp:332
virtual const Operator * GetRestriction() const
Restriction operator from input vectors for the operator to linear algebra (linear system) vectors...
Definition: operator.hpp:118
This is the most general type, no assumptions on F and G.
Definition: operator.hpp:291
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
Definition: operator.hpp:707
virtual const Operator * GetOutputRestrictionTranspose() const
Transpose of GetOutputRestriction, directly available in this form to facilitate matrix-free RAP-type...
Definition: operator.hpp:129
ID for class PetscParMatrix, MATSHELL format.
Definition: operator.hpp:262
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: operator.hpp:822
Operator(int s=0)
Construct a square Operator with given size s (default 0).
Definition: operator.hpp:59
IdentityOperator(int n)
Create an identity operator of size n.
Definition: operator.hpp:678
bool IsIdentityProlongation(const Operator *P)
Definition: operator.hpp:689
Type
Enumeration defining IDs for some classes derived from Operator.
Definition: operator.hpp:255
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:276
Scaled Operator B: x -&gt; a A(x).
Definition: operator.hpp:695
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:316
Vector w
Auxiliary vectors.
Definition: operator.hpp:847
int NumRows() const
Get the number of rows (size of output) of the Operator. Synonym with Height().
Definition: operator.hpp:69
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:241
Base abstract class for second order time dependent operators.
Definition: operator.hpp:603
int NumCols() const
Get the number of columns (size of input) of the Operator. Synonym with Width().
Definition: operator.hpp:75
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:703
The transpose of a given operator. Switches the roles of the methods Mult() and MultTranspose().
Definition: operator.hpp:718
Type GetType() const
Return the type ID of the Operator class.
Definition: operator.hpp:276
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:755
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:292
The TimedependentAdjointOperator extends the TimeDependentOperator class to use features in SUNDIALS ...
Definition: operator.hpp:507
virtual void QuadratureIntegration(const Vector &y, Vector &qdot) const
Provide the operator integration of a quadrature equation.
Definition: operator.hpp:522
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:224
Rectangular Operator for imposing essential boundary conditions on the input space using only the act...
Definition: operator.hpp:905
Host memory; using new[] and delete[].
virtual void AssembleDiagonal(Vector &diag) const
Computes the diagonal entries into diag. Typically, this operation only makes sense for linear Operat...
Definition: operator.hpp:107
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:457
Keep the diagonal value.
Definition: operator.hpp:50
ID for class PetscParMatrix, MATAIJ format.
Definition: operator.hpp:260
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
DiagonalPolicy
Defines operator diagonal policy upon elimination of rows and/or columns.
Definition: operator.hpp:46
General triple product operator x -&gt; A*B*C*x, with ownership of the factors.
Definition: operator.hpp:807
int dim
Definition: ex24.cpp:53
virtual void Mult(const Vector &x, Vector &y) const
Operator application.
Definition: operator.hpp:681
EvalMode GetEvalMode() const
Return the current evaluation mode. See SetEvalMode() for details.
Definition: operator.hpp:339
Operator(int h, int w)
Construct an Operator with the given height (output size) and width (input size). ...
Definition: operator.hpp:63
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:256
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:230
TripleProductOperator(const Operator *A, const Operator *B, const Operator *C, bool ownA, bool ownB, bool ownC)
Definition: operator.cpp:366
virtual int SUNMassMult(const Vector &x, Vector &v)
Compute the mass matrix-vector product .
Definition: operator.cpp:282
void PrintMatlab(std::ostream &out, int n, int m=0) const
Prints operator with input size n and output size m in Matlab format.
Definition: operator.cpp:188
ID for class ComplexSparseMatrix.
Definition: operator.hpp:267
Operator * SetupRAP(const Operator *Pi, const Operator *Po)
Returns RAP Operator of this, using input/output Prolongation matrices Pi corresponds to &quot;P&quot;...
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:752
EvalMode
Evaluation mode. See SetEvalMode() for details.
Definition: operator.hpp:296
Vector data type.
Definition: vector.hpp:60
Identity Operator I: x -&gt; x.
Definition: operator.hpp:674
ID for class HypreParMatrix.
Definition: operator.hpp:259
virtual void MultTranspose(const Vector &x, Vector &y) const
Application of the transpose.
Definition: operator.hpp:711
TransposeOperator(const Operator *a)
Construct the transpose of a given operator *a.
Definition: operator.hpp:725
DiagonalPolicy diag_policy
Diagonal policy for constrained dofs.
Definition: operator.hpp:849
RefCoord s
Base class for solvers.
Definition: operator.hpp:651
Operator * A
The unconstrained Operator.
Definition: operator.hpp:845
ID for class ComplexOperator.
Definition: operator.hpp:266
Square Operator for imposing essential boundary conditions using only the action, Mult()...
Definition: operator.hpp:841
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
double EstimateLargestEigenvalue(Operator &opr, Vector &v0, int numSteps=10, double tolerance=1e-8, int seed=12345)
Returns an estimate of the largest eigenvalue of the operator opr using the iterative power method...
Definition: operator.cpp:632
Abstract operator.
Definition: operator.hpp:24
virtual ~Operator()
Virtual destructor.
Definition: operator.hpp:251
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:825
ID for class PetscParMatrix, MATIS format.
Definition: operator.hpp:261
MemoryClass
Memory classes identify sets of memory types.
Definition: mem_manager.hpp:73
Set the diagonal value to zero.
Definition: operator.hpp:48
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
virtual int SUNImplicitSetupB(const double t, const Vector &x, const Vector &xB, const Vector &fxB, int jokB, int *jcurB, double gammaB)
Setup the ODE linear system or , where .
Definition: operator.hpp:562
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:531
virtual const Operator * GetOutputRestriction() const
Restriction operator from output vectors for the operator to linear algebra (linear system) vectors...
Definition: operator.hpp:132
virtual MemoryClass GetMemoryClass() const
Return the MemoryClass preferred by the Operator.
Definition: operator.hpp:777