MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
bilinearform.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_BILINEARFORM
13#define MFEM_BILINEARFORM
14
15#include "../config/config.hpp"
16#include "../linalg/linalg.hpp"
17#include "fespace.hpp"
18#include "gridfunc.hpp"
19#include "linearform.hpp"
20#include "bilininteg.hpp"
21#include "bilinearform_ext.hpp"
22#include "staticcond.hpp"
23#include "hybridization.hpp"
24
25namespace mfem
26{
27
28/** @brief Enumeration defining the assembly level for bilinear and nonlinear
29 form classes derived from Operator. For more details, see
30 https://mfem.org/howto/assembly_levels */
31enum class AssemblyLevel
32{
33 /// In the case of a BilinearForm LEGACY corresponds to a fully assembled
34 /// form, i.e. a global sparse matrix in MFEM, Hypre or PETSC format.
35 /// In the case of a NonlinearForm LEGACY corresponds to an operator that
36 /// is fully evaluated on the fly.
37 /// This assembly level is ALWAYS performed on the host.
38 LEGACY = 0,
39 /// @deprecated Use LEGACY instead.
40 LEGACYFULL = 0,
41 /// Fully assembled form, i.e. a global sparse matrix in MFEM format. This
42 /// assembly is compatible with device execution.
43 FULL,
44 /// Form assembled at element level, which computes and stores dense element
45 /// matrices.
46 ELEMENT,
47 /// Partially-assembled form, which computes and stores data only at
48 /// quadrature points.
49 PARTIAL,
50 /// "Matrix-free" form that computes all of its action on-the-fly without any
51 /// substantial storage.
52 NONE,
53};
54
55
56/** @brief A "square matrix" operator for the associated FE space and
57 BLFIntegrators The sum of all the BLFIntegrators can be used form the matrix
58 M. This class also supports other assembly levels specified via the
59 SetAssemblyLevel() function. */
60class BilinearForm : public Matrix
61{
63
64protected:
65 /// Sparse matrix $ M $ to be associated with the form. Owned.
67
68 /** @brief Sparse Matrix $ M_e $ used to store the eliminations
69 from the b.c. Owned.
70 $ M + M_e = M_{original} $ */
72
73 /// FE space on which the form lives. Not owned.
75
76 /** @brief The ::AssemblyLevel of the form (AssemblyLevel::LEGACY,
77 AssemblyLevel::FULL, AssemblyLevel::ELEMENT, AssemblyLevel::PARTIAL) */
79
80 /// Element batch size used in the form action (1, 8, num_elems, etc.)
81 int batch;
82
83 /** @brief Extension for supporting Full Assembly (FA),
84 Element Assembly (EA),Partial Assembly (PA),
85 or Matrix Free assembly (MF). */
87
88 /** Indicates if the sparse matrix is sorted after assembly when using
89 Full Assembly (FA). */
90 bool sort_sparse_matrix = false;
91
92 /** @brief Indicates the Mesh::sequence corresponding to the current state of
93 the BilinearForm. */
95
96 /** @brief Indicates the BilinearFormIntegrator%s stored in #domain_integs,
97 #boundary_integs, #interior_face_integs, and #boundary_face_integs are
98 owned by another BilinearForm. */
100
101 /// Set of Domain Integrators to be applied.
103
104 /// Element attribute marker (should be of length mesh->attributes.Max() or
105 /// 0 if mesh->attributes is empty)
106 /// Includes all by default.
107 /// 0 - ignore attribute
108 /// 1 - include attribute
109 Array<Array<int>*> domain_integs_marker; ///< Entries are not owned.
110
111 /// Set of Boundary Integrators to be applied.
113 Array<Array<int>*> boundary_integs_marker; ///< Entries are not owned.
114
115 /// Set of interior face Integrators to be applied.
117
118 /// Set of boundary face Integrators to be applied.
120 Array<Array<int>*> boundary_face_integs_marker; ///< Entries are not owned.
121
124
126
129
130 /** @brief This data member allows one to specify what should be done to the
131 diagonal matrix entries and corresponding RHS values upon elimination of
132 the constrained DoFs. */
134
136
137 /// Allocate appropriate SparseMatrix and assign it to #mat
138 void AllocMat();
139
140 /** @brief For partially conforming trial and/or test FE spaces, complete the
141 assembly process by performing $ P^t A P $ where $ A $ is the
142 internal sparse matrix and $ P $ is the conforming prolongation
143 matrix of the trial/test FE space. After this call the
144 BilinearForm becomes an operator on the conforming FE space. */
145 void ConformingAssemble();
146
147 /// may be used in the construction of derived classes
149 {
150 fes = NULL; sequence = -1;
151 mat = mat_e = NULL; extern_bfs = 0; element_matrices = NULL;
152 static_cond = NULL; hybridization = NULL;
156 batch = 1;
157 ext = NULL;
158 }
159
160private:
161 /// Copy construction is not supported; body is undefined.
162 BilinearForm(const BilinearForm &);
163
164 /// Copy assignment is not supported; body is undefined.
165 BilinearForm &operator=(const BilinearForm &);
166
167public:
168 /// Creates bilinear form associated with FE space @a *f.
169 /** The pointer @a f is not owned by the newly constructed object. */
171
172 /** @brief Create a BilinearForm on the FiniteElementSpace @a f, using the
173 same integrators as the BilinearForm @a bf.
174
175 The pointer @a f is not owned by the newly constructed object.
176
177 The integrators in @a bf are copied as pointers and they are not owned by
178 the newly constructed BilinearForm.
179
180 The optional parameter @a ps is used to initialize the internal flag
181 #precompute_sparsity, see UsePrecomputedSparsity() for details. */
183
184 /// Get the size of the BilinearForm as a square matrix.
185 int Size() const { return height; }
186
187 /// Set the desired assembly level.
188 /** Valid choices are:
189
190 - AssemblyLevel::LEGACY (default)
191 - AssemblyLevel::FULL
192 - AssemblyLevel::PARTIAL
193 - AssemblyLevel::ELEMENT
194 - AssemblyLevel::NONE
195
196 If used, this method must be called before assembly. */
197 void SetAssemblyLevel(AssemblyLevel assembly_level);
198
199 /** @brief Force the sparse matrix column indices to be sorted when using
200 AssemblyLevel::FULL.
201
202 When assembling on device the assembly algorithm uses atomic operations
203 to insert values in the sparse matrix, which can result in different
204 column index orderings across runs. Calling this method with @a enable_it
205 set to @a true forces a sorting algorithm to be called at the end of the
206 assembly procedure to ensure sorted column indices (and therefore
207 deterministic results).
208 */
209 void EnableSparseMatrixSorting(bool enable_it)
210 {
211 sort_sparse_matrix = enable_it;
212 }
213
214 /// Returns the assembly level
216
218
219 /** @brief Enable the use of static condensation. For details see the
220 description for class StaticCondensation in fem/staticcond.hpp This
221 method should be called before assembly. If the number of unknowns after
222 static condensation is not reduced, it is not enabled. */
224
225 /** @brief Check if static condensation was actually enabled by a previous
226 call to EnableStaticCondensation(). */
228
229 /// Return the trace FE space associated with static condensation.
231 { return static_cond ? static_cond->GetTraceFESpace() : NULL; }
232
233 /// Enable hybridization.
234 /** For details see the description for class
235 Hybridization in fem/hybridization.hpp. This method should be called
236 before assembly. */
237 void EnableHybridization(FiniteElementSpace *constr_space,
238 BilinearFormIntegrator *constr_integ,
239 const Array<int> &ess_tdof_list);
240
241 /** @brief For scalar FE spaces, precompute the sparsity pattern of the
242 matrix (assuming dense element matrices) based on the types of
243 integrators present in the bilinear form. */
245
246 /** @brief Use the given CSR sparsity pattern to allocate the internal
247 SparseMatrix.
248
249 - The @a I and @a J arrays must define a square graph with size equal to
250 GetVSize() of the associated FiniteElementSpace.
251 - This method should be called after enabling static condensation or
252 hybridization, if used.
253 - In the case of static condensation, @a I and @a J are not used.
254 - The ownership of the arrays @a I and @a J remains with the caller. */
255 void UseSparsity(int *I, int *J, bool isSorted);
256
257 /// Use the sparsity of @a A to allocate the internal SparseMatrix.
258 void UseSparsity(SparseMatrix &A);
259
260 /** @brief Pre-allocate the internal SparseMatrix before assembly.
261 If the internal flag #precompute_sparsity
262 is set, the matrix is allocated in CSR format (i.e.
263 finalized) and the entries are initialized with zeros. */
264 void AllocateMatrix() { if (mat == NULL) { AllocMat(); } }
265
266 /// Access all the integrators added with AddDomainIntegrator().
268
269 /** @brief Access all boundary markers added with AddDomainIntegrator().
270 If no marker was specified when the integrator was added, the
271 corresponding pointer (to Array<int>) will be NULL. */
273
274 /// Access all the integrators added with AddBoundaryIntegrator().
276 /** @brief Access all boundary markers added with AddBoundaryIntegrator().
277 If no marker was specified when the integrator was added, the
278 corresponding pointer (to Array<int>) will be NULL. */
280
281 /// Access all integrators added with AddInteriorFaceIntegrator().
283
284 /// Access all integrators added with AddBdrFaceIntegrator().
286
287 /** @brief Access all boundary markers added with AddBdrFaceIntegrator().
288 If no marker was specified when the integrator was added, the
289 corresponding pointer (to Array<int>) will be NULL. */
291 { return &boundary_face_integs_marker; }
292
293 /// Returns a reference to: $ M_{ij} $
294 const real_t &operator()(int i, int j) { return (*mat)(i,j); }
295
296 /// Returns a reference to: $ M_{ij} $
297 virtual real_t &Elem(int i, int j);
298
299 /// Returns constant reference to: $ M_{ij} $
300 virtual const real_t &Elem(int i, int j) const;
301
302 /// Matrix vector multiplication: $ y = M x $
303 virtual void Mult(const Vector &x, Vector &y) const;
304
305 /** @brief Matrix vector multiplication with the original uneliminated
306 matrix. The original matrix is $ M + M_e $ so we have:
307 $ y = M x + M_e x $ */
308 void FullMult(const Vector &x, Vector &y) const
309 { mat->Mult(x, y); mat_e->AddMult(x, y); }
310
311 /// Add the matrix vector multiple to a vector: $ y += a M x $
312 virtual void AddMult(const Vector &x, Vector &y, const real_t a = 1.0) const
313 { mat -> AddMult (x, y, a); }
314
315 /** @brief Add the original uneliminated matrix vector multiple to a vector.
316 The original matrix is $ M + Me $ so we have:
317 $ y += M x + M_e x $ */
318 void FullAddMult(const Vector &x, Vector &y) const
319 { mat->AddMult(x, y); mat_e->AddMult(x, y); }
320
321 /// Add the matrix transpose vector multiplication: $ y += a M^T x $
322 virtual void AddMultTranspose(const Vector & x, Vector & y,
323 const real_t a = 1.0) const
324 { mat->AddMultTranspose(x, y, a); }
325
326 /** @brief Add the original uneliminated matrix transpose vector
327 multiple to a vector. The original matrix is $ M + M_e $
328 so we have: $ y += M^T x + {M_e}^T x $ */
329 void FullAddMultTranspose(const Vector & x, Vector & y) const
330 { mat->AddMultTranspose(x, y); mat_e->AddMultTranspose(x, y); }
331
332 /// Matrix transpose vector multiplication: $ y = M^T x $
333 virtual void MultTranspose(const Vector & x, Vector & y) const;
334
335 /// Compute $ y^T M x $
336 real_t InnerProduct(const Vector &x, const Vector &y) const
337 { return mat->InnerProduct (x, y); }
338
339 /** @brief Returns a pointer to (approximation) of the matrix inverse:
340 $ M^{-1} $ (currently returns NULL) */
341 virtual MatrixInverse *Inverse() const;
342
343 /** @brief Finalizes the matrix initialization if the ::AssemblyLevel is
344 AssemblyLevel::LEGACY.
345 THe matrix that gets finalized is different if you are using static
346 condensation or hybridization.*/
347 virtual void Finalize(int skip_zeros = 1);
348
349 /** @brief Returns a const reference to the sparse matrix: $ M $
350 *
351 This will fail if HasSpMat() is false. */
352 const SparseMatrix &SpMat() const
353 {
354 MFEM_VERIFY(mat, "mat is NULL and can't be dereferenced");
355 return *mat;
356 }
357
358 /** @brief Returns a reference to the sparse matrix: $ M $
359 *
360 This will fail if HasSpMat() is false. */
362 {
363 MFEM_VERIFY(mat, "mat is NULL and can't be dereferenced");
364 return *mat;
365 }
366
367 /** @brief Returns true if the sparse matrix is not null, false otherwise.
368 *
369 @sa SpMat(). */
370 bool HasSpMat()
371 {
372 return mat != nullptr;
373 }
374
375
376 /** @brief Nullifies the internal matrix $ M $ and returns a pointer
377 to it. Used for transferring ownership. */
378 SparseMatrix *LoseMat() { SparseMatrix *tmp = mat; mat = NULL; return tmp; }
379
380 /** @brief Returns a const reference to the sparse matrix of eliminated b.c.:
381 $ M_e $
382
383 This will fail if HasSpMatElim() is false. */
384 const SparseMatrix &SpMatElim() const
385 {
386 MFEM_VERIFY(mat_e, "mat_e is NULL and can't be dereferenced");
387 return *mat_e;
388 }
389
390 /** @brief Returns a reference to the sparse matrix of eliminated b.c.:
391 $ M_e $
392
393 This will fail if HasSpMatElim() is false. */
395 {
396 MFEM_VERIFY(mat_e, "mat_e is NULL and can't be dereferenced");
397 return *mat_e;
398 }
399
400 /** @brief Returns true if the sparse matrix of eliminated b.c.s is not
401 null, false otherwise.
402
403 @sa SpMatElim(). */
405 {
406 return mat_e != nullptr;
407 }
408
409 /// Adds new Domain Integrator. Assumes ownership of @a bfi.
411
412 /// Adds new Domain Integrator restricted to certain elements specified by
413 /// the @a elem_attr_marker.
415 Array<int> &elem_marker);
416
417 /// Adds new Boundary Integrator. Assumes ownership of @a bfi.
419
420 /** @brief Adds new Boundary Integrator, restricted to specific boundary
421 attributes.
422
423 Assumes ownership of @a bfi. The array @a bdr_marker is stored internally
424 as a pointer to the given Array<int> object. */
426 Array<int> &bdr_marker);
427
428 /// Adds new interior Face Integrator. Assumes ownership of @a bfi.
430
431 /// Adds new boundary Face Integrator. Assumes ownership of @a bfi.
433
434 /** @brief Adds new boundary Face Integrator, restricted to specific boundary
435 attributes.
436
437 Assumes ownership of @a bfi. The array @a bdr_marker is stored internally
438 as a pointer to the given Array<int> object. */
440 Array<int> &bdr_marker);
441
442 /// Sets all sparse values of $ M $ and $ M_e $ to 'a'.
443 void operator=(const real_t a)
444 {
445 if (mat != NULL) { *mat = a; }
446 if (mat_e != NULL) { *mat_e = a; }
447 }
448
449 /// Assembles the form i.e. sums over all domain/bdr integrators.
450 void Assemble(int skip_zeros = 1);
451
452 /** @brief Assemble the diagonal of the bilinear form into @a diag. Note that
453 @a diag is a tdof Vector.
454
455 When the AssemblyLevel is not LEGACY, and the mesh has hanging nodes,
456 this method returns |P^T| d_l, where d_l is the diagonal of the form
457 before applying conforming assembly, P^T is the transpose of the
458 conforming prolongation, and |.| denotes the entry-wise absolute value.
459 In general, this is just an approximation of the exact diagonal for this
460 case. */
461 virtual void AssembleDiagonal(Vector &diag) const;
462
463 /// Get the finite element space prolongation operator.
464 virtual const Operator *GetProlongation() const
465 { return fes->GetConformingProlongation(); }
466
467 /// Get the finite element space restriction operator
468 virtual const Operator *GetRestriction() const
469 { return fes->GetConformingRestriction(); }
470
471 /// Get the output finite element space prolongation matrix
472 virtual const Operator *GetOutputProlongation() const
473 { return GetProlongation(); }
474
475 /** @brief Returns the output fe space restriction matrix, transposed
476
477 Logically, this is the transpose of GetOutputRestriction, but in
478 practice it is convenient to have it in transposed form for
479 construction of RAP operators in matrix-free methods. */
482
483 /// Get the output finite element space restriction matrix
484 virtual const Operator *GetOutputRestriction() const
485 { return GetRestriction(); }
486
487 /// Compute serial RAP operator and store it in @a A as a SparseMatrix.
489 {
490 MFEM_ASSERT(mat, "SerialRAP requires the SparseMatrix to be assembled.");
492 A.Reset(mat, false);
493 }
494
495 /** @brief Form the linear system A X = B, corresponding to this bilinear
496 form and the linear form @a b(.). */
497 /** This method applies any necessary transformations to the linear system
498 such as: eliminating boundary conditions; applying conforming constraints
499 for non-conforming AMR; parallel assembly; static condensation;
500 hybridization.
501
502 The GridFunction-size vector @a x must contain the essential b.c. The
503 BilinearForm and the LinearForm-size vector @a b must be assembled.
504
505 The vector @a X is initialized with a suitable initial guess: when using
506 hybridization, the vector @a X is set to zero; otherwise, the essential
507 entries of @a X are set to the corresponding b.c. and all other entries
508 are set to zero (@a copy_interior == 0) or copied from @a x
509 (@a copy_interior != 0).
510
511 This method can be called multiple times (with the same @a ess_tdof_list
512 array) to initialize different right-hand sides and boundary condition
513 values.
514
515 After solving the linear system, the finite element solution @a x can be
516 recovered by calling RecoverFEMSolution() (with the same vectors @a X,
517 @a b, and @a x).
518
519 NOTE: If there are no transformations, @a X simply reuses the data of
520 @a x. */
521 virtual void FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x,
522 Vector &b, OperatorHandle &A, Vector &X,
523 Vector &B, int copy_interior = 0);
524
525 /** @brief Form the linear system A X = B, corresponding to this bilinear
526 form and the linear form @a b(.). */
527 /** Version of the method FormLinearSystem() where the system matrix is
528 returned in the variable @a A, of type OpType, holding a *reference* to
529 the system matrix (created with the method OpType::MakeRef()). The
530 reference will be invalidated when SetOperatorType(), Update(), or the
531 destructor is called. */
532 template <typename OpType>
533 void FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x, Vector &b,
534 OpType &A, Vector &X, Vector &B,
535 int copy_interior = 0)
536 {
538 FormLinearSystem(ess_tdof_list, x, b, Ah, X, B, copy_interior);
539 OpType *A_ptr = Ah.Is<OpType>();
540 MFEM_VERIFY(A_ptr, "invalid OpType used");
541 A.MakeRef(*A_ptr);
542 }
543
544 /// Form the linear system matrix @a A, see FormLinearSystem() for details.
545 virtual void FormSystemMatrix(const Array<int> &ess_tdof_list,
546 OperatorHandle &A);
547
548 /// Form the linear system matrix A, see FormLinearSystem() for details.
549 /** Version of the method FormSystemMatrix() where the system matrix is
550 returned in the variable @a A, of type OpType, holding a *reference* to
551 the system matrix (created with the method OpType::MakeRef()). The
552 reference will be invalidated when SetOperatorType(), Update(), or the
553 destructor is called. */
554 template <typename OpType>
555 void FormSystemMatrix(const Array<int> &ess_tdof_list, OpType &A)
556 {
558 FormSystemMatrix(ess_tdof_list, Ah);
559 OpType *A_ptr = Ah.Is<OpType>();
560 MFEM_VERIFY(A_ptr, "invalid OpType used");
561 A.MakeRef(*A_ptr);
562 }
563
564 /// Recover the solution of a linear system formed with FormLinearSystem().
565 /** Call this method after solving a linear system constructed using the
566 FormLinearSystem() method to recover the solution as a GridFunction-size
567 vector in @a x. Use the same arguments as in the FormLinearSystem() call.
568 */
569 virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x);
570
571 /// Compute and store internally all element matrices.
573
574 /// Free the memory used by the element matrices.
576 { delete element_matrices; element_matrices = NULL; }
577
578 /// Compute the element matrix of the given element
579 /** The element matrix is computed by calling the domain integrators
580 or the one stored internally by a prior call of ComputeElementMatrices()
581 is returned when available.
582 */
583 void ComputeElementMatrix(int i, DenseMatrix &elmat);
584
585 /// Compute the boundary element matrix of the given boundary element
586 void ComputeBdrElementMatrix(int i, DenseMatrix &elmat);
587
588 /// Assemble the given element matrix
589 /** The element matrix @a elmat is assembled for the element @a i, i.e.
590 added to the system matrix. The flag @a skip_zeros skips the zero
591 elements of the matrix, unless they are breaking the symmetry of
592 the system matrix.
593 */
594 void AssembleElementMatrix(int i, const DenseMatrix &elmat,
595 int skip_zeros = 1);
596
597 /// Assemble the given element matrix
598 /** The element matrix @a elmat is assembled for the element @a i, i.e.
599 added to the system matrix. The vdofs of the element are returned
600 in @a vdofs. The flag @a skip_zeros skips the zero elements of the
601 matrix, unless they are breaking the symmetry of the system matrix.
602 */
603 void AssembleElementMatrix(int i, const DenseMatrix &elmat,
604 Array<int> &vdofs, int skip_zeros = 1);
605
606 /// Assemble the given boundary element matrix
607 /** The boundary element matrix @a elmat is assembled for the boundary
608 element @a i, i.e. added to the system matrix. The flag @a skip_zeros
609 skips the zero elements of the matrix, unless they are breaking the
610 symmetry of the system matrix.
611 */
612 void AssembleBdrElementMatrix(int i, const DenseMatrix &elmat,
613 int skip_zeros = 1);
614
615 /// Assemble the given boundary element matrix
616 /** The boundary element matrix @a elmat is assembled for the boundary
617 element @a i, i.e. added to the system matrix. The vdofs of the element
618 are returned in @a vdofs. The flag @a skip_zeros skips the zero elements
619 of the matrix, unless they are breaking the symmetry of the system
620 matrix. */
621 void AssembleBdrElementMatrix(int i, const DenseMatrix &elmat,
622 Array<int> &vdofs, int skip_zeros = 1);
623
624 /// Eliminate essential boundary DOFs from the system.
625 /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
626 the essential part of the boundary. By default, the diagonal at the
627 essential DOFs is set to 1.0. This behavior is controlled by the argument
628 @a dpolicy. */
629 void EliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
630 const Vector &sol, Vector &rhs,
631 DiagonalPolicy dpolicy = DIAG_ONE);
632
633 /// Eliminate essential boundary DOFs from the system matrix.
634 void EliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
635 DiagonalPolicy dpolicy = DIAG_ONE);
636 /// Perform elimination and set the diagonal entry to the given value
637 void EliminateEssentialBCDiag(const Array<int> &bdr_attr_is_ess,
638 real_t value);
639
640 /// Eliminate the given @a vdofs. NOTE: here, @a vdofs is a list of DOFs.
641 /** In this case the eliminations are applied to the internal $ M $
642 and @a rhs without storing the elimination matrix $ M_e $. */
643 void EliminateVDofs(const Array<int> &vdofs, const Vector &sol, Vector &rhs,
644 DiagonalPolicy dpolicy = DIAG_ONE);
645
646 /** @brief Eliminate the given @a vdofs, storing the eliminated part
647 internally in $ M_e $.
648
649 This method works in conjunction with EliminateVDofsInRHS() and allows
650 elimination of boundary conditions in multiple right-hand sides. In this
651 method, @a vdofs is a list of DOFs. */
652 void EliminateVDofs(const Array<int> &vdofs,
653 DiagonalPolicy dpolicy = DIAG_ONE);
654
655 /** @brief Similar to
656 EliminateVDofs(const Array<int> &, const Vector &,
657 Vector &, DiagonalPolicy)
658 but here @a ess_dofs is a marker (boolean) array on all vector-dofs
659 (@a ess_dofs[i] < 0 is true). */
660 void EliminateEssentialBCFromDofs(const Array<int> &ess_dofs,
661 const Vector &sol,
662 Vector &rhs,
663 DiagonalPolicy dpolicy = DIAG_ONE);
664
665 /** @brief Similar to EliminateVDofs(const Array<int> &, DiagonalPolicy) but
666 here @a ess_dofs is a marker (boolean) array on all vector-dofs
667 (@a ess_dofs[i] < 0 is true). */
668 void EliminateEssentialBCFromDofs(const Array<int> &ess_dofs,
669 DiagonalPolicy dpolicy = DIAG_ONE);
670 /// Perform elimination and set the diagonal entry to the given value
671 void EliminateEssentialBCFromDofsDiag(const Array<int> &ess_dofs,
672 real_t value);
673
674 /** @brief Use the stored eliminated part of the matrix (see
675 EliminateVDofs(const Array<int> &, DiagonalPolicy)) to modify the r.h.s.
676 @a b; @a vdofs is a list of DOFs (non-directional, i.e. >= 0). */
677 void EliminateVDofsInRHS(const Array<int> &vdofs, const Vector &x,
678 Vector &b);
679
680 /** @brief Compute inner product for full uneliminated matrix:
681 $ y^T M x + y^T M_e x $ */
682 real_t FullInnerProduct(const Vector &x, const Vector &y) const
683 { return mat->InnerProduct(x, y) + mat_e->InnerProduct(x, y); }
684
685 /** @brief Update the @a FiniteElementSpace and delete all data associated
686 with the old one. */
687 virtual void Update(FiniteElementSpace *nfes = NULL);
688
689 /// (DEPRECATED) Return the FE space associated with the BilinearForm.
690 /** @deprecated Use FESpace() instead. */
691 MFEM_DEPRECATED FiniteElementSpace *GetFES() { return fes; }
692
693 /// Return the FE space associated with the BilinearForm.
695
696 /// Read-only access to the associated FiniteElementSpace.
697 const FiniteElementSpace *FESpace() const { return fes; }
698
699 /** @brief Sets Operator::DiagonalPolicy used upon construction of the
700 linear system.
701 Policies include:
702
703 - DIAG_ZERO (Set the diagonal values to zero)
704 - DIAG_ONE (Set the diagonal values to one)
705 - DIAG_KEEP (Keep the diagonal values)
706 */
708
709 /// Indicate that integrators are not owned by the BilinearForm
711
712 /** @brief Deletes internal matrices, bilinear integrators, and the
713 BilinearFormExtension */
714 virtual ~BilinearForm();
715};
716
717
718/**
719 Class for assembling of bilinear forms `a(u,v)` defined on different
720 trial and test spaces. The assembled matrix `M` is such that
721
722 a(u,v) = V^t M U
723
724 where `U` and `V` are the vectors representing the functions `u` and `v`,
725 respectively. The first argument, `u`, of `a(,)` is in the trial space
726 and the second argument, `v`, is in the test space. Thus,
727
728 # of rows of M = dimension of the test space and
729 # of cols of M = dimension of the trial space.
730
731 Both trial and test spaces should be defined on the same mesh.
732*/
734{
735protected:
736 SparseMatrix *mat; ///< Owned.
737 SparseMatrix *mat_e; ///< Owned.
738
740 *test_fes; ///< Not owned
741
742 /// The form assembly level (full, partial, etc.)
744
745 /** Extension for supporting Full Assembly (FA), Element Assembly (EA),
746 Partial Assembly (PA), or Matrix Free assembly (MF). */
748
749 /** @brief Indicates the BilinearFormIntegrator%s stored in
750 MixedBilinearForm#domain_integs, MixedBilinearForm#boundary_integs,
751 MixedBilinearForm#trace_face_integs and
752 MixedBilinearForm#boundary_trace_face_integs
753 are owned by another MixedBilinearForm. */
755
756 /// Domain integrators.
758 /// Entries are not owned.
760
761 /// Boundary integrators.
763 /// Entries are not owned.
765
766 /// Trace face (skeleton) integrators.
768
769 /// Boundary trace face (skeleton) integrators.
771 /// Entries are not owned.
773
776
777private:
778 /// Copy construction is not supported; body is undefined.
780
781 /// Copy assignment is not supported; body is undefined.
782 MixedBilinearForm &operator=(const MixedBilinearForm &);
783
784public:
785 /** @brief Construct a MixedBilinearForm on the given trial, @a tr_fes, and
786 test, @a te_fes, FiniteElementSpace%s. */
787 /** The pointers @a tr_fes and @a te_fes are not owned by the newly
788 constructed object. */
790 FiniteElementSpace *te_fes);
791
792 /** @brief Create a MixedBilinearForm on the given trial, @a tr_fes, and
793 test, @a te_fes, FiniteElementSpace%s, using the same integrators as the
794 MixedBilinearForm @a mbf.
795
796 The pointers @a tr_fes and @a te_fes are not owned by the newly
797 constructed object.
798
799 The integrators in @a mbf are copied as pointers and they are not owned
800 by the newly constructed MixedBilinearForm. */
802 FiniteElementSpace *te_fes,
803 MixedBilinearForm *mbf);
804
805 /// Returns a reference to: $ M_{ij} $
806 virtual real_t &Elem(int i, int j);
807
808 /// Returns a reference to: $ M_{ij} $
809 virtual const real_t &Elem(int i, int j) const;
810
811 /// Matrix multiplication: $ y = M x $
812 virtual void Mult(const Vector & x, Vector & y) const;
813
814 /// Add the matrix vector multiple to a vector: $ y += a M x $
815 virtual void AddMult(const Vector & x, Vector & y,
816 const real_t a = 1.0) const;
817
818 /// Matrix transpose vector multiplication: $ y = M^T x $
819 virtual void MultTranspose(const Vector & x, Vector & y) const;
820
821 /// Add the matrix transpose vector multiplication: $ y += a M^T x $
822 virtual void AddMultTranspose(const Vector & x, Vector & y,
823 const real_t a = 1.0) const;
824
825 /** @brief Returns a pointer to (approximation) of the matrix inverse:
826 $ M^{-1} $ (currently unimplemented and returns NULL)*/
827 virtual MatrixInverse *Inverse() const;
828
829 /** @brief Finalizes the matrix initialization if the ::AssemblyLevel is
830 AssemblyLevel::LEGACY.*/
831 virtual void Finalize(int skip_zeros = 1);
832
833 /** @brief Extract the associated matrix as SparseMatrix blocks. The number
834 of block rows and columns is given by the vector dimensions (vdim) of the
835 test and trial spaces, respectively. */
836 void GetBlocks(Array2D<SparseMatrix *> &blocks) const;
837
838 /// Returns a const reference to the sparse matrix: $ M $
839 /** This will segfault if the usual sparse mat is not defined
840 like when static condensation is being used or AllocMat() has
841 not yet been called. */
842 const SparseMatrix &SpMat() const { return *mat; }
843
844 /// Returns a reference to the sparse matrix: $ M $
845 SparseMatrix &SpMat() { return *mat; }
846
847 /** @brief Nullifies the internal matrix $ M $ and returns a pointer
848 to it. Used for transferring ownership. */
849 SparseMatrix *LoseMat() { SparseMatrix *tmp = mat; mat = NULL; return tmp; }
850
851 /// Adds a domain integrator. Assumes ownership of @a bfi.
853
854 /// Adds a domain integrator. Assumes ownership of @a bfi.
856 Array<int> &elem_marker);
857
858 /// Adds a boundary integrator. Assumes ownership of @a bfi.
860
861 /// Adds a boundary integrator. Assumes ownership of @a bfi.
863 Array<int> &bdr_marker);
864
865 /** @brief Add a trace face integrator. Assumes ownership of @a bfi.
866
867 This type of integrator assembles terms over all faces of the mesh using
868 the face FE from the trial space and the two adjacent volume FEs from
869 the test space. */
871
872 /// Adds a boundary trace face integrator. Assumes ownership of @a bfi.
874
875 /// Adds a boundary trace face integrator. Assumes ownership of @a bfi.
877 Array<int> &bdr_marker);
878
879 /// Access all integrators added with AddDomainIntegrator().
881 /** @brief Access all domain markers added with AddDomainIntegrator().
882 If no marker was specified when the integrator was added, the
883 corresponding pointer (to Array<int>) will be NULL. */
885
886 /// Access all integrators added with AddBoundaryIntegrator().
888
889 /** @brief Access all boundary markers added with AddBoundaryIntegrator().
890
891 If no marker was specified when the integrator was added, the
892 corresponding pointer (to Array<int>) will be NULL. */
894
895 /// Access all integrators added with AddTraceFaceIntegrator().
897
898 /// Access all integrators added with AddBdrTraceFaceIntegrator().
900 { return &boundary_trace_face_integs; }
901
902 /** @brief Access all boundary markers added with AddBdrTraceFaceIntegrator()
903
904 If no marker was specified when the integrator was added, the
905 corresponding pointer (to Array<int>) will be NULL. */
908
909 /// Sets all sparse values of $ M $ to @a a.
910 void operator=(const real_t a) { *mat = a; }
911
912 /// Set the desired assembly level. The default is AssemblyLevel::LEGACY.
913 /** This method must be called before assembly. See ::AssemblyLevel*/
914 void SetAssemblyLevel(AssemblyLevel assembly_level);
915
916 void Assemble(int skip_zeros = 1);
917
918 /** @brief Assemble the diagonal of ADA^T into diag, where A is this mixed
919 bilinear form and D is a diagonal. */
920 void AssembleDiagonal_ADAt(const Vector &D, Vector &diag) const;
921
922 /// Get the input finite element space prolongation matrix
923 virtual const Operator *GetProlongation() const
924 { return trial_fes->GetProlongationMatrix(); }
925
926 /// Get the input finite element space restriction matrix
927 virtual const Operator *GetRestriction() const
928 { return trial_fes->GetRestrictionMatrix(); }
929
930 /// Get the test finite element space prolongation matrix
931 virtual const Operator *GetOutputProlongation() const
932 { return test_fes->GetProlongationMatrix(); }
933
934 /// Get the test finite element space restriction matrix
935 virtual const Operator *GetOutputRestriction() const
936 { return test_fes->GetRestrictionMatrix(); }
937
938 /** @brief For partially conforming trial and/or test FE spaces, complete the
939 assembly process by performing $ P2^t A P1 $ where $ A $ is the
940 internal sparse matrix; $ P1 $ and $ P2 $ are the conforming
941 prolongation matrices of the trial and test FE spaces, respectively.
942 After this call the MixedBilinearForm becomes an operator on the
943 conforming FE spaces. */
944 void ConformingAssemble();
945
946 /// Compute the element matrix of the given element
947 void ComputeElementMatrix(int i, DenseMatrix &elmat);
948
949 /// Compute the boundary element matrix of the given boundary element
950 void ComputeBdrElementMatrix(int i, DenseMatrix &elmat);
951
952 /// Assemble the given element matrix
953 /** The element matrix @a elmat is assembled for the element @a i, i.e.
954 added to the system matrix. The flag @a skip_zeros skips the zero
955 elements of the matrix, unless they are breaking the symmetry of
956 the system matrix.
957 */
958 void AssembleElementMatrix(int i, const DenseMatrix &elmat,
959 int skip_zeros = 1);
960
961 /// Assemble the given element matrix
962 /** The element matrix @a elmat is assembled for the element @a i, i.e.
963 added to the system matrix. The vdofs of the element are returned
964 in @a trial_vdofs and @a test_vdofs. The flag @a skip_zeros skips
965 the zero elements of the matrix, unless they are breaking the symmetry
966 of the system matrix.
967 */
968 void AssembleElementMatrix(int i, const DenseMatrix &elmat,
970 int skip_zeros = 1);
971
972 /// Assemble the given boundary element matrix
973 /** The boundary element matrix @a elmat is assembled for the boundary
974 element @a i, i.e. added to the system matrix. The flag @a skip_zeros
975 skips the zero elements of the matrix, unless they are breaking the
976 symmetry of the system matrix.
977 */
978 void AssembleBdrElementMatrix(int i, const DenseMatrix &elmat,
979 int skip_zeros = 1);
980
981 /// Assemble the given boundary element matrix
982 /** The boundary element matrix @a elmat is assembled for the boundary
983 element @a i, i.e. added to the system matrix. The vdofs of the element
984 are returned in @a trial_vdofs and @a test_vdofs. The flag @a skip_zeros
985 skips the zero elements of the matrix, unless they are breaking the
986 symmetry of the system matrix.*/
987 void AssembleBdrElementMatrix(int i, const DenseMatrix &elmat,
990 int skip_zeros = 1);
991
992 /// Eliminate essential boundary DOFs from the columns of the system.
993 /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
994 the essential part of the boundary. All entries in the columns will be
995 set to 0.0 through elimination.*/
996 void EliminateTrialDofs(const Array<int> &bdr_attr_is_ess,
997 const Vector &sol, Vector &rhs);
998
999 /// Eliminate the list of DOFs from the columns of the system.
1000 /** @a marked_vdofs is the of colunm numbers that will be eliminated. All
1001 entries in the columns will be set to 0.0 through elimination.*/
1002 void EliminateEssentialBCFromTrialDofs(const Array<int> &marked_vdofs,
1003 const Vector &sol, Vector &rhs);
1004
1005 /// Eliminate essential boundary DOFs from the rows of the system.
1006 /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
1007 the essential part of the boundary. All entries in the rows will be
1008 set to 0.0 through elimination.*/
1009 virtual void EliminateTestDofs(const Array<int> &bdr_attr_is_ess);
1010
1011 /** @brief Return in @a A that is column-constrained.
1012
1013 This returns the same operator as FormRectangularLinearSystem(), but does
1014 without the transformations of the right-hand side. */
1015 virtual void FormRectangularSystemMatrix(const Array<int> &trial_tdof_list,
1016 const Array<int> &test_tdof_list,
1017 OperatorHandle &A);
1018
1019 /** @brief Form the column-constrained linear system matrix A.
1020
1021 Version of the method FormRectangularSystemMatrix() where the system
1022 matrix is returned in the variable @a A, of type OpType, holding a
1023 *reference* to the system matrix (created with the method
1024 OpType::MakeRef()). The reference will be invalidated when
1025 SetOperatorType(), Update(), or the destructor is called. */
1026 template <typename OpType>
1027 void FormRectangularSystemMatrix(const Array<int> &trial_tdof_list,
1028 const Array<int> &test_tdof_list, OpType &A)
1029 {
1030 OperatorHandle Ah;
1031 FormRectangularSystemMatrix(trial_tdof_list, test_tdof_list, Ah);
1032 OpType *A_ptr = Ah.Is<OpType>();
1033 MFEM_VERIFY(A_ptr, "invalid OpType used");
1034 A.MakeRef(*A_ptr);
1035 }
1036
1037 /** @brief Form the linear system A X = B, corresponding to this mixed
1038 bilinear form and the linear form @a b(.).
1039
1040 Return in @a A a *reference* to the system matrix that is
1041 column-constrained. The reference will be invalidated when
1042 SetOperatorType(), Update(), or the destructor is called. */
1043 virtual void FormRectangularLinearSystem(const Array<int> &trial_tdof_list,
1044 const Array<int> &test_tdof_list,
1045 Vector &x, Vector &b,
1046 OperatorHandle &A, Vector &X,
1047 Vector &B);
1048
1049 /** @brief Form the linear system A X = B, corresponding to this bilinear
1050 form and the linear form @a b(.).
1051
1052 Version of the method FormRectangularLinearSystem() where the system
1053 matrix is returned in the variable @a A, of type OpType, holding a
1054 *reference* to the system matrix (created with the method
1055 OpType::MakeRef()). The reference will be invalidated when
1056 SetOperatorType(), Update(), or the destructor is called. */
1057 template <typename OpType>
1058 void FormRectangularLinearSystem(const Array<int> &trial_tdof_list,
1059 const Array<int> &test_tdof_list,
1060 Vector &x, Vector &b,
1061 OpType &A, Vector &X, Vector &B)
1062 {
1063 OperatorHandle Ah;
1064 FormRectangularLinearSystem(trial_tdof_list, test_tdof_list, x, b,
1065 Ah, X, B);
1066 OpType *A_ptr = Ah.Is<OpType>();
1067 MFEM_VERIFY(A_ptr, "invalid OpType used");
1068 A.MakeRef(*A_ptr);
1069 }
1070
1071 /// Must be called after making changes to #trial_fes or #test_fes.
1072 void Update();
1073
1074 /// Return the trial FE space associated with the BilinearForm.
1076
1077 /// Read-only access to the associated trial FiniteElementSpace.
1078 const FiniteElementSpace *TrialFESpace() const { return trial_fes; }
1079
1080 /// Return the test FE space associated with the BilinearForm.
1082
1083 /// Read-only access to the associated test FiniteElementSpace.
1084 const FiniteElementSpace *TestFESpace() const { return test_fes; }
1085
1086 /** @brief Deletes internal matrices, bilinear integrators, and the
1087 BilinearFormExtension */
1088 virtual ~MixedBilinearForm();
1089};
1090
1091
1092/**
1093 Class for constructing the matrix representation of a linear operator,
1094 `v = L u`, from one FiniteElementSpace (domain) to another FiniteElementSpace
1095 (range). The constructed matrix `A` is such that
1096
1097 V = A U
1098
1099 where `U` and `V` are the vectors of degrees of freedom representing the
1100 functions `u` and `v`, respectively. The dimensions of `A` are
1101
1102 number of rows of A = dimension of the range space and
1103 number of cols of A = dimension of the domain space.
1104
1105 This class is very similar to MixedBilinearForm. One difference is that
1106 the linear operator `L` is defined using a special kind of
1107 BilinearFormIntegrator (we reuse its functionality instead of defining a
1108 new class). The other difference with the MixedBilinearForm class is that
1109 the "assembly" process overwrites the global matrix entries using the
1110 local element matrices instead of adding them.
1111
1112 Note that if we define the bilinear form `b(u,v) := (Lu,v)` using an inner
1113 product in the range space, then its matrix representation, `B`, is
1114
1115 B = M A, (since V^t B U = b(u,v) = (Lu,v) = V^t M A U)
1116
1117 where `M` denotes the mass matrix for the inner product in the range space:
1118 `V1^t M V2 = (v1,v2)`. Similarly, if `c(u,w) := (Lu,Lw)` then
1119
1120 C = A^t M A.
1121*/
1123{
1124private:
1125 /// Copy construction is not supported; body is undefined.
1127
1128 /// Copy assignment is not supported; body is undefined.
1130
1131public:
1132 /** @brief Construct a DiscreteLinearOperator on the given
1133 FiniteElementSpace%s @a domain_fes and @a range_fes. */
1134 /** The pointers @a domain_fes and @a range_fes are not owned by the newly
1135 constructed object. */
1137 FiniteElementSpace *range_fes)
1138 : MixedBilinearForm(domain_fes, range_fes) { }
1139
1140 /// Adds a domain interpolator. Assumes ownership of @a di.
1142 { AddDomainIntegrator(di); }
1144 Array<int> &elem_marker)
1145 { AddDomainIntegrator(di, elem_marker); }
1146
1147 /// Adds a trace face interpolator. Assumes ownership of @a di.
1149 { AddTraceFaceIntegrator(di); }
1150
1151 /// Access all interpolators added with AddDomainInterpolator().
1154
1155 /// Set the desired assembly level. The default is AssemblyLevel::FULL.
1156 /** This method must be called before assembly. */
1157 void SetAssemblyLevel(AssemblyLevel assembly_level);
1158
1159 /** @brief Construct the internal matrix representation of the discrete
1160 linear operator. */
1161 virtual void Assemble(int skip_zeros = 1);
1162
1163 /** @brief Get the output finite element space restriction matrix in
1164 transposed form. */
1167};
1168
1169}
1170
1171#endif
Dynamic 2D array using row-major layout.
Definition: array.hpp:372
Class extending the BilinearForm class to support different AssemblyLevels.
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:27
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level.
void AllocateMatrix()
Pre-allocate the internal SparseMatrix before assembly. If the internal flag precompute_sparsity is s...
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization if the AssemblyLevel is AssemblyLevel::LEGACY....
virtual const Operator * GetProlongation() const
Get the finite element space prolongation operator.
void SetDiagonalPolicy(DiagonalPolicy policy)
Sets Operator::DiagonalPolicy used upon construction of the linear system. Policies include:
Array< BilinearFormIntegrator * > domain_integs
Set of Domain Integrators to be applied.
Array< Array< int > * > domain_integs_marker
Entries are not owned.
BilinearFormExtension * ext
Extension for supporting Full Assembly (FA), Element Assembly (EA),Partial Assembly (PA),...
void EnableStaticCondensation()
Enable the use of static condensation. For details see the description for class StaticCondensation i...
virtual void AssembleDiagonal(Vector &diag) const
Assemble the diagonal of the bilinear form into diag. Note that diag is a tdof Vector.
void UsePrecomputedSparsity(int ps=1)
For scalar FE spaces, precompute the sparsity pattern of the matrix (assuming dense element matrices)...
void AssembleBdrElementMatrix(int i, const DenseMatrix &elmat, int skip_zeros=1)
Assemble the given boundary element matrix.
Array< Array< int > * > boundary_face_integs_marker
Entries are not owned.
bool HasSpMatElim()
Returns true if the sparse matrix of eliminated b.c.s is not null, false otherwise.
void UseSparsity(int *I, int *J, bool isSorted)
Use the given CSR sparsity pattern to allocate the internal SparseMatrix.
void EliminateVDofsInRHS(const Array< int > &vdofs, const Vector &x, Vector &b)
Use the stored eliminated part of the matrix (see EliminateVDofs(const Array<int> &,...
void FullAddMult(const Vector &x, Vector &y) const
Add the original uneliminated matrix vector multiple to a vector. The original matrix is so we have:...
void FullMult(const Vector &x, Vector &y) const
Matrix vector multiplication with the original uneliminated matrix. The original matrix is so we hav...
virtual void MultTranspose(const Vector &x, Vector &y) const
Matrix transpose vector multiplication: .
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication: .
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void ComputeElementMatrices()
Compute and store internally all element matrices.
Array< BilinearFormIntegrator * > boundary_face_integs
Set of boundary face Integrators to be applied.
Array< BilinearFormIntegrator * > * GetFBFI()
Access all integrators added with AddInteriorFaceIntegrator().
virtual void Update(FiniteElementSpace *nfes=NULL)
Update the FiniteElementSpace and delete all data associated with the old one.
void EliminateVDofs(const Array< int > &vdofs, const Vector &sol, Vector &rhs, DiagonalPolicy dpolicy=DIAG_ONE)
Eliminate the given vdofs. NOTE: here, vdofs is a list of DOFs.
Array< BilinearFormIntegrator * > * GetDBFI()
Access all the integrators added with AddDomainIntegrator().
void FreeElementMatrices()
Free the memory used by the element matrices.
virtual ~BilinearForm()
Deletes internal matrices, bilinear integrators, and the BilinearFormExtension.
BilinearForm()
may be used in the construction of derived classes
Hybridization * hybridization
Owned.
DenseTensor * element_matrices
Owned.
void EliminateEssentialBCFromDofsDiag(const Array< int > &ess_dofs, real_t value)
Perform elimination and set the diagonal entry to the given value.
DiagonalPolicy diag_policy
This data member allows one to specify what should be done to the diagonal matrix entries and corresp...
void EnableSparseMatrixSorting(bool enable_it)
Force the sparse matrix column indices to be sorted when using AssemblyLevel::FULL.
FiniteElementSpace * FESpace()
Return the FE space associated with the BilinearForm.
void AddBoundaryIntegrator(BilinearFormIntegrator *bfi)
Adds new Boundary Integrator. Assumes ownership of bfi.
const FiniteElementSpace * FESpace() const
Read-only access to the associated FiniteElementSpace.
void FullAddMultTranspose(const Vector &x, Vector &y) const
Add the original uneliminated matrix transpose vector multiple to a vector. The original matrix is s...
void AllocMat()
Allocate appropriate SparseMatrix and assign it to mat.
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
AssemblyLevel GetAssemblyLevel() const
Returns the assembly level.
virtual const Operator * GetOutputRestriction() const
Get the output finite element space restriction matrix.
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Recover the solution of a linear system formed with FormLinearSystem().
virtual MatrixInverse * Inverse() const
Returns a pointer to (approximation) of the matrix inverse: (currently returns NULL)
real_t FullInnerProduct(const Vector &x, const Vector &y) const
Compute inner product for full uneliminated matrix: .
real_t InnerProduct(const Vector &x, const Vector &y) const
Compute .
const real_t & operator()(int i, int j)
Returns a reference to: .
void operator=(const real_t a)
Sets all sparse values of and to 'a'.
SparseMatrix & SpMatElim()
Returns a reference to the sparse matrix of eliminated b.c.: .
SparseMatrix * LoseMat()
Nullifies the internal matrix and returns a pointer to it. Used for transferring ownership.
virtual const Operator * GetRestriction() const
Get the finite element space restriction operator.
bool HasSpMat()
Returns true if the sparse matrix is not null, false otherwise.
void EnableHybridization(FiniteElementSpace *constr_space, BilinearFormIntegrator *constr_integ, const Array< int > &ess_tdof_list)
Enable hybridization.
int extern_bfs
Indicates the BilinearFormIntegrators stored in domain_integs, boundary_integs, interior_face_integs,...
void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new boundary Face Integrator. Assumes ownership of bfi.
const SparseMatrix & SpMatElim() const
Returns a const reference to the sparse matrix of eliminated b.c.: .
void AssembleElementMatrix(int i, const DenseMatrix &elmat, int skip_zeros=1)
Assemble the given element matrix.
void FormSystemMatrix(const Array< int > &ess_tdof_list, OpType &A)
Form the linear system matrix A, see FormLinearSystem() for details.
SparseMatrix & SpMat()
Returns a reference to the sparse matrix: .
Array< BilinearFormIntegrator * > * GetBBFI()
Access all the integrators added with AddBoundaryIntegrator().
Array< Array< int > * > boundary_integs_marker
Entries are not owned.
virtual real_t & Elem(int i, int j)
Returns a reference to: .
long sequence
Indicates the Mesh::sequence corresponding to the current state of the BilinearForm.
Array< BilinearFormIntegrator * > interior_face_integs
Set of interior face Integrators to be applied.
virtual void AddMultTranspose(const Vector &x, Vector &y, const real_t a=1.0) const
Add the matrix transpose vector multiplication: .
void ComputeBdrElementMatrix(int i, DenseMatrix &elmat)
Compute the boundary element matrix of the given boundary element.
Array< Array< int > * > * GetBBFI_Marker()
Access all boundary markers added with AddBoundaryIntegrator(). If no marker was specified when the i...
void EliminateEssentialBCFromDofs(const Array< int > &ess_dofs, const Vector &sol, Vector &rhs, DiagonalPolicy dpolicy=DIAG_ONE)
Similar to EliminateVDofs(const Array<int> &, const Vector &, Vector &,...
SparseMatrix * mat
Sparse matrix to be associated with the form. Owned.
AssemblyLevel assembly
The AssemblyLevel of the form (AssemblyLevel::LEGACY, AssemblyLevel::FULL, AssemblyLevel::ELEMENT,...
void SerialRAP(OperatorHandle &A)
Compute serial RAP operator and store it in A as a SparseMatrix.
FiniteElementSpace * SCFESpace() const
Return the trace FE space associated with static condensation.
virtual const Operator * GetOutputProlongation() const
Get the output finite element space prolongation matrix.
Array< BilinearFormIntegrator * > * GetBFBFI()
Access all integrators added with AddBdrFaceIntegrator().
FiniteElementSpace * fes
FE space on which the form lives. Not owned.
int batch
Element batch size used in the form action (1, 8, num_elems, etc.)
virtual void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
Form the linear system A X = B, corresponding to this bilinear form and the linear form b(....
bool StaticCondensationIsEnabled() const
Check if static condensation was actually enabled by a previous call to EnableStaticCondensation().
void UseExternalIntegrators()
Indicate that integrators are not owned by the BilinearForm.
int Size() const
Get the size of the BilinearForm as a square matrix.
MFEM_DEPRECATED FiniteElementSpace * GetFES()
(DEPRECATED) Return the FE space associated with the BilinearForm.
virtual void AddMult(const Vector &x, Vector &y, const real_t a=1.0) const
Add the matrix vector multiple to a vector: .
Hybridization * GetHybridization() const
Array< Array< int > * > * GetBFBFI_Marker()
Access all boundary markers added with AddBdrFaceIntegrator(). If no marker was specified when the in...
void ComputeElementMatrix(int i, DenseMatrix &elmat)
Compute the element matrix of the given element.
Array< int > vdofs
void EliminateEssentialBC(const Array< int > &bdr_attr_is_ess, const Vector &sol, Vector &rhs, DiagonalPolicy dpolicy=DIAG_ONE)
Eliminate essential boundary DOFs from the system.
void ConformingAssemble()
For partially conforming trial and/or test FE spaces, complete the assembly process by performing wh...
Array< BilinearFormIntegrator * > boundary_integs
Set of Boundary Integrators to be applied.
virtual const Operator * GetOutputRestrictionTranspose() const
Returns the output fe space restriction matrix, transposed.
void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OpType &A, Vector &X, Vector &B, int copy_interior=0)
Form the linear system A X = B, corresponding to this bilinear form and the linear form b(....
void EliminateEssentialBCDiag(const Array< int > &bdr_attr_is_ess, real_t value)
Perform elimination and set the diagonal entry to the given value.
StaticCondensation * static_cond
Owned.
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
Array< Array< int > * > * GetDBFI_Marker()
Access all boundary markers added with AddDomainIntegrator(). If no marker was specified when the int...
const SparseMatrix & SpMat() const
Returns a const reference to the sparse matrix: .
SparseMatrix * mat_e
Sparse Matrix used to store the eliminations from the b.c. Owned. .
void AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new interior Face Integrator. Assumes ownership of bfi.
Data type dense matrix using column-major storage.
Definition: densemat.hpp:24
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:1104
void AddDomainInterpolator(DiscreteInterpolator *di, Array< int > &elem_marker)
void AddDomainInterpolator(DiscreteInterpolator *di)
Adds a domain interpolator. Assumes ownership of di.
void AddTraceFaceInterpolator(DiscreteInterpolator *di)
Adds a trace face interpolator. Assumes ownership of di.
virtual void Assemble(int skip_zeros=1)
Construct the internal matrix representation of the discrete linear operator.
Array< Array< int > * > * GetDI_Marker()
DiscreteLinearOperator(FiniteElementSpace *domain_fes, FiniteElementSpace *range_fes)
Construct a DiscreteLinearOperator on the given FiniteElementSpaces domain_fes and range_fes.
virtual const Operator * GetOutputRestrictionTranspose() const
Get the output finite element space restriction matrix in transposed form.
Array< BilinearFormIntegrator * > * GetDI()
Access all interpolators added with AddDomainInterpolator().
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level. The default is AssemblyLevel::FULL.
Data and methods for fully-assembled bilinear forms.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:220
const SparseMatrix * GetConformingRestriction() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.cpp:1309
virtual const SparseMatrix * GetRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.hpp:620
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
Definition: fespace.hpp:597
const Operator * GetRestrictionTransposeOperator() const
Return an operator that performs the transpose of GetRestrictionOperator.
Definition: fespace.cpp:1324
const SparseMatrix * GetConformingProlongation() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.cpp:1302
Auxiliary class Hybridization, used to implement BilinearForm hybridization.
Abstract data type for matrix inverse.
Definition: matrix.hpp:63
Abstract data type matrix.
Definition: matrix.hpp:28
Class extending the MixedBilinearForm class to support different AssemblyLevels.
virtual const Operator * GetOutputRestriction() const
Get the test finite element space restriction matrix.
virtual void MultTranspose(const Vector &x, Vector &y) const
Matrix transpose vector multiplication: .
const FiniteElementSpace * TrialFESpace() const
Read-only access to the associated trial FiniteElementSpace.
void ConformingAssemble()
For partially conforming trial and/or test FE spaces, complete the assembly process by performing wh...
MixedBilinearFormExtension * ext
void AssembleBdrElementMatrix(int i, const DenseMatrix &elmat, int skip_zeros=1)
Assemble the given boundary element matrix.
const FiniteElementSpace * TestFESpace() const
Read-only access to the associated test FiniteElementSpace.
SparseMatrix & SpMat()
Returns a reference to the sparse matrix: .
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization if the AssemblyLevel is AssemblyLevel::LEGACY.
void Assemble(int skip_zeros=1)
virtual const Operator * GetProlongation() const
Get the input finite element space prolongation matrix.
void AddBoundaryIntegrator(BilinearFormIntegrator *bfi)
Adds a boundary integrator. Assumes ownership of bfi.
virtual void EliminateTestDofs(const Array< int > &bdr_attr_is_ess)
Eliminate essential boundary DOFs from the rows of the system.
virtual void FormRectangularSystemMatrix(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, OperatorHandle &A)
Return in A that is column-constrained.
Array< BilinearFormIntegrator * > * GetTFBFI()
Access all integrators added with AddTraceFaceIntegrator().
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level. The default is AssemblyLevel::LEGACY.
int extern_bfs
Indicates the BilinearFormIntegrators stored in MixedBilinearForm::domain_integs, MixedBilinearForm::...
void AssembleDiagonal_ADAt(const Vector &D, Vector &diag) const
Assemble the diagonal of ADA^T into diag, where A is this mixed bilinear form and D is a diagonal.
Array< BilinearFormIntegrator * > * GetBBFI()
Access all integrators added with AddBoundaryIntegrator().
void AddBdrTraceFaceIntegrator(BilinearFormIntegrator *bfi)
Adds a boundary trace face integrator. Assumes ownership of bfi.
Array< Array< int > * > boundary_integs_marker
Entries are not owned.
Array< BilinearFormIntegrator * > trace_face_integs
Trace face (skeleton) integrators.
virtual const Operator * GetRestriction() const
Get the input finite element space restriction matrix.
FiniteElementSpace * TestFESpace()
Return the test FE space associated with the BilinearForm.
Array< BilinearFormIntegrator * > boundary_trace_face_integs
Boundary trace face (skeleton) integrators.
virtual void Mult(const Vector &x, Vector &y) const
Matrix multiplication: .
Array< BilinearFormIntegrator * > domain_integs
Domain integrators.
SparseMatrix * LoseMat()
Nullifies the internal matrix and returns a pointer to it. Used for transferring ownership.
SparseMatrix * mat
Owned.
void Update()
Must be called after making changes to trial_fes or test_fes.
Array< Array< int > * > * GetBTFBFI_Marker()
Access all boundary markers added with AddBdrTraceFaceIntegrator()
virtual void AddMult(const Vector &x, Vector &y, const real_t a=1.0) const
Add the matrix vector multiple to a vector: .
void ComputeElementMatrix(int i, DenseMatrix &elmat)
Compute the element matrix of the given element.
Array< Array< int > * > domain_integs_marker
Entries are not owned.
Array< BilinearFormIntegrator * > boundary_integs
Boundary integrators.
void EliminateEssentialBCFromTrialDofs(const Array< int > &marked_vdofs, const Vector &sol, Vector &rhs)
Eliminate the list of DOFs from the columns of the system.
void EliminateTrialDofs(const Array< int > &bdr_attr_is_ess, const Vector &sol, Vector &rhs)
Eliminate essential boundary DOFs from the columns of the system.
void AssembleElementMatrix(int i, const DenseMatrix &elmat, int skip_zeros=1)
Assemble the given element matrix.
Array< Array< int > * > boundary_trace_face_integs_marker
Entries are not owned.
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)
Form the linear system A X = B, corresponding to this mixed bilinear form and the linear form b(....
const SparseMatrix & SpMat() const
Returns a const reference to the sparse matrix: .
virtual real_t & Elem(int i, int j)
Returns a reference to: .
void ComputeBdrElementMatrix(int i, DenseMatrix &elmat)
Compute the boundary element matrix of the given boundary element.
void FormRectangularSystemMatrix(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, OpType &A)
Form the column-constrained linear system matrix A.
SparseMatrix * mat_e
Owned.
FiniteElementSpace * trial_fes
Not owned.
void GetBlocks(Array2D< SparseMatrix * > &blocks) const
Extract the associated matrix as SparseMatrix blocks. The number of block rows and columns is given b...
FiniteElementSpace * TrialFESpace()
Return the trial FE space associated with the BilinearForm.
Array< BilinearFormIntegrator * > * GetDBFI()
Access all integrators added with AddDomainIntegrator().
Array< Array< int > * > * GetBBFI_Marker()
Access all boundary markers added with AddBoundaryIntegrator().
void AddTraceFaceIntegrator(BilinearFormIntegrator *bfi)
Add a trace face integrator. Assumes ownership of bfi.
void operator=(const real_t a)
Sets all sparse values of to a.
virtual void AddMultTranspose(const Vector &x, Vector &y, const real_t a=1.0) const
Add the matrix transpose vector multiplication: .
Array< Array< int > * > * GetDBFI_Marker()
Access all domain markers added with AddDomainIntegrator(). If no marker was specified when the integ...
Array< BilinearFormIntegrator * > * GetBTFBFI()
Access all integrators added with AddBdrTraceFaceIntegrator().
virtual ~MixedBilinearForm()
Deletes internal matrices, bilinear integrators, and the BilinearFormExtension.
FiniteElementSpace * test_fes
Not owned.
virtual MatrixInverse * Inverse() const
Returns a pointer to (approximation) of the matrix inverse: (currently unimplemented and returns NUL...
virtual const Operator * GetOutputProlongation() const
Get the test finite element space prolongation matrix.
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds a domain integrator. Assumes ownership of bfi.
AssemblyLevel assembly
The form assembly level (full, partial, etc.)
void FormRectangularLinearSystem(const Array< int > &trial_tdof_list, const Array< int > &test_tdof_list, Vector &x, Vector &b, OpType &A, Vector &X, Vector &B)
Form the linear system A X = B, corresponding to this bilinear form and the linear form b(....
Pointer to an Operator of a specified type.
Definition: handle.hpp:34
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
Definition: handle.hpp:145
OpType * Is() const
Return the Operator pointer dynamically cast to a specified OpType.
Definition: handle.hpp:108
Abstract operator.
Definition: operator.hpp:25
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
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
@ DIAG_KEEP
Keep the diagonal value.
Definition: operator.hpp:51
Data type sparse matrix.
Definition: sparsemat.hpp:51
virtual void AddMultTranspose(const Vector &x, Vector &y, const real_t a=1.0) const
y += At * x (default) or y += a * At * x
Definition: sparsemat.cpp:961
real_t InnerProduct(const Vector &x, const Vector &y) const
Compute y^t A x.
Definition: sparsemat.cpp:1227
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:755
virtual void AddMult(const Vector &x, Vector &y, const real_t a=1.0) const
y += A * x (default) or y += a * A * x
Definition: sparsemat.cpp:762
FiniteElementSpace * GetTraceFESpace()
Return a pointer to the reduced/trace FE space.
Definition: staticcond.hpp:109
Vector data type.
Definition: vector.hpp:80
real_t b
Definition: lissajous.cpp:42
real_t a
Definition: lissajous.cpp:41
AssemblyLevel
Enumeration defining the assembly level for bilinear and nonlinear form classes derived from Operator...
float real_t
Definition: config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition: lor_mms.hpp:30