MFEM  v4.6.0 Finite element discretization library
weakform.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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_DPGWEAKFORM
13 #define MFEM_DPGWEAKFORM
14
15 #include "mfem.hpp"
16 #include "blockstaticcond.hpp"
17
18 namespace mfem
19 {
20
21 /** @brief Class representing the DPG weak formulation.
22  Given the variational formulation
23  a(u,v) = b(v), (or A u = b, where <Au,v> = a(u,v))
24  this class forms the DPG linear system
25  A^T G^-1 A u = A^T G^-1 b
26  This system results from the minimum residual formulation
27  u = argmin_w ||G^-1(b - Aw)||.
28  Here G is a symmetric positive definite matrix resulting from the discretization of
29  the Riesz operator on the test space. Since the test space is broken
30  (discontinuous), G is defined and inverted element-wise and the assembly
31  of the global system is performed in the same manner as the standard FEM method.
32  Note that DPGWeakForm can handle multiple Finite Element spaces.*/
34 {
35
36 protected:
37
39
40  bool initialized = false;
41
42  Mesh * mesh = nullptr;
43  int height, width;
44  int nblocks;
47
48  /// Block matrix \f$M \f$ to be associated with the Block bilinear form. Owned.
49  BlockMatrix *mat = nullptr;
50
51  /// Block vector \f$y \f$ to be associated with the Block linear form
52  BlockVector * y = nullptr;
53
54  /** @brief Block Matrix \f$M_e \f$ used to store the eliminations
55  from the b.c. Owned.
56  \f$M + M_e = M_{original} \f$ */
57  BlockMatrix *mat_e = nullptr;
58
59  /// Trial FE spaces
61
62  /// Flags to determine if a FiniteElementSpace is Trace
64
65  /// Test FE Collections (Broken)
68
69  /// Set of Trial Integrators to be applied for matrix B.
71
72  /// Set of Test Space (broken) Integrators to be applied for matrix G
74
75  /// Set of Linear Form Integrators to be applied.
77
78  /// Block Prolongation
79  BlockMatrix * P = nullptr;
80  /// Block Restriction
81  BlockMatrix * R = nullptr;
82
84
85  void Init();
86  void ReleaseInitMemory();
87
88  /// Allocate appropriate BlockMatrix and assign it to mat
89  void AllocMat();
90
91  void ConformingAssemble();
92
93  void ComputeOffsets();
94
95  virtual void BuildProlongation();
96
97  bool store_matrices = false;
98
99  /** Store the matrix L^-1 B and Vector L^-1 l
100  where G = L L^t */
104
105 public:
106  /// Default constructor. User must call SetSpaces to setup the FE spaces
108  {
109  height = 0;
110  width = 0;
111  }
112
113  /// Creates bilinear form associated with FE spaces @a fes_.
116  {
117  SetSpaces(fes_,fecol_);
118  }
119
120  void SetTestFECollVdim(int test_fec, int vdim)
121  {
122  test_fecols_vdims[test_fec] = vdim;
123  }
124
127  {
128  trial_fes = fes_;
129  test_fecols = fecol_;
131  test_fecols_vdims = 1;
132  nblocks = trial_fes.Size();
133  mesh = trial_fes[0]->GetMesh();
134
136  for (int i = 0; i < nblocks; i++)
137  {
138  IsTraceFes[i] =
139  (dynamic_cast<const H1_Trace_FECollection*>(trial_fes[i]->FEColl()) ||
140  dynamic_cast<const ND_Trace_FECollection*>(trial_fes[i]->FEColl()) ||
141  dynamic_cast<const RT_Trace_FECollection*>(trial_fes[i]->FEColl()));
142  }
143  Init();
144  }
145
146  /// Get the size of the bilinear form of the DPGWeakForm
147  int Size() const { return height; }
148
149  /// Pre-allocate the internal BlockMatrix before assembly.
150  void AllocateMatrix() { if (mat == nullptr) { AllocMat(); } }
151
152  /// Finalizes the matrix initialization.
153  void Finalize(int skip_zeros = 1);
154
155  /// Returns a reference to the BlockMatrix: \f$M \f$
157  {
158  MFEM_VERIFY(mat, "mat is NULL and can't be dereferenced");
159  return *mat;
160  }
161
162  /// Returns a reference to the sparse matrix of eliminated b.c.: \f$M_e \f$
164  {
165  MFEM_VERIFY(mat_e, "mat_e is NULL and can't be dereferenced");
166  return *mat_e;
167  }
168
169  /// Adds new Trial Integrator. Assumes ownership of @a bfi.
170  void AddTrialIntegrator(BilinearFormIntegrator *bfi, int n, int m);
171
172  /// Adds new Test Integrator. Assumes ownership of @a bfi.
173  void AddTestIntegrator(BilinearFormIntegrator *bfi, int n, int m);
174
175  /// Adds new Domain LF Integrator. Assumes ownership of @a bfi.
177
178  /// Assembles the form i.e. sums over all integrators.
179  void Assemble(int skip_zeros = 1);
180
181  /** @brief Form the linear system A X = B, corresponding to this DPG weak
182  form */
183  /** This method applies any necessary transformations to the linear system
184  such as: eliminating boundary conditions; applying conforming constraints
185  for non-conforming AMR; static condensation;
186
187  The GridFunction-size vector @a x must contain the essential b.c. The
188  DPGWeakForm must be assembled.
189
190  The vector @a X is initialized with a suitable initial guess: the essential
191  entries of @a X are set to the corresponding b.c. and all other entries
192  are set to zero (@a copy_interior == 0) or copied from @a x
193  (@a copy_interior != 0).
194
195  After solving the linear system, the finite element solution @a x can be
196  recovered by calling RecoverFEMSolution() (with the same vectors @a X,
197  and @a x).
198
199  NOTE: If there are no transformations, @a X simply reuses the data of
200  @a x. */
201  virtual void FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x,
202  OperatorHandle &A, Vector &X,
203  Vector &B, int copy_interior = 0);
204
205  /** @brief Form the linear system A X = B, corresponding to this DPG weak form
206  Version of the method FormLinearSystem() where the system matrix is
207  returned in the variable @a A, of type OpType, holding a *reference* to
208  the system matrix (created with the method OpType::MakeRef()). */
209  template <typename OpType>
210  void FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x,
211  OpType &A, Vector &X, Vector &B,
212  int copy_interior = 0)
213  {
214  OperatorHandle Ah;
215  FormLinearSystem(ess_tdof_list, x, Ah, X, B, copy_interior);
216  OpType *A_ptr = Ah.Is<OpType>();
217  MFEM_VERIFY(A_ptr, "invalid OpType used");
218  A.MakeRef(*A_ptr);
219  }
220
221  /// Form the linear system matrix @a A, see FormLinearSystem() for details.
222  virtual void FormSystemMatrix(const Array<int> &ess_tdof_list,
223  OperatorHandle &A);
224
225  /// Form the linear system matrix A, see FormLinearSystem() for details.
226  /** Version of the method FormSystemMatrix() where the system matrix is
227  returned in the variable @a A, of type OpType, holding a *reference* to
228  the system matrix (created with the method OpType::MakeRef()). */
229  template <typename OpType>
230  void FormSystemMatrix(const Array<int> &ess_tdof_list, OpType &A)
231  {
232  OperatorHandle Ah;
233  FormSystemMatrix(ess_tdof_list, Ah);
234  OpType *A_ptr = Ah.Is<OpType>();
235  MFEM_VERIFY(A_ptr, "invalid OpType used");
236  A.MakeRef(*A_ptr);
237  }
238
239  /// Eliminate the given @a vdofs, storing the eliminated part internally in \f$M_e \f$.
240  /** This method works in conjunction with EliminateVDofsInRHS() and allows
241  elimination of boundary conditions in multiple right-hand sides. In this
242  method, @a vdofs is a list of DOFs. */
243  void EliminateVDofs(const Array<int> &vdofs,
245
246  /** @brief Use the stored eliminated part of the matrix (see
247  EliminateVDofs(const Array<int> &, DiagonalPolicy)) to modify the r.h.s.
248  @a b; @a vdofs is a list of DOFs. */
249  void EliminateVDofsInRHS(const Array<int> &vdofs, const Vector &x, Vector &b);
250
251  /// Recover the solution of a linear system formed with FormLinearSystem().
252  /** Call this method after solving a linear system constructed using the
253  FormLinearSystem() method to recover the solution as a GridFunction-size
254  vector in @a x. Use the same arguments as in the FormLinearSystem() call.
255  */
256  virtual void RecoverFEMSolution(const Vector &X,Vector &x);
257
258  /// Sets diagonal policy used upon construction of the linear system.
259  /** Policies include:
260  - DIAG_ZERO (Set the diagonal values to zero)
261  - DIAG_ONE (Set the diagonal values to one)
262  - DIAG_KEEP (Keep the diagonal values)
263  */
265  {
266  diag_policy = policy;
267  }
268
269  /// Update the DPGWeakForm after mesh modifications (AMR)
270  virtual void Update();
271
272  /// Store internal element matrices used for computation of residual after solve
273  void StoreMatrices(bool store_matrices_ = true)
274  {
275  store_matrices = store_matrices_;
276  if (Bmat.Size() == 0)
277  {
278  Bmat.SetSize(mesh->GetNE());
279  fvec.SetSize(mesh->GetNE());
280  for (int i =0; i<mesh->GetNE(); i++)
281  {
282  Bmat[i] = nullptr;
283  fvec[i] = nullptr;
284  }
285  }
286  }
287
289
290  /// Compute DPG residual based error estimator
291  Vector & ComputeResidual(const BlockVector & x);
292
293  virtual ~DPGWeakForm();
294
295 };
296
297 } // namespace mfem
298
299 #endif
void AllocMat()
Allocate appropriate BlockMatrix and assign it to mat.
Definition: weakform.cpp:77
void ComputeOffsets()
Definition: weakform.cpp:61
void EliminateVDofsInRHS(const Array< int > &vdofs, const Vector &x, Vector &b)
Use the stored eliminated part of the matrix (see EliminateVDofs(const Array<int> &...
Definition: weakform.cpp:567
BlockMatrix & BlockMat()
Returns a reference to the BlockMatrix: .
Definition: weakform.hpp:156
DPGWeakForm(Array< FiniteElementSpace * > &fes_, Array< FiniteElementCollection *> &fecol_)
Creates bilinear form associated with FE spaces fes_.
Definition: weakform.hpp:114
A class to handle Vectors in a block fashion.
Definition: blockvector.hpp:30
virtual void RecoverFEMSolution(const Vector &X, Vector &x)
Recover the solution of a linear system formed with FormLinearSystem().
Definition: weakform.cpp:598
void AddTestIntegrator(BilinearFormIntegrator *bfi, int n, int m)
Adds new Test Integrator. Assumes ownership of bfi.
Definition: weakform.cpp:117
BlockMatrix * P
Block Prolongation.
Definition: weakform.hpp:79
void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
Definition: weakform.cpp:97
BlockMatrix * R
Block Restriction.
Definition: weakform.hpp:81
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
OpType * Is() const
Return the Operator pointer dynamically cast to a specified OpType.
Definition: handle.hpp:108
void FormSystemMatrix(const Array< int > &ess_tdof_list, OpType &A)
Form the linear system matrix A, see FormLinearSystem() for details.
Definition: weakform.hpp:230
void SetSpaces(Array< FiniteElementSpace * > &fes_, Array< FiniteElementCollection *> &fecol_)
Definition: weakform.hpp:125
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
Definition: weakform.cpp:531
Array< int > IsTraceFes
Flags to determine if a FiniteElementSpace is Trace.
Definition: weakform.hpp:63
void AllocateMatrix()
Pre-allocate the internal BlockMatrix before assembly.
Definition: weakform.hpp:150
Array< FiniteElementSpace *> trial_fes
Trial FE spaces.
Definition: weakform.hpp:60
Class that performs static condensation of interior dofs for multiple FE spaces. This class is used i...
Abstract base class LinearFormIntegrator.
Definition: lininteg.hpp:24
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all integrators.
Definition: weakform.cpp:247
void StoreMatrices(bool store_matrices_=true)
Store internal element matrices used for computation of residual after solve.
Definition: weakform.hpp:273
virtual void BuildProlongation()
Definition: weakform.cpp:133
BlockMatrix * mat_e
Block Matrix used to store the eliminations from the b.c. Owned. .
Definition: weakform.hpp:57
double b
Definition: lissajous.cpp:42
Arbitrary order "H^{1/2}-conforming" trace finite elements defined on the interface between mesh elem...
Definition: fe_coll.hpp:319
void ConformingAssemble()
Definition: weakform.cpp:151
virtual void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
Form the linear system A X = B, corresponding to this DPG weak form.
Definition: weakform.cpp:476
BlockMatrix & BlockMatElim()
Returns a reference to the sparse matrix of eliminated b.c.: .
Definition: weakform.hpp:163
void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, OpType &A, Vector &X, Vector &B, int copy_interior=0)
Form the linear system A X = B, corresponding to this DPG weak form Version of the method FormLinearS...
Definition: weakform.hpp:210
Set the diagonal value to one.
Definition: operator.hpp:50
virtual void Update()
Update the DPGWeakForm after mesh modifications (AMR)
Definition: weakform.cpp:671
Adds new Domain LF Integrator. Assumes ownership of bfi.
Definition: weakform.cpp:125
Dynamic 2D array using row-major layout.
Definition: array.hpp:354
Array< DenseMatrix *> Bmat
Definition: weakform.hpp:101
BlockStaticCondensation * static_cond
Owned.
Definition: weakform.hpp:38
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:24
int Size() const
Get the size of the bilinear form of the DPGWeakForm.
Definition: weakform.hpp:147
BlockMatrix * mat
Block matrix to be associated with the Block bilinear form. Owned.
Definition: weakform.hpp:49
Array2D< Array< BilinearFormIntegrator *> *> test_integs
Set of Test Space (broken) Integrators to be applied for matrix G.
Definition: weakform.hpp:73
void SetTestFECollVdim(int test_fec, int vdim)
Definition: weakform.hpp:120
Array< Array< LinearFormIntegrator *> *> lfis
Set of Linear Form Integrators to be applied.
Definition: weakform.hpp:76
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
void EliminateVDofs(const Array< int > &vdofs, Operator::DiagonalPolicy dpolicy=Operator::DIAG_ONE)
Eliminate the given vdofs, storing the eliminated part internally in .
Definition: weakform.cpp:574
Array2D< Array< BilinearFormIntegrator *> *> trial_integs
Set of Trial Integrators to be applied for matrix B.
Definition: weakform.hpp:70
void ReleaseInitMemory()
Definition: weakform.cpp:629
void SetDiagonalPolicy(Operator::DiagonalPolicy policy)
Sets diagonal policy used upon construction of the linear system.
Definition: weakform.hpp:264
virtual ~DPGWeakForm()
Definition: weakform.cpp:811
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
Array< int > test_fecols_vdims
Definition: weakform.hpp:67
DPGWeakForm()
Default constructor. User must call SetSpaces to setup the FE spaces.
Definition: weakform.hpp:107
DiagonalPolicy
Defines operator diagonal policy upon elimination of rows and/or columns.
Definition: operator.hpp:47
void EnableStaticCondensation()
Definition: weakform.cpp:717
Array< int > tdof_offsets
Definition: weakform.hpp:46
Array< int > dof_offsets
Definition: weakform.hpp:45
Vector & ComputeResidual(const BlockVector &x)
Compute DPG residual based error estimator.
Definition: weakform.cpp:723
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
Vector data type.
Definition: vector.hpp:58
Class representing the DPG weak formulation. Given the variational formulation a(u,v) = b(v), (or A u = b, where <Au,v> = a(u,v)) this class forms the DPG linear system A^T G^-1 A u = A^T G^-1 b This system results from the minimum residual formulation u = argmin_w ||G^-1(b - Aw)||. Here G is a symmetric positive definite matrix resulting from the discretization of the Riesz operator on the test space. Since the test space is broken (discontinuous), G is defined and inverted element-wise and the assembly of the global system is performed in the same manner as the standard FEM method. Note that DPGWeakForm can handle multiple Finite Element spaces.
Definition: weakform.hpp:33
BlockVector * y
Block vector to be associated with the Block linear form.
Definition: weakform.hpp:52
mfem::Operator::DiagonalPolicy diag_policy
Definition: weakform.hpp:83
Array< FiniteElementCollection * > test_fecols
Test FE Collections (Broken)
Definition: weakform.hpp:66
void AddTrialIntegrator(BilinearFormIntegrator *bfi, int n, int m)
Adds new Trial Integrator. Assumes ownership of bfi.
Definition: weakform.cpp:105
Array< Vector *> fvec
Definition: weakform.hpp:102