NeoPZ
pzeulerconslaw.h
Go to the documentation of this file.
1 
6 #ifndef TPZEULERCONSLAW_H
7 #define TPZEULERCONSLAW_H
8 
9 #include <iostream>
10 #include "TPZMaterial.h"
11 #include "tpzoutofrange.h"
12 #include "pzfmatrix.h"
13 #include "pzvec.h"
14 #include "pzconslaw.h"
15 #include "pzartdiff.h"
16 
17 #include "pzlog.h"
18 
19 #ifdef LOG4CXX
20 extern LoggerPtr fluxroe;
21 extern LoggerPtr fluxappr;
22 #endif
23 
24 #ifdef _AUTODIFF
25 #include "fadType.h"
26 #define _TFAD
27 #endif
28 
33 class TPZEulerConsLaw : public TPZConservationLaw
34 {
35  public :
36 
37  TPZEulerConsLaw(int nummat,REAL timeStep,
38  REAL gamma,int dim,
39  TPZArtDiffType artdiff);
40 
41  ~TPZEulerConsLaw();
42 
43  TPZEulerConsLaw();
44 
45  TPZEulerConsLaw(const TPZEulerConsLaw &cp) :
46  TPZRegisterClassId(&TPZEulerConsLaw::ClassId),TPZConservationLaw(cp), fArtDiff(cp.fArtDiff),fDiff(cp.fDiff),
47  fConvVol(cp.fConvVol),fConvFace(cp.fConvFace)
48  {
49  }
50 
51  virtual TPZMaterial * NewMaterial() override
52  {
53  return new TPZEulerConsLaw(*this);
54  }
55 
62  void SetTimeDiscr(TPZTimeDiscr Diff, TPZTimeDiscr ConvVol, TPZTimeDiscr ConvFace);
63 
65  TPZArtDiff & ArtDiff(){ return fArtDiff; };
66 
68  REAL OptimalCFL(int degree);
69 
71  virtual REAL SetTimeStep(REAL maxveloc,REAL deltax,int degree) override;
72 
74  static int NStateVariables(int dim);
75 
77  virtual int NStateVariables() const override;
78 
85  REAL DeltaX(REAL detJac);
86 
91  REAL Det(TPZFMatrix<REAL> & Mat);
92 
100  template< class T >
101  static void Pressure(REAL gamma, int dim, T& press, TPZVec<T> &U);
102 
104  template <class T>
105  static void cSpeed(TPZVec<T> & sol, REAL gamma, T & c);
106 
108  template <class T>
109  static void uRes(TPZVec<T> & sol, T & us);
110 
111  virtual STATE Pressure(TPZVec<STATE> &U) override;
112 
113  virtual void Print(std::ostream & out) override;
114 
115  virtual std::string Name() override;
116 
117  virtual int VariableIndex(const std::string &name) override;
118 
119  virtual int NSolutionVariables(int var) override;
120 
121  virtual int NFluxes() override;
122 
123 
128  template <class T>
129  void ComputeGhostState(TPZVec<T> &solL, TPZVec<T> &solR, TPZVec<REAL> &normal, TPZBndCond &bc, int & entropyFix);
130 
131 protected:
132  virtual void Solution(TPZVec<STATE> &Sol,TPZFMatrix<STATE> &DSol,TPZFMatrix<REAL> &axes,int var,TPZVec<STATE> &Solout) override;
133 public:
134 
136  virtual void Solution(TPZMaterialData &data, int var, TPZVec<STATE> &Solout) override
137  {
138  TPZConservationLaw::Solution(data,var,Solout);
139  }
146  template < class T >
147  void Flux(TPZVec<T> &U,TPZVec<T> &Fx,TPZVec<T> &Fy,TPZVec<T> &Fz);
148 
158  template <class T>
159  static void JacobFlux(REAL gamma, int dim, TPZVec<T> & U,TPZVec<TPZDiffMatrix<T> > &Ai);
160 
162  template <class T>
163  void Test_Flux(TPZVec<T> &solL, TPZVec<T> &solR, TPZVec<REAL> & normal, REAL gamma, TPZVec<T> & flux);
164 
175  template <class T>
176  static void Roe_Flux(TPZVec<T> &solL, TPZVec<T> &solR,
177  TPZVec<REAL> & normal, REAL gamma,
178  TPZVec<T> & flux, int entropyFix = 1);
179 
190  template <class T>
191  static void ApproxRoe_Flux(TPZVec<T> &solL, TPZVec<T> &solR,
192  TPZVec<REAL> & normal, REAL gamma,
193  TPZVec<T> & flux, int entropyFix = 1);
194 
195 
197  void SetDelta(REAL delta);
198 
199 public:
201  template <class T>
202  static void Roe_Flux(const T & rho_f,
203  const T & rhou_f,
204  const T & rhov_f,
205  const T & rhow_f,
206  const T & rhoE_f,
207  const T & rho_t,
208  const T & rhou_t,
209  const T & rhov_t,
210  const T & rhow_t,
211  const T & rhoE_t,
212  const REAL nx,
213  const REAL ny,
214  const REAL nz,
215  const REAL gam,
216  T & flux_rho,
217  T & flux_rhou,
218  T & flux_rhov,
219  T & flux_rhow,
220  T & flux_rhoE, int entropyFix = 1);
221 
222  template <class T>
223  static void Roe_Flux(const T & rho_f,
224  const T & rhou_f,
225  const T & rhov_f,
226  const T & rhoE_f,
227  const T & rho_t,
228  const T & rhou_t,
229  const T & rhov_t,
230  const T & rhoE_t,
231  const REAL nx,
232  const REAL ny,
233  const REAL gam,
234  T &flux_rho,
235  T &flux_rhou,
236  T &flux_rhov,
237  T &flux_rhoE, int entropyFix = 1);
239  template <class T>
240  static void ApproxRoe_Flux(const T & rho_f,
241  const T & rhou_f,
242  const T & rhov_f,
243  const T & rhow_f,
244  const T & rhoE_f,
245  const T & rho_t,
246  const T & rhou_t,
247  const T & rhov_t,
248  const T & rhow_t,
249  const T & rhoE_t,
250  const REAL nx,
251  const REAL ny,
252  const REAL nz,
253  const REAL gam,
254  T & flux_rho,
255  T & flux_rhou,
256  T & flux_rhov,
257  T & flux_rhow,
258  T & flux_rhoE, int entropyFix = 1);
259 
260  template <class T>
261  static void ApproxRoe_Flux(const T & rho_f,
262  const T & rhou_f,
263  const T & rhov_f,
264  const T & rhoE_f,
265  const T & rho_t,
266  const T & rhou_t,
267  const T & rhov_t,
268  const T & rhoE_t,
269  const REAL nx,
270  const REAL ny,
271  const REAL gam,
272  T &flux_rho,
273  T &flux_rhou,
274  T &flux_rhov,
275  T &flux_rhoE, int entropyFix = 1);
278 #ifdef _AUTODIFF
279 
291  void PrepareFAD(TPZVec<STATE> & sol, TPZFMatrix<STATE> & dsol,
292  TPZFMatrix<REAL> & phi, TPZFMatrix<REAL> & dphi,
293  TPZVec<FADREAL> & FADsol,
294  TPZVec<FADREAL> & FADdsol);
295 
306  void PrepareInterfaceFAD(
307  TPZVec<STATE> &solL,TPZVec<STATE> &solR,
308  TPZFMatrix<REAL> &phiL,TPZFMatrix<REAL> &phiR,
309  TPZVec<FADREAL> & FADsolL,
310  TPZVec<FADREAL> & FADsolR);
311 
320  template <class T>
321  void PrepareFastestInterfaceFAD(
322  TPZVec<STATE> &solL,TPZVec<STATE> &solR,
323  TPZVec<T> & FADsolL,
324  TPZVec<T> & FADsolR);
327 #endif
328 
332  virtual void Contribute(TPZMaterialData &data, REAL weight, TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef) override;
333 
334  virtual void Contribute(TPZMaterialData &data, REAL weight, TPZFMatrix<STATE> &ef) override;
335 
336  virtual void ContributeInterface(TPZMaterialData &data, TPZMaterialData &dataleft, TPZMaterialData &dataright, REAL weight, TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef) override;
337 
338  virtual void ContributeInterface(TPZMaterialData &data, TPZMaterialData &dataleft, TPZMaterialData &dataright, REAL weight, TPZFMatrix<STATE> &ef) override;
339 
340  virtual void ContributeBC(TPZMaterialData &data,
341  REAL weight,
342  TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef,
343  TPZBndCond &bc) override;
344 
345  virtual void ContributeBC(TPZMaterialData &data,
346  REAL weight,
347  TPZFMatrix<STATE> &ef,
348  TPZBndCond &bc) override
349  {
350  TPZDiscontinuousGalerkin::ContributeBC(data,weight,ef,bc);
351  }
352 
353  virtual void ContributeBCInterface(TPZMaterialData &data, TPZMaterialData &dataleft,
354  REAL weight,
355  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,
356  TPZBndCond &bc) override;
357 
358  virtual void ContributeBCInterface(TPZMaterialData &data, TPZMaterialData &dataleft,
359  REAL weight,
360  TPZFMatrix<STATE> &ef,
361  TPZBndCond &bc) override;
362 
366  void ContributeFastestBCInterface(int dim,
367  TPZVec<REAL> &x,
368  TPZVec<STATE> &solL, TPZFMatrix<STATE> &dsolL,
369  REAL weight, TPZVec<REAL> &normal,
370  TPZFMatrix<REAL> &phiL, TPZFMatrix<REAL> &dphiL,
371  TPZFMatrix<REAL> &axesleft,
372  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef, TPZBndCond &bc);
373 
374  template <int dim>
375  void ContributeFastestBCInterface_dim(
376  TPZVec<REAL> &x,
377  TPZVec<STATE> &solL, TPZFMatrix<STATE> &dsolL,
378  REAL weight, TPZVec<REAL> &normal,
379  TPZFMatrix<REAL> &phiL,TPZFMatrix<REAL> &dphiL,
380  TPZFMatrix<REAL> &axesleft,
381  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc);
382 
383  virtual void ContributeLast(TPZVec<REAL> &x,TPZFMatrix<REAL> &jacinv,
384  TPZVec<STATE> &sol,TPZFMatrix<STATE> &dsol,
385  REAL weight,
386  TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi,
387  TPZFMatrix<STATE> &ef);
388 
389  virtual void ContributeAdv(TPZVec<REAL> &x,TPZFMatrix<REAL> &jacinv,
390  TPZVec<STATE> &sol,TPZFMatrix<STATE> &dsol,
391  REAL weight,
392  TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi,
393  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef);
394 
395  virtual void ContributeAdv(TPZVec<REAL> &x,TPZFMatrix<REAL> &jacinv,
396  TPZVec<STATE> &sol,TPZFMatrix<STATE> &dsol,
397  REAL weight,
398  TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi,
399  TPZFMatrix<STATE> &ef);
400 
401  void ContributeApproxImplDiff(TPZVec<REAL> &x,
402  TPZFMatrix<REAL> &jacinv,
403  TPZVec<STATE> &sol,TPZFMatrix<STATE> &dsol,
404  REAL weight,
405  TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi,
406  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef);
407 
408  void ContributeExplDiff(TPZVec<REAL> &x,
409  TPZFMatrix<REAL> &jacinv,
410  TPZVec<STATE> &sol,TPZFMatrix<STATE> &dsol,
411  REAL weight,
412  TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi,
413  TPZFMatrix<STATE> &ef);
414 
415 #ifdef _AUTODIFF
416  void ContributeImplDiff(TPZVec<REAL> &x,
417  TPZFMatrix<REAL> &jacinv,
418  TPZVec<FADREAL> &sol,TPZVec<FADREAL> &dsol,
419  REAL weight,
420  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef);
421 
422  void ContributeFastestImplDiff(int dim,
423  TPZVec<REAL> &x,
424  TPZFMatrix<REAL> &jacinv,
425  TPZVec<STATE> &sol,TPZFMatrix<STATE> &dsol,
426  TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi,
427  REAL weight,
428  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef);
429 
430 #endif
431 
432  void ContributeExplConvFace(TPZVec<REAL> &x,
433  TPZVec<STATE> &solL,TPZVec<STATE> &solR,
434  REAL weight,TPZVec<REAL> &normal,
435  TPZFMatrix<REAL> &phiL,TPZFMatrix<REAL> &phiR,
436  TPZFMatrix<STATE> &ef, int entropyFix = 1);
437 
438  void ContributeApproxImplConvFace(TPZVec<REAL> &x, REAL faceSize,
439  TPZVec<STATE> &solL,TPZVec<STATE> &solR,
440  REAL weight,TPZVec<REAL> &normal,
441  TPZFMatrix<REAL> &phiL,TPZFMatrix<REAL> &phiR,
442  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef, int entropyFix = 1
443  );
444 
445 #ifdef _AUTODIFF
446 
447  void ContributeApproxImplConvFace(TPZVec<REAL> &x, REAL faceSize,
448  TPZVec<FADREAL> &solL,TPZVec<FADREAL> &solR,
449  REAL weight,TPZVec<REAL> &normal,
450  TPZFMatrix<REAL> &phiL,TPZFMatrix<REAL> &phiR,
451  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef, int entropyFix = 1
452  );
453 
454  void ContributeImplConvFace(TPZVec<REAL> &x,
455  TPZVec<FADREAL> &solL,TPZVec<FADREAL> &solR,
456  REAL weight,TPZVec<REAL> &normal,
457  TPZFMatrix<REAL> &phiL,TPZFMatrix<REAL> &phiR,
458  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef, int entropyFix = 1);
459 
460  void ContributeFastestImplConvFace(int dim,
461  TPZVec<REAL> &x,
462  TPZVec<STATE> &solL,TPZVec<STATE> &solR,
463  REAL weight,TPZVec<REAL> &normal,
464  TPZFMatrix<REAL> &phiL,TPZFMatrix<REAL> &phiR,
465  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef, int entropyFix = 1);
466 
467  template <int dim>
468  void ContributeFastestImplConvFace_dim(
469  TPZVec<REAL> &x,
470  TPZVec<STATE> &solL,TPZVec<STATE> &solR,
471  REAL weight,TPZVec<REAL> &normal,
472  TPZFMatrix<REAL> &phiL,TPZFMatrix<REAL> &phiR,
473  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef, int entropyFix = 1);
474 
475  template <class T>
476  void ContributeFastestImplConvFace_T(TPZVec<REAL> &x,
477  TPZVec<T> &FADsolL,TPZVec<T> &FADsolR,
478  REAL weight,TPZVec<REAL> &normal,
479  TPZFMatrix<REAL> &phiL,TPZFMatrix<REAL> &phiR,
480  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,int entropyFix = 1);
481 
482 #endif
483 
484  void ContributeImplConvVol(TPZVec<REAL> &x,
485  TPZVec<STATE> &sol,TPZFMatrix<STATE> &dsol,
486  REAL weight,
487  TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi,
488  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef);
489 
490  void ContributeExplConvVol(TPZVec<REAL> &x,
491  TPZVec<STATE> &sol,
492  REAL weight,
493  TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi,
494  TPZFMatrix<STATE> &ef);
495 
496 #ifdef _AUTODIFF
497  void ContributeImplConvVol(TPZVec<REAL> &x,
498  TPZVec<FADREAL> &sol,TPZVec<FADREAL> &dsol,
499  REAL weight,
500  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef);
501 #endif
502 
503  void ContributeExplT1(TPZVec<REAL> &x,
504  TPZVec<STATE> &sol,TPZFMatrix<STATE> &dsol,
505  REAL weight,
506  TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi,
507  TPZFMatrix<STATE> &ef);
508 
509  void ContributeImplT1(TPZVec<REAL> &x,
510  TPZVec<STATE> &sol,TPZFMatrix<STATE> &dsol,
511  REAL weight,
512  TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi,
513  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef);
514 
515  void ContributeExplT2(TPZVec<REAL> &x,
516  TPZVec<STATE> &sol,
517  REAL weight,
518  TPZFMatrix<REAL> &phi,
519  TPZFMatrix<STATE> &ef);
525  virtual void Write(TPZStream &buf, int withclassid) const override;
526 
528  void Read(TPZStream &buf, void *context) override;
529 
531 public:
532 virtual int ClassId() const override;
533 
534 
538 protected:
539 
541  TPZArtDiff fArtDiff;
542 
544  TPZTimeDiscr fDiff, fConvVol, fConvFace;
545 
547 };
548 
549 
550 inline std::string TPZEulerConsLaw::Name()
551 {
552  return "TPZEulerConsLaw";
553 }
554 
555 template < class T >
556 inline void TPZEulerConsLaw::Flux(TPZVec<T> &U,TPZVec<T> &Fx,TPZVec<T> &Fy,TPZVec<T> &Fz) {
557 
558  T press;
559  Pressure(fGamma, fDim, press, U);
560  int nstate = NStateVariables();
561  if(nstate < 3 && nstate > 5){
562  PZError << "TPZEulerConsLaw::Flux case not implemented\n";
563  Fx.Resize(0);
564  Fy.Resize(0);
565  Fz.Resize(0);
566  exit(-1);
567  }
568 
569  Fx.Resize(5,0.0);//v�ido
570  Fy.Resize(5,0.0);//para
571  Fz.Resize(5,0.0);//R , R , R
572 
573  if(nstate == 5){
574  Fx[0] = U[1];//ro u
575  Fx[1] = (U[1]/U[0])*U[1] + press;//ro u2 + p
576  Fx[2] = U[1]*(U[2]/U[0]);//ro u v
577  Fx[3] = U[1]*(U[3]/U[0]);//ro u w
578  Fx[4] = (U[4]+press)*(U[1]/U[0]);//(ro e + p) u
579 
580  Fy[0] = U[2];//ro v
581  Fy[1] = U[2]*(U[1]/U[0]);//ro v u
582  Fy[2] = (U[2]/U[0])*U[2] + press;//ro v2 + p
583  Fy[3] = U[2]*(U[3]/U[0]);//ro v w
584  Fy[4] = (U[4] + press)*(U[2]/U[0]);//(ro e + p) v
585 
586  Fz[0] = U[3];//ro w
587  Fz[1] = U[3]*(U[1]/U[0]);//ro w u
588  Fz[2] = U[3]*(U[2]/U[0]);//ro w v
589  Fz[3] = (U[3]/U[0])*U[3] + press;//ro w2 + p
590  Fz[4] = (U[4] + press)*(U[3]/U[0]);//(ro e + p) w
591  return;
592  }
593 
594  if(nstate == 4){
595  Fx[0] = U[1];//ro u
596  Fx[1] = U[1]*U[1] / U[0] + press;//ro u2 + p
597  Fx[2] = U[1]*U[2] / U[0];//ro u v
598  Fx[3] = (U[3]+press)*U[1] / U[0];//(E + p) u
599 
600  Fy[0] = U[2];//ro v
601  Fy[1] = U[1]*U[2] / U[0];//ro u v
602  Fy[2] = U[2]*U[2] / U[0] + press;//ro v2 + p
603  Fy[3] = (U[3] + press)*U[2] / U[0];//(E + p) v
604  return;
605  }
606 
607  if(nstate == 3){
608  Fx[0] = U[1];//ro u
609  Fx[1] = (U[1]/U[0])*U[1] + press;//ro u2 + p
610  Fx[2] = (U[2]+press)*(U[1]/U[0]);//(ro e + p) u
611  }
612 }
613 
614 template <class T>
615 inline void TPZEulerConsLaw::JacobFlux(REAL gamma, int dim, TPZVec<T> & U,TPZVec<TPZDiffMatrix<T> > &Ai)
616 {
617 
618  Ai.Resize(dim);
619  int i;
620  for(i=0;i<dim;i++)Ai[i].Redim(TPZEulerConsLaw::NStateVariables(dim), TPZEulerConsLaw::NStateVariables(dim));
621 
622  if(U[0] < REAL(1.e-6)) {
623  PZError << "TPZEulerConsLaw::JacobFlux: Density negative " << U[0] << std::endl;
624  TPZOutofRange obj;
625  throw(obj);
626  }
627 
628  T u,v,w,e;
629  REAL gamma1 = gamma-1.;
630  REAL gamma2 = gamma1/2.;
631  REAL gamma3 = gamma-3;
632 
633  if(dim == 3){
634 
635  u = U[1]/U[0];
636  v = U[2]/U[0];
637  w = U[3]/U[0];
638  e = U[4]/U[0];
639 
640  T u2 = u*u;
641  T v2 = v*v;
642  T w2 = w*w;
643  T vel = u2+v2+w2;
644 
645  Ai[0](0,0) = 0.;
646  Ai[0](0,1) = 1.;
647  Ai[0](0,2) = 0.;
648  Ai[0](0,3) = 0.;
649  Ai[0](0,4) = 0.;
650 
651  Ai[0](1,0) = gamma2*vel-u2;
652  Ai[0](1,1) = -gamma3*u;
653  Ai[0](1,2) = -gamma1*v;
654  Ai[0](1,3) = -gamma1*w;
655  Ai[0](1,4) = gamma1;
656 
657  Ai[0](2,0) = -u*v;
658  Ai[0](2,1) = v;
659  Ai[0](2,2) = u;
660  Ai[0](2,3) = 0.;
661  Ai[0](2,4) = 0.;
662 
663  Ai[0](3,0) = -u*w;
664  Ai[0](3,1) = w;
665  Ai[0](3,2) = 0.;
666  Ai[0](3,3) = u;
667  Ai[0](3,4) = 0.;
668 
669  Ai[0](4,0) = -gamma*e*u + gamma1*u*vel;
670  Ai[0](4,1) = gamma*e - gamma1*u2 - gamma2*vel;
671  Ai[0](4,2) = -gamma1*u*v;
672  Ai[0](4,3) = -gamma1*u*w;
673  Ai[0](4,4) = gamma*u;
674 
675  Ai[1](0,0) = 0.;
676  Ai[1](0,1) = 0.;
677  Ai[1](0,2) = 1.;
678  Ai[1](0,3) = 0.;
679  Ai[1](0,4) = 0.;
680 
681  Ai[1](1,0) = -u*v;
682  Ai[1](1,1) = v;
683  Ai[1](1,2) = u;
684  Ai[1](1,3) = 0.;
685  Ai[1](1,4) = 0.;
686 
687  Ai[1](2,0) = gamma2*vel-v2;
688  Ai[1](2,1) = -gamma1*u;
689  Ai[1](2,2) = -gamma3*v;
690  Ai[1](2,3) = -gamma1*w;
691  Ai[1](2,4) = gamma1;
692 
693  Ai[1](3,0) = -v*w;
694  Ai[1](3,1) = 0.;
695  Ai[1](3,2) = w;
696  Ai[1](3,3) = v;
697  Ai[1](3,4) = 0.;
698 
699  Ai[1](4,0) = -gamma*e*v + gamma1*v*vel;
700  Ai[1](4,1) = -gamma1*u*v;
701  Ai[1](4,2) = gamma*e - gamma1*v2 - gamma2*vel;
702  Ai[1](4,3) = -gamma1*v*w;
703  Ai[1](4,4) = gamma*v;
704 
705  Ai[2](0,0) = 0.;
706  Ai[2](0,1) = 0.;
707  Ai[2](0,2) = 0.;
708  Ai[2](0,3) = 1.;
709  Ai[2](0,4) = 0.;
710 
711  Ai[2](1,0) = -u*w;
712  Ai[2](1,1) = w;
713  Ai[2](1,2) = 0.;
714  Ai[2](1,3) = u;
715  Ai[2](1,4) = 0.;
716 
717  Ai[2](2,0) = -v*w;
718  Ai[2](2,1) = 0.;
719  Ai[2](2,2) = w;
720  Ai[2](2,3) = v;
721  Ai[2](2,4) = 0.;
722 
723  Ai[2](3,0) = gamma2*vel-w2;
724  Ai[2](3,1) = -gamma1*u;
725  Ai[2](3,2) = -gamma1*v;
726  Ai[2](3,3) = -gamma3*w;
727  Ai[2](3,4) = gamma1;
728 
729  Ai[2](4,0) = -gamma*e*w + gamma1*w*vel;
730  Ai[2](4,1) = -gamma1*u*w;
731  Ai[2](4,2) = -gamma1*v*w;
732  Ai[2](4,3) = gamma*e - gamma1*w2 - gamma2*vel;
733  Ai[2](4,4) = gamma*w;
734 
735  } else if(dim == 2){
736 
737  u = U[1]/U[0];
738  v = U[2]/U[0];
739  e = U[3]/U[0];
740 
741  T u2 = u*u;
742  T v2 = v*v;
743  T vel = u2+v2;
744 
745  Ai[0](0,0) = 0.;
746  Ai[0](0,1) = 1.;
747  Ai[0](0,2) = 0.;
748  Ai[0](0,3) = 0.;
749 
750  Ai[0](1,0) = gamma2*vel-u2;
751  Ai[0](1,1) = -gamma3*u;
752  Ai[0](1,2) = -gamma1*v;
753  Ai[0](1,3) = gamma1;
754 
755  Ai[0](2,0) = -u*v;
756  Ai[0](2,1) = v;
757  Ai[0](2,2) = u;
758  Ai[0](2,3) = 0.;
759 
760  Ai[0](3,0) = -gamma*e*u + gamma1*u*vel;
761  Ai[0](3,1) = gamma*e - gamma1*u2 - gamma2*vel;
762  Ai[0](3,2) = -gamma1*u*v;
763  Ai[0](3,3) = gamma*u;
764 
765  Ai[1](0,0) = 0.;
766  Ai[1](0,1) = 0.;
767  Ai[1](0,2) = 1.;
768  Ai[1](0,3) = 0.;
769 
770  Ai[1](1,0) = -u*v;
771  Ai[1](1,1) = v;
772  Ai[1](1,2) = u;
773  Ai[1](1,3) = 0.;
774 
775  Ai[1](2,0) = gamma2*vel-v2;
776  Ai[1](2,1) = -gamma1*u;
777  Ai[1](2,2) = -gamma3*v;
778  Ai[1](2,3) = gamma1;
779 
780  Ai[1](3,0) = -gamma*e*v + gamma1*v*vel;
781  Ai[1](3,1) = -gamma1*u*v;
782  Ai[1](3,2) = gamma*e - gamma1*v2 - gamma2*vel;
783  Ai[1](3,3) = gamma*v;
784 
785  } else if(dim == 1){
786 
787  u = U[1]/U[0];
788  e = U[2]/U[0];
789 
790  T u2 = u*u;
791  T vel = u2;
792 
793  Ai[0](0,0) = 0.;
794  Ai[0](0,1) = 1.;
795  Ai[0](0,2) = 0.;
796 
797  Ai[0](1,0) = gamma2*vel-u2;
798  Ai[0](1,1) = -gamma3*u;
799  Ai[0](1,2) = gamma1;
800 
801  Ai[0](2,0) = -gamma*e*u + gamma1*u*vel;
802  Ai[0](2,1) = gamma*e - gamma1*u2 - gamma2*vel;
803  Ai[0](2,2) = gamma*u;
804  }
805 }
806 
807 
808 #ifdef _AUTODIFF
809 template <class T>
810 inline REAL val(T & number)
811 {
812  return number.val();
813 }
814 #endif
815 
816 template< class T >
817 inline void TPZEulerConsLaw::Pressure(REAL gamma, int dim, T & press, TPZVec<T> &U)
818 {
819  if(fabs(val(U[0])) < 1.e-6) {
820  PZError << "\nTPZEulerConsLaw::Pressure> Density negative "
821  << U[0] << std::endl;
822  TPZOutofRange obj;
823  throw(obj);
824  }
825  // Press� = (gam-1)*(E - ro*||(u,v,w)||/2)
826  // onde aqui ro_e = E (nota�)
827 
828  int nstate = NStateVariables(dim);
829  press = 0.0;
830  if(U.NElements() != nstate) U.Resize(nstate);
831  if(nstate == 5){
832  //U = (U0,U1,U2,U3,U4) = (ro , ro u , ro v , ro w , ro e)
833  T rho_velocity = ( U[1]*U[1] + U[2]*U[2] + U[3]*U[3] )/U[0];
834  press = ((gamma-1.)*( U[4] - REAL(0.5) * rho_velocity ));
835  } else
836  if(nstate == 4){
837  //U = (U0,U1,U2,U3,U4) = (ro , ro u , ro v , ro e)
838  T rho_velocity = ( U[1]*U[1] + U[2]*U[2] )/U[0];
839  press = ((gamma-1.)*( U[3] - REAL(0.5) * rho_velocity ));
840  } else
841  if(nstate == 3){
842  //U = (U0,U1,U2,U3,U4) = (ro , ro u , ro e)
843  T rho_velocity = ( U[1]*U[1] )/U[0];
844  press = ((gamma-1.)*( U[2] - REAL(0.5) * rho_velocity ));
845  } else {
846  std::cout << "\nTPZEulerConsLaw::Pressure> Unknown case - returning zero\n";
847  press = 0.0;
848  return;
849  }
850  if(val(press) < 0){
851  T temp = ((T)(gamma-1.))*U[nstate-1];
852  PZError << "TPZEulerConsLaw::Pressure> Negative pressure: " << press << " (gama-1)*E = " << temp << std::endl;
853  TPZOutofRange obj;
854  throw(obj);
855  }
856 }
857 
858 //----------------Test Flux
859 template <class T>
860 void TPZEulerConsLaw::Test_Flux(TPZVec<T> &solL, TPZVec<T> &solR, TPZVec<REAL> & normal, REAL gamma, TPZVec<T> & flux)
861 {
862  int nState = flux.NElements();
863  for(int i = 0; i < nState; i++)
864  {
865  flux[i] = (solR[i]-solL[i]);
866  }
867 }
868 
869 
870 //----------------Roe Flux
871 template <class T>
872 void TPZEulerConsLaw::Roe_Flux(TPZVec<T> &solL, TPZVec<T> &solR, TPZVec<REAL> & normal, REAL gamma, TPZVec<T> & flux, int entropyFix)
873 {
874  // Normals outgoing from the BC elements into the
875  // mesh elements -> all the normals are opposited to
876  // the common convention -> changing the left/right
877  // elements and normals.
878  int nState = solL.NElements();
879  if(nState == 5)
880  {
881  Roe_Flux(solL[0], solL[1], solL[2], solL[3], solL[4],
882  solR[0], solR[1], solR[2], solR[3], solR[4],
883  normal[0], normal[1], normal[2],
884  gamma,
885  flux[0], flux[1], flux[2], flux[3], flux[4], entropyFix);
886 
887  }else if(nState == 4)
888  {
889  Roe_Flux(solL[0], solL[1], solL[2], solL[3],
890  solR[0], solR[1], solR[2], solR[3],
891  normal[0], normal[1],
892  gamma,
893  flux[0], flux[1], flux[2], flux[3], entropyFix);
894  }else if(nState == 3)
895  {
896  //using the 2D expression for 1d problem
897  T auxL = REAL(0.),
898  auxR = REAL(0.),
899  fluxaux = REAL(0.);
900  auxL = flux[0];
901  Roe_Flux(solL[0], solL[1], auxL, solL[2],
902  solR[0], solR[1], auxR, solR[2],
903  normal[0], 0,
904  gamma,
905  flux[0], flux[1], fluxaux, flux[2], entropyFix);
906  }else
907  {
908  PZError << "No flux on " << nState << " state variables.\n";
909  }
910 
911 }
912 
913 //left = **_f right = **_t
914 template <class T>
915 inline void TPZEulerConsLaw::Roe_Flux(
916  const T & rho_f, const T & rhou_f, const T & rhov_f, const T & rhow_f,
917  const T & rhoE_f, const T & rho_t, const T & rhou_t, const T & rhov_t, const T & rhow_t,
918  const T & rhoE_t, const REAL nx, const REAL ny, const REAL nz, const REAL gam,
919  T & flux_rho, T &flux_rhou, T &flux_rhov,
920  T & flux_rhow, T &flux_rhoE, int entropyFix){
921 
922  T alpha1,alpha2,alpha3,alpha4,alpha5,alpha;
923  T a1,a2,a3,a4,a5,b1,b2,b3,b4,b5;
924  T ep_t, ep_f, p_t, p_f;
925  T rhouv_t, rhouv_f, rhouw_t, rhouw_f, rhovw_t, rhovw_f;
926  T lambda_f, lambda_t;
927  T delta_rho, delta_rhou, delta_rhov, delta_rhow, delta_rhoE;
928  T hnx, hny, hnz;
929  T tempo11, usc;
930 
931  flux_rho = 0;
932  flux_rhou = 0;
933  flux_rhov = 0;
934  flux_rhow = 0;
935  flux_rhoE = 0;
936 
937  REAL gam1 = gam - 1.0;
938  T irho_f = REAL(1.0)/rho_f;
939  T irho_t = REAL(1.0)/rho_t;
940 
941  //
942  //.. Compute the ROE Averages
943  //
944  //.... some useful quantities
945  T coef1 = sqrt(rho_f);
946  T coef2 = sqrt(rho_t);
947  T somme_coef = coef1 + coef2;
948  T isomme_coef = REAL(1.0)/somme_coef;
949  T u_f = rhou_f*irho_f;
950  T v_f = rhov_f*irho_f;
951  T w_f = rhow_f*irho_f;
952  T h_f = (gam * rhoE_f*irho_f) - (.5*gam1) * (u_f * u_f + v_f * v_f + w_f * w_f);
953  T u_t = rhou_t*irho_t;
954  T v_t = rhov_t*irho_t;
955  T w_t = rhow_t*irho_t;
956  T h_t = (gam * rhoE_t*irho_t) - (.5*gam1) * (u_t * u_t + v_t * v_t + w_t * w_t);
957 
958  //.... averages
959  //REAL rho_ave = coef1 * coef2;
960  T u_ave = (coef1 * u_f + coef2 * u_t) * isomme_coef;
961  T v_ave = (coef1 * v_f + coef2 * v_t) * isomme_coef;
962  T w_ave = (coef1 * w_f + coef2 * w_t) * isomme_coef;
963  T h_ave = (coef1 * h_f + coef2 * h_t) * isomme_coef;
964  //
965  //.. Compute Speed of sound
966  T scal = u_ave * nx + v_ave * ny + w_ave * nz;
967  T norme = sqrt(nx * nx + ny * ny + nz * nz);
968  T inorme = REAL(1.0)/norme;
969  T u2pv2pw2 = u_ave * u_ave + v_ave * v_ave + w_ave * w_ave;
970  T c_speed = gam1 * (h_ave - REAL(0.5) * u2pv2pw2);
971  if(c_speed < REAL(1e-6)) c_speed = 1e-6;// <!> zeroes the derivatives? // avoid division by 0 if critical
972  c_speed = sqrt(c_speed);
973  T c_speed2 = c_speed * norme;
974  //
975  //.. Compute the eigenvalues of the Jacobian matrix
976  T eig_val1 = scal - c_speed2;
977  T eig_val2 = scal;
978  T eig_val3 = scal + c_speed2;
979  //
980  //.. Compute the ROE flux
981  //.... In this part many tests upon the eigenvalues
982  //.... are done to simplify calculations
983  //.... Here we use the two formes of the ROE flux :
984  //.... phi(Wl,Wr) = F(Wl) + A-(Wroe)(Wr - Wl)
985  //.... phi(Wl,Wr) = F(Wr) - A+(Wroe)(Wr - Wl)
986  //
987  if(eig_val2 <= REAL(0.0)) {
988  p_t = gam1 * (rhoE_t - REAL(0.5) * (rhou_t * rhou_t +
989  rhov_t * rhov_t + rhow_t * rhow_t) * irho_t);
990  ep_t = rhoE_t + p_t;
991  rhouv_t = rhou_t * v_t;
992  rhouw_t = rhou_t * w_t;
993  rhovw_t = rhov_t * w_t;
994  flux_rho = rhou_t * nx + rhov_t * ny + rhow_t * nz;
995  flux_rhou = (rhou_t * u_t + p_t) * nx + rhouv_t * ny + rhouw_t * nz;
996  flux_rhov = rhouv_t * nx + (rhov_t * v_t + p_t) * ny + rhovw_t * nz;
997  flux_rhow = rhouw_t * nx + rhovw_t * ny + (rhow_t * w_t + p_t) * nz;
998  flux_rhoE = ep_t * (u_t * nx + v_t * ny + w_t * nz);
999  //
1000  //.... A Entropic modification
1001  //
1002  p_f = gam1 * (rhoE_f - REAL(0.5) * (rhou_f * rhou_f + rhov_f * rhov_f
1003  + rhow_f * rhow_f) * irho_f);
1004  lambda_f = u_f * nx + v_f * ny + w_f * nz + norme
1005  * sqrt(gam * p_f * irho_f);
1006  lambda_t = u_t * nx + v_t * ny + w_t * nz + norme
1007  * sqrt(gam * p_t * irho_t);
1008  if (entropyFix && (lambda_f < REAL(0.)) && (lambda_t > REAL(0.))) {
1009  eig_val3 = lambda_t * (eig_val3 - lambda_f) / (lambda_t - lambda_f);
1010  }
1011  //
1012  if (eig_val3 > REAL(0.0)) {
1013  //.. In this case A+ is obtained by multiplying the last
1014  //.. colomne of T-1 with the last row of T with eig_val3 //Cedric
1015  delta_rho = rho_t - rho_f; //right - left
1016  delta_rhou = rhou_t - rhou_f; //**_t - **_f
1017  delta_rhov = rhov_t - rhov_f;
1018  delta_rhow = rhow_t - rhow_f;
1019  delta_rhoE = rhoE_t - rhoE_f;
1020  //
1021  scal = scal * inorme;
1022  hnx = nx * inorme;
1023  hny = ny * inorme;
1024  hnz = nz * inorme;
1025  usc = REAL(1.0)/c_speed;
1026  tempo11 = gam1 * usc;
1027  //.. Last columne of the matrix T-1
1028  a1 = usc;
1029  a2 = u_ave * usc + hnx;
1030  a3 = v_ave * usc + hny;
1031  a4 = w_ave * usc + hnz;
1032  a5 = REAL(0.5) * u2pv2pw2 * usc + REAL(2.5) * c_speed + scal;
1033  //.. Last row of the matrix T * eig_val3
1034  b1 = REAL(0.5) * (REAL(0.5) * tempo11 * u2pv2pw2 - scal);
1035  b2 = REAL(0.5) * (hnx - tempo11 * u_ave);
1036  b3 = REAL(0.5) * (hny - tempo11 * v_ave);
1037  b4 = REAL(0.5) * (hnz - tempo11 * w_ave);
1038  b5 = REAL(0.5) * tempo11;
1039  //
1040  alpha1 = b1 * delta_rho;
1041  alpha2 = b2 * delta_rhou;
1042  alpha3 = b3 * delta_rhov;
1043  alpha4 = b4 * delta_rhow;
1044  alpha5 = b5 * delta_rhoE;
1045  alpha = eig_val3 * (alpha1 + alpha2 + alpha3 + alpha4 + alpha5);
1046  //
1047  flux_rho -= a1 * alpha;
1048  flux_rhou -= a2 * alpha;
1049  flux_rhov -= a3 * alpha;
1050  flux_rhow -= a4 * alpha;
1051  flux_rhoE -= a5 * alpha;
1052  }
1053  }
1054  //
1055  if(eig_val2 > REAL(0.0)) {
1056  p_f = gam1 * (rhoE_f - REAL(0.5) * (rhou_f * rhou_f +
1057  rhov_f * rhov_f + rhow_f * rhow_f) * irho_f);
1058  ep_f = rhoE_f + p_f;
1059  rhouv_f = rhou_f * v_f;
1060  rhouw_f = rhou_f * w_f;
1061  rhovw_f = rhov_f * w_f;
1062  flux_rho = rhou_f * nx + rhov_f * ny + rhow_f * nz;
1063  flux_rhou = (rhou_f * u_f + p_f) * nx + rhouv_f * ny + rhouw_f * nz;
1064  flux_rhov = rhouv_f * nx + (rhov_f * v_f + p_f) * ny + rhovw_f * nz;
1065  flux_rhow = rhouw_f * nx + rhovw_f * ny + (rhow_f * w_f + p_f) * nz;
1066  flux_rhoE = ep_f * (u_f * nx + v_f * ny + w_f * nz);
1067  //
1068  // A Entropic modification
1069  //
1070  p_t = gam1 * (rhoE_t - REAL(0.5) * (rhou_t * rhou_t +
1071  rhov_t * rhov_t + rhow_t * rhow_t) * irho_t);
1072  lambda_f = u_f * nx + v_f * ny + w_f * nz - norme
1073  * sqrt(gam * p_f * irho_f);
1074  lambda_t = u_t * nx + v_t * ny + w_t * nz - norme
1075  * sqrt(gam * p_t * irho_t);
1076  if (entropyFix && (lambda_f < REAL(0.)) && (lambda_t > REAL(0.))) {
1077  eig_val1 = lambda_f * (lambda_t - eig_val1) / (lambda_t - lambda_f);
1078  }
1079  //
1080  if (eig_val1 < REAL(0.0)) {
1081  //.. In this case A+ is obtained by multiplying the first
1082  //.. columne of T-1 with the first row of T with eig_val1
1083  delta_rho = rho_t - rho_f;
1084  delta_rhou = rhou_t - rhou_f;
1085  delta_rhov = rhov_t - rhov_f;
1086  delta_rhow = rhow_t - rhow_f;
1087  delta_rhoE = rhoE_t - rhoE_f;
1088  //
1089  scal = scal * inorme;
1090  hnx = nx * inorme;
1091  hny = ny * inorme;
1092  hnz = nz * inorme;
1093  usc = REAL(1.0)/c_speed;
1094  tempo11 = gam1 * usc;
1095  //.. First colomne of the matrix T-1
1096  a1 = usc;
1097  a2 = u_ave * usc - hnx;
1098  a3 = v_ave * usc - hny;
1099  a4 = w_ave * usc - hnz;
1100  a5 = REAL(0.5) * u2pv2pw2 * usc + REAL(2.5) * c_speed - scal;
1101  //.. First row of the matrix T * eig_val1
1102  b1 = REAL(0.5) * (REAL(0.5) * tempo11 * u2pv2pw2 + scal);
1103  b2 = -REAL(0.5) * (hnx + tempo11 * u_ave);
1104  b3 = -REAL(0.5) * (hny + tempo11 * v_ave);
1105  b4 = -REAL(0.5) * (hnz + tempo11 * w_ave);
1106  b5 = REAL(0.5) * tempo11;
1107  //
1108  alpha1 = b1 * delta_rho;
1109  alpha2 = b2 * delta_rhou;
1110  alpha3 = b3 * delta_rhov;
1111  alpha4 = b4 * delta_rhow;
1112  alpha5 = b5 * delta_rhoE;
1113  alpha = eig_val1 * (alpha1 + alpha2 + alpha3 + alpha4 + alpha5);
1114  //
1115  flux_rho += a1 * alpha;
1116  flux_rhou += a2 * alpha;
1117  flux_rhov += a3 * alpha;
1118  flux_rhow += a4 * alpha;
1119  flux_rhoE += a5 * alpha;
1120  }
1121  }
1122 }
1123 
1124 //left = **_f right = **_t
1125 template <class T>
1126 inline void TPZEulerConsLaw::Roe_Flux(const T & rho_f, const T & rhou_f, const T & rhov_f, const T & rhoE_f,
1127  const T & rho_t, const T & rhou_t, const T & rhov_t, const T & rhoE_t,
1128  const REAL nx, const REAL ny, const REAL gam,
1129  T & flux_rho, T & flux_rhou,T & flux_rhov, T & flux_rhoE, int entropyFix){
1130 
1131  T alpha1,alpha2,alpha3,alpha4,a1,a2,a3,a4,b1,b2,b3,b4,alpha;
1132  T ep_t, ep_f, p_t, p_f;
1133  T rhouv_t, rhouv_f;
1134  T lambda_f, lambda_t;
1135  T delta_rho, delta_rhou,delta_rhov, delta_rhoE;
1136  REAL hnx, hny;
1137  T tempo11, usc;
1138 
1139  flux_rho = 0;
1140  flux_rhou = 0;
1141  flux_rhov = 0;
1142  flux_rhoE = 0;
1143 
1144  REAL gam1 = gam - REAL(1.0);
1145  //REAL gam2 = gam * (gam - 1.0);
1146  //REAL igam = 1.0 / (gam - 1.0);
1147 
1148  //
1149  //.. Compute the ROE Averages
1150  //
1151  //.... some useful quantities
1152  T coef1 = sqrt(rho_f);
1153  T coef2 = sqrt(rho_t);
1154  T somme_coef = coef1 + coef2;
1155  T u_f = rhou_f/rho_f;
1156  T v_f = rhov_f/rho_f;
1157  T h_f = (gam * rhoE_f/rho_f) - (u_f * u_f + v_f * v_f) * (gam1 / REAL(2.0));
1158  T u_t = rhou_t/rho_t;
1159  T v_t = rhov_t/rho_t;
1160  T h_t = (gam * rhoE_t/rho_t) - (u_t * u_t + v_t * v_t) * (gam1 / REAL(2.0));
1161 
1162  //.... averages
1163  //REAL rho_ave = coef1 * coef2;
1164  T u_ave = (coef1 * u_f + coef2 * u_t) / somme_coef;
1165  T v_ave = (coef1 * v_f + coef2 * v_t) / somme_coef;
1166  T h_ave = (coef1 * h_f + coef2 * h_t) / somme_coef;
1167  //
1168  // cout << "Correct coef1 " << coef1 << "coef2 " << coef2 << "h_f " << h_f << "h_t " << h_t << endl;
1169 
1170  //.. Compute Speed of sound
1171  T scal = u_ave * nx + v_ave * ny;
1172  REAL norme = sqrt(nx * nx + ny * ny);
1173  T u2pv2 = u_ave * u_ave + v_ave * v_ave;
1174  T c_speed = gam1 * (h_ave - REAL(0.5) * u2pv2);
1175  // cout << "c_speed " << c_speed << endl << "h_ave " << h_ave << endl << "u2pv2 " << u2pv2 << endl;
1176 
1177  if(c_speed < REAL(1e-6)) c_speed = REAL(1e-6); // avoid division by 0 if critical
1178  c_speed = sqrt(c_speed);
1179  T c_speed2 = c_speed * norme;
1180  //
1181  //.. Compute the eigenvalues of the Jacobian matrix
1182  T eig_val1 = scal - c_speed2;
1183  T eig_val2 = scal;
1184  T eig_val3 = scal + c_speed2;
1185  // cout << "Eigenvalues correct" << eig_val1 << endl << eig_val2 << endl <<
1186  // eig_val3 << endl;
1187 
1188  //
1189  //.. Compute the ROE flux
1190  //.... In this part many tests upon the eigenvalues
1191  //.... are done to simplify calculations
1192  //.... Here we use the two formes of the ROE flux :
1193  //.... phi(Wl,Wr) = F(Wl) + A-(Wroe)(Wr - Wl)
1194  //.... phi(Wl,Wr) = F(Wr) - A+(Wroe)(Wr - Wl)
1195  //
1196  if(eig_val2 <= REAL(0.0)) {
1197  p_t = gam1 * (rhoE_t - REAL(0.5) * (rhou_t * rhou_t + rhov_t * rhov_t) / rho_t);
1198  ep_t = rhoE_t + p_t;
1199  rhouv_t = rhou_t * v_t;
1200  flux_rho = rhou_t * nx + rhov_t * ny;
1201  flux_rhou = (rhou_t * u_t + p_t) * nx + rhouv_t * ny;
1202  flux_rhov = rhouv_t * nx + (rhov_t * v_t + p_t) * ny;
1203  flux_rhoE = ep_t * (u_t * nx + v_t * ny);
1204  //
1205  //.... A Entropic modification
1206  //
1207  p_f = gam1 * (rhoE_f - REAL(0.5) * (rhou_f * rhou_f + rhov_f * rhov_f) / rho_f);
1208  lambda_f = u_f * nx + v_f * ny + norme * sqrt(gam * p_f / rho_f);
1209  lambda_t = u_t * nx + v_t * ny + norme
1210  * sqrt(gam * p_t / rho_t);
1211  if (entropyFix && (lambda_f < REAL(0.)) && (lambda_t > REAL(0.))) {
1212  eig_val3 = lambda_t * (eig_val3 - lambda_f) / (lambda_t - lambda_f);
1213  }
1214  //
1215  if (eig_val3 > REAL(0.0)) {
1216  //.. In this case A+ is obtained by multiplying the last
1217  //.. colomne of T-1 with the last row of T with eig_val3
1218  delta_rho = rho_t - rho_f;
1219  delta_rhou = rhou_t - rhou_f;
1220  delta_rhov = rhov_t - rhov_f;
1221  delta_rhoE = rhoE_t - rhoE_f;
1222  //
1223  scal = scal / norme;
1224  hnx = nx / norme;
1225  hny = ny / norme;
1226  usc = REAL(1.0)/c_speed;
1227  tempo11 = gam1 * usc;
1228  //.. Last columne of the matrix T-1
1229  a1 = usc;
1230  a2 = u_ave * usc + hnx;
1231  a3 = v_ave * usc + hny;
1232  a4 = REAL(0.5) * u2pv2 * usc + REAL(2.5) * c_speed + scal;
1233  //.. Last row of the matrix T * eig_val3
1234  b1 = REAL(0.5) * eig_val3 * (REAL(0.5) * tempo11 * u2pv2 - scal);
1235  b2 = REAL(0.5) * eig_val3 * (hnx - tempo11 * u_ave);
1236  b3 = REAL(0.5) * eig_val3 * (hny - tempo11 * v_ave);
1237  b4 = REAL(0.5) * eig_val3 * tempo11;
1238  //
1239  alpha1 = a1 * b1 * delta_rho;
1240  alpha2 = a1 * b2 * delta_rhou;
1241  alpha3 = a1 * b3 * delta_rhov;
1242  alpha4 = a1 * b4 * delta_rhoE;
1243  alpha = alpha1 + alpha2 + alpha3 + alpha4;
1244  //
1245  flux_rho -= alpha;
1246  flux_rhou -= a2 * b1 * delta_rho + a2 * b2 * delta_rhou +
1247  a2 * b3 * delta_rhov + a2 * b4 * delta_rhoE;
1248  flux_rhov -= a3 * b1 * delta_rho + a3 * b2 * delta_rhou +
1249  a3 * b3 * delta_rhov + a3 * b4 * delta_rhoE;
1250  flux_rhoE -= a4 * b1 * delta_rho + a4 * b2 * delta_rhou +
1251  a4 * b3 * delta_rhov + a4 * b4 * delta_rhoE;
1252  }
1253  }
1254  //
1255  if(eig_val2 > REAL(0.0)) {
1256  p_f = gam1 * (rhoE_f - REAL(0.5) * (rhou_f * rhou_f +
1257  rhov_f * rhov_f) / rho_f);
1258  ep_f = rhoE_f + p_f;
1259  rhouv_f = rhou_f * v_f;
1260  flux_rho = rhou_f * nx + rhov_f * ny;
1261  flux_rhou = (rhou_f * u_f + p_f) * nx + rhouv_f * ny;
1262  flux_rhov = rhouv_f * nx + (rhov_f * v_f + p_f) * ny;
1263  flux_rhoE = ep_f * (u_f * nx + v_f * ny);
1264  //
1265  // A Entropic modification
1266  //
1267  p_t = gam1 * (rhoE_t - REAL(0.5) * (rhou_t * rhou_t +
1268  rhov_t * rhov_t) / rho_t);
1269  lambda_f = u_f * nx + v_f * ny - norme * sqrt(gam * p_f / rho_f);
1270  lambda_t = u_t * nx + v_t * ny - norme * sqrt(gam * p_t / rho_t);
1271  if (entropyFix && (lambda_f < REAL(0.)) && (lambda_t > REAL(0.))) {
1272  eig_val1 = lambda_f * (lambda_t - eig_val1) / (lambda_t - lambda_f);
1273  }
1274  //
1275  if (eig_val1 < REAL(0.0)) {
1276  //.. In this case A+ is obtained by multiplying the first
1277  //.. columne of T-1 with the first row of T with eig_val1
1278  delta_rho = rho_t - rho_f;
1279  delta_rhou = rhou_t - rhou_f;
1280  delta_rhov = rhov_t - rhov_f;
1281  delta_rhoE = rhoE_t - rhoE_f;
1282  //
1283  scal = scal / norme;
1284  hnx = nx / norme;
1285  hny = ny / norme;
1286  usc = REAL(1.0)/c_speed;
1287  tempo11 = gam1 * usc;
1288  //.. First colomne of the matrix T-1
1289  a1 = usc;
1290  a2 = u_ave * usc - hnx;
1291  a3 = v_ave * usc - hny;
1292  a4 = REAL(0.5) * u2pv2 * usc + REAL(2.5) * c_speed - scal;
1293  //.. First row of the matrix T * eig_val1
1294  b1 = REAL(0.5) * eig_val1 * (REAL(0.5) * tempo11 * u2pv2 + scal);
1295  b2 = -REAL(0.5) * eig_val1 * (hnx + tempo11 * u_ave);
1296  b3 = -REAL(0.5) * eig_val1 * (hny + tempo11 * v_ave);
1297  b4 = REAL(0.5) * eig_val1 * tempo11;
1298  //
1299  alpha1 = a1 * b1 * delta_rho;
1300  alpha2 = a1 * b2 * delta_rhou;
1301  alpha3 = a1 * b3 * delta_rhov;
1302  alpha4 = a1 * b4 * delta_rhoE;
1303  alpha = alpha1 + alpha2 + alpha3 + alpha4;
1304  //
1305  flux_rho += alpha;
1306  flux_rhou += a2 * b1 * delta_rho + a2 * b2 * delta_rhou +
1307  a2 * b3 * delta_rhov + a2 * b4 * delta_rhoE;
1308  flux_rhov += a3 * b1 * delta_rho + a3 * b2 * delta_rhou +
1309  a3 * b3 * delta_rhov + a3 * b4 * delta_rhoE;
1310  flux_rhoE += a4 * b1 * delta_rho + a4 * b2 * delta_rhou +
1311  a4 * b3 * delta_rhov + a4 * b4 * delta_rhoE;
1312  }
1313  }
1314 }
1315 
1316 
1317 template <class T>
1318 void TPZEulerConsLaw::cSpeed(TPZVec<T> & sol, REAL gamma, T & c)
1319 {
1320  if(sol[0] < REAL(1e-10))
1321  {
1322  PZError << "TPZEulerConsLaw::cSpeed Too low or negative density\n";
1323  TPZOutofRange obj;
1324  throw(obj);
1325  }
1326 
1327  int dim = sol.NElements() - 2;
1328  T press, temp;
1329  TPZEulerConsLaw::Pressure(gamma, dim, press, sol);
1330  temp = gamma * press;
1331 
1332  if(temp < REAL(1e-10)) // too low or negative
1333  {
1334  PZError << "TPZEulerConsLaw::cSpeed Too low or negative numerator\n";
1335  }
1336  c = sqrt(gamma * press/ sol[0]);
1337 }
1338 
1339 template <class T>
1340 inline void TPZEulerConsLaw::uRes(TPZVec<T> & sol, T & us)
1341 {
1342  if(sol[0] < REAL(1e-10))
1343  {
1344  PZError << "TPZEulerConsLaw::cSpeed Too low or negative density\n";
1345  TPZOutofRange obj;
1346  throw(obj);
1347  }
1348 
1349  T temp;
1350  switch(sol.NElements())
1351  {
1352  case(3):
1353  us = sol[1]/sol[0];
1354  case(4):
1355  temp = sol[1]*sol[1] + sol[2]*sol[2];
1356  if(temp < REAL(1e-40))
1357  {
1358  PZError << "TPZEulerConsLaw::uRes Zero Velocity\n";
1359  TPZOutofRange obj;
1360  throw(obj);
1361  // exit(-1);
1362  }
1363  us = sqrt(temp)/sol[0];
1364  break;
1365  case(5):
1366  temp = sol[1]*sol[1] + sol[2]*sol[2] + sol[3]*sol[3];
1367  if(temp < REAL(1e-40))
1368  {
1369  PZError << "TPZEulerConsLaw::uRes Zero Velocity\n";
1370  TPZOutofRange obj;
1371  throw(obj);
1372  // exit(-1);
1373  }
1374  us = sqrt(temp)/sol[0];
1375  break;
1376  default:
1377  PZError << "TPZEulerConsLaw::uRes Error: invalid Dimension\n";
1378  }
1379 }
1380 
1381 #endif
TPZEulerConsLaw::ContributeAdv
virtual void ContributeAdv(TPZVec< REAL > &x, TPZFMatrix< REAL > &jacinv, TPZVec< STATE > &sol, TPZFMatrix< STATE > &dsol, REAL weight, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphi, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef)
Definition: pzeulerconslaw.cpp:502
TPZEulerConsLaw
This material implements the weak statement of the compressible euler equations.
Definition: pzeulerconslaw.h:33
fabs
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ fabs
Definition: tfadfunc.h:140
TPZEulerConsLaw::JacobFlux
static void JacobFlux(REAL gamma, int dim, TPZVec< T > &U, TPZVec< TPZDiffMatrix< T > > &Ai)
Jacobian of the tensor flux of Euler.
Definition: pzeulerconslaw.h:615
TPZEulerConsLaw::VariableIndex
virtual int VariableIndex(const std::string &name) override
Returns the relative index of a variable according to its name.
Definition: pzeulerconslaw.cpp:157
TPZEulerConsLaw::Solution
virtual void Solution(TPZMaterialData &data, int var, TPZVec< STATE > &Solout) override
returns the solution associated with the var index based on the finite element approximation ...
Definition: pzeulerconslaw.h:136
pzlog.h
Contains definitions to LOGPZ_DEBUG, LOGPZ_INFO, LOGPZ_WARN, LOGPZ_ERROR and LOGPZ_FATAL, and the implementation of the inline InitializePZLOG(string) function using log4cxx library or not. It must to be called out of "#ifdef LOG4CXX" scope.
TPZEulerConsLaw::Test_Flux
void Test_Flux(TPZVec< T > &solL, TPZVec< T > &solR, TPZVec< REAL > &normal, REAL gamma, TPZVec< T > &flux)
Test flux -> returns the averaged state variables across an interface.
Definition: pzeulerconslaw.h:860
TPZEulerConsLaw::ContributeBCInterface
virtual void ContributeBCInterface(TPZMaterialData &data, TPZMaterialData &dataleft, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc) override
It computes a contribution to stiffness matrix and load vector at one BC integration point...
Definition: pzeulerconslaw.cpp:702
bc
clarg::argBool bc("-bc", "binary checkpoints", false)
TPZMaterial.h
TPZOutofRange
This class is used as an exception thrown on an outofrange condition.
Definition: tpzoutofrange.h:24
TPZEulerConsLaw::ClassId
virtual int ClassId() const override
Class identificator.
Definition: pzeulerconslaw.cpp:2515
TPZEulerConsLaw::fConvFace
TPZTimeDiscr fConvFace
Definition: pzeulerconslaw.h:544
pzvec.h
Templated vector implementation.
TPZEulerConsLaw::Name
virtual std::string Name() override
Returns the material name.
Definition: pzeulerconslaw.h:550
TPZEulerConsLaw::ContributeExplConvFace
void ContributeExplConvFace(TPZVec< REAL > &x, TPZVec< STATE > &solL, TPZVec< STATE > &solR, REAL weight, TPZVec< REAL > &normal, TPZFMatrix< REAL > &phiL, TPZFMatrix< REAL > &phiR, TPZFMatrix< STATE > &ef, int entropyFix=1)
Definition: pzeulerconslaw.cpp:1347
TPZEulerConsLaw::ContributeBC
virtual void ContributeBC(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ef, TPZBndCond &bc) override
Contributes to the residual vector the boundary conditions.
Definition: pzeulerconslaw.h:345
TPZEulerConsLaw::ComputeGhostState
void ComputeGhostState(TPZVec< T > &solL, TPZVec< T > &solR, TPZVec< REAL > &normal, TPZBndCond &bc, int &entropyFix)
Computes the ghost state variables bsed on the BC type.
Definition: pzeulerconslaw.cpp:937
TPZEulerConsLaw::TPZEulerConsLaw
TPZEulerConsLaw(const TPZEulerConsLaw &cp)
Definition: pzeulerconslaw.h:45
TPZEulerConsLaw::SetDelta
void SetDelta(REAL delta)
Sets the delta parameter inside the artifficial diffusion term.
Definition: pzeulerconslaw.cpp:272
TPZEulerConsLaw::fConvVol
TPZTimeDiscr fConvVol
Definition: pzeulerconslaw.h:544
degree
void degree(int root, int adj_num, int adj_row[], int adj[], int mask[], int deg[], int *iccsze, int ls[], int node_num)
Definition: rcm.cpp:875
TPZEulerConsLaw::ContributeBC
virtual void ContributeBC(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc) override
Contributes to the residual vector the boundary conditions.
Definition: pzeulerconslaw.cpp:662
TPZEulerConsLaw::Flux
void Flux(TPZVec< T > &U, TPZVec< T > &Fx, TPZVec< T > &Fy, TPZVec< T > &Fz)
tensor of the three-dimensional flux of Euler
Definition: pzeulerconslaw.h:556
TPZVec
This class implements a simple vector storage scheme for a templated class T. Utility.
Definition: pzgeopoint.h:19
val
REAL val(STATE &number)
Returns value of the variable.
Definition: pzartdiff.h:23
TPZTimeDiscr
TPZTimeDiscr
Indicates the type of time discretization.
Definition: pzconslaw.h:34
TPZArtDiff
This class adds to the term of diffusion to the variacional formulation of the differential equation...
Definition: pzartdiff.h:62
TPZEulerConsLaw::uRes
static void uRes(TPZVec< T > &sol, T &us)
Definition: pzeulerconslaw.h:1340
TPZEulerConsLaw::ContributeApproxImplDiff
void ContributeApproxImplDiff(TPZVec< REAL > &x, TPZFMatrix< REAL > &jacinv, TPZVec< STATE > &sol, TPZFMatrix< STATE > &dsol, REAL weight, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphi, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef)
Definition: pzeulerconslaw.cpp:1263
TPZEulerConsLaw::ContributeApproxImplConvFace
void ContributeApproxImplConvFace(TPZVec< REAL > &x, REAL faceSize, TPZVec< STATE > &solL, TPZVec< STATE > &solR, REAL weight, TPZVec< REAL > &normal, TPZFMatrix< REAL > &phiL, TPZFMatrix< REAL > &phiR, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, int entropyFix=1)
Definition: pzeulerconslaw.cpp:1451
TPZEulerConsLaw::Write
virtual void Write(TPZStream &buf, int withclassid) const override
Saves the element data to a stream.
Definition: pzeulerconslaw.cpp:1946
TPZEulerConsLaw::NFluxes
virtual int NFluxes() override
Returns the number of fluxes associated to this material.
Definition: pzeulerconslaw.cpp:211
TPZEulerConsLaw::ApproxRoe_Flux
static void ApproxRoe_Flux(TPZVec< T > &solL, TPZVec< T > &solR, TPZVec< REAL > &normal, REAL gamma, TPZVec< T > &flux, int entropyFix=1)
This flux encapsulates the two and three dimensional fluxes acquired from the Mouse program...
Definition: pzeulerconslaw.cpp:1984
TPZEulerConsLaw::Print
virtual void Print(std::ostream &out) override
Prints the state of internal variables.
Definition: pzeulerconslaw.cpp:102
TPZEulerConsLaw::Read
void Read(TPZStream &buf, void *context) override
Reads the element data from a stream.
Definition: pzeulerconslaw.cpp:1958
TPZEulerConsLaw::Det
REAL Det(TPZFMatrix< REAL > &Mat)
Computes the determinant of a 2d or 3d matrix. Used by recompute the element size.
Definition: pzeulerconslaw.cpp:184
TPZEulerConsLaw::ContributeImplT1
void ContributeImplT1(TPZVec< REAL > &x, TPZVec< STATE > &sol, TPZFMatrix< STATE > &dsol, REAL weight, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphi, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef)
Definition: pzeulerconslaw.cpp:1873
TPZMaterial
This abstract class defines the behaviour which each derived class needs to implement.
Definition: TPZMaterial.h:39
TPZConservationLaw::fDim
int fDim
Dimension of the problem.
Definition: pzconslaw.h:243
TPZEulerConsLaw::ContributeExplT2
void ContributeExplT2(TPZVec< REAL > &x, TPZVec< STATE > &sol, REAL weight, TPZFMatrix< REAL > &phi, TPZFMatrix< STATE > &ef)
Definition: pzeulerconslaw.cpp:1912
TPZVec::Resize
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object reallocating the necessary storage, copying the existing objects to the new...
Definition: pzvec.h:373
TPZEulerConsLaw::ContributeFastestBCInterface
void ContributeFastestBCInterface(int dim, TPZVec< REAL > &x, TPZVec< STATE > &solL, TPZFMatrix< STATE > &dsolL, REAL weight, TPZVec< REAL > &normal, TPZFMatrix< REAL > &phiL, TPZFMatrix< REAL > &dphiL, TPZFMatrix< REAL > &axesleft, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc)
TPZEulerConsLaw::DeltaX
REAL DeltaX(REAL detJac)
Estimates the deltax (element diameter) based on the inverse of the jacobian.
Definition: pzeulerconslaw.cpp:179
TPZEulerConsLaw::fArtDiff
TPZArtDiff fArtDiff
diffusive term
Definition: pzeulerconslaw.h:541
TPZEulerConsLaw::SetTimeDiscr
void SetTimeDiscr(TPZTimeDiscr Diff, TPZTimeDiscr ConvVol, TPZTimeDiscr ConvFace)
Configures the time discretization of some contributions.
Definition: pzeulerconslaw.cpp:61
TPZDiscontinuousGalerkin::ContributeBC
virtual void ContributeBC(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc) override=0
It computes a contribution to the stiffness matrix and load vector at one BC integration point...
pzfmatrix.h
Contains TPZMatrixclass which implements full matrix (using column major representation).
TPZBndCond
This class defines the boundary condition for TPZMaterial objects.
Definition: pzbndcond.h:29
TPZEulerConsLaw::ContributeLast
virtual void ContributeLast(TPZVec< REAL > &x, TPZFMatrix< REAL > &jacinv, TPZVec< STATE > &sol, TPZFMatrix< STATE > &dsol, REAL weight, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphi, TPZFMatrix< STATE > &ef)
Definition: pzeulerconslaw.cpp:474
sqrt
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ sqrt
Definition: tfadfunc.h:120
TPZEulerConsLaw::ContributeFastestBCInterface_dim
void ContributeFastestBCInterface_dim(TPZVec< REAL > &x, TPZVec< STATE > &solL, TPZFMatrix< STATE > &dsolL, REAL weight, TPZVec< REAL > &normal, TPZFMatrix< REAL > &phiL, TPZFMatrix< REAL > &dphiL, TPZFMatrix< REAL > &axesleft, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc)
TPZEulerConsLaw::ContributeInterface
virtual void ContributeInterface(TPZMaterialData &data, TPZMaterialData &dataleft, TPZMaterialData &dataright, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef) override
Contributes to the residual vector and tangent matrix the face-based quantities.
Definition: pzeulerconslaw.cpp:577
TPZEulerConsLaw::ContributeExplT1
void ContributeExplT1(TPZVec< REAL > &x, TPZVec< STATE > &sol, TPZFMatrix< STATE > &dsol, REAL weight, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphi, TPZFMatrix< STATE > &ef)
Definition: pzeulerconslaw.cpp:1840
TPZEulerConsLaw::NStateVariables
virtual int NStateVariables() const override
Object-based overload.
Definition: pzeulerconslaw.cpp:91
TPZEulerConsLaw::OptimalCFL
REAL OptimalCFL(int degree)
returns the best value for the CFL number based on the interpolation degree.
Definition: pzeulerconslaw.cpp:68
gamma
long double gamma(unsigned int n)
Evaluate the factorial of a integer.
Definition: tpzgaussrule.cpp:567
TPZEulerConsLaw::ContributeExplDiff
void ContributeExplDiff(TPZVec< REAL > &x, TPZFMatrix< REAL > &jacinv, TPZVec< STATE > &sol, TPZFMatrix< STATE > &dsol, REAL weight, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphi, TPZFMatrix< STATE > &ef)
Definition: pzeulerconslaw.cpp:1284
TPZConservationLaw
Implements the interface for conservation laws, keeping track of the timestep as well.
Definition: pzconslaw.h:71
pzartdiff.h
Contains the TPZArtDiff class which implements a numerical diffusivity coefficient.
TPZEulerConsLaw::cSpeed
static void cSpeed(TPZVec< T > &sol, REAL gamma, T &c)
Evaluates the speed of sound in the fluid.
Definition: pzeulerconslaw.h:1318
TPZConservationLaw::Solution
virtual void Solution(TPZVec< STATE > &Sol, TPZFMatrix< STATE > &DSol, TPZFMatrix< REAL > &axes, int var, TPZVec< STATE > &Solout) override=0
TPZDiffMatrix
Matrix class to hold the flux derivatives A B C and diffusive matrix coefficients. Matrix.
Definition: pzdiffmatrix.h:27
TPZEulerConsLaw::Solution
virtual void Solution(TPZVec< STATE > &Sol, TPZFMatrix< STATE > &DSol, TPZFMatrix< REAL > &axes, int var, TPZVec< STATE > &Solout) override
Definition: pzeulerconslaw.cpp:218
TPZEulerConsLaw::ContributeExplConvVol
void ContributeExplConvVol(TPZVec< REAL > &x, TPZVec< STATE > &sol, REAL weight, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphi, TPZFMatrix< STATE > &ef)
Definition: pzeulerconslaw.cpp:1763
TPZEulerConsLaw::fDiff
TPZTimeDiscr fDiff
variables indication whether the following terms are implicit
Definition: pzeulerconslaw.h:544
TPZEulerConsLaw::Roe_Flux
static void Roe_Flux(TPZVec< T > &solL, TPZVec< T > &solR, TPZVec< REAL > &normal, REAL gamma, TPZVec< T > &flux, int entropyFix=1)
This flux encapsulates the two and three dimensional fluxes acquired from the Mouse program...
Definition: pzeulerconslaw.h:872
TPZConservationLaw::fGamma
REAL fGamma
Ratio between specific heat is constant and the specific heat the constant volume of a polytropic gas...
Definition: pzconslaw.h:255
TPZEulerConsLaw::NewMaterial
virtual TPZMaterial * NewMaterial() override
To create another material of the same type.
Definition: pzeulerconslaw.h:51
TPZStream
Defines the interface for saving and reading data. Persistency.
Definition: TPZStream.h:50
tpzoutofrange.h
Contains the TPZOutofRange class.
TPZVec::NElements
int64_t NElements() const
Returns the number of elements of the vector.
Definition: pzvec.h:190
pzconslaw.h
Contains the TPZConservationLaw class which implements the interface for conservation laws...
TPZEulerConsLaw::ContributeImplConvVol
void ContributeImplConvVol(TPZVec< REAL > &x, TPZVec< STATE > &sol, TPZFMatrix< STATE > &dsol, REAL weight, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphi, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef)
Definition: pzeulerconslaw.cpp:1794
TPZEulerConsLaw::NSolutionVariables
virtual int NSolutionVariables(int var) override
Returns the number of variables associated with the variable indexed by var.
Definition: pzeulerconslaw.cpp:169
TPZEulerConsLaw::ArtDiff
TPZArtDiff & ArtDiff()
Returns a reference to the artificial diffusion term.
Definition: pzeulerconslaw.h:65
TPZEulerConsLaw::Pressure
static void Pressure(REAL gamma, int dim, T &press, TPZVec< T > &U)
Thermodynamic pressure determined by the law of an ideal gas.
Definition: pzeulerconslaw.h:817
TPZArtDiffType
TPZArtDiffType
Enumerate to define the possible types of artificial diffusion term to stabilize the numerical scheme...
Definition: pzartdiff.h:43
PZError
#define PZError
Defines the output device to error messages and the DebugStop() function.
Definition: pzerror.h:15
TPZEulerConsLaw::SetTimeStep
virtual REAL SetTimeStep(REAL maxveloc, REAL deltax, int degree) override
See declaration in base class.
Definition: pzeulerconslaw.cpp:73
test.obj
obj
Definition: test.py:225
TPZEulerConsLaw::Contribute
virtual void Contribute(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef) override
Contributes to the residual vector and tangent matrix the volume-based quantities.
Definition: pzeulerconslaw.cpp:392
Generated on Sun Aug 16 2020 11:17:03 for NeoPZ by   doxygen 1.8.13