MFEM  v4.5.2
Finite element discretization library
bilinearform_ext.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_BILINEARFORM_EXT
13 #define MFEM_BILINEARFORM_EXT
14 
15 #include "../config/config.hpp"
16 #include "fespace.hpp"
17 #include "../general/device.hpp"
18 
19 namespace mfem
20 {
21 
22 class BilinearForm;
23 class MixedBilinearForm;
24 class DiscreteLinearOperator;
25 
26 /// Class extending the BilinearForm class to support different AssemblyLevels.
27 /** FA - Full Assembly
28  PA - Partial Assembly
29  EA - Element Assembly
30  MF - Matrix Free
31 */
33 {
34 protected:
35  BilinearForm *a; ///< Not owned
36 
37 public:
39 
40  virtual MemoryClass GetMemoryClass() const
41  { return Device::GetDeviceMemoryClass(); }
42 
43  /// Get the finite element space prolongation matrix
44  virtual const Operator *GetProlongation() const;
45 
46  /// Get the finite element space restriction matrix
47  virtual const Operator *GetRestriction() const;
48 
49  /// Assemble at the level given for the BilinearFormExtension subclass
50  virtual void Assemble() = 0;
51 
52  virtual void AssembleDiagonal(Vector &diag) const
53  {
54  MFEM_ABORT("AssembleDiagonal not implemented for this assembly level!");
55  }
56 
57  virtual void FormSystemMatrix(const Array<int> &ess_tdof_list,
58  OperatorHandle &A) = 0;
59  virtual void FormLinearSystem(const Array<int> &ess_tdof_list,
60  Vector &x, Vector &b,
61  OperatorHandle &A, Vector &X, Vector &B,
62  int copy_interior = 0) = 0;
63  virtual void Update() = 0;
64 };
65 
66 /// Data and methods for partially-assembled bilinear forms
68 {
69 protected:
70  const FiniteElementSpace *trial_fes, *test_fes; // Not owned
71  mutable Vector localX, localY;
74  const Operator *elem_restrict; // Not owned
77 
78 public:
80 
81  void Assemble();
82  void AssembleDiagonal(Vector &diag) const;
83  void FormSystemMatrix(const Array<int> &ess_tdof_list, OperatorHandle &A);
84  void FormLinearSystem(const Array<int> &ess_tdof_list,
85  Vector &x, Vector &b,
86  OperatorHandle &A, Vector &X, Vector &B,
87  int copy_interior = 0);
88  void Mult(const Vector &x, Vector &y) const;
89  void MultTranspose(const Vector &x, Vector &y) const;
90  void Update();
91 
92 protected:
94 };
95 
96 /// Data and methods for element-assembled bilinear forms
98 {
99 protected:
100  int ne;
101  int elemDofs;
102  // The element matrices are stored row major
105  int faceDofs;
108 
109 public:
111 
112  void Assemble();
113  void Mult(const Vector &x, Vector &y) const;
114  void MultTranspose(const Vector &x, Vector &y) const;
115 };
116 
117 /// Data and methods for fully-assembled bilinear forms
119 {
120 private:
121  SparseMatrix *mat;
122  mutable Vector dg_x, dg_y;
123 
124 public:
126 
127  void Assemble();
128  void RAP(OperatorHandle &A);
129  /** @note Always does `DIAG_ONE` policy to be consistent with
130  `Operator::FormConstrainedSystemOperator`. */
131  void EliminateBC(const Array<int> &ess_dofs, OperatorHandle &A);
132  void FormSystemMatrix(const Array<int> &ess_tdof_list, OperatorHandle &A);
133  void FormLinearSystem(const Array<int> &ess_tdof_list,
134  Vector &x, Vector &b,
135  OperatorHandle &A, Vector &X, Vector &B,
136  int copy_interior = 0);
137  void Mult(const Vector &x, Vector &y) const;
138  void MultTranspose(const Vector &x, Vector &y) const;
139 
140  /** DGMult and DGMultTranspose use the extended L-vector to perform the
141  computation. */
142  void DGMult(const Vector &x, Vector &y) const;
143  void DGMultTranspose(const Vector &x, Vector &y) const;
144 };
145 
146 /// Data and methods for matrix-free bilinear forms
148 {
149 protected:
150  const FiniteElementSpace *trial_fes, *test_fes; // Not owned
151  mutable Vector localX, localY;
154  const Operator *elem_restrict; // Not owned
157 
158 public:
160 
161  void Assemble();
162  void AssembleDiagonal(Vector &diag) const;
163  void FormSystemMatrix(const Array<int> &ess_tdof_list, OperatorHandle &A);
164  void FormLinearSystem(const Array<int> &ess_tdof_list,
165  Vector &x, Vector &b,
166  OperatorHandle &A, Vector &X, Vector &B,
167  int copy_interior = 0);
168  void Mult(const Vector &x, Vector &y) const;
169  void MultTranspose(const Vector &x, Vector &y) const;
170  void Update();
171 };
172 
173 /// Class extending the MixedBilinearForm class to support different AssemblyLevels.
174 /** FA - Full Assembly
175  PA - Partial Assembly
176  EA - Element Assembly
177  MF - Matrix Free
178 */
180 {
181 protected:
182  MixedBilinearForm *a; ///< Not owned
183 
184 public:
186 
187  virtual MemoryClass GetMemoryClass() const
188  { return Device::GetMemoryClass(); }
189 
190  /// Get the finite element space prolongation matrix
191  virtual const Operator *GetProlongation() const;
192 
193  /// Get the finite element space restriction matrix
194  virtual const Operator *GetRestriction() const;
195 
196  /// Get the output finite element space restriction matrix
197  virtual const Operator *GetOutputProlongation() const;
198 
199  /// Get the output finite element space restriction matrix
200  virtual const Operator *GetOutputRestriction() const;
201 
202  virtual void Assemble() = 0;
203  virtual void FormRectangularSystemOperator(const Array<int> &trial_tdof_list,
204  const Array<int> &test_tdof_list,
205  OperatorHandle &A) = 0;
206  virtual void FormRectangularLinearSystem(const Array<int> &trial_tdof_list,
207  const Array<int> &test_tdof_list,
208  Vector &x, Vector &b,
209  OperatorHandle &A, Vector &X, Vector &B) = 0;
210 
211  virtual void AssembleDiagonal_ADAt(const Vector &D, Vector &diag) const = 0;
212 
213  virtual void Update() = 0;
214 };
215 
216 /// Data and methods for partially-assembled mixed bilinear forms
218 {
219 protected:
220  const FiniteElementSpace *trial_fes, *test_fes; // Not owned
222  const Operator *elem_restrict_trial; // Not owned
223  const Operator *elem_restrict_test; // Not owned
224 
225  /// Helper function to set up inputs/outputs for Mult or MultTranspose
226  void SetupMultInputs(const Operator *elem_restrict_x,
227  const Vector &x, Vector &localX,
228  const Operator *elem_restrict_y,
229  Vector &y, Vector &localY, const double c) const;
230 
231 public:
233 
234  /// Partial assembly of all internal integrators
235  void Assemble();
236  /**
237  @brief Setup OperatorHandle A to contain constrained linear operator
238 
239  OperatorHandle A contains matrix-free constrained operator formed for RAP
240  system where ess_tdof_list are in trial space and eliminated from
241  "columns" of A.
242  */
243  void FormRectangularSystemOperator(const Array<int> &trial_tdof_list,
244  const Array<int> &test_tdof_list,
245  OperatorHandle &A);
246  /**
247  Setup OperatorHandle A to contain constrained linear operator and
248  eliminate columns corresponding to essential dofs from system,
249  updating RHS B vector with the results.
250  */
251  void FormRectangularLinearSystem(const Array<int> &trial_tdof_list,
252  const Array<int> &test_tdof_list,
253  Vector &x, Vector &b,
254  OperatorHandle &A, Vector &X, Vector &B);
255  /// y = A*x
256  void Mult(const Vector &x, Vector &y) const;
257  /// y += c*A*x
258  void AddMult(const Vector &x, Vector &y, const double c=1.0) const;
259  /// y = A^T*x
260  void MultTranspose(const Vector &x, Vector &y) const;
261  /// y += c*A^T*x
262  void AddMultTranspose(const Vector &x, Vector &y, const double c=1.0) const;
263  /// Assemble the diagonal of ADA^T for a diagonal vector D.
264  void AssembleDiagonal_ADAt(const Vector &D, Vector &diag) const;
265 
266  /// Update internals for when a new MixedBilinearForm is given to this class
267  void Update();
268 };
269 
270 
271 /**
272  @brief Partial assembly extension for DiscreteLinearOperator
273 
274  This acts very much like PAMixedBilinearFormExtension, but its
275  FormRectangularSystemOperator implementation emulates 'Set' rather than
276  'Add' in the assembly case.
277 */
279 {
280 public:
282 
283  /// Partial assembly of all internal integrators
284  void Assemble();
285 
286  void AddMult(const Vector &x, Vector &y, const double c=1.0) const;
287 
288  void AddMultTranspose(const Vector &x, Vector &y, const double c=1.0) const;
289 
291  OperatorHandle& A);
292 
293  const Operator * GetOutputRestrictionTranspose() const;
294 
295 private:
296  Vector test_multiplicity;
297 };
298 
299 }
300 
301 #endif
static MemoryClass GetMemoryClass()
(DEPRECATED) Equivalent to GetDeviceMemoryClass().
Definition: device.hpp:285
void AddMultTranspose(const Vector &x, Vector &y, const double c=1.0) const
y += c*A^T*x
void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
virtual const Operator * GetRestriction() const
Get the finite element space restriction matrix.
void FormRectangularSystemOperator(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, OperatorHandle &A)
Setup OperatorHandle A to contain constrained linear operator.
Data and methods for matrix-free bilinear forms.
void SetupMultInputs(const Operator *elem_restrict_x, const Vector &x, Vector &localX, const Operator *elem_restrict_y, Vector &y, Vector &localY, const double c) const
Helper function to set up inputs/outputs for Mult or MultTranspose.
const FiniteElementSpace * test_fes
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
MFBilinearFormExtension(BilinearForm *form)
void Assemble()
Assemble at the level given for the BilinearFormExtension subclass.
void Assemble()
Assemble at the level given for the BilinearFormExtension subclass.
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
PAMixedBilinearFormExtension(MixedBilinearForm *form)
void DGMultTranspose(const Vector &x, Vector &y) const
BilinearFormExtension(BilinearForm *form)
void AssembleDiagonal(Vector &diag) const
Computes the diagonal entries into diag. Typically, this operation only makes sense for linear Operat...
virtual MemoryClass GetMemoryClass() const
Return the MemoryClass preferred by the Operator.
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 ...
virtual const Operator * GetOutputProlongation() const
Get the output finite element space restriction matrix.
virtual void AssembleDiagonal(Vector &diag) const
Computes the diagonal entries into diag. Typically, this operation only makes sense for linear Operat...
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 ...
Data and methods for partially-assembled bilinear forms.
virtual void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)=0
Class extending the MixedBilinearForm class to support different AssemblyLevels.
void AssembleDiagonal_ADAt(const Vector &D, Vector &diag) const
Assemble the diagonal of ADA^T for a diagonal vector D.
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
virtual const Operator * GetRestriction() const
Get the finite element space restriction matrix.
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 ...
PADiscreteLinearOperatorExtension(DiscreteLinearOperator *linop)
void Assemble()
Partial assembly of all internal integrators.
virtual void FormRectangularSystemOperator(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, OperatorHandle &A)=0
Data type sparse matrix.
Definition: sparsemat.hpp:50
static MemoryClass GetDeviceMemoryClass()
Get the current Device MemoryClass. This is the MemoryClass used by most MFEM device kernels to acces...
Definition: device.hpp:281
Data and methods for fully-assembled bilinear forms.
double b
Definition: lissajous.cpp:42
BilinearForm * a
Not owned.
virtual void Assemble()=0
Assemble at the level given for the BilinearFormExtension subclass.
void AssembleDiagonal(Vector &diag) const
Computes the diagonal entries into diag. Typically, this operation only makes sense for linear Operat...
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
void EliminateBC(const Array< int > &ess_dofs, OperatorHandle &A)
Class extending the BilinearForm class to support different AssemblyLevels.
void FormRectangularSystemOperator(const Array< int > &, const Array< int > &, OperatorHandle &A)
Setup OperatorHandle A to contain constrained linear operator.
FABilinearFormExtension(BilinearForm *form)
const FaceRestriction * bdr_face_restrict_lex
void Assemble()
Partial assembly of all internal integrators.
void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
void FormRectangularLinearSystem(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B)
const Operator * GetOutputRestrictionTranspose() const
Transpose of GetOutputRestriction, directly available in this form to facilitate matrix-free RAP-type...
void Assemble()
Assemble at the level given for the BilinearFormExtension subclass.
void SetupRestrictionOperators(const L2FaceValues m)
Data and methods for element-assembled bilinear forms.
void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
const FiniteElementSpace * test_fes
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
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 ...
virtual const Operator * GetOutputRestriction() const
Get the output finite element space restriction matrix.
void Assemble()
Assemble at the level given for the BilinearFormExtension subclass.
void Mult(const Vector &x, Vector &y) const
y = A*x
const FaceRestriction * int_face_restrict_lex
void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
EABilinearFormExtension(BilinearForm *form)
MixedBilinearForm * a
Not owned.
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
const FaceRestriction * bdr_face_restrict_lex
const FaceRestriction * int_face_restrict_lex
void RAP(OperatorHandle &A)
const FiniteElementSpace * test_fes
void Update()
Update internals for when a new MixedBilinearForm is given to this class.
void AddMult(const Vector &x, Vector &y, const double c=1.0) const
y += c*A*x
void AddMult(const Vector &x, Vector &y, const double c=1.0) const
y += c*A*x
virtual const Operator * GetProlongation() const
Get the finite element space prolongation matrix.
virtual const Operator * GetProlongation() const
Get the finite element space prolongation matrix.
Data and methods for partially-assembled mixed bilinear forms.
Vector data type.
Definition: vector.hpp:60
void DGMult(const Vector &x, Vector &y) const
virtual MemoryClass GetMemoryClass() const
Return the MemoryClass preferred by the Operator.
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)=0
void MultTranspose(const Vector &x, Vector &y) const
y = A^T*x
const FiniteElementSpace * trial_fes
Base class for operators that extracts Face degrees of freedom.
Abstract operator.
Definition: operator.hpp:24
const FiniteElementSpace * trial_fes
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
MixedBilinearFormExtension(MixedBilinearForm *form)
void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
MemoryClass
Memory classes identify sets of memory types.
Definition: mem_manager.hpp:73
const FiniteElementSpace * trial_fes
void AddMultTranspose(const Vector &x, Vector &y, const double c=1.0) const
y += c*A^T*x
virtual void FormRectangularLinearSystem(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B)=0
Partial assembly extension for DiscreteLinearOperator.
virtual void AssembleDiagonal_ADAt(const Vector &D, Vector &diag) const =0