/
stokes2.py
50 lines (42 loc) · 1.22 KB
/
stokes2.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
from dolfin import *
from mshr import *
b = 0.7
embryo = Ellipse(Point(0.0, 0.0), 1, b)
mesh = generate_mesh(embryo, 32)
# Define function spaces
P2 = VectorElement("CG", triangle, 2)
P1 = FiniteElement("CG", triangle, 1)
TH = MixedElement([P2, P1])
W = FunctionSpace(mesh, TH)
g = Constant(0.0)
mu = Constant(1.0)
force = Constant((0.0, 0.0))
# Specify Boundary Conditions
flow_profile = ("-sin(atan2(x[1]/b, x[0]))*sin(nharmonic*atan2(x[1]/b, x[0]))",
"+cos(atan2(x[1]/b, x[0]))*sin(nharmonic*atan2(x[1]/b, x[0]))")
bc = DirichletBC(W.sub(0),
Expression(flow_profile, degree=2, b=b, nharmonic=2),
"on_boundary")
# Define trial and test functions
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
def epsilon(u):
return grad(u) + nabla_grad(u)
a = inner(mu*epsilon(u) + p*Identity(2), epsilon(v))*dx -div(u)*q*dx -1e-10*p*q*dx
L = dot(force, v)*dx + g*q*dx
# Solve system
U = Function(W)
solve(a == L, U, bc)
# Get sub-functions
u, p = U.split()
from vedo.dolfin import plot
plot(
u,
mode='mesh and arrows',
scale=0.1,
warp_zfactor=-0.1,
lw=0,
scalarbar='horizontal',
axes={'xlabel_size':0.01,'ylabel_size':0.01, 'ztitle':''},
title="Velocity",
)