from dolfin import *

""" This script solves the Poisson example 3."""

# Create mesh and 
mesh = Mesh("dolfin-channel.xml")

# Define function space
V = FunctionSpace(mesh, "Lagrange", 1)

# Define Dirichlet boundary x = 0
def left_boundary(x):
    return near(x[0],0.0)

# Define Dirichlet boundary x = 1
def right_boundary(x):
    return near(x[0],1.0)

# Define boundary condition
g_D0 = Expression("0.5*cos(2*DOLFIN_PI*x[1])")
bc0 = DirichletBC(V, g_D0, left_boundary)
g_D1 = Constant(1.0)
bc1 = DirichletBC(V, g_D1, right_boundary)
bcs = [bc0,bc1]

# Define trial an test functions
u = TrialFunction(V)
v = TestFunction(V)

# Define source term
f = Expression("2*cos(2*DOLFIN_PI*x[0])*cos(2*DOLFIN_PI*x[1])")
g_N = Expression("sin(DOLFIN_PI*x[0])*sin(DOLFIN_PI*x[1])")

# Define variational forms
a = inner(grad(u), grad(v))*dx
# Mind the minus sign here (opposed to the code snippet in the lecture 2"!
L = f*v*dx - g_N*v*ds

# Compute solution
u = Function(V)
solve(a == L, u, bcs)

# Save solution in VTK format
u_file = File("poisson_3.pvd")
u_file << u

# Plot solution
plot(u, interactive=True)