gamma = 1.4
class Compressible2DEuler:
    def pressure(U):
        return (U[3]-(U[1]**2+U[2]**2)/2)*(gamma-1)
    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):
        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]
        res = as_matrix( [ [rho*u, rho*v],
                           [rho*u**2 + p, rho*u*v],
                           [rho*u*v, rho*v**2 + p],
                           [(rE+p)*u, (rE+p)*v]
                         ])
        return res
    def alpha(U,n):
        rho, u,v, p = Compressible2DEuler.toPrim(U)
        return abs(u*n[0]+v*n[1]) + sqrt(gamma*p/rho)
class LLFFlux:
    def flux(Descr,U,n,dt):
        return jump(U)
        UL = U('-')
        UR = U('+')
        return UL-UR
        FL = Descr.F_c(UL)*n
        FR = Descr.F_c(UR)*n
        flux = (FL+FR)
        max_value = lambda a,b: (a+b+abs(a-b))/2.
        visc = max_value( Descr.alpha(UL,n), Descr.alpha(UR,n) )
        return flux + visc/dt*(UR-UL)

Model = Compressible2DEuler
NumFLux = LLFFlux

space = Space(2,4)
u = TrialFunction(space)
v = TestFunction(space)
dt = NamedConstant(triangle, "dt")    # time step

# just a starting point...
def femdgModel(description):
    return inner(description.F_c(u),grad(v))*dx

F = -dt*femdgModel(Model)