MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
weakform.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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_DPGWEAKFORM
13#define MFEM_DPGWEAKFORM
14
15#include "mfem.hpp"
16#include "blockstaticcond.hpp"
17
18namespace 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
36protected:
37
39
40 bool initialized = false;
41
42 Mesh * mesh = nullptr;
47
48 /// Block matrix $ M $ to be associated with the Block bilinear form. Owned.
49 BlockMatrix *mat = nullptr;
50
51 /// Block vector $ y $ to be associated with the Block linear form
52 BlockVector * y = nullptr;
53
54 /** @brief Block Matrix $ M_e $ used to store the eliminations
55 from the b.c. Owned.
56 $ M + M_e = M_{original} $ */
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
105public:
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_;
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: $ M $
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.: $ M_e $
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 {
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 {
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 $ M_e $.
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
292
293 virtual ~DPGWeakForm();
294
295};
296
297} // namespace mfem
298
299#endif
Dynamic 2D array using row-major layout.
Definition array.hpp:372
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:697
int Size() const
Return the logical size of the array.
Definition array.hpp:144
Abstract base class BilinearFormIntegrator.
Class that performs static condensation of interior dofs for multiple FE spaces. This class is used i...
A class to handle Vectors in a block fashion.
Class representing the DPG weak formulation. Given the variational formulation a(u,...
Definition weakform.hpp:34
Array2D< Array< BilinearFormIntegrator * > * > test_integs
Set of Test Space (broken) Integrators to be applied for matrix G.
Definition weakform.hpp:73
BlockStaticCondensation * static_cond
Owned.
Definition weakform.hpp:38
BlockMatrix * mat
Block matrix to be associated with the Block bilinear form. Owned.
Definition weakform.hpp:49
void StoreMatrices(bool store_matrices_=true)
Store internal element matrices used for computation of residual after solve.
Definition weakform.hpp:273
Array< int > tdof_offsets
Definition weakform.hpp:46
Array< FiniteElementSpace * > trial_fes
Trial FE spaces.
Definition weakform.hpp:60
virtual void Update()
Update the DPGWeakForm after mesh modifications (AMR)
Definition weakform.cpp:671
Array< Array< LinearFormIntegrator * > * > lfis
Set of Linear Form Integrators to be applied.
Definition weakform.hpp:76
void AddTrialIntegrator(BilinearFormIntegrator *bfi, int n, int m)
Adds new Trial Integrator. Assumes ownership of bfi.
Definition weakform.cpp:105
Array< FiniteElementCollection * > test_fecols
Test FE Collections (Broken)
Definition weakform.hpp:66
mfem::Operator::DiagonalPolicy diag_policy
Definition weakform.hpp:83
DPGWeakForm()
Default constructor. User must call SetSpaces to setup the FE spaces.
Definition weakform.hpp:107
Array< DenseMatrix * > Bmat
Definition weakform.hpp:101
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
Vector & ComputeResidual(const BlockVector &x)
Compute DPG residual based error estimator.
Definition weakform.cpp:723
Array2D< Array< BilinearFormIntegrator * > * > trial_integs
Set of Trial Integrators to be applied for matrix B.
Definition weakform.hpp:70
virtual void BuildProlongation()
Definition weakform.cpp:133
Array< int > dof_offsets
Definition weakform.hpp:45
void AddDomainLFIntegrator(LinearFormIntegrator *lfi, int n)
Adds new Domain LF Integrator. Assumes ownership of bfi.
Definition weakform.cpp:125
void SetDiagonalPolicy(Operator::DiagonalPolicy policy)
Sets diagonal policy used upon construction of the linear system.
Definition weakform.hpp:264
void AllocMat()
Allocate appropriate BlockMatrix and assign it to mat.
Definition weakform.cpp:77
void ConformingAssemble()
Definition weakform.cpp:151
BlockMatrix & BlockMatElim()
Returns a reference to the sparse matrix of eliminated b.c.: .
Definition weakform.hpp:163
int Size() const
Get the size of the bilinear form of the DPGWeakForm.
Definition weakform.hpp:147
Array< Vector * > fvec
Definition weakform.hpp:102
void EnableStaticCondensation()
Definition weakform.cpp:717
BlockMatrix * R
Block Restriction.
Definition weakform.hpp:81
virtual void RecoverFEMSolution(const Vector &X, Vector &x)
Recover the solution of a linear system formed with FormLinearSystem().
Definition weakform.cpp:598
Array< int > test_fecols_vdims
Definition weakform.hpp:67
virtual ~DPGWeakForm()
Definition weakform.cpp:811
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 SetTestFECollVdim(int test_fec, int vdim)
Definition weakform.hpp:120
void AddTestIntegrator(BilinearFormIntegrator *bfi, int n, int m)
Adds new Test Integrator. Assumes ownership of bfi.
Definition weakform.cpp:117
void ComputeOffsets()
Definition weakform.cpp:61
BlockVector * y
Block vector to be associated with the Block linear form.
Definition weakform.hpp:52
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
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
void ReleaseInitMemory()
Definition weakform.cpp:629
BlockMatrix & BlockMat()
Returns a reference to the BlockMatrix: .
Definition weakform.hpp:156
void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
Definition weakform.cpp:97
void AllocateMatrix()
Pre-allocate the internal BlockMatrix before assembly.
Definition weakform.hpp:150
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all integrators.
Definition weakform.cpp:247
BlockMatrix * mat_e
Block Matrix used to store the eliminations from the b.c. Owned. .
Definition weakform.hpp:57
DPGWeakForm(Array< FiniteElementSpace * > &fes_, Array< FiniteElementCollection * > &fecol_)
Creates bilinear form associated with FE spaces fes_.
Definition weakform.hpp:114
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
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
void SetSpaces(Array< FiniteElementSpace * > &fes_, Array< FiniteElementCollection * > &fecol_)
Definition weakform.hpp:125
Array< int > IsTraceFes
Flags to determine if a FiniteElementSpace is Trace.
Definition weakform.hpp:63
BlockMatrix * P
Block Prolongation.
Definition weakform.hpp:79
Arbitrary order "H^{1/2}-conforming" trace finite elements defined on the interface between mesh elem...
Definition fe_coll.hpp:322
Abstract base class LinearFormIntegrator.
Definition lininteg.hpp:25
Mesh data type.
Definition mesh.hpp:56
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
Arbitrary order H(curl)-trace finite elements defined on the interface between mesh elements (faces,...
Definition fe_coll.hpp:511
Pointer to an Operator of a specified type.
Definition handle.hpp:34
OpType * Is() const
Return the Operator pointer dynamically cast to a specified OpType.
Definition handle.hpp:108
DiagonalPolicy
Defines operator diagonal policy upon elimination of rows and/or columns.
Definition operator.hpp:48
@ DIAG_ONE
Set the diagonal value to one.
Definition operator.hpp:50
Arbitrary order "H^{-1/2}-conforming" face finite elements defined on the interface between mesh elem...
Definition fe_coll.hpp:445
Vector data type.
Definition vector.hpp:80
real_t b
Definition lissajous.cpp:42