diff --git a/dune/fem-dg/models/modelwrapper.hh b/dune/fem-dg/models/modelwrapper.hh index 38cc7a7d4c0edaf14d824e25c76f3748c5333425..65f3d2343799b2b18dd2e3202f6a59568947c4ca 100644 --- a/dune/fem-dg/models/modelwrapper.hh +++ b/dune/fem-dg/models/modelwrapper.hh @@ -245,6 +245,7 @@ namespace Fem template <class LocalEvaluation> inline bool hasBoundaryValue( const LocalEvaluation& local ) const { + return true; Dune::FieldVector< int, dimRange > bndIds; return impl_.isDirichletIntersection( local.intersection(), bndIds ); } @@ -255,6 +256,10 @@ namespace Fem const RangeType& uLeft, RangeType& uRight ) const { + uRight = uLeft; + return; + + Dune::FieldVector< int, dimRange > bndIds; #ifndef NDEBUG const bool isDirichlet = @@ -332,6 +337,8 @@ namespace Fem const RangeType& u, DomainType& velocity) const { + velocity = AdditionalType :: velocity( en, x, u ); + /* for(int i=0; i<dimDomain; ++i) { // U = (rho, rho v_0,...,rho v_(d-1), e ) @@ -339,6 +346,7 @@ namespace Fem // sign is needed. velocity[i] = u[i+1]; } + */ } // we have physical check for this model @@ -350,17 +358,10 @@ namespace Fem // calculate jump between left and right value template< class Entity > inline bool physical(const Entity& entity, - const DomainType& xGlobal, + const DomainType& x, const RangeType& u) const { - if (u[0]<1e-8) - return false; - else - { - //std::cout << eulerFlux_.rhoeps(u) << std::endl; - // return (eulerFlux_.rhoeps(u) > 1e-8); - return true ; - } + return AdditionalType :: physical( entity, x, u ) > 0; } // adjust average value if necessary diff --git a/dune/fem-dg/solver/dg.hh b/dune/fem-dg/solver/dg.hh index 5f151ab65239492ae962b83719cbec926c2f753c..e9d4c9f59b3f83a72219234cbf8a466dfa4157ec 100644 --- a/dune/fem-dg/solver/dg.hh +++ b/dune/fem-dg/solver/dg.hh @@ -30,7 +30,7 @@ namespace Fem class Additional, Solver::Enum solverId = Solver::Enum::fem, Formulation::Enum formId = Formulation::Enum::primal, - AdvectionLimiter::Enum limiterId = AdvectionLimiter::Enum::limited, + AdvectionLimiter::Enum limiterId = AdvectionLimiter::Enum::unlimited, AdvectionFlux::Enum advFluxId = AdvectionFlux::Enum::llf, DiffusionFlux::Enum diffFluxId = DiffusionFlux::Enum::primal > diff --git a/pydemo/euler/euler.ufl b/pydemo/euler/euler.ufl index 588738be267b80e2e1f9acdd897a37f11d7c4d01..4efd3b885f440daa351f10df981fefe4bdc43e13 100644 --- a/pydemo/euler/euler.ufl +++ b/pydemo/euler/euler.ufl @@ -1,7 +1,11 @@ gamma = 1.4 class Compressible2DEuler: + def rhoeps(U): + v = 0.5*(U[1]**2 + U[2]**2)/U[0] + rE = U[-1] + return rE - v def pressure(U): - return (U[3]-(U[1]**2+U[2]**2)/2)*(gamma-1) + return (gamma-1)*Compressible2DEuler.rhoeps(U) def toCons(U): return [U[0],U[0]*U[1],U[0]*U[2],U[3]/(gamma-1)+U[0]/2*(U[1]**2+U[2]**2)] def toPrim(U): @@ -20,8 +24,8 @@ class Compressible2DEuler: return abs(u*n[0]+v*n[1]) + sqrt(gamma*p/rho) def velocity(U): return as_vector( [U[1],U[2]] ) -# def physical(U) -# return conditional( (U[0]>1e-8) ) + def physical(U): + return conditional( (U[0]>1e-8), conditional( Compressible2DEuler.rhoeps(U) > 1e-8 , 1, 0 ), 0 ) Model = Compressible2DEuler