"""
Compute the 
"""

from dolfin import *

#
# Mesh
# 
mesh = Mesh("airfoil_exterior.xml")
V    = FunctionSpace(mesh,"CG",1)

#
# Boundary Conditions
#
potential_left = Constant(1.0)
def left_boundary(x,on_boundary):
  return abs(x[0]+2.0) < DOLFIN_EPS and on_boundary
bc_left = DirichletBC(V,potential_left,left_boundary)
  
potential_right = Constant(2.0)
def right_boundary(x,on_boundary):
  return abs(x[0]-12.0) < DOLFIN_EPS and on_boundary
bc_right = DirichletBC(V,potential_right,right_boundary)

bc = [bc_left,bc_right]

#
# Set up the bilinear form
# 
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u),grad(v))*dx
L = Constant(0.0)*v*dx

#
# Solve the system
# 
u = Function(V)
solve(a == L,u,bc)

file = File("02_potential_flow_airfoil.pvd")
file << u

u_vec = project(grad(u),VectorFunctionSpace(mesh,"CG",1))
file = File("02_potential_flow_airfoil_flux.pvd")
file << u_vec
plot(u_vec,interactive=True)