MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
bilinearform_ext.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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 AddMult(const Vector &x, Vector &y, const double c=1.0) const = 0;
212  virtual void AddMultTranspose(const Vector &x, Vector &y,
213  const double c=1.0) const = 0;
214 
215  virtual void AssembleDiagonal_ADAt(const Vector &D, Vector &diag) const = 0;
216 
217  virtual void Update() = 0;
218 };
219 
220 /// Data and methods for partially-assembled mixed bilinear forms
222 {
223 protected:
224  const FiniteElementSpace *trial_fes, *test_fes; // Not owned
226  const Operator *elem_restrict_trial; // Not owned
227  const Operator *elem_restrict_test; // Not owned
228 
229  /// Helper function to set up inputs/outputs for Mult or MultTranspose
230  void SetupMultInputs(const Operator *elem_restrict_x,
231  const Vector &x, Vector &localX,
232  const Operator *elem_restrict_y,
233  Vector &y, Vector &localY, const double c) const;
234 
235 public:
237 
238  /// Partial assembly of all internal integrators
239  void Assemble();
240  /**
241  @brief Setup OperatorHandle A to contain constrained linear operator
242 
243  OperatorHandle A contains matrix-free constrained operator formed for RAP
244  system where ess_tdof_list are in trial space and eliminated from
245  "columns" of A.
246  */
247  void FormRectangularSystemOperator(const Array<int> &trial_tdof_list,
248  const Array<int> &test_tdof_list,
249  OperatorHandle &A);
250  /**
251  Setup OperatorHandle A to contain constrained linear operator and
252  eliminate columns corresponding to essential dofs from system,
253  updating RHS B vector with the results.
254  */
255  void FormRectangularLinearSystem(const Array<int> &trial_tdof_list,
256  const Array<int> &test_tdof_list,
257  Vector &x, Vector &b,
258  OperatorHandle &A, Vector &X, Vector &B);
259  /// y = A*x
260  void Mult(const Vector &x, Vector &y) const;
261  /// y += c*A*x
262  void AddMult(const Vector &x, Vector &y, const double c=1.0) const;
263  /// y = A^T*x
264  void MultTranspose(const Vector &x, Vector &y) const;
265  /// y += c*A^T*x
266  void AddMultTranspose(const Vector &x, Vector &y, const double c=1.0) const;
267  /// Assemble the diagonal of ADA^T for a diagonal vector D.
268  void AssembleDiagonal_ADAt(const Vector &D, Vector &diag) const;
269 
270  /// Update internals for when a new MixedBilinearForm is given to this class
271  void Update();
272 };
273 
274 
275 /**
276  @brief Partial assembly extension for DiscreteLinearOperator
277 
278  This acts very much like PAMixedBilinearFormExtension, but its
279  FormRectangularSystemOperator implementation emulates 'Set' rather than
280  'Add' in the assembly case.
281 */
283 {
284 public:
286 
287  /// Partial assembly of all internal integrators
288  void Assemble();
289 
290  void AddMult(const Vector &x, Vector &y, const double c) const;
291 
292  void AddMultTranspose(const Vector &x, Vector &y, const double c=1.0) const;
293 
295  OperatorHandle& A);
296 
297  const Operator * GetOutputRestrictionTranspose() const;
298 
299 private:
300  Vector test_multiplicity;
301 };
302 
303 }
304 
305 #endif
static MemoryClass GetMemoryClass()
(DEPRECATED) Equivalent to GetDeviceMemoryClass().
Definition: device.hpp:285
void AssembleDiagonal_ADAt(const Vector &D, Vector &diag) const
Assemble the diagonal of ADA^T for a diagonal vector D.
void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
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.
const FiniteElementSpace * test_fes
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
void DGMultTranspose(const Vector &x, Vector &y) const
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.
void AddMultTranspose(const Vector &x, Vector &y, const double c=1.0) const
y += c*A^T*x
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
PAMixedBilinearFormExtension(MixedBilinearForm *form)
void MultTranspose(const Vector &x, Vector &y) const
y = A^T*x
virtual void AddMult(const Vector &x, Vector &y, const double c=1.0) const =0
virtual const Operator * GetRestriction() const
Get the finite element space restriction matrix.
virtual const Operator * GetRestriction() const
Get the finite element space restriction matrix.
BilinearFormExtension(BilinearForm *form)
virtual void AddMultTranspose(const Vector &x, Vector &y, const double c=1.0) const =0
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.
const Operator * GetOutputRestrictionTranspose() const
Transpose of GetOutputRestriction, directly available in this form to facilitate matrix-free RAP-type...
virtual const Operator * GetOutputProlongation() const
Get the output 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 ...
virtual const Operator * GetProlongation() const
Get the finite element space prolongation matrix.
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
virtual void AssembleDiagonal(Vector &diag) const
Computes the diagonal entries into diag. Typically, this operation only makes sense for linear Operat...
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.
void AssembleDiagonal(Vector &diag) const
Computes the diagonal entries into diag. Typically, this operation only makes sense for linear Operat...
virtual void Assemble()=0
Assemble at the level given for the BilinearFormExtension subclass.
void AddMult(const Vector &x, Vector &y, const double c=1.0) const
y += c*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.
virtual const Operator * GetProlongation() const
Get the finite element space prolongation matrix.
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)
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 ...
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 ...
void Assemble()
Assemble at the level given for the BilinearFormExtension subclass.
void SetupRestrictionOperators(const L2FaceValues m)
void Mult(const Vector &x, Vector &y) const
y = A*x
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
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 Assemble()
Assemble at the level given for the BilinearFormExtension subclass.
const FaceRestriction * int_face_restrict_lex
void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
EABilinearFormExtension(BilinearForm *form)
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
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.
MixedBilinearForm * a
Not owned.
A &quot;square matrix&quot; 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 DGMult(const Vector &x, Vector &y) const
void AddMultTranspose(const Vector &x, Vector &y, const double c=1.0) const
y += c*A^T*x
virtual const Operator * GetOutputRestriction() const
Get the output finite element space restriction matrix.
void RAP(OperatorHandle &A)
const FiniteElementSpace * test_fes
void Update()
Update internals for when a new MixedBilinearForm is given to this class.
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 ...
Data and methods for partially-assembled mixed bilinear forms.
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Vector data type.
Definition: vector.hpp:60
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)=0
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.
const FiniteElementSpace * trial_fes
Base class for operators that extracts Face degrees of freedom.
Abstract operator.
Definition: operator.hpp:24
const FiniteElementSpace * trial_fes
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 AddMult(const Vector &x, Vector &y, const double c) const
y += c*A*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