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