# Solve Poisson problem on unit square using "primal" discontinuous Petrov-Galerkin method

from dolfin import *

# Test case
u_exact = Expression("sin(pi*x[0])*sin(2*pi*x[1])")
f_exact = Expression("sin(pi*x[0])*sin(2*pi*x[1])*5*pi*pi")


# Create mesh and define function space
degree  = 1
delta_p = 2
N = 40
mesh = UnitSquareMesh(N,N)

u0 = Constant(0.0)

# define spaces
U = FunctionSpace(mesh, "CG", degree)					    # trial functions
S = FunctionSpace(mesh, "RT", degree, restriction="facet") 	# fluxes
V = FunctionSpace(mesh, "DG", degree+delta_p)               # discontinuous test functions
E = MixedFunctionSpace([U,S,V])
(u,sig,v)    = TrialFunctions(E)
(du,dsig,dv) = TestFunctions(E)
n = FacetNormal(mesh)

def inner_v(v1,v2):
	return inner(grad(v1),grad(v2))*dx + v1*v2*dx

def b(u,s,v):
	return inner(grad(u),grad(v))*dx - inner(dot(s('+'),n('+')),jump(v))*dS - inner(dot(s,n),v)*ds

a = inner_v(v,dv) + b(u,sig,dv) + b(du,dsig,v)
L = inner(f_exact,dv)*dx

bcs = []
bcs.append(DirichletBC(E.sub(0), u_exact, DomainBoundary())) # boundary conditions on u

uSol = Function(E)
solve(a==L, uSol, bcs)
(u_,sig_,v_) = uSol.split()

# Estimate the error
L2_err = sqrt(assemble((u_-u_exact)*(u_-u_exact)*dx))
V_err  = sqrt(assemble(inner_v(v_,v_)))
print 'L2 error = ', L2_err, ' V-error = ', V_err



# Save the solution
u_file = File("u.pvd",'compressed')
u_file << u_
s_file = File("sigma.pvd",'compressed')
s_file << sig_
v_file = File("v.pvd",'compressed')
v_file << v_