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