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