__author__ = "Anders Logg <logg@simula.no>"
__date__ = "2012-01-17"
__copyright__ = "Copyright (C) 2012 Anders Logg"
__license__  = "GNU LGPL version 3 or any later version"

# Last changed: 2012-01-25

from cbc.flow import *

eps = DOLFIN_EPS

class Inflow(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] < eps and x[1] > eps and x[1] < 1.0 - eps
        return on_boundary

class Outflow(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] > 1.0 - eps and x[1] > eps and x[1] < 1.0 - eps

class Noslip(SubDomain):
    def inside(self, x, on_boundary):
        inflow = x[0] < eps and x[1] > eps and x[1] < 1.0 - eps
        outflow = x[0] > 1.0 - eps and x[1] > eps and x[1] < 1.0 - eps
        return on_boundary and (not inflow) and (not outflow)

class Channel(NavierStokes):

    def mesh(self):
        return Mesh("dolfin_channel.xml")

    def density(self):
        return 1000.0

    def viscosity(self):
        # Viscosity of water / (kg / (m*s))
        return 0.001002

    def velocity_dirichlet_values(self):
        return [(0, 0)]

    def velocity_dirichlet_boundaries(self):
        return [Noslip()]

    def pressure_dirichlet_values(self):
        return [1000.0, 0]

    def pressure_dirichlet_boundaries(self):
        return [Inflow(), Outflow()]

    def velocity_initial_condition(self):
        return (0, 0)

    def pressure_initial_condition(self):
        return "0.0"

    def end_time(self):
        return 0.1

    def time_step(self):
        return 5e-4

# Solve problem
problem = Channel()
problem.parameters["solver_parameters"]["plot_solution"] = False
problem.parameters["solver_parameters"]["save_solution"] = True
u, p = problem.solve()

# Compute average velocity in x-direction
mesh = u.function_space().mesh()
velocity = assemble(u[0]*dx) / assemble(Constant(1)*dx, mesh=mesh)
print
print "Average x-velocity: %.16g" % velocity