# Copyright (C) 2005 Johan Hoffman and Anders Logg.
# Licensed under the GNU GPL Version 2.
#
# Modified by Garth N. Wells 2005
#
# First added:  2005-04-04
# Last changed: 2005-12-31
#
# The bilinear form a(u,v) and linear form L(v) for
# convection-diffusion using cG(1)cG(1) with SUPG stabilisation.
#
# Compile this form with FFC: ffc ConvectionDiffusion.form.

scalar   = FiniteElement("Lagrange", "triangle", 1)
vector   = FiniteElement("Vector Lagrange", "triangle", 1)
constant = FiniteElement("Discontinuous Lagrange", "triangle", 0)

v    = BasisFunction(scalar) # test function
u1   = BasisFunction(scalar) # value at next time step
u0   = Function(scalar)      # value at previous time step
w    = Function(vector)      # convection
stab = Function(scalar)      # tau/2|w|
f    = Function(scalar)      # source term
h    = Function(constant)    # element size
k    = Constant()            # time step
c    = Constant()            # diffusion

beta = mult(h,stab)  # stabiliation term
r1   = dot(w, grad(u1)) - c*div(grad(u1)) # "steady" residual at next time step
r0   = dot(w, grad(u0)) - c*div(grad(u0)) # "steady" residual at previous time step
SUPG = beta*dot(w, grad(v))               # SUPG weighting

a_galerkin  = v*u1*dx + 0.5*k*(v*dot(w, grad(u1))*dx + c*dot( grad(v), grad(u1) )*dx) 
a_stabilise = SUPG*u1*dx + 0.5*k*SUPG*r1*dx

L_galerkin  = v*u0*dx - 0.5*k*(v*dot(w, grad(u0))*dx + c*dot(grad(v), grad(u0))*dx) + v*f*dx 
L_stabilise = SUPG*u0*dx - 0.5*k*SUPG*r0*dx + SUPG*f*dx

a = a_galerkin + a_stabilise
L = L_galerkin + L_stabilise