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_BILININTEG 13 #define MFEM_BILININTEG 14 15 #include "../config/config.hpp" 16 #include "nonlininteg.hpp" 17 #include "fespace.hpp" 18 19 namespace mfem 20 { 21 22 // Local maximum size of dofs and quads in 1D 23 constexpr int HCURL_MAX_D1D = 5; 24 #ifdef MFEM_USE_HIP 25 constexpr int HCURL_MAX_Q1D = 5; 26 #else 27 constexpr int HCURL_MAX_Q1D = 6; 28 #endif 29 30 constexpr int HDIV_MAX_D1D = 5; 31 constexpr int HDIV_MAX_Q1D = 6; 32 33 /// Abstract base class BilinearFormIntegrator 34 class BilinearFormIntegrator : public NonlinearFormIntegrator 35 { 36 protected: BilinearFormIntegrator(const IntegrationRule * ir=NULL)37 BilinearFormIntegrator(const IntegrationRule *ir = NULL) 38 : NonlinearFormIntegrator(ir) { } 39 40 public: 41 // TODO: add support for other assembly levels (in addition to PA) and their 42 // actions. 43 44 // TODO: for mixed meshes the quadrature rules to be used by methods like 45 // AssemblePA() can be given as a QuadratureSpace, e.g. using a new method: 46 // SetQuadratureSpace(). 47 48 // TODO: the methods for the various assembly levels make sense even in the 49 // base class NonlinearFormIntegrator, except that not all assembly levels 50 // make sense for the action of the nonlinear operator (but they all make 51 // sense for its Jacobian). 52 53 using NonlinearFormIntegrator::AssemblePA; 54 55 /// Method defining partial assembly. 56 /** The result of the partial assembly is stored internally so that it can be 57 used later in the methods AddMultPA() and AddMultTransposePA(). */ 58 virtual void AssemblePA(const FiniteElementSpace &fes); 59 /** Used with BilinearFormIntegrators that have different spaces. */ 60 virtual void AssemblePA(const FiniteElementSpace &trial_fes, 61 const FiniteElementSpace &test_fes); 62 63 virtual void AssemblePAInteriorFaces(const FiniteElementSpace &fes); 64 65 virtual void AssemblePABoundaryFaces(const FiniteElementSpace &fes); 66 67 /// Assemble diagonal and add it to Vector @a diag. 68 virtual void AssembleDiagonalPA(Vector &diag); 69 70 /// Assemble diagonal of ADA^T (A is this integrator) and add it to @a diag. 71 virtual void AssembleDiagonalPA_ADAt(const Vector &D, Vector &diag); 72 73 /// Method for partially assembled action. 74 /** Perform the action of integrator on the input @a x and add the result to 75 the output @a y. Both @a x and @a y are E-vectors, i.e. they represent 76 the element-wise discontinuous version of the FE space. 77 78 This method can be called only after the method AssemblePA() has been 79 called. */ 80 virtual void AddMultPA(const Vector &x, Vector &y) const; 81 82 /// Method for partially assembled transposed action. 83 /** Perform the transpose action of integrator on the input @a x and add the 84 result to the output @a y. Both @a x and @a y are E-vectors, i.e. they 85 represent the element-wise discontinuous version of the FE space. 86 87 This method can be called only after the method AssemblePA() has been 88 called. */ 89 virtual void AddMultTransposePA(const Vector &x, Vector &y) const; 90 91 /// Method defining element assembly. 92 /** The result of the element assembly is added to the @a emat Vector if 93 @a add is true. Otherwise, if @a add is false, we set @a emat. */ 94 virtual void AssembleEA(const FiniteElementSpace &fes, Vector &emat, 95 const bool add = true); 96 /** Used with BilinearFormIntegrators that have different spaces. */ 97 // virtual void AssembleEA(const FiniteElementSpace &trial_fes, 98 // const FiniteElementSpace &test_fes, 99 // Vector &emat); 100 101 /// Method defining matrix-free assembly. 102 /** The result of fully matrix-free assembly is stored internally so that it 103 can be used later in the methods AddMultMF() and AddMultTransposeMF(). */ 104 virtual void AssembleMF(const FiniteElementSpace &fes); 105 106 /** Perform the action of integrator on the input @a x and add the result to 107 the output @a y. Both @a x and @a y are E-vectors, i.e. they represent 108 the element-wise discontinuous version of the FE space. 109 110 This method can be called only after the method AssembleMF() has been 111 called. */ 112 virtual void AddMultMF(const Vector &x, Vector &y) const; 113 114 /** Perform the transpose action of integrator on the input @a x and add the 115 result to the output @a y. Both @a x and @a y are E-vectors, i.e. they 116 represent the element-wise discontinuous version of the FE space. 117 118 This method can be called only after the method AssemblePA() has been 119 called. */ 120 virtual void AddMultTransposeMF(const Vector &x, Vector &y) const; 121 122 /// Assemble diagonal and add it to Vector @a diag. 123 virtual void AssembleDiagonalMF(Vector &diag); 124 125 virtual void AssembleEAInteriorFaces(const FiniteElementSpace &fes, 126 Vector &ea_data_int, 127 Vector &ea_data_ext, 128 const bool add = true); 129 130 virtual void AssembleEABoundaryFaces(const FiniteElementSpace &fes, 131 Vector &ea_data_bdr, 132 const bool add = true); 133 134 /// Given a particular Finite Element computes the element matrix elmat. 135 virtual void AssembleElementMatrix(const FiniteElement &el, 136 ElementTransformation &Trans, 137 DenseMatrix &elmat); 138 139 /** Compute the local matrix representation of a bilinear form 140 a(u,v) defined on different trial (given by u) and test 141 (given by v) spaces. The rows in the local matrix correspond 142 to the test dofs and the columns -- to the trial dofs. */ 143 virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, 144 const FiniteElement &test_fe, 145 ElementTransformation &Trans, 146 DenseMatrix &elmat); 147 148 virtual void AssembleFaceMatrix(const FiniteElement &el1, 149 const FiniteElement &el2, 150 FaceElementTransformations &Trans, 151 DenseMatrix &elmat); 152 153 /** Abstract method used for assembling TraceFaceIntegrators in a 154 MixedBilinearForm. */ 155 virtual void AssembleFaceMatrix(const FiniteElement &trial_face_fe, 156 const FiniteElement &test_fe1, 157 const FiniteElement &test_fe2, 158 FaceElementTransformations &Trans, 159 DenseMatrix &elmat); 160 161 /// @brief Perform the local action of the BilinearFormIntegrator. 162 /// Note that the default implementation in the base class is general but not 163 /// efficient. 164 virtual void AssembleElementVector(const FiniteElement &el, 165 ElementTransformation &Tr, 166 const Vector &elfun, Vector &elvect); 167 168 /// @brief Perform the local action of the BilinearFormIntegrator resulting 169 /// from a face integral term. 170 /// Note that the default implementation in the base class is general but not 171 /// efficient. 172 virtual void AssembleFaceVector(const FiniteElement &el1, 173 const FiniteElement &el2, 174 FaceElementTransformations &Tr, 175 const Vector &elfun, Vector &elvect); 176 AssembleElementGrad(const FiniteElement & el,ElementTransformation & Tr,const Vector & elfun,DenseMatrix & elmat)177 virtual void AssembleElementGrad(const FiniteElement &el, 178 ElementTransformation &Tr, 179 const Vector &elfun, DenseMatrix &elmat) 180 { AssembleElementMatrix(el, Tr, elmat); } 181 AssembleFaceGrad(const FiniteElement & el1,const FiniteElement & el2,FaceElementTransformations & Tr,const Vector & elfun,DenseMatrix & elmat)182 virtual void AssembleFaceGrad(const FiniteElement &el1, 183 const FiniteElement &el2, 184 FaceElementTransformations &Tr, 185 const Vector &elfun, DenseMatrix &elmat) 186 { AssembleFaceMatrix(el1, el2, Tr, elmat); } 187 188 /** @brief Virtual method required for Zienkiewicz-Zhu type error estimators. 189 190 The purpose of the method is to compute a local "flux" finite element 191 function given a local finite element solution. The "flux" function has 192 to be computed in terms of its coefficients (represented by the Vector 193 @a flux) which multiply the basis functions defined by the FiniteElement 194 @a fluxelem. Typically, the "flux" function will have more than one 195 component and consequently @a flux should be store the coefficients of 196 all components: first all coefficient for component 0, then all 197 coefficients for component 1, etc. What the "flux" function represents 198 depends on the specific integrator. For example, in the case of 199 DiffusionIntegrator, the flux is the gradient of the solution multiplied 200 by the diffusion coefficient. 201 202 @param[in] el FiniteElement of the solution. 203 @param[in] Trans The ElementTransformation describing the physical 204 position of the mesh element. 205 @param[in] u Solution coefficients representing the expansion of the 206 solution function in the basis of @a el. 207 @param[in] fluxelem FiniteElement of the "flux". 208 @param[out] flux "Flux" coefficients representing the expansion of the 209 "flux" function in the basis of @a fluxelem. The size 210 of @a flux as a Vector has to be set by this method, 211 e.g. using Vector::SetSize(). 212 @param[in] with_coef If zero (the default value is 1) the implementation 213 of the method may choose not to scale the "flux" 214 function by any coefficients describing the 215 integrator. 216 */ ComputeElementFlux(const FiniteElement & el,ElementTransformation & Trans,Vector & u,const FiniteElement & fluxelem,Vector & flux,bool with_coef=true)217 virtual void ComputeElementFlux(const FiniteElement &el, 218 ElementTransformation &Trans, 219 Vector &u, 220 const FiniteElement &fluxelem, 221 Vector &flux, bool with_coef = true) { } 222 223 /** @brief Virtual method required for Zienkiewicz-Zhu type error estimators. 224 225 The purpose of this method is to compute a local number that measures the 226 energy of a given "flux" function (see ComputeElementFlux() for a 227 description of the "flux" function). Typically, the energy of a "flux" 228 function should be equal to a_local(u,u), if the "flux" is defined from 229 a solution u; here a_local(.,.) denotes the element-local bilinear 230 form represented by the integrator. 231 232 @param[in] fluxelem FiniteElement of the "flux". 233 @param[in] Trans The ElementTransformation describing the physical 234 position of the mesh element. 235 @param[in] flux "Flux" coefficients representing the expansion of the 236 "flux" function in the basis of @a fluxelem. 237 @param[out] d_energy If not NULL, the given Vector should be set to 238 represent directional energy split that can be used 239 for anisotropic error estimation. 240 @returns The computed energy. 241 */ ComputeFluxEnergy(const FiniteElement & fluxelem,ElementTransformation & Trans,Vector & flux,Vector * d_energy=NULL)242 virtual double ComputeFluxEnergy(const FiniteElement &fluxelem, 243 ElementTransformation &Trans, 244 Vector &flux, Vector *d_energy = NULL) 245 { return 0.0; } 246 ~BilinearFormIntegrator()247 virtual ~BilinearFormIntegrator() { } 248 }; 249 250 /** Wraps a given @a BilinearFormIntegrator and transposes the resulting element 251 matrices. See for example ex9, ex9p. */ 252 class TransposeIntegrator : public BilinearFormIntegrator 253 { 254 private: 255 int own_bfi; 256 BilinearFormIntegrator *bfi; 257 258 DenseMatrix bfi_elmat; 259 260 public: TransposeIntegrator(BilinearFormIntegrator * bfi_,int own_bfi_=1)261 TransposeIntegrator (BilinearFormIntegrator *bfi_, int own_bfi_ = 1) 262 { bfi = bfi_; own_bfi = own_bfi_; } 263 264 virtual void SetIntRule(const IntegrationRule *ir); 265 266 virtual void AssembleElementMatrix(const FiniteElement &el, 267 ElementTransformation &Trans, 268 DenseMatrix &elmat); 269 270 virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, 271 const FiniteElement &test_fe, 272 ElementTransformation &Trans, 273 DenseMatrix &elmat); 274 275 using BilinearFormIntegrator::AssembleFaceMatrix; 276 virtual void AssembleFaceMatrix(const FiniteElement &el1, 277 const FiniteElement &el2, 278 FaceElementTransformations &Trans, 279 DenseMatrix &elmat); 280 281 using BilinearFormIntegrator::AssemblePA; 282 AssemblePA(const FiniteElementSpace & fes)283 virtual void AssemblePA(const FiniteElementSpace& fes) 284 { 285 bfi->AssemblePA(fes); 286 } 287 AssemblePAInteriorFaces(const FiniteElementSpace & fes)288 virtual void AssemblePAInteriorFaces(const FiniteElementSpace &fes) 289 { 290 bfi->AssemblePAInteriorFaces(fes); 291 } 292 AssemblePABoundaryFaces(const FiniteElementSpace & fes)293 virtual void AssemblePABoundaryFaces(const FiniteElementSpace &fes) 294 { 295 bfi->AssemblePABoundaryFaces(fes); 296 } 297 AddMultTransposePA(const Vector & x,Vector & y) const298 virtual void AddMultTransposePA(const Vector &x, Vector &y) const 299 { 300 bfi->AddMultPA(x, y); 301 } 302 AddMultPA(const Vector & x,Vector & y) const303 virtual void AddMultPA(const Vector& x, Vector& y) const 304 { 305 bfi->AddMultTransposePA(x, y); 306 } 307 308 virtual void AssembleEA(const FiniteElementSpace &fes, Vector &emat, 309 const bool add); 310 311 virtual void AssembleEAInteriorFaces(const FiniteElementSpace &fes, 312 Vector &ea_data_int, 313 Vector &ea_data_ext, 314 const bool add); 315 316 virtual void AssembleEABoundaryFaces(const FiniteElementSpace &fes, 317 Vector &ea_data_bdr, 318 const bool add); 319 ~TransposeIntegrator()320 virtual ~TransposeIntegrator() { if (own_bfi) { delete bfi; } } 321 }; 322 323 class LumpedIntegrator : public BilinearFormIntegrator 324 { 325 private: 326 int own_bfi; 327 BilinearFormIntegrator *bfi; 328 329 public: LumpedIntegrator(BilinearFormIntegrator * bfi_,int own_bfi_=1)330 LumpedIntegrator (BilinearFormIntegrator *bfi_, int own_bfi_ = 1) 331 { bfi = bfi_; own_bfi = own_bfi_; } 332 333 virtual void SetIntRule(const IntegrationRule *ir); 334 335 virtual void AssembleElementMatrix(const FiniteElement &el, 336 ElementTransformation &Trans, 337 DenseMatrix &elmat); 338 ~LumpedIntegrator()339 virtual ~LumpedIntegrator() { if (own_bfi) { delete bfi; } } 340 }; 341 342 /// Integrator that inverts the matrix assembled by another integrator. 343 class InverseIntegrator : public BilinearFormIntegrator 344 { 345 private: 346 int own_integrator; 347 BilinearFormIntegrator *integrator; 348 349 public: InverseIntegrator(BilinearFormIntegrator * integ,int own_integ=1)350 InverseIntegrator(BilinearFormIntegrator *integ, int own_integ = 1) 351 { integrator = integ; own_integrator = own_integ; } 352 353 virtual void SetIntRule(const IntegrationRule *ir); 354 355 virtual void AssembleElementMatrix(const FiniteElement &el, 356 ElementTransformation &Trans, 357 DenseMatrix &elmat); 358 ~InverseIntegrator()359 virtual ~InverseIntegrator() { if (own_integrator) { delete integrator; } } 360 }; 361 362 /// Integrator defining a sum of multiple Integrators. 363 class SumIntegrator : public BilinearFormIntegrator 364 { 365 private: 366 int own_integrators; 367 mutable DenseMatrix elem_mat; 368 Array<BilinearFormIntegrator*> integrators; 369 370 public: SumIntegrator(int own_integs=1)371 SumIntegrator(int own_integs = 1) { own_integrators = own_integs; } 372 373 virtual void SetIntRule(const IntegrationRule *ir); 374 AddIntegrator(BilinearFormIntegrator * integ)375 void AddIntegrator(BilinearFormIntegrator *integ) 376 { integrators.Append(integ); } 377 378 virtual void AssembleElementMatrix(const FiniteElement &el, 379 ElementTransformation &Trans, 380 DenseMatrix &elmat); 381 virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, 382 const FiniteElement &test_fe, 383 ElementTransformation &Trans, 384 DenseMatrix &elmat); 385 386 using BilinearFormIntegrator::AssembleFaceMatrix; 387 virtual void AssembleFaceMatrix(const FiniteElement &el1, 388 const FiniteElement &el2, 389 FaceElementTransformations &Trans, 390 DenseMatrix &elmat); 391 392 virtual void AssembleFaceMatrix(const FiniteElement &trial_face_fe, 393 const FiniteElement &test_fe1, 394 const FiniteElement &test_fe2, 395 FaceElementTransformations &Trans, 396 DenseMatrix &elmat); 397 398 using BilinearFormIntegrator::AssemblePA; 399 virtual void AssemblePA(const FiniteElementSpace& fes); 400 401 virtual void AssembleDiagonalPA(Vector &diag); 402 403 virtual void AssemblePAInteriorFaces(const FiniteElementSpace &fes); 404 405 virtual void AssemblePABoundaryFaces(const FiniteElementSpace &fes); 406 407 virtual void AddMultTransposePA(const Vector &x, Vector &y) const; 408 409 virtual void AddMultPA(const Vector& x, Vector& y) const; 410 411 virtual void AssembleMF(const FiniteElementSpace &fes); 412 413 virtual void AddMultMF(const Vector &x, Vector &y) const; 414 415 virtual void AddMultTransposeMF(const Vector &x, Vector &y) const; 416 417 virtual void AssembleDiagonalMF(Vector &diag); 418 419 virtual void AssembleEA(const FiniteElementSpace &fes, Vector &emat, 420 const bool add); 421 422 virtual void AssembleEAInteriorFaces(const FiniteElementSpace &fes, 423 Vector &ea_data_int, 424 Vector &ea_data_ext, 425 const bool add); 426 427 virtual void AssembleEABoundaryFaces(const FiniteElementSpace &fes, 428 Vector &ea_data_bdr, 429 const bool add); 430 431 virtual ~SumIntegrator(); 432 }; 433 434 /** An abstract class for integrating the product of two scalar basis functions 435 with an optional scalar coefficient. */ 436 class MixedScalarIntegrator: public BilinearFormIntegrator 437 { 438 public: 439 440 virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, 441 const FiniteElement &test_fe, 442 ElementTransformation &Trans, 443 DenseMatrix &elmat); 444 445 /// Support for use in BilinearForm. Can be used only when appropriate. AssembleElementMatrix(const FiniteElement & fe,ElementTransformation & Trans,DenseMatrix & elmat)446 virtual void AssembleElementMatrix(const FiniteElement &fe, 447 ElementTransformation &Trans, 448 DenseMatrix &elmat) 449 { AssembleElementMatrix2(fe, fe, Trans, elmat); } 450 451 protected: 452 /// This parameter can be set by derived methods to enable single shape 453 /// evaluation in case CalcTestShape() and CalcTrialShape() return the same 454 /// result if given the same FiniteElement. The default is false. 455 bool same_calc_shape; 456 MixedScalarIntegrator()457 MixedScalarIntegrator() : same_calc_shape(false), Q(NULL) {} MixedScalarIntegrator(Coefficient & q)458 MixedScalarIntegrator(Coefficient &q) : same_calc_shape(false), Q(&q) {} 459 VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const460 inline virtual bool VerifyFiniteElementTypes( 461 const FiniteElement & trial_fe, 462 const FiniteElement & test_fe) const 463 { 464 return (trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR && 465 test_fe.GetRangeType() == mfem::FiniteElement::SCALAR ); 466 } 467 FiniteElementTypeFailureMessage() const468 inline virtual const char * FiniteElementTypeFailureMessage() const 469 { 470 return "MixedScalarIntegrator: " 471 "Trial and test spaces must both be scalar fields."; 472 } 473 GetIntegrationOrder(const FiniteElement & trial_fe,const FiniteElement & test_fe,ElementTransformation & Trans)474 inline virtual int GetIntegrationOrder(const FiniteElement & trial_fe, 475 const FiniteElement & test_fe, 476 ElementTransformation &Trans) 477 { return trial_fe.GetOrder() + test_fe.GetOrder() + Trans.OrderW(); } 478 479 CalcTestShape(const FiniteElement & test_fe,ElementTransformation & Trans,Vector & shape)480 inline virtual void CalcTestShape(const FiniteElement & test_fe, 481 ElementTransformation &Trans, 482 Vector & shape) 483 { test_fe.CalcPhysShape(Trans, shape); } 484 CalcTrialShape(const FiniteElement & trial_fe,ElementTransformation & Trans,Vector & shape)485 inline virtual void CalcTrialShape(const FiniteElement & trial_fe, 486 ElementTransformation &Trans, 487 Vector & shape) 488 { trial_fe.CalcPhysShape(Trans, shape); } 489 490 Coefficient *Q; 491 492 private: 493 494 #ifndef MFEM_THREAD_SAFE 495 Vector test_shape; 496 Vector trial_shape; 497 #endif 498 499 }; 500 501 /** An abstract class for integrating the inner product of two vector basis 502 functions with an optional scalar, vector, or matrix coefficient. */ 503 class MixedVectorIntegrator: public BilinearFormIntegrator 504 { 505 public: 506 507 virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, 508 const FiniteElement &test_fe, 509 ElementTransformation &Trans, 510 DenseMatrix &elmat); 511 512 /// Support for use in BilinearForm. Can be used only when appropriate. AssembleElementMatrix(const FiniteElement & fe,ElementTransformation & Trans,DenseMatrix & elmat)513 virtual void AssembleElementMatrix(const FiniteElement &fe, 514 ElementTransformation &Trans, 515 DenseMatrix &elmat) 516 { AssembleElementMatrix2(fe, fe, Trans, elmat); } 517 518 protected: 519 /// This parameter can be set by derived methods to enable single shape 520 /// evaluation in case CalcTestShape() and CalcTrialShape() return the same 521 /// result if given the same FiniteElement. The default is false. 522 bool same_calc_shape; 523 MixedVectorIntegrator()524 MixedVectorIntegrator() 525 : same_calc_shape(false), Q(NULL), VQ(NULL), DQ(NULL), MQ(NULL) {} MixedVectorIntegrator(Coefficient & q)526 MixedVectorIntegrator(Coefficient &q) 527 : same_calc_shape(false), Q(&q), VQ(NULL), DQ(NULL), MQ(NULL) {} MixedVectorIntegrator(VectorCoefficient & vq,bool diag=true)528 MixedVectorIntegrator(VectorCoefficient &vq, bool diag = true) 529 : same_calc_shape(false), Q(NULL), VQ(diag?NULL:&vq), DQ(diag?&vq:NULL), 530 MQ(NULL) {} MixedVectorIntegrator(MatrixCoefficient & mq)531 MixedVectorIntegrator(MatrixCoefficient &mq) 532 : same_calc_shape(false), Q(NULL), VQ(NULL), DQ(NULL), MQ(&mq) {} 533 VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const534 inline virtual bool VerifyFiniteElementTypes( 535 const FiniteElement & trial_fe, 536 const FiniteElement & test_fe) const 537 { 538 return (trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR && 539 test_fe.GetRangeType() == mfem::FiniteElement::VECTOR ); 540 } 541 FiniteElementTypeFailureMessage() const542 inline virtual const char * FiniteElementTypeFailureMessage() const 543 { 544 return "MixedVectorIntegrator: " 545 "Trial and test spaces must both be vector fields"; 546 } 547 GetIntegrationOrder(const FiniteElement & trial_fe,const FiniteElement & test_fe,ElementTransformation & Trans)548 inline virtual int GetIntegrationOrder(const FiniteElement & trial_fe, 549 const FiniteElement & test_fe, 550 ElementTransformation &Trans) 551 { return trial_fe.GetOrder() + test_fe.GetOrder() + Trans.OrderW(); } 552 553 CalcTestShape(const FiniteElement & test_fe,ElementTransformation & Trans,DenseMatrix & shape)554 inline virtual void CalcTestShape(const FiniteElement & test_fe, 555 ElementTransformation &Trans, 556 DenseMatrix & shape) 557 { test_fe.CalcVShape(Trans, shape); } 558 CalcTrialShape(const FiniteElement & trial_fe,ElementTransformation & Trans,DenseMatrix & shape)559 inline virtual void CalcTrialShape(const FiniteElement & trial_fe, 560 ElementTransformation &Trans, 561 DenseMatrix & shape) 562 { trial_fe.CalcVShape(Trans, shape); } 563 564 Coefficient *Q; 565 VectorCoefficient *VQ; 566 DiagonalMatrixCoefficient *DQ; 567 MatrixCoefficient *MQ; 568 569 private: 570 571 #ifndef MFEM_THREAD_SAFE 572 Vector V; 573 Vector D; 574 DenseMatrix M; 575 DenseMatrix test_shape; 576 DenseMatrix trial_shape; 577 DenseMatrix test_shape_tmp; 578 #endif 579 580 }; 581 582 /** An abstract class for integrating the product of a scalar basis function and 583 the inner product of a vector basis function with a vector coefficient. In 584 2D the inner product can be replaced with a cross product. */ 585 class MixedScalarVectorIntegrator: public BilinearFormIntegrator 586 { 587 public: 588 589 virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, 590 const FiniteElement &test_fe, 591 ElementTransformation &Trans, 592 DenseMatrix &elmat); 593 594 /// Support for use in BilinearForm. Can be used only when appropriate. 595 /** Appropriate use cases are classes derived from 596 MixedScalarVectorIntegrator where the trial and test spaces can be the 597 same. Examples of such classes are: MixedVectorDivergenceIntegrator, 598 MixedScalarWeakDivergenceIntegrator, etc. */ AssembleElementMatrix(const FiniteElement & fe,ElementTransformation & Trans,DenseMatrix & elmat)599 virtual void AssembleElementMatrix(const FiniteElement &fe, 600 ElementTransformation &Trans, 601 DenseMatrix &elmat) 602 { AssembleElementMatrix2(fe, fe, Trans, elmat); } 603 604 protected: 605 MixedScalarVectorIntegrator(VectorCoefficient & vq,bool transpose_=false,bool cross_2d_=false)606 MixedScalarVectorIntegrator(VectorCoefficient &vq, bool transpose_ = false, 607 bool cross_2d_ = false) 608 : VQ(&vq), transpose(transpose_), cross_2d(cross_2d_) {} 609 VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const610 inline virtual bool VerifyFiniteElementTypes( 611 const FiniteElement & trial_fe, 612 const FiniteElement & test_fe) const 613 { 614 return ((transpose && 615 trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR && 616 test_fe.GetRangeType() == mfem::FiniteElement::SCALAR ) || 617 (!transpose && 618 trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR && 619 test_fe.GetRangeType() == mfem::FiniteElement::VECTOR ) 620 ); 621 } 622 FiniteElementTypeFailureMessage() const623 inline virtual const char * FiniteElementTypeFailureMessage() const 624 { 625 if ( transpose ) 626 { 627 return "MixedScalarVectorIntegrator: " 628 "Trial space must be a vector field " 629 "and the test space must be a scalar field"; 630 } 631 else 632 { 633 return "MixedScalarVectorIntegrator: " 634 "Trial space must be a scalar field " 635 "and the test space must be a vector field"; 636 } 637 } 638 GetIntegrationOrder(const FiniteElement & trial_fe,const FiniteElement & test_fe,ElementTransformation & Trans)639 inline virtual int GetIntegrationOrder(const FiniteElement & trial_fe, 640 const FiniteElement & test_fe, 641 ElementTransformation &Trans) 642 { return trial_fe.GetOrder() + test_fe.GetOrder() + Trans.OrderW(); } 643 644 CalcVShape(const FiniteElement & vector_fe,ElementTransformation & Trans,DenseMatrix & shape)645 inline virtual void CalcVShape(const FiniteElement & vector_fe, 646 ElementTransformation &Trans, 647 DenseMatrix & shape) 648 { vector_fe.CalcVShape(Trans, shape); } 649 CalcShape(const FiniteElement & scalar_fe,ElementTransformation & Trans,Vector & shape)650 inline virtual void CalcShape(const FiniteElement & scalar_fe, 651 ElementTransformation &Trans, 652 Vector & shape) 653 { scalar_fe.CalcPhysShape(Trans, shape); } 654 655 VectorCoefficient *VQ; 656 bool transpose; 657 bool cross_2d; // In 2D use a cross product rather than a dot product 658 659 private: 660 661 #ifndef MFEM_THREAD_SAFE 662 Vector V; 663 DenseMatrix vshape; 664 Vector shape; 665 Vector vshape_tmp; 666 #endif 667 668 }; 669 670 /** Class for integrating the bilinear form a(u,v) := (Q u, v) in either 1D, 2D, 671 or 3D and where Q is an optional scalar coefficient, u and v are each in H1 672 or L2. */ 673 class MixedScalarMassIntegrator : public MixedScalarIntegrator 674 { 675 public: MixedScalarMassIntegrator()676 MixedScalarMassIntegrator() { same_calc_shape = true; } MixedScalarMassIntegrator(Coefficient & q)677 MixedScalarMassIntegrator(Coefficient &q) 678 : MixedScalarIntegrator(q) { same_calc_shape = true; } 679 }; 680 681 /** Class for integrating the bilinear form a(u,v) := (Q u, v) in either 2D, or 682 3D and where Q is a vector coefficient, u is in H1 or L2 and v is in H(Curl) 683 or H(Div). */ 684 class MixedVectorProductIntegrator : public MixedScalarVectorIntegrator 685 { 686 public: MixedVectorProductIntegrator(VectorCoefficient & vq)687 MixedVectorProductIntegrator(VectorCoefficient &vq) 688 : MixedScalarVectorIntegrator(vq) {} 689 }; 690 691 /** Class for integrating the bilinear form a(u,v) := (Q D u, v) in 1D where Q 692 is an optional scalar coefficient, u is in H1, and v is in L2. */ 693 class MixedScalarDerivativeIntegrator : public MixedScalarIntegrator 694 { 695 public: MixedScalarDerivativeIntegrator()696 MixedScalarDerivativeIntegrator() {} MixedScalarDerivativeIntegrator(Coefficient & q)697 MixedScalarDerivativeIntegrator(Coefficient &q) 698 : MixedScalarIntegrator(q) {} 699 700 protected: VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const701 inline virtual bool VerifyFiniteElementTypes( 702 const FiniteElement & trial_fe, 703 const FiniteElement & test_fe) const 704 { 705 return (trial_fe.GetDim() == 1 && test_fe.GetDim() == 1 && 706 trial_fe.GetDerivType() == mfem::FiniteElement::GRAD && 707 test_fe.GetRangeType() == mfem::FiniteElement::SCALAR ); 708 } 709 FiniteElementTypeFailureMessage() const710 inline virtual const char * FiniteElementTypeFailureMessage() const 711 { 712 return "MixedScalarDerivativeIntegrator: " 713 "Trial and test spaces must both be scalar fields in 1D " 714 "and the trial space must implement CalcDShape."; 715 } 716 CalcTrialShape(const FiniteElement & trial_fe,ElementTransformation & Trans,Vector & shape)717 inline virtual void CalcTrialShape(const FiniteElement & trial_fe, 718 ElementTransformation &Trans, 719 Vector & shape) 720 { 721 DenseMatrix dshape(shape.GetData(), shape.Size(), 1); 722 trial_fe.CalcPhysDShape(Trans, dshape); 723 } 724 }; 725 726 /** Class for integrating the bilinear form a(u,v) := -(Q u, D v) in 1D where Q 727 is an optional scalar coefficient, u is in L2, and v is in H1. */ 728 class MixedScalarWeakDerivativeIntegrator : public MixedScalarIntegrator 729 { 730 public: MixedScalarWeakDerivativeIntegrator()731 MixedScalarWeakDerivativeIntegrator() {} MixedScalarWeakDerivativeIntegrator(Coefficient & q)732 MixedScalarWeakDerivativeIntegrator(Coefficient &q) 733 : MixedScalarIntegrator(q) {} 734 735 protected: VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const736 inline virtual bool VerifyFiniteElementTypes( 737 const FiniteElement & trial_fe, 738 const FiniteElement & test_fe) const 739 { 740 return (trial_fe.GetDim() == 1 && test_fe.GetDim() == 1 && 741 trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR && 742 test_fe.GetDerivType() == mfem::FiniteElement::GRAD ); 743 } 744 FiniteElementTypeFailureMessage() const745 inline virtual const char * FiniteElementTypeFailureMessage() const 746 { 747 return "MixedScalarWeakDerivativeIntegrator: " 748 "Trial and test spaces must both be scalar fields in 1D " 749 "and the test space must implement CalcDShape with " 750 "map type \"VALUE\"."; 751 } 752 CalcTestShape(const FiniteElement & test_fe,ElementTransformation & Trans,Vector & shape)753 inline virtual void CalcTestShape(const FiniteElement & test_fe, 754 ElementTransformation &Trans, 755 Vector & shape) 756 { 757 DenseMatrix dshape(shape.GetData(), shape.Size(), 1); 758 test_fe.CalcPhysDShape(Trans, dshape); 759 shape *= -1.0; 760 } 761 }; 762 763 /** Class for integrating the bilinear form a(u,v) := (Q div u, v) in either 2D 764 or 3D where Q is an optional scalar coefficient, u is in H(Div), and v is a 765 scalar field. */ 766 class MixedScalarDivergenceIntegrator : public MixedScalarIntegrator 767 { 768 public: MixedScalarDivergenceIntegrator()769 MixedScalarDivergenceIntegrator() {} MixedScalarDivergenceIntegrator(Coefficient & q)770 MixedScalarDivergenceIntegrator(Coefficient &q) 771 : MixedScalarIntegrator(q) {} 772 773 protected: VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const774 inline virtual bool VerifyFiniteElementTypes( 775 const FiniteElement & trial_fe, 776 const FiniteElement & test_fe) const 777 { 778 return (trial_fe.GetDerivType() == mfem::FiniteElement::DIV && 779 test_fe.GetRangeType() == mfem::FiniteElement::SCALAR ); 780 } 781 FiniteElementTypeFailureMessage() const782 inline virtual const char * FiniteElementTypeFailureMessage() const 783 { 784 return "MixedScalarDivergenceIntegrator: " 785 "Trial must be H(Div) and the test space must be a " 786 "scalar field"; 787 } 788 GetIntegrationOrder(const FiniteElement & trial_fe,const FiniteElement & test_fe,ElementTransformation & Trans)789 inline virtual int GetIntegrationOrder(const FiniteElement & trial_fe, 790 const FiniteElement & test_fe, 791 ElementTransformation &Trans) 792 { return trial_fe.GetOrder() + test_fe.GetOrder() + Trans.OrderW() - 1; } 793 CalcTrialShape(const FiniteElement & trial_fe,ElementTransformation & Trans,Vector & shape)794 inline virtual void CalcTrialShape(const FiniteElement & trial_fe, 795 ElementTransformation &Trans, 796 Vector & shape) 797 { trial_fe.CalcPhysDivShape(Trans, shape); } 798 }; 799 800 /** Class for integrating the bilinear form a(u,v) := (V div u, v) in either 2D 801 or 3D where V is a vector coefficient, u is in H(Div), and v is a vector 802 field. */ 803 class MixedVectorDivergenceIntegrator : public MixedScalarVectorIntegrator 804 { 805 public: MixedVectorDivergenceIntegrator(VectorCoefficient & vq)806 MixedVectorDivergenceIntegrator(VectorCoefficient &vq) 807 : MixedScalarVectorIntegrator(vq) {} 808 809 protected: VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const810 inline virtual bool VerifyFiniteElementTypes( 811 const FiniteElement & trial_fe, 812 const FiniteElement & test_fe) const 813 { 814 return (trial_fe.GetDerivType() == mfem::FiniteElement::DIV && 815 test_fe.GetRangeType() == mfem::FiniteElement::VECTOR ); 816 } 817 FiniteElementTypeFailureMessage() const818 inline virtual const char * FiniteElementTypeFailureMessage() const 819 { 820 return "MixedVectorDivergenceIntegrator: " 821 "Trial must be H(Div) and the test space must be a " 822 "vector field"; 823 } 824 825 // Subtract one due to the divergence and add one for the coefficient 826 // which is assumed to be at least linear. GetIntegrationOrder(const FiniteElement & trial_fe,const FiniteElement & test_fe,ElementTransformation & Trans)827 inline virtual int GetIntegrationOrder(const FiniteElement & trial_fe, 828 const FiniteElement & test_fe, 829 ElementTransformation &Trans) 830 { return trial_fe.GetOrder() + test_fe.GetOrder() + Trans.OrderW() - 1 + 1; } 831 CalcShape(const FiniteElement & scalar_fe,ElementTransformation & Trans,Vector & shape)832 inline virtual void CalcShape(const FiniteElement & scalar_fe, 833 ElementTransformation &Trans, 834 Vector & shape) 835 { scalar_fe.CalcPhysDivShape(Trans, shape); } 836 }; 837 838 /** Class for integrating the bilinear form a(u,v) := -(Q u, div v) in either 2D 839 or 3D where Q is an optional scalar coefficient, u is in L2 or H1, and v is 840 in H(Div). */ 841 class MixedScalarWeakGradientIntegrator : public MixedScalarIntegrator 842 { 843 public: MixedScalarWeakGradientIntegrator()844 MixedScalarWeakGradientIntegrator() {} MixedScalarWeakGradientIntegrator(Coefficient & q)845 MixedScalarWeakGradientIntegrator(Coefficient &q) 846 : MixedScalarIntegrator(q) {} 847 848 protected: VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const849 inline virtual bool VerifyFiniteElementTypes( 850 const FiniteElement & trial_fe, 851 const FiniteElement & test_fe) const 852 { 853 return (trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR && 854 test_fe.GetDerivType() == mfem::FiniteElement::DIV ); 855 } 856 FiniteElementTypeFailureMessage() const857 inline virtual const char * FiniteElementTypeFailureMessage() const 858 { 859 return "MixedScalarWeakGradientIntegrator: " 860 "Trial space must be a scalar field " 861 "and the test space must be H(Div)"; 862 } 863 GetIntegrationOrder(const FiniteElement & trial_fe,const FiniteElement & test_fe,ElementTransformation & Trans)864 inline virtual int GetIntegrationOrder(const FiniteElement & trial_fe, 865 const FiniteElement & test_fe, 866 ElementTransformation &Trans) 867 { return trial_fe.GetOrder() + test_fe.GetOrder() + Trans.OrderW() - 1; } 868 CalcTestShape(const FiniteElement & test_fe,ElementTransformation & Trans,Vector & shape)869 virtual void CalcTestShape(const FiniteElement & test_fe, 870 ElementTransformation &Trans, 871 Vector & shape) 872 { 873 test_fe.CalcPhysDivShape(Trans, shape); 874 shape *= -1.0; 875 } 876 }; 877 878 /** Class for integrating the bilinear form a(u,v) := (Q curl u, v) in 2D where 879 Q is an optional scalar coefficient, u is in H(Curl), and v is in L2 or 880 H1. */ 881 class MixedScalarCurlIntegrator : public MixedScalarIntegrator 882 { 883 public: MixedScalarCurlIntegrator()884 MixedScalarCurlIntegrator() {} MixedScalarCurlIntegrator(Coefficient & q)885 MixedScalarCurlIntegrator(Coefficient &q) 886 : MixedScalarIntegrator(q) {} 887 888 protected: VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const889 inline virtual bool VerifyFiniteElementTypes( 890 const FiniteElement & trial_fe, 891 const FiniteElement & test_fe) const 892 { 893 return (trial_fe.GetDim() == 2 && test_fe.GetDim() == 2 && 894 trial_fe.GetDerivType() == mfem::FiniteElement::CURL && 895 test_fe.GetRangeType() == mfem::FiniteElement::SCALAR ); 896 } 897 FiniteElementTypeFailureMessage() const898 inline virtual const char * FiniteElementTypeFailureMessage() const 899 { 900 return "MixedScalarCurlIntegrator: " 901 "Trial must be H(Curl) and the test space must be a " 902 "scalar field"; 903 } 904 GetIntegrationOrder(const FiniteElement & trial_fe,const FiniteElement & test_fe,ElementTransformation & Trans)905 inline virtual int GetIntegrationOrder(const FiniteElement & trial_fe, 906 const FiniteElement & test_fe, 907 ElementTransformation &Trans) 908 { return trial_fe.GetOrder() + test_fe.GetOrder() + Trans.OrderW() - 1; } 909 CalcTrialShape(const FiniteElement & trial_fe,ElementTransformation & Trans,Vector & shape)910 inline virtual void CalcTrialShape(const FiniteElement & trial_fe, 911 ElementTransformation &Trans, 912 Vector & shape) 913 { 914 DenseMatrix dshape(shape.GetData(), shape.Size(), 1); 915 trial_fe.CalcPhysCurlShape(Trans, dshape); 916 } 917 }; 918 919 /** Class for integrating the bilinear form a(u,v) := (Q u, curl v) in 2D where 920 Q is an optional scalar coefficient, u is in L2 or H1, and v is in 921 H(Curl). */ 922 class MixedScalarWeakCurlIntegrator : public MixedScalarIntegrator 923 { 924 public: MixedScalarWeakCurlIntegrator()925 MixedScalarWeakCurlIntegrator() {} MixedScalarWeakCurlIntegrator(Coefficient & q)926 MixedScalarWeakCurlIntegrator(Coefficient &q) 927 : MixedScalarIntegrator(q) {} 928 929 protected: VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const930 inline virtual bool VerifyFiniteElementTypes( 931 const FiniteElement & trial_fe, 932 const FiniteElement & test_fe) const 933 { 934 return (trial_fe.GetDim() == 2 && test_fe.GetDim() == 2 && 935 trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR && 936 test_fe.GetDerivType() == mfem::FiniteElement::CURL ); 937 } 938 FiniteElementTypeFailureMessage() const939 inline virtual const char * FiniteElementTypeFailureMessage() const 940 { 941 return "MixedScalarWeakCurlIntegrator: " 942 "Trial space must be a scalar field " 943 "and the test space must be H(Curl)"; 944 } 945 CalcTestShape(const FiniteElement & test_fe,ElementTransformation & Trans,Vector & shape)946 inline virtual void CalcTestShape(const FiniteElement & test_fe, 947 ElementTransformation &Trans, 948 Vector & shape) 949 { 950 DenseMatrix dshape(shape.GetData(), shape.Size(), 1); 951 test_fe.CalcPhysCurlShape(Trans, dshape); 952 } 953 }; 954 955 /** Class for integrating the bilinear form a(u,v) := (Q u, v) in either 2D or 956 3D and where Q is an optional coefficient (of type scalar, matrix, or 957 diagonal matrix) u and v are each in H(Curl) or H(Div). */ 958 class MixedVectorMassIntegrator : public MixedVectorIntegrator 959 { 960 public: MixedVectorMassIntegrator()961 MixedVectorMassIntegrator() { same_calc_shape = true; } MixedVectorMassIntegrator(Coefficient & q)962 MixedVectorMassIntegrator(Coefficient &q) 963 : MixedVectorIntegrator(q) { same_calc_shape = true; } MixedVectorMassIntegrator(DiagonalMatrixCoefficient & dq)964 MixedVectorMassIntegrator(DiagonalMatrixCoefficient &dq) 965 : MixedVectorIntegrator(dq, true) { same_calc_shape = true; } MixedVectorMassIntegrator(MatrixCoefficient & mq)966 MixedVectorMassIntegrator(MatrixCoefficient &mq) 967 : MixedVectorIntegrator(mq) { same_calc_shape = true; } 968 }; 969 970 /** Class for integrating the bilinear form a(u,v) := (V x u, v) in 3D and where 971 V is a vector coefficient u and v are each in H(Curl) or H(Div). */ 972 class MixedCrossProductIntegrator : public MixedVectorIntegrator 973 { 974 public: MixedCrossProductIntegrator(VectorCoefficient & vq)975 MixedCrossProductIntegrator(VectorCoefficient &vq) 976 : MixedVectorIntegrator(vq, false) { same_calc_shape = true; } 977 }; 978 979 /** Class for integrating the bilinear form a(u,v) := (V . u, v) in 2D or 3D and 980 where V is a vector coefficient u is in H(Curl) or H(Div) and v is in H1 or 981 L2. */ 982 class MixedDotProductIntegrator : public MixedScalarVectorIntegrator 983 { 984 public: MixedDotProductIntegrator(VectorCoefficient & vq)985 MixedDotProductIntegrator(VectorCoefficient &vq) 986 : MixedScalarVectorIntegrator(vq, true) {} 987 VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const988 inline virtual bool VerifyFiniteElementTypes( 989 const FiniteElement & trial_fe, 990 const FiniteElement & test_fe) const 991 { 992 return (trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR && 993 test_fe.GetRangeType() == mfem::FiniteElement::SCALAR ); 994 } 995 FiniteElementTypeFailureMessage() const996 inline virtual const char * FiniteElementTypeFailureMessage() const 997 { 998 return "MixedDotProductIntegrator: " 999 "Trial space must be a vector field " 1000 "and the test space must be a scalar field"; 1001 } 1002 }; 1003 1004 /** Class for integrating the bilinear form a(u,v) := (-V . u, Div v) in 2D or 1005 3D and where V is a vector coefficient u is in H(Curl) or H(Div) and v is in 1006 RT. */ 1007 class MixedWeakGradDotIntegrator : public MixedScalarVectorIntegrator 1008 { 1009 public: MixedWeakGradDotIntegrator(VectorCoefficient & vq)1010 MixedWeakGradDotIntegrator(VectorCoefficient &vq) 1011 : MixedScalarVectorIntegrator(vq, true) {} 1012 VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const1013 inline virtual bool VerifyFiniteElementTypes( 1014 const FiniteElement & trial_fe, 1015 const FiniteElement & test_fe) const 1016 { 1017 return (trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR && 1018 test_fe.GetRangeType() == mfem::FiniteElement::VECTOR && 1019 test_fe.GetDerivType() == mfem::FiniteElement::DIV ); 1020 } 1021 FiniteElementTypeFailureMessage() const1022 inline virtual const char * FiniteElementTypeFailureMessage() const 1023 { 1024 return "MixedWeakGradDotIntegrator: " 1025 "Trial space must be a vector field " 1026 "and the test space must be a vector field with a divergence"; 1027 } 1028 1029 // Subtract one due to the gradient and add one for the coefficient 1030 // which is assumed to be at least linear. GetIntegrationOrder(const FiniteElement & trial_fe,const FiniteElement & test_fe,ElementTransformation & Trans)1031 inline virtual int GetIntegrationOrder(const FiniteElement & trial_fe, 1032 const FiniteElement & test_fe, 1033 ElementTransformation &Trans) 1034 { return trial_fe.GetOrder() + test_fe.GetOrder() + Trans.OrderW() - 1 + 1; } 1035 CalcShape(const FiniteElement & scalar_fe,ElementTransformation & Trans,Vector & shape)1036 inline virtual void CalcShape(const FiniteElement & scalar_fe, 1037 ElementTransformation &Trans, 1038 Vector & shape) 1039 { scalar_fe.CalcPhysDivShape(Trans, shape); shape *= -1.0; } 1040 }; 1041 1042 /** Class for integrating the bilinear form a(u,v) := (V x u, Grad v) in 3D and 1043 where V is a vector coefficient u is in H(Curl) or H(Div) and v is in H1. */ 1044 class MixedWeakDivCrossIntegrator : public MixedVectorIntegrator 1045 { 1046 public: MixedWeakDivCrossIntegrator(VectorCoefficient & vq)1047 MixedWeakDivCrossIntegrator(VectorCoefficient &vq) 1048 : MixedVectorIntegrator(vq, false) {} 1049 VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const1050 inline virtual bool VerifyFiniteElementTypes( 1051 const FiniteElement & trial_fe, 1052 const FiniteElement & test_fe) const 1053 { 1054 return (trial_fe.GetDim() == 3 && test_fe.GetDim() == 3 && 1055 trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR && 1056 test_fe.GetRangeType() == mfem::FiniteElement::SCALAR && 1057 test_fe.GetDerivType() == mfem::FiniteElement::GRAD ); 1058 } 1059 FiniteElementTypeFailureMessage() const1060 inline virtual const char * FiniteElementTypeFailureMessage() const 1061 { 1062 return "MixedWeakDivCrossIntegrator: " 1063 "Trial space must be a vector field in 3D " 1064 "and the test space must be a scalar field with a gradient"; 1065 } 1066 CalcTestShape(const FiniteElement & test_fe,ElementTransformation & Trans,DenseMatrix & shape)1067 inline virtual void CalcTestShape(const FiniteElement & test_fe, 1068 ElementTransformation &Trans, 1069 DenseMatrix & shape) 1070 { test_fe.CalcPhysDShape(Trans, shape); shape *= -1.0; } 1071 }; 1072 1073 /** Class for integrating the bilinear form a(u,v) := (Q Grad u, Grad v) in 3D 1074 or in 2D and where Q is a scalar or matrix coefficient u and v are both in 1075 H1. */ 1076 class MixedGradGradIntegrator : public MixedVectorIntegrator 1077 { 1078 public: MixedGradGradIntegrator()1079 MixedGradGradIntegrator() { same_calc_shape = true; } MixedGradGradIntegrator(Coefficient & q)1080 MixedGradGradIntegrator(Coefficient &q) 1081 : MixedVectorIntegrator(q) { same_calc_shape = true; } MixedGradGradIntegrator(DiagonalMatrixCoefficient & dq)1082 MixedGradGradIntegrator(DiagonalMatrixCoefficient &dq) 1083 : MixedVectorIntegrator(dq, true) { same_calc_shape = true; } MixedGradGradIntegrator(MatrixCoefficient & mq)1084 MixedGradGradIntegrator(MatrixCoefficient &mq) 1085 : MixedVectorIntegrator(mq) { same_calc_shape = true; } 1086 VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const1087 inline virtual bool VerifyFiniteElementTypes( 1088 const FiniteElement & trial_fe, 1089 const FiniteElement & test_fe) const 1090 { 1091 return (trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR && 1092 trial_fe.GetDerivType() == mfem::FiniteElement::GRAD && 1093 test_fe.GetRangeType() == mfem::FiniteElement::SCALAR && 1094 test_fe.GetDerivType() == mfem::FiniteElement::GRAD ); 1095 } 1096 FiniteElementTypeFailureMessage() const1097 inline virtual const char * FiniteElementTypeFailureMessage() const 1098 { 1099 return "MixedGradGradIntegrator: " 1100 "Trial and test spaces must both be scalar fields " 1101 "with a gradient operator."; 1102 } 1103 GetIntegrationOrder(const FiniteElement & trial_fe,const FiniteElement & test_fe,ElementTransformation & Trans)1104 inline virtual int GetIntegrationOrder(const FiniteElement & trial_fe, 1105 const FiniteElement & test_fe, 1106 ElementTransformation &Trans) 1107 { 1108 // Same as DiffusionIntegrator 1109 return test_fe.Space() == FunctionSpace::Pk ? 1110 trial_fe.GetOrder() + test_fe.GetOrder() - 2 : 1111 trial_fe.GetOrder() + test_fe.GetOrder() + test_fe.GetDim() - 1; 1112 } 1113 CalcTrialShape(const FiniteElement & trial_fe,ElementTransformation & Trans,DenseMatrix & shape)1114 inline virtual void CalcTrialShape(const FiniteElement & trial_fe, 1115 ElementTransformation &Trans, 1116 DenseMatrix & shape) 1117 { trial_fe.CalcPhysDShape(Trans, shape); } 1118 CalcTestShape(const FiniteElement & test_fe,ElementTransformation & Trans,DenseMatrix & shape)1119 inline virtual void CalcTestShape(const FiniteElement & test_fe, 1120 ElementTransformation &Trans, 1121 DenseMatrix & shape) 1122 { test_fe.CalcPhysDShape(Trans, shape); } 1123 }; 1124 1125 /** Class for integrating the bilinear form a(u,v) := (V x Grad u, Grad v) in 3D 1126 or in 2D and where V is a vector coefficient u and v are both in H1. */ 1127 class MixedCrossGradGradIntegrator : public MixedVectorIntegrator 1128 { 1129 public: MixedCrossGradGradIntegrator(VectorCoefficient & vq)1130 MixedCrossGradGradIntegrator(VectorCoefficient &vq) 1131 : MixedVectorIntegrator(vq, false) { same_calc_shape = true; } 1132 VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const1133 inline virtual bool VerifyFiniteElementTypes( 1134 const FiniteElement & trial_fe, 1135 const FiniteElement & test_fe) const 1136 { 1137 return (trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR && 1138 trial_fe.GetDerivType() == mfem::FiniteElement::GRAD && 1139 test_fe.GetRangeType() == mfem::FiniteElement::SCALAR && 1140 test_fe.GetDerivType() == mfem::FiniteElement::GRAD ); 1141 } 1142 FiniteElementTypeFailureMessage() const1143 inline virtual const char * FiniteElementTypeFailureMessage() const 1144 { 1145 return "MixedCrossGradGradIntegrator: " 1146 "Trial and test spaces must both be scalar fields " 1147 "with a gradient operator."; 1148 } 1149 CalcTrialShape(const FiniteElement & trial_fe,ElementTransformation & Trans,DenseMatrix & shape)1150 inline virtual void CalcTrialShape(const FiniteElement & trial_fe, 1151 ElementTransformation &Trans, 1152 DenseMatrix & shape) 1153 { trial_fe.CalcPhysDShape(Trans, shape); } 1154 CalcTestShape(const FiniteElement & test_fe,ElementTransformation & Trans,DenseMatrix & shape)1155 inline virtual void CalcTestShape(const FiniteElement & test_fe, 1156 ElementTransformation &Trans, 1157 DenseMatrix & shape) 1158 { test_fe.CalcPhysDShape(Trans, shape); } 1159 }; 1160 1161 /** Class for integrating the bilinear form a(u,v) := (Q Curl u, Curl v) in 3D 1162 and where Q is a scalar or matrix coefficient u and v are both in 1163 H(Curl). */ 1164 class MixedCurlCurlIntegrator : public MixedVectorIntegrator 1165 { 1166 public: MixedCurlCurlIntegrator()1167 MixedCurlCurlIntegrator() { same_calc_shape = true; } MixedCurlCurlIntegrator(Coefficient & q)1168 MixedCurlCurlIntegrator(Coefficient &q) 1169 : MixedVectorIntegrator(q) { same_calc_shape = true; } MixedCurlCurlIntegrator(DiagonalMatrixCoefficient & dq)1170 MixedCurlCurlIntegrator(DiagonalMatrixCoefficient &dq) 1171 : MixedVectorIntegrator(dq, true) { same_calc_shape = true; } MixedCurlCurlIntegrator(MatrixCoefficient & mq)1172 MixedCurlCurlIntegrator(MatrixCoefficient &mq) 1173 : MixedVectorIntegrator(mq) { same_calc_shape = true; } 1174 VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const1175 inline virtual bool VerifyFiniteElementTypes( 1176 const FiniteElement & trial_fe, 1177 const FiniteElement & test_fe) const 1178 { 1179 return (trial_fe.GetDim() == 3 && test_fe.GetDim() == 3 && 1180 trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR && 1181 trial_fe.GetDerivType() == mfem::FiniteElement::CURL && 1182 test_fe.GetRangeType() == mfem::FiniteElement::VECTOR && 1183 test_fe.GetDerivType() == mfem::FiniteElement::CURL ); 1184 } 1185 FiniteElementTypeFailureMessage() const1186 inline virtual const char * FiniteElementTypeFailureMessage() const 1187 { 1188 return "MixedCurlCurlIntegrator" 1189 "Trial and test spaces must both be vector fields in 3D " 1190 "with a curl."; 1191 } 1192 CalcTrialShape(const FiniteElement & trial_fe,ElementTransformation & Trans,DenseMatrix & shape)1193 inline virtual void CalcTrialShape(const FiniteElement & trial_fe, 1194 ElementTransformation &Trans, 1195 DenseMatrix & shape) 1196 { trial_fe.CalcPhysCurlShape(Trans, shape); } 1197 CalcTestShape(const FiniteElement & test_fe,ElementTransformation & Trans,DenseMatrix & shape)1198 inline virtual void CalcTestShape(const FiniteElement & test_fe, 1199 ElementTransformation &Trans, 1200 DenseMatrix & shape) 1201 { test_fe.CalcPhysCurlShape(Trans, shape); } 1202 }; 1203 1204 /** Class for integrating the bilinear form a(u,v) := (V x Curl u, Curl v) in 3D 1205 and where V is a vector coefficient u and v are both in H(Curl). */ 1206 class MixedCrossCurlCurlIntegrator : public MixedVectorIntegrator 1207 { 1208 public: MixedCrossCurlCurlIntegrator(VectorCoefficient & vq)1209 MixedCrossCurlCurlIntegrator(VectorCoefficient &vq) 1210 : MixedVectorIntegrator(vq, false) { same_calc_shape = true; } 1211 VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const1212 inline virtual bool VerifyFiniteElementTypes( 1213 const FiniteElement & trial_fe, 1214 const FiniteElement & test_fe) const 1215 { 1216 return (trial_fe.GetDim() == 3 && test_fe.GetDim() == 3 && 1217 trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR && 1218 trial_fe.GetDerivType() == mfem::FiniteElement::CURL && 1219 test_fe.GetRangeType() == mfem::FiniteElement::VECTOR && 1220 test_fe.GetDerivType() == mfem::FiniteElement::CURL ); 1221 } 1222 FiniteElementTypeFailureMessage() const1223 inline virtual const char * FiniteElementTypeFailureMessage() const 1224 { 1225 return "MixedCrossCurlCurlIntegrator: " 1226 "Trial and test spaces must both be vector fields in 3D " 1227 "with a curl."; 1228 } 1229 CalcTrialShape(const FiniteElement & trial_fe,ElementTransformation & Trans,DenseMatrix & shape)1230 inline virtual void CalcTrialShape(const FiniteElement & trial_fe, 1231 ElementTransformation &Trans, 1232 DenseMatrix & shape) 1233 { trial_fe.CalcPhysCurlShape(Trans, shape); } 1234 CalcTestShape(const FiniteElement & test_fe,ElementTransformation & Trans,DenseMatrix & shape)1235 inline virtual void CalcTestShape(const FiniteElement & test_fe, 1236 ElementTransformation &Trans, 1237 DenseMatrix & shape) 1238 { test_fe.CalcPhysCurlShape(Trans, shape); } 1239 }; 1240 1241 /** Class for integrating the bilinear form a(u,v) := (V x Curl u, Grad v) in 3D 1242 and where V is a vector coefficient u is in H(Curl) and v is in H1. */ 1243 class MixedCrossCurlGradIntegrator : public MixedVectorIntegrator 1244 { 1245 public: MixedCrossCurlGradIntegrator(VectorCoefficient & vq)1246 MixedCrossCurlGradIntegrator(VectorCoefficient &vq) 1247 : MixedVectorIntegrator(vq, false) {} 1248 VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const1249 inline virtual bool VerifyFiniteElementTypes( 1250 const FiniteElement & trial_fe, 1251 const FiniteElement & test_fe) const 1252 { 1253 return (trial_fe.GetDim() == 3 && test_fe.GetDim() == 3 && 1254 trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR && 1255 trial_fe.GetDerivType() == mfem::FiniteElement::CURL && 1256 test_fe.GetRangeType() == mfem::FiniteElement::SCALAR && 1257 test_fe.GetDerivType() == mfem::FiniteElement::GRAD ); 1258 } 1259 FiniteElementTypeFailureMessage() const1260 inline virtual const char * FiniteElementTypeFailureMessage() const 1261 { 1262 return "MixedCrossCurlGradIntegrator" 1263 "Trial space must be a vector field in 3D with a curl" 1264 "and the test space must be a scalar field with a gradient"; 1265 } 1266 CalcTrialShape(const FiniteElement & trial_fe,ElementTransformation & Trans,DenseMatrix & shape)1267 inline virtual void CalcTrialShape(const FiniteElement & trial_fe, 1268 ElementTransformation &Trans, 1269 DenseMatrix & shape) 1270 { trial_fe.CalcPhysCurlShape(Trans, shape); } 1271 CalcTestShape(const FiniteElement & test_fe,ElementTransformation & Trans,DenseMatrix & shape)1272 inline virtual void CalcTestShape(const FiniteElement & test_fe, 1273 ElementTransformation &Trans, 1274 DenseMatrix & shape) 1275 { test_fe.CalcPhysDShape(Trans, shape); } 1276 }; 1277 1278 /** Class for integrating the bilinear form a(u,v) := (V x Grad u, Curl v) in 3D 1279 and where V is a scalar coefficient u is in H1 and v is in H(Curl). */ 1280 class MixedCrossGradCurlIntegrator : public MixedVectorIntegrator 1281 { 1282 public: MixedCrossGradCurlIntegrator(VectorCoefficient & vq)1283 MixedCrossGradCurlIntegrator(VectorCoefficient &vq) 1284 : MixedVectorIntegrator(vq, false) {} 1285 VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const1286 inline virtual bool VerifyFiniteElementTypes( 1287 const FiniteElement & trial_fe, 1288 const FiniteElement & test_fe) const 1289 { 1290 return (trial_fe.GetDim() == 3 && test_fe.GetDim() == 3 && 1291 trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR && 1292 trial_fe.GetDerivType() == mfem::FiniteElement::GRAD && 1293 test_fe.GetRangeType() == mfem::FiniteElement::VECTOR && 1294 test_fe.GetDerivType() == mfem::FiniteElement::CURL ); 1295 } 1296 FiniteElementTypeFailureMessage() const1297 inline virtual const char * FiniteElementTypeFailureMessage() const 1298 { 1299 return "MixedCrossGradCurlIntegrator" 1300 "Trial space must be a scalar field in 3D with a gradient" 1301 "and the test space must be a vector field with a curl"; 1302 } 1303 CalcTrialShape(const FiniteElement & trial_fe,ElementTransformation & Trans,DenseMatrix & shape)1304 inline virtual void CalcTrialShape(const FiniteElement & trial_fe, 1305 ElementTransformation &Trans, 1306 DenseMatrix & shape) 1307 { trial_fe.CalcPhysDShape(Trans, shape); } 1308 CalcTestShape(const FiniteElement & test_fe,ElementTransformation & Trans,DenseMatrix & shape)1309 inline virtual void CalcTestShape(const FiniteElement & test_fe, 1310 ElementTransformation &Trans, 1311 DenseMatrix & shape) 1312 { test_fe.CalcPhysCurlShape(Trans, shape); } 1313 }; 1314 1315 /** Class for integrating the bilinear form a(u,v) := (V x u, Curl v) in 3D and 1316 where V is a vector coefficient u is in H(Curl) or H(Div) and v is in 1317 H(Curl). */ 1318 class MixedWeakCurlCrossIntegrator : public MixedVectorIntegrator 1319 { 1320 public: MixedWeakCurlCrossIntegrator(VectorCoefficient & vq)1321 MixedWeakCurlCrossIntegrator(VectorCoefficient &vq) 1322 : MixedVectorIntegrator(vq, false) {} 1323 VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const1324 inline virtual bool VerifyFiniteElementTypes( 1325 const FiniteElement & trial_fe, 1326 const FiniteElement & test_fe) const 1327 { 1328 return (trial_fe.GetDim() == 3 && test_fe.GetDim() == 3 && 1329 trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR && 1330 test_fe.GetRangeType() == mfem::FiniteElement::VECTOR && 1331 test_fe.GetDerivType() == mfem::FiniteElement::CURL ); 1332 } 1333 FiniteElementTypeFailureMessage() const1334 inline virtual const char * FiniteElementTypeFailureMessage() const 1335 { 1336 return "MixedWeakCurlCrossIntegrator: " 1337 "Trial space must be a vector field in 3D " 1338 "and the test space must be a vector field with a curl"; 1339 } 1340 CalcTestShape(const FiniteElement & test_fe,ElementTransformation & Trans,DenseMatrix & shape)1341 inline virtual void CalcTestShape(const FiniteElement & test_fe, 1342 ElementTransformation &Trans, 1343 DenseMatrix & shape) 1344 { test_fe.CalcPhysCurlShape(Trans, shape); } 1345 }; 1346 1347 /** Class for integrating the bilinear form a(u,v) := (V x u, Curl v) in 2D and 1348 where V is a vector coefficient u is in H(Curl) or H(Div) and v is in 1349 H(Curl). */ 1350 class MixedScalarWeakCurlCrossIntegrator : public MixedScalarVectorIntegrator 1351 { 1352 public: MixedScalarWeakCurlCrossIntegrator(VectorCoefficient & vq)1353 MixedScalarWeakCurlCrossIntegrator(VectorCoefficient &vq) 1354 : MixedScalarVectorIntegrator(vq, true, true) {} 1355 VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const1356 inline virtual bool VerifyFiniteElementTypes( 1357 const FiniteElement & trial_fe, 1358 const FiniteElement & test_fe) const 1359 { 1360 return (trial_fe.GetDim() == 2 && test_fe.GetDim() == 2 && 1361 trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR && 1362 test_fe.GetRangeType() == mfem::FiniteElement::VECTOR && 1363 test_fe.GetDerivType() == mfem::FiniteElement::CURL ); 1364 } 1365 FiniteElementTypeFailureMessage() const1366 inline virtual const char * FiniteElementTypeFailureMessage() const 1367 { 1368 return "MixedScalarWeakCurlCrossIntegrator: " 1369 "Trial space must be a vector field in 2D " 1370 "and the test space must be a vector field with a curl"; 1371 } 1372 CalcShape(const FiniteElement & scalar_fe,ElementTransformation & Trans,Vector & shape)1373 inline virtual void CalcShape(const FiniteElement & scalar_fe, 1374 ElementTransformation &Trans, 1375 Vector & shape) 1376 { 1377 DenseMatrix dshape(shape.GetData(), shape.Size(), 1); 1378 scalar_fe.CalcPhysCurlShape(Trans, dshape); 1379 } 1380 }; 1381 1382 /** Class for integrating the bilinear form a(u,v) := (V x Grad u, v) in 3D or 1383 in 2D and where V is a vector coefficient u is in H1 and v is in H(Curl) or 1384 H(Div). */ 1385 class MixedCrossGradIntegrator : public MixedVectorIntegrator 1386 { 1387 public: MixedCrossGradIntegrator(VectorCoefficient & vq)1388 MixedCrossGradIntegrator(VectorCoefficient &vq) 1389 : MixedVectorIntegrator(vq, false) {} 1390 VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const1391 inline virtual bool VerifyFiniteElementTypes( 1392 const FiniteElement & trial_fe, 1393 const FiniteElement & test_fe) const 1394 { 1395 return (trial_fe.GetDim() == 3 && test_fe.GetDim() == 3 && 1396 trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR && 1397 trial_fe.GetDerivType() == mfem::FiniteElement::GRAD && 1398 test_fe.GetRangeType() == mfem::FiniteElement::VECTOR ); 1399 } 1400 FiniteElementTypeFailureMessage() const1401 inline virtual const char * FiniteElementTypeFailureMessage() const 1402 { 1403 return "MixedCrossGradIntegrator: " 1404 "Trial space must be a scalar field with a gradient operator" 1405 " and the test space must be a vector field both in 3D."; 1406 } 1407 CalcTrialShape(const FiniteElement & trial_fe,ElementTransformation & Trans,DenseMatrix & shape)1408 inline virtual void CalcTrialShape(const FiniteElement & trial_fe, 1409 ElementTransformation &Trans, 1410 DenseMatrix & shape) 1411 { trial_fe.CalcPhysDShape(Trans, shape); } 1412 CalcTestShape(const FiniteElement & test_fe,ElementTransformation & Trans,DenseMatrix & shape)1413 inline virtual void CalcTestShape(const FiniteElement & test_fe, 1414 ElementTransformation &Trans, 1415 DenseMatrix & shape) 1416 { test_fe.CalcVShape(Trans, shape); } 1417 }; 1418 1419 /** Class for integrating the bilinear form a(u,v) := (V x Curl u, v) in 3D and 1420 where V is a vector coefficient u is in H(Curl) and v is in H(Curl) or 1421 H(Div). */ 1422 class MixedCrossCurlIntegrator : public MixedVectorIntegrator 1423 { 1424 public: MixedCrossCurlIntegrator(VectorCoefficient & vq)1425 MixedCrossCurlIntegrator(VectorCoefficient &vq) 1426 : MixedVectorIntegrator(vq, false) {} 1427 VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const1428 inline virtual bool VerifyFiniteElementTypes( 1429 const FiniteElement & trial_fe, 1430 const FiniteElement & test_fe) const 1431 { 1432 return (trial_fe.GetDim() == 3 && test_fe.GetDim() == 3 && 1433 trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR && 1434 trial_fe.GetDerivType() == mfem::FiniteElement::CURL && 1435 test_fe.GetRangeType() == mfem::FiniteElement::VECTOR ); 1436 } 1437 FiniteElementTypeFailureMessage() const1438 inline virtual const char * FiniteElementTypeFailureMessage() const 1439 { 1440 return "MixedCrossCurlIntegrator: " 1441 "Trial space must be a vector field in 3D with a curl " 1442 "and the test space must be a vector field"; 1443 } 1444 CalcTrialShape(const FiniteElement & trial_fe,ElementTransformation & Trans,DenseMatrix & shape)1445 inline virtual void CalcTrialShape(const FiniteElement & trial_fe, 1446 ElementTransformation &Trans, 1447 DenseMatrix & shape) 1448 { trial_fe.CalcPhysCurlShape(Trans, shape); } 1449 }; 1450 1451 /** Class for integrating the bilinear form a(u,v) := (V x Curl u, v) in 2D and 1452 where V is a vector coefficient u is in H(Curl) and v is in H(Curl) or 1453 H(Div). */ 1454 class MixedScalarCrossCurlIntegrator : public MixedScalarVectorIntegrator 1455 { 1456 public: MixedScalarCrossCurlIntegrator(VectorCoefficient & vq)1457 MixedScalarCrossCurlIntegrator(VectorCoefficient &vq) 1458 : MixedScalarVectorIntegrator(vq, false, true) {} 1459 VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const1460 inline virtual bool VerifyFiniteElementTypes( 1461 const FiniteElement & trial_fe, 1462 const FiniteElement & test_fe) const 1463 { 1464 return (trial_fe.GetDim() == 2 && test_fe.GetDim() == 2 && 1465 trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR && 1466 trial_fe.GetDerivType() == mfem::FiniteElement::CURL && 1467 test_fe.GetRangeType() == mfem::FiniteElement::VECTOR ); 1468 } 1469 FiniteElementTypeFailureMessage() const1470 inline virtual const char * FiniteElementTypeFailureMessage() const 1471 { 1472 return "MixedCrossCurlIntegrator: " 1473 "Trial space must be a vector field in 2D with a curl " 1474 "and the test space must be a vector field"; 1475 } 1476 CalcShape(const FiniteElement & scalar_fe,ElementTransformation & Trans,Vector & shape)1477 inline virtual void CalcShape(const FiniteElement & scalar_fe, 1478 ElementTransformation &Trans, 1479 Vector & shape) 1480 { 1481 DenseMatrix dshape(shape.GetData(), shape.Size(), 1); 1482 scalar_fe.CalcPhysCurlShape(Trans, dshape); shape *= -1.0; 1483 } 1484 }; 1485 1486 /** Class for integrating the bilinear form a(u,v) := (V x Grad u, v) in 2D and 1487 where V is a vector coefficient u is in H1 and v is in H1 or L2. */ 1488 class MixedScalarCrossGradIntegrator : public MixedScalarVectorIntegrator 1489 { 1490 public: MixedScalarCrossGradIntegrator(VectorCoefficient & vq)1491 MixedScalarCrossGradIntegrator(VectorCoefficient &vq) 1492 : MixedScalarVectorIntegrator(vq, true, true) {} 1493 VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const1494 inline virtual bool VerifyFiniteElementTypes( 1495 const FiniteElement & trial_fe, 1496 const FiniteElement & test_fe) const 1497 { 1498 return (trial_fe.GetDim() == 2 && test_fe.GetDim() == 2 && 1499 trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR && 1500 trial_fe.GetDerivType() == mfem::FiniteElement::GRAD && 1501 test_fe.GetRangeType() == mfem::FiniteElement::SCALAR ); 1502 } 1503 FiniteElementTypeFailureMessage() const1504 inline virtual const char * FiniteElementTypeFailureMessage() const 1505 { 1506 return "MixedScalarCrossGradIntegrator: " 1507 "Trial space must be a scalar field in 2D with a gradient " 1508 "and the test space must be a scalar field"; 1509 } 1510 CalcVShape(const FiniteElement & vector_fe,ElementTransformation & Trans,DenseMatrix & shape)1511 inline virtual void CalcVShape(const FiniteElement & vector_fe, 1512 ElementTransformation &Trans, 1513 DenseMatrix & shape) 1514 { vector_fe.CalcPhysDShape(Trans, shape); } 1515 }; 1516 1517 /** Class for integrating the bilinear form a(u,v) := (V x u, v) in 2D and where 1518 V is a vector coefficient u is in ND or RT and v is in H1 or L2. */ 1519 class MixedScalarCrossProductIntegrator : public MixedScalarVectorIntegrator 1520 { 1521 public: MixedScalarCrossProductIntegrator(VectorCoefficient & vq)1522 MixedScalarCrossProductIntegrator(VectorCoefficient &vq) 1523 : MixedScalarVectorIntegrator(vq, true, true) {} 1524 VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const1525 inline virtual bool VerifyFiniteElementTypes( 1526 const FiniteElement & trial_fe, 1527 const FiniteElement & test_fe) const 1528 { 1529 return (trial_fe.GetDim() == 2 && test_fe.GetDim() == 2 && 1530 trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR && 1531 test_fe.GetRangeType() == mfem::FiniteElement::SCALAR ); 1532 } 1533 FiniteElementTypeFailureMessage() const1534 inline virtual const char * FiniteElementTypeFailureMessage() const 1535 { 1536 return "MixedScalarCrossProductIntegrator: " 1537 "Trial space must be a vector field in 2D " 1538 "and the test space must be a scalar field"; 1539 } 1540 }; 1541 1542 /** Class for integrating the bilinear form a(u,v) := (V x z u, v) in 2D and 1543 where V is a vector coefficient u is in H1 or L2 and v is in ND or RT. */ 1544 class MixedScalarWeakCrossProductIntegrator : public MixedScalarVectorIntegrator 1545 { 1546 public: MixedScalarWeakCrossProductIntegrator(VectorCoefficient & vq)1547 MixedScalarWeakCrossProductIntegrator(VectorCoefficient &vq) 1548 : MixedScalarVectorIntegrator(vq, false, true) {} 1549 VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const1550 inline virtual bool VerifyFiniteElementTypes( 1551 const FiniteElement & trial_fe, 1552 const FiniteElement & test_fe) const 1553 { 1554 return (trial_fe.GetDim() == 2 && test_fe.GetDim() == 2 && 1555 trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR && 1556 test_fe.GetRangeType() == mfem::FiniteElement::VECTOR ); 1557 } 1558 FiniteElementTypeFailureMessage() const1559 inline virtual const char * FiniteElementTypeFailureMessage() const 1560 { 1561 return "MixedScalarWeakCrossProductIntegrator: " 1562 "Trial space must be a scalar field in 2D " 1563 "and the test space must be a vector field"; 1564 } 1565 CalcShape(const FiniteElement & scalar_fe,ElementTransformation & Trans,Vector & shape)1566 inline virtual void CalcShape(const FiniteElement & scalar_fe, 1567 ElementTransformation &Trans, 1568 Vector & shape) 1569 { scalar_fe.CalcPhysShape(Trans, shape); shape *= -1.0; } 1570 }; 1571 1572 /** Class for integrating the bilinear form a(u,v) := (V . Grad u, v) in 2D or 1573 3D and where V is a vector coefficient, u is in H1 and v is in H1 or L2. */ 1574 class MixedDirectionalDerivativeIntegrator : public MixedScalarVectorIntegrator 1575 { 1576 public: MixedDirectionalDerivativeIntegrator(VectorCoefficient & vq)1577 MixedDirectionalDerivativeIntegrator(VectorCoefficient &vq) 1578 : MixedScalarVectorIntegrator(vq, true) {} 1579 VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const1580 inline virtual bool VerifyFiniteElementTypes( 1581 const FiniteElement & trial_fe, 1582 const FiniteElement & test_fe) const 1583 { 1584 return (trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR && 1585 trial_fe.GetDerivType() == mfem::FiniteElement::GRAD && 1586 test_fe.GetRangeType() == mfem::FiniteElement::SCALAR ); 1587 } 1588 FiniteElementTypeFailureMessage() const1589 inline virtual const char * FiniteElementTypeFailureMessage() const 1590 { 1591 return "MixedDirectionalDerivativeIntegrator: " 1592 "Trial space must be a scalar field with a gradient " 1593 "and the test space must be a scalar field"; 1594 } 1595 CalcVShape(const FiniteElement & vector_fe,ElementTransformation & Trans,DenseMatrix & shape)1596 inline virtual void CalcVShape(const FiniteElement & vector_fe, 1597 ElementTransformation &Trans, 1598 DenseMatrix & shape) 1599 { vector_fe.CalcPhysDShape(Trans, shape); } 1600 }; 1601 1602 /** Class for integrating the bilinear form a(u,v) := (-V . Grad u, Div v) in 2D 1603 or 3D and where V is a vector coefficient, u is in H1 and v is in RT. */ 1604 class MixedGradDivIntegrator : public MixedScalarVectorIntegrator 1605 { 1606 public: MixedGradDivIntegrator(VectorCoefficient & vq)1607 MixedGradDivIntegrator(VectorCoefficient &vq) 1608 : MixedScalarVectorIntegrator(vq, true) {} 1609 VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const1610 inline virtual bool VerifyFiniteElementTypes( 1611 const FiniteElement & trial_fe, 1612 const FiniteElement & test_fe) const 1613 { 1614 return (trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR && 1615 trial_fe.GetDerivType() == mfem::FiniteElement::GRAD && 1616 test_fe.GetRangeType() == mfem::FiniteElement::VECTOR && 1617 test_fe.GetDerivType() == mfem::FiniteElement::DIV ); 1618 } 1619 FiniteElementTypeFailureMessage() const1620 inline virtual const char * FiniteElementTypeFailureMessage() const 1621 { 1622 return "MixedGradDivIntegrator: " 1623 "Trial space must be a scalar field with a gradient" 1624 "and the test space must be a vector field with a divergence"; 1625 } 1626 CalcVShape(const FiniteElement & vector_fe,ElementTransformation & Trans,DenseMatrix & shape)1627 inline virtual void CalcVShape(const FiniteElement & vector_fe, 1628 ElementTransformation &Trans, 1629 DenseMatrix & shape) 1630 { vector_fe.CalcPhysDShape(Trans, shape); shape *= -1.0; } 1631 CalcShape(const FiniteElement & scalar_fe,ElementTransformation & Trans,Vector & shape)1632 inline virtual void CalcShape(const FiniteElement & scalar_fe, 1633 ElementTransformation &Trans, 1634 Vector & shape) 1635 { scalar_fe.CalcPhysDivShape(Trans, shape); } 1636 }; 1637 1638 /** Class for integrating the bilinear form a(u,v) := (-V Div u, Grad v) in 2D 1639 or 3D and where V is a vector coefficient, u is in RT and v is in H1. */ 1640 class MixedDivGradIntegrator : public MixedScalarVectorIntegrator 1641 { 1642 public: MixedDivGradIntegrator(VectorCoefficient & vq)1643 MixedDivGradIntegrator(VectorCoefficient &vq) 1644 : MixedScalarVectorIntegrator(vq, false) {} 1645 VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const1646 inline virtual bool VerifyFiniteElementTypes( 1647 const FiniteElement & trial_fe, 1648 const FiniteElement & test_fe) const 1649 { 1650 return (trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR && 1651 trial_fe.GetDerivType() == mfem::FiniteElement::DIV && 1652 test_fe.GetRangeType() == mfem::FiniteElement::SCALAR && 1653 test_fe.GetDerivType() == mfem::FiniteElement::GRAD 1654 ); 1655 } 1656 FiniteElementTypeFailureMessage() const1657 inline virtual const char * FiniteElementTypeFailureMessage() const 1658 { 1659 return "MixedDivGradIntegrator: " 1660 "Trial space must be a vector field with a divergence" 1661 "and the test space must be a scalar field with a gradient"; 1662 } 1663 CalcVShape(const FiniteElement & vector_fe,ElementTransformation & Trans,DenseMatrix & shape)1664 inline virtual void CalcVShape(const FiniteElement & vector_fe, 1665 ElementTransformation &Trans, 1666 DenseMatrix & shape) 1667 { vector_fe.CalcPhysDShape(Trans, shape); shape *= -1.0; } 1668 CalcShape(const FiniteElement & scalar_fe,ElementTransformation & Trans,Vector & shape)1669 inline virtual void CalcShape(const FiniteElement & scalar_fe, 1670 ElementTransformation &Trans, 1671 Vector & shape) 1672 { scalar_fe.CalcPhysDivShape(Trans, shape); } 1673 }; 1674 1675 /** Class for integrating the bilinear form a(u,v) := (-V u, Grad v) in 2D or 3D 1676 and where V is a vector coefficient, u is in H1 or L2 and v is in H1. */ 1677 class MixedScalarWeakDivergenceIntegrator : public MixedScalarVectorIntegrator 1678 { 1679 public: MixedScalarWeakDivergenceIntegrator(VectorCoefficient & vq)1680 MixedScalarWeakDivergenceIntegrator(VectorCoefficient &vq) 1681 : MixedScalarVectorIntegrator(vq, false) {} 1682 VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const1683 inline virtual bool VerifyFiniteElementTypes( 1684 const FiniteElement & trial_fe, 1685 const FiniteElement & test_fe) const 1686 { 1687 return (trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR && 1688 test_fe.GetRangeType() == mfem::FiniteElement::SCALAR && 1689 test_fe.GetDerivType() == mfem::FiniteElement::GRAD ); 1690 } 1691 FiniteElementTypeFailureMessage() const1692 inline virtual const char * FiniteElementTypeFailureMessage() const 1693 { 1694 return "MixedScalarWeakDivergenceIntegrator: " 1695 "Trial space must be a scalar field " 1696 "and the test space must be a scalar field with a gradient"; 1697 } 1698 CalcVShape(const FiniteElement & vector_fe,ElementTransformation & Trans,DenseMatrix & shape)1699 inline virtual void CalcVShape(const FiniteElement & vector_fe, 1700 ElementTransformation &Trans, 1701 DenseMatrix & shape) 1702 { vector_fe.CalcPhysDShape(Trans, shape); shape *= -1.0; } 1703 }; 1704 1705 /** Class for integrating the bilinear form a(u,v) := (Q grad u, v) in either 2D 1706 or 3D and where Q is an optional coefficient (of type scalar, matrix, or 1707 diagonal matrix) u is in H1 and v is in H(Curl) or H(Div). */ 1708 class MixedVectorGradientIntegrator : public MixedVectorIntegrator 1709 { 1710 public: MixedVectorGradientIntegrator()1711 MixedVectorGradientIntegrator() {} MixedVectorGradientIntegrator(Coefficient & q)1712 MixedVectorGradientIntegrator(Coefficient &q) 1713 : MixedVectorIntegrator(q) {} MixedVectorGradientIntegrator(DiagonalMatrixCoefficient & dq)1714 MixedVectorGradientIntegrator(DiagonalMatrixCoefficient &dq) 1715 : MixedVectorIntegrator(dq, true) {} MixedVectorGradientIntegrator(MatrixCoefficient & mq)1716 MixedVectorGradientIntegrator(MatrixCoefficient &mq) 1717 : MixedVectorIntegrator(mq) {} 1718 1719 protected: VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const1720 inline virtual bool VerifyFiniteElementTypes( 1721 const FiniteElement & trial_fe, 1722 const FiniteElement & test_fe) const 1723 { 1724 return (trial_fe.GetDerivType() == mfem::FiniteElement::GRAD && 1725 test_fe.GetRangeType() == mfem::FiniteElement::VECTOR ); 1726 } 1727 FiniteElementTypeFailureMessage() const1728 inline virtual const char * FiniteElementTypeFailureMessage() const 1729 { 1730 return "MixedVectorGradientIntegrator: " 1731 "Trial spaces must be H1 and the test space must be a " 1732 "vector field in 2D or 3D"; 1733 } 1734 CalcTrialShape(const FiniteElement & trial_fe,ElementTransformation & Trans,DenseMatrix & shape)1735 inline virtual void CalcTrialShape(const FiniteElement & trial_fe, 1736 ElementTransformation &Trans, 1737 DenseMatrix & shape) 1738 { 1739 trial_fe.CalcPhysDShape(Trans, shape); 1740 } 1741 1742 using BilinearFormIntegrator::AssemblePA; 1743 virtual void AssemblePA(const FiniteElementSpace &trial_fes, 1744 const FiniteElementSpace &test_fes); 1745 1746 virtual void AddMultPA(const Vector&, Vector&) const; 1747 1748 private: 1749 DenseMatrix Jinv; 1750 1751 // PA extension 1752 Vector pa_data; 1753 const DofToQuad *mapsO; ///< Not owned. DOF-to-quad map, open. 1754 const DofToQuad *mapsC; ///< Not owned. DOF-to-quad map, closed. 1755 const GeometricFactors *geom; ///< Not owned 1756 int dim, ne, dofs1D, quad1D; 1757 }; 1758 1759 /** Class for integrating the bilinear form a(u,v) := (Q curl u, v) in 3D and 1760 where Q is an optional coefficient (of type scalar, matrix, or diagonal 1761 matrix) u is in H(Curl) and v is in H(Div) or H(Curl). */ 1762 class MixedVectorCurlIntegrator : public MixedVectorIntegrator 1763 { 1764 public: MixedVectorCurlIntegrator()1765 MixedVectorCurlIntegrator() {} MixedVectorCurlIntegrator(Coefficient & q)1766 MixedVectorCurlIntegrator(Coefficient &q) 1767 : MixedVectorIntegrator(q) {} MixedVectorCurlIntegrator(DiagonalMatrixCoefficient & dq)1768 MixedVectorCurlIntegrator(DiagonalMatrixCoefficient &dq) 1769 : MixedVectorIntegrator(dq, true) {} MixedVectorCurlIntegrator(MatrixCoefficient & mq)1770 MixedVectorCurlIntegrator(MatrixCoefficient &mq) 1771 : MixedVectorIntegrator(mq) {} 1772 1773 protected: VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const1774 inline virtual bool VerifyFiniteElementTypes( 1775 const FiniteElement & trial_fe, 1776 const FiniteElement & test_fe) const 1777 { 1778 return (trial_fe.GetDim() == 3 && test_fe.GetDim() == 3 && 1779 trial_fe.GetDerivType() == mfem::FiniteElement::CURL && 1780 test_fe.GetRangeType() == mfem::FiniteElement::VECTOR ); 1781 } 1782 FiniteElementTypeFailureMessage() const1783 inline virtual const char * FiniteElementTypeFailureMessage() const 1784 { 1785 return "MixedVectorCurlIntegrator: " 1786 "Trial space must be H(Curl) and the test space must be a " 1787 "vector field in 3D"; 1788 } 1789 CalcTrialShape(const FiniteElement & trial_fe,ElementTransformation & Trans,DenseMatrix & shape)1790 inline virtual void CalcTrialShape(const FiniteElement & trial_fe, 1791 ElementTransformation &Trans, 1792 DenseMatrix & shape) 1793 { 1794 trial_fe.CalcPhysCurlShape(Trans, shape); 1795 } 1796 1797 using BilinearFormIntegrator::AssemblePA; 1798 virtual void AssemblePA(const FiniteElementSpace &trial_fes, 1799 const FiniteElementSpace &test_fes); 1800 1801 virtual void AddMultPA(const Vector&, Vector&) const; 1802 1803 private: 1804 // PA extension 1805 Vector pa_data; 1806 const DofToQuad *mapsO; ///< Not owned. DOF-to-quad map, open. 1807 const DofToQuad *mapsC; ///< Not owned. DOF-to-quad map, closed. 1808 const DofToQuad *mapsOtest; ///< Not owned. DOF-to-quad map, open. 1809 const DofToQuad *mapsCtest; ///< Not owned. DOF-to-quad map, closed. 1810 const GeometricFactors *geom; ///< Not owned 1811 int dim, ne, dofs1D, dofs1Dtest,quad1D, testType, trialType, coeffDim; 1812 }; 1813 1814 /** Class for integrating the bilinear form a(u,v) := (Q u, curl v) in 3D and 1815 where Q is an optional coefficient (of type scalar, matrix, or diagonal 1816 matrix) u is in H(Div) or H(Curl) and v is in H(Curl). */ 1817 class MixedVectorWeakCurlIntegrator : public MixedVectorIntegrator 1818 { 1819 public: MixedVectorWeakCurlIntegrator()1820 MixedVectorWeakCurlIntegrator() {} MixedVectorWeakCurlIntegrator(Coefficient & q)1821 MixedVectorWeakCurlIntegrator(Coefficient &q) 1822 : MixedVectorIntegrator(q) {} MixedVectorWeakCurlIntegrator(DiagonalMatrixCoefficient & dq)1823 MixedVectorWeakCurlIntegrator(DiagonalMatrixCoefficient &dq) 1824 : MixedVectorIntegrator(dq, true) {} MixedVectorWeakCurlIntegrator(MatrixCoefficient & mq)1825 MixedVectorWeakCurlIntegrator(MatrixCoefficient &mq) 1826 : MixedVectorIntegrator(mq) {} 1827 1828 protected: VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const1829 inline virtual bool VerifyFiniteElementTypes( 1830 const FiniteElement & trial_fe, 1831 const FiniteElement & test_fe) const 1832 { 1833 return (trial_fe.GetDim() == 3 && test_fe.GetDim() == 3 && 1834 trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR && 1835 test_fe.GetDerivType() == mfem::FiniteElement::CURL ); 1836 } 1837 FiniteElementTypeFailureMessage() const1838 inline virtual const char * FiniteElementTypeFailureMessage() const 1839 { 1840 return "MixedVectorWeakCurlIntegrator: " 1841 "Trial space must be vector field in 3D and the " 1842 "test space must be H(Curl)"; 1843 } 1844 CalcTestShape(const FiniteElement & test_fe,ElementTransformation & Trans,DenseMatrix & shape)1845 inline virtual void CalcTestShape(const FiniteElement & test_fe, 1846 ElementTransformation &Trans, 1847 DenseMatrix & shape) 1848 { 1849 test_fe.CalcPhysCurlShape(Trans, shape); 1850 } 1851 1852 using BilinearFormIntegrator::AssemblePA; 1853 virtual void AssemblePA(const FiniteElementSpace &trial_fes, 1854 const FiniteElementSpace &test_fes); 1855 1856 virtual void AddMultPA(const Vector&, Vector&) const; 1857 1858 private: 1859 // PA extension 1860 Vector pa_data; 1861 const DofToQuad *mapsO; ///< Not owned. DOF-to-quad map, open. 1862 const DofToQuad *mapsC; ///< Not owned. DOF-to-quad map, closed. 1863 const GeometricFactors *geom; ///< Not owned 1864 int dim, ne, dofs1D, quad1D, testType, trialType, coeffDim; 1865 }; 1866 1867 /** Class for integrating the bilinear form a(u,v) := - (Q u, grad v) in either 1868 2D or 3D and where Q is an optional coefficient (of type scalar, matrix, or 1869 diagonal matrix) u is in H(Div) or H(Curl) and v is in H1. */ 1870 class MixedVectorWeakDivergenceIntegrator : public MixedVectorIntegrator 1871 { 1872 public: MixedVectorWeakDivergenceIntegrator()1873 MixedVectorWeakDivergenceIntegrator() {} MixedVectorWeakDivergenceIntegrator(Coefficient & q)1874 MixedVectorWeakDivergenceIntegrator(Coefficient &q) 1875 : MixedVectorIntegrator(q) {} MixedVectorWeakDivergenceIntegrator(DiagonalMatrixCoefficient & dq)1876 MixedVectorWeakDivergenceIntegrator(DiagonalMatrixCoefficient &dq) 1877 : MixedVectorIntegrator(dq, true) {} MixedVectorWeakDivergenceIntegrator(MatrixCoefficient & mq)1878 MixedVectorWeakDivergenceIntegrator(MatrixCoefficient &mq) 1879 : MixedVectorIntegrator(mq) {} 1880 1881 protected: VerifyFiniteElementTypes(const FiniteElement & trial_fe,const FiniteElement & test_fe) const1882 inline virtual bool VerifyFiniteElementTypes( 1883 const FiniteElement & trial_fe, 1884 const FiniteElement & test_fe) const 1885 { 1886 return (trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR && 1887 test_fe.GetDerivType() == mfem::FiniteElement::GRAD ); 1888 } 1889 FiniteElementTypeFailureMessage() const1890 inline virtual const char * FiniteElementTypeFailureMessage() const 1891 { 1892 return "MixedVectorWeakDivergenceIntegrator: " 1893 "Trial space must be vector field and the " 1894 "test space must be H1"; 1895 } 1896 CalcTestShape(const FiniteElement & test_fe,ElementTransformation & Trans,DenseMatrix & shape)1897 inline virtual void CalcTestShape(const FiniteElement & test_fe, 1898 ElementTransformation &Trans, 1899 DenseMatrix & shape) 1900 { 1901 test_fe.CalcPhysDShape(Trans, shape); 1902 shape *= -1.0; 1903 } 1904 }; 1905 1906 /** Class for integrating the bilinear form a(u,v) := (Q grad u, v) where Q is a 1907 scalar coefficient, and v is a vector with components v_i in the same (H1) space 1908 as u. 1909 1910 See also MixedVectorGradientIntegrator when v is in H(curl). */ 1911 class GradientIntegrator : public BilinearFormIntegrator 1912 { 1913 protected: 1914 Coefficient *Q; 1915 1916 private: 1917 Vector shape; 1918 DenseMatrix dshape; 1919 DenseMatrix gshape; 1920 DenseMatrix Jadj; 1921 DenseMatrix elmat_comp; 1922 // PA extension 1923 Vector pa_data; 1924 const DofToQuad *trial_maps, *test_maps; ///< Not owned 1925 const GeometricFactors *geom; ///< Not owned 1926 int dim, ne, nq; 1927 int trial_dofs1D, test_dofs1D, quad1D; 1928 1929 public: GradientIntegrator()1930 GradientIntegrator() : 1931 Q{NULL}, trial_maps{NULL}, test_maps{NULL}, geom{NULL} 1932 { } GradientIntegrator(Coefficient * q_)1933 GradientIntegrator(Coefficient *q_) : 1934 Q{q_}, trial_maps{NULL}, test_maps{NULL}, geom{NULL} 1935 { } GradientIntegrator(Coefficient & q)1936 GradientIntegrator(Coefficient &q) : 1937 Q{&q}, trial_maps{NULL}, test_maps{NULL}, geom{NULL} 1938 { } 1939 1940 virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, 1941 const FiniteElement &test_fe, 1942 ElementTransformation &Trans, 1943 DenseMatrix &elmat); 1944 1945 using BilinearFormIntegrator::AssemblePA; 1946 virtual void AssemblePA(const FiniteElementSpace &trial_fes, 1947 const FiniteElementSpace &test_fes); 1948 1949 virtual void AddMultPA(const Vector &x, Vector &y) const; 1950 virtual void AddMultTransposePA(const Vector &x, Vector &y) const; 1951 1952 static const IntegrationRule &GetRule(const FiniteElement &trial_fe, 1953 const FiniteElement &test_fe, 1954 ElementTransformation &Trans); 1955 }; 1956 1957 /** Class for integrating the bilinear form a(u,v) := (Q grad u, grad v) where Q 1958 can be a scalar or a matrix coefficient. */ 1959 class DiffusionIntegrator: public BilinearFormIntegrator 1960 { 1961 protected: 1962 Coefficient *Q; 1963 VectorCoefficient *VQ; 1964 MatrixCoefficient *MQ; 1965 SymmetricMatrixCoefficient *SMQ; 1966 1967 private: 1968 Vector vec, vecdxt, pointflux, shape; 1969 #ifndef MFEM_THREAD_SAFE 1970 DenseMatrix dshape, dshapedxt, invdfdx, M, dshapedxt_m; 1971 DenseMatrix te_dshape, te_dshapedxt; 1972 Vector D; 1973 #endif 1974 1975 // PA extension 1976 const FiniteElementSpace *fespace; 1977 const DofToQuad *maps; ///< Not owned 1978 const GeometricFactors *geom; ///< Not owned 1979 int dim, ne, dofs1D, quad1D; 1980 Vector pa_data; 1981 bool symmetric = true; ///< False if using a nonsymmetric matrix coefficient 1982 1983 public: 1984 /// Construct a diffusion integrator with coefficient Q = 1 DiffusionIntegrator()1985 DiffusionIntegrator() 1986 : Q(NULL), VQ(NULL), MQ(NULL), SMQ(NULL), maps(NULL), geom(NULL) { } 1987 1988 /// Construct a diffusion integrator with a scalar coefficient q DiffusionIntegrator(Coefficient & q)1989 DiffusionIntegrator(Coefficient &q) 1990 : Q(&q), VQ(NULL), MQ(NULL), SMQ(NULL), maps(NULL), geom(NULL) { } 1991 1992 /// Construct a diffusion integrator with a vector coefficient q DiffusionIntegrator(VectorCoefficient & q)1993 DiffusionIntegrator(VectorCoefficient &q) 1994 : Q(NULL), VQ(&q), MQ(NULL), SMQ(NULL), maps(NULL), geom(NULL) { } 1995 1996 /// Construct a diffusion integrator with a matrix coefficient q DiffusionIntegrator(MatrixCoefficient & q)1997 DiffusionIntegrator(MatrixCoefficient &q) 1998 : Q(NULL), VQ(NULL), MQ(&q), SMQ(NULL), maps(NULL), geom(NULL) { } 1999 2000 /// Construct a diffusion integrator with a symmetric matrix coefficient q DiffusionIntegrator(SymmetricMatrixCoefficient & q)2001 DiffusionIntegrator(SymmetricMatrixCoefficient &q) 2002 : Q(NULL), VQ(NULL), MQ(NULL), SMQ(&q), maps(NULL), geom(NULL) { } 2003 2004 /** Given a particular Finite Element computes the element stiffness matrix 2005 elmat. */ 2006 virtual void AssembleElementMatrix(const FiniteElement &el, 2007 ElementTransformation &Trans, 2008 DenseMatrix &elmat); 2009 /** Given a trial and test Finite Element computes the element stiffness 2010 matrix elmat. */ 2011 virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, 2012 const FiniteElement &test_fe, 2013 ElementTransformation &Trans, 2014 DenseMatrix &elmat); 2015 2016 /// Perform the local action of the BilinearFormIntegrator 2017 virtual void AssembleElementVector(const FiniteElement &el, 2018 ElementTransformation &Tr, 2019 const Vector &elfun, Vector &elvect); 2020 2021 virtual void ComputeElementFlux(const FiniteElement &el, 2022 ElementTransformation &Trans, 2023 Vector &u, const FiniteElement &fluxelem, 2024 Vector &flux, bool with_coef = true); 2025 2026 virtual double ComputeFluxEnergy(const FiniteElement &fluxelem, 2027 ElementTransformation &Trans, 2028 Vector &flux, Vector *d_energy = NULL); 2029 2030 using BilinearFormIntegrator::AssemblePA; 2031 2032 virtual void AssembleMF(const FiniteElementSpace &fes); 2033 2034 virtual void AssemblePA(const FiniteElementSpace &fes); 2035 2036 virtual void AssembleEA(const FiniteElementSpace &fes, Vector &emat, 2037 const bool add); 2038 2039 virtual void AssembleDiagonalPA(Vector &diag); 2040 2041 virtual void AssembleDiagonalMF(Vector &diag); 2042 2043 virtual void AddMultMF(const Vector&, Vector&) const; 2044 2045 virtual void AddMultPA(const Vector&, Vector&) const; 2046 2047 virtual void AddMultTransposePA(const Vector&, Vector&) const; 2048 2049 static const IntegrationRule &GetRule(const FiniteElement &trial_fe, 2050 const FiniteElement &test_fe); 2051 SupportsCeed() const2052 bool SupportsCeed() const { return DeviceCanUseCeed(); } 2053 }; 2054 2055 /** Class for local mass matrix assembling a(u,v) := (Q u, v) */ 2056 class MassIntegrator: public BilinearFormIntegrator 2057 { 2058 protected: 2059 #ifndef MFEM_THREAD_SAFE 2060 Vector shape, te_shape; 2061 #endif 2062 Coefficient *Q; 2063 // PA extension 2064 const FiniteElementSpace *fespace; 2065 Vector pa_data; 2066 const DofToQuad *maps; ///< Not owned 2067 const GeometricFactors *geom; ///< Not owned 2068 int dim, ne, nq, dofs1D, quad1D; 2069 2070 public: MassIntegrator(const IntegrationRule * ir=NULL)2071 MassIntegrator(const IntegrationRule *ir = NULL) 2072 : BilinearFormIntegrator(ir), Q(NULL), maps(NULL), geom(NULL) { } 2073 2074 /// Construct a mass integrator with coefficient q MassIntegrator(Coefficient & q,const IntegrationRule * ir=NULL)2075 MassIntegrator(Coefficient &q, const IntegrationRule *ir = NULL) 2076 : BilinearFormIntegrator(ir), Q(&q), maps(NULL), geom(NULL) { } 2077 2078 /** Given a particular Finite Element computes the element mass matrix 2079 elmat. */ 2080 virtual void AssembleElementMatrix(const FiniteElement &el, 2081 ElementTransformation &Trans, 2082 DenseMatrix &elmat); 2083 virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, 2084 const FiniteElement &test_fe, 2085 ElementTransformation &Trans, 2086 DenseMatrix &elmat); 2087 2088 using BilinearFormIntegrator::AssemblePA; 2089 2090 virtual void AssembleMF(const FiniteElementSpace &fes); 2091 2092 virtual void AssemblePA(const FiniteElementSpace &fes); 2093 2094 virtual void AssembleEA(const FiniteElementSpace &fes, Vector &emat, 2095 const bool add); 2096 2097 virtual void AssembleDiagonalPA(Vector &diag); 2098 2099 virtual void AssembleDiagonalMF(Vector &diag); 2100 2101 virtual void AddMultMF(const Vector&, Vector&) const; 2102 2103 virtual void AddMultPA(const Vector&, Vector&) const; 2104 2105 virtual void AddMultTransposePA(const Vector&, Vector&) const; 2106 2107 static const IntegrationRule &GetRule(const FiniteElement &trial_fe, 2108 const FiniteElement &test_fe, 2109 ElementTransformation &Trans); 2110 SupportsCeed() const2111 bool SupportsCeed() const { return DeviceCanUseCeed(); } 2112 }; 2113 2114 /** Mass integrator (u, v) restricted to the boundary of a domain */ 2115 class BoundaryMassIntegrator : public MassIntegrator 2116 { 2117 public: BoundaryMassIntegrator(Coefficient & q)2118 BoundaryMassIntegrator(Coefficient &q) : MassIntegrator(q) { } 2119 2120 using BilinearFormIntegrator::AssembleFaceMatrix; 2121 2122 virtual void AssembleFaceMatrix(const FiniteElement &el1, 2123 const FiniteElement &el2, 2124 FaceElementTransformations &Trans, 2125 DenseMatrix &elmat); 2126 }; 2127 2128 /// alpha (q . grad u, v) 2129 class ConvectionIntegrator : public BilinearFormIntegrator 2130 { 2131 protected: 2132 VectorCoefficient *Q; 2133 double alpha; 2134 // PA extension 2135 Vector pa_data; 2136 const DofToQuad *maps; ///< Not owned 2137 const GeometricFactors *geom; ///< Not owned 2138 int dim, ne, nq, dofs1D, quad1D; 2139 2140 private: 2141 #ifndef MFEM_THREAD_SAFE 2142 DenseMatrix dshape, adjJ, Q_ir; 2143 Vector shape, vec2, BdFidxT; 2144 #endif 2145 2146 public: ConvectionIntegrator(VectorCoefficient & q,double a=1.0)2147 ConvectionIntegrator(VectorCoefficient &q, double a = 1.0) 2148 : Q(&q) { alpha = a; } 2149 2150 virtual void AssembleElementMatrix(const FiniteElement &, 2151 ElementTransformation &, 2152 DenseMatrix &); 2153 2154 using BilinearFormIntegrator::AssemblePA; 2155 2156 virtual void AssembleMF(const FiniteElementSpace &fes); 2157 2158 virtual void AssemblePA(const FiniteElementSpace&); 2159 2160 virtual void AssembleEA(const FiniteElementSpace &fes, Vector &emat, 2161 const bool add); 2162 2163 virtual void AssembleDiagonalPA(Vector &diag); 2164 2165 virtual void AssembleDiagonalMF(Vector &diag); 2166 2167 virtual void AddMultMF(const Vector&, Vector&) const; 2168 2169 virtual void AddMultPA(const Vector&, Vector&) const; 2170 2171 static const IntegrationRule &GetRule(const FiniteElement &el, 2172 ElementTransformation &Trans); 2173 2174 static const IntegrationRule &GetRule(const FiniteElement &trial_fe, 2175 const FiniteElement &test_fe, 2176 ElementTransformation &Trans); 2177 SupportsCeed() const2178 bool SupportsCeed() const { return DeviceCanUseCeed(); } 2179 }; 2180 2181 // Alias for @ConvectionIntegrator. 2182 using NonconservativeConvectionIntegrator = ConvectionIntegrator; 2183 2184 /// -alpha (u, q . grad v), negative transpose of ConvectionIntegrator 2185 class ConservativeConvectionIntegrator : public TransposeIntegrator 2186 { 2187 public: ConservativeConvectionIntegrator(VectorCoefficient & q,double a=1.0)2188 ConservativeConvectionIntegrator(VectorCoefficient &q, double a = 1.0) 2189 : TransposeIntegrator(new ConvectionIntegrator(q, -a)) { } 2190 }; 2191 2192 /// alpha (q . grad u, v) using the "group" FE discretization 2193 class GroupConvectionIntegrator : public BilinearFormIntegrator 2194 { 2195 protected: 2196 VectorCoefficient *Q; 2197 double alpha; 2198 2199 private: 2200 DenseMatrix dshape, adjJ, Q_nodal, grad; 2201 Vector shape; 2202 2203 public: GroupConvectionIntegrator(VectorCoefficient & q,double a=1.0)2204 GroupConvectionIntegrator(VectorCoefficient &q, double a = 1.0) 2205 : Q(&q) { alpha = a; } 2206 virtual void AssembleElementMatrix(const FiniteElement &, 2207 ElementTransformation &, 2208 DenseMatrix &); 2209 }; 2210 2211 /** Class for integrating the bilinear form a(u,v) := (Q u, v), 2212 where u=(u1,...,un) and v=(v1,...,vn); ui and vi are defined 2213 by scalar FE through standard transformation. */ 2214 class VectorMassIntegrator: public BilinearFormIntegrator 2215 { 2216 private: 2217 int vdim; 2218 Vector shape, te_shape, vec; 2219 DenseMatrix partelmat; 2220 DenseMatrix mcoeff; 2221 int Q_order; 2222 2223 protected: 2224 Coefficient *Q; 2225 VectorCoefficient *VQ; 2226 MatrixCoefficient *MQ; 2227 // PA extension 2228 Vector pa_data; 2229 const DofToQuad *maps; ///< Not owned 2230 const GeometricFactors *geom; ///< Not owned 2231 int dim, ne, nq, dofs1D, quad1D; 2232 2233 public: 2234 /// Construct an integrator with coefficient 1.0 VectorMassIntegrator()2235 VectorMassIntegrator() 2236 : vdim(-1), Q_order(0), Q(NULL), VQ(NULL), MQ(NULL) { } 2237 /** Construct an integrator with scalar coefficient q. If possible, save 2238 memory by using a scalar integrator since the resulting matrix is block 2239 diagonal with the same diagonal block repeated. */ VectorMassIntegrator(Coefficient & q,int qo=0)2240 VectorMassIntegrator(Coefficient &q, int qo = 0) 2241 : vdim(-1), Q_order(qo), Q(&q), VQ(NULL), MQ(NULL) { } VectorMassIntegrator(Coefficient & q,const IntegrationRule * ir)2242 VectorMassIntegrator(Coefficient &q, const IntegrationRule *ir) 2243 : BilinearFormIntegrator(ir), vdim(-1), Q_order(0), Q(&q), VQ(NULL), 2244 MQ(NULL) { } 2245 /// Construct an integrator with diagonal coefficient q VectorMassIntegrator(VectorCoefficient & q,int qo=0)2246 VectorMassIntegrator(VectorCoefficient &q, int qo = 0) 2247 : vdim(q.GetVDim()), Q_order(qo), Q(NULL), VQ(&q), MQ(NULL) { } 2248 /// Construct an integrator with matrix coefficient q VectorMassIntegrator(MatrixCoefficient & q,int qo=0)2249 VectorMassIntegrator(MatrixCoefficient &q, int qo = 0) 2250 : vdim(q.GetVDim()), Q_order(qo), Q(NULL), VQ(NULL), MQ(&q) { } 2251 GetVDim() const2252 int GetVDim() const { return vdim; } SetVDim(int vdim)2253 void SetVDim(int vdim) { this->vdim = vdim; } 2254 2255 virtual void AssembleElementMatrix(const FiniteElement &el, 2256 ElementTransformation &Trans, 2257 DenseMatrix &elmat); 2258 virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, 2259 const FiniteElement &test_fe, 2260 ElementTransformation &Trans, 2261 DenseMatrix &elmat); 2262 using BilinearFormIntegrator::AssemblePA; 2263 virtual void AssemblePA(const FiniteElementSpace &fes); 2264 virtual void AssembleMF(const FiniteElementSpace &fes); 2265 virtual void AssembleDiagonalPA(Vector &diag); 2266 virtual void AssembleDiagonalMF(Vector &diag); 2267 virtual void AddMultPA(const Vector &x, Vector &y) const; 2268 virtual void AddMultMF(const Vector &x, Vector &y) const; SupportsCeed() const2269 bool SupportsCeed() const { return DeviceCanUseCeed(); } 2270 }; 2271 2272 2273 /** Class for integrating (div u, p) where u is a vector field given by 2274 VectorFiniteElement through Piola transformation (for RT elements); p is 2275 scalar function given by FiniteElement through standard transformation. 2276 Here, u is the trial function and p is the test function. 2277 2278 Note: the element matrix returned by AssembleElementMatrix2 does NOT depend 2279 on the ElementTransformation Trans. */ 2280 class VectorFEDivergenceIntegrator : public BilinearFormIntegrator 2281 { 2282 protected: 2283 Coefficient *Q; 2284 2285 using BilinearFormIntegrator::AssemblePA; 2286 virtual void AssemblePA(const FiniteElementSpace &trial_fes, 2287 const FiniteElementSpace &test_fes); 2288 2289 virtual void AddMultPA(const Vector&, Vector&) const; 2290 virtual void AddMultTransposePA(const Vector&, Vector&) const; 2291 2292 private: 2293 #ifndef MFEM_THREAD_SAFE 2294 Vector divshape, shape; 2295 #endif 2296 2297 // PA extension 2298 Vector pa_data; 2299 const DofToQuad *mapsO; ///< Not owned. DOF-to-quad map, open. 2300 const DofToQuad *L2mapsO; ///< Not owned. DOF-to-quad map, open. 2301 const DofToQuad *mapsC; ///< Not owned. DOF-to-quad map, closed. 2302 int dim, ne, dofs1D, L2dofs1D, quad1D; 2303 2304 public: VectorFEDivergenceIntegrator()2305 VectorFEDivergenceIntegrator() { Q = NULL; } VectorFEDivergenceIntegrator(Coefficient & q)2306 VectorFEDivergenceIntegrator(Coefficient &q) { Q = &q; } AssembleElementMatrix(const FiniteElement & el,ElementTransformation & Trans,DenseMatrix & elmat)2307 virtual void AssembleElementMatrix(const FiniteElement &el, 2308 ElementTransformation &Trans, 2309 DenseMatrix &elmat) { } 2310 virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, 2311 const FiniteElement &test_fe, 2312 ElementTransformation &Trans, 2313 DenseMatrix &elmat); 2314 2315 virtual void AssembleDiagonalPA_ADAt(const Vector &D, Vector &diag); 2316 }; 2317 2318 2319 /** Integrator for `(-Q u, grad v)` for Nedelec (`u`) and H1 (`v`) elements. 2320 This is equivalent to a weak divergence of the Nedelec basis functions. */ 2321 class VectorFEWeakDivergenceIntegrator: public BilinearFormIntegrator 2322 { 2323 protected: 2324 Coefficient *Q; 2325 2326 private: 2327 #ifndef MFEM_THREAD_SAFE 2328 DenseMatrix dshape; 2329 DenseMatrix dshapedxt; 2330 DenseMatrix vshape; 2331 DenseMatrix invdfdx; 2332 #endif 2333 2334 public: VectorFEWeakDivergenceIntegrator()2335 VectorFEWeakDivergenceIntegrator() { Q = NULL; } VectorFEWeakDivergenceIntegrator(Coefficient & q)2336 VectorFEWeakDivergenceIntegrator(Coefficient &q) { Q = &q; } AssembleElementMatrix(const FiniteElement & el,ElementTransformation & Trans,DenseMatrix & elmat)2337 virtual void AssembleElementMatrix(const FiniteElement &el, 2338 ElementTransformation &Trans, 2339 DenseMatrix &elmat) { } 2340 virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, 2341 const FiniteElement &test_fe, 2342 ElementTransformation &Trans, 2343 DenseMatrix &elmat); 2344 }; 2345 2346 /** Integrator for (curl u, v) for Nedelec and RT elements. If the trial and 2347 test spaces are switched, assembles the form (u, curl v). */ 2348 class VectorFECurlIntegrator: public BilinearFormIntegrator 2349 { 2350 protected: 2351 Coefficient *Q; 2352 2353 private: 2354 #ifndef MFEM_THREAD_SAFE 2355 DenseMatrix curlshapeTrial; 2356 DenseMatrix vshapeTest; 2357 DenseMatrix curlshapeTrial_dFT; 2358 #endif 2359 2360 public: VectorFECurlIntegrator()2361 VectorFECurlIntegrator() { Q = NULL; } VectorFECurlIntegrator(Coefficient & q)2362 VectorFECurlIntegrator(Coefficient &q) { Q = &q; } AssembleElementMatrix(const FiniteElement & el,ElementTransformation & Trans,DenseMatrix & elmat)2363 virtual void AssembleElementMatrix(const FiniteElement &el, 2364 ElementTransformation &Trans, 2365 DenseMatrix &elmat) { } 2366 virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, 2367 const FiniteElement &test_fe, 2368 ElementTransformation &Trans, 2369 DenseMatrix &elmat); 2370 }; 2371 2372 /// Class for integrating (Q D_i(u), v); u and v are scalars 2373 class DerivativeIntegrator : public BilinearFormIntegrator 2374 { 2375 protected: 2376 Coefficient* Q; 2377 2378 private: 2379 int xi; 2380 DenseMatrix dshape, dshapedxt, invdfdx; 2381 Vector shape, dshapedxi; 2382 2383 public: DerivativeIntegrator(Coefficient & q,int i)2384 DerivativeIntegrator(Coefficient &q, int i) : Q(&q), xi(i) { } AssembleElementMatrix(const FiniteElement & el,ElementTransformation & Trans,DenseMatrix & elmat)2385 virtual void AssembleElementMatrix(const FiniteElement &el, 2386 ElementTransformation &Trans, 2387 DenseMatrix &elmat) 2388 { AssembleElementMatrix2(el,el,Trans,elmat); } 2389 virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, 2390 const FiniteElement &test_fe, 2391 ElementTransformation &Trans, 2392 DenseMatrix &elmat); 2393 }; 2394 2395 /// Integrator for (curl u, curl v) for Nedelec elements 2396 class CurlCurlIntegrator: public BilinearFormIntegrator 2397 { 2398 private: 2399 Vector vec, pointflux; 2400 #ifndef MFEM_THREAD_SAFE 2401 Vector D; 2402 DenseMatrix curlshape, curlshape_dFt, M; 2403 DenseMatrix vshape, projcurl; 2404 #endif 2405 2406 protected: 2407 Coefficient *Q; 2408 DiagonalMatrixCoefficient *DQ; 2409 MatrixCoefficient *MQ; 2410 SymmetricMatrixCoefficient *SMQ; 2411 2412 // PA extension 2413 Vector pa_data; 2414 const DofToQuad *mapsO; ///< Not owned. DOF-to-quad map, open. 2415 const DofToQuad *mapsC; ///< Not owned. DOF-to-quad map, closed. 2416 const GeometricFactors *geom; ///< Not owned 2417 int dim, ne, nq, dofs1D, quad1D; 2418 bool symmetric = true; ///< False if using a nonsymmetric matrix coefficient 2419 2420 public: CurlCurlIntegrator()2421 CurlCurlIntegrator() { Q = NULL; DQ = NULL; MQ = NULL; SMQ = NULL; } 2422 /// Construct a bilinear form integrator for Nedelec elements CurlCurlIntegrator(Coefficient & q,const IntegrationRule * ir=NULL)2423 CurlCurlIntegrator(Coefficient &q, const IntegrationRule *ir = NULL) : 2424 BilinearFormIntegrator(ir), Q(&q), DQ(NULL), MQ(NULL), SMQ(NULL) { } CurlCurlIntegrator(DiagonalMatrixCoefficient & dq,const IntegrationRule * ir=NULL)2425 CurlCurlIntegrator(DiagonalMatrixCoefficient &dq, 2426 const IntegrationRule *ir = NULL) : 2427 BilinearFormIntegrator(ir), Q(NULL), DQ(&dq), MQ(NULL), SMQ(NULL) { } CurlCurlIntegrator(MatrixCoefficient & mq,const IntegrationRule * ir=NULL)2428 CurlCurlIntegrator(MatrixCoefficient &mq, const IntegrationRule *ir = NULL) : 2429 BilinearFormIntegrator(ir), Q(NULL), DQ(NULL), MQ(&mq), SMQ(NULL) { } CurlCurlIntegrator(SymmetricMatrixCoefficient & smq,const IntegrationRule * ir=NULL)2430 CurlCurlIntegrator(SymmetricMatrixCoefficient &smq, 2431 const IntegrationRule *ir = NULL) : 2432 BilinearFormIntegrator(ir), Q(NULL), DQ(NULL), MQ(NULL), SMQ(&smq) { } 2433 2434 /* Given a particular Finite Element, compute the 2435 element curl-curl matrix elmat */ 2436 virtual void AssembleElementMatrix(const FiniteElement &el, 2437 ElementTransformation &Trans, 2438 DenseMatrix &elmat); 2439 2440 virtual void ComputeElementFlux(const FiniteElement &el, 2441 ElementTransformation &Trans, 2442 Vector &u, const FiniteElement &fluxelem, 2443 Vector &flux, bool with_coef); 2444 2445 virtual double ComputeFluxEnergy(const FiniteElement &fluxelem, 2446 ElementTransformation &Trans, 2447 Vector &flux, Vector *d_energy = NULL); 2448 2449 using BilinearFormIntegrator::AssemblePA; 2450 virtual void AssemblePA(const FiniteElementSpace &fes); 2451 virtual void AddMultPA(const Vector &x, Vector &y) const; 2452 virtual void AssembleDiagonalPA(Vector& diag); 2453 }; 2454 2455 /** Integrator for (curl u, curl v) for FE spaces defined by 'dim' copies of a 2456 scalar FE space. */ 2457 class VectorCurlCurlIntegrator: public BilinearFormIntegrator 2458 { 2459 private: 2460 #ifndef MFEM_THREAD_SAFE 2461 DenseMatrix dshape_hat, dshape, curlshape, Jadj, grad_hat, grad; 2462 #endif 2463 2464 protected: 2465 Coefficient *Q; 2466 2467 public: VectorCurlCurlIntegrator()2468 VectorCurlCurlIntegrator() { Q = NULL; } 2469 VectorCurlCurlIntegrator(Coefficient & q)2470 VectorCurlCurlIntegrator(Coefficient &q) : Q(&q) { } 2471 2472 /// Assemble an element matrix 2473 virtual void AssembleElementMatrix(const FiniteElement &el, 2474 ElementTransformation &Trans, 2475 DenseMatrix &elmat); 2476 /// Compute element energy: (1/2) (curl u, curl u)_E 2477 virtual double GetElementEnergy(const FiniteElement &el, 2478 ElementTransformation &Tr, 2479 const Vector &elfun); 2480 }; 2481 2482 /** Integrator for (Q u, v), where Q is an optional coefficient (of type scalar, 2483 vector (diagonal matrix), or matrix), trial function u is in H(Curl) or 2484 H(Div), and test function v is in H(Curl), H(Div), or v=(v1,...,vn), where 2485 vi are in H1. */ 2486 class VectorFEMassIntegrator: public BilinearFormIntegrator 2487 { 2488 private: Init(Coefficient * q,DiagonalMatrixCoefficient * dq,MatrixCoefficient * mq,SymmetricMatrixCoefficient * smq)2489 void Init(Coefficient *q, DiagonalMatrixCoefficient *dq, MatrixCoefficient *mq, 2490 SymmetricMatrixCoefficient *smq) 2491 { Q = q; DQ = dq; MQ = mq; SMQ = smq; } 2492 2493 #ifndef MFEM_THREAD_SAFE 2494 Vector shape; 2495 Vector D; 2496 DenseMatrix K; 2497 DenseMatrix partelmat; 2498 DenseMatrix test_vshape; 2499 DenseMatrix trial_vshape; 2500 #endif 2501 2502 protected: 2503 Coefficient *Q; 2504 DiagonalMatrixCoefficient *DQ; 2505 MatrixCoefficient *MQ; 2506 SymmetricMatrixCoefficient *SMQ; 2507 2508 // PA extension 2509 Vector pa_data; 2510 const DofToQuad *mapsO; ///< Not owned. DOF-to-quad map, open. 2511 const DofToQuad *mapsC; ///< Not owned. DOF-to-quad map, closed. 2512 const DofToQuad *mapsOtest; ///< Not owned. DOF-to-quad map, open. 2513 const DofToQuad *mapsCtest; ///< Not owned. DOF-to-quad map, closed. 2514 const GeometricFactors *geom; ///< Not owned 2515 int dim, ne, nq, dofs1D, dofs1Dtest, quad1D, trial_fetype, test_fetype; 2516 bool symmetric = true; ///< False if using a nonsymmetric matrix coefficient 2517 2518 public: VectorFEMassIntegrator()2519 VectorFEMassIntegrator() { Init(NULL, NULL, NULL, NULL); } VectorFEMassIntegrator(Coefficient * q_)2520 VectorFEMassIntegrator(Coefficient *q_) { Init(q_, NULL, NULL, NULL); } VectorFEMassIntegrator(Coefficient & q)2521 VectorFEMassIntegrator(Coefficient &q) { Init(&q, NULL, NULL, NULL); } VectorFEMassIntegrator(DiagonalMatrixCoefficient * dq_)2522 VectorFEMassIntegrator(DiagonalMatrixCoefficient *dq_) { Init(NULL, dq_, NULL, NULL); } VectorFEMassIntegrator(DiagonalMatrixCoefficient & dq)2523 VectorFEMassIntegrator(DiagonalMatrixCoefficient &dq) { Init(NULL, &dq, NULL, NULL); } VectorFEMassIntegrator(MatrixCoefficient * mq_)2524 VectorFEMassIntegrator(MatrixCoefficient *mq_) { Init(NULL, NULL, mq_, NULL); } VectorFEMassIntegrator(MatrixCoefficient & mq)2525 VectorFEMassIntegrator(MatrixCoefficient &mq) { Init(NULL, NULL, &mq, NULL); } VectorFEMassIntegrator(SymmetricMatrixCoefficient & smq)2526 VectorFEMassIntegrator(SymmetricMatrixCoefficient &smq) { Init(NULL, NULL, NULL, &smq); } VectorFEMassIntegrator(SymmetricMatrixCoefficient * smq)2527 VectorFEMassIntegrator(SymmetricMatrixCoefficient *smq) { Init(NULL, NULL, NULL, smq); } 2528 2529 virtual void AssembleElementMatrix(const FiniteElement &el, 2530 ElementTransformation &Trans, 2531 DenseMatrix &elmat); 2532 virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, 2533 const FiniteElement &test_fe, 2534 ElementTransformation &Trans, 2535 DenseMatrix &elmat); 2536 2537 using BilinearFormIntegrator::AssemblePA; 2538 virtual void AssemblePA(const FiniteElementSpace &fes); 2539 virtual void AssemblePA(const FiniteElementSpace &trial_fes, 2540 const FiniteElementSpace &test_fes); 2541 virtual void AddMultPA(const Vector &x, Vector &y) const; 2542 virtual void AssembleDiagonalPA(Vector& diag); 2543 }; 2544 2545 /** Integrator for (Q div u, p) where u=(v1,...,vn) and all vi are in the same 2546 scalar FE space; p is also in a (different) scalar FE space. */ 2547 class VectorDivergenceIntegrator : public BilinearFormIntegrator 2548 { 2549 protected: 2550 Coefficient *Q; 2551 2552 private: 2553 Vector shape; 2554 Vector divshape; 2555 DenseMatrix dshape; 2556 DenseMatrix gshape; 2557 DenseMatrix Jadj; 2558 // PA extension 2559 Vector pa_data; 2560 const DofToQuad *trial_maps, *test_maps; ///< Not owned 2561 const GeometricFactors *geom; ///< Not owned 2562 int dim, ne, nq; 2563 int trial_dofs1D, test_dofs1D, quad1D; 2564 2565 public: VectorDivergenceIntegrator()2566 VectorDivergenceIntegrator() : 2567 Q(NULL), trial_maps(NULL), test_maps(NULL), geom(NULL) 2568 { } VectorDivergenceIntegrator(Coefficient * q_)2569 VectorDivergenceIntegrator(Coefficient *q_) : 2570 Q(q_), trial_maps(NULL), test_maps(NULL), geom(NULL) 2571 { } VectorDivergenceIntegrator(Coefficient & q)2572 VectorDivergenceIntegrator(Coefficient &q) : 2573 Q(&q), trial_maps(NULL), test_maps(NULL), geom(NULL) 2574 { } 2575 2576 virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, 2577 const FiniteElement &test_fe, 2578 ElementTransformation &Trans, 2579 DenseMatrix &elmat); 2580 2581 using BilinearFormIntegrator::AssemblePA; 2582 virtual void AssemblePA(const FiniteElementSpace &trial_fes, 2583 const FiniteElementSpace &test_fes); 2584 2585 virtual void AddMultPA(const Vector &x, Vector &y) const; 2586 virtual void AddMultTransposePA(const Vector &x, Vector &y) const; 2587 2588 static const IntegrationRule &GetRule(const FiniteElement &trial_fe, 2589 const FiniteElement &test_fe, 2590 ElementTransformation &Trans); 2591 }; 2592 2593 /// (Q div u, div v) for RT elements 2594 class DivDivIntegrator: public BilinearFormIntegrator 2595 { 2596 protected: 2597 Coefficient *Q; 2598 2599 using BilinearFormIntegrator::AssemblePA; 2600 virtual void AssemblePA(const FiniteElementSpace &fes); 2601 virtual void AddMultPA(const Vector &x, Vector &y) const; 2602 virtual void AssembleDiagonalPA(Vector& diag); 2603 2604 private: 2605 #ifndef MFEM_THREAD_SAFE 2606 Vector divshape; 2607 #endif 2608 2609 // PA extension 2610 Vector pa_data; 2611 const DofToQuad *mapsO; ///< Not owned. DOF-to-quad map, open. 2612 const DofToQuad *mapsC; ///< Not owned. DOF-to-quad map, closed. 2613 const GeometricFactors *geom; ///< Not owned 2614 int dim, ne, dofs1D, quad1D; 2615 2616 public: DivDivIntegrator()2617 DivDivIntegrator() { Q = NULL; } DivDivIntegrator(Coefficient & q)2618 DivDivIntegrator(Coefficient &q) : Q(&q) { } 2619 2620 virtual void AssembleElementMatrix(const FiniteElement &el, 2621 ElementTransformation &Trans, 2622 DenseMatrix &elmat); 2623 }; 2624 2625 /** Integrator for 2626 2627 (Q grad u, grad v) = sum_i (Q grad u_i, grad v_i) e_i e_i^T 2628 2629 for vector FE spaces, where e_i is the unit vector in the i-th direction. 2630 The resulting local element matrix is square, of size <tt> vdim*dof </tt>, 2631 where \c vdim is the vector dimension space and \c dof is the local degrees 2632 of freedom. The integrator is not aware of the true vector dimension and 2633 must use \c VectorCoefficient, \c MatrixCoefficient, or a caller-specified 2634 value to determine the vector space. For a scalar coefficient, the caller 2635 may manually specify the vector dimension or the vector dimension is assumed 2636 to be the spatial dimension (i.e. 2-dimension or 3-dimension). 2637 */ 2638 class VectorDiffusionIntegrator : public BilinearFormIntegrator 2639 { 2640 protected: 2641 Coefficient *Q = NULL; 2642 VectorCoefficient *VQ = NULL; 2643 MatrixCoefficient *MQ = NULL; 2644 2645 // PA extension 2646 const DofToQuad *maps; ///< Not owned 2647 const GeometricFactors *geom; ///< Not owned 2648 int dim, sdim, ne, dofs1D, quad1D; 2649 Vector pa_data; 2650 2651 private: 2652 DenseMatrix dshape, dshapedxt, pelmat; 2653 int vdim = -1; 2654 DenseMatrix mcoeff; 2655 Vector vcoeff; 2656 2657 public: VectorDiffusionIntegrator()2658 VectorDiffusionIntegrator() { } 2659 2660 /** \brief Integrator with unit coefficient for caller-specified vector 2661 dimension. 2662 2663 If the vector dimension does not match the true dimension of the space, 2664 the resulting element matrix will be mathematically invalid. */ VectorDiffusionIntegrator(int vector_dimension)2665 VectorDiffusionIntegrator(int vector_dimension) 2666 : vdim(vector_dimension) { } 2667 VectorDiffusionIntegrator(Coefficient & q)2668 VectorDiffusionIntegrator(Coefficient &q) 2669 : Q(&q) { } 2670 2671 /** \brief Integrator with scalar coefficient for caller-specified vector 2672 dimension. 2673 2674 The element matrix is block-diagonal with \c vdim copies of the element 2675 matrix integrated with the \c Coefficient. 2676 2677 If the vector dimension does not match the true dimension of the space, 2678 the resulting element matrix will be mathematically invalid. */ VectorDiffusionIntegrator(Coefficient & q,int vector_dimension)2679 VectorDiffusionIntegrator(Coefficient &q, int vector_dimension) 2680 : Q(&q), vdim(vector_dimension) { } 2681 2682 /** \brief Integrator with \c VectorCoefficient. The vector dimension of the 2683 \c FiniteElementSpace is assumed to be the same as the dimension of the 2684 \c Vector. 2685 2686 The element matrix is block-diagonal and each block is integrated with 2687 coefficient q_i. 2688 2689 If the vector dimension does not match the true dimension of the space, 2690 the resulting element matrix will be mathematically invalid. */ VectorDiffusionIntegrator(VectorCoefficient & vq)2691 VectorDiffusionIntegrator(VectorCoefficient &vq) 2692 : VQ(&vq), vdim(vq.GetVDim()) { } 2693 2694 /** \brief Integrator with \c MatrixCoefficient. The vector dimension of the 2695 \c FiniteElementSpace is assumed to be the same as the dimension of the 2696 \c Matrix. 2697 2698 The element matrix is populated in each block. Each block is integrated 2699 with coefficient q_ij. 2700 2701 If the vector dimension does not match the true dimension of the space, 2702 the resulting element matrix will be mathematically invalid. */ VectorDiffusionIntegrator(MatrixCoefficient & mq)2703 VectorDiffusionIntegrator(MatrixCoefficient& mq) 2704 : MQ(&mq), vdim(mq.GetVDim()) { } 2705 2706 virtual void AssembleElementMatrix(const FiniteElement &el, 2707 ElementTransformation &Trans, 2708 DenseMatrix &elmat); 2709 virtual void AssembleElementVector(const FiniteElement &el, 2710 ElementTransformation &Tr, 2711 const Vector &elfun, Vector &elvect); 2712 using BilinearFormIntegrator::AssemblePA; 2713 virtual void AssemblePA(const FiniteElementSpace &fes); 2714 virtual void AssembleMF(const FiniteElementSpace &fes); 2715 virtual void AssembleDiagonalPA(Vector &diag); 2716 virtual void AssembleDiagonalMF(Vector &diag); 2717 virtual void AddMultPA(const Vector &x, Vector &y) const; 2718 virtual void AddMultMF(const Vector &x, Vector &y) const; SupportsCeed() const2719 bool SupportsCeed() const { return DeviceCanUseCeed(); } 2720 }; 2721 2722 /** Integrator for the linear elasticity form: 2723 a(u,v) = (lambda div(u), div(v)) + (2 mu e(u), e(v)), 2724 where e(v) = (1/2) (grad(v) + grad(v)^T). 2725 This is a 'Vector' integrator, i.e. defined for FE spaces 2726 using multiple copies of a scalar FE space. */ 2727 class ElasticityIntegrator : public BilinearFormIntegrator 2728 { 2729 protected: 2730 double q_lambda, q_mu; 2731 Coefficient *lambda, *mu; 2732 2733 private: 2734 #ifndef MFEM_THREAD_SAFE 2735 Vector shape; 2736 DenseMatrix dshape, gshape, pelmat; 2737 Vector divshape; 2738 #endif 2739 2740 public: ElasticityIntegrator(Coefficient & l,Coefficient & m)2741 ElasticityIntegrator(Coefficient &l, Coefficient &m) 2742 { lambda = &l; mu = &m; } 2743 /** With this constructor lambda = q_l * m and mu = q_m * m; 2744 if dim * q_l + 2 * q_m = 0 then trace(sigma) = 0. */ ElasticityIntegrator(Coefficient & m,double q_l,double q_m)2745 ElasticityIntegrator(Coefficient &m, double q_l, double q_m) 2746 { lambda = NULL; mu = &m; q_lambda = q_l; q_mu = q_m; } 2747 2748 virtual void AssembleElementMatrix(const FiniteElement &, 2749 ElementTransformation &, 2750 DenseMatrix &); 2751 2752 /** Compute the stress corresponding to the local displacement @a u and 2753 interpolate it at the nodes of the given @a fluxelem. Only the symmetric 2754 part of the stress is stored, so that the size of @a flux is equal to 2755 the number of DOFs in @a fluxelem times dim*(dim+1)/2. In 2D, the order 2756 of the stress components is: s_xx, s_yy, s_xy. In 3D, it is: s_xx, s_yy, 2757 s_zz, s_xy, s_xz, s_yz. In other words, @a flux is the local vector for 2758 a FE space with dim*(dim+1)/2 vector components, based on the finite 2759 element @a fluxelem. */ 2760 virtual void ComputeElementFlux(const FiniteElement &el, 2761 ElementTransformation &Trans, 2762 Vector &u, 2763 const FiniteElement &fluxelem, 2764 Vector &flux, bool with_coef = true); 2765 2766 /** Compute the element energy (integral of the strain energy density) 2767 corresponding to the stress represented by @a flux which is a vector of 2768 coefficients multiplying the basis functions defined by @a fluxelem. In 2769 other words, @a flux is the local vector for a FE space with 2770 dim*(dim+1)/2 vector components, based on the finite element @a fluxelem. 2771 The number of components, dim*(dim+1)/2 is such that it represents the 2772 symmetric part of the (symmetric) stress tensor. The order of the 2773 components is: s_xx, s_yy, s_xy in 2D, and s_xx, s_yy, s_zz, s_xy, s_xz, 2774 s_yz in 3D. */ 2775 virtual double ComputeFluxEnergy(const FiniteElement &fluxelem, 2776 ElementTransformation &Trans, 2777 Vector &flux, Vector *d_energy = NULL); 2778 }; 2779 2780 /** Integrator for the DG form: 2781 alpha < rho_u (u.n) {v},[w] > + beta < rho_u |u.n| [v],[w] >, 2782 where v and w are the trial and test variables, respectively, and rho/u are 2783 given scalar/vector coefficients. {v} represents the average value of v on 2784 the face and [v] is the jump such that {v}=(v1+v2)/2 and [v]=(v1-v2) for the 2785 face between elements 1 and 2. For boundary elements, v2=0. The vector 2786 coefficient, u, is assumed to be continuous across the faces and when given 2787 the scalar coefficient, rho, is assumed to be discontinuous. The integrator 2788 uses the upwind value of rho, rho_u, which is value from the side into which 2789 the vector coefficient, u, points. 2790 2791 One use case for this integrator is to discretize the operator -u.grad(v) 2792 with a DG formulation. The resulting formulation uses the 2793 ConvectionIntegrator (with coefficient u, and parameter alpha = -1) and the 2794 transpose of the DGTraceIntegrator (with coefficient u, and parameters alpha 2795 = 1, beta = -1/2 to use the upwind face flux, see also 2796 NonconservativeDGTraceIntegrator). This discretization and the handling of 2797 the inflow and outflow boundaries is illustrated in Example 9/9p. 2798 2799 Another use case for this integrator is to discretize the operator -div(u v) 2800 with a DG formulation. The resulting formulation is conservative and 2801 consists of the ConservativeConvectionIntegrator (with coefficient u, and 2802 parameter alpha = -1) plus the DGTraceIntegrator (with coefficient u, and 2803 parameters alpha = -1, beta = -1/2 to use the upwind face flux). 2804 */ 2805 class DGTraceIntegrator : public BilinearFormIntegrator 2806 { 2807 protected: 2808 Coefficient *rho; 2809 VectorCoefficient *u; 2810 double alpha, beta; 2811 // PA extension 2812 Vector pa_data; 2813 const DofToQuad *maps; ///< Not owned 2814 const FaceGeometricFactors *geom; ///< Not owned 2815 int dim, nf, nq, dofs1D, quad1D; 2816 2817 private: 2818 Vector shape1, shape2; 2819 2820 public: 2821 /// Construct integrator with rho = 1, b = 0.5*a. DGTraceIntegrator(VectorCoefficient & u_,double a)2822 DGTraceIntegrator(VectorCoefficient &u_, double a) 2823 { rho = NULL; u = &u_; alpha = a; beta = 0.5*a; } 2824 2825 /// Construct integrator with rho = 1. DGTraceIntegrator(VectorCoefficient & u_,double a,double b)2826 DGTraceIntegrator(VectorCoefficient &u_, double a, double b) 2827 { rho = NULL; u = &u_; alpha = a; beta = b; } 2828 DGTraceIntegrator(Coefficient & rho_,VectorCoefficient & u_,double a,double b)2829 DGTraceIntegrator(Coefficient &rho_, VectorCoefficient &u_, 2830 double a, double b) 2831 { rho = &rho_; u = &u_; alpha = a; beta = b; } 2832 2833 using BilinearFormIntegrator::AssembleFaceMatrix; 2834 virtual void AssembleFaceMatrix(const FiniteElement &el1, 2835 const FiniteElement &el2, 2836 FaceElementTransformations &Trans, 2837 DenseMatrix &elmat); 2838 2839 using BilinearFormIntegrator::AssemblePA; 2840 2841 virtual void AssemblePAInteriorFaces(const FiniteElementSpace &fes); 2842 2843 virtual void AssemblePABoundaryFaces(const FiniteElementSpace &fes); 2844 2845 virtual void AddMultTransposePA(const Vector &x, Vector &y) const; 2846 2847 virtual void AddMultPA(const Vector&, Vector&) const; 2848 2849 virtual void AssembleEAInteriorFaces(const FiniteElementSpace& fes, 2850 Vector &ea_data_int, 2851 Vector &ea_data_ext, 2852 const bool add); 2853 2854 virtual void AssembleEABoundaryFaces(const FiniteElementSpace& fes, 2855 Vector &ea_data_bdr, 2856 const bool add); 2857 2858 static const IntegrationRule &GetRule(Geometry::Type geom, int order, 2859 FaceElementTransformations &T); 2860 2861 private: 2862 void SetupPA(const FiniteElementSpace &fes, FaceType type); 2863 }; 2864 2865 // Alias for @a DGTraceIntegrator. 2866 using ConservativeDGTraceIntegrator = DGTraceIntegrator; 2867 2868 /** Integrator that represents the face terms used for the non-conservative 2869 DG discretization of the convection equation: 2870 -alpha < rho_u (u.n) {v},[w] > + beta < rho_u |u.n| [v],[w] >. 2871 2872 This integrator can be used with together with ConvectionIntegrator to 2873 implement an upwind DG discretization in non-conservative form, see ex9 and 2874 ex9p. */ 2875 class NonconservativeDGTraceIntegrator : public TransposeIntegrator 2876 { 2877 public: NonconservativeDGTraceIntegrator(VectorCoefficient & u,double a)2878 NonconservativeDGTraceIntegrator(VectorCoefficient &u, double a) 2879 : TransposeIntegrator(new DGTraceIntegrator(u, -a, 0.5*a)) { } 2880 NonconservativeDGTraceIntegrator(VectorCoefficient & u,double a,double b)2881 NonconservativeDGTraceIntegrator(VectorCoefficient &u, double a, double b) 2882 : TransposeIntegrator(new DGTraceIntegrator(u, -a, b)) { } 2883 NonconservativeDGTraceIntegrator(Coefficient & rho,VectorCoefficient & u,double a,double b)2884 NonconservativeDGTraceIntegrator(Coefficient &rho, VectorCoefficient &u, 2885 double a, double b) 2886 : TransposeIntegrator(new DGTraceIntegrator(rho, u, -a, b)) { } 2887 }; 2888 2889 /** Integrator for the DG form: 2890 2891 - < {(Q grad(u)).n}, [v] > + sigma < [u], {(Q grad(v)).n} > 2892 + kappa < {h^{-1} Q} [u], [v] >, 2893 2894 where Q is a scalar or matrix diffusion coefficient and u, v are the trial 2895 and test spaces, respectively. The parameters sigma and kappa determine the 2896 DG method to be used (when this integrator is added to the "broken" 2897 DiffusionIntegrator): 2898 * sigma = -1, kappa >= kappa0: symm. interior penalty (IP or SIPG) method, 2899 * sigma = +1, kappa > 0: non-symmetric interior penalty (NIPG) method, 2900 * sigma = +1, kappa = 0: the method of Baumann and Oden. */ 2901 class DGDiffusionIntegrator : public BilinearFormIntegrator 2902 { 2903 protected: 2904 Coefficient *Q; 2905 MatrixCoefficient *MQ; 2906 double sigma, kappa; 2907 2908 // these are not thread-safe! 2909 Vector shape1, shape2, dshape1dn, dshape2dn, nor, nh, ni; 2910 DenseMatrix jmat, dshape1, dshape2, mq, adjJ; 2911 2912 public: DGDiffusionIntegrator(const double s,const double k)2913 DGDiffusionIntegrator(const double s, const double k) 2914 : Q(NULL), MQ(NULL), sigma(s), kappa(k) { } DGDiffusionIntegrator(Coefficient & q,const double s,const double k)2915 DGDiffusionIntegrator(Coefficient &q, const double s, const double k) 2916 : Q(&q), MQ(NULL), sigma(s), kappa(k) { } DGDiffusionIntegrator(MatrixCoefficient & q,const double s,const double k)2917 DGDiffusionIntegrator(MatrixCoefficient &q, const double s, const double k) 2918 : Q(NULL), MQ(&q), sigma(s), kappa(k) { } 2919 using BilinearFormIntegrator::AssembleFaceMatrix; 2920 virtual void AssembleFaceMatrix(const FiniteElement &el1, 2921 const FiniteElement &el2, 2922 FaceElementTransformations &Trans, 2923 DenseMatrix &elmat); 2924 }; 2925 2926 /** Integrator for the "BR2" diffusion stabilization term 2927 2928 sum_e eta (r_e([u]), r_e([v])) 2929 2930 where r_e is the lifting operator defined on each edge e. The parameter eta 2931 can be chosen to be one to obtain a stable discretization. The constructor 2932 for this integrator requires the finite element space because the lifting 2933 operator depends on the element-wise inverse mass matrix. 2934 2935 BR2 stands for the second method of Bassi and Rebay: 2936 2937 - F. Bassi and S. Rebay. A high order discontinuous Galerkin method for 2938 compressible turbulent flows. In B. Cockburn, G. E. Karniadakis, and 2939 C.-W. Shu, editors, Discontinuous Galerkin Methods, pages 77-88. Springer 2940 Berlin Heidelberg, 2000. 2941 - D. N. Arnold, F. Brezzi, B. Cockburn, and L. D. Marini. Unified analysis 2942 of discontinuous Galerkin methods for elliptic problems. SIAM Journal on 2943 Numerical Analysis, 39(5):1749-1779, 2002. 2944 */ 2945 class DGDiffusionBR2Integrator : public BilinearFormIntegrator 2946 { 2947 protected: 2948 double eta; 2949 2950 // Block factorizations of local mass matrices, with offsets for the case of 2951 // not equally sized blocks (mixed meshes, p-refinement) 2952 Array<double> Minv; 2953 Array<int> ipiv; 2954 Array<int> ipiv_offsets, Minv_offsets; 2955 2956 Vector shape1, shape2; 2957 2958 DenseMatrix R11, R12, R21, R22; 2959 DenseMatrix MinvR11, MinvR12, MinvR21, MinvR22; 2960 DenseMatrix Re, MinvRe; 2961 2962 public: 2963 DGDiffusionBR2Integrator(class FiniteElementSpace *fes, double e = 1.0); 2964 using BilinearFormIntegrator::AssembleFaceMatrix; 2965 virtual void AssembleFaceMatrix(const FiniteElement &el1, 2966 const FiniteElement &el2, 2967 FaceElementTransformations &Trans, 2968 DenseMatrix &elmat); 2969 }; 2970 2971 /** Integrator for the DG elasticity form, for the formulations see: 2972 - PhD Thesis of Jonas De Basabe, High-Order Finite %Element Methods for 2973 Seismic Wave Propagation, UT Austin, 2009, p. 23, and references therein 2974 - Peter Hansbo and Mats G. Larson, Discontinuous Galerkin and the 2975 Crouzeix-Raviart %Element: Application to Elasticity, PREPRINT 2000-09, 2976 p.3 2977 2978 \f[ 2979 - \left< \{ \tau(u) \}, [v] \right> + \alpha \left< \{ \tau(v) \}, [u] 2980 \right> + \kappa \left< h^{-1} \{ \lambda + 2 \mu \} [u], [v] \right> 2981 \f] 2982 2983 where \f$ \left<u, v\right> = \int_{F} u \cdot v \f$, and \f$ F \f$ is a 2984 face which is either a boundary face \f$ F_b \f$ of an element \f$ K \f$ or 2985 an interior face \f$ F_i \f$ separating elements \f$ K_1 \f$ and \f$ K_2 \f$. 2986 2987 In the bilinear form above \f$ \tau(u) \f$ is traction, and it's also 2988 \f$ \tau(u) = \sigma(u) \cdot \vec{n} \f$, where \f$ \sigma(u) \f$ is 2989 stress, and \f$ \vec{n} \f$ is the unit normal vector w.r.t. to \f$ F \f$. 2990 2991 In other words, we have 2992 \f[ 2993 - \left< \{ \sigma(u) \cdot \vec{n} \}, [v] \right> + \alpha \left< \{ 2994 \sigma(v) \cdot \vec{n} \}, [u] \right> + \kappa \left< h^{-1} \{ 2995 \lambda + 2 \mu \} [u], [v] \right> 2996 \f] 2997 2998 For isotropic media 2999 \f[ 3000 \begin{split} 3001 \sigma(u) &= \lambda \nabla \cdot u I + 2 \mu \varepsilon(u) \\ 3002 &= \lambda \nabla \cdot u I + 2 \mu \frac{1}{2} (\nabla u + \nabla 3003 u^T) \\ 3004 &= \lambda \nabla \cdot u I + \mu (\nabla u + \nabla u^T) 3005 \end{split} 3006 \f] 3007 3008 where \f$ I \f$ is identity matrix, \f$ \lambda \f$ and \f$ \mu \f$ are Lame 3009 coefficients (see ElasticityIntegrator), \f$ u, v \f$ are the trial and test 3010 functions, respectively. 3011 3012 The parameters \f$ \alpha \f$ and \f$ \kappa \f$ determine the DG method to 3013 use (when this integrator is added to the "broken" ElasticityIntegrator): 3014 3015 - IIPG, \f$\alpha = 0\f$, 3016 C. Dawson, S. Sun, M. Wheeler, Compatible algorithms for coupled flow and 3017 transport, Comp. Meth. Appl. Mech. Eng., 193(23-26), 2565-2580, 2004. 3018 3019 - SIPG, \f$\alpha = -1\f$, 3020 M. Grote, A. Schneebeli, D. Schotzau, Discontinuous Galerkin Finite 3021 %Element Method for the Wave Equation, SINUM, 44(6), 2408-2431, 2006. 3022 3023 - NIPG, \f$\alpha = 1\f$, 3024 B. Riviere, M. Wheeler, V. Girault, A Priori Error Estimates for Finite 3025 %Element Methods Based on Discontinuous Approximation Spaces for Elliptic 3026 Problems, SINUM, 39(3), 902-931, 2001. 3027 3028 This is a '%Vector' integrator, i.e. defined for FE spaces using multiple 3029 copies of a scalar FE space. 3030 */ 3031 class DGElasticityIntegrator : public BilinearFormIntegrator 3032 { 3033 public: DGElasticityIntegrator(double alpha_,double kappa_)3034 DGElasticityIntegrator(double alpha_, double kappa_) 3035 : lambda(NULL), mu(NULL), alpha(alpha_), kappa(kappa_) { } 3036 DGElasticityIntegrator(Coefficient & lambda_,Coefficient & mu_,double alpha_,double kappa_)3037 DGElasticityIntegrator(Coefficient &lambda_, Coefficient &mu_, 3038 double alpha_, double kappa_) 3039 : lambda(&lambda_), mu(&mu_), alpha(alpha_), kappa(kappa_) { } 3040 3041 using BilinearFormIntegrator::AssembleFaceMatrix; 3042 virtual void AssembleFaceMatrix(const FiniteElement &el1, 3043 const FiniteElement &el2, 3044 FaceElementTransformations &Trans, 3045 DenseMatrix &elmat); 3046 3047 protected: 3048 Coefficient *lambda, *mu; 3049 double alpha, kappa; 3050 3051 #ifndef MFEM_THREAD_SAFE 3052 // values of all scalar basis functions for one component of u (which is a 3053 // vector) at the integration point in the reference space 3054 Vector shape1, shape2; 3055 // values of derivatives of all scalar basis functions for one component 3056 // of u (which is a vector) at the integration point in the reference space 3057 DenseMatrix dshape1, dshape2; 3058 // Adjugate of the Jacobian of the transformation: adjJ = det(J) J^{-1} 3059 DenseMatrix adjJ; 3060 // gradient of shape functions in the real (physical, not reference) 3061 // coordinates, scaled by det(J): 3062 // dshape_ps(jdof,jm) = sum_{t} adjJ(t,jm)*dshape(jdof,t) 3063 DenseMatrix dshape1_ps, dshape2_ps; 3064 Vector nor; // nor = |weight(J_face)| n 3065 Vector nL1, nL2; // nL1 = (lambda1 * ip.weight / detJ1) nor 3066 Vector nM1, nM2; // nM1 = (mu1 * ip.weight / detJ1) nor 3067 Vector dshape1_dnM, dshape2_dnM; // dshape1_dnM = dshape1_ps . nM1 3068 // 'jmat' corresponds to the term: kappa <h^{-1} {lambda + 2 mu} [u], [v]> 3069 DenseMatrix jmat; 3070 #endif 3071 3072 static void AssembleBlock( 3073 const int dim, const int row_ndofs, const int col_ndofs, 3074 const int row_offset, const int col_offset, 3075 const double jmatcoef, const Vector &col_nL, const Vector &col_nM, 3076 const Vector &row_shape, const Vector &col_shape, 3077 const Vector &col_dshape_dnM, const DenseMatrix &col_dshape, 3078 DenseMatrix &elmat, DenseMatrix &jmat); 3079 }; 3080 3081 /** Integrator for the DPG form: < v, [w] > over all faces (the interface) where 3082 the trial variable v is defined on the interface and the test variable w is 3083 defined inside the elements, generally in a DG space. */ 3084 class TraceJumpIntegrator : public BilinearFormIntegrator 3085 { 3086 private: 3087 Vector face_shape, shape1, shape2; 3088 3089 public: TraceJumpIntegrator()3090 TraceJumpIntegrator() { } 3091 using BilinearFormIntegrator::AssembleFaceMatrix; 3092 virtual void AssembleFaceMatrix(const FiniteElement &trial_face_fe, 3093 const FiniteElement &test_fe1, 3094 const FiniteElement &test_fe2, 3095 FaceElementTransformations &Trans, 3096 DenseMatrix &elmat); 3097 }; 3098 3099 /** Integrator for the form: < v, [w.n] > over all faces (the interface) where 3100 the trial variable v is defined on the interface and the test variable w is 3101 in an H(div)-conforming space. */ 3102 class NormalTraceJumpIntegrator : public BilinearFormIntegrator 3103 { 3104 private: 3105 Vector face_shape, normal, shape1_n, shape2_n; 3106 DenseMatrix shape1, shape2; 3107 3108 public: NormalTraceJumpIntegrator()3109 NormalTraceJumpIntegrator() { } 3110 using BilinearFormIntegrator::AssembleFaceMatrix; 3111 virtual void AssembleFaceMatrix(const FiniteElement &trial_face_fe, 3112 const FiniteElement &test_fe1, 3113 const FiniteElement &test_fe2, 3114 FaceElementTransformations &Trans, 3115 DenseMatrix &elmat); 3116 }; 3117 3118 /** Abstract class to serve as a base for local interpolators to be used in the 3119 DiscreteLinearOperator class. */ 3120 class DiscreteInterpolator : public BilinearFormIntegrator { }; 3121 3122 3123 /** Class for constructing the gradient as a DiscreteLinearOperator from an 3124 H1-conforming space to an H(curl)-conforming space. The range space can be 3125 vector L2 space as well. */ 3126 class GradientInterpolator : public DiscreteInterpolator 3127 { 3128 public: GradientInterpolator()3129 GradientInterpolator() : dofquad_fe(NULL) { } ~GradientInterpolator()3130 virtual ~GradientInterpolator() { delete dofquad_fe; } 3131 AssembleElementMatrix2(const FiniteElement & h1_fe,const FiniteElement & nd_fe,ElementTransformation & Trans,DenseMatrix & elmat)3132 virtual void AssembleElementMatrix2(const FiniteElement &h1_fe, 3133 const FiniteElement &nd_fe, 3134 ElementTransformation &Trans, 3135 DenseMatrix &elmat) 3136 { nd_fe.ProjectGrad(h1_fe, Trans, elmat); } 3137 3138 using BilinearFormIntegrator::AssemblePA; 3139 3140 /** @brief Setup method for PA data. 3141 3142 @param[in] trial_fes H1 Lagrange space 3143 @param[in] test_fes H(curl) Nedelec space 3144 */ 3145 virtual void AssemblePA(const FiniteElementSpace &trial_fes, 3146 const FiniteElementSpace &test_fes); 3147 3148 virtual void AddMultPA(const Vector &x, Vector &y) const; 3149 virtual void AddMultTransposePA(const Vector &x, Vector &y) const; 3150 3151 private: 3152 /// 1D finite element that generates and owns the 1D DofToQuad maps below 3153 FiniteElement * dofquad_fe; 3154 3155 bool B_id; // is the B basis operator (maps_C_C) the identity? 3156 const DofToQuad *maps_C_C; // one-d map with Lobatto rows, Lobatto columns 3157 const DofToQuad *maps_O_C; // one-d map with Legendre rows, Lobatto columns 3158 int dim, ne, o_dofs1D, c_dofs1D; 3159 }; 3160 3161 3162 /** Class for constructing the identity map as a DiscreteLinearOperator. This 3163 is the discrete embedding matrix when the domain space is a subspace of 3164 the range space. Otherwise, a dof projection matrix is constructed. */ 3165 class IdentityInterpolator : public DiscreteInterpolator 3166 { 3167 public: AssembleElementMatrix2(const FiniteElement & dom_fe,const FiniteElement & ran_fe,ElementTransformation & Trans,DenseMatrix & elmat)3168 virtual void AssembleElementMatrix2(const FiniteElement &dom_fe, 3169 const FiniteElement &ran_fe, 3170 ElementTransformation &Trans, 3171 DenseMatrix &elmat) 3172 { ran_fe.Project(dom_fe, Trans, elmat); } 3173 3174 using BilinearFormIntegrator::AssemblePA; 3175 3176 virtual void AssemblePA(const FiniteElementSpace &trial_fes, 3177 const FiniteElementSpace &test_fes); 3178 3179 virtual void AddMultPA(const Vector &x, Vector &y) const; 3180 virtual void AddMultTransposePA(const Vector &x, Vector &y) const; 3181 3182 private: 3183 /// 1D finite element that generates and owns the 1D DofToQuad maps below 3184 FiniteElement * dofquad_fe; 3185 3186 const DofToQuad *maps_C_C; // one-d map with Lobatto rows, Lobatto columns 3187 const DofToQuad *maps_O_C; // one-d map with Legendre rows, Lobatto columns 3188 int dim, ne, o_dofs1D, c_dofs1D; 3189 3190 Vector pa_data; 3191 }; 3192 3193 3194 /** Class for constructing the (local) discrete curl matrix which can be used 3195 as an integrator in a DiscreteLinearOperator object to assemble the global 3196 discrete curl matrix. */ 3197 class CurlInterpolator : public DiscreteInterpolator 3198 { 3199 public: AssembleElementMatrix2(const FiniteElement & dom_fe,const FiniteElement & ran_fe,ElementTransformation & Trans,DenseMatrix & elmat)3200 virtual void AssembleElementMatrix2(const FiniteElement &dom_fe, 3201 const FiniteElement &ran_fe, 3202 ElementTransformation &Trans, 3203 DenseMatrix &elmat) 3204 { ran_fe.ProjectCurl(dom_fe, Trans, elmat); } 3205 }; 3206 3207 3208 /** Class for constructing the (local) discrete divergence matrix which can 3209 be used as an integrator in a DiscreteLinearOperator object to assemble 3210 the global discrete divergence matrix. 3211 3212 Note: Since the dofs in the L2_FECollection are nodal values, the local 3213 discrete divergence matrix (with an RT-type domain space) will depend on 3214 the transformation. On the other hand, the local matrix returned by 3215 VectorFEDivergenceIntegrator is independent of the transformation. */ 3216 class DivergenceInterpolator : public DiscreteInterpolator 3217 { 3218 public: AssembleElementMatrix2(const FiniteElement & dom_fe,const FiniteElement & ran_fe,ElementTransformation & Trans,DenseMatrix & elmat)3219 virtual void AssembleElementMatrix2(const FiniteElement &dom_fe, 3220 const FiniteElement &ran_fe, 3221 ElementTransformation &Trans, 3222 DenseMatrix &elmat) 3223 { ran_fe.ProjectDiv(dom_fe, Trans, elmat); } 3224 }; 3225 3226 3227 /** A trace face interpolator class for interpolating the normal component of 3228 the domain space, e.g. vector H1, into the range space, e.g. the trace of 3229 RT which uses FiniteElement::INTEGRAL map type. */ 3230 class NormalInterpolator : public DiscreteInterpolator 3231 { 3232 public: 3233 virtual void AssembleElementMatrix2(const FiniteElement &dom_fe, 3234 const FiniteElement &ran_fe, 3235 ElementTransformation &Trans, 3236 DenseMatrix &elmat); 3237 }; 3238 3239 /** Interpolator of a scalar coefficient multiplied by a scalar field onto 3240 another scalar field. Note that this can produce inaccurate fields unless 3241 the target is sufficiently high order. */ 3242 class ScalarProductInterpolator : public DiscreteInterpolator 3243 { 3244 public: ScalarProductInterpolator(Coefficient & sc)3245 ScalarProductInterpolator(Coefficient & sc) : Q(&sc) { } 3246 3247 virtual void AssembleElementMatrix2(const FiniteElement &dom_fe, 3248 const FiniteElement &ran_fe, 3249 ElementTransformation &Trans, 3250 DenseMatrix &elmat); 3251 3252 protected: 3253 Coefficient *Q; 3254 }; 3255 3256 /** Interpolator of a scalar coefficient multiplied by a vector field onto 3257 another vector field. Note that this can produce inaccurate fields unless 3258 the target is sufficiently high order. */ 3259 class ScalarVectorProductInterpolator : public DiscreteInterpolator 3260 { 3261 public: ScalarVectorProductInterpolator(Coefficient & sc)3262 ScalarVectorProductInterpolator(Coefficient & sc) 3263 : Q(&sc) { } 3264 3265 virtual void AssembleElementMatrix2(const FiniteElement &dom_fe, 3266 const FiniteElement &ran_fe, 3267 ElementTransformation &Trans, 3268 DenseMatrix &elmat); 3269 protected: 3270 Coefficient *Q; 3271 }; 3272 3273 /** Interpolator of a vector coefficient multiplied by a scalar field onto 3274 another vector field. Note that this can produce inaccurate fields unless 3275 the target is sufficiently high order. */ 3276 class VectorScalarProductInterpolator : public DiscreteInterpolator 3277 { 3278 public: VectorScalarProductInterpolator(VectorCoefficient & vc)3279 VectorScalarProductInterpolator(VectorCoefficient & vc) 3280 : VQ(&vc) { } 3281 3282 virtual void AssembleElementMatrix2(const FiniteElement &dom_fe, 3283 const FiniteElement &ran_fe, 3284 ElementTransformation &Trans, 3285 DenseMatrix &elmat); 3286 protected: 3287 VectorCoefficient *VQ; 3288 }; 3289 3290 /** Interpolator of the 2D cross product between a vector coefficient and an 3291 H(curl)-conforming field onto an L2-conforming field. */ 3292 class ScalarCrossProductInterpolator : public DiscreteInterpolator 3293 { 3294 public: ScalarCrossProductInterpolator(VectorCoefficient & vc)3295 ScalarCrossProductInterpolator(VectorCoefficient & vc) 3296 : VQ(&vc) { } 3297 3298 virtual void AssembleElementMatrix2(const FiniteElement &nd_fe, 3299 const FiniteElement &l2_fe, 3300 ElementTransformation &Trans, 3301 DenseMatrix &elmat); 3302 protected: 3303 VectorCoefficient *VQ; 3304 }; 3305 3306 /** Interpolator of the cross product between a vector coefficient and an 3307 H(curl)-conforming field onto an H(div)-conforming field. The range space 3308 can also be vector L2. */ 3309 class VectorCrossProductInterpolator : public DiscreteInterpolator 3310 { 3311 public: VectorCrossProductInterpolator(VectorCoefficient & vc)3312 VectorCrossProductInterpolator(VectorCoefficient & vc) 3313 : VQ(&vc) { } 3314 3315 virtual void AssembleElementMatrix2(const FiniteElement &nd_fe, 3316 const FiniteElement &rt_fe, 3317 ElementTransformation &Trans, 3318 DenseMatrix &elmat); 3319 protected: 3320 VectorCoefficient *VQ; 3321 }; 3322 3323 /** Interpolator of the inner product between a vector coefficient and an 3324 H(div)-conforming field onto an L2-conforming field. The range space can 3325 also be H1. */ 3326 class VectorInnerProductInterpolator : public DiscreteInterpolator 3327 { 3328 public: VectorInnerProductInterpolator(VectorCoefficient & vc)3329 VectorInnerProductInterpolator(VectorCoefficient & vc) : VQ(&vc) { } 3330 3331 virtual void AssembleElementMatrix2(const FiniteElement &rt_fe, 3332 const FiniteElement &l2_fe, 3333 ElementTransformation &Trans, 3334 DenseMatrix &elmat); 3335 protected: 3336 VectorCoefficient *VQ; 3337 }; 3338 3339 3340 3341 // PA Diffusion Assemble 2D kernel 3342 template<const int T_SDIM> 3343 void PADiffusionSetup2D(const int Q1D, 3344 const int coeffDim, 3345 const int NE, 3346 const Array<double> &w, 3347 const Vector &j, 3348 const Vector &c, 3349 Vector &d); 3350 3351 } 3352 #endif 3353