Dear all,
As an introduction to blood flow modelling I tried to start with Poisson equation with different boundary conditions and subdomains.
The code compiles but I have a feeling that it doesn't work properly.
Especially, that from FEniCS I get:
*** Warning: Found no facets matching domain for boundary condition.
I would be grateful if you could look at the code and tell me if you see any mistakes
from fenics import*
from mshr import*
from math import pi, sin, cos, tan, pi, fabs, sqrt
Angle of the bypass
alpha = pi/4
Radius of the blood vessel
r1 = 1.5
Radiuss of the bypass
r2 = 0.5
Length of the blood vessel
l1 = 10
Length of the bypass
l2 = 6
Check if dimensions are correct
condition0 = (alpha < pi/2)
condition1 = (cos(alpha)r2 < r1)
condition2 = (sin(alpha)r1 < l1/2)
if (condition0) and (condition1) and (condition2):
print 'Dimensions are correct'
else:
print 'Dimensions are incorrect'
Generate main blood vessel
W1 = Cylinder(Point(-l1/2, 0, 0), Point(l1/2, 0, 0), r1, r1, 64)
Generate bypass
W2 = Cylinder(Point(0, 0, 0), Point(-cos(alpha)l2, sin(alpha)l2, 0), r2, r2, 64)
Define domain
domain = W1 + W2
Define left blood vessel boundary
class Left(SubDomain):
def inside(self, x, on_boundary):
return ((x[0] + l1/2) < DOLFIN_EPS) and on_boundary
Define right blood vessel boundary
class Right(SubDomain):
def inside(self, x, on_boundary):
return ((x[0] - l1/2) < DOLFIN_EPS) and on_boundary
Define upper bypass boundary
class Bypass(SubDomain):
def inside(self, x, on_boundary):
return (((fabs(x[0]/tan(alpha) - x[1] + (cos(alpha)l2)/tan(alpha) + \
sin(alpha)l2))/sqrt((1/tan(alpha))**2 + 1)) < DOLFIN_EPS) \
and on_boundary
Initialize sub-domain instances
left = Left()
right = Right()
bypass = Bypass()
Generate mesh
mesh = generate_mesh(domain, 64)
Mark domains
domains = CellFunction("size_t", mesh)
domains.set_all(0)
Mark boundaries
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)
bypass.mark(boundaries, 3)
Define function space
V = FunctionSpace(mesh, 'P', 2)
Define Dirichlet boundary
bcs = [DirichletBC(V, 2.0, boundaries, 0), DirichletBC(V, 3.0, boundaries, 1) ]
Definiowanie measures
dx = Measure('dx', domain = mesh, subdomain_data = domains)
ds = Measure('ds', domain = mesh, subdomain_data = boundaries)
Define trial function
u = TrialFunction(V)
Define test function
v = TestFunction(V)
Define value of Laplace operator
f = Constant(1.0)
Define value of Neumann boundary condition
g = Constant(0.0)
Define bilinear form
a = inner(grad(u), grad(v))*dx(0)
Define linear operator
L = fvdx(0) + gvds(2) + gvds(3)
Compute solution
u = Function(V)
solve(a == L, u, bcs)
Save solution in VTK format
vtkfile = File('Solutions/Poisson_subdomains.pvd')
vtkfile << u
