MFEM  v4.6.0
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 /// Abstract base class BilinearFormIntegrator
25 {
26 protected:
28  : NonlinearFormIntegrator(ir) { }
29 
30 public:
31  // TODO: add support for other assembly levels (in addition to PA) and their
32  // actions.
33 
34  // TODO: for mixed meshes the quadrature rules to be used by methods like
35  // AssemblePA() can be given as a QuadratureSpace, e.g. using a new method:
36  // SetQuadratureSpace().
37 
38  // TODO: the methods for the various assembly levels make sense even in the
39  // base class NonlinearFormIntegrator, except that not all assembly levels
40  // make sense for the action of the nonlinear operator (but they all make
41  // sense for its Jacobian).
42 
44 
45  /// Method defining partial assembly.
46  /** The result of the partial assembly is stored internally so that it can be
47  used later in the methods AddMultPA() and AddMultTransposePA(). */
48  virtual void AssemblePA(const FiniteElementSpace &fes);
49  /** Used with BilinearFormIntegrators that have different spaces. */
50  virtual void AssemblePA(const FiniteElementSpace &trial_fes,
51  const FiniteElementSpace &test_fes);
52 
53  /// Method defining partial assembly on NURBS patches.
54  /** The result of the partial assembly is stored internally so that it can be
55  used later in the method AddMultNURBSPA(). */
56  virtual void AssembleNURBSPA(const FiniteElementSpace &fes);
57 
58  virtual void AssemblePABoundary(const FiniteElementSpace &fes);
59 
60  virtual void AssemblePAInteriorFaces(const FiniteElementSpace &fes);
61 
62  virtual void AssemblePABoundaryFaces(const FiniteElementSpace &fes);
63 
64  /// Assemble diagonal and add it to Vector @a diag.
65  virtual void AssembleDiagonalPA(Vector &diag);
66 
67  /// Assemble diagonal of ADA^T (A is this integrator) and add it to @a diag.
68  virtual void AssembleDiagonalPA_ADAt(const Vector &D, Vector &diag);
69 
70  /// Method for partially assembled action.
71  /** Perform the action of integrator on the input @a x and add the result to
72  the output @a y. Both @a x and @a y are E-vectors, i.e. they represent
73  the element-wise discontinuous version of the FE space.
74 
75  This method can be called only after the method AssemblePA() has been
76  called. */
77  virtual void AddMultPA(const Vector &x, Vector &y) const;
78 
79  /// Method for partially assembled action on NURBS patches.
80  virtual void AddMultNURBSPA(const Vector&x, Vector&y) const;
81 
82  /// Method for partially assembled transposed action.
83  /** Perform the transpose action of integrator on the input @a x and add the
84  result to the output @a y. Both @a x and @a y are E-vectors, i.e. they
85  represent the element-wise discontinuous version of the FE space.
86 
87  This method can be called only after the method AssemblePA() has been
88  called. */
89  virtual void AddMultTransposePA(const Vector &x, Vector &y) const;
90 
91  /// Method defining element assembly.
92  /** The result of the element assembly is added to the @a emat Vector if
93  @a add is true. Otherwise, if @a add is false, we set @a emat. */
94  virtual void AssembleEA(const FiniteElementSpace &fes, Vector &emat,
95  const bool add = true);
96  /** Used with BilinearFormIntegrators that have different spaces. */
97  // virtual void AssembleEA(const FiniteElementSpace &trial_fes,
98  // const FiniteElementSpace &test_fes,
99  // Vector &emat);
100 
101  /// Method defining matrix-free assembly.
102  /** The result of fully matrix-free assembly is stored internally so that it
103  can be used later in the methods AddMultMF() and AddMultTransposeMF(). */
104  virtual void AssembleMF(const FiniteElementSpace &fes);
105 
106  /** Perform the action of integrator on the input @a x and add the result to
107  the output @a y. Both @a x and @a y are E-vectors, i.e. they represent
108  the element-wise discontinuous version of the FE space.
109 
110  This method can be called only after the method AssembleMF() has been
111  called. */
112  virtual void AddMultMF(const Vector &x, Vector &y) const;
113 
114  /** Perform the transpose action of integrator on the input @a x and add the
115  result to the output @a y. Both @a x and @a y are E-vectors, i.e. they
116  represent the element-wise discontinuous version of the FE space.
117 
118  This method can be called only after the method AssemblePA() has been
119  called. */
120  virtual void AddMultTransposeMF(const Vector &x, Vector &y) const;
121 
122  /// Assemble diagonal and add it to Vector @a diag.
123  virtual void AssembleDiagonalMF(Vector &diag);
124 
125  virtual void AssembleEAInteriorFaces(const FiniteElementSpace &fes,
126  Vector &ea_data_int,
127  Vector &ea_data_ext,
128  const bool add = true);
129 
130  virtual void AssembleEABoundaryFaces(const FiniteElementSpace &fes,
131  Vector &ea_data_bdr,
132  const bool add = true);
133 
134  /// Given a particular Finite Element computes the element matrix elmat.
135  virtual void AssembleElementMatrix(const FiniteElement &el,
136  ElementTransformation &Trans,
137  DenseMatrix &elmat);
138 
139  /** Compute the local matrix representation of a bilinear form
140  a(u,v) defined on different trial (given by u) and test
141  (given by v) spaces. The rows in the local matrix correspond
142  to the test dofs and the columns -- to the trial dofs. */
143  virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
144  const FiniteElement &test_fe,
145  ElementTransformation &Trans,
146  DenseMatrix &elmat);
147 
148  /** Given a particular NURBS patch, computes the patch matrix as a
149  SparseMatrix @a smat.
150  */
151  virtual void AssemblePatchMatrix(const int patch,
152  const FiniteElementSpace &fes,
153  SparseMatrix*& smat);
154 
155  virtual void AssembleFaceMatrix(const FiniteElement &el1,
156  const FiniteElement &el2,
158  DenseMatrix &elmat);
159 
160  /** Abstract method used for assembling TraceFaceIntegrators in a
161  MixedBilinearForm. */
162  virtual void AssembleFaceMatrix(const FiniteElement &trial_face_fe,
163  const FiniteElement &test_fe1,
164  const FiniteElement &test_fe2,
166  DenseMatrix &elmat);
167 
168  /** Abstract method used for assembling TraceFaceIntegrators for
169  DPG weak formulations. */
170  virtual void AssembleTraceFaceMatrix(int elem,
171  const FiniteElement &trial_face_fe,
172  const FiniteElement &test_fe,
174  DenseMatrix &elmat);
175 
176 
177  /// @brief Perform the local action of the BilinearFormIntegrator.
178  /// Note that the default implementation in the base class is general but not
179  /// efficient.
180  virtual void AssembleElementVector(const FiniteElement &el,
182  const Vector &elfun, Vector &elvect);
183 
184  /// @brief Perform the local action of the BilinearFormIntegrator resulting
185  /// from a face integral term.
186  /// Note that the default implementation in the base class is general but not
187  /// efficient.
188  virtual void AssembleFaceVector(const FiniteElement &el1,
189  const FiniteElement &el2,
191  const Vector &elfun, Vector &elvect);
192 
193  virtual void AssembleElementGrad(const FiniteElement &el,
195  const Vector &elfun, DenseMatrix &elmat)
196  { AssembleElementMatrix(el, Tr, elmat); }
197 
198  virtual void AssembleFaceGrad(const FiniteElement &el1,
199  const FiniteElement &el2,
201  const Vector &elfun, DenseMatrix &elmat)
202  { AssembleFaceMatrix(el1, el2, Tr, elmat); }
203 
204  /** @brief Virtual method required for Zienkiewicz-Zhu type error estimators.
205 
206  The purpose of the method is to compute a local "flux" finite element
207  function given a local finite element solution. The "flux" function has
208  to be computed in terms of its coefficients (represented by the Vector
209  @a flux) which multiply the basis functions defined by the FiniteElement
210  @a fluxelem. Typically, the "flux" function will have more than one
211  component and consequently @a flux should be store the coefficients of
212  all components: first all coefficient for component 0, then all
213  coefficients for component 1, etc. What the "flux" function represents
214  depends on the specific integrator. For example, in the case of
215  DiffusionIntegrator, the flux is the gradient of the solution multiplied
216  by the diffusion coefficient.
217 
218  @param[in] el FiniteElement of the solution.
219  @param[in] Trans The ElementTransformation describing the physical
220  position of the mesh element.
221  @param[in] u Solution coefficients representing the expansion of the
222  solution function in the basis of @a el.
223  @param[in] fluxelem FiniteElement of the "flux".
224  @param[out] flux "Flux" coefficients representing the expansion of the
225  "flux" function in the basis of @a fluxelem. The size
226  of @a flux as a Vector has to be set by this method,
227  e.g. using Vector::SetSize().
228  @param[in] with_coef If zero (the default value is 1) the implementation
229  of the method may choose not to scale the "flux"
230  function by any coefficients describing the
231  integrator.
232  @param[in] ir If passed (the default value is NULL), the implementation
233  of the method will ignore the integration rule provided
234  by the @a fluxelem parameter and, instead, compute the
235  discrete flux at the points specified by the integration
236  rule @a ir.
237  */
238  virtual void ComputeElementFlux(const FiniteElement &el,
239  ElementTransformation &Trans,
240  Vector &u,
241  const FiniteElement &fluxelem,
242  Vector &flux, bool with_coef = true,
243  const IntegrationRule *ir = NULL) { }
244 
245  /** @brief Virtual method required for Zienkiewicz-Zhu type error estimators.
246 
247  The purpose of this method is to compute a local number that measures the
248  energy of a given "flux" function (see ComputeElementFlux() for a
249  description of the "flux" function). Typically, the energy of a "flux"
250  function should be equal to a_local(u,u), if the "flux" is defined from
251  a solution u; here a_local(.,.) denotes the element-local bilinear
252  form represented by the integrator.
253 
254  @param[in] fluxelem FiniteElement of the "flux".
255  @param[in] Trans The ElementTransformation describing the physical
256  position of the mesh element.
257  @param[in] flux "Flux" coefficients representing the expansion of the
258  "flux" function in the basis of @a fluxelem.
259  @param[out] d_energy If not NULL, the given Vector should be set to
260  represent directional energy split that can be used
261  for anisotropic error estimation.
262  @returns The computed energy.
263  */
264  virtual double ComputeFluxEnergy(const FiniteElement &fluxelem,
265  ElementTransformation &Trans,
266  Vector &flux, Vector *d_energy = NULL)
267  { return 0.0; }
268 
270 };
271 
272 /** Wraps a given @a BilinearFormIntegrator and transposes the resulting element
273  matrices. See for example ex9, ex9p. */
275 {
276 private:
277  int own_bfi;
279 
280  DenseMatrix bfi_elmat;
281 
282 public:
283  TransposeIntegrator (BilinearFormIntegrator *bfi_, int own_bfi_ = 1)
284  { bfi = bfi_; own_bfi = own_bfi_; }
285 
286  virtual void SetIntRule(const IntegrationRule *ir);
287 
288  virtual void AssembleElementMatrix(const FiniteElement &el,
289  ElementTransformation &Trans,
290  DenseMatrix &elmat);
291 
292  virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
293  const FiniteElement &test_fe,
294  ElementTransformation &Trans,
295  DenseMatrix &elmat);
296 
298  virtual void AssembleFaceMatrix(const FiniteElement &el1,
299  const FiniteElement &el2,
301  DenseMatrix &elmat);
302 
304 
305  virtual void AssemblePA(const FiniteElementSpace& fes)
306  {
307  bfi->AssemblePA(fes);
308  }
309 
310  virtual void AssemblePA(const FiniteElementSpace &trial_fes,
311  const FiniteElementSpace &test_fes)
312  {
313  bfi->AssemblePA(test_fes, trial_fes); // Reverse test and trial
314  }
315 
317  {
318  bfi->AssemblePAInteriorFaces(fes);
319  }
320 
322  {
323  bfi->AssemblePABoundaryFaces(fes);
324  }
325 
326  virtual void AddMultTransposePA(const Vector &x, Vector &y) const
327  {
328  bfi->AddMultPA(x, y);
329  }
330 
331  virtual void AddMultPA(const Vector& x, Vector& y) const
332  {
333  bfi->AddMultTransposePA(x, y);
334  }
335 
336  virtual void AssembleEA(const FiniteElementSpace &fes, Vector &emat,
337  const bool add);
338 
339  virtual void AssembleEAInteriorFaces(const FiniteElementSpace &fes,
340  Vector &ea_data_int,
341  Vector &ea_data_ext,
342  const bool add);
343 
344  virtual void AssembleEABoundaryFaces(const FiniteElementSpace &fes,
345  Vector &ea_data_bdr,
346  const bool add);
347 
348  virtual ~TransposeIntegrator() { if (own_bfi) { delete bfi; } }
349 };
350 
352 {
353 private:
354  int own_bfi;
356 
357 public:
358  LumpedIntegrator (BilinearFormIntegrator *bfi_, int own_bfi_ = 1)
359  { bfi = bfi_; own_bfi = own_bfi_; }
360 
361  virtual void SetIntRule(const IntegrationRule *ir);
362 
363  virtual void AssembleElementMatrix(const FiniteElement &el,
364  ElementTransformation &Trans,
365  DenseMatrix &elmat);
366 
367  virtual ~LumpedIntegrator() { if (own_bfi) { delete bfi; } }
368 };
369 
370 /// Integrator that inverts the matrix assembled by another integrator.
372 {
373 private:
374  int own_integrator;
375  BilinearFormIntegrator *integrator;
376 
377 public:
378  InverseIntegrator(BilinearFormIntegrator *integ, int own_integ = 1)
379  { integrator = integ; own_integrator = own_integ; }
380 
381  virtual void SetIntRule(const IntegrationRule *ir);
382 
383  virtual void AssembleElementMatrix(const FiniteElement &el,
384  ElementTransformation &Trans,
385  DenseMatrix &elmat);
386 
387  virtual ~InverseIntegrator() { if (own_integrator) { delete integrator; } }
388 };
389 
390 /// Integrator defining a sum of multiple Integrators.
392 {
393 private:
394  int own_integrators;
395  mutable DenseMatrix elem_mat;
396  Array<BilinearFormIntegrator*> integrators;
397 
398 public:
399  SumIntegrator(int own_integs = 1) { own_integrators = own_integs; }
400 
401  virtual void SetIntRule(const IntegrationRule *ir);
402 
404  { integrators.Append(integ); }
405 
406  virtual void AssembleElementMatrix(const FiniteElement &el,
407  ElementTransformation &Trans,
408  DenseMatrix &elmat);
409  virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
410  const FiniteElement &test_fe,
411  ElementTransformation &Trans,
412  DenseMatrix &elmat);
413 
415  virtual void AssembleFaceMatrix(const FiniteElement &el1,
416  const FiniteElement &el2,
418  DenseMatrix &elmat);
419 
420  virtual void AssembleFaceMatrix(const FiniteElement &trial_face_fe,
421  const FiniteElement &test_fe1,
422  const FiniteElement &test_fe2,
424  DenseMatrix &elmat);
425 
427  virtual void AssemblePA(const FiniteElementSpace& fes);
428 
429  virtual void AssembleDiagonalPA(Vector &diag);
430 
431  virtual void AssemblePAInteriorFaces(const FiniteElementSpace &fes);
432 
433  virtual void AssemblePABoundaryFaces(const FiniteElementSpace &fes);
434 
435  virtual void AddMultTransposePA(const Vector &x, Vector &y) const;
436 
437  virtual void AddMultPA(const Vector& x, Vector& y) const;
438 
439  virtual void AssembleMF(const FiniteElementSpace &fes);
440 
441  virtual void AddMultMF(const Vector &x, Vector &y) const;
442 
443  virtual void AddMultTransposeMF(const Vector &x, Vector &y) const;
444 
445  virtual void AssembleDiagonalMF(Vector &diag);
446 
447  virtual void AssembleEA(const FiniteElementSpace &fes, Vector &emat,
448  const bool add);
449 
450  virtual void AssembleEAInteriorFaces(const FiniteElementSpace &fes,
451  Vector &ea_data_int,
452  Vector &ea_data_ext,
453  const bool add);
454 
455  virtual void AssembleEABoundaryFaces(const FiniteElementSpace &fes,
456  Vector &ea_data_bdr,
457  const bool add);
458 
459  virtual ~SumIntegrator();
460 };
461 
462 /** An abstract class for integrating the product of two scalar basis functions
463  with an optional scalar coefficient. */
465 {
466 public:
467 
468  virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
469  const FiniteElement &test_fe,
470  ElementTransformation &Trans,
471  DenseMatrix &elmat);
472 
473  /// Support for use in BilinearForm. Can be used only when appropriate.
474  virtual void AssembleElementMatrix(const FiniteElement &fe,
475  ElementTransformation &Trans,
476  DenseMatrix &elmat)
477  { AssembleElementMatrix2(fe, fe, Trans, elmat); }
478 
479 protected:
480  /// This parameter can be set by derived methods to enable single shape
481  /// evaluation in case CalcTestShape() and CalcTrialShape() return the same
482  /// result if given the same FiniteElement. The default is false.
484 
487 
488  inline virtual bool VerifyFiniteElementTypes(
489  const FiniteElement & trial_fe,
490  const FiniteElement & test_fe) const
491  {
492  return (trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR &&
494  }
495 
496  inline virtual const char * FiniteElementTypeFailureMessage() const
497  {
498  return "MixedScalarIntegrator: "
499  "Trial and test spaces must both be scalar fields.";
500  }
501 
502  inline virtual int GetIntegrationOrder(const FiniteElement & trial_fe,
503  const FiniteElement & test_fe,
504  ElementTransformation &Trans)
505  { return trial_fe.GetOrder() + test_fe.GetOrder() + Trans.OrderW(); }
506 
507 
508  inline virtual void CalcTestShape(const FiniteElement & test_fe,
509  ElementTransformation &Trans,
510  Vector & shape)
511  { test_fe.CalcPhysShape(Trans, shape); }
512 
513  inline virtual void CalcTrialShape(const FiniteElement & trial_fe,
514  ElementTransformation &Trans,
515  Vector & shape)
516  { trial_fe.CalcPhysShape(Trans, shape); }
517 
519 
520 private:
521 
522 #ifndef MFEM_THREAD_SAFE
523  Vector test_shape;
524  Vector trial_shape;
525 #endif
526 
527 };
528 
529 /** An abstract class for integrating the inner product of two vector basis
530  functions with an optional scalar, vector, or matrix coefficient. */
532 {
533 public:
534 
535  virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
536  const FiniteElement &test_fe,
537  ElementTransformation &Trans,
538  DenseMatrix &elmat);
539 
540  /// Support for use in BilinearForm. Can be used only when appropriate.
541  virtual void AssembleElementMatrix(const FiniteElement &fe,
542  ElementTransformation &Trans,
543  DenseMatrix &elmat)
544  { AssembleElementMatrix2(fe, fe, Trans, elmat); }
545 
546 protected:
547  /// This parameter can be set by derived methods to enable single shape
548  /// evaluation in case CalcTestShape() and CalcTrialShape() return the same
549  /// result if given the same FiniteElement. The default is false.
551 
553  : same_calc_shape(false), Q(NULL), VQ(NULL), DQ(NULL), MQ(NULL) {}
555  : same_calc_shape(false), Q(&q), VQ(NULL), DQ(NULL), MQ(NULL) {}
557  : same_calc_shape(false), Q(NULL), VQ(diag?NULL:&vq), DQ(diag?&vq:NULL),
558  MQ(NULL) {}
560  : same_calc_shape(false), Q(NULL), VQ(NULL), DQ(NULL), MQ(&mq) {}
561 
562  inline virtual bool VerifyFiniteElementTypes(
563  const FiniteElement & trial_fe,
564  const FiniteElement & test_fe) const
565  {
566  return (trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
568  }
569 
570  inline virtual const char * FiniteElementTypeFailureMessage() const
571  {
572  return "MixedVectorIntegrator: "
573  "Trial and test spaces must both be vector fields";
574  }
575 
576  inline virtual int GetIntegrationOrder(const FiniteElement & trial_fe,
577  const FiniteElement & test_fe,
578  ElementTransformation &Trans)
579  { return trial_fe.GetOrder() + test_fe.GetOrder() + Trans.OrderW(); }
580 
581 
582  inline virtual int GetTestVDim(const FiniteElement & test_fe)
583  { return std::max(space_dim, test_fe.GetVDim()); }
584 
585  inline virtual void CalcTestShape(const FiniteElement & test_fe,
586  ElementTransformation &Trans,
587  DenseMatrix & shape)
588  { test_fe.CalcVShape(Trans, shape); }
589 
590  inline virtual int GetTrialVDim(const FiniteElement & trial_fe)
591  { return std::max(space_dim, trial_fe.GetVDim()); }
592 
593  inline virtual void CalcTrialShape(const FiniteElement & trial_fe,
594  ElementTransformation &Trans,
595  DenseMatrix & shape)
596  { trial_fe.CalcVShape(Trans, shape); }
597 
603 
604 private:
605 
606 #ifndef MFEM_THREAD_SAFE
607  Vector V;
608  Vector D;
609  DenseMatrix M;
610  DenseMatrix test_shape;
611  DenseMatrix trial_shape;
612  DenseMatrix shape_tmp;
613 #endif
614 
615 };
616 
617 /** An abstract class for integrating the product of a scalar basis function and
618  the inner product of a vector basis function with a vector coefficient. In
619  2D the inner product can be replaced with a cross product. */
621 {
622 public:
623 
624  virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
625  const FiniteElement &test_fe,
626  ElementTransformation &Trans,
627  DenseMatrix &elmat);
628 
629  /// Support for use in BilinearForm. Can be used only when appropriate.
630  /** Appropriate use cases are classes derived from
631  MixedScalarVectorIntegrator where the trial and test spaces can be the
632  same. Examples of such classes are: MixedVectorDivergenceIntegrator,
633  MixedScalarWeakDivergenceIntegrator, etc. */
634  virtual void AssembleElementMatrix(const FiniteElement &fe,
635  ElementTransformation &Trans,
636  DenseMatrix &elmat)
637  { AssembleElementMatrix2(fe, fe, Trans, elmat); }
638 
639 protected:
640 
641  MixedScalarVectorIntegrator(VectorCoefficient &vq, bool transpose_ = false,
642  bool cross_2d_ = false)
643  : VQ(&vq), transpose(transpose_), cross_2d(cross_2d_) {}
644 
645  inline virtual bool VerifyFiniteElementTypes(
646  const FiniteElement & trial_fe,
647  const FiniteElement & test_fe) const
648  {
649  return ((transpose &&
651  test_fe.GetRangeType() == mfem::FiniteElement::SCALAR ) ||
652  (!transpose &&
655  );
656  }
657 
658  inline virtual const char * FiniteElementTypeFailureMessage() const
659  {
660  if ( transpose )
661  {
662  return "MixedScalarVectorIntegrator: "
663  "Trial space must be a vector field "
664  "and the test space must be a scalar field";
665  }
666  else
667  {
668  return "MixedScalarVectorIntegrator: "
669  "Trial space must be a scalar field "
670  "and the test space must be a vector field";
671  }
672  }
673 
674  inline virtual int GetIntegrationOrder(const FiniteElement & trial_fe,
675  const FiniteElement & test_fe,
676  ElementTransformation &Trans)
677  { return trial_fe.GetOrder() + test_fe.GetOrder() + Trans.OrderW(); }
678 
679 
680  inline virtual int GetVDim(const FiniteElement & vector_fe)
681  { return std::max(space_dim, vector_fe.GetVDim()); }
682 
683  inline virtual void CalcVShape(const FiniteElement & vector_fe,
684  ElementTransformation &Trans,
685  DenseMatrix & shape_)
686  { vector_fe.CalcVShape(Trans, shape_); }
687 
688  inline virtual void CalcShape(const FiniteElement & scalar_fe,
689  ElementTransformation &Trans,
690  Vector & shape_)
691  { scalar_fe.CalcPhysShape(Trans, shape_); }
692 
695  bool transpose;
696  bool cross_2d; // In 2D use a cross product rather than a dot product
697 
698 private:
699 
700 #ifndef MFEM_THREAD_SAFE
701  Vector V;
702  DenseMatrix vshape;
703  Vector shape;
704  Vector vshape_tmp;
705 #endif
706 
707 };
708 
709 /** Class for integrating the bilinear form a(u,v) := (Q u, v) in either 1D, 2D,
710  or 3D and where Q is an optional scalar coefficient, u and v are each in H1
711  or L2. */
713 {
714 public:
717  : MixedScalarIntegrator(q) { same_calc_shape = true; }
718 };
719 
720 /** Class for integrating the bilinear form a(u,v) := (Q u, v) in either 2D, or
721  3D and where Q is a vector coefficient, u is in H1 or L2 and v is in H(Curl)
722  or H(Div). */
724 {
725 public:
728 };
729 
730 /** Class for integrating the bilinear form a(u,v) := (Q D u, v) in 1D where Q
731  is an optional scalar coefficient, u is in H1, and v is in L2. */
733 {
734 public:
737  : MixedScalarIntegrator(q) {}
738 
739 protected:
740  inline virtual bool VerifyFiniteElementTypes(
741  const FiniteElement & trial_fe,
742  const FiniteElement & test_fe) const
743  {
744  return (trial_fe.GetDim() == 1 && test_fe.GetDim() == 1 &&
745  trial_fe.GetDerivType() == mfem::FiniteElement::GRAD &&
747  }
748 
749  inline virtual const char * FiniteElementTypeFailureMessage() const
750  {
751  return "MixedScalarDerivativeIntegrator: "
752  "Trial and test spaces must both be scalar fields in 1D "
753  "and the trial space must implement CalcDShape.";
754  }
755 
756  inline virtual void CalcTrialShape(const FiniteElement & trial_fe,
757  ElementTransformation &Trans,
758  Vector & shape)
759  {
760  DenseMatrix dshape(shape.GetData(), shape.Size(), 1);
761  trial_fe.CalcPhysDShape(Trans, dshape);
762  }
763 };
764 
765 /** Class for integrating the bilinear form a(u,v) := -(Q u, D v) in 1D where Q
766  is an optional scalar coefficient, u is in L2, and v is in H1. */
768 {
769 public:
772  : MixedScalarIntegrator(q) {}
773 
774 protected:
775  inline virtual bool VerifyFiniteElementTypes(
776  const FiniteElement & trial_fe,
777  const FiniteElement & test_fe) const
778  {
779  return (trial_fe.GetDim() == 1 && test_fe.GetDim() == 1 &&
782  }
783 
784  inline virtual const char * FiniteElementTypeFailureMessage() const
785  {
786  return "MixedScalarWeakDerivativeIntegrator: "
787  "Trial and test spaces must both be scalar fields in 1D "
788  "and the test space must implement CalcDShape with "
789  "map type \"VALUE\".";
790  }
791 
792  inline virtual void CalcTestShape(const FiniteElement & test_fe,
793  ElementTransformation &Trans,
794  Vector & shape)
795  {
796  DenseMatrix dshape(shape.GetData(), shape.Size(), 1);
797  test_fe.CalcPhysDShape(Trans, dshape);
798  shape *= -1.0;
799  }
800 };
801 
802 /** Class for integrating the bilinear form a(u,v) := (Q div u, v) in either 2D
803  or 3D where Q is an optional scalar coefficient, u is in H(Div), and v is a
804  scalar field. */
806 {
807 public:
810  : MixedScalarIntegrator(q) {}
811 
812 protected:
813  inline virtual bool VerifyFiniteElementTypes(
814  const FiniteElement & trial_fe,
815  const FiniteElement & test_fe) const
816  {
817  return (trial_fe.GetDerivType() == mfem::FiniteElement::DIV &&
819  }
820 
821  inline virtual const char * FiniteElementTypeFailureMessage() const
822  {
823  return "MixedScalarDivergenceIntegrator: "
824  "Trial must be H(Div) and the test space must be a "
825  "scalar field";
826  }
827 
828  inline virtual int GetIntegrationOrder(const FiniteElement & trial_fe,
829  const FiniteElement & test_fe,
830  ElementTransformation &Trans)
831  { return trial_fe.GetOrder() + test_fe.GetOrder() + Trans.OrderW() - 1; }
832 
833  inline virtual void CalcTrialShape(const FiniteElement & trial_fe,
834  ElementTransformation &Trans,
835  Vector & shape)
836  { trial_fe.CalcPhysDivShape(Trans, shape); }
837 };
838 
839 /** Class for integrating the bilinear form a(u,v) := (V div u, v) in either 2D
840  or 3D where V is a vector coefficient, u is in H(Div), and v is a vector
841  field. */
843 {
844 public:
847 
848 protected:
849  inline virtual bool VerifyFiniteElementTypes(
850  const FiniteElement & trial_fe,
851  const FiniteElement & test_fe) const
852  {
853  return (trial_fe.GetDerivType() == mfem::FiniteElement::DIV &&
855  }
856 
857  inline virtual const char * FiniteElementTypeFailureMessage() const
858  {
859  return "MixedVectorDivergenceIntegrator: "
860  "Trial must be H(Div) and the test space must be a "
861  "vector field";
862  }
863 
864  // Subtract one due to the divergence and add one for the coefficient
865  // which is assumed to be at least linear.
866  inline virtual int GetIntegrationOrder(const FiniteElement & trial_fe,
867  const FiniteElement & test_fe,
868  ElementTransformation &Trans)
869  { return trial_fe.GetOrder() + test_fe.GetOrder() + Trans.OrderW() - 1 + 1; }
870 
871  inline virtual void CalcShape(const FiniteElement & scalar_fe,
872  ElementTransformation &Trans,
873  Vector & shape)
874  { scalar_fe.CalcPhysDivShape(Trans, shape); }
875 };
876 
877 /** Class for integrating the bilinear form a(u,v) := -(Q u, div v) in either 2D
878  or 3D where Q is an optional scalar coefficient, u is in L2 or H1, and v is
879  in H(Div). */
881 {
882 public:
885  : MixedScalarIntegrator(q) {}
886 
887 protected:
888  inline virtual bool VerifyFiniteElementTypes(
889  const FiniteElement & trial_fe,
890  const FiniteElement & test_fe) const
891  {
892  return (trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR &&
893  test_fe.GetDerivType() == mfem::FiniteElement::DIV );
894  }
895 
896  inline virtual const char * FiniteElementTypeFailureMessage() const
897  {
898  return "MixedScalarWeakGradientIntegrator: "
899  "Trial space must be a scalar field "
900  "and the test space must be H(Div)";
901  }
902 
903  inline virtual int GetIntegrationOrder(const FiniteElement & trial_fe,
904  const FiniteElement & test_fe,
905  ElementTransformation &Trans)
906  { return trial_fe.GetOrder() + test_fe.GetOrder() + Trans.OrderW() - 1; }
907 
908  virtual void CalcTestShape(const FiniteElement & test_fe,
909  ElementTransformation &Trans,
910  Vector & shape)
911  {
912  test_fe.CalcPhysDivShape(Trans, shape);
913  shape *= -1.0;
914  }
915 };
916 
917 /** Class for integrating the bilinear form a(u,v) := (Q curl u, v) in 2D where
918  Q is an optional scalar coefficient, u is in H(Curl), and v is in L2 or
919  H1. */
921 {
922 public:
925  : MixedScalarIntegrator(q) {}
926 
927 protected:
928  inline virtual bool VerifyFiniteElementTypes(
929  const FiniteElement & trial_fe,
930  const FiniteElement & test_fe) const
931  {
932  return (trial_fe.GetDim() == 2 && test_fe.GetDim() == 2 &&
933  trial_fe.GetDerivType() == mfem::FiniteElement::CURL &&
935  }
936 
937  inline virtual const char * FiniteElementTypeFailureMessage() const
938  {
939  return "MixedScalarCurlIntegrator: "
940  "Trial must be H(Curl) and the test space must be a "
941  "scalar field";
942  }
943 
944  inline virtual int GetIntegrationOrder(const FiniteElement & trial_fe,
945  const FiniteElement & test_fe,
946  ElementTransformation &Trans)
947  { return trial_fe.GetOrder() + test_fe.GetOrder() + Trans.OrderW() - 1; }
948 
949  inline virtual void CalcTrialShape(const FiniteElement & trial_fe,
950  ElementTransformation &Trans,
951  Vector & shape)
952  {
953  DenseMatrix dshape(shape.GetData(), shape.Size(), 1);
954  trial_fe.CalcPhysCurlShape(Trans, dshape);
955  }
956 
958  virtual void AssemblePA(const FiniteElementSpace &trial_fes,
959  const FiniteElementSpace &test_fes);
960 
961  virtual void AddMultPA(const Vector&, Vector&) const;
962  virtual void AddMultTransposePA(const Vector &x, Vector &y) const;
963 
964  // PA extension
966  const DofToQuad *mapsO; ///< Not owned. DOF-to-quad map, open.
967  const DofToQuad *mapsC; ///< Not owned. DOF-to-quad map, closed.
969 };
970 
971 /** Class for integrating the bilinear form a(u,v) := (Q u, curl v) in 2D where
972  Q is an optional scalar coefficient, u is in L2 or H1, and v is in
973  H(Curl). Partial assembly (PA) is supported but could be further optimized
974  by using more efficient threading and shared memory.
975 */
977 {
978 public:
981  : MixedScalarIntegrator(q) {}
982 
983 protected:
984  inline virtual bool VerifyFiniteElementTypes(
985  const FiniteElement & trial_fe,
986  const FiniteElement & test_fe) const
987  {
988  return (trial_fe.GetDim() == 2 && test_fe.GetDim() == 2 &&
991  }
992 
993  inline virtual const char * FiniteElementTypeFailureMessage() const
994  {
995  return "MixedScalarWeakCurlIntegrator: "
996  "Trial space must be a scalar field "
997  "and the test space must be H(Curl)";
998  }
999 
1000  inline virtual void CalcTestShape(const FiniteElement & test_fe,
1001  ElementTransformation &Trans,
1002  Vector & shape)
1003  {
1004  DenseMatrix dshape(shape.GetData(), shape.Size(), 1);
1005  test_fe.CalcPhysCurlShape(Trans, dshape);
1006  }
1007 };
1008 
1009 /** Class for integrating the bilinear form a(u,v) := (Q u, v) in either 2D or
1010  3D and where Q is an optional coefficient (of type scalar, matrix, or
1011  diagonal matrix) u and v are each in H(Curl) or H(Div). */
1013 {
1014 public:
1017  : MixedVectorIntegrator(q) { same_calc_shape = true; }
1019  : MixedVectorIntegrator(dq, true) { same_calc_shape = true; }
1021  : MixedVectorIntegrator(mq) { same_calc_shape = true; }
1022 };
1023 
1024 /** Class for integrating the bilinear form a(u,v) := (V x u, v) in 3D and where
1025  V is a vector coefficient u and v are each in H(Curl) or H(Div). */
1027 {
1028 public:
1030  : MixedVectorIntegrator(vq, false) { same_calc_shape = true; }
1031 };
1032 
1033 /** Class for integrating the bilinear form a(u,v) := (V . u, v) in 2D or 3D and
1034  where V is a vector coefficient u is in H(Curl) or H(Div) and v is in H1 or
1035  L2. */
1037 {
1038 public:
1040  : MixedScalarVectorIntegrator(vq, true) {}
1041 
1042  inline virtual bool VerifyFiniteElementTypes(
1043  const FiniteElement & trial_fe,
1044  const FiniteElement & test_fe) const
1045  {
1046  return (trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
1048  }
1049 
1050  inline virtual const char * FiniteElementTypeFailureMessage() const
1051  {
1052  return "MixedDotProductIntegrator: "
1053  "Trial space must be a vector field "
1054  "and the test space must be a scalar field";
1055  }
1056 };
1057 
1058 /** Class for integrating the bilinear form a(u,v) := (-V . u, Div v) in 2D or
1059  3D and where V is a vector coefficient u is in H(Curl) or H(Div) and v is in
1060  RT. */
1062 {
1063 public:
1065  : MixedScalarVectorIntegrator(vq, true) {}
1066 
1067  inline virtual bool VerifyFiniteElementTypes(
1068  const FiniteElement & trial_fe,
1069  const FiniteElement & test_fe) const
1070  {
1071  return (trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
1073  test_fe.GetDerivType() == mfem::FiniteElement::DIV );
1074  }
1075 
1076  inline virtual const char * FiniteElementTypeFailureMessage() const
1077  {
1078  return "MixedWeakGradDotIntegrator: "
1079  "Trial space must be a vector field "
1080  "and the test space must be a vector field with a divergence";
1081  }
1082 
1083  // Subtract one due to the gradient and add one for the coefficient
1084  // which is assumed to be at least linear.
1085  inline virtual int GetIntegrationOrder(const FiniteElement & trial_fe,
1086  const FiniteElement & test_fe,
1087  ElementTransformation &Trans)
1088  { return trial_fe.GetOrder() + test_fe.GetOrder() + Trans.OrderW() - 1 + 1; }
1089 
1090  inline virtual void CalcShape(const FiniteElement & scalar_fe,
1091  ElementTransformation &Trans,
1092  Vector & shape)
1093  { scalar_fe.CalcPhysDivShape(Trans, shape); shape *= -1.0; }
1094 };
1095 
1096 /** Class for integrating the bilinear form a(u,v) := (V x u, Grad v) in 3D and
1097  where V is a vector coefficient u is in H(Curl) or H(Div) and v is in H1. */
1099 {
1100 public:
1102  : MixedVectorIntegrator(vq, false) {}
1103 
1104  inline virtual bool VerifyFiniteElementTypes(
1105  const FiniteElement & trial_fe,
1106  const FiniteElement & test_fe) const
1107  {
1108  return (trial_fe.GetVDim() == 3 &&
1109  trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
1111  test_fe.GetDerivType() == mfem::FiniteElement::GRAD );
1112  }
1113 
1114  inline virtual const char * FiniteElementTypeFailureMessage() const
1115  {
1116  return "MixedWeakDivCrossIntegrator: "
1117  "Trial space must be a vector field in 3D "
1118  "and the test space must be a scalar field with a gradient";
1119  }
1120 
1121  inline virtual int GetTestVDim(const FiniteElement & test_fe)
1122  { return space_dim; }
1123 
1124  inline virtual void CalcTestShape(const FiniteElement & test_fe,
1125  ElementTransformation &Trans,
1126  DenseMatrix & shape)
1127  { test_fe.CalcPhysDShape(Trans, shape); shape *= -1.0; }
1128 };
1129 
1130 /** Class for integrating the bilinear form a(u,v) := (Q Grad u, Grad v) in 3D
1131  or in 2D and where Q is a scalar or matrix coefficient u and v are both in
1132  H1. */
1134 {
1135 public:
1138  : MixedVectorIntegrator(q) { same_calc_shape = true; }
1140  : MixedVectorIntegrator(dq, true) { same_calc_shape = true; }
1142  : MixedVectorIntegrator(mq) { same_calc_shape = true; }
1143 
1144  inline virtual bool VerifyFiniteElementTypes(
1145  const FiniteElement & trial_fe,
1146  const FiniteElement & test_fe) const
1147  {
1148  return (trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR &&
1149  trial_fe.GetDerivType() == mfem::FiniteElement::GRAD &&
1151  test_fe.GetDerivType() == mfem::FiniteElement::GRAD );
1152  }
1153 
1154  inline virtual const char * FiniteElementTypeFailureMessage() const
1155  {
1156  return "MixedGradGradIntegrator: "
1157  "Trial and test spaces must both be scalar fields "
1158  "with a gradient operator.";
1159  }
1160 
1161  inline virtual int GetIntegrationOrder(const FiniteElement & trial_fe,
1162  const FiniteElement & test_fe,
1163  ElementTransformation &Trans)
1164  {
1165  // Same as DiffusionIntegrator
1166  return test_fe.Space() == FunctionSpace::Pk ?
1167  trial_fe.GetOrder() + test_fe.GetOrder() - 2 :
1168  trial_fe.GetOrder() + test_fe.GetOrder() + test_fe.GetDim() - 1;
1169  }
1170 
1171  inline virtual int GetTrialVDim(const FiniteElement & trial_fe)
1172  { return space_dim; }
1173 
1174  inline virtual void CalcTrialShape(const FiniteElement & trial_fe,
1175  ElementTransformation &Trans,
1176  DenseMatrix & shape)
1177  { trial_fe.CalcPhysDShape(Trans, shape); }
1178 
1179  inline virtual int GetTestVDim(const FiniteElement & test_fe)
1180  { return space_dim; }
1181 
1182  inline virtual void CalcTestShape(const FiniteElement & test_fe,
1183  ElementTransformation &Trans,
1184  DenseMatrix & shape)
1185  { test_fe.CalcPhysDShape(Trans, shape); }
1186 };
1187 
1188 /** Class for integrating the bilinear form a(u,v) := (V x Grad u, Grad v) in 3D
1189  or in 2D and where V is a vector coefficient u and v are both in H1. */
1191 {
1192 public:
1194  : MixedVectorIntegrator(vq, false) { same_calc_shape = true; }
1195 
1196  inline virtual bool VerifyFiniteElementTypes(
1197  const FiniteElement & trial_fe,
1198  const FiniteElement & test_fe) const
1199  {
1200  return (trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR &&
1201  trial_fe.GetDerivType() == mfem::FiniteElement::GRAD &&
1203  test_fe.GetDerivType() == mfem::FiniteElement::GRAD );
1204  }
1205 
1206  inline virtual const char * FiniteElementTypeFailureMessage() const
1207  {
1208  return "MixedCrossGradGradIntegrator: "
1209  "Trial and test spaces must both be scalar fields "
1210  "with a gradient operator.";
1211  }
1212 
1213  inline virtual int GetTrialVDim(const FiniteElement & trial_fe)
1214  { return space_dim; }
1215 
1216  inline virtual void CalcTrialShape(const FiniteElement & trial_fe,
1217  ElementTransformation &Trans,
1218  DenseMatrix & shape)
1219  { trial_fe.CalcPhysDShape(Trans, shape); }
1220 
1221  inline virtual int GetTestVDim(const FiniteElement & test_fe)
1222  { return space_dim; }
1223 
1224  inline virtual void CalcTestShape(const FiniteElement & test_fe,
1225  ElementTransformation &Trans,
1226  DenseMatrix & shape)
1227  { test_fe.CalcPhysDShape(Trans, shape); }
1228 };
1229 
1230 /** Class for integrating the bilinear form a(u,v) := (Q Curl u, Curl v) in 3D
1231  and where Q is a scalar or matrix coefficient u and v are both in
1232  H(Curl). */
1234 {
1235 public:
1238  : MixedVectorIntegrator(q) { same_calc_shape = true; }
1240  : MixedVectorIntegrator(dq, true) { same_calc_shape = true; }
1242  : MixedVectorIntegrator(mq) { same_calc_shape = true; }
1243 
1244  inline virtual bool VerifyFiniteElementTypes(
1245  const FiniteElement & trial_fe,
1246  const FiniteElement & test_fe) const
1247  {
1248  return (trial_fe.GetCurlDim() == 3 && test_fe.GetCurlDim() == 3 &&
1249  trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
1250  trial_fe.GetDerivType() == mfem::FiniteElement::CURL &&
1252  test_fe.GetDerivType() == mfem::FiniteElement::CURL );
1253  }
1254 
1255  inline virtual const char * FiniteElementTypeFailureMessage() const
1256  {
1257  return "MixedCurlCurlIntegrator"
1258  "Trial and test spaces must both be vector fields in 3D "
1259  "with a curl.";
1260  }
1261 
1262  inline virtual int GetTrialVDim(const FiniteElement & trial_fe)
1263  { return trial_fe.GetCurlDim(); }
1264 
1265  inline virtual void CalcTrialShape(const FiniteElement & trial_fe,
1266  ElementTransformation &Trans,
1267  DenseMatrix & shape)
1268  { trial_fe.CalcPhysCurlShape(Trans, shape); }
1269 
1270  inline virtual int GetTestVDim(const FiniteElement & test_fe)
1271  { return test_fe.GetCurlDim(); }
1272 
1273  inline virtual void CalcTestShape(const FiniteElement & test_fe,
1274  ElementTransformation &Trans,
1275  DenseMatrix & shape)
1276  { test_fe.CalcPhysCurlShape(Trans, shape); }
1277 };
1278 
1279 /** Class for integrating the bilinear form a(u,v) := (V x Curl u, Curl v) in 3D
1280  and where V is a vector coefficient u and v are both in H(Curl). */
1282 {
1283 public:
1285  : MixedVectorIntegrator(vq, false) { same_calc_shape = true; }
1286 
1287  inline virtual bool VerifyFiniteElementTypes(
1288  const FiniteElement & trial_fe,
1289  const FiniteElement & test_fe) const
1290  {
1291  return (trial_fe.GetCurlDim() == 3 && trial_fe.GetVDim() == 3 &&
1292  test_fe.GetCurlDim() == 3 && test_fe.GetVDim() == 3 &&
1293  trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
1294  trial_fe.GetDerivType() == mfem::FiniteElement::CURL &&
1296  test_fe.GetDerivType() == mfem::FiniteElement::CURL );
1297  }
1298 
1299  inline virtual const char * FiniteElementTypeFailureMessage() const
1300  {
1301  return "MixedCrossCurlCurlIntegrator: "
1302  "Trial and test spaces must both be vector fields in 3D "
1303  "with a curl.";
1304  }
1305 
1306  inline virtual int GetTrialVDim(const FiniteElement & trial_fe)
1307  { return trial_fe.GetCurlDim(); }
1308 
1309  inline virtual void CalcTrialShape(const FiniteElement & trial_fe,
1310  ElementTransformation &Trans,
1311  DenseMatrix & shape)
1312  { trial_fe.CalcPhysCurlShape(Trans, shape); }
1313 
1314  inline virtual int GetTestVDim(const FiniteElement & test_fe)
1315  { return test_fe.GetCurlDim(); }
1316 
1317  inline virtual void CalcTestShape(const FiniteElement & test_fe,
1318  ElementTransformation &Trans,
1319  DenseMatrix & shape)
1320  { test_fe.CalcPhysCurlShape(Trans, shape); }
1321 };
1322 
1323 /** Class for integrating the bilinear form a(u,v) := (V x Curl u, Grad v) in 3D
1324  and where V is a vector coefficient u is in H(Curl) and v is in H1. */
1326 {
1327 public:
1329  : MixedVectorIntegrator(vq, false) {}
1330 
1331  inline virtual bool VerifyFiniteElementTypes(
1332  const FiniteElement & trial_fe,
1333  const FiniteElement & test_fe) const
1334  {
1335  return (trial_fe.GetCurlDim() == 3 &&
1336  trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
1337  trial_fe.GetDerivType() == mfem::FiniteElement::CURL &&
1339  test_fe.GetDerivType() == mfem::FiniteElement::GRAD );
1340  }
1341 
1342  inline virtual const char * FiniteElementTypeFailureMessage() const
1343  {
1344  return "MixedCrossCurlGradIntegrator"
1345  "Trial space must be a vector field in 3D with a curl"
1346  "and the test space must be a scalar field with a gradient";
1347  }
1348 
1349  inline virtual int GetTrialVDim(const FiniteElement & trial_fe)
1350  { return trial_fe.GetCurlDim(); }
1351 
1352  inline virtual void CalcTrialShape(const FiniteElement & trial_fe,
1353  ElementTransformation &Trans,
1354  DenseMatrix & shape)
1355  { trial_fe.CalcPhysCurlShape(Trans, shape); }
1356 
1357  inline virtual int GetTestVDim(const FiniteElement & test_fe)
1358  { return space_dim; }
1359 
1360  inline virtual void CalcTestShape(const FiniteElement & test_fe,
1361  ElementTransformation &Trans,
1362  DenseMatrix & shape)
1363  { test_fe.CalcPhysDShape(Trans, shape); }
1364 };
1365 
1366 /** Class for integrating the bilinear form a(u,v) := (V x Grad u, Curl v) in 3D
1367  and where V is a scalar coefficient u is in H1 and v is in H(Curl). */
1369 {
1370 public:
1372  : MixedVectorIntegrator(vq, false) {}
1373 
1374  inline virtual bool VerifyFiniteElementTypes(
1375  const FiniteElement & trial_fe,
1376  const FiniteElement & test_fe) const
1377  {
1378  return (test_fe.GetCurlDim() == 3 &&
1379  trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR &&
1380  trial_fe.GetDerivType() == mfem::FiniteElement::GRAD &&
1382  test_fe.GetDerivType() == mfem::FiniteElement::CURL );
1383  }
1384 
1385  inline virtual const char * FiniteElementTypeFailureMessage() const
1386  {
1387  return "MixedCrossGradCurlIntegrator"
1388  "Trial space must be a scalar field in 3D with a gradient"
1389  "and the test space must be a vector field with a curl";
1390  }
1391 
1392  inline virtual int GetTrialVDim(const FiniteElement & trial_fe)
1393  { return space_dim; }
1394 
1395  inline virtual void CalcTrialShape(const FiniteElement & trial_fe,
1396  ElementTransformation &Trans,
1397  DenseMatrix & shape)
1398  { trial_fe.CalcPhysDShape(Trans, shape); }
1399 
1400  inline virtual int GetTestVDim(const FiniteElement & test_fe)
1401  { return test_fe.GetCurlDim(); }
1402 
1403  inline virtual void CalcTestShape(const FiniteElement & test_fe,
1404  ElementTransformation &Trans,
1405  DenseMatrix & shape)
1406  { test_fe.CalcPhysCurlShape(Trans, shape); }
1407 };
1408 
1409 /** Class for integrating the bilinear form a(u,v) := (V x u, Curl v) in 3D and
1410  where V is a vector coefficient u is in H(Curl) or H(Div) and v is in
1411  H(Curl). */
1413 {
1414 public:
1416  : MixedVectorIntegrator(vq, false) {}
1417 
1418  inline virtual bool VerifyFiniteElementTypes(
1419  const FiniteElement & trial_fe,
1420  const FiniteElement & test_fe) const
1421  {
1422  return (trial_fe.GetVDim() == 3 && test_fe.GetCurlDim() == 3 &&
1423  trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
1425  test_fe.GetDerivType() == mfem::FiniteElement::CURL );
1426  }
1427 
1428  inline virtual const char * FiniteElementTypeFailureMessage() const
1429  {
1430  return "MixedWeakCurlCrossIntegrator: "
1431  "Trial space must be a vector field in 3D "
1432  "and the test space must be a vector field with a curl";
1433  }
1434 
1435  inline virtual int GetTestVDim(const FiniteElement & test_fe)
1436  { return test_fe.GetCurlDim(); }
1437 
1438  inline virtual void CalcTestShape(const FiniteElement & test_fe,
1439  ElementTransformation &Trans,
1440  DenseMatrix & shape)
1441  { test_fe.CalcPhysCurlShape(Trans, shape); }
1442 };
1443 
1444 /** Class for integrating the bilinear form a(u,v) := (V x u, Curl v) in 2D and
1445  where V is a vector coefficient u is in H(Curl) or H(Div) and v is in
1446  H(Curl). */
1448 {
1449 public:
1451  : MixedScalarVectorIntegrator(vq, true, true) {}
1452 
1453  inline virtual bool VerifyFiniteElementTypes(
1454  const FiniteElement & trial_fe,
1455  const FiniteElement & test_fe) const
1456  {
1457  return (trial_fe.GetDim() == 2 && test_fe.GetDim() == 2 &&
1458  trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
1460  test_fe.GetDerivType() == mfem::FiniteElement::CURL );
1461  }
1462 
1463  inline virtual const char * FiniteElementTypeFailureMessage() const
1464  {
1465  return "MixedScalarWeakCurlCrossIntegrator: "
1466  "Trial space must be a vector field in 2D "
1467  "and the test space must be a vector field with a curl";
1468  }
1469 
1470  inline virtual void CalcShape(const FiniteElement & scalar_fe,
1471  ElementTransformation &Trans,
1472  Vector & shape)
1473  {
1474  DenseMatrix dshape(shape.GetData(), shape.Size(), 1);
1475  scalar_fe.CalcPhysCurlShape(Trans, dshape);
1476  }
1477 };
1478 
1479 /** Class for integrating the bilinear form a(u,v) := (V x Grad u, v) in 3D or
1480  in 2D and where V is a vector coefficient u is in H1 and v is in H(Curl) or
1481  H(Div). */
1483 {
1484 public:
1486  : MixedVectorIntegrator(vq, false) {}
1487 
1488  inline virtual bool VerifyFiniteElementTypes(
1489  const FiniteElement & trial_fe,
1490  const FiniteElement & test_fe) const
1491  {
1492  return (test_fe.GetVDim() == 3 &&
1493  trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR &&
1494  trial_fe.GetDerivType() == mfem::FiniteElement::GRAD &&
1496  }
1497 
1498  inline virtual const char * FiniteElementTypeFailureMessage() const
1499  {
1500  return "MixedCrossGradIntegrator: "
1501  "Trial space must be a scalar field with a gradient operator"
1502  " and the test space must be a vector field both in 3D.";
1503  }
1504 
1505  inline virtual int GetTrialVDim(const FiniteElement & trial_fe)
1506  { return space_dim; }
1507 
1508  inline virtual void CalcTrialShape(const FiniteElement & trial_fe,
1509  ElementTransformation &Trans,
1510  DenseMatrix & shape)
1511  { trial_fe.CalcPhysDShape(Trans, shape); }
1512 
1513  inline virtual void CalcTestShape(const FiniteElement & test_fe,
1514  ElementTransformation &Trans,
1515  DenseMatrix & shape)
1516  { test_fe.CalcVShape(Trans, shape); }
1517 };
1518 
1519 /** Class for integrating the bilinear form a(u,v) := (V x Curl u, v) in 3D and
1520  where V is a vector coefficient u is in H(Curl) and v is in H(Curl) or
1521  H(Div). */
1523 {
1524 public:
1526  : MixedVectorIntegrator(vq, false) {}
1527 
1528  inline virtual bool VerifyFiniteElementTypes(
1529  const FiniteElement & trial_fe,
1530  const FiniteElement & test_fe) const
1531  {
1532  return (trial_fe.GetCurlDim() == 3 && test_fe.GetVDim() == 3 &&
1533  trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
1534  trial_fe.GetDerivType() == mfem::FiniteElement::CURL &&
1536  }
1537 
1538  inline virtual const char * FiniteElementTypeFailureMessage() const
1539  {
1540  return "MixedCrossCurlIntegrator: "
1541  "Trial space must be a vector field in 3D with a curl "
1542  "and the test space must be a vector field";
1543  }
1544 
1545  inline virtual int GetTrialVDim(const FiniteElement & trial_fe)
1546  { return trial_fe.GetCurlDim(); }
1547 
1548  inline virtual void CalcTrialShape(const FiniteElement & trial_fe,
1549  ElementTransformation &Trans,
1550  DenseMatrix & shape)
1551  { trial_fe.CalcPhysCurlShape(Trans, shape); }
1552 };
1553 
1554 /** Class for integrating the bilinear form a(u,v) := (V x Curl u, v) in 2D and
1555  where V is a vector coefficient u is in H(Curl) and v is in H(Curl) or
1556  H(Div). */
1558 {
1559 public:
1561  : MixedScalarVectorIntegrator(vq, false, true) {}
1562 
1563  inline virtual bool VerifyFiniteElementTypes(
1564  const FiniteElement & trial_fe,
1565  const FiniteElement & test_fe) const
1566  {
1567  return (trial_fe.GetDim() == 2 && test_fe.GetDim() == 2 &&
1568  trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
1569  trial_fe.GetDerivType() == mfem::FiniteElement::CURL &&
1571  }
1572 
1573  inline virtual const char * FiniteElementTypeFailureMessage() const
1574  {
1575  return "MixedCrossCurlIntegrator: "
1576  "Trial space must be a vector field in 2D with a curl "
1577  "and the test space must be a vector field";
1578  }
1579 
1580  inline virtual void CalcShape(const FiniteElement & scalar_fe,
1581  ElementTransformation &Trans,
1582  Vector & shape)
1583  {
1584  DenseMatrix dshape(shape.GetData(), shape.Size(), 1);
1585  scalar_fe.CalcPhysCurlShape(Trans, dshape); shape *= -1.0;
1586  }
1587 };
1588 
1589 /** Class for integrating the bilinear form a(u,v) := (V x Grad u, v) in 2D and
1590  where V is a vector coefficient u is in H1 and v is in H1 or L2. */
1592 {
1593 public:
1595  : MixedScalarVectorIntegrator(vq, true, true) {}
1596 
1597  inline virtual bool VerifyFiniteElementTypes(
1598  const FiniteElement & trial_fe,
1599  const FiniteElement & test_fe) const
1600  {
1601  return (trial_fe.GetDim() == 2 && test_fe.GetDim() == 2 &&
1602  trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR &&
1603  trial_fe.GetDerivType() == mfem::FiniteElement::GRAD &&
1605  }
1606 
1607  inline virtual const char * FiniteElementTypeFailureMessage() const
1608  {
1609  return "MixedScalarCrossGradIntegrator: "
1610  "Trial space must be a scalar field in 2D with a gradient "
1611  "and the test space must be a scalar field";
1612  }
1613 
1614  inline int GetVDim(const FiniteElement & vector_fe)
1615  { return space_dim; }
1616 
1617  inline virtual void CalcVShape(const FiniteElement & vector_fe,
1618  ElementTransformation &Trans,
1619  DenseMatrix & shape)
1620  { vector_fe.CalcPhysDShape(Trans, shape); }
1621 };
1622 
1623 /** Class for integrating the bilinear form a(u,v) := (V x u, v) in 2D and where
1624  V is a vector coefficient u is in ND or RT and v is in H1 or L2. */
1626 {
1627 public:
1629  : MixedScalarVectorIntegrator(vq, true, true) {}
1630 
1631  inline virtual bool VerifyFiniteElementTypes(
1632  const FiniteElement & trial_fe,
1633  const FiniteElement & test_fe) const
1634  {
1635  return (trial_fe.GetDim() == 2 && test_fe.GetDim() == 2 &&
1636  trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
1638  }
1639 
1640  inline virtual const char * FiniteElementTypeFailureMessage() const
1641  {
1642  return "MixedScalarCrossProductIntegrator: "
1643  "Trial space must be a vector field in 2D "
1644  "and the test space must be a scalar field";
1645  }
1646 };
1647 
1648 /** Class for integrating the bilinear form a(u,v) := (V x z u, v) in 2D and
1649  where V is a vector coefficient u is in H1 or L2 and v is in ND or RT. */
1651 {
1652 public:
1654  : MixedScalarVectorIntegrator(vq, false, true) {}
1655 
1656  inline virtual bool VerifyFiniteElementTypes(
1657  const FiniteElement & trial_fe,
1658  const FiniteElement & test_fe) const
1659  {
1660  return (trial_fe.GetDim() == 2 && test_fe.GetDim() == 2 &&
1661  trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR &&
1663  }
1664 
1665  inline virtual const char * FiniteElementTypeFailureMessage() const
1666  {
1667  return "MixedScalarWeakCrossProductIntegrator: "
1668  "Trial space must be a scalar field in 2D "
1669  "and the test space must be a vector field";
1670  }
1671 
1672  inline virtual void CalcShape(const FiniteElement & scalar_fe,
1673  ElementTransformation &Trans,
1674  Vector & shape)
1675  { scalar_fe.CalcPhysShape(Trans, shape); shape *= -1.0; }
1676 };
1677 
1678 /** Class for integrating the bilinear form a(u,v) := (V . Grad u, v) in 2D or
1679  3D and where V is a vector coefficient, u is in H1 and v is in H1 or L2. */
1681 {
1682 public:
1684  : MixedScalarVectorIntegrator(vq, true) {}
1685 
1686  inline virtual bool VerifyFiniteElementTypes(
1687  const FiniteElement & trial_fe,
1688  const FiniteElement & test_fe) const
1689  {
1690  return (trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR &&
1691  trial_fe.GetDerivType() == mfem::FiniteElement::GRAD &&
1693  }
1694 
1695  inline virtual const char * FiniteElementTypeFailureMessage() const
1696  {
1697  return "MixedDirectionalDerivativeIntegrator: "
1698  "Trial space must be a scalar field with a gradient "
1699  "and the test space must be a scalar field";
1700  }
1701 
1702  inline virtual int GetVDim(const FiniteElement & vector_fe)
1703  { return space_dim; }
1704 
1705  inline virtual void CalcVShape(const FiniteElement & vector_fe,
1706  ElementTransformation &Trans,
1707  DenseMatrix & shape)
1708  { vector_fe.CalcPhysDShape(Trans, shape); }
1709 };
1710 
1711 /** Class for integrating the bilinear form a(u,v) := (-V . Grad u, Div v) in 2D
1712  or 3D and where V is a vector coefficient, u is in H1 and v is in RT. */
1714 {
1715 public:
1717  : MixedScalarVectorIntegrator(vq, true) {}
1718 
1719  inline virtual bool VerifyFiniteElementTypes(
1720  const FiniteElement & trial_fe,
1721  const FiniteElement & test_fe) const
1722  {
1723  return (trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR &&
1724  trial_fe.GetDerivType() == mfem::FiniteElement::GRAD &&
1726  test_fe.GetDerivType() == mfem::FiniteElement::DIV );
1727  }
1728 
1729  inline virtual const char * FiniteElementTypeFailureMessage() const
1730  {
1731  return "MixedGradDivIntegrator: "
1732  "Trial space must be a scalar field with a gradient"
1733  "and the test space must be a vector field with a divergence";
1734  }
1735 
1736  inline virtual int GetVDim(const FiniteElement & vector_fe)
1737  { return space_dim; }
1738 
1739  inline virtual void CalcVShape(const FiniteElement & vector_fe,
1740  ElementTransformation &Trans,
1741  DenseMatrix & shape)
1742  { vector_fe.CalcPhysDShape(Trans, shape); shape *= -1.0; }
1743 
1744  inline virtual void CalcShape(const FiniteElement & scalar_fe,
1745  ElementTransformation &Trans,
1746  Vector & shape)
1747  { scalar_fe.CalcPhysDivShape(Trans, shape); }
1748 };
1749 
1750 /** Class for integrating the bilinear form a(u,v) := (-V Div u, Grad v) in 2D
1751  or 3D and where V is a vector coefficient, u is in RT and v is in H1. */
1753 {
1754 public:
1756  : MixedScalarVectorIntegrator(vq, false) {}
1757 
1758  inline virtual bool VerifyFiniteElementTypes(
1759  const FiniteElement & trial_fe,
1760  const FiniteElement & test_fe) const
1761  {
1762  return (trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
1763  trial_fe.GetDerivType() == mfem::FiniteElement::DIV &&
1766  );
1767  }
1768 
1769  inline virtual const char * FiniteElementTypeFailureMessage() const
1770  {
1771  return "MixedDivGradIntegrator: "
1772  "Trial space must be a vector field with a divergence"
1773  "and the test space must be a scalar field with a gradient";
1774  }
1775 
1776  inline virtual int GetVDim(const FiniteElement & vector_fe)
1777  { return space_dim; }
1778 
1779  inline virtual void CalcVShape(const FiniteElement & vector_fe,
1780  ElementTransformation &Trans,
1781  DenseMatrix & shape)
1782  { vector_fe.CalcPhysDShape(Trans, shape); shape *= -1.0; }
1783 
1784  inline virtual void CalcShape(const FiniteElement & scalar_fe,
1785  ElementTransformation &Trans,
1786  Vector & shape)
1787  { scalar_fe.CalcPhysDivShape(Trans, shape); }
1788 };
1789 
1790 /** Class for integrating the bilinear form a(u,v) := (-V u, Grad v) in 2D or 3D
1791  and where V is a vector coefficient, u is in H1 or L2 and v is in H1. */
1793 {
1794 public:
1796  : MixedScalarVectorIntegrator(vq, false) {}
1797 
1798  inline virtual bool VerifyFiniteElementTypes(
1799  const FiniteElement & trial_fe,
1800  const FiniteElement & test_fe) const
1801  {
1802  return (trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR &&
1804  test_fe.GetDerivType() == mfem::FiniteElement::GRAD );
1805  }
1806 
1807  inline virtual const char * FiniteElementTypeFailureMessage() const
1808  {
1809  return "MixedScalarWeakDivergenceIntegrator: "
1810  "Trial space must be a scalar field "
1811  "and the test space must be a scalar field with a gradient";
1812  }
1813 
1814  inline int GetVDim(const FiniteElement & vector_fe)
1815  { return space_dim; }
1816 
1817  inline virtual void CalcVShape(const FiniteElement & vector_fe,
1818  ElementTransformation &Trans,
1819  DenseMatrix & shape)
1820  { vector_fe.CalcPhysDShape(Trans, shape); shape *= -1.0; }
1821 };
1822 
1823 /** Class for integrating the bilinear form a(u,v) := (Q grad u, v) in either 2D
1824  or 3D and where Q is an optional coefficient (of type scalar, matrix, or
1825  diagonal matrix) u is in H1 and v is in H(Curl) or H(Div). Partial assembly
1826  (PA) is supported but could be further optimized by using more efficient
1827  threading and shared memory.
1828 */
1830 {
1831 public:
1834  : MixedVectorIntegrator(q) {}
1836  : MixedVectorIntegrator(dq, true) {}
1838  : MixedVectorIntegrator(mq) {}
1839 
1840 protected:
1841  inline virtual bool VerifyFiniteElementTypes(
1842  const FiniteElement & trial_fe,
1843  const FiniteElement & test_fe) const
1844  {
1845  return (trial_fe.GetDerivType() == mfem::FiniteElement::GRAD &&
1847  }
1848 
1849  inline virtual const char * FiniteElementTypeFailureMessage() const
1850  {
1851  return "MixedVectorGradientIntegrator: "
1852  "Trial spaces must be H1 and the test space must be a "
1853  "vector field in 2D or 3D";
1854  }
1855 
1856  inline virtual int GetTrialVDim(const FiniteElement & trial_fe)
1857  { return space_dim; }
1858 
1859  inline virtual void CalcTrialShape(const FiniteElement & trial_fe,
1860  ElementTransformation &Trans,
1861  DenseMatrix & shape)
1862  {
1863  trial_fe.CalcPhysDShape(Trans, shape);
1864  }
1865 
1867  virtual void AssemblePA(const FiniteElementSpace &trial_fes,
1868  const FiniteElementSpace &test_fes);
1869 
1870  virtual void AddMultPA(const Vector&, Vector&) const;
1871  virtual void AddMultTransposePA(const Vector&, Vector&) const;
1872 
1873 private:
1874  DenseMatrix Jinv;
1875 
1876  // PA extension
1877  Vector pa_data;
1878  const DofToQuad *mapsO; ///< Not owned. DOF-to-quad map, open.
1879  const DofToQuad *mapsC; ///< Not owned. DOF-to-quad map, closed.
1880  const GeometricFactors *geom; ///< Not owned
1881  int dim, ne, dofs1D, quad1D;
1882 };
1883 
1884 /** Class for integrating the bilinear form a(u,v) := (Q curl u, v) in 3D and
1885  where Q is an optional coefficient (of type scalar, matrix, or diagonal
1886  matrix) u is in H(Curl) and v is in H(Div) or H(Curl). */
1888 {
1889 public:
1892  : MixedVectorIntegrator(q) {}
1894  : MixedVectorIntegrator(dq, true) {}
1896  : MixedVectorIntegrator(mq) {}
1897 
1898 protected:
1899  inline virtual bool VerifyFiniteElementTypes(
1900  const FiniteElement & trial_fe,
1901  const FiniteElement & test_fe) const
1902  {
1903  return (trial_fe.GetCurlDim() == 3 && test_fe.GetVDim() == 3 &&
1904  trial_fe.GetDerivType() == mfem::FiniteElement::CURL &&
1906  }
1907 
1908  inline virtual const char * FiniteElementTypeFailureMessage() const
1909  {
1910  return "MixedVectorCurlIntegrator: "
1911  "Trial space must be H(Curl) and the test space must be a "
1912  "vector field in 3D";
1913  }
1914 
1915  inline virtual int GetTrialVDim(const FiniteElement & trial_fe)
1916  { return trial_fe.GetCurlDim(); }
1917 
1918  inline virtual void CalcTrialShape(const FiniteElement & trial_fe,
1919  ElementTransformation &Trans,
1920  DenseMatrix & shape)
1921  {
1922  trial_fe.CalcPhysCurlShape(Trans, shape);
1923  }
1924 
1926  virtual void AssemblePA(const FiniteElementSpace &trial_fes,
1927  const FiniteElementSpace &test_fes);
1928 
1929  virtual void AddMultPA(const Vector&, Vector&) const;
1930  virtual void AddMultTransposePA(const Vector&, Vector&) const;
1931 
1932 private:
1933  // PA extension
1934  Vector pa_data;
1935  const DofToQuad *mapsO; ///< Not owned. DOF-to-quad map, open.
1936  const DofToQuad *mapsC; ///< Not owned. DOF-to-quad map, closed.
1937  const DofToQuad *mapsOtest; ///< Not owned. DOF-to-quad map, open.
1938  const DofToQuad *mapsCtest; ///< Not owned. DOF-to-quad map, closed.
1939  const GeometricFactors *geom; ///< Not owned
1940  int dim, ne, dofs1D, dofs1Dtest,quad1D, testType, trialType, coeffDim;
1941 };
1942 
1943 /** Class for integrating the bilinear form a(u,v) := (Q u, curl v) in 3D and
1944  where Q is an optional coefficient (of type scalar, matrix, or diagonal
1945  matrix) u is in H(Div) or H(Curl) and v is in H(Curl). */
1947 {
1948 public:
1951  : MixedVectorIntegrator(q) {}
1953  : MixedVectorIntegrator(dq, true) {}
1955  : MixedVectorIntegrator(mq) {}
1956 
1957 protected:
1958  inline virtual bool VerifyFiniteElementTypes(
1959  const FiniteElement & trial_fe,
1960  const FiniteElement & test_fe) const
1961  {
1962  return (trial_fe.GetVDim() == 3 && test_fe.GetCurlDim() == 3 &&
1963  trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
1964  test_fe.GetDerivType() == mfem::FiniteElement::CURL );
1965  }
1966 
1967  inline virtual const char * FiniteElementTypeFailureMessage() const
1968  {
1969  return "MixedVectorWeakCurlIntegrator: "
1970  "Trial space must be vector field in 3D and the "
1971  "test space must be H(Curl)";
1972  }
1973 
1974  inline virtual int GetTestVDim(const FiniteElement & test_fe)
1975  { return test_fe.GetCurlDim(); }
1976 
1977  inline virtual void CalcTestShape(const FiniteElement & test_fe,
1978  ElementTransformation &Trans,
1979  DenseMatrix & shape)
1980  {
1981  test_fe.CalcPhysCurlShape(Trans, shape);
1982  }
1983 
1985  virtual void AssemblePA(const FiniteElementSpace &trial_fes,
1986  const FiniteElementSpace &test_fes);
1987 
1988  virtual void AddMultPA(const Vector&, Vector&) const;
1989  virtual void AddMultTransposePA(const Vector&, Vector&) const;
1990 
1991 private:
1992  // PA extension
1993  Vector pa_data;
1994  const DofToQuad *mapsO; ///< Not owned. DOF-to-quad map, open.
1995  const DofToQuad *mapsC; ///< Not owned. DOF-to-quad map, closed.
1996  const GeometricFactors *geom; ///< Not owned
1997  int dim, ne, dofs1D, quad1D, testType, trialType, coeffDim;
1998 };
1999 
2000 /** Class for integrating the bilinear form a(u,v) := - (Q u, grad v) in either
2001  2D or 3D and where Q is an optional coefficient (of type scalar, matrix, or
2002  diagonal matrix) u is in H(Div) or H(Curl) and v is in H1. */
2004 {
2005 public:
2008  : MixedVectorIntegrator(q) {}
2010  : MixedVectorIntegrator(dq, true) {}
2012  : MixedVectorIntegrator(mq) {}
2013 
2014 protected:
2015  inline virtual bool VerifyFiniteElementTypes(
2016  const FiniteElement & trial_fe,
2017  const FiniteElement & test_fe) const
2018  {
2019  return (trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
2020  test_fe.GetDerivType() == mfem::FiniteElement::GRAD );
2021  }
2022 
2023  inline virtual const char * FiniteElementTypeFailureMessage() const
2024  {
2025  return "MixedVectorWeakDivergenceIntegrator: "
2026  "Trial space must be vector field and the "
2027  "test space must be H1";
2028  }
2029 
2030  inline virtual int GetTestVDim(const FiniteElement & test_fe)
2031  { return space_dim; }
2032 
2033  inline virtual void CalcTestShape(const FiniteElement & test_fe,
2034  ElementTransformation &Trans,
2035  DenseMatrix & shape)
2036  {
2037  test_fe.CalcPhysDShape(Trans, shape);
2038  shape *= -1.0;
2039  }
2040 };
2041 
2042 /** Class for integrating the bilinear form a(u,v) := (Q grad u, v) where Q is a
2043  scalar coefficient, and v is a vector with components v_i in the same (H1) space
2044  as u.
2045 
2046  See also MixedVectorGradientIntegrator when v is in H(curl). */
2048 {
2049 protected:
2051 
2052 private:
2053  Vector shape;
2054  DenseMatrix dshape;
2055  DenseMatrix gshape;
2056  DenseMatrix Jadj;
2057  DenseMatrix elmat_comp;
2058  // PA extension
2059  Vector pa_data;
2060  const DofToQuad *trial_maps, *test_maps; ///< Not owned
2061  const GeometricFactors *geom; ///< Not owned
2062  int dim, ne, nq;
2063  int trial_dofs1D, test_dofs1D, quad1D;
2064 
2065 public:
2067  Q{NULL}, trial_maps{NULL}, test_maps{NULL}, geom{NULL}
2068  { }
2070  Q{q_}, trial_maps{NULL}, test_maps{NULL}, geom{NULL}
2071  { }
2073  Q{&q}, trial_maps{NULL}, test_maps{NULL}, geom{NULL}
2074  { }
2075 
2076  virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
2077  const FiniteElement &test_fe,
2078  ElementTransformation &Trans,
2079  DenseMatrix &elmat);
2080 
2082  virtual void AssemblePA(const FiniteElementSpace &trial_fes,
2083  const FiniteElementSpace &test_fes);
2084 
2085  virtual void AddMultPA(const Vector &x, Vector &y) const;
2086  virtual void AddMultTransposePA(const Vector &x, Vector &y) const;
2087 
2088  static const IntegrationRule &GetRule(const FiniteElement &trial_fe,
2089  const FiniteElement &test_fe,
2090  ElementTransformation &Trans);
2091 };
2092 
2093 /** Class for integrating the bilinear form a(u,v) := (Q grad u, grad v) where Q
2094  can be a scalar or a matrix coefficient. */
2096 {
2097 protected:
2101 
2102 private:
2103  Vector vec, vecdxt, pointflux, shape;
2104 #ifndef MFEM_THREAD_SAFE
2105  DenseMatrix dshape, dshapedxt, invdfdx, M, dshapedxt_m;
2106  DenseMatrix te_dshape, te_dshapedxt;
2107  Vector D;
2108 #endif
2109 
2110  // PA extension
2111  const FiniteElementSpace *fespace;
2112  const DofToQuad *maps; ///< Not owned
2113  const GeometricFactors *geom; ///< Not owned
2114  int dim, ne, dofs1D, quad1D;
2115  Vector pa_data;
2116  bool symmetric = true; ///< False if using a nonsymmetric matrix coefficient
2117 
2118  // Data for NURBS patch PA
2119 
2120  // Type for a variable-row-length 2D array, used for data related to 1D
2121  // quadrature rules in each dimension.
2122  typedef std::vector<std::vector<int>> IntArrayVar2D;
2123 
2124  int numPatches = 0;
2125  static constexpr int numTypes = 2; // Number of rule types
2126 
2127  // In the case integrationMode == Mode::PATCHWISE_REDUCED, an approximate
2128  // integration rule with sparse nonzero weights is computed by NNLSSolver,
2129  // for each 1D basis function on each patch, in each spatial dimension. For a
2130  // fixed 1D basis function b_i with DOF index i, in the tensor product basis
2131  // of patch p, the prescribed exact 1D rule is of the form
2132  // \sum_k a_{i,j,k} w_k for some integration points indexed by k, with
2133  // weights w_k and coefficients a_{i,j,k} depending on Q(x), an element
2134  // transformation, b_i, and b_j, for all 1D basis functions b_j whose support
2135  // overlaps that of b_i. Define the constraint matrix G = [g_{j,k}] with
2136  // g_{j,k} = a_{i,j,k} and the vector of exact weights w = [w_k]. A reduced
2137  // rule should have different weights w_r, many of them zero, and should
2138  // approximately satisfy Gw_r = Gw. A sparse approximate solution to this
2139  // underdetermined system is computed by NNLSSolver, and its data is stored
2140  // in the following members.
2141 
2142  // For each patch p, spatial dimension d (total dim), and rule type t (total
2143  // numTypes), an std::vector<Vector> of reduced quadrature weights for all
2144  // basis functions is stored in reducedWeights[t + numTypes * (d + dim * p)],
2145  // reshaped as rw(t,d,p). Note that nd may vary with respect to the patch and
2146  // spatial dimension. Array reducedIDs is treated similarly.
2147  std::vector<std::vector<Vector>> reducedWeights;
2148  std::vector<IntArrayVar2D> reducedIDs;
2149  std::vector<Array<int>> pQ1D, pD1D;
2150  std::vector<std::vector<Array2D<double>>> pB, pG;
2151  std::vector<IntArrayVar2D> pminD, pmaxD, pminQ, pmaxQ, pminDD, pmaxDD;
2152 
2153  std::vector<Array<const IntegrationRule*>> pir1d;
2154 
2155  void SetupPatchPA(const int patch, Mesh *mesh, bool unitWeights=false);
2156 
2157  void SetupPatchBasisData(Mesh *mesh, unsigned int patch);
2158 
2159  /** Called by AssemblePatchMatrix for sparse matrix assembly on a NURBS patch
2160  with full 1D quadrature rules. */
2161  void AssemblePatchMatrix_fullQuadrature(const int patch,
2162  const FiniteElementSpace &fes,
2163  SparseMatrix*& smat);
2164 
2165  /** Called by AssemblePatchMatrix for sparse matrix assembly on a NURBS patch
2166  with reduced 1D quadrature rules. */
2167  void AssemblePatchMatrix_reducedQuadrature(const int patch,
2168  const FiniteElementSpace &fes,
2169  SparseMatrix*& smat);
2170 
2171 public:
2172  /// Construct a diffusion integrator with coefficient Q = 1
2174  : BilinearFormIntegrator(ir),
2175  Q(NULL), VQ(NULL), MQ(NULL), maps(NULL), geom(NULL) { }
2176 
2177  /// Construct a diffusion integrator with a scalar coefficient q
2179  : BilinearFormIntegrator(ir),
2180  Q(&q), VQ(NULL), MQ(NULL), maps(NULL), geom(NULL) { }
2181 
2182  /// Construct a diffusion integrator with a vector coefficient q
2184  const IntegrationRule *ir = nullptr)
2185  : BilinearFormIntegrator(ir),
2186  Q(NULL), VQ(&q), MQ(NULL), maps(NULL), geom(NULL) { }
2187 
2188  /// Construct a diffusion integrator with a matrix coefficient q
2190  const IntegrationRule *ir = nullptr)
2191  : BilinearFormIntegrator(ir),
2192  Q(NULL), VQ(NULL), MQ(&q), maps(NULL), geom(NULL) { }
2193 
2194  /** Given a particular Finite Element computes the element stiffness matrix
2195  elmat. */
2196  virtual void AssembleElementMatrix(const FiniteElement &el,
2197  ElementTransformation &Trans,
2198  DenseMatrix &elmat);
2199  /** Given a trial and test Finite Element computes the element stiffness
2200  matrix elmat. */
2201  virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
2202  const FiniteElement &test_fe,
2203  ElementTransformation &Trans,
2204  DenseMatrix &elmat);
2205 
2206  virtual void AssemblePatchMatrix(const int patch,
2207  const FiniteElementSpace &fes,
2208  SparseMatrix*& smat);
2209 
2210  virtual void AssembleNURBSPA(const FiniteElementSpace &fes);
2211 
2212  void AssemblePatchPA(const int patch, const FiniteElementSpace &fes);
2213 
2214  /// Perform the local action of the BilinearFormIntegrator
2215  virtual void AssembleElementVector(const FiniteElement &el,
2217  const Vector &elfun, Vector &elvect);
2218 
2219  virtual void ComputeElementFlux(const FiniteElement &el,
2220  ElementTransformation &Trans,
2221  Vector &u, const FiniteElement &fluxelem,
2222  Vector &flux, bool with_coef = true,
2223  const IntegrationRule *ir = NULL);
2224 
2225  virtual double ComputeFluxEnergy(const FiniteElement &fluxelem,
2226  ElementTransformation &Trans,
2227  Vector &flux, Vector *d_energy = NULL);
2228 
2230 
2231  virtual void AssembleMF(const FiniteElementSpace &fes);
2232 
2233  virtual void AssemblePA(const FiniteElementSpace &fes);
2234 
2235  virtual void AssembleEA(const FiniteElementSpace &fes, Vector &emat,
2236  const bool add);
2237 
2238  virtual void AssembleDiagonalPA(Vector &diag);
2239 
2240  virtual void AssembleDiagonalMF(Vector &diag);
2241 
2242  virtual void AddMultMF(const Vector&, Vector&) const;
2243 
2244  virtual void AddMultPA(const Vector&, Vector&) const;
2245 
2246  virtual void AddMultTransposePA(const Vector&, Vector&) const;
2247 
2248  virtual void AddMultNURBSPA(const Vector&, Vector&) const;
2249 
2250  void AddMultPatchPA(const int patch, const Vector &x, Vector &y) const;
2251 
2252  static const IntegrationRule &GetRule(const FiniteElement &trial_fe,
2253  const FiniteElement &test_fe);
2254 
2255  bool SupportsCeed() const { return DeviceCanUseCeed(); }
2256 
2257  Coefficient *GetCoefficient() const { return Q; }
2258 };
2259 
2260 /** Class for local mass matrix assembling a(u,v) := (Q u, v) */
2262 {
2263  friend class DGMassInverse;
2264 protected:
2265 #ifndef MFEM_THREAD_SAFE
2267 #endif
2269  // PA extension
2272  const DofToQuad *maps; ///< Not owned
2273  const GeometricFactors *geom; ///< Not owned
2274  const FaceGeometricFactors *face_geom; ///< Not owned
2275  int dim, ne, nq, dofs1D, quad1D;
2276 
2277 public:
2279  : BilinearFormIntegrator(ir), Q(NULL), maps(NULL), geom(NULL) { }
2280 
2281  /// Construct a mass integrator with coefficient q
2283  : BilinearFormIntegrator(ir), Q(&q), maps(NULL), geom(NULL) { }
2284 
2285  /** Given a particular Finite Element computes the element mass matrix
2286  elmat. */
2287  virtual void AssembleElementMatrix(const FiniteElement &el,
2288  ElementTransformation &Trans,
2289  DenseMatrix &elmat);
2290  virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
2291  const FiniteElement &test_fe,
2292  ElementTransformation &Trans,
2293  DenseMatrix &elmat);
2294 
2296 
2297  virtual void AssembleMF(const FiniteElementSpace &fes);
2298 
2299  virtual void AssemblePA(const FiniteElementSpace &fes);
2300 
2301  virtual void AssemblePABoundary(const FiniteElementSpace &fes);
2302 
2303  virtual void AssembleEA(const FiniteElementSpace &fes, Vector &emat,
2304  const bool add);
2305 
2306  virtual void AssembleDiagonalPA(Vector &diag);
2307 
2308  virtual void AssembleDiagonalMF(Vector &diag);
2309 
2310  virtual void AddMultMF(const Vector&, Vector&) const;
2311 
2312  virtual void AddMultPA(const Vector&, Vector&) const;
2313 
2314  virtual void AddMultTransposePA(const Vector&, Vector&) const;
2315 
2316  static const IntegrationRule &GetRule(const FiniteElement &trial_fe,
2317  const FiniteElement &test_fe,
2318  ElementTransformation &Trans);
2319 
2320  bool SupportsCeed() const { return DeviceCanUseCeed(); }
2321 
2322  const Coefficient *GetCoefficient() const { return Q; }
2323 };
2324 
2325 /** Mass integrator (u, v) restricted to the boundary of a domain */
2327 {
2328 public:
2330 
2332 
2333  virtual void AssembleFaceMatrix(const FiniteElement &el1,
2334  const FiniteElement &el2,
2336  DenseMatrix &elmat);
2337 };
2338 
2339 /// alpha (q . grad u, v)
2341 {
2342 protected:
2344  double alpha;
2345  // PA extension
2347  const DofToQuad *maps; ///< Not owned
2348  const GeometricFactors *geom; ///< Not owned
2349  int dim, ne, nq, dofs1D, quad1D;
2350 
2351 private:
2352 #ifndef MFEM_THREAD_SAFE
2353  DenseMatrix dshape, adjJ, Q_ir;
2354  Vector shape, vec2, BdFidxT;
2355 #endif
2356 
2357 public:
2359  : Q(&q) { alpha = a; }
2360 
2361  virtual void AssembleElementMatrix(const FiniteElement &,
2363  DenseMatrix &);
2364 
2366 
2367  virtual void AssembleMF(const FiniteElementSpace &fes);
2368 
2369  virtual void AssemblePA(const FiniteElementSpace&);
2370 
2371  virtual void AssembleEA(const FiniteElementSpace &fes, Vector &emat,
2372  const bool add);
2373 
2374  virtual void AssembleDiagonalPA(Vector &diag);
2375 
2376  virtual void AssembleDiagonalMF(Vector &diag);
2377 
2378  virtual void AddMultMF(const Vector&, Vector&) const;
2379 
2380  virtual void AddMultPA(const Vector&, Vector&) const;
2381 
2382  virtual void AddMultTransposePA(const Vector &x, Vector &y) const;
2383 
2384  static const IntegrationRule &GetRule(const FiniteElement &el,
2385  ElementTransformation &Trans);
2386 
2387  static const IntegrationRule &GetRule(const FiniteElement &trial_fe,
2388  const FiniteElement &test_fe,
2389  ElementTransformation &Trans);
2390 
2391  bool SupportsCeed() const { return DeviceCanUseCeed(); }
2392 };
2393 
2394 // Alias for @ConvectionIntegrator.
2396 
2397 /// -alpha (u, q . grad v), negative transpose of ConvectionIntegrator
2399 {
2400 public:
2403 };
2404 
2405 /// alpha (q . grad u, v) using the "group" FE discretization
2407 {
2408 protected:
2410  double alpha;
2411 
2412 private:
2413  DenseMatrix dshape, adjJ, Q_nodal, grad;
2414  Vector shape;
2415 
2416 public:
2418  : Q(&q) { alpha = a; }
2419  virtual void AssembleElementMatrix(const FiniteElement &,
2421  DenseMatrix &);
2422 };
2423 
2424 /** Class for integrating the bilinear form a(u,v) := (Q u, v),
2425  where u=(u1,...,un) and v=(v1,...,vn); ui and vi are defined
2426  by scalar FE through standard transformation. */
2428 {
2429 private:
2430  int vdim;
2431  Vector shape, te_shape, vec;
2432  DenseMatrix partelmat;
2433  DenseMatrix mcoeff;
2434  int Q_order;
2435 
2436 protected:
2440  // PA extension
2442  const DofToQuad *maps; ///< Not owned
2443  const GeometricFactors *geom; ///< Not owned
2444  int dim, ne, nq, dofs1D, quad1D;
2445 
2446 public:
2447  /// Construct an integrator with coefficient 1.0
2449  : vdim(-1), Q_order(0), Q(NULL), VQ(NULL), MQ(NULL) { }
2450  /** Construct an integrator with scalar coefficient q. If possible, save
2451  memory by using a scalar integrator since the resulting matrix is block
2452  diagonal with the same diagonal block repeated. */
2454  : vdim(-1), Q_order(qo), Q(&q), VQ(NULL), MQ(NULL) { }
2456  : BilinearFormIntegrator(ir), vdim(-1), Q_order(0), Q(&q), VQ(NULL),
2457  MQ(NULL) { }
2458  /// Construct an integrator with diagonal coefficient q
2460  : vdim(q.GetVDim()), Q_order(qo), Q(NULL), VQ(&q), MQ(NULL) { }
2461  /// Construct an integrator with matrix coefficient q
2463  : vdim(q.GetVDim()), Q_order(qo), Q(NULL), VQ(NULL), MQ(&q) { }
2464 
2465  int GetVDim() const { return vdim; }
2466  void SetVDim(int vdim_) { vdim = vdim_; }
2467 
2468  virtual void AssembleElementMatrix(const FiniteElement &el,
2469  ElementTransformation &Trans,
2470  DenseMatrix &elmat);
2471  virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
2472  const FiniteElement &test_fe,
2473  ElementTransformation &Trans,
2474  DenseMatrix &elmat);
2476  virtual void AssemblePA(const FiniteElementSpace &fes);
2477  virtual void AssembleMF(const FiniteElementSpace &fes);
2478  virtual void AssembleDiagonalPA(Vector &diag);
2479  virtual void AssembleDiagonalMF(Vector &diag);
2480  virtual void AddMultPA(const Vector &x, Vector &y) const;
2481  virtual void AddMultMF(const Vector &x, Vector &y) const;
2482  bool SupportsCeed() const { return DeviceCanUseCeed(); }
2483 };
2484 
2485 
2486 /** Class for integrating (div u, p) where u is a vector field given by
2487  VectorFiniteElement through Piola transformation (for RT elements); p is
2488  scalar function given by FiniteElement through standard transformation.
2489  Here, u is the trial function and p is the test function.
2490 
2491  Note: if the test space does not have map type INTEGRAL, then the element
2492  matrix returned by AssembleElementMatrix2 will not depend on the
2493  ElementTransformation Trans. */
2495 {
2496 protected:
2498 
2500  virtual void AssemblePA(const FiniteElementSpace &trial_fes,
2501  const FiniteElementSpace &test_fes);
2502 
2503  virtual void AddMultPA(const Vector&, Vector&) const;
2504  virtual void AddMultTransposePA(const Vector&, Vector&) const;
2505 
2506 private:
2507 #ifndef MFEM_THREAD_SAFE
2508  Vector divshape, shape;
2509 #endif
2510 
2511  // PA extension
2512  Vector pa_data;
2513  const DofToQuad *mapsO; ///< Not owned. DOF-to-quad map, open.
2514  const DofToQuad *L2mapsO; ///< Not owned. DOF-to-quad map, open.
2515  const DofToQuad *mapsC; ///< Not owned. DOF-to-quad map, closed.
2516  int dim, ne, dofs1D, L2dofs1D, quad1D;
2517 
2518 public:
2521  virtual void AssembleElementMatrix(const FiniteElement &el,
2522  ElementTransformation &Trans,
2523  DenseMatrix &elmat) { }
2524  virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
2525  const FiniteElement &test_fe,
2526  ElementTransformation &Trans,
2527  DenseMatrix &elmat);
2528 
2529  virtual void AssembleDiagonalPA_ADAt(const Vector &D, Vector &diag);
2530 };
2531 
2532 
2533 /** Integrator for `(-Q u, grad v)` for Nedelec (`u`) and H1 (`v`) elements.
2534  This is equivalent to a weak divergence of the Nedelec basis functions. */
2536 {
2537 protected:
2539 
2540 private:
2541 #ifndef MFEM_THREAD_SAFE
2542  DenseMatrix dshape;
2543  DenseMatrix dshapedxt;
2544  DenseMatrix vshape;
2545  DenseMatrix invdfdx;
2546 #endif
2547 
2548 public:
2551  virtual void AssembleElementMatrix(const FiniteElement &el,
2552  ElementTransformation &Trans,
2553  DenseMatrix &elmat) { }
2554  virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
2555  const FiniteElement &test_fe,
2556  ElementTransformation &Trans,
2557  DenseMatrix &elmat);
2558 };
2559 
2560 /** Integrator for (curl u, v) for Nedelec and RT elements. If the trial and
2561  test spaces are switched, assembles the form (u, curl v). */
2563 {
2564 protected:
2566 
2567 private:
2568 #ifndef MFEM_THREAD_SAFE
2569  DenseMatrix curlshapeTrial;
2570  DenseMatrix vshapeTest;
2571  DenseMatrix curlshapeTrial_dFT;
2572 #endif
2573 
2574 public:
2577  virtual void AssembleElementMatrix(const FiniteElement &el,
2578  ElementTransformation &Trans,
2579  DenseMatrix &elmat) { }
2580  virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
2581  const FiniteElement &test_fe,
2582  ElementTransformation &Trans,
2583  DenseMatrix &elmat);
2584 };
2585 
2586 /// Class for integrating (Q D_i(u), v); u and v are scalars
2588 {
2589 protected:
2591 
2592 private:
2593  int xi;
2594  DenseMatrix dshape, dshapedxt, invdfdx;
2595  Vector shape, dshapedxi;
2596 
2597 public:
2598  DerivativeIntegrator(Coefficient &q, int i) : Q(&q), xi(i) { }
2599  virtual void AssembleElementMatrix(const FiniteElement &el,
2600  ElementTransformation &Trans,
2601  DenseMatrix &elmat)
2602  { AssembleElementMatrix2(el,el,Trans,elmat); }
2603  virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
2604  const FiniteElement &test_fe,
2605  ElementTransformation &Trans,
2606  DenseMatrix &elmat);
2607 };
2608 
2609 /// Integrator for (curl u, curl v) for Nedelec elements
2611 {
2612 private:
2613  Vector vec, pointflux;
2614 #ifndef MFEM_THREAD_SAFE
2615  Vector D;
2616  DenseMatrix curlshape, curlshape_dFt, M;
2617  DenseMatrix te_curlshape, te_curlshape_dFt;
2618  DenseMatrix vshape, projcurl;
2619 #endif
2620 
2621 protected:
2625 
2626  // PA extension
2628  const DofToQuad *mapsO; ///< Not owned. DOF-to-quad map, open.
2629  const DofToQuad *mapsC; ///< Not owned. DOF-to-quad map, closed.
2630  const GeometricFactors *geom; ///< Not owned
2631  int dim, ne, nq, dofs1D, quad1D;
2632  bool symmetric = true; ///< False if using a nonsymmetric matrix coefficient
2633 
2634 public:
2635  CurlCurlIntegrator() { Q = NULL; DQ = NULL; MQ = NULL; }
2636  /// Construct a bilinear form integrator for Nedelec elements
2638  BilinearFormIntegrator(ir), Q(&q), DQ(NULL), MQ(NULL) { }
2640  const IntegrationRule *ir = NULL) :
2641  BilinearFormIntegrator(ir), Q(NULL), DQ(&dq), MQ(NULL) { }
2643  BilinearFormIntegrator(ir), Q(NULL), DQ(NULL), MQ(&mq) { }
2644 
2645  /* Given a particular Finite Element, compute the
2646  element curl-curl matrix elmat */
2647  virtual void AssembleElementMatrix(const FiniteElement &el,
2648  ElementTransformation &Trans,
2649  DenseMatrix &elmat);
2650 
2651  virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
2652  const FiniteElement &test_fe,
2653  ElementTransformation &Trans,
2654  DenseMatrix &elmat);
2655 
2656  virtual void ComputeElementFlux(const FiniteElement &el,
2657  ElementTransformation &Trans,
2658  Vector &u, const FiniteElement &fluxelem,
2659  Vector &flux, bool with_coef,
2660  const IntegrationRule *ir = NULL);
2661 
2662  virtual double ComputeFluxEnergy(const FiniteElement &fluxelem,
2663  ElementTransformation &Trans,
2664  Vector &flux, Vector *d_energy = NULL);
2665 
2667  virtual void AssemblePA(const FiniteElementSpace &fes);
2668  virtual void AddMultPA(const Vector &x, Vector &y) const;
2669  virtual void AssembleDiagonalPA(Vector& diag);
2670 
2671  const Coefficient *GetCoefficient() const { return Q; }
2672 };
2673 
2674 /** Integrator for (curl u, curl v) for FE spaces defined by 'dim' copies of a
2675  scalar FE space. */
2677 {
2678 private:
2679 #ifndef MFEM_THREAD_SAFE
2680  DenseMatrix dshape_hat, dshape, curlshape, Jadj, grad_hat, grad;
2681 #endif
2682 
2683 protected:
2685 
2686 public:
2688 
2690 
2691  /// Assemble an element matrix
2692  virtual void AssembleElementMatrix(const FiniteElement &el,
2693  ElementTransformation &Trans,
2694  DenseMatrix &elmat);
2695  /// Compute element energy: (1/2) (curl u, curl u)_E
2696  virtual double GetElementEnergy(const FiniteElement &el,
2698  const Vector &elfun);
2699 };
2700 
2701 /** Class for integrating the bilinear form a(u,v) := (Q curl u, v) where Q is
2702  an optional scalar coefficient, and v is a vector with components v_i in
2703  the L2 or H1 space. This integrator handles 3 cases:
2704  (a) u āˆˆ H(curl) in 3D, v is a 3D vector with components v_i in L^2 or H^1
2705  (b) u āˆˆ H(curl) in 2D, v is a scalar field in L^2 or H^1
2706  (c) u is a scalar field in H^1, i.e, curl u := [0 1;-1 0]grad u and v is a
2707  2D vector field with components v_i in L^2 or H^1 space.
2708  Note: Case (b) can also be handled by MixedScalarCurlIntegrator */
2710 {
2711 protected:
2713 
2714 private:
2715  Vector shape;
2716  DenseMatrix dshape;
2717  DenseMatrix curlshape;
2718  DenseMatrix elmat_comp;
2719 public:
2720  MixedCurlIntegrator() : Q{NULL} { }
2723 
2724  virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
2725  const FiniteElement &test_fe,
2726  ElementTransformation &Trans,
2727  DenseMatrix &elmat);
2728 };
2729 
2730 /** Integrator for (Q u, v), where Q is an optional coefficient (of type scalar,
2731  vector (diagonal matrix), or matrix), trial function u is in H(Curl) or
2732  H(Div), and test function v is in H(Curl), H(Div), or v=(v1,...,vn), where
2733  vi are in H1. */
2735 {
2736 private:
2738  { Q = q; DQ = dq; MQ = mq; }
2739 
2740 #ifndef MFEM_THREAD_SAFE
2741  Vector shape;
2742  Vector D;
2743  DenseMatrix K;
2744  DenseMatrix partelmat;
2745  DenseMatrix test_vshape;
2746  DenseMatrix trial_vshape;
2747 #endif
2748 
2749 protected:
2753 
2754  // PA extension
2756  const DofToQuad *mapsO; ///< Not owned. DOF-to-quad map, open.
2757  const DofToQuad *mapsC; ///< Not owned. DOF-to-quad map, closed.
2758  const DofToQuad *mapsOtest; ///< Not owned. DOF-to-quad map, open.
2759  const DofToQuad *mapsCtest; ///< Not owned. DOF-to-quad map, closed.
2760  const GeometricFactors *geom; ///< Not owned
2762  bool symmetric = true; ///< False if using a nonsymmetric matrix coefficient
2763 
2764 public:
2765  VectorFEMassIntegrator() { Init(NULL, NULL, NULL); }
2766  VectorFEMassIntegrator(Coefficient *q_) { Init(q_, NULL, NULL); }
2767  VectorFEMassIntegrator(Coefficient &q) { Init(&q, NULL, NULL); }
2770  VectorFEMassIntegrator(MatrixCoefficient *mq_) { Init(NULL, NULL, mq_); }
2771  VectorFEMassIntegrator(MatrixCoefficient &mq) { Init(NULL, NULL, &mq); }
2772 
2773  virtual void AssembleElementMatrix(const FiniteElement &el,
2774  ElementTransformation &Trans,
2775  DenseMatrix &elmat);
2776  virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
2777  const FiniteElement &test_fe,
2778  ElementTransformation &Trans,
2779  DenseMatrix &elmat);
2780 
2782  virtual void AssemblePA(const FiniteElementSpace &fes);
2783  virtual void AssemblePA(const FiniteElementSpace &trial_fes,
2784  const FiniteElementSpace &test_fes);
2785  virtual void AddMultPA(const Vector &x, Vector &y) const;
2786  virtual void AddMultTransposePA(const Vector &x, Vector &y) const;
2787  virtual void AssembleDiagonalPA(Vector& diag);
2788 
2789  const Coefficient *GetCoefficient() const { return Q; }
2790 };
2791 
2792 /** Integrator for (Q div u, p) where u=(v1,...,vn) and all vi are in the same
2793  scalar FE space; p is also in a (different) scalar FE space. */
2795 {
2796 protected:
2798 
2799 private:
2800  Vector shape;
2801  Vector divshape;
2802  DenseMatrix dshape;
2803  DenseMatrix gshape;
2804  DenseMatrix Jadj;
2805  // PA extension
2806  Vector pa_data;
2807  const DofToQuad *trial_maps, *test_maps; ///< Not owned
2808  const GeometricFactors *geom; ///< Not owned
2809  int dim, ne, nq;
2810  int trial_dofs1D, test_dofs1D, quad1D;
2811 
2812 public:
2814  Q(NULL), trial_maps(NULL), test_maps(NULL), geom(NULL)
2815  { }
2817  Q(q_), trial_maps(NULL), test_maps(NULL), geom(NULL)
2818  { }
2820  Q(&q), trial_maps(NULL), test_maps(NULL), geom(NULL)
2821  { }
2822 
2823  virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
2824  const FiniteElement &test_fe,
2825  ElementTransformation &Trans,
2826  DenseMatrix &elmat);
2827 
2829  virtual void AssemblePA(const FiniteElementSpace &trial_fes,
2830  const FiniteElementSpace &test_fes);
2831 
2832  virtual void AddMultPA(const Vector &x, Vector &y) const;
2833  virtual void AddMultTransposePA(const Vector &x, Vector &y) const;
2834 
2835  static const IntegrationRule &GetRule(const FiniteElement &trial_fe,
2836  const FiniteElement &test_fe,
2837  ElementTransformation &Trans);
2838 };
2839 
2840 /// (Q div u, div v) for RT elements
2842 {
2843 protected:
2845 
2847  virtual void AssemblePA(const FiniteElementSpace &fes);
2848  virtual void AddMultPA(const Vector &x, Vector &y) const;
2849  virtual void AssembleDiagonalPA(Vector& diag);
2850 
2851 private:
2852 #ifndef MFEM_THREAD_SAFE
2853  Vector divshape, te_divshape;
2854 #endif
2855 
2856  // PA extension
2857  Vector pa_data;
2858  const DofToQuad *mapsO; ///< Not owned. DOF-to-quad map, open.
2859  const DofToQuad *mapsC; ///< Not owned. DOF-to-quad map, closed.
2860  const GeometricFactors *geom; ///< Not owned
2861  int dim, ne, dofs1D, quad1D;
2862 
2863 public:
2864  DivDivIntegrator() { Q = NULL; }
2866  BilinearFormIntegrator(ir), Q(&q) { }
2867 
2868  virtual void AssembleElementMatrix(const FiniteElement &el,
2869  ElementTransformation &Trans,
2870  DenseMatrix &elmat);
2871 
2872  virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
2873  const FiniteElement &test_fe,
2874  ElementTransformation &Trans,
2875  DenseMatrix &elmat);
2876 
2877  const Coefficient *GetCoefficient() const { return Q; }
2878 };
2879 
2880 /** Integrator for
2881 
2882  (Q grad u, grad v) = sum_i (Q grad u_i, grad v_i) e_i e_i^T
2883 
2884  for vector FE spaces, where e_i is the unit vector in the i-th direction.
2885  The resulting local element matrix is square, of size <tt> vdim*dof </tt>,
2886  where \c vdim is the vector dimension space and \c dof is the local degrees
2887  of freedom. The integrator is not aware of the true vector dimension and
2888  must use \c VectorCoefficient, \c MatrixCoefficient, or a caller-specified
2889  value to determine the vector space. For a scalar coefficient, the caller
2890  may manually specify the vector dimension or the vector dimension is assumed
2891  to be the spatial dimension (i.e. 2-dimension or 3-dimension).
2892 */
2894 {
2895 protected:
2896  Coefficient *Q = NULL;
2899 
2900  // PA extension
2901  const DofToQuad *maps; ///< Not owned
2902  const GeometricFactors *geom; ///< Not owned
2905 
2906 private:
2907  DenseMatrix dshape, dshapedxt, pelmat;
2908  int vdim = -1;
2909  DenseMatrix mcoeff;
2910  Vector vcoeff;
2911 
2912 public:
2914 
2915  /** \brief Integrator with unit coefficient for caller-specified vector
2916  dimension.
2917 
2918  If the vector dimension does not match the true dimension of the space,
2919  the resulting element matrix will be mathematically invalid. */
2920  VectorDiffusionIntegrator(int vector_dimension)
2921  : vdim(vector_dimension) { }
2922 
2924  : Q(&q) { }
2925 
2927  : BilinearFormIntegrator(ir), Q(&q) { }
2928 
2929  /** \brief Integrator with scalar coefficient for caller-specified vector
2930  dimension.
2931 
2932  The element matrix is block-diagonal with \c vdim copies of the element
2933  matrix integrated with the \c Coefficient.
2934 
2935  If the vector dimension does not match the true dimension of the space,
2936  the resulting element matrix will be mathematically invalid. */
2937  VectorDiffusionIntegrator(Coefficient &q, int vector_dimension)
2938  : Q(&q), vdim(vector_dimension) { }
2939 
2940  /** \brief Integrator with \c VectorCoefficient. The vector dimension of the
2941  \c FiniteElementSpace is assumed to be the same as the dimension of the
2942  \c Vector.
2943 
2944  The element matrix is block-diagonal and each block is integrated with
2945  coefficient q_i.
2946 
2947  If the vector dimension does not match the true dimension of the space,
2948  the resulting element matrix will be mathematically invalid. */
2950  : VQ(&vq), vdim(vq.GetVDim()) { }
2951 
2952  /** \brief Integrator with \c MatrixCoefficient. The vector dimension of the
2953  \c FiniteElementSpace is assumed to be the same as the dimension of the
2954  \c Matrix.
2955 
2956  The element matrix is populated in each block. Each block is integrated
2957  with coefficient q_ij.
2958 
2959  If the vector dimension does not match the true dimension of the space,
2960  the resulting element matrix will be mathematically invalid. */
2962  : MQ(&mq), vdim(mq.GetVDim()) { }
2963 
2964  virtual void AssembleElementMatrix(const FiniteElement &el,
2965  ElementTransformation &Trans,
2966  DenseMatrix &elmat);
2967  virtual void AssembleElementVector(const FiniteElement &el,
2969  const Vector &elfun, Vector &elvect);
2971  virtual void AssemblePA(const FiniteElementSpace &fes);
2972  virtual void AssembleMF(const FiniteElementSpace &fes);
2973  virtual void AssembleDiagonalPA(Vector &diag);
2974  virtual void AssembleDiagonalMF(Vector &diag);
2975  virtual void AddMultPA(const Vector &x, Vector &y) const;
2976  virtual void AddMultMF(const Vector &x, Vector &y) const;
2977  bool SupportsCeed() const { return DeviceCanUseCeed(); }
2978 };
2979 
2980 /** Integrator for the linear elasticity form:
2981  a(u,v) = (lambda div(u), div(v)) + (2 mu e(u), e(v)),
2982  where e(v) = (1/2) (grad(v) + grad(v)^T).
2983  This is a 'Vector' integrator, i.e. defined for FE spaces
2984  using multiple copies of a scalar FE space. */
2986 {
2987 protected:
2988  double q_lambda, q_mu;
2990 
2991 private:
2992 #ifndef MFEM_THREAD_SAFE
2993  Vector shape;
2994  DenseMatrix dshape, gshape, pelmat;
2995  Vector divshape;
2996 #endif
2997 
2998 public:
3000  { lambda = &l; mu = &m; }
3001  /** With this constructor lambda = q_l * m and mu = q_m * m;
3002  if dim * q_l + 2 * q_m = 0 then trace(sigma) = 0. */
3003  ElasticityIntegrator(Coefficient &m, double q_l, double q_m)
3004  { lambda = NULL; mu = &m; q_lambda = q_l; q_mu = q_m; }
3005 
3006  virtual void AssembleElementMatrix(const FiniteElement &,
3008  DenseMatrix &);
3009 
3010  /** Compute the stress corresponding to the local displacement @a u and
3011  interpolate it at the nodes of the given @a fluxelem. Only the symmetric
3012  part of the stress is stored, so that the size of @a flux is equal to
3013  the number of DOFs in @a fluxelem times dim*(dim+1)/2. In 2D, the order
3014  of the stress components is: s_xx, s_yy, s_xy. In 3D, it is: s_xx, s_yy,
3015  s_zz, s_xy, s_xz, s_yz. In other words, @a flux is the local vector for
3016  a FE space with dim*(dim+1)/2 vector components, based on the finite
3017  element @a fluxelem. The integration rule is taken from @a fluxelem.
3018  @a ir exists to specific an alternative integration rule. */
3019  virtual void ComputeElementFlux(const FiniteElement &el,
3020  ElementTransformation &Trans,
3021  Vector &u,
3022  const FiniteElement &fluxelem,
3023  Vector &flux, bool with_coef = true,
3024  const IntegrationRule *ir = NULL);
3025 
3026  /** Compute the element energy (integral of the strain energy density)
3027  corresponding to the stress represented by @a flux which is a vector of
3028  coefficients multiplying the basis functions defined by @a fluxelem. In
3029  other words, @a flux is the local vector for a FE space with
3030  dim*(dim+1)/2 vector components, based on the finite element @a fluxelem.
3031  The number of components, dim*(dim+1)/2 is such that it represents the
3032  symmetric part of the (symmetric) stress tensor. The order of the
3033  components is: s_xx, s_yy, s_xy in 2D, and s_xx, s_yy, s_zz, s_xy, s_xz,
3034  s_yz in 3D. */
3035  virtual double ComputeFluxEnergy(const FiniteElement &fluxelem,
3036  ElementTransformation &Trans,
3037  Vector &flux, Vector *d_energy = NULL);
3038 };
3039 
3040 /** Integrator for the DG form:
3041  alpha < rho_u (u.n) {v},[w] > + beta < rho_u |u.n| [v],[w] >,
3042  where v and w are the trial and test variables, respectively, and rho/u are
3043  given scalar/vector coefficients. {v} represents the average value of v on
3044  the face and [v] is the jump such that {v}=(v1+v2)/2 and [v]=(v1-v2) for the
3045  face between elements 1 and 2. For boundary elements, v2=0. The vector
3046  coefficient, u, is assumed to be continuous across the faces and when given
3047  the scalar coefficient, rho, is assumed to be discontinuous. The integrator
3048  uses the upwind value of rho, rho_u, which is value from the side into which
3049  the vector coefficient, u, points.
3050 
3051  One use case for this integrator is to discretize the operator -u.grad(v)
3052  with a DG formulation. The resulting formulation uses the
3053  ConvectionIntegrator (with coefficient u, and parameter alpha = -1) and the
3054  transpose of the DGTraceIntegrator (with coefficient u, and parameters alpha
3055  = 1, beta = -1/2 to use the upwind face flux, see also
3056  NonconservativeDGTraceIntegrator). This discretization and the handling of
3057  the inflow and outflow boundaries is illustrated in Example 9/9p.
3058 
3059  Another use case for this integrator is to discretize the operator -div(u v)
3060  with a DG formulation. The resulting formulation is conservative and
3061  consists of the ConservativeConvectionIntegrator (with coefficient u, and
3062  parameter alpha = -1) plus the DGTraceIntegrator (with coefficient u, and
3063  parameters alpha = -1, beta = -1/2 to use the upwind face flux).
3064  */
3066 {
3067 protected:
3070  double alpha, beta;
3071  // PA extension
3073  const DofToQuad *maps; ///< Not owned
3074  const FaceGeometricFactors *geom; ///< Not owned
3075  int dim, nf, nq, dofs1D, quad1D;
3076 
3077 private:
3078  Vector shape1, shape2;
3079 
3080 public:
3081  /// Construct integrator with rho = 1, b = 0.5*a.
3083  { rho = NULL; u = &u_; alpha = a; beta = 0.5*a; }
3084 
3085  /// Construct integrator with rho = 1.
3087  { rho = NULL; u = &u_; alpha = a; beta = b; }
3088 
3090  double a, double b)
3091  { rho = &rho_; u = &u_; alpha = a; beta = b; }
3092 
3094  virtual void AssembleFaceMatrix(const FiniteElement &el1,
3095  const FiniteElement &el2,
3097  DenseMatrix &elmat);
3098 
3100 
3101  virtual void AssemblePAInteriorFaces(const FiniteElementSpace &fes);
3102 
3103  virtual void AssemblePABoundaryFaces(const FiniteElementSpace &fes);
3104 
3105  virtual void AddMultTransposePA(const Vector &x, Vector &y) const;
3106 
3107  virtual void AddMultPA(const Vector&, Vector&) const;
3108 
3109  virtual void AssembleEAInteriorFaces(const FiniteElementSpace& fes,
3110  Vector &ea_data_int,
3111  Vector &ea_data_ext,
3112  const bool add);
3113 
3114  virtual void AssembleEABoundaryFaces(const FiniteElementSpace& fes,
3115  Vector &ea_data_bdr,
3116  const bool add);
3117 
3118  static const IntegrationRule &GetRule(Geometry::Type geom, int order,
3120 
3121 private:
3122  void SetupPA(const FiniteElementSpace &fes, FaceType type);
3123 };
3124 
3125 // Alias for @a DGTraceIntegrator.
3127 
3128 /** Integrator that represents the face terms used for the non-conservative
3129  DG discretization of the convection equation:
3130  -alpha < rho_u (u.n) {v},[w] > + beta < rho_u |u.n| [v],[w] >.
3131 
3132  This integrator can be used with together with ConvectionIntegrator to
3133  implement an upwind DG discretization in non-conservative form, see ex9 and
3134  ex9p. */
3136 {
3137 public:
3139  : TransposeIntegrator(new DGTraceIntegrator(u, -a, 0.5*a)) { }
3140 
3142  : TransposeIntegrator(new DGTraceIntegrator(u, -a, b)) { }
3143 
3145  double a, double b)
3146  : TransposeIntegrator(new DGTraceIntegrator(rho, u, -a, b)) { }
3147 };
3148 
3149 /** Integrator for the DG form:
3150 
3151  - < {(Q grad(u)).n}, [v] > + sigma < [u], {(Q grad(v)).n} >
3152  + kappa < {h^{-1} Q} [u], [v] >
3153 
3154  where Q is a scalar or matrix diffusion coefficient and u, v are the trial
3155  and test spaces, respectively. The parameters sigma and kappa determine the
3156  DG method to be used (when this integrator is added to the "broken"
3157  DiffusionIntegrator):
3158  * sigma = -1, kappa >= kappa0: symm. interior penalty (IP or SIPG) method,
3159  * sigma = +1, kappa > 0: non-symmetric interior penalty (NIPG) method,
3160  * sigma = +1, kappa = 0: the method of Baumann and Oden. */
3162 {
3163 protected:
3166  double sigma, kappa;
3167 
3168  // these are not thread-safe!
3171 
3172 public:
3173  DGDiffusionIntegrator(const double s, const double k)
3174  : Q(NULL), MQ(NULL), sigma(s), kappa(k) { }
3175  DGDiffusionIntegrator(Coefficient &q, const double s, const double k)
3176  : Q(&q), MQ(NULL), sigma(s), kappa(k) { }
3177  DGDiffusionIntegrator(MatrixCoefficient &q, const double s, const double k)
3178  : Q(NULL), MQ(&q), sigma(s), kappa(k) { }
3180  virtual void AssembleFaceMatrix(const FiniteElement &el1,
3181  const FiniteElement &el2,
3183  DenseMatrix &elmat);
3184 };
3185 
3186 /** Integrator for the "BR2" diffusion stabilization term
3187 
3188  sum_e eta (r_e([u]), r_e([v]))
3189 
3190  where r_e is the lifting operator defined on each edge e (potentially
3191  weighted by a coefficient Q). The parameter eta can be chosen to be one to
3192  obtain a stable discretization. The constructor for this integrator requires
3193  the finite element space because the lifting operator depends on the
3194  element-wise inverse mass matrix.
3195 
3196  BR2 stands for the second method of Bassi and Rebay:
3197 
3198  - F. Bassi and S. Rebay. A high order discontinuous Galerkin method for
3199  compressible turbulent flows. In B. Cockburn, G. E. Karniadakis, and
3200  C.-W. Shu, editors, Discontinuous Galerkin Methods, pages 77-88. Springer
3201  Berlin Heidelberg, 2000.
3202  - D. N. Arnold, F. Brezzi, B. Cockburn, and L. D. Marini. Unified analysis
3203  of discontinuous Galerkin methods for elliptic problems. SIAM Journal on
3204  Numerical Analysis, 39(5):1749-1779, 2002.
3205 */
3207 {
3208 protected:
3209  double eta;
3210 
3211  // Block factorizations of local mass matrices, with offsets for the case of
3212  // not equally sized blocks (mixed meshes, p-refinement)
3216 
3218 
3220 
3224 
3225  /// Precomputes the inverses (LU factorizations) of the local mass matrices.
3226  /** @a fes must be a DG space, so the mass matrix is block diagonal, and its
3227  inverse can be computed locally. This is required for the computation of
3228  the lifting operators @a r_e.
3229  */
3230  void PrecomputeMassInverse(class FiniteElementSpace &fes);
3231 
3232 public:
3233  DGDiffusionBR2Integrator(class FiniteElementSpace &fes, double e = 1.0);
3235  double e = 1.0);
3236  MFEM_DEPRECATED DGDiffusionBR2Integrator(class FiniteElementSpace *fes,
3237  double e = 1.0);
3238 
3240  virtual void AssembleFaceMatrix(const FiniteElement &el1,
3241  const FiniteElement &el2,
3243  DenseMatrix &elmat);
3244 };
3245 
3246 /** Integrator for the DG elasticity form, for the formulations see:
3247  - PhD Thesis of Jonas De Basabe, High-Order Finite %Element Methods for
3248  Seismic Wave Propagation, UT Austin, 2009, p. 23, and references therein
3249  - Peter Hansbo and Mats G. Larson, Discontinuous Galerkin and the
3250  Crouzeix-Raviart %Element: Application to Elasticity, PREPRINT 2000-09,
3251  p.3
3252 
3253  \f[
3254  - \left< \{ \tau(u) \}, [v] \right> + \alpha \left< \{ \tau(v) \}, [u]
3255  \right> + \kappa \left< h^{-1} \{ \lambda + 2 \mu \} [u], [v] \right>
3256  \f]
3257 
3258  where \f$ \left<u, v\right> = \int_{F} u \cdot v \f$, and \f$ F \f$ is a
3259  face which is either a boundary face \f$ F_b \f$ of an element \f$ K \f$ or
3260  an interior face \f$ F_i \f$ separating elements \f$ K_1 \f$ and \f$ K_2 \f$.
3261 
3262  In the bilinear form above \f$ \tau(u) \f$ is traction, and it's also
3263  \f$ \tau(u) = \sigma(u) \cdot \vec{n} \f$, where \f$ \sigma(u) \f$ is
3264  stress, and \f$ \vec{n} \f$ is the unit normal vector w.r.t. to \f$ F \f$.
3265 
3266  In other words, we have
3267  \f[
3268  - \left< \{ \sigma(u) \cdot \vec{n} \}, [v] \right> + \alpha \left< \{
3269  \sigma(v) \cdot \vec{n} \}, [u] \right> + \kappa \left< h^{-1} \{
3270  \lambda + 2 \mu \} [u], [v] \right>
3271  \f]
3272 
3273  For isotropic media
3274  \f[
3275  \begin{split}
3276  \sigma(u) &= \lambda \nabla \cdot u I + 2 \mu \varepsilon(u) \\
3277  &= \lambda \nabla \cdot u I + 2 \mu \frac{1}{2} (\nabla u + \nabla
3278  u^T) \\
3279  &= \lambda \nabla \cdot u I + \mu (\nabla u + \nabla u^T)
3280  \end{split}
3281  \f]
3282 
3283  where \f$ I \f$ is identity matrix, \f$ \lambda \f$ and \f$ \mu \f$ are Lame
3284  coefficients (see ElasticityIntegrator), \f$ u, v \f$ are the trial and test
3285  functions, respectively.
3286 
3287  The parameters \f$ \alpha \f$ and \f$ \kappa \f$ determine the DG method to
3288  use (when this integrator is added to the "broken" ElasticityIntegrator):
3289 
3290  - IIPG, \f$\alpha = 0\f$,
3291  C. Dawson, S. Sun, M. Wheeler, Compatible algorithms for coupled flow and
3292  transport, Comp. Meth. Appl. Mech. Eng., 193(23-26), 2565-2580, 2004.
3293 
3294  - SIPG, \f$\alpha = -1\f$,
3295  M. Grote, A. Schneebeli, D. Schotzau, Discontinuous Galerkin Finite
3296  %Element Method for the Wave Equation, SINUM, 44(6), 2408-2431, 2006.
3297 
3298  - NIPG, \f$\alpha = 1\f$,
3299  B. Riviere, M. Wheeler, V. Girault, A Priori Error Estimates for Finite
3300  %Element Methods Based on Discontinuous Approximation Spaces for Elliptic
3301  Problems, SINUM, 39(3), 902-931, 2001.
3302 
3303  This is a '%Vector' integrator, i.e. defined for FE spaces using multiple
3304  copies of a scalar FE space.
3305  */
3307 {
3308 public:
3309  DGElasticityIntegrator(double alpha_, double kappa_)
3310  : lambda(NULL), mu(NULL), alpha(alpha_), kappa(kappa_) { }
3311 
3313  double alpha_, double kappa_)
3314  : lambda(&lambda_), mu(&mu_), alpha(alpha_), kappa(kappa_) { }
3315 
3317  virtual void AssembleFaceMatrix(const FiniteElement &el1,
3318  const FiniteElement &el2,
3320  DenseMatrix &elmat);
3321 
3322 protected:
3324  double alpha, kappa;
3325 
3326 #ifndef MFEM_THREAD_SAFE
3327  // values of all scalar basis functions for one component of u (which is a
3328  // vector) at the integration point in the reference space
3330  // values of derivatives of all scalar basis functions for one component
3331  // of u (which is a vector) at the integration point in the reference space
3333  // Adjugate of the Jacobian of the transformation: adjJ = det(J) J^{-1}
3335  // gradient of shape functions in the real (physical, not reference)
3336  // coordinates, scaled by det(J):
3337  // dshape_ps(jdof,jm) = sum_{t} adjJ(t,jm)*dshape(jdof,t)
3339  Vector nor; // nor = |weight(J_face)| n
3340  Vector nL1, nL2; // nL1 = (lambda1 * ip.weight / detJ1) nor
3341  Vector nM1, nM2; // nM1 = (mu1 * ip.weight / detJ1) nor
3342  Vector dshape1_dnM, dshape2_dnM; // dshape1_dnM = dshape1_ps . nM1
3343  // 'jmat' corresponds to the term: kappa <h^{-1} {lambda + 2 mu} [u], [v]>
3345 #endif
3346 
3347  static void AssembleBlock(
3348  const int dim, const int row_ndofs, const int col_ndofs,
3349  const int row_offset, const int col_offset,
3350  const double jmatcoef, const Vector &col_nL, const Vector &col_nM,
3351  const Vector &row_shape, const Vector &col_shape,
3352  const Vector &col_dshape_dnM, const DenseMatrix &col_dshape,
3353  DenseMatrix &elmat, DenseMatrix &jmat);
3354 };
3355 
3356 /** Integrator for the DPG form: < v, [w] > over all faces (the interface) where
3357  the trial variable v is defined on the interface and the test variable w is
3358  defined inside the elements, generally in a DG space. */
3360 {
3361 private:
3362  Vector face_shape, shape1, shape2;
3363 
3364 public:
3367  virtual void AssembleFaceMatrix(const FiniteElement &trial_face_fe,
3368  const FiniteElement &test_fe1,
3369  const FiniteElement &test_fe2,
3371  DenseMatrix &elmat);
3372 };
3373 
3374 /** Integrator for the form: < v, [w.n] > over all faces (the interface) where
3375  the trial variable v is defined on the interface and the test variable w is
3376  in an H(div)-conforming space. */
3378 {
3379 private:
3380  Vector face_shape, normal, shape1_n, shape2_n;
3381  DenseMatrix shape1, shape2;
3382 
3383 public:
3386  virtual void AssembleFaceMatrix(const FiniteElement &trial_face_fe,
3387  const FiniteElement &test_fe1,
3388  const FiniteElement &test_fe2,
3390  DenseMatrix &elmat);
3391 };
3392 
3393 /** Integrator for the DPG form: < v, w > over a face (the interface) where
3394  the trial variable v is defined on the interface
3395  (H^-1/2 i.e., v:=uā‹…n normal trace of H(div))
3396  and the test variable w is in an H1-conforming space. */
3398 {
3399 private:
3400  Vector face_shape, shape;
3401 public:
3403  void AssembleTraceFaceMatrix(int elem,
3404  const FiniteElement &trial_face_fe,
3405  const FiniteElement &test_fe,
3407  DenseMatrix &elmat);
3408 };
3409 
3410 /** Integrator for the form: < v, w.n > over a face (the interface) where
3411  the trial variable v is defined on the interface (H^1/2, i.e., trace of H1)
3412  and the test variable w is in an H(div)-conforming space. */
3414 {
3415 private:
3416  Vector face_shape, normal, shape_n;
3417  DenseMatrix shape;
3418 
3419 public:
3421  virtual void AssembleTraceFaceMatrix(int ielem,
3422  const FiniteElement &trial_face_fe,
3423  const FiniteElement &test_fe,
3425  DenseMatrix &elmat);
3426 };
3427 
3428 
3429 /** Integrator for the form: < v, w Ɨ n > over a face (the interface)
3430  * In 3D the trial variable v is defined on the interface (H^-1/2(curl), trace of H(curl))
3431  * In 2D it's defined on the interface (H^1/2, trace of H1)
3432  * The test variable w is in an H(curl)-conforming space. */
3434 {
3435 private:
3436  DenseMatrix face_shape, shape, shape_n;
3437  Vector normal;
3438  Vector temp;
3439 
3440  void cross_product(const Vector & x, const DenseMatrix & Y, DenseMatrix & Z)
3441  {
3442  int dim = x.Size();
3443  MFEM_VERIFY(Y.Width() == dim, "Size mismatch");
3444  int dimc = dim == 3 ? dim : 1;
3445  int h = Y.Height();
3446  Z.SetSize(h,dimc);
3447  if (dim == 3)
3448  {
3449  for (int i = 0; i<h; i++)
3450  {
3451  Z(i,0) = x(2) * Y(i,1) - x(1) * Y(i,2);
3452  Z(i,1) = x(0) * Y(i,2) - x(2) * Y(i,0);
3453  Z(i,2) = x(1) * Y(i,0) - x(0) * Y(i,1);
3454  }
3455  }
3456  else
3457  {
3458  for (int i = 0; i<h; i++)
3459  {
3460  Z(i,0) = x(1) * Y(i,0) - x(0) * Y(i,1);
3461  }
3462  }
3463  }
3464 
3465 public:
3467  void AssembleTraceFaceMatrix(int elem,
3468  const FiniteElement &trial_face_fe,
3469  const FiniteElement &test_fe,
3471  DenseMatrix &elmat);
3472 };
3473 
3474 /** Abstract class to serve as a base for local interpolators to be used in the
3475  DiscreteLinearOperator class. */
3477 
3478 
3479 /** Class for constructing the gradient as a DiscreteLinearOperator from an
3480  H1-conforming space to an H(curl)-conforming space. The range space can be
3481  vector L2 space as well. */
3483 {
3484 public:
3485  GradientInterpolator() : dofquad_fe(NULL) { }
3486  virtual ~GradientInterpolator() { delete dofquad_fe; }
3487 
3488  virtual void AssembleElementMatrix2(const FiniteElement &h1_fe,
3489  const FiniteElement &nd_fe,
3490  ElementTransformation &Trans,
3491  DenseMatrix &elmat)
3492  { nd_fe.ProjectGrad(h1_fe, Trans, elmat); }
3493 
3495 
3496  /** @brief Setup method for PA data.
3497 
3498  @param[in] trial_fes H1 Lagrange space
3499  @param[in] test_fes H(curl) Nedelec space
3500  */
3501  virtual void AssemblePA(const FiniteElementSpace &trial_fes,
3502  const FiniteElementSpace &test_fes);
3503 
3504  virtual void AddMultPA(const Vector &x, Vector &y) const;
3505  virtual void AddMultTransposePA(const Vector &x, Vector &y) const;
3506 
3507 private:
3508  /// 1D finite element that generates and owns the 1D DofToQuad maps below
3509  FiniteElement *dofquad_fe;
3510 
3511  bool B_id; // is the B basis operator (maps_C_C) the identity?
3512  const DofToQuad *maps_C_C; // one-d map with Lobatto rows, Lobatto columns
3513  const DofToQuad *maps_O_C; // one-d map with Legendre rows, Lobatto columns
3514  int dim, ne, o_dofs1D, c_dofs1D;
3515 };
3516 
3517 
3518 /** Class for constructing the identity map as a DiscreteLinearOperator. This
3519  is the discrete embedding matrix when the domain space is a subspace of
3520  the range space. Otherwise, a dof projection matrix is constructed. */
3522 {
3523 public:
3524  IdentityInterpolator(): dofquad_fe(NULL) { }
3525 
3526  virtual void AssembleElementMatrix2(const FiniteElement &dom_fe,
3527  const FiniteElement &ran_fe,
3528  ElementTransformation &Trans,
3529  DenseMatrix &elmat)
3530  { ran_fe.Project(dom_fe, Trans, elmat); }
3531 
3533 
3534  virtual void AssemblePA(const FiniteElementSpace &trial_fes,
3535  const FiniteElementSpace &test_fes);
3536 
3537  virtual void AddMultPA(const Vector &x, Vector &y) const;
3538  virtual void AddMultTransposePA(const Vector &x, Vector &y) const;
3539 
3540  virtual ~IdentityInterpolator() { delete dofquad_fe; }
3541 
3542 private:
3543  /// 1D finite element that generates and owns the 1D DofToQuad maps below
3544  FiniteElement *dofquad_fe;
3545 
3546  const DofToQuad *maps_C_C; // one-d map with Lobatto rows, Lobatto columns
3547  const DofToQuad *maps_O_C; // one-d map with Legendre rows, Lobatto columns
3548  int dim, ne, o_dofs1D, c_dofs1D;
3549 
3550  Vector pa_data;
3551 };
3552 
3553 
3554 /** Class for constructing the (local) discrete curl matrix which can be used
3555  as an integrator in a DiscreteLinearOperator object to assemble the global
3556  discrete curl matrix. */
3558 {
3559 public:
3560  virtual void AssembleElementMatrix2(const FiniteElement &dom_fe,
3561  const FiniteElement &ran_fe,
3562  ElementTransformation &Trans,
3563  DenseMatrix &elmat)
3564  { ran_fe.ProjectCurl(dom_fe, Trans, elmat); }
3565 };
3566 
3567 
3568 /** Class for constructing the (local) discrete divergence matrix which can
3569  be used as an integrator in a DiscreteLinearOperator object to assemble
3570  the global discrete divergence matrix.
3571 
3572  Note: Since the dofs in the L2_FECollection are nodal values, the local
3573  discrete divergence matrix (with an RT-type domain space) will depend on
3574  the transformation. On the other hand, the local matrix returned by
3575  VectorFEDivergenceIntegrator is independent of the transformation. */
3577 {
3578 public:
3579  virtual void AssembleElementMatrix2(const FiniteElement &dom_fe,
3580  const FiniteElement &ran_fe,
3581  ElementTransformation &Trans,
3582  DenseMatrix &elmat)
3583  { ran_fe.ProjectDiv(dom_fe, Trans, elmat); }
3584 };
3585 
3586 
3587 /** A trace face interpolator class for interpolating the normal component of
3588  the domain space, e.g. vector H1, into the range space, e.g. the trace of
3589  RT which uses FiniteElement::INTEGRAL map type. */
3591 {
3592 public:
3593  virtual void AssembleElementMatrix2(const FiniteElement &dom_fe,
3594  const FiniteElement &ran_fe,
3595  ElementTransformation &Trans,
3596  DenseMatrix &elmat);
3597 };
3598 
3599 /** Interpolator of a scalar coefficient multiplied by a scalar field onto
3600  another scalar field. Note that this can produce inaccurate fields unless
3601  the target is sufficiently high order. */
3603 {
3604 public:
3606 
3607  virtual void AssembleElementMatrix2(const FiniteElement &dom_fe,
3608  const FiniteElement &ran_fe,
3609  ElementTransformation &Trans,
3610  DenseMatrix &elmat);
3611 
3612 protected:
3614 };
3615 
3616 /** Interpolator of a scalar coefficient multiplied by a vector field onto
3617  another vector field. Note that this can produce inaccurate fields unless
3618  the target is sufficiently high order. */
3620 {
3621 public:
3623  : Q(&sc) { }
3624 
3625  virtual void AssembleElementMatrix2(const FiniteElement &dom_fe,
3626  const FiniteElement &ran_fe,
3627  ElementTransformation &Trans,
3628  DenseMatrix &elmat);
3629 protected:
3631 };
3632 
3633 /** Interpolator of a vector coefficient multiplied by a scalar field onto
3634  another vector field. Note that this can produce inaccurate fields unless
3635  the target is sufficiently high order. */
3637 {
3638 public:
3640  : VQ(&vc) { }
3641 
3642  virtual void AssembleElementMatrix2(const FiniteElement &dom_fe,
3643  const FiniteElement &ran_fe,
3644  ElementTransformation &Trans,
3645  DenseMatrix &elmat);
3646 protected:
3648 };
3649 
3650 /** Interpolator of the 2D cross product between a vector coefficient and an
3651  H(curl)-conforming field onto an L2-conforming field. */
3653 {
3654 public:
3656  : VQ(&vc) { }
3657 
3658  virtual void AssembleElementMatrix2(const FiniteElement &nd_fe,
3659  const FiniteElement &l2_fe,
3660  ElementTransformation &Trans,
3661  DenseMatrix &elmat);
3662 protected:
3664 };
3665 
3666 /** Interpolator of the cross product between a vector coefficient and an
3667  H(curl)-conforming field onto an H(div)-conforming field. The range space
3668  can also be vector L2. */
3670 {
3671 public:
3673  : VQ(&vc) { }
3674 
3675  virtual void AssembleElementMatrix2(const FiniteElement &nd_fe,
3676  const FiniteElement &rt_fe,
3677  ElementTransformation &Trans,
3678  DenseMatrix &elmat);
3679 protected:
3681 };
3682 
3683 /** Interpolator of the inner product between a vector coefficient and an
3684  H(div)-conforming field onto an L2-conforming field. The range space can
3685  also be H1. */
3687 {
3688 public:
3690 
3691  virtual void AssembleElementMatrix2(const FiniteElement &rt_fe,
3692  const FiniteElement &l2_fe,
3693  ElementTransformation &Trans,
3694  DenseMatrix &elmat);
3695 protected:
3697 };
3698 
3699 }
3700 #endif
Abstract class for all finite elements.
Definition: fe_base.hpp:233
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 void AssemblePA(const FiniteElementSpace &trial_fes, const FiniteElementSpace &test_fes)
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:474
virtual const char * FiniteElementTypeFailureMessage() const
DiagonalMatrixCoefficient * DQ
Definition: bilininteg.hpp:601
virtual void AssembleElementMatrix(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &elmat)
Support for use in BilinearForm. Can be used only when appropriate.
Definition: bilininteg.hpp:541
virtual void AssembleMF(const FiniteElementSpace &fes)
Method defining matrix-free assembly.
virtual void AssembleEA(const FiniteElementSpace &fes, Vector &emat, const bool add)
Method defining element assembly.
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:219
virtual int GetTrialVDim(const FiniteElement &trial_fe)
virtual bool VerifyFiniteElementTypes(const FiniteElement &trial_fe, const FiniteElement &test_fe) const
Definition: bilininteg.hpp:488
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:96
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:84
virtual void CalcTestShape(const FiniteElement &test_fe, ElementTransformation &Trans, Vector &shape)
Definition: bilininteg.hpp:908
virtual void AddMultTransposePA(const Vector &, Vector &) const
Method for partially assembled transposed action.
MixedCrossCurlGradIntegrator(VectorCoefficient &vq)
virtual int GetTrialVDim(const FiniteElement &trial_fe)
bool symmetric
False if using a nonsymmetric matrix coefficient.
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:245
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:123
Integrator defining a sum of multiple Integrators.
Definition: bilininteg.hpp:391
SumIntegrator(int own_integs=1)
Definition: bilininteg.hpp:399
virtual void AssembleFaceMatrix(const FiniteElement &el1, const FiniteElement &el2, FaceElementTransformations &Trans, DenseMatrix &elmat)
Definition: bilininteg.cpp:307
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
virtual int OrderW() const =0
Return the order of the determinant of the Jacobian (weight) of the transformation.
MixedScalarCrossProductIntegrator(VectorCoefficient &vq)
bool SupportsCeed() const
Indicates whether this integrator can use a Ceed backend.
virtual void AssemblePABoundary(const FiniteElementSpace &fes)
Definition: bilininteg.cpp:42
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:238
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 AssembleElementMatrix2(const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &Trans, DenseMatrix &elmat)
virtual void CalcTrialShape(const FiniteElement &trial_fe, ElementTransformation &Trans, DenseMatrix &shape)
VectorCoefficient * u
virtual int GetIntegrationOrder(const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &Trans)
Definition: bilininteg.hpp:674
virtual int GetTestVDim(const FiniteElement &test_fe)
Definition: bilininteg.hpp:582
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:316
virtual void CalcShape(const FiniteElement &scalar_fe, ElementTransformation &Trans, Vector &shape)
ElasticityIntegrator(Coefficient &l, Coefficient &m)
virtual void AddMultPA(const Vector &, Vector &) const
Method for partially assembled action.
virtual void AddMultTransposePA(const Vector &x, Vector &y) const
Method for partially assembled transposed action.
Definition: bilininteg.cpp:376
virtual int GetTestVDim(const FiniteElement &test_fe)
bool SupportsCeed() const
Indicates whether this integrator can use a Ceed backend.
virtual void AssemblePA(const FiniteElementSpace &fes)
Method defining partial assembly.
virtual void AssembleMF(const FiniteElementSpace &fes)
Method defining matrix-free assembly.
VectorCoefficient * VQ
Definition: bilininteg.hpp:600
VectorDiffusionIntegrator(Coefficient &q)
virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &Trans, DenseMatrix &elmat)
virtual const char * FiniteElementTypeFailureMessage() const
virtual bool VerifyFiniteElementTypes(const FiniteElement &trial_fe, const FiniteElement &test_fe) const
Definition: bilininteg.hpp:775
int Space() const
Returns the type of FunctionSpace on the element.
Definition: fe_base.hpp:337
virtual const char * FiniteElementTypeFailureMessage() const
virtual void CalcTestShape(const FiniteElement &test_fe, ElementTransformation &Trans, Vector &shape)
Definition: bilininteg.hpp:508
virtual bool VerifyFiniteElementTypes(const FiniteElement &trial_fe, const FiniteElement &test_fe) const
const GeometricFactors * geom
Not owned.
virtual void AddMultTransposePA(const Vector &x, Vector &y) const
Method for partially assembled transposed action.
virtual void CalcTrialShape(const FiniteElement &trial_fe, ElementTransformation &Trans, Vector &shape)
Definition: bilininteg.hpp:756
virtual void AssembleDiagonalMF(Vector &diag)
Assemble diagonal and add it to Vector diag.
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:813
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:792
virtual void AssembleFaceMatrix(const FiniteElement &el1, const FiniteElement &el2, FaceElementTransformations &Trans, DenseMatrix &elmat)
MixedScalarVectorIntegrator(VectorCoefficient &vq, bool transpose_=false, bool cross_2d_=false)
Definition: bilininteg.hpp:641
MixedScalarMassIntegrator(Coefficient &q)
Definition: bilininteg.hpp:716
int GetCurlDim() const
Returns the dimension of the curl for vector-valued finite elements.
Definition: fe_base.hpp:317
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 &x, Vector &y) const
Definition: bilininteg.cpp:392
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
virtual void AssemblePA(const FiniteElementSpace &fes)
Method defining partial assembly.
Definition: bilininteg.cpp:336
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:258
virtual void AssembleElementMatrix(const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
Given a particular Finite Element computes the element matrix elmat.
ScalarCrossProductInterpolator(VectorCoefficient &vc)
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &Trans, DenseMatrix &elmat)
virtual void AssemblePA(const FiniteElementSpace &fes)
Method defining partial assembly.
MixedGradGradIntegrator(DiagonalMatrixCoefficient &dq)
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
virtual void AssemblePAInteriorFaces(const FiniteElementSpace &fes)
Definition: bilininteg.cpp:352
virtual int GetIntegrationOrder(const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &Trans)
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:129
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:236
virtual void AssembleElementMatrix(const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
Given a particular Finite Element computes the element matrix elmat.
Definition: bilininteg.cpp:251
virtual void AssembleDiagonalPA(Vector &diag)
Assemble diagonal and add it to Vector diag.
const DofToQuad * maps
Not owned.
virtual bool VerifyFiniteElementTypes(const FiniteElement &trial_fe, const FiniteElement &test_fe) const
BilinearFormIntegrator(const IntegrationRule *ir=NULL)
Definition: bilininteg.hpp:27
virtual void AssemblePA(const FiniteElementSpace &)
Method defining partial assembly.
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:305
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 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 AddMultMF(const Vector &x, Vector &y) const
virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &Trans, DenseMatrix &elmat)
Definition: bilininteg.cpp:765
void AddIntegrator(BilinearFormIntegrator *integ)
Definition: bilininteg.hpp:403
const DofToQuad * mapsC
Not owned. DOF-to-quad map, closed.
virtual void CalcShape(const FiniteElement &scalar_fe, ElementTransformation &Trans, Vector &shape_)
Definition: bilininteg.hpp:688
Structure for storing mesh geometric factors: coordinates, Jacobians, and determinants of the Jacobia...
Definition: mesh.hpp:2183
virtual const char * FiniteElementTypeFailureMessage() const
Definition: bilininteg.hpp:857
DiffusionIntegrator(VectorCoefficient &q, const IntegrationRule *ir=nullptr)
Construct a diffusion integrator with a vector coefficient q.
virtual void AddMultPA(const Vector &x, Vector &y) const
Method for partially assembled action.
virtual void AddMultPA(const Vector &x, Vector &y) const
Method for partially assembled action.
virtual const char * FiniteElementTypeFailureMessage() const
const FaceGeometricFactors * face_geom
Not owned.
virtual void AssembleTraceFaceMatrix(int ielem, const FiniteElement &trial_face_fe, const FiniteElement &test_fe, FaceElementTransformations &Trans, DenseMatrix &elmat)
virtual int GetTestVDim(const FiniteElement &test_fe)
virtual int GetVDim(const FiniteElement &vector_fe)
virtual void AddMultNURBSPA(const Vector &, Vector &) const
Method for partially assembled action on NURBS patches.
virtual void CalcTrialShape(const FiniteElement &trial_fe, ElementTransformation &Trans, Vector &shape)
Definition: bilininteg.hpp:513
virtual const char * FiniteElementTypeFailureMessage() const
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:264
MixedCrossProductIntegrator(VectorCoefficient &vq)
virtual void AssemblePA(const FiniteElementSpace &fes)
Method defining partial assembly.
virtual void AddMultNURBSPA(const Vector &x, Vector &y) const
Method for partially assembled action on NURBS patches.
Definition: bilininteg.cpp:105
VectorDiffusionIntegrator(int vector_dimension)
Integrator with unit coefficient for caller-specified vector dimension.
DGDiffusionBR2Integrator(class FiniteElementSpace &fes, double e=1.0)
virtual void AddMultPA(const Vector &x, Vector &y) const
Method for partially assembled action.
virtual const char * FiniteElementTypeFailureMessage() const
Definition: bilininteg.hpp:993
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)
virtual void AddMultPA(const Vector &, Vector &) const
Method for partially assembled action.
MixedVectorWeakDivergenceIntegrator(DiagonalMatrixCoefficient &dq)
MixedVectorGradientIntegrator(Coefficient &q)
(Q div u, div v) for RT elements
virtual void AddMultPA(const Vector &x, Vector &y) const
Method for partially assembled action.
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:317
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 AddMultTransposePA(const Vector &, Vector &) const
Method for partially assembled transposed action.
virtual void AssembleDiagonalPA(Vector &diag)
Assemble diagonal and add it to Vector diag.
Definition: bilininteg.cpp:344
Coefficient * GetCoefficient() const
virtual bool VerifyFiniteElementTypes(const FiniteElement &trial_fe, const FiniteElement &test_fe) const
virtual void AssemblePA(const FiniteElementSpace &trial_fes, const FiniteElementSpace &test_fes)
Setup method for PA data.
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:457
virtual void CalcTestShape(const FiniteElement &test_fe, ElementTransformation &Trans, DenseMatrix &shape)
const GeometricFactors * geom
Not owned.
MixedCurlCurlIntegrator(MatrixCoefficient &mq)
virtual void AssembleNURBSPA(const FiniteElementSpace &fes)
Method defining partial assembly on NURBS patches.
MatrixCoefficient * MQ
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 AddMultPA(const Vector &x, Vector &y) const
Method for partially assembled action.
virtual void AssembleDiagonalPA(Vector &diag)
Assemble diagonal and add it to Vector diag.
virtual void AssembleEA(const FiniteElementSpace &fes, Vector &emat, const bool add=true)
Method defining element assembly.
Definition: bilininteg.cpp:66
VectorMassIntegrator(VectorCoefficient &q, int qo=0)
Construct an integrator with diagonal coefficient q.
virtual void AssembleMF(const FiniteElementSpace &fes)
Method defining matrix-free assembly.
virtual void AddMultPA(const Vector &x, Vector &y) const
Method for partially assembled action.
Definition: bilininteg.cpp:99
virtual void AssembleFaceMatrix(const FiniteElement &el1, const FiniteElement &el2, FaceElementTransformations &Trans, DenseMatrix &elmat)
Structure for storing face geometric factors: coordinates, Jacobians, determinants of the Jacobians...
Definition: mesh.hpp:2237
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 bool VerifyFiniteElementTypes(const FiniteElement &trial_fe, const FiniteElement &test_fe) const
Definition: bilininteg.hpp:645
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:126
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:221
virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &Trans, DenseMatrix &elmat)
Definition: bilininteg.cpp:690
virtual void AssembleMF(const FiniteElementSpace &fes)
Method defining matrix-free assembly.
virtual void AssemblePA(const FiniteElementSpace &trial_fes, const FiniteElementSpace &test_fes)
virtual void AssembleElementMatrix2(const FiniteElement &dom_fe, const FiniteElement &ran_fe, ElementTransformation &Trans, DenseMatrix &elmat)
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:227
Integrator that inverts the matrix assembled by another integrator.
Definition: bilininteg.hpp:371
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 void AssemblePA(const FiniteElementSpace &fes)
Method defining partial assembly.
virtual bool VerifyFiniteElementTypes(const FiniteElement &trial_fe, const FiniteElement &test_fe) const
const Coefficient * GetCoefficient() const
virtual void AssembleNURBSPA(const FiniteElementSpace &fes)
Method defining partial assembly on NURBS patches.
Definition: bilininteg.cpp:29
virtual void AssembleDiagonalMF(Vector &diag)
Assemble diagonal and add it to Vector diag.
Definition: bilininteg.cpp:408
virtual void AddMultTransposePA(const Vector &x, Vector &y) const
Method for partially assembled transposed action.
virtual void AssemblePA(const FiniteElementSpace &trial_fes, const FiniteElementSpace &test_fes)
Definition: bilininteg.hpp:310
-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:749
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:93
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
virtual void AssembleDiagonalPA(Vector &diag)
Assemble diagonal and add it to Vector diag.
virtual void AssemblePA(const FiniteElementSpace &trial_fes, const FiniteElementSpace &test_fes)
VectorDiffusionIntegrator(MatrixCoefficient &mq)
Integrator with MatrixCoefficient. The vector dimension of the FiniteElementSpace is assumed to be th...
VectorDivergenceIntegrator(Coefficient &q)
virtual void AssemblePA(const FiniteElementSpace &fes)
Method defining partial assembly.
virtual void AssembleEA(const FiniteElementSpace &fes, Vector &emat, const bool add)
Method defining element assembly.
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:264
virtual int GetTestVDim(const FiniteElement &test_fe)
virtual const char * FiniteElementTypeFailureMessage() const
Definition: bilininteg.hpp:658
virtual const char * FiniteElementTypeFailureMessage() const
virtual void AddMultPA(const Vector &, Vector &) const
Method for partially assembled action.
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
Data type sparse matrix.
Definition: sparsemat.hpp:50
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:759
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:271
MixedScalarWeakCurlIntegrator(Coefficient &q)
Definition: bilininteg.hpp:980
virtual void AssemblePA(const FiniteElementSpace &fes)
Method defining partial assembly.
MixedCurlCurlIntegrator(Coefficient &q)
GradientIntegrator(Coefficient &q)
virtual void AddMultMF(const Vector &, Vector &) 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 void AssemblePA(const FiniteElementSpace &trial_fes, const FiniteElementSpace &test_fes)
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:280
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:384
TransposeIntegrator(BilinearFormIntegrator *bfi_, int own_bfi_=1)
Definition: bilininteg.hpp:283
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun, DenseMatrix &elmat)
Assemble the local gradient matrix.
Definition: bilininteg.hpp:193
static const IntegrationRule & GetRule(const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &Trans)
MixedVectorWeakCurlIntegrator(DiagonalMatrixCoefficient &dq)
VectorFEDivergenceIntegrator(Coefficient &q)
virtual void AssembleEAInteriorFaces(const FiniteElementSpace &fes, Vector &ea_data_int, Vector &ea_data_ext, const bool add)
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
virtual void AddMultPA(const Vector &x, Vector &y) const
Method for partially assembled action.
VectorFEMassIntegrator(MatrixCoefficient *mq_)
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:71
VectorFECurlIntegrator(Coefficient &q)
virtual void CalcVShape(const FiniteElement &vector_fe, ElementTransformation &Trans, DenseMatrix &shape_)
Definition: bilininteg.hpp:683
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:368
virtual void AssemblePatchMatrix(const int patch, const FiniteElementSpace &fes, SparseMatrix *&smat)
Definition: bilininteg.cpp:157
MixedVectorMassIntegrator(MatrixCoefficient &mq)
virtual void AssemblePatchMatrix(const int patch, const FiniteElementSpace &fes, SparseMatrix *&smat)
virtual void AssembleElementMatrix(const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
Given a particular Finite Element computes the element matrix elmat.
MixedScalarDivergenceIntegrator(Coefficient &q)
Defini