1
# This demo solves the Stokes equations, using stabilized
2
# first order elements for the velocity and pressure. The
3
# sub domains for the different boundary conditions used
4
# in this simulation are computed by the demo program in
5
# src/demo/mesh/subdomains.
7
# Original implementation: ../cpp/main.cpp by Anders Logg
10
__author__ = "Kristian B. Oelgaard (k.b.oelgaard@tudelft.nl)"
11
__date__ = "2007-11-15 -- 2007-12-03"
12
__copyright__ = "Copyright (C) 2007 Kristian B. Oelgaard"
13
__license__ = "GNU LGPL Version 2.1"
17
# Load mesh and create finite elements
18
mesh = Mesh("../../../../../../data/meshes/dolfin-2.xml.gz")
19
scalar = FiniteElement("Lagrange", "triangle", 1)
20
vector = VectorElement("Lagrange", "triangle", 1)
21
system = vector + scalar
24
sub_domains = MeshFunction("uint", mesh, "../subdomains.xml.gz");
26
# Function for no-slip boundary condition for velocity
27
class Noslip(Function):
28
def __init__(self, mesh):
29
Function.__init__(self, mesh)
31
def eval(self, values, x):
41
# Function for inflow boundary condition for velocity
42
class Inflow(Function):
43
def __init__(self, mesh):
44
Function.__init__(self, mesh)
46
def eval(self, values, x):
56
# Create functions for boundary conditions
59
zero = Function(mesh, 0.0)
61
# Define sub systems for boundary conditions
62
velocity = SubSystem(0)
63
pressure = SubSystem(1)
65
# No-slip boundary condition for velocity
66
bc0 = DirichletBC(noslip, sub_domains, 0, velocity)
68
# Inflow boundary condition for velocity
69
bc1 = DirichletBC(inflow, sub_domains, 1, velocity)
71
# Boundary condition for pressure at outflow
72
bc2 = DirichletBC(zero, sub_domains, 2, pressure)
74
# Collect boundary conditions
78
# Define variational problem
79
(v, q) = TestFunctions(system)
80
(u, p) = TrialFunctions(system)
82
f = Function(vector, mesh, 0.0)
83
h = MeshSize("triangle", mesh)
88
a = (dot(grad(v), grad(u)) - div(v)*p + q*div(u) + delta*dot(grad(q), grad(p)))*dx
89
L = dot(v + mult(delta, grad(q)), f)*dx
92
pde = LinearPDE(a, L, mesh, bcs)
95
(U, P) = pde.solve().split()
102
ufile = File("velocity.xml")
104
pfile = File("pressure.xml")
107
# Save solution in VTK format
108
ufile_pvd = File("velocity.pvd")
110
pfile_pvd = File("pressure.pvd")