from fenics import *
from mshr import *
import numpy as np

finalt = 5 #Final time
steps = 5000 #number of steps
dt = finalt/steps #time steps

# model parameters
mu = 0.001 #dynamic viscosity
rho = 1 #density

#creating mesh
channel = Rectangle(Point(0,0), Point(2.2, 0.41))
cylinder = Circle(Point(0.2,0.2), 0.05)
domain = channel - cylinder
mesh = generate_mesh(domain,64)

#Defining spaces
V = VectorFunctionSpace(mesh,'P', 2)
Q = FunctionSpace(mesh, 'P', 1)

#defining boundaries:
inflow = 'near(x[0], 0)'
outflow = 'near(x[0], 2.2)'
walls = 'near(x[1], 0) || near(x[1], 0.41)'
cylinder = 'on_boundary && x[0]>0.1 && x[0]<0.3 &&x[1]>0.1 && x[1]<0.3'

#defining inflow profile
inflowProfile('4.0*1.5*x[1]*(0.41-x[1])/pow(0.41,2)', '0')

#defining boundary conditions:

bcuInflow = DirichletBC(V,Expression(inflowProfile, degree=2), inflow)
bcuWalls = DirichletBC(V, Constant((0, 0)), walls)
bcuCylinder = DirichletBC(V, Constant((0,0)), cylinder)
bcpOutflow = DirichletBC(Q, Constant(0), outflow)
bcu = [bcuInflow, bcuWalls, bcuCylinder]
bcp = [bcpOutflow]