MFEM  v4.2.0 Finite element discretization library
ex18.hpp
Go to the documentation of this file.
1 // MFEM Example 18 - Serial/Parallel Shared Code
2
3 #include "mfem.hpp"
4
5 using namespace std;
6 using namespace mfem;
7
8 // Problem definition
9 extern int problem;
10
11 // Maximum characteristic speed (updated by integrators)
12 extern double max_char_speed;
13
14 extern const int num_equation;
15 extern const double specific_heat_ratio;
16 extern const double gas_constant;
17
18 // Time-dependent operator for the right-hand side of the ODE representing the
19 // DG weak form.
21 {
22 private:
23  const int dim;
24
25  FiniteElementSpace &vfes;
26  Operator &A;
27  SparseMatrix &Aflux;
28  DenseTensor Me_inv;
29
30  mutable Vector state;
31  mutable DenseMatrix f;
32  mutable DenseTensor flux;
33  mutable Vector z;
34
35  void GetFlux(const DenseMatrix &state, DenseTensor &flux) const;
36
37 public:
39  Operator &_A, SparseMatrix &_Aflux);
40
41  virtual void Mult(const Vector &x, Vector &y) const;
42
43  virtual ~FE_Evolution() { }
44 };
45
46 // Implements a simple Rusanov flux
48 {
49 private:
50  Vector flux1;
51  Vector flux2;
52
53 public:
54  RiemannSolver();
55  double Eval(const Vector &state1, const Vector &state2,
56  const Vector &nor, Vector &flux);
57 };
58
59
60 // Constant (in time) mixed bilinear form multiplying the flux grid function.
61 // The form is (vec(v), grad(w)) where the trial space = vector L2 space (mesh
62 // dim) and test space = scalar L2 space.
64 {
65 private:
66  Vector shape;
67  DenseMatrix flux;
68  DenseMatrix dshapedr;
69  DenseMatrix dshapedx;
70
71 public:
72  DomainIntegrator(const int dim);
73
74  virtual void AssembleElementMatrix2(const FiniteElement &trial_fe,
75  const FiniteElement &test_fe,
77  DenseMatrix &elmat);
78 };
79
80 // Interior face term: <F.n(u),[w]>
82 {
83 private:
84  RiemannSolver rsolver;
85  Vector shape1;
86  Vector shape2;
87  Vector funval1;
88  Vector funval2;
89  Vector nor;
90  Vector fluxN;
91
92 public:
93  FaceIntegrator(RiemannSolver &rsolver_, const int dim);
94
95  virtual void AssembleFaceVector(const FiniteElement &el1,
96  const FiniteElement &el2,
98  const Vector &elfun, Vector &elvect);
99 };
100
101 // Implementation of class FE_Evolution
103  Operator &_A, SparseMatrix &_Aflux)
104  : TimeDependentOperator(_A.Height()),
105  dim(_vfes.GetFE(0)->GetDim()),
106  vfes(_vfes),
107  A(_A),
108  Aflux(_Aflux),
109  Me_inv(vfes.GetFE(0)->GetDof(), vfes.GetFE(0)->GetDof(), vfes.GetNE()),
110  state(num_equation),
111  f(num_equation, dim),
112  flux(vfes.GetNDofs(), dim, num_equation),
113  z(A.Height())
114 {
115  // Standard local assembly and inversion for energy mass matrices.
116  const int dof = vfes.GetFE(0)->GetDof();
117  DenseMatrix Me(dof);
118  DenseMatrixInverse inv(&Me);
119  MassIntegrator mi;
120  for (int i = 0; i < vfes.GetNE(); i++)
121  {
122  mi.AssembleElementMatrix(*vfes.GetFE(i), *vfes.GetElementTransformation(i), Me);
123  inv.Factor();
124  inv.GetInverseMatrix(Me_inv(i));
125  }
126 }
127
128 void FE_Evolution::Mult(const Vector &x, Vector &y) const
129 {
130  // 0. Reset wavespeed computation before operator application.
131  max_char_speed = 0.;
132
133  // 1. Create the vector z with the face terms -<F.n(u), [w]>.
134  A.Mult(x, z);
135
136  // 2. Add the element terms.
137  // i. computing the flux approximately as a grid function by interpolating
138  // at the solution nodes.
139  // ii. multiplying this grid function by a (constant) mixed bilinear form for
140  // each of the num_equation, computing (F(u), grad(w)) for each equation.
141
142  DenseMatrix xmat(x.GetData(), vfes.GetNDofs(), num_equation);
143  GetFlux(xmat, flux);
144
145  for (int k = 0; k < num_equation; k++)
146  {
147  Vector fk(flux(k).GetData(), dim * vfes.GetNDofs());
148  Vector zk(z.GetData() + k * vfes.GetNDofs(), vfes.GetNDofs());
150  }
151
152  // 3. Multiply element-wise by the inverse mass matrices.
153  Vector zval;
154  Array<int> vdofs;
155  const int dof = vfes.GetFE(0)->GetDof();
156  DenseMatrix zmat, ymat(dof, num_equation);
157
158  for (int i = 0; i < vfes.GetNE(); i++)
159  {
160  // Return the vdofs ordered byNODES
161  vfes.GetElementVDofs(i, vdofs);
162  z.GetSubVector(vdofs, zval);
163  zmat.UseExternalData(zval.GetData(), dof, num_equation);
164  mfem::Mult(Me_inv(i), zmat, ymat);
165  y.SetSubVector(vdofs, ymat.GetData());
166  }
167 }
168
169 // Physicality check (at end)
170 bool StateIsPhysical(const Vector &state, const int dim);
171
172 // Pressure (EOS) computation
173 inline double ComputePressure(const Vector &state, int dim)
174 {
175  const double den = state(0);
176  const Vector den_vel(state.GetData() + 1, dim);
177  const double den_energy = state(1 + dim);
178
179  double den_vel2 = 0;
180  for (int d = 0; d < dim; d++) { den_vel2 += den_vel(d) * den_vel(d); }
181  den_vel2 /= den;
182
183  return (specific_heat_ratio - 1.0) * (den_energy - 0.5 * den_vel2);
184 }
185
186 // Compute the vector flux F(u)
187 void ComputeFlux(const Vector &state, int dim, DenseMatrix &flux)
188 {
189  const double den = state(0);
190  const Vector den_vel(state.GetData() + 1, dim);
191  const double den_energy = state(1 + dim);
192
193  MFEM_ASSERT(StateIsPhysical(state, dim), "");
194
195  const double pres = ComputePressure(state, dim);
196
197  for (int d = 0; d < dim; d++)
198  {
199  flux(0, d) = den_vel(d);
200  for (int i = 0; i < dim; i++)
201  {
202  flux(1+i, d) = den_vel(i) * den_vel(d) / den;
203  }
204  flux(1+d, d) += pres;
205  }
206
207  const double H = (den_energy + pres) / den;
208  for (int d = 0; d < dim; d++)
209  {
210  flux(1+dim, d) = den_vel(d) * H;
211  }
212 }
213
214 // Compute the scalar F(u).n
215 void ComputeFluxDotN(const Vector &state, const Vector &nor,
216  Vector &fluxN)
217 {
218  // NOTE: nor in general is not a unit normal
219  const int dim = nor.Size();
220  const double den = state(0);
221  const Vector den_vel(state.GetData() + 1, dim);
222  const double den_energy = state(1 + dim);
223
224  MFEM_ASSERT(StateIsPhysical(state, dim), "");
225
226  const double pres = ComputePressure(state, dim);
227
228  double den_velN = 0;
229  for (int d = 0; d < dim; d++) { den_velN += den_vel(d) * nor(d); }
230
231  fluxN(0) = den_velN;
232  for (int d = 0; d < dim; d++)
233  {
234  fluxN(1+d) = den_velN * den_vel(d) / den + pres * nor(d);
235  }
236
237  const double H = (den_energy + pres) / den;
238  fluxN(1 + dim) = den_velN * H;
239 }
240
241 // Compute the maximum characteristic speed.
242 inline double ComputeMaxCharSpeed(const Vector &state, const int dim)
243 {
244  const double den = state(0);
245  const Vector den_vel(state.GetData() + 1, dim);
246
247  double den_vel2 = 0;
248  for (int d = 0; d < dim; d++) { den_vel2 += den_vel(d) * den_vel(d); }
249  den_vel2 /= den;
250
251  const double pres = ComputePressure(state, dim);
252  const double sound = sqrt(specific_heat_ratio * pres / den);
253  const double vel = sqrt(den_vel2 / den);
254
255  return vel + sound;
256 }
257
258 // Compute the flux at solution nodes.
259 void FE_Evolution::GetFlux(const DenseMatrix &x, DenseTensor &flux) const
260 {
261  const int dof = flux.SizeI();
262  const int dim = flux.SizeJ();
263
264  for (int i = 0; i < dof; i++)
265  {
266  for (int k = 0; k < num_equation; k++) { state(k) = x(i, k); }
267  ComputeFlux(state, dim, f);
268
269  for (int d = 0; d < dim; d++)
270  {
271  for (int k = 0; k < num_equation; k++)
272  {
273  flux(i, d, k) = f(k, d);
274  }
275  }
276
277  // Update max char speed
278  const double mcs = ComputeMaxCharSpeed(state, dim);
279  if (mcs > max_char_speed) { max_char_speed = mcs; }
280  }
281 }
282
283 // Implementation of class RiemannSolver
285  flux1(num_equation),
286  flux2(num_equation) { }
287
288 double RiemannSolver::Eval(const Vector &state1, const Vector &state2,
289  const Vector &nor, Vector &flux)
290 {
291  // NOTE: nor in general is not a unit normal
292  const int dim = nor.Size();
293
294  MFEM_ASSERT(StateIsPhysical(state1, dim), "");
295  MFEM_ASSERT(StateIsPhysical(state2, dim), "");
296
297  const double maxE1 = ComputeMaxCharSpeed(state1, dim);
298  const double maxE2 = ComputeMaxCharSpeed(state2, dim);
299
300  const double maxE = max(maxE1, maxE2);
301
302  ComputeFluxDotN(state1, nor, flux1);
303  ComputeFluxDotN(state2, nor, flux2);
304
305  double normag = 0;
306  for (int i = 0; i < dim; i++)
307  {
308  normag += nor(i) * nor(i);
309  }
310  normag = sqrt(normag);
311
312  for (int i = 0; i < num_equation; i++)
313  {
314  flux(i) = 0.5 * (flux1(i) + flux2(i))
315  - 0.5 * maxE * (state2(i) - state1(i)) * normag;
316  }
317
318  return maxE;
319 }
320
321 // Implementation of class DomainIntegrator
322 DomainIntegrator::DomainIntegrator(const int dim) : flux(num_equation, dim) { }
323
325  const FiniteElement &test_fe,
327  DenseMatrix &elmat)
328 {
329  // Assemble the form (vec(v), grad(w))
330
331  // Trial space = vector L2 space (mesh dim)
332  // Test space = scalar L2 space
333
334  const int dof_trial = trial_fe.GetDof();
335  const int dof_test = test_fe.GetDof();
336  const int dim = trial_fe.GetDim();
337
338  shape.SetSize(dof_trial);
339  dshapedr.SetSize(dof_test, dim);
340  dshapedx.SetSize(dof_test, dim);
341
342  elmat.SetSize(dof_test, dof_trial * dim);
343  elmat = 0.0;
344
345  const int maxorder = max(trial_fe.GetOrder(), test_fe.GetOrder());
346  const int intorder = 2 * maxorder;
347  const IntegrationRule *ir = &IntRules.Get(trial_fe.GetGeomType(), intorder);
348
349  for (int i = 0; i < ir->GetNPoints(); i++)
350  {
351  const IntegrationPoint &ip = ir->IntPoint(i);
352
353  // Calculate the shape functions
354  trial_fe.CalcShape(ip, shape);
355  shape *= ip.weight;
356
357  // Compute the physical gradients of the test functions
358  Tr.SetIntPoint(&ip);
359  test_fe.CalcDShape(ip, dshapedr);
361
362  for (int d = 0; d < dim; d++)
363  {
364  for (int j = 0; j < dof_test; j++)
365  {
366  for (int k = 0; k < dof_trial; k++)
367  {
368  elmat(j, k + d * dof_trial) += shape(k) * dshapedx(j, d);
369  }
370  }
371  }
372  }
373 }
374
375 // Implementation of class FaceIntegrator
377  rsolver(rsolver_),
378  funval1(num_equation),
379  funval2(num_equation),
380  nor(dim),
381  fluxN(num_equation) { }
382
384  const FiniteElement &el2,
386  const Vector &elfun, Vector &elvect)
387 {
388  // Compute the term <F.n(u),[w]> on the interior faces.
389  const int dof1 = el1.GetDof();
390  const int dof2 = el2.GetDof();
391
392  shape1.SetSize(dof1);
393  shape2.SetSize(dof2);
394
395  elvect.SetSize((dof1 + dof2) * num_equation);
396  elvect = 0.0;
397
398  DenseMatrix elfun1_mat(elfun.GetData(), dof1, num_equation);
399  DenseMatrix elfun2_mat(elfun.GetData() + dof1 * num_equation, dof2,
400  num_equation);
401
402  DenseMatrix elvect1_mat(elvect.GetData(), dof1, num_equation);
403  DenseMatrix elvect2_mat(elvect.GetData() + dof1 * num_equation, dof2,
404  num_equation);
405
406  // Integration order calculation from DGTraceIntegrator
407  int intorder;
408  if (Tr.Elem2No >= 0)
409  intorder = (min(Tr.Elem1->OrderW(), Tr.Elem2->OrderW()) +
410  2*max(el1.GetOrder(), el2.GetOrder()));
411  else
412  {
413  intorder = Tr.Elem1->OrderW() + 2*el1.GetOrder();
414  }
415  if (el1.Space() == FunctionSpace::Pk)
416  {
417  intorder++;
418  }
419  const IntegrationRule *ir = &IntRules.Get(Tr.GetGeometryType(), intorder);
420
421  for (int i = 0; i < ir->GetNPoints(); i++)
422  {
423  const IntegrationPoint &ip = ir->IntPoint(i);
424
425  Tr.SetAllIntPoints(&ip); // set face and element int. points
426
427  // Calculate basis functions on both elements at the face
428  el1.CalcShape(Tr.GetElement1IntPoint(), shape1);
429  el2.CalcShape(Tr.GetElement2IntPoint(), shape2);
430
431  // Interpolate elfun at the point
432  elfun1_mat.MultTranspose(shape1, funval1);
433  elfun2_mat.MultTranspose(shape2, funval2);
434
435  // Get the normal vector and the flux on the face
436  CalcOrtho(Tr.Jacobian(), nor);
437  const double mcs = rsolver.Eval(funval1, funval2, nor, fluxN);
438
439  // Update max char speed
440  if (mcs > max_char_speed) { max_char_speed = mcs; }
441
442  fluxN *= ip.weight;
443  for (int k = 0; k < num_equation; k++)
444  {
445  for (int s = 0; s < dof1; s++)
446  {
447  elvect1_mat(s, k) -= fluxN(k) * shape1(s);
448  }
449  for (int s = 0; s < dof2; s++)
450  {
451  elvect2_mat(s, k) += fluxN(k) * shape2(s);
452  }
453  }
454  }
455 }
456
457 // Check that the state is physical - enabled in debug mode
458 bool StateIsPhysical(const Vector &state, const int dim)
459 {
460  const double den = state(0);
461  const Vector den_vel(state.GetData() + 1, dim);
462  const double den_energy = state(1 + dim);
463
464  if (den < 0)
465  {
466  cout << "Negative density: ";
467  for (int i = 0; i < state.Size(); i++)
468  {
469  cout << state(i) << " ";
470  }
471  cout << endl;
472  return false;
473  }
474  if (den_energy <= 0)
475  {
476  cout << "Negative energy: ";
477  for (int i = 0; i < state.Size(); i++)
478  {
479  cout << state(i) << " ";
480  }
481  cout << endl;
482  return false;
483  }
484
485  double den_vel2 = 0;
486  for (int i = 0; i < dim; i++) { den_vel2 += den_vel(i) * den_vel(i); }
487  den_vel2 /= den;
488
489  const double pres = (specific_heat_ratio - 1.0) * (den_energy - 0.5 * den_vel2);
490
491  if (pres <= 0)
492  {
493  cout << "Negative pressure: " << pres << ", state: ";
494  for (int i = 0; i < state.Size(); i++)
495  {
496  cout << state(i) << " ";
497  }
498  cout << endl;
499  return false;
500  }
501  return true;
502 }
503
504 // Initial condition
505 void InitialCondition(const Vector &x, Vector &y)
506 {
507  MFEM_ASSERT(x.Size() == 2, "");
508
509  double radius = 0, Minf = 0, beta = 0;
510  if (problem == 1)
511  {
512  // "Fast vortex"
514  Minf = 0.5;
515  beta = 1. / 5.;
516  }
517  else if (problem == 2)
518  {
519  // "Slow vortex"
521  Minf = 0.05;
522  beta = 1. / 50.;
523  }
524  else
525  {
526  mfem_error("Cannot recognize problem."
527  "Options are: 1 - fast vortex, 2 - slow vortex");
528  }
529
530  const double xc = 0.0, yc = 0.0;
531
532  // Nice units
533  const double vel_inf = 1.;
534  const double den_inf = 1.;
535
536  // Derive remainder of background state from this and Minf
537  const double pres_inf = (den_inf / specific_heat_ratio) * (vel_inf / Minf) *
538  (vel_inf / Minf);
539  const double temp_inf = pres_inf / (den_inf * gas_constant);
540
542  r2rad += (x(0) - xc) * (x(0) - xc);
543  r2rad += (x(1) - yc) * (x(1) - yc);
545
546  const double shrinv1 = 1.0 / (specific_heat_ratio - 1.);
547
548  const double velX = vel_inf * (1 - beta * (x(1) - yc) / radius * exp(
550  const double velY = vel_inf * beta * (x(0) - xc) / radius * exp(-0.5 * r2rad);
551  const double vel2 = velX * velX + velY * velY;
552
553  const double specific_heat = gas_constant * specific_heat_ratio * shrinv1;
554  const double temp = temp_inf - 0.5 * (vel_inf * beta) *
555  (vel_inf * beta) / specific_heat * exp(-r2rad);
556
557  const double den = den_inf * pow(temp/temp_inf, shrinv1);
558  const double pres = den * gas_constant * temp;
559  const double energy = shrinv1 * pres / den + 0.5 * vel2;
560
561  y(0) = den;
562  y(1) = den * velX;
563  y(2) = den * velY;
564  y(3) = den * energy;
565 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:245
Abstract class for all finite elements.
Definition: fe.hpp:235
Return the adjugate of the Jacobian matrix of the transformation at the currently set IntegrationPoin...
Definition: eltrans.hpp:127
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:521
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe.hpp:309
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:397
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
virtual int OrderW() const =0
Return the order of the determinant of the Jacobian (weight) of the transformation.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:915
Geometry::Type GetGeometryType() const
Return the Geometry::Type of the reference element.
Definition: eltrans.hpp:149
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:459
void GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:173
Base abstract class for first order time dependent operators.
Definition: operator.hpp:267
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:476
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:85
A specialized ElementTransformation class representing a face and its two neighboring elements...
Definition: eltrans.hpp:467
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:495
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe.hpp:319
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
void AddMult(const Vector &x, Vector &y, const double a=1.0) const
y += A * x (default) or y += a * A * x
Definition: sparsemat.cpp:601
int Size() const
Returns the size of the vector.
Definition: vector.hpp:160
FaceIntegrator(RiemannSolver &rsolver_, const int dim)
Definition: ex18.hpp:376
int Space() const
Returns the type of FunctionSpace on the element.
Definition: fe.hpp:329
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
const IntegrationPoint & GetElement2IntPoint()
Get a const reference to the integration point in neighboring element 2 corresponding to the currentl...
Definition: eltrans.hpp:574
void GetInverseMatrix(DenseMatrix &Ainv) const
Compute and return the inverse matrix in Ainv.
Definition: densemat.cpp:3249
void Factor()
Factor the current DenseMatrix, *a.
Definition: densemat.cpp:3237
double * GetData() const
Returns the matrix data array.
Definition: densemat.hpp:100
ElementTransformation * Elem2
Definition: eltrans.hpp:509
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:169
void CalcOrtho(const DenseMatrix &J, Vector &n)
Definition: densemat.cpp:2305
double max_char_speed
Definition: ex18.cpp:59
RiemannSolver()
Definition: ex18.hpp:284
FE_Evolution(FiniteElementSpace &_vfes, Operator &_A, SparseMatrix &_Aflux)
Definition: ex18.hpp:102
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe.hpp:312
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:427
double Eval(const Vector &state1, const Vector &state2, const Vector &nor, Vector &flux)
Definition: ex18.hpp:288
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:248
void ComputeFluxDotN(const Vector &state, const Vector &nor, Vector &fluxN)
Definition: ex18.hpp:215
virtual void AssembleElementMatrix(const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
Definition: bilininteg.cpp:891
virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &Tr, DenseMatrix &elmat)
Definition: ex18.hpp:324
Data type sparse matrix.
Definition: sparsemat.hpp:46
const int num_equation
Definition: ex18.cpp:54
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:153
void SetAllIntPoints(const IntegrationPoint *face_ip)
Set the integration point in the Face and the two neighboring elements, if present.
Definition: eltrans.hpp:553
const double specific_heat_ratio
Definition: ex18.cpp:55
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint, using the method SetIntPoint().
Definition: eltrans.hpp:111
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual void AssembleFaceVector(const FiniteElement &el1, const FiniteElement &el2, FaceElementTransformations &Tr, const Vector &elfun, Vector &elvect)
Perform the local action of the NonlinearFormIntegrator resulting from a face integral term...
Definition: ex18.hpp:383
void ComputeFlux(const Vector &state, int dim, DenseMatrix &flux)
Definition: ex18.hpp:187
int SizeI() const
Definition: densemat.hpp:765
double ComputeMaxCharSpeed(const Vector &state, const int dim)
Definition: ex18.hpp:242
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition: fespace.hpp:460
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:31
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:315
int problem
Definition: ex15.cpp:54
DomainIntegrator(const int dim)
Definition: ex18.hpp:322
const IntegrationPoint & GetElement1IntPoint()
Get a const reference to the integration point in neighboring element 1 corresponding to the currentl...
Definition: eltrans.hpp:564
virtual void Mult(const Vector &x, Vector &y) const
Perform the action of the operator: y = k = f(x, t), where k solves the algebraic equation F(x...
Definition: ex18.hpp:128
int SizeJ() const
Definition: densemat.hpp:766
Class for integration point with weight.
Definition: intrules.hpp:25
This class is used to express the local action of a general nonlinear finite element operator...
Definition: nonlininteg.hpp:26
int dim
Definition: ex24.cpp:53
double ComputePressure(const Vector &state, int dim)
Definition: ex18.hpp:173
void InitialCondition(const Vector &x, Vector &y)
Definition: ex18.hpp:505
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:1798
ElementTransformation * Elem1
Definition: eltrans.hpp:509
virtual ~FE_Evolution()
Definition: ex18.hpp:43
Vector data type.
Definition: vector.hpp:51
bool StateIsPhysical(const Vector &state, const int dim)
Definition: ex18.hpp:458
void UseExternalData(double *d, int h, int w)
Change the data array and the size of the DenseMatrix.
Definition: densemat.hpp:65
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:90
Abstract operator.
Definition: operator.hpp:24
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
const double gas_constant
Definition: ex18.cpp:56
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:728