MFEM  v3.3.2 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  IntegrationPoint eip1;
92  IntegrationPoint eip2;
93
94 public:
95  FaceIntegrator(RiemannSolver &rsolver_, const int dim);
96
97  virtual void AssembleFaceVector(const FiniteElement &el1,
98  const FiniteElement &el2,
100  const Vector &elfun, Vector &elvect);
101 };
102
103 // Implementation of class FE_Evolution
105  Operator &_A, SparseMatrix &_Aflux)
106  : TimeDependentOperator(_A.Height()),
107  dim(_vfes.GetFE(0)->GetDim()),
108  vfes(_vfes),
109  A(_A),
110  Aflux(_Aflux),
111  Me_inv(vfes.GetFE(0)->GetDof(), vfes.GetFE(0)->GetDof(), vfes.GetNE()),
112  state(num_equation),
113  f(num_equation, dim),
114  flux(vfes.GetNDofs(), dim, num_equation),
115  z(A.Height())
116 {
117  // Standard local assembly and inversion for energy mass matrices.
118  const int dof = vfes.GetFE(0)->GetDof();
119  DenseMatrix Me(dof);
120  DenseMatrixInverse inv(&Me);
121  MassIntegrator mi;
122  for (int i = 0; i < vfes.GetNE(); i++)
123  {
124  mi.AssembleElementMatrix(*vfes.GetFE(i), *vfes.GetElementTransformation(i), Me);
125  inv.Factor();
126  inv.GetInverseMatrix(Me_inv(i));
127  }
128 }
129
130 void FE_Evolution::Mult(const Vector &x, Vector &y) const
131 {
132  // 0. Reset wavespeed computation before operator application.
133  max_char_speed = 0.;
134
135  // 1. Create the vector z with the face terms -<F.n(u), [w]>.
136  A.Mult(x, z);
137
138  // 2. Add the element terms.
139  // i. computing the flux approximately as a grid function by interpolating
140  // at the solution nodes.
141  // ii. multiplying this grid function by a (constant) mixed bilinear form for
142  // each of the num_equation, computing (F(u), grad(w)) for each equation.
143
144  DenseMatrix xmat(x.GetData(), vfes.GetNDofs(), num_equation);
145  GetFlux(xmat, flux);
146
147  for (int k = 0; k < num_equation; k++)
148  {
149  Vector fk(flux(k).GetData(), dim * vfes.GetNDofs());
150  Vector zk(z.GetData() + k * vfes.GetNDofs(), vfes.GetNDofs());
152  }
153
154  // 3. Multiply element-wise by the inverse mass matrices.
155  Vector zval;
156  Array<int> vdofs;
157  const int dof = vfes.GetFE(0)->GetDof();
158  DenseMatrix zmat, ymat(dof, num_equation);
159
160  for (int i = 0; i < vfes.GetNE(); i++)
161  {
162  // Return the vdofs ordered byNODES
163  vfes.GetElementVDofs(i, vdofs);
164  z.GetSubVector(vdofs, zval);
165  zmat.UseExternalData(zval.GetData(), dof, num_equation);
166  mfem::Mult(Me_inv(i), zmat, ymat);
167  y.SetSubVector(vdofs, ymat.GetData());
168  }
169 }
170
171 // Physicality check (at end)
172 bool StateIsPhysical(const Vector &state, const int dim);
173
174 // Pressure (EOS) computation
175 inline double ComputePressure(const Vector &state, int dim)
176 {
177  const double den = state(0);
178  const Vector den_vel(state.GetData() + 1, dim);
179  const double den_energy = state(1 + dim);
180
181  double den_vel2 = 0;
182  for (int d = 0; d < dim; d++) { den_vel2 += den_vel(d) * den_vel(d); }
183  den_vel2 /= den;
184
185  return (specific_heat_ratio - 1.0) * (den_energy - 0.5 * den_vel2);
186 }
187
188 // Compute the vector flux F(u)
189 void ComputeFlux(const Vector &state, int dim, DenseMatrix &flux)
190 {
191  const double den = state(0);
192  const Vector den_vel(state.GetData() + 1, dim);
193  const double den_energy = state(1 + dim);
194
195  MFEM_ASSERT(StateIsPhysical(state, dim), "");
196
197  const double pres = ComputePressure(state, dim);
198
199  for (int d = 0; d < dim; d++)
200  {
201  flux(0, d) = den_vel(d);
202  for (int i = 0; i < dim; i++)
203  {
204  flux(1+i, d) = den_vel(i) * den_vel(d) / den;
205  }
206  flux(1+d, d) += pres;
207  }
208
209  const double H = (den_energy + pres) / den;
210  for (int d = 0; d < dim; d++)
211  {
212  flux(1+dim, d) = den_vel(d) * H;
213  }
214 }
215
216 // Compute the scalar F(u).n
217 void ComputeFluxDotN(const Vector &state, const Vector &nor,
218  Vector &fluxN)
219 {
220  // NOTE: nor in general is not a unit normal
221  const int dim = nor.Size();
222  const double den = state(0);
223  const Vector den_vel(state.GetData() + 1, dim);
224  const double den_energy = state(1 + dim);
225
226  MFEM_ASSERT(StateIsPhysical(state, dim), "");
227
228  const double pres = ComputePressure(state, dim);
229
230  double den_velN = 0;
231  for (int d = 0; d < dim; d++) { den_velN += den_vel(d) * nor(d); }
232
233  fluxN(0) = den_velN;
234  for (int d = 0; d < dim; d++)
235  {
236  fluxN(1+d) = den_velN * den_vel(d) / den + pres * nor(d);
237  }
238
239  const double H = (den_energy + pres) / den;
240  fluxN(1 + dim) = den_velN * H;
241 }
242
243 // Compute the maximum characteristic speed.
244 inline double ComputeMaxCharSpeed(const Vector &state, const int dim)
245 {
246  const double den = state(0);
247  const Vector den_vel(state.GetData() + 1, dim);
248
249  double den_vel2 = 0;
250  for (int d = 0; d < dim; d++) { den_vel2 += den_vel(d) * den_vel(d); }
251  den_vel2 /= den;
252
253  const double pres = ComputePressure(state, dim);
254  const double sound = sqrt(specific_heat_ratio * pres / den);
255  const double vel = sqrt(den_vel2 / den);
256
257  return vel + sound;
258 }
259
260 // Compute the flux at solution nodes.
261 void FE_Evolution::GetFlux(const DenseMatrix &x, DenseTensor &flux) const
262 {
263  const int dof = flux.SizeI();
264  const int dim = flux.SizeJ();
265
266  for (int i = 0; i < dof; i++)
267  {
268  for (int k = 0; k < num_equation; k++) { state(k) = x(i, k); }
269  ComputeFlux(state, dim, f);
270
271  for (int d = 0; d < dim; d++)
272  {
273  for (int k = 0; k < num_equation; k++)
274  {
275  flux(i, d, k) = f(k, d);
276  }
277  }
278
279  // Update max char speed
280  const double mcs = ComputeMaxCharSpeed(state, dim);
281  if (mcs > max_char_speed) { max_char_speed = mcs; }
282  }
283 }
284
285 // Implementation of class RiemannSolver
287  flux1(num_equation),
288  flux2(num_equation) { }
289
290 double RiemannSolver::Eval(const Vector &state1, const Vector &state2,
291  const Vector &nor, Vector &flux)
292 {
293  // NOTE: nor in general is not a unit normal
294  const int dim = nor.Size();
295
296  MFEM_ASSERT(StateIsPhysical(state1, dim), "");
297  MFEM_ASSERT(StateIsPhysical(state2, dim), "");
298
299  const double maxE1 = ComputeMaxCharSpeed(state1, dim);
300  const double maxE2 = ComputeMaxCharSpeed(state2, dim);
301
302  const double maxE = max(maxE1, maxE2);
303
304  ComputeFluxDotN(state1, nor, flux1);
305  ComputeFluxDotN(state2, nor, flux2);
306
307  double normag = 0;
308  for (int i = 0; i < dim; i++)
309  {
310  normag += nor(i) * nor(i);
311  }
312  normag = sqrt(normag);
313
314  for (int i = 0; i < num_equation; i++)
315  {
316  flux(i) = 0.5 * (flux1(i) + flux2(i))
317  - 0.5 * maxE * (state2(i) - state1(i)) * normag;
318  }
319
320  return maxE;
321 }
322
323 // Implementation of class DomainIntegrator
324 DomainIntegrator::DomainIntegrator(const int dim) : flux(num_equation, dim) { }
325
327  const FiniteElement &test_fe,
329  DenseMatrix &elmat)
330 {
331  // Assemble the form (vec(v), grad(w))
332
333  // Trial space = vector L2 space (mesh dim)
334  // Test space = scalar L2 space
335
336  const int dof_trial = trial_fe.GetDof();
337  const int dof_test = test_fe.GetDof();
338  const int dim = trial_fe.GetDim();
339
340  shape.SetSize(dof_trial);
341  dshapedr.SetSize(dof_test, dim);
342  dshapedx.SetSize(dof_test, dim);
343
344  elmat.SetSize(dof_test, dof_trial * dim);
345  elmat = 0.0;
346
347  const int maxorder = max(trial_fe.GetOrder(), test_fe.GetOrder());
348  const int intorder = 2 * maxorder;
349  const IntegrationRule *ir = &IntRules.Get(trial_fe.GetGeomType(), intorder);
350
351  for (int i = 0; i < ir->GetNPoints(); i++)
352  {
353  const IntegrationPoint &ip = ir->IntPoint(i);
354
355  // Calculate the shape functions
356  trial_fe.CalcShape(ip, shape);
357  shape *= ip.weight;
358
359  // Compute the physical gradients of the test functions
360  Tr.SetIntPoint(&ip);
361  test_fe.CalcDShape(ip, dshapedr);
363
364  for (int d = 0; d < dim; d++)
365  {
366  for (int j = 0; j < dof_test; j++)
367  {
368  for (int k = 0; k < dof_trial; k++)
369  {
370  elmat(j, k + d * dof_trial) += shape(k) * dshapedx(j, d);
371  }
372  }
373  }
374  }
375 }
376
377 // Implementation of class FaceIntegrator
379  rsolver(rsolver_),
380  funval1(num_equation),
381  funval2(num_equation),
382  nor(dim),
383  fluxN(num_equation) { }
384
386  const FiniteElement &el2,
388  const Vector &elfun, Vector &elvect)
389 {
390  // Compute the term <F.n(u),[w]> on the interior faces.
391  const int dof1 = el1.GetDof();
392  const int dof2 = el2.GetDof();
393
394  shape1.SetSize(dof1);
395  shape2.SetSize(dof2);
396
397  elvect.SetSize((dof1 + dof2) * num_equation);
398  elvect = 0.0;
399
400  DenseMatrix elfun1_mat(elfun.GetData(), dof1, num_equation);
401  DenseMatrix elfun2_mat(elfun.GetData() + dof1 * num_equation, dof2,
402  num_equation);
403
404  DenseMatrix elvect1_mat(elvect.GetData(), dof1, num_equation);
405  DenseMatrix elvect2_mat(elvect.GetData() + dof1 * num_equation, dof2,
406  num_equation);
407
408  // Integration order calculation from DGTraceIntegrator
409  int intorder;
410  if (Tr.Elem2No >= 0)
411  intorder = (min(Tr.Elem1->OrderW(), Tr.Elem2->OrderW()) +
412  2*max(el1.GetOrder(), el2.GetOrder()));
413  else
414  {
415  intorder = Tr.Elem1->OrderW() + 2*el1.GetOrder();
416  }
417  if (el1.Space() == FunctionSpace::Pk)
418  {
419  intorder++;
420  }
421  const IntegrationRule *ir = &IntRules.Get(Tr.FaceGeom, intorder);
422
423  for (int i = 0; i < ir->GetNPoints(); i++)
424  {
425  const IntegrationPoint &ip = ir->IntPoint(i);
426
427  Tr.Loc1.Transform(ip, eip1);
428  Tr.Loc2.Transform(ip, eip2);
429
430  // Calculate basis functions on both elements at the face
431  el1.CalcShape(eip1, shape1);
432  el2.CalcShape(eip2, shape2);
433
434  // Interpolate elfun at the point
435  elfun1_mat.MultTranspose(shape1, funval1);
436  elfun2_mat.MultTranspose(shape2, funval2);
437
438  Tr.Face->SetIntPoint(&ip);
439
440  // Get the normal vector and the flux on the face
441  CalcOrtho(Tr.Face->Jacobian(), nor);
442  const double mcs = rsolver.Eval(funval1, funval2, nor, fluxN);
443
444  // Update max char speed
445  if (mcs > max_char_speed) { max_char_speed = mcs; }
446
447  fluxN *= ip.weight;
448  for (int k = 0; k < num_equation; k++)
449  {
450  for (int s = 0; s < dof1; s++)
451  {
452  elvect1_mat(s, k) -= fluxN(k) * shape1(s);
453  }
454  for (int s = 0; s < dof2; s++)
455  {
456  elvect2_mat(s, k) += fluxN(k) * shape2(s);
457  }
458  }
459  }
460 }
461
462 // Check that the state is physical - enabled in debug mode
463 bool StateIsPhysical(const Vector &state, const int dim)
464 {
465  const double den = state(0);
466  const Vector den_vel(state.GetData() + 1, dim);
467  const double den_energy = state(1 + dim);
468
469  if (den < 0)
470  {
471  cout << "Negative density: ";
472  for (int i = 0; i < state.Size(); i++)
473  {
474  cout << state(i) << " ";
475  }
476  cout << endl;
477  return false;
478  }
479  if (den_energy <= 0)
480  {
481  cout << "Negative energy: ";
482  for (int i = 0; i < state.Size(); i++)
483  {
484  cout << state(i) << " ";
485  }
486  cout << endl;
487  return false;
488  }
489
490  double den_vel2 = 0;
491  for (int i = 0; i < dim; i++) { den_vel2 += den_vel(i) * den_vel(i); }
492  den_vel2 /= den;
493
494  const double pres = (specific_heat_ratio - 1.0) * (den_energy - 0.5 * den_vel2);
495
496  if (pres <= 0)
497  {
498  cout << "Negative pressure: " << pres << ", state: ";
499  for (int i = 0; i < state.Size(); i++)
500  {
501  cout << state(i) << " ";
502  }
503  cout << endl;
504  return false;
505  }
506  return true;
507 }
508
509 // Initial condition
510 void InitialCondition(const Vector &x, Vector &y)
511 {
512  const int dim = x.Size();
513  MFEM_ASSERT(dim == 2, "");
514
515  double radius = 0, Minf = 0, beta = 0;
516  if (problem == 1)
517  {
518  // "Fast vortex"
520  Minf = 0.5;
521  beta = 1. / 5.;
522  }
523  else if (problem == 2)
524  {
525  // "Slow vortex"
527  Minf = 0.05;
528  beta = 1. / 50.;
529  }
530  else
531  {
532  mfem_error("Cannot recognize problem."
533  "Options are: 1 - fast vortex, 2 - slow vortex");
534  }
535
536  const double xc = 0.0, yc = 0.0;
537
538  // Nice units
539  const double vel_inf = 1.;
540  const double den_inf = 1.;
541
542  // Derive remainder of background state from this and Minf
543  const double pres_inf = (den_inf / specific_heat_ratio) * (vel_inf / Minf) *
544  (vel_inf / Minf);
545  const double temp_inf = pres_inf / (den_inf * gas_constant);
546
548  r2rad += (x(0) - xc) * (x(0) - xc);
549  r2rad += (x(1) - yc) * (x(1) - yc);
551
552  const double shrinv1 = 1.0 / (specific_heat_ratio - 1.);
553
554  const double velX = vel_inf * (1 - beta * (x(1) - yc) / radius * exp(
556  const double velY = vel_inf * beta * (x(0) - xc) / radius * exp(-0.5 * r2rad);
557  const double vel2 = velX * velX + velY * velY;
558
559  const double specific_heat = gas_constant * specific_heat_ratio * shrinv1;
560  const double temp = temp_inf - 0.5 * (vel_inf * beta) *
561  (vel_inf * beta) / specific_heat * exp(-r2rad);
562
563  const double den = den_inf * pow(temp/temp_inf, shrinv1);
564  const double pres = den * gas_constant * temp;
565  const double energy = shrinv1 * pres / den + 0.5 * vel2;
566
567  y(0) = den;
568  y(1) = den * velX;
569  y(2) = den * velY;
570  y(3) = den * energy;
571 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:222
Abstract class for Finite Elements.
Definition: fe.hpp:47
Definition: eltrans.hpp:72
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:498
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe.hpp:116
ElementTransformation * Face
Definition: eltrans.hpp:347
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:161
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:83
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:861
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:320
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:133
Base abstract class for time dependent operators.
Definition: operator.hpp:142
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:468
void SetIntPoint(const IntegrationPoint *ip)
Definition: eltrans.hpp:52
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Definition: vector.cpp:462
int GetOrder() const
Returns the order of the finite element.
Definition: fe.hpp:125
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:481
int Size() const
Returns the size of the vector.
Definition: vector.hpp:113
FaceIntegrator(RiemannSolver &rsolver_, const int dim)
Definition: ex18.hpp:378
IntegrationPointTransformation Loc2
Definition: eltrans.hpp:348
int Space() const
Returns the type of space on each element.
Definition: fe.hpp:128
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
void GetInverseMatrix(DenseMatrix &Ainv) const
Compute and return the inverse matrix in Ainv.
Definition: densemat.hpp:562
void Factor()
Factor the current DenseMatrix, *a.
Definition: densemat.cpp:4142
double * GetData() const
Returns the matrix data array.
Definition: densemat.hpp:94
ElementTransformation * Elem2
Definition: eltrans.hpp:347
double * GetData() const
Definition: vector.hpp:121
void CalcOrtho(const DenseMatrix &J, Vector &n)
Definition: densemat.cpp:3278
double max_char_speed
Definition: ex18.cpp:59
RiemannSolver()
Definition: ex18.hpp:286
FE_Evolution(FiniteElementSpace &_vfes, Operator &_A, SparseMatrix &_Aflux)
Definition: ex18.hpp:104
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:184
double Eval(const Vector &state1, const Vector &state2, const Vector &nor, Vector &flux)
Definition: ex18.hpp:290
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:225
void ComputeFluxDotN(const Vector &state, const Vector &nor, Vector &fluxN)
Definition: ex18.hpp:217
virtual void AssembleElementMatrix(const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
Definition: bilininteg.cpp:721
virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &Tr, DenseMatrix &elmat)
Definition: ex18.hpp:326
int dim
Definition: ex3.cpp:47
Data type sparse matrix.
Definition: sparsemat.hpp:38
IntegrationPointTransformation Loc1
Definition: eltrans.hpp:348
const int num_equation
Definition: ex18.cpp:54
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:67
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:385
void ComputeFlux(const Vector &state, int dim, DenseMatrix &flux)
Definition: ex18.hpp:189
int SizeI() const
Definition: densemat.hpp:658
double ComputeMaxCharSpeed(const Vector &state, const int dim)
Definition: ex18.hpp:244
int GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe.hpp:119
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i&#39;th element.
Definition: fespace.hpp:205
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:22
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:122
int problem
Definition: ex15.cpp:57
void mfem_error(const char *msg)
Definition: error.cpp:107
DomainIntegrator(const int dim)
Definition: ex18.hpp:324
void Transform(const IntegrationPoint &, IntegrationPoint &)
Definition: eltrans.cpp:515
virtual ~FE_Evolution()
Definition: ex18.hpp:43
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:130
int SizeJ() const
Definition: densemat.hpp:659
Class for integration point with weight.
Definition: intrules.hpp:25
double ComputePressure(const Vector &state, int dim)
Definition: ex18.hpp:175
void InitialCondition(const Vector &x, Vector &y)
Definition: ex18.hpp:510
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement associated with i&#39;th element.
Definition: fespace.cpp:1171
ElementTransformation * Elem1
Definition: eltrans.hpp:347
Vector data type.
Definition: vector.hpp:41
bool StateIsPhysical(const Vector &state, const int dim)
Definition: ex18.hpp:463
void UseExternalData(double *d, int h, int w)
Change the data array and the size of the DenseMatrix.
Definition: densemat.hpp:64
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:86
Abstract operator.
Definition: operator.hpp:21
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:634