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