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