some clean up

pydemo/euler/euler.py

@@ -16,17 +16,19 @@ def CompressibleEuler(dim, gamma):
             v = as_vector( [U[i] for i in range(1,dim+1)] )
             rhoEps = U[dim+1]/(gamma-1.)
             kin = dot(v,v) * 0.5*U[0]
-            return [U[0], *(U[0]*v), rhoEps+kin]
+            return as_vector( [U[0], *(U[0]*v), rhoEps+kin] )
         def toPrim(U):
-            return [U[0], Model.velocity(U), Model.pressure(U)]
+            return U[0], Model.velocity(U), Model.pressure(U)
         def F_c(U):
             assert dim==2
             rho, v, p = Model.toPrim(U)
             rE = U[dim+1]
-            res = [ [rho*v[0], rho*v[1]],
+            # TODO make indpendent of dim
+            res = as_matrix([
+                    [rho*v[0], rho*v[1]],
                     [rho*v[0]*v[0] + p, rho*v[0]*v[1]],
                     [rho*v[0]*v[1], rho*v[1]*v[1] + p],
-                    [(rE+p)*v[0], (rE+p)*v[1]] ]
+                    [(rE+p)*v[0], (rE+p)*v[1]] ] )
             return res
         def maxLambda(U,n):
             rho, v, p = Model.toPrim(U)
@@ -48,7 +50,7 @@ def CompressibleEuler(dim, gamma):
     return Model
 
 def riemanProblem(x,x0,UL,UR):
-    return conditional(x[0]<x0,as_vector(UL),as_vector(UR))
+    return conditional(x[0]<x0,UL,UR)
 
 def sod(dim,gamma):
     space = Space(dim,dim+2)

pydemo/scalar/scalar.py

@@ -4,7 +4,7 @@ class Transport1D:
     name = 'transport'
 
     def F_c(U):
-        return [ [U[0], 0] ]
+        return as_matrix( [ [U[0], 0] ] )
 
     def outflowValue(x,u):
         return u
@@ -13,38 +13,42 @@ class Transport1D:
     def maxLambda(U,n):
         return abs(n[0])
     def velocity(U):
-        return [1,0]
+        return as_vector( [1,0] )
     def physical(U):
         return 1
     def jump(U,V):
-        return U[0]-V[0]
+        return U-V
 
 class LinearAdvectionDiffusion1D:
     name = 'linearAD'
 
-    def F_c(U):
-        return [ [U[0], 0] ]
+    # def F_c(U):
+    #     return as_matrix( [ [U[0], 0] ] )
     def F_v(U,DU):
         return 0.1*DU
-
-    def dirichletValue(x,u):
-        return as_vector( [cos(2*pi*x[0])] )
-    boundaryValue = {1: dirichletValue}
-
     def maxLambda(U,n):
         return abs(n[0])
     def velocity(U):
-        return [1,0]
+        return as_vector( [0,0] )
     def physical(U):
         return 1
     def jump(U,V):
-        return U[0]-V[0]
+        return U-V
+class LinearAdvectionDiffusion1DDirichlet(LinearAdvectionDiffusion1D):
+    def dirichletValue(x,u):
+        return as_vector( [cos(2*pi*x[0])] )
+    boundaryValue = {1: dirichletValue}
+class LinearAdvectionDiffusion1DNeuman(LinearAdvectionDiffusion1D):
+    def outflowFlux(x,u,n):
+        return 0
+        # return LinearAdvectionDiffusion1D.F_c(u)*n
+    boundaryFlux = {1: outflowFlux}
 
 class Burgers1D:
     name = 'burgers'
 
     def F_c(U):
-        return [ [U[0]*U[0]/2, 0] ]
+        return as_matrix( [ [U[0]*U[0]/2, 0] ] )
 
     def outflowValue(x,u):
         return u
@@ -53,11 +57,11 @@ class Burgers1D:
     def maxLambda(U,n):
         return abs(U[0]*n[0])
     def velocity(U):
-        return [U[0],0]
+        return as_vector( [U[0],0] )
     def physical(U):
         return 1
     def jump(U,V):
-        return U[0]-V[0]
+        return U-V
 
 space = Space(2,1)
 x = SpatialCoordinate(space.cell())
@@ -67,8 +71,11 @@ def riemanProblem(x,x0,UL,UR):
 def constantTransport():
     return Transport1D, as_vector( [0.1] )
 
+def sinProblem():
+    return LinearAdvectionDiffusion1DDirichlet, as_vector( [sin(2*pi*x[0])] )
+
 def cosProblem():
-    return LinearAdvectionDiffusion1D, as_vector( [cos(2*pi*x[0])] )
+    return LinearAdvectionDiffusion1DNeuman, as_vector( [cos(2*pi*x[0])] )
 
 def burgersShock():
     return Burgers1D, riemanProblem(x[0],-0.5,[1],[0])

pydemo/scalar/shock_tube.py

@@ -8,8 +8,9 @@ from dune.femdg import createFemDGSolver
 # from scalar import riemanProblem as problem
 # from scalar import constantTransport as probelm
 # from scalar import cosProblem as problem
+from scalar import cosProblem as problem
 # from scalar import burgersShock as problem
-from scalar import burgersVW as problem
+# from scalar import burgersVW as problem
 # from scalar import burgersStationaryShock as problem
 
 Model, initial = problem()
@@ -18,7 +19,7 @@ Model, initial = problem()
 parameter.append("parameter")
 parameter.append({"fem.verboserank": -1})
 
-x0,x1,N = [-1, 0], [1, 0.1], [20, 5]
+x0,x1,N = [-1, 0], [1, 0.1], [50, 7]
 grid = structuredGrid(x0,x1,N)
 # grid = create.grid("ALUSimplex", cartesianDomain(x0,x1,N))
 dimR     = 1

python/dune/femdg/_operators.py

@@ -153,7 +153,7 @@ def createFemDGSolver(Model, space):
 
     hasAdvection = hasattr(Model,"F_c")
     if hasAdvection:
-        advModel = inner(as_matrix(Model.F_c(u)),grad(v))*dx
+        advModel = inner(Model.F_c(u),grad(v))*dx
     else:
         advModel = inner(grad(u-u),grad(v))*dx    # TODO: make a better empty model
     if hasattr(Model,"S_ns"):