Make sod work

dune/fem-dg/models/modelwrapper.hh

@@ -338,15 +338,6 @@ namespace Fem
                           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 )
-        // we store \rho u but do not need to divide by \rho here since only
-        // sign is needed.
-        velocity[i] = u[i+1];
-      }
-      */
     }
 
     // we have physical check for this model
@@ -382,6 +373,7 @@ namespace Fem
                      const RangeType& uRight,
                      RangeType& jump) const
     {
+      jump = AdditionalType :: jump( it, x, uLeft, uRight );
       // take pressure as shock detection values
       //const RangeFieldType pl = pressure( uLeft );
       //const RangeFieldType pr = pressure( uRight );

pydemo/euler/euler.ufl

@@ -2,7 +2,7 @@ gamma = 1.4
 class Compressible2DEuler:
     def rhoeps(U):
         v = 0.5*(U[1]**2 + U[2]**2)/U[0]
-        rE = U[-1]
+        rE = U[3]
         return rE - v
     def pressure(U):
         return (gamma-1)*Compressible2DEuler.rhoeps(U)
@@ -12,7 +12,7 @@ class Compressible2DEuler:
         return [U[0],U[1]/U[0],U[2]/U[0],Compressible2DEuler.pressure(U)]
     def F_c(U):
         rho, u,v, p = Compressible2DEuler.toPrim(U)
-        rE = U[-1]
+        rE = U[3]
         res = as_matrix( [ [rho*u, rho*v],
                            [rho*u**2 + p, rho*u*v],
                            [rho*u*v, rho*v**2 + p],
@@ -26,6 +26,10 @@ class Compressible2DEuler:
         return as_vector( [U[1],U[2]] )
     def physical(U):
         return conditional( (U[0]>1e-8), conditional( Compressible2DEuler.rhoeps(U) > 1e-8 , 1, 0 ), 0 )
+    def jump(U,V):
+        pL = Compressible2DEuler.pressure(U)
+        pR = Compressible2DEuler.pressure(V)
+        return (pL - pR)/(0.5*(pL + pR))
 
 Model = Compressible2DEuler
 

python/dune/femdg/_operators.py

@@ -234,6 +234,16 @@ def createFemDGSolver(Model, space):
             targs=['class Entity, class Point'], static=True,
             predefined=predefined)
 
+    jump = getattr(Model,"jump",None)
+    jump = jump(u,v)
+    generateMethod(struct, jump,
+            'double', 'jump',
+            args=['const Intersection& it', 'const Point &x',
+                  'const RangeType &u',
+                  'const RangeType &v'],
+            targs=['class Intersection, class Point'], static=True,
+            predefined=predefined)
+
     writer = SourceWriter(StringWriter())
     writer.emit([struct])