MFEM  v4.6.0
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  /// Attributes of all mesh elements.
73  mutable Vector tmp_evec; // Work array
74  mutable Vector localX, localY;
77  const Operator *elem_restrict; // Not owned
80 
81 public:
83 
84  void Assemble();
85  void AssembleDiagonal(Vector &diag) const;
86  void FormSystemMatrix(const Array<int> &ess_tdof_list, OperatorHandle &A);
87  void FormLinearSystem(const Array<int> &ess_tdof_list,
88  Vector &x, Vector &b,
89  OperatorHandle &A, Vector &X, Vector &B,
90  int copy_interior = 0);
91  void Mult(const Vector &x, Vector &y) const;
92  void MultTranspose(const Vector &x, Vector &y) const;
93  void Update();
94 
95 protected:
97 
98  /// @brief Accumulate the action (or transpose) of the integrator on @a x
99  /// into @a y, taking into account the (possibly null) @a markers array.
100  ///
101  /// If @a markers is non-null, then only those elements or boundary elements
102  /// whose attribute is marked in the markers array will be added to @a y.
103  ///
104  /// @param integ The integrator (domain, boundary, or boundary face).
105  /// @param x Input E-vector.
106  /// @param markers Marked attributes (possibly null, meaning all attributes).
107  /// @param attributes Array of element or boundary element attributes.
108  /// @param transpose Compute the action or transpose of the integrator .
109  /// @param y Output E-vector
110  void AddMultWithMarkers(const BilinearFormIntegrator &integ,
111  const Vector &x,
112  const Array<int> *markers,
113  const Array<int> &attributes,
114  const bool transpose,
115  Vector &y) const;
116 };
117 
118 /// Data and methods for element-assembled bilinear forms
120 {
121 protected:
122  int ne;
123  int elemDofs;
124  // The element matrices are stored row major
127  int faceDofs;
130 
131 public:
133 
134  void Assemble();
135  void Mult(const Vector &x, Vector &y) const;
136  void MultTranspose(const Vector &x, Vector &y) const;
137 };
138 
139 /// Data and methods for fully-assembled bilinear forms
141 {
142 private:
143  SparseMatrix *mat;
144  mutable Vector dg_x, dg_y;
145 
146 public:
148 
149  void Assemble();
150  void RAP(OperatorHandle &A);
151  /** @note Always does `DIAG_ONE` policy to be consistent with
152  `Operator::FormConstrainedSystemOperator`. */
153  void EliminateBC(const Array<int> &ess_dofs, OperatorHandle &A);
154  void FormSystemMatrix(const Array<int> &ess_tdof_list, OperatorHandle &A);
155  void FormLinearSystem(const Array<int> &ess_tdof_list,
156  Vector &x, Vector &b,
157  OperatorHandle &A, Vector &X, Vector &B,
158  int copy_interior = 0);
159  void Mult(const Vector &x, Vector &y) const;
160  void MultTranspose(const Vector &x, Vector &y) const;
161 
162  /** DGMult and DGMultTranspose use the extended L-vector to perform the
163  computation. */
164  void DGMult(const Vector &x, Vector &y) const;
165  void DGMultTranspose(const Vector &x, Vector &y) const;
166 };
167 
168 /// Data and methods for matrix-free bilinear forms
170 {
171 protected:
172  const FiniteElementSpace *trial_fes, *test_fes; // Not owned
173  mutable Vector localX, localY;
176  const Operator *elem_restrict; // Not owned
179 
180 public:
182 
183  void Assemble();
184  void AssembleDiagonal(Vector &diag) const;
185  void FormSystemMatrix(const Array<int> &ess_tdof_list, OperatorHandle &A);
186  void FormLinearSystem(const Array<int> &ess_tdof_list,
187  Vector &x, Vector &b,
188  OperatorHandle &A, Vector &X, Vector &B,
189  int copy_interior = 0);
190  void Mult(const Vector &x, Vector &y) const;
191  void MultTranspose(const Vector &x, Vector &y) const;
192  void Update();
193 };
194 
195 /// Class extending the MixedBilinearForm class to support different AssemblyLevels.
196 /** FA - Full Assembly
197  PA - Partial Assembly
198  EA - Element Assembly
199  MF - Matrix Free
200 */
202 {
203 protected:
204  MixedBilinearForm *a; ///< Not owned
205 
206 public:
208 
209  virtual MemoryClass GetMemoryClass() const
210  { return Device::GetMemoryClass(); }
211 
212  /// Get the finite element space prolongation matrix
213  virtual const Operator *GetProlongation() const;
214 
215  /// Get the finite element space restriction matrix
216  virtual const Operator *GetRestriction() const;
217 
218  /// Get the output finite element space restriction matrix
219  virtual const Operator *GetOutputProlongation() const;
220 
221  /// Get the output finite element space restriction matrix
222  virtual const Operator *GetOutputRestriction() const;
223 
224  virtual void Assemble() = 0;
225  virtual void FormRectangularSystemOperator(const Array<int> &trial_tdof_list,
226  const Array<int> &test_tdof_list,
227  OperatorHandle &A) = 0;
228  virtual void FormRectangularLinearSystem(const Array<int> &trial_tdof_list,
229  const Array<int> &test_tdof_list,
230  Vector &x, Vector &b,
231  OperatorHandle &A, Vector &X, Vector &B) = 0;
232 
233  virtual void AssembleDiagonal_ADAt(const Vector &D, Vector &diag) const = 0;
234 
235  virtual void Update() = 0;
236 };
237 
238 /// Data and methods for partially-assembled mixed bilinear forms
240 {
241 protected:
242  const FiniteElementSpace *trial_fes, *test_fes; // Not owned
244  const Operator *elem_restrict_trial; // Not owned
245  const Operator *elem_restrict_test; // Not owned
246 
247  /// Helper function to set up inputs/outputs for Mult or MultTranspose
248  void SetupMultInputs(const Operator *elem_restrict_x,
249  const Vector &x, Vector &localX,
250  const Operator *elem_restrict_y,
251  Vector &y, Vector &localY, const double c) const;
252 
253 public:
255 
256  /// Partial assembly of all internal integrators
257  void Assemble();
258  /**
259  @brief Setup OperatorHandle A to contain constrained linear operator
260 
261  OperatorHandle A contains matrix-free constrained operator formed for RAP
262  system where ess_tdof_list are in trial space and eliminated from
263  "columns" of A.
264  */
265  void FormRectangularSystemOperator(const Array<int> &trial_tdof_list,
266  const Array<int> &test_tdof_list,
267  OperatorHandle &A);
268  /**
269  Setup OperatorHandle A to contain constrained linear operator and
270  eliminate columns corresponding to essential dofs from system,
271  updating RHS B vector with the results.
272  */
273  void FormRectangularLinearSystem(const Array<int> &trial_tdof_list,
274  const Array<int> &test_tdof_list,
275  Vector &x, Vector &b,
276  OperatorHandle &A, Vector &X, Vector &B);
277  /// y = A*x
278  void Mult(const Vector &x, Vector &y) const;
279  /// y += c*A*x
280  void AddMult(const Vector &x, Vector &y, const double c=1.0) const;
281  /// y = A^T*x
282  void MultTranspose(const Vector &x, Vector &y) const;
283  /// y += c*A^T*x
284  void AddMultTranspose(const Vector &x, Vector &y, const double c=1.0) const;
285  /// Assemble the diagonal of ADA^T for a diagonal vector D.
286  void AssembleDiagonal_ADAt(const Vector &D, Vector &diag) const;
287 
288  /// Update internals for when a new MixedBilinearForm is given to this class
289  void Update();
290 };
291 
292 
293 /**
294  @brief Partial assembly extension for DiscreteLinearOperator
295 
296  This acts very much like PAMixedBilinearFormExtension, but its
297  FormRectangularSystemOperator implementation emulates 'Set' rather than
298  'Add' in the assembly case.
299 */
301 {
302 public:
304 
305  /// Partial assembly of all internal integrators
306  void Assemble();
307 
308  void AddMult(const Vector &x, Vector &y, const double c=1.0) const;
309 
310  void AddMultTranspose(const Vector &x, Vector &y, const double c=1.0) const;
311 
313  OperatorHandle& A);
314 
315  const Operator * GetOutputRestrictionTranspose() const;
316 
317 private:
318  Vector test_multiplicity;
319 };
320 
321 }
322 
323 #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
Array< int > elem_attributes
Attributes of all mesh elements.
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.
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:24
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:219
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
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:58
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
void AddMultWithMarkers(const BilinearFormIntegrator &integ, const Vector &x, const Array< int > *markers, const Array< int > &attributes, const bool transpose, Vector &y) const
Accumulate the action (or transpose) of the integrator on x into y, taking into account the (possibly...
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