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_PBILINEARFORM 13 #define MFEM_PBILINEARFORM 14 15 #include "../config/config.hpp" 16 17 #ifdef MFEM_USE_MPI 18 19 #include <mpi.h> 20 #include "pfespace.hpp" 21 #include "pgridfunc.hpp" 22 #include "bilinearform.hpp" 23 24 namespace mfem 25 { 26 27 /// Class for parallel bilinear form 28 class ParBilinearForm : public BilinearForm 29 { 30 friend FABilinearFormExtension; 31 protected: 32 ParFiniteElementSpace *pfes; ///< Points to the same object as #fes 33 34 /// Auxiliary objects used in TrueAddMult(). 35 mutable ParGridFunction X, Y; 36 37 OperatorHandle p_mat, p_mat_e; 38 39 bool keep_nbr_block; 40 41 // Allocate mat - called when (mat == NULL && fbfi.Size() > 0) 42 void pAllocMat(); 43 44 void AssembleSharedFaces(int skip_zeros = 1); 45 46 private: 47 /// Copy construction is not supported; body is undefined. 48 ParBilinearForm(const ParBilinearForm &); 49 50 /// Copy assignment is not supported; body is undefined. 51 ParBilinearForm &operator=(const ParBilinearForm &); 52 53 public: 54 /// Creates parallel bilinear form associated with the FE space @a *pf. 55 /** The pointer @a pf is not owned by the newly constructed object. */ ParBilinearForm(ParFiniteElementSpace * pf)56 ParBilinearForm(ParFiniteElementSpace *pf) 57 : BilinearForm(pf), pfes(pf), 58 p_mat(Operator::Hypre_ParCSR), p_mat_e(Operator::Hypre_ParCSR) 59 { keep_nbr_block = false; } 60 61 /** @brief Create a ParBilinearForm on the ParFiniteElementSpace @a *pf, 62 using the same integrators as the ParBilinearForm @a *bf. 63 64 The pointer @a pf is not owned by the newly constructed object. 65 66 The integrators in @a bf are copied as pointers and they are not owned by 67 the newly constructed ParBilinearForm. */ ParBilinearForm(ParFiniteElementSpace * pf,ParBilinearForm * bf)68 ParBilinearForm(ParFiniteElementSpace *pf, ParBilinearForm *bf) 69 : BilinearForm(pf, bf), pfes(pf), 70 p_mat(Operator::Hypre_ParCSR), p_mat_e(Operator::Hypre_ParCSR) 71 { keep_nbr_block = false; } 72 73 /** When set to true and the ParBilinearForm has interior face integrators, 74 the local SparseMatrix will include the rows (in addition to the columns) 75 corresponding to face-neighbor dofs. The default behavior is to disregard 76 those rows. Must be called before the first Assemble call. */ KeepNbrBlock(bool knb=true)77 void KeepNbrBlock(bool knb = true) { keep_nbr_block = knb; } 78 79 /** @brief Set the operator type id for the parallel matrix/operator when 80 using AssemblyLevel::LEGACY. */ 81 /** If using static condensation or hybridization, call this method *after* 82 enabling it. */ SetOperatorType(Operator::Type tid)83 void SetOperatorType(Operator::Type tid) 84 { 85 p_mat.SetType(tid); p_mat_e.SetType(tid); 86 if (hybridization) { hybridization->SetOperatorType(tid); } 87 if (static_cond) { static_cond->SetOperatorType(tid); } 88 } 89 90 /// Assemble the local matrix 91 void Assemble(int skip_zeros = 1); 92 93 /** @brief Assemble the diagonal of the bilinear form into @a diag. Note that 94 @a diag is a true-dof Vector. 95 96 When the AssemblyLevel is not LEGACY, and the mesh is nonconforming, 97 this method returns |P^T| d_l, where d_l is the local diagonal of the 98 form before applying parallel/conforming assembly, P^T is the transpose 99 of the parallel/conforming prolongation, and |.| denotes the entry-wise 100 absolute value. In general, this is just an approximation of the exact 101 diagonal for this case. */ 102 virtual void AssembleDiagonal(Vector &diag) const; 103 104 /// Returns the matrix assembled on the true dofs, i.e. P^t A P. 105 /** The returned matrix has to be deleted by the caller. */ ParallelAssemble()106 HypreParMatrix *ParallelAssemble() { return ParallelAssemble(mat); } 107 108 /// Returns the eliminated matrix assembled on the true dofs, i.e. P^t A_e P. 109 /** The returned matrix has to be deleted by the caller. */ ParallelAssembleElim()110 HypreParMatrix *ParallelAssembleElim() { return ParallelAssemble(mat_e); } 111 112 /// Return the matrix @a m assembled on the true dofs, i.e. P^t A P. 113 /** The returned matrix has to be deleted by the caller. */ 114 HypreParMatrix *ParallelAssemble(SparseMatrix *m); 115 116 /** @brief Returns the matrix assembled on the true dofs, i.e. 117 @a A = P^t A_local P, in the format (type id) specified by @a A. */ ParallelAssemble(OperatorHandle & A)118 void ParallelAssemble(OperatorHandle &A) { ParallelAssemble(A, mat); } 119 120 /** Returns the eliminated matrix assembled on the true dofs, i.e. 121 @a A_elim = P^t A_elim_local P in the format (type id) specified by @a A. 122 */ ParallelAssembleElim(OperatorHandle & A_elim)123 void ParallelAssembleElim(OperatorHandle &A_elim) 124 { ParallelAssemble(A_elim, mat_e); } 125 126 /** Returns the matrix @a A_local assembled on the true dofs, i.e. 127 @a A = P^t A_local P in the format (type id) specified by @a A. */ 128 void ParallelAssemble(OperatorHandle &A, SparseMatrix *A_local); 129 130 /// Eliminate essential boundary DOFs from a parallel assembled system. 131 /** The array @a bdr_attr_is_ess marks boundary attributes that constitute 132 the essential part of the boundary. */ 133 void ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess, 134 HypreParMatrix &A, 135 const HypreParVector &X, 136 HypreParVector &B) const; 137 138 /// Eliminate essential boundary DOFs from a parallel assembled matrix @a A. 139 /** The array @a bdr_attr_is_ess marks boundary attributes that constitute 140 the essential part of the boundary. The eliminated part is stored in a 141 matrix A_elim such that A_original = A_new + A_elim. Returns a pointer to 142 the newly allocated matrix A_elim which should be deleted by the caller. 143 The matrices @a A and A_elim can be used to eliminate boundary conditions 144 in multiple right-hand sides, by calling the function EliminateBC() (from 145 hypre.hpp). */ 146 HypreParMatrix *ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess, 147 HypreParMatrix &A) const; 148 149 /// Eliminate essential true DOFs from a parallel assembled matrix @a A. 150 /** Given a list of essential true dofs and the parallel assembled matrix 151 @a A, eliminate the true dofs from the matrix, storing the eliminated 152 part in a matrix A_elim such that A_original = A_new + A_elim. Returns a 153 pointer to the newly allocated matrix A_elim which should be deleted by 154 the caller. The matrices @a A and A_elim can be used to eliminate 155 boundary conditions in multiple right-hand sides, by calling the function 156 EliminateBC() (from hypre.hpp). */ ParallelEliminateTDofs(const Array<int> & tdofs_list,HypreParMatrix & A) const157 HypreParMatrix *ParallelEliminateTDofs(const Array<int> &tdofs_list, 158 HypreParMatrix &A) const 159 { return A.EliminateRowsCols(tdofs_list); } 160 161 /** @brief Compute @a y += @a a (P^t A P) @a x, where @a x and @a y are 162 vectors on the true dofs. */ 163 void TrueAddMult(const Vector &x, Vector &y, const double a = 1.0) const; 164 165 /// Return the parallel FE space associated with the ParBilinearForm. ParFESpace() const166 ParFiniteElementSpace *ParFESpace() const { return pfes; } 167 168 /// Return the parallel trace FE space associated with static condensation. SCParFESpace() const169 ParFiniteElementSpace *SCParFESpace() const 170 { return static_cond ? static_cond->GetParTraceFESpace() : NULL; } 171 172 /// Get the parallel finite element space prolongation matrix GetProlongation() const173 virtual const Operator *GetProlongation() const 174 { return pfes->GetProlongationMatrix(); } 175 /// Get the transpose of GetRestriction, useful for matrix-free RAP GetRestrictionTranspose() const176 virtual const Operator *GetRestrictionTranspose() const 177 { return pfes->GetRestrictionTransposeOperator(); } 178 /// Get the parallel finite element space restriction matrix GetRestriction() const179 virtual const Operator *GetRestriction() const 180 { return pfes->GetRestrictionMatrix(); } 181 182 using BilinearForm::FormLinearSystem; 183 using BilinearForm::FormSystemMatrix; 184 185 virtual void FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x, 186 Vector &b, OperatorHandle &A, Vector &X, 187 Vector &B, int copy_interior = 0); 188 189 virtual void FormSystemMatrix(const Array<int> &ess_tdof_list, 190 OperatorHandle &A); 191 192 /** Call this method after solving a linear system constructed using the 193 FormLinearSystem method to recover the solution as a ParGridFunction-size 194 vector in x. Use the same arguments as in the FormLinearSystem call. */ 195 virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x); 196 197 virtual void Update(FiniteElementSpace *nfes = NULL); 198 199 void EliminateVDofsInRHS(const Array<int> &vdofs, const Vector &x, Vector &b); 200 ~ParBilinearForm()201 virtual ~ParBilinearForm() { } 202 }; 203 204 /// Class for parallel bilinear form using different test and trial FE spaces. 205 class ParMixedBilinearForm : public MixedBilinearForm 206 { 207 protected: 208 /// Points to the same object as #trial_fes 209 ParFiniteElementSpace *trial_pfes; 210 /// Points to the same object as #test_fes 211 ParFiniteElementSpace *test_pfes; 212 /// Auxiliary objects used in TrueAddMult(). 213 mutable ParGridFunction X, Y; 214 215 /// Matrix and eliminated matrix 216 OperatorHandle p_mat, p_mat_e; 217 218 private: 219 /// Copy construction is not supported; body is undefined. 220 ParMixedBilinearForm(const ParMixedBilinearForm &); 221 222 /// Copy assignment is not supported; body is undefined. 223 ParMixedBilinearForm &operator=(const ParMixedBilinearForm &); 224 225 public: 226 /** @brief Construct a ParMixedBilinearForm on the given FiniteElementSpace%s 227 @a trial_fes and @a test_fes. */ 228 /** The pointers @a trial_fes and @a test_fes are not owned by the newly 229 constructed object. */ ParMixedBilinearForm(ParFiniteElementSpace * trial_fes,ParFiniteElementSpace * test_fes)230 ParMixedBilinearForm(ParFiniteElementSpace *trial_fes, 231 ParFiniteElementSpace *test_fes) 232 : MixedBilinearForm(trial_fes, test_fes), 233 p_mat(Operator::Hypre_ParCSR), p_mat_e(Operator::Hypre_ParCSR) 234 { 235 trial_pfes = trial_fes; 236 test_pfes = test_fes; 237 } 238 239 /** @brief Create a ParMixedBilinearForm on the given FiniteElementSpace%s 240 @a trial_fes and @a test_fes, using the same integrators as the 241 ParMixedBilinearForm @a mbf. 242 243 The pointers @a trial_fes and @a test_fes are not owned by the newly 244 constructed object. 245 246 The integrators in @a mbf are copied as pointers and they are not owned 247 by the newly constructed ParMixedBilinearForm. */ ParMixedBilinearForm(ParFiniteElementSpace * trial_fes,ParFiniteElementSpace * test_fes,ParMixedBilinearForm * mbf)248 ParMixedBilinearForm(ParFiniteElementSpace *trial_fes, 249 ParFiniteElementSpace *test_fes, 250 ParMixedBilinearForm * mbf) 251 : MixedBilinearForm(trial_fes, test_fes, mbf), 252 p_mat(Operator::Hypre_ParCSR), p_mat_e(Operator::Hypre_ParCSR) 253 { 254 trial_pfes = trial_fes; 255 test_pfes = test_fes; 256 } 257 258 /// Returns the matrix assembled on the true dofs, i.e. P_test^t A P_trial. 259 HypreParMatrix *ParallelAssemble(); 260 261 /** @brief Returns the matrix assembled on the true dofs, i.e. 262 @a A = P_test^t A_local P_trial, in the format (type id) specified by 263 @a A. */ 264 void ParallelAssemble(OperatorHandle &A); 265 266 using MixedBilinearForm::FormRectangularSystemMatrix; 267 using MixedBilinearForm::FormRectangularLinearSystem; 268 269 /** @brief Return in @a A a parallel (on truedofs) version of this operator. 270 271 This returns the same operator as FormRectangularLinearSystem(), but does 272 without the transformations of the right-hand side. */ 273 virtual void FormRectangularSystemMatrix(const Array<int> &trial_tdof_list, 274 const Array<int> &test_tdof_list, 275 OperatorHandle &A); 276 277 /** @brief Form the parallel linear system A X = B, corresponding to this mixed 278 bilinear form and the linear form @a b(.). 279 280 Return in @a A a *reference* to the system matrix that is column-constrained. 281 The reference will be invalidated when SetOperatorType(), Update(), or the 282 destructor is called. */ 283 virtual void FormRectangularLinearSystem(const Array<int> &trial_tdof_list, 284 const Array<int> &test_tdof_list, Vector &x, 285 Vector &b, OperatorHandle &A, Vector &X, 286 Vector &B); 287 288 /// Compute y += a (P^t A P) x, where x and y are vectors on the true dofs 289 void TrueAddMult(const Vector &x, Vector &y, const double a = 1.0) const; 290 ~ParMixedBilinearForm()291 virtual ~ParMixedBilinearForm() { } 292 }; 293 294 /** The parallel matrix representation a linear operator between parallel finite 295 element spaces */ 296 class ParDiscreteLinearOperator : public DiscreteLinearOperator 297 { 298 protected: 299 /// Points to the same object as #trial_fes 300 ParFiniteElementSpace *domain_fes; 301 /// Points to the same object as #test_fes 302 ParFiniteElementSpace *range_fes; 303 304 private: 305 /// Copy construction is not supported; body is undefined. 306 ParDiscreteLinearOperator(const ParDiscreteLinearOperator &); 307 308 /// Copy assignment is not supported; body is undefined. 309 ParDiscreteLinearOperator &operator=(const ParDiscreteLinearOperator &); 310 311 public: 312 /** @brief Construct a ParDiscreteLinearOperator on the given 313 FiniteElementSpace%s @a dfes (domain FE space) and @a rfes (range FE 314 space). */ 315 /** The pointers @a dfes and @a rfes are not owned by the newly constructed 316 object. */ ParDiscreteLinearOperator(ParFiniteElementSpace * dfes,ParFiniteElementSpace * rfes)317 ParDiscreteLinearOperator(ParFiniteElementSpace *dfes, 318 ParFiniteElementSpace *rfes) 319 : DiscreteLinearOperator(dfes, rfes) { domain_fes=dfes; range_fes=rfes; } 320 321 /// Returns the matrix "assembled" on the true dofs 322 HypreParMatrix *ParallelAssemble() const; 323 324 /** @brief Returns the matrix assembled on the true dofs, i.e. 325 @a A = R_test A_local P_trial, in the format (type id) specified by 326 @a A. */ 327 void ParallelAssemble(OperatorHandle &A); 328 329 /** Extract the parallel blocks corresponding to the vector dimensions of the 330 domain and range parallel finite element spaces */ 331 void GetParBlocks(Array2D<HypreParMatrix *> &blocks) const; 332 333 using MixedBilinearForm::FormRectangularSystemMatrix; 334 335 /** @brief Return in @a A a parallel (on truedofs) version of this operator. */ 336 virtual void FormRectangularSystemMatrix(OperatorHandle &A); 337 ~ParDiscreteLinearOperator()338 virtual ~ParDiscreteLinearOperator() { } 339 }; 340 341 } 342 343 #endif // MFEM_USE_MPI 344 345 #endif 346