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