"""
Solve Laplace eqn using FENICs
"""

from dolfin import *

N = 64  # number of elements in the mesh
mesh = UnitSquareMesh(N,N)

# state space
Y = FunctionSpace(mesh,"Lagrange",1)
y_ = TrialFunction(Y)
v_ = TestFunction(Y)

# define the bilinear form
a = inner(grad(y_),grad(v_))*dx
# define the right hand side (degree=2 of interpolation)
f = Expression("sin(pi*x[0]) + cos(pi*x[1])", degree=2)
l = inner(f,v_)*dx

# define the Dirichlet boundary conditions
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        # in this example we mark all boundary as Dirichlet
        return on_boundary

y0 = Constant(0.0)
bc = DirichletBC(Y,y0,DirichletBoundary())

# assemble the system and apply boundary conditions
A = assemble(a)
b = assemble(l)
bc.apply(A,b)

# solve the system
y = Function(Y)
solve(A,y.vector(),b)

# plot the solution - does not work in Docker...
"""
plot(y,  title="Solution",  rescale=True)
# hold plot
interactive()
"""

# alternatively, one can also save the solution as VTK file
file_y = File("y.pvd")
file_y << y