Latest changes.

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

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
              >

pydemo/euler/euler.ufl

 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