1 // Copyright (c) 2010-2021, 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 { 61 friend FABilinearFormExtension; 62 63 protected: 64 /// Sparse matrix \f$ M \f$ to be associated with the form. Owned. 65 SparseMatrix *mat; 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$ */ 70 SparseMatrix *mat_e; 71 72 /// FE space on which the form lives. Not owned. 73 FiniteElementSpace *fes; 74 75 /// The assembly level of the form (full, partial, etc.) 76 AssemblyLevel assembly; 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). */ 81 BilinearFormExtension *ext; 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. */ 90 int extern_bfs; 91 92 /// Set of Domain Integrators to be applied. 93 Array<BilinearFormIntegrator*> domain_integs; 94 /// Element attribute marker (should be of length mesh->attributes) 95 /// Includes all by default. 96 /// 0 - ignore attribute 97 /// 1 - include attribute 98 Array<Array<int>*> domain_integs_marker; 99 100 /// Set of Boundary Integrators to be applied. 101 Array<BilinearFormIntegrator*> boundary_integs; 102 Array<Array<int>*> boundary_integs_marker; ///< Entries are not owned. 103 104 /// Set of interior face Integrators to be applied. 105 Array<BilinearFormIntegrator*> interior_face_integs; 106 107 /// Set of boundary face Integrators to be applied. 108 Array<BilinearFormIntegrator*> boundary_face_integs; 109 Array<Array<int>*> boundary_face_integs_marker; ///< Entries are not owned. 110 111 DenseMatrix elemmat; 112 Array<int> vdofs; 113 114 DenseTensor *element_matrices; ///< Owned. 115 116 StaticCondensation *static_cond; ///< Owned. 117 Hybridization *hybridization; ///< Owned. 118 119 /** This data member allows one to specify what should be done to the 120 diagonal matrix entries and corresponding RHS values upon elimination of 121 the constrained DoFs. */ 122 DiagonalPolicy diag_policy; 123 124 int precompute_sparsity; 125 // Allocate appropriate SparseMatrix and assign it to mat 126 void AllocMat(); 127 128 void ConformingAssemble(); 129 130 // may be used in the construction of derived classes BilinearForm()131 BilinearForm() : Matrix (0) 132 { 133 fes = NULL; sequence = -1; 134 mat = mat_e = NULL; extern_bfs = 0; element_matrices = NULL; 135 static_cond = NULL; hybridization = NULL; 136 precompute_sparsity = 0; 137 diag_policy = DIAG_KEEP; 138 assembly = AssemblyLevel::LEGACY; 139 batch = 1; 140 ext = NULL; 141 } 142 143 private: 144 /// Copy construction is not supported; body is undefined. 145 BilinearForm(const BilinearForm &); 146 147 /// Copy assignment is not supported; body is undefined. 148 BilinearForm &operator=(const BilinearForm &); 149 150 public: 151 /// Creates bilinear form associated with FE space @a *f. 152 /** The pointer @a f is not owned by the newly constructed object. */ 153 BilinearForm(FiniteElementSpace *f); 154 155 /** @brief Create a BilinearForm on the FiniteElementSpace @a f, using the 156 same integrators as the BilinearForm @a bf. 157 158 The pointer @a f is not owned by the newly constructed object. 159 160 The integrators in @a bf are copied as pointers and they are not owned by 161 the newly constructed BilinearForm. 162 163 The optional parameter @a ps is used to initialize the internal flag 164 #precompute_sparsity, see UsePrecomputedSparsity() for details. */ 165 BilinearForm(FiniteElementSpace *f, BilinearForm *bf, int ps = 0); 166 167 /// Get the size of the BilinearForm as a square matrix. Size() const168 int Size() const { return height; } 169 170 /// Set the desired assembly level. 171 /** Valid choices are: 172 173 - AssemblyLevel::LEGACY (default) 174 - AssemblyLevel::FULL 175 - AssemblyLevel::PARTIAL 176 - AssemblyLevel::ELEMENT 177 - AssemblyLevel::NONE 178 179 This method must be called before assembly. */ 180 void SetAssemblyLevel(AssemblyLevel assembly_level); 181 182 /// Returns the assembly level GetAssemblyLevel() const183 AssemblyLevel GetAssemblyLevel() const { return assembly; } 184 GetHybridization() const185 Hybridization *GetHybridization() const { return hybridization; } 186 187 /** @brief Enable the use of static condensation. For details see the 188 description for class StaticCondensation in fem/staticcond.hpp This method 189 should be called before assembly. If the number of unknowns after static 190 condensation is not reduced, it is not enabled. */ 191 void EnableStaticCondensation(); 192 193 /** @brief Check if static condensation was actually enabled by a previous 194 call to EnableStaticCondensation(). */ StaticCondensationIsEnabled() const195 bool StaticCondensationIsEnabled() const { return static_cond; } 196 197 /// Return the trace FE space associated with static condensation. SCFESpace() const198 FiniteElementSpace *SCFESpace() const 199 { return static_cond ? static_cond->GetTraceFESpace() : NULL; } 200 201 /// Enable hybridization. 202 /** For details see the description for class 203 Hybridization in fem/hybridization.hpp. This method should be called 204 before assembly. */ 205 void EnableHybridization(FiniteElementSpace *constr_space, 206 BilinearFormIntegrator *constr_integ, 207 const Array<int> &ess_tdof_list); 208 209 /** @brief For scalar FE spaces, precompute the sparsity pattern of the matrix 210 (assuming dense element matrices) based on the types of integrators 211 present in the bilinear form. */ UsePrecomputedSparsity(int ps=1)212 void UsePrecomputedSparsity(int ps = 1) { precompute_sparsity = ps; } 213 214 /** @brief Use the given CSR sparsity pattern to allocate the internal 215 SparseMatrix. 216 217 - The @a I and @a J arrays must define a square graph with size equal to 218 GetVSize() of the associated FiniteElementSpace. 219 - This method should be called after enabling static condensation or 220 hybridization, if used. 221 - In the case of static condensation, @a I and @a J are not used. 222 - The ownership of the arrays @a I and @a J remains with the caller. */ 223 void UseSparsity(int *I, int *J, bool isSorted); 224 225 /// Use the sparsity of @a A to allocate the internal SparseMatrix. 226 void UseSparsity(SparseMatrix &A); 227 228 /// Pre-allocate the internal SparseMatrix before assembly. 229 /** If the flag 'precompute sparsity' 230 is set, the matrix is allocated in CSR format (i.e. 231 finalized) and the entries are initialized with zeros. */ AllocateMatrix()232 void AllocateMatrix() { if (mat == NULL) { AllocMat(); } } 233 234 /// Access all the integrators added with AddDomainIntegrator(). GetDBFI()235 Array<BilinearFormIntegrator*> *GetDBFI() { return &domain_integs; } 236 237 /// Access all the integrators added with AddBoundaryIntegrator(). GetBBFI()238 Array<BilinearFormIntegrator*> *GetBBFI() { return &boundary_integs; } 239 /** @brief Access all boundary markers added with AddBoundaryIntegrator(). 240 If no marker was specified when the integrator was added, the 241 corresponding pointer (to Array<int>) will be NULL. */ GetBBFI_Marker()242 Array<Array<int>*> *GetBBFI_Marker() { return &boundary_integs_marker; } 243 244 /// Access all integrators added with AddInteriorFaceIntegrator(). GetFBFI()245 Array<BilinearFormIntegrator*> *GetFBFI() { return &interior_face_integs; } 246 247 /// Access all integrators added with AddBdrFaceIntegrator(). GetBFBFI()248 Array<BilinearFormIntegrator*> *GetBFBFI() { return &boundary_face_integs; } 249 /** @brief Access all boundary markers added with AddBdrFaceIntegrator(). 250 If no marker was specified when the integrator was added, the 251 corresponding pointer (to Array<int>) will be NULL. */ GetBFBFI_Marker()252 Array<Array<int>*> *GetBFBFI_Marker() 253 { return &boundary_face_integs_marker; } 254 255 /// Returns a reference to: \f$ M_{ij} \f$ operator ()(int i,int j)256 const double &operator()(int i, int j) { return (*mat)(i,j); } 257 258 /// Returns a reference to: \f$ M_{ij} \f$ 259 virtual double &Elem(int i, int j); 260 261 /// Returns constant reference to: \f$ M_{ij} \f$ 262 virtual const double &Elem(int i, int j) const; 263 264 /// Matrix vector multiplication: \f$ y = M x \f$ 265 virtual void Mult(const Vector &x, Vector &y) const; 266 267 /** @brief Matrix vector multiplication with the original uneliminated 268 matrix. The original matrix is \f$ M + M_e \f$ so we have: 269 \f$ y = M x + M_e x \f$ */ FullMult(const Vector & x,Vector & y) const270 void FullMult(const Vector &x, Vector &y) const 271 { mat->Mult(x, y); mat_e->AddMult(x, y); } 272 273 /// Add the matrix vector multiple to a vector: \f$ y += a M x \f$ AddMult(const Vector & x,Vector & y,const double a=1.0) const274 virtual void AddMult(const Vector &x, Vector &y, const double a = 1.0) const 275 { mat -> AddMult (x, y, a); } 276 277 /** @brief Add the original uneliminated matrix vector multiple to a vector. 278 The original matrix is \f$ M + Me \f$ so we have: 279 \f$ y += M x + M_e x \f$ */ FullAddMult(const Vector & x,Vector & y) const280 void FullAddMult(const Vector &x, Vector &y) const 281 { mat->AddMult(x, y); mat_e->AddMult(x, y); } 282 283 /// Add the matrix transpose vector multiplication: \f$ y += a M^T x \f$ AddMultTranspose(const Vector & x,Vector & y,const double a=1.0) const284 virtual void AddMultTranspose(const Vector & x, Vector & y, 285 const double a = 1.0) const 286 { mat->AddMultTranspose(x, y, a); } 287 288 /** @brief Add the original uneliminated matrix transpose vector 289 multiple to a vector. The original matrix is \f$ M + M_e \f$ 290 so we have: \f$ y += M^T x + {M_e}^T x \f$ */ FullAddMultTranspose(const Vector & x,Vector & y) const291 void FullAddMultTranspose(const Vector & x, Vector & y) const 292 { mat->AddMultTranspose(x, y); mat_e->AddMultTranspose(x, y); } 293 294 /// Matrix transpose vector multiplication: \f$ y = M^T x \f$ MultTranspose(const Vector & x,Vector & y) const295 virtual void MultTranspose(const Vector & x, Vector & y) const 296 { y = 0.0; AddMultTranspose (x, y); } 297 298 /// Compute \f$ y^T M x \f$ InnerProduct(const Vector & x,const Vector & y) const299 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. SpMat() const309 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$ SpMat()316 SparseMatrix &SpMat() 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. */ LoseMat()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$ SpMatElim() const327 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$ SpMatElim()334 SparseMatrix &SpMatElim() 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. 341 void AddDomainIntegrator(BilinearFormIntegrator *bfi); 342 /// Adds new Domain Integrator restricted to certain elements specified by 343 /// the @a elem_attr_marker. 344 void AddDomainIntegrator(BilinearFormIntegrator *bfi, 345 Array<int> &elem_marker); 346 347 /// Adds new Boundary Integrator. Assumes ownership of @a bfi. 348 void AddBoundaryIntegrator(BilinearFormIntegrator *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. */ 355 void AddBoundaryIntegrator(BilinearFormIntegrator *bfi, 356 Array<int> &bdr_marker); 357 358 /// Adds new interior Face Integrator. Assumes ownership of @a bfi. 359 void AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi); 360 361 /// Adds new boundary Face Integrator. Assumes ownership of @a bfi. 362 void AddBdrFaceIntegrator(BilinearFormIntegrator *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. */ 369 void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi, 370 Array<int> &bdr_marker); 371 372 /// Sets all sparse values of \f$ M \f$ and \f$ M_e \f$ to 'a'. operator =(const double 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. GetProlongation() const394 virtual const Operator *GetProlongation() const 395 { return fes->GetConformingProlongation(); } 396 /// Get the finite element space restriction operator GetRestriction() const397 virtual const Operator *GetRestriction() const 398 { return fes->GetConformingRestriction(); } 399 /// Get the output finite element space prolongation matrix GetOutputProlongation() const400 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. */ GetOutputRestrictionTranspose() const407 virtual const Operator *GetOutputRestrictionTranspose() const 408 { return GetOutputProlongation(); } 409 /// Get the output finite element space restriction matrix GetOutputRestriction() const410 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> FormLinearSystem(const Array<int> & ess_tdof_list,Vector & x,Vector & b,OpType & A,Vector & X,Vector & B,int copy_interior=0)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> FormSystemMatrix(const Array<int> & ess_tdof_list,OpType & A)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. FreeElementMatrices()493 void FreeElementMatrices() 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$ FullInnerProduct(const Vector & x,const Vector & y) const594 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. */ GetFES()602 MFEM_DEPRECATED FiniteElementSpace *GetFES() { return fes; } 603 604 /// Return the FE space associated with the BilinearForm. FESpace()605 FiniteElementSpace *FESpace() { return fes; } 606 /// Read-only access to the associated FiniteElementSpace. FESpace() const607 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 UseExternalIntegrators()619 void UseExternalIntegrators() { extern_bfs = 1; } 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.) 651 AssemblyLevel assembly; 652 653 /** Extension for supporting Full Assembly (FA), Element Assembly (EA), 654 Partial Assembly (PA), or Matrix Free assembly (MF). */ 655 MixedBilinearFormExtension *ext; 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. */ 660 int extern_bfs; 661 662 /// Domain integrators. 663 Array<BilinearFormIntegrator*> domain_integs; 664 665 /// Boundary integrators. 666 Array<BilinearFormIntegrator*> boundary_integs; 667 Array<Array<int>*> boundary_integs_marker; ///< Entries are not owned. 668 669 /// Trace face (skeleton) integrators. 670 Array<BilinearFormIntegrator*> trace_face_integs; 671 672 /// Boundary trace face (skeleton) integrators. 673 Array<BilinearFormIntegrator*> boundary_trace_face_integs; 674 /// Entries are not owned. 675 Array<Array<int>*> boundary_trace_face_integs_marker; 676 677 DenseMatrix elemmat; 678 Array<int> trial_vdofs, test_vdofs; 679 680 private: 681 /// Copy construction is not supported; body is undefined. 682 MixedBilinearForm(const MixedBilinearForm &); 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. */ 692 MixedBilinearForm(FiniteElementSpace *tr_fes, 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. */ 704 MixedBilinearForm(FiniteElementSpace *tr_fes, 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$ SpMat() const735 const SparseMatrix &SpMat() const { return *mat; } 736 737 /// Returns a reference to the sparse matrix: \f$ M \f$ SpMat()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. */ LoseMat()742 SparseMatrix *LoseMat() { SparseMatrix *tmp = mat; mat = NULL; return tmp; } 743 744 /// Adds a domain integrator. Assumes ownership of @a bfi. 745 void AddDomainIntegrator(BilinearFormIntegrator *bfi); 746 747 /// Adds a boundary integrator. Assumes ownership of @a bfi. 748 void AddBoundaryIntegrator(BilinearFormIntegrator *bfi); 749 750 /// Adds a boundary integrator. Assumes ownership of @a bfi. 751 void AddBoundaryIntegrator (BilinearFormIntegrator * 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. */ 759 void AddTraceFaceIntegrator(BilinearFormIntegrator *bfi); 760 761 /// Adds a boundary trace face integrator. Assumes ownership of @a bfi. 762 void AddBdrTraceFaceIntegrator (BilinearFormIntegrator * bfi); 763 764 /// Adds a boundary trace face integrator. Assumes ownership of @a bfi. 765 void AddBdrTraceFaceIntegrator (BilinearFormIntegrator * bfi, 766 Array<int> &bdr_marker); 767 768 /// Access all integrators added with AddDomainIntegrator(). GetDBFI()769 Array<BilinearFormIntegrator*> *GetDBFI() { return &domain_integs; } 770 771 /// Access all integrators added with AddBoundaryIntegrator(). GetBBFI()772 Array<BilinearFormIntegrator*> *GetBBFI() { return &boundary_integs; } 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. */ GetBBFI_Marker()776 Array<Array<int>*> *GetBBFI_Marker() { return &boundary_integs_marker; } 777 778 /// Access all integrators added with AddTraceFaceIntegrator(). GetTFBFI()779 Array<BilinearFormIntegrator*> *GetTFBFI() { return &trace_face_integs; } 780 781 /// Access all integrators added with AddBdrTraceFaceIntegrator(). GetBTFBFI()782 Array<BilinearFormIntegrator*> *GetBTFBFI() 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. */ GetBTFBFI_Marker()787 Array<Array<int>*> *GetBTFBFI_Marker() 788 { return &boundary_trace_face_integs_marker; } 789 790 /// Sets all sparse values of \f$ M \f$ to @a a. operator =(const double 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 GetProlongation() const804 virtual const Operator *GetProlongation() const 805 { return trial_fes->GetProlongationMatrix(); } 806 807 /// Get the input finite element space restriction matrix GetRestriction() const808 virtual const Operator *GetRestriction() const 809 { return trial_fes->GetRestrictionMatrix(); } 810 811 /// Get the test finite element space prolongation matrix GetOutputProlongation() const812 virtual const Operator *GetOutputProlongation() const 813 { return test_fes->GetProlongationMatrix(); } 814 815 /// Get the test finite element space restriction matrix GetOutputRestriction() const816 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, 849 Array<int> &trial_vdofs, Array<int> &test_vdofs, 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, 869 Array<int> &trial_vdofs, Array<int> &test_vdofs, 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> FormRectangularSystemMatrix(const Array<int> & trial_tdof_list,const Array<int> & test_tdof_list,OpType & A)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> FormRectangularLinearSystem(const Array<int> & trial_tdof_list,const Array<int> & test_tdof_list,Vector & x,Vector & b,OpType & A,Vector & X,Vector & B)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. TrialFESpace()943 FiniteElementSpace *TrialFESpace() { return trial_fes; } 944 /// Read-only access to the associated trial FiniteElementSpace. TrialFESpace() const945 const FiniteElementSpace *TrialFESpace() const { return trial_fes; } 946 947 /// Return the test FE space associated with the BilinearForm. TestFESpace()948 FiniteElementSpace *TestFESpace() { return test_fes; } 949 /// Read-only access to the associated test FiniteElementSpace. TestFESpace() const950 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 */ 986 class DiscreteLinearOperator : public MixedBilinearForm 987 { 988 private: 989 /// Copy construction is not supported; body is undefined. 990 DiscreteLinearOperator(const DiscreteLinearOperator &); 991 992 /// Copy assignment is not supported; body is undefined. 993 DiscreteLinearOperator &operator=(const DiscreteLinearOperator &); 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. */ DiscreteLinearOperator(FiniteElementSpace * domain_fes,FiniteElementSpace * range_fes)1000 DiscreteLinearOperator(FiniteElementSpace *domain_fes, 1001 FiniteElementSpace *range_fes) 1002 : MixedBilinearForm(domain_fes, range_fes) { } 1003 1004 /// Adds a domain interpolator. Assumes ownership of @a di. AddDomainInterpolator(DiscreteInterpolator * di)1005 void AddDomainInterpolator(DiscreteInterpolator *di) 1006 { AddDomainIntegrator(di); } 1007 1008 /// Adds a trace face interpolator. Assumes ownership of @a di. AddTraceFaceInterpolator(DiscreteInterpolator * di)1009 void AddTraceFaceInterpolator(DiscreteInterpolator *di) 1010 { AddTraceFaceIntegrator(di); } 1011 1012 /// Access all interpolators added with AddDomainInterpolator(). GetDI()1013 Array<BilinearFormIntegrator*> *GetDI() { return &domain_integs; } 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. */ GetOutputRestrictionTranspose() const1025 virtual const Operator *GetOutputRestrictionTranspose() const 1026 { return test_fes->GetRestrictionTransposeOperator(); } 1027 }; 1028 1029 } 1030 1031 #endif 1032