← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to demo/pde/stokes/stabilized/python/demo.py

  • Committer: Johannes Ring
  • Date: 2008-03-05 22:43:06 UTC
  • Revision ID: johannr@simula.no-20080305224306-2npsdyhfdpl2esji
The BIG commit!

Show diffs side-by-side

added added

removed removed

Lines of Context:
demo/pde/stokes/stabilized/python/demo.py
 
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.
 
6
#
 
7
# Original implementation: ../cpp/main.cpp by Anders Logg
 
8
#
 
9
 
 
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"
 
14
 
 
15
from dolfin import *
 
16
 
 
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
 
22
 
 
23
# Load subdomains
 
24
sub_domains = MeshFunction("uint", mesh, "../subdomains.xml.gz");
 
25
 
 
26
# Function for no-slip boundary condition for velocity
 
27
class Noslip(Function):
 
28
    def __init__(self, mesh):
 
29
        Function.__init__(self, mesh)
 
30
 
 
31
    def eval(self, values, x):
 
32
        values[0] = 0.0
 
33
        values[1] = 0.0
 
34
 
 
35
    def rank(self):
 
36
        return 1
 
37
 
 
38
    def dim(self, i):
 
39
        return 2
 
40
 
 
41
# Function for inflow boundary condition for velocity
 
42
class Inflow(Function):
 
43
    def __init__(self, mesh):
 
44
        Function.__init__(self, mesh)
 
45
 
 
46
    def eval(self, values, x):
 
47
        values[0] = -1.0
 
48
        values[1] = 0.0
 
49
 
 
50
    def rank(self):
 
51
        return 1
 
52
 
 
53
    def dim(self, i):
 
54
        return 2
 
55
 
 
56
# Create functions for boundary conditions
 
57
noslip = Noslip(mesh)
 
58
inflow = Inflow(mesh)
 
59
zero = Function(mesh, 0.0)
 
60
 
 
61
# Define sub systems for boundary conditions
 
62
velocity = SubSystem(0)
 
63
pressure = SubSystem(1)
 
64
 
 
65
# No-slip boundary condition for velocity
 
66
bc0 = DirichletBC(noslip, sub_domains, 0, velocity)
 
67
 
 
68
# Inflow boundary condition for velocity
 
69
bc1 = DirichletBC(inflow, sub_domains, 1, velocity)
 
70
 
 
71
# Boundary condition for pressure at outflow
 
72
bc2 = DirichletBC(zero, sub_domains, 2, pressure)
 
73
 
 
74
# Collect boundary conditions
 
75
bcs = [bc0, bc1, bc2]
 
76
 
 
77
 
 
78
# Define variational problem
 
79
(v, q) = TestFunctions(system)
 
80
(u, p) = TrialFunctions(system)
 
81
 
 
82
f = Function(vector, mesh, 0.0)
 
83
h = MeshSize("triangle", mesh)
 
84
 
 
85
beta  = 0.2
 
86
delta = beta*h*h
 
87
 
 
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
 
90
 
 
91
# Set up PDE
 
92
pde = LinearPDE(a, L, mesh, bcs)
 
93
 
 
94
# Solve PDE
 
95
(U, P) = pde.solve().split()
 
96
 
 
97
# Plot solution
 
98
plot(U)
 
99
plot(P)
 
100
 
 
101
# Save solution
 
102
ufile = File("velocity.xml")
 
103
ufile << U
 
104
pfile = File("pressure.xml")
 
105
pfile << P
 
106
 
 
107
# Save solution in VTK format
 
108
ufile_pvd = File("velocity.pvd")
 
109
ufile_pvd << U
 
110
pfile_pvd = File("pressure.pvd")
 
111
pfile_pvd << P
 
112