
from dolfin import *

mesh = Mesh("model_liten_rett.xml.gz")

V = FunctionSpace(mesh, "CG", 1)

u = TrialFunction(V)
v = TestFunction(V)


f = Constant(0)

a = u*v*dx 
L = f*v*dx 

inflow_boundary = 'liten_inflow_bd_r.xml'
outflow_boundary = 'liten_outflow_bd_r.xml'
noslip_boundary = 'liten_no_slip_bd_r.xml'

	
sub_domain_no_slip =  MeshFunction("uint", self.mesh, noslip_boundary)
sub_domain_inflow =  MeshFunction("uint", self.mesh, inflow_boundary)
sub_domain_outflow =  MeshFunction("uint", self.mesh, outflow_boundary)

# Create no-slip boundary condition for velocity
#bv = DirichletBC(V, Constant((0.0, 0.0, 0.0)), sub_domain_no_slip, 3)

pressure_bc = Expression("(100-x[2])*b")
pressure_bc.b = 0.1/20.0

# Create boundary conditions for pressure
bp0 = DirichletBC(V, pressure_bc, sub_domain_inflow, 3)
bp1 = DirichletBC(V, pressure_bc, sub_domain_outflow, 3)

bc = [bp0, bp1]

problem = VariationalProblem(a, L, bc)
u = problem.solve()

plot(u)
interactive()
