##########################################################
#FEniCS tutorial demo program: Poisson equation with Dirichlet conditions.
#Test problem is chosen to give an exact solution at all nodes of the mesh.
#  -Laplace(u) = f    in the unit square
#            u = u_D  on the boundary
#  u_D = 1 + x^2 + 2y^2
#    f = -6
###########################################################

module ft01

using FenicsPy

# Create mesh and define function space
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, "P", 1)

# Define boundary condition
u_D = Expression("1 + x[0]*x[0] + 2*x[1]*x[1]", degree=2)

bc = DirichletBC(V, u_D, "on_boundary")

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
a = dot(grad(u), grad(v))*dx
L = f*v*dx

# Compute solution
u = FeFunction(V)
solve(a == L, u, bc)

# Save solution to file in VTK format
vtkfile = File("poisson/solution.pvd")
vtkfile << u

# Compute error in L2 norm
error_L2 = errornorm(u_D, u, "L2")

# Compute maximum error at vertices
vertex_values_u_D = u_D.compute_vertex_values(mesh)
vertex_values_u = u.compute_vertex_values(mesh)
error_max = max(abs.(vertex_values_u_D - vertex_values_u)...)

# Print errors
println("error_L2  =", error_L2)
println("error_max =", error_max)

# Plot solution and mesh
#plot(u)
#plot(mesh)

end # module ft01
