"""
Solve
  - div(grad(p)) = -div(rho g(x,y)),
subject to the boundary conditions
  p = 0                for x in \Gamma_1.
  
Assumes: constant density; variable gravity vector.
"""

from dolfin import *

# Define radial gravity vector
gravity = Expression(("-3.3*x[0]/sqrt(x[0]*x[0]+x[1]*x[1])", "-3.3*x[1]/sqrt(x[0]*x[0]+x[1]*x[1])"), degree=2)

# Create mesh and define function space
mesh = Mesh("circle.xml")
boundaries = MeshFunction('size_t', mesh, "circle_facet_region.xml")

# Define constant density
rho = Constant(1.0)

# Define spaces
V = FunctionSpace(mesh, "Lagrange", 1)
q, p = TestFunction(V), TrialFunction(V)

# Define Dirichlet boundaries
pressure_reference_value = 0.0
bc = DirichletBC(V, pressure_reference_value, boundaries, 1)

# Variational problem
bform_A = inner(grad(q), grad(p)) * dx
lform_L = inner(grad(q), rho*gravity) * dx

# Compute solution
p = Function(V)
solve(bform_A == lform_L, p, bc)

# Save solution
file = File("circle_p.pvd")
file << p

