from firedrake import *
m=20
tmesh = Mesh("trapezium.msh")


V = FunctionSpace(tmesh, "DG", 0)
W = VectorFunctionSpace(tmesh, "CG", 1)

timestep = 0.05




# Advection velocity
velocity = Expression(("1.0", " "1.0"))
u0 = project(velocity, W)
#Inflow BC
inflow = Expression("(x[1] < 0.2) && (x[0] > 0.5) ? 1.0 : 0")

ic = Expression("(x[1] = 0.0 ) && (x[0] > 0.5) ? 1.0 : 0")
D0 = Function(V)
D0.interpolate(inflow)

#Set up element normal
n = FacetNormal(mesh)
#
un = 0.5*(dot(u0, n) + abs(dot(u0, n)))

#Previous timestep for BE timestep
D_=Function(V)

D_.interpolate(ic)

D = TrialFunction(V)
phi = TestFunction(V)


a1 = ((D/timestep)*phi - D*dot(u0, grad(phi)))*dx
a2 = dot(jump(phi), un('+')*D('+') - un('-')*D('-'))*dS
a3 = dot(phi, un*D)*(ds("Right")+ds("Top"))   # outflow at top wall
a = a1 + a2 + a3

L = (D_/timestep)*phi*dx - D0*phi*dot(u0, n)*(ds("Bottom")+ds("Left"))  # inflow at bottom wall


out = Function(V)




outfile = File("upwind_unsteady.pvd")

outfile << D_
t=0.0
end=1
while (t <= end):
       solve(a == L, out)
       D_.assign(out)
       t+=timestep
       outfile << out