~fenics-core/dolfin/remove-init-mesh

« back to all changes in this revision

Viewing changes to demo/undocumented/lift-drag2/python/demo_lift-drag.py

  • Committer: Garth N. Wells
  • Date: 2013-01-08 17:02:41 UTC
  • Revision ID: gnw20@cam.ac.uk-20130108170241-lgl9nr9ppq9v4y1w
Remove a demo

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
"""This demo demonstrates how to compute functionals (or forms in
 
2
general) over subsets of the mesh. The two functionals lift and
 
3
drag are computed for the pressure field around a dolphin. Here, we
 
4
use the pressure field obtained from solving the Stokes equations
 
5
(see demo program in the sub directory
 
6
src/demo/pde/stokes/taylor-hood).
 
7
 
 
8
The calculation only includes the pressure contribution (not shear
 
9
forces).
 
10
"""
 
11
 
 
12
# Copyright (C) 2007 Kristian B. Oelgaard
 
13
#
 
14
# This file is part of DOLFIN.
 
15
#
 
16
# DOLFIN is free software: you can redistribute it and/or modify
 
17
# it under the terms of the GNU Lesser General Public License as published by
 
18
# the Free Software Foundation, either version 3 of the License, or
 
19
# (at your option) any later version.
 
20
#
 
21
# DOLFIN is distributed in the hope that it will be useful,
 
22
# but WITHOUT ANY WARRANTY; without even the implied warranty of
 
23
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 
24
# GNU Lesser General Public License for more details.
 
25
#
 
26
# You should have received a copy of the GNU Lesser General Public License
 
27
# along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
 
28
#
 
29
# Modified by Anders Logg, 2008.
 
30
# Modified by Garth N. Wells, 2009.
 
31
#
 
32
# First added:  2007-11-14
 
33
# Last changed: 2008-12-27
 
34
 
 
35
from dolfin import *
 
36
 
 
37
# Read the mesh from file
 
38
mesh =  Mesh("../mesh.xml.gz")
 
39
 
 
40
# Create FunctionSpace for pressure field
 
41
Vp = FunctionSpace(mesh, "CG", 1)
 
42
p = Function(Vp, "../pressure.xml.gz")
 
43
 
 
44
# Define sub domain for the dolphin
 
45
class Fish(SubDomain):
 
46
    def inside(self, x, on_boundary):
 
47
        return x[0] > DOLFIN_EPS and x[0] < (1.0 - DOLFIN_EPS) and \
 
48
               x[1] > DOLFIN_EPS and x[1] < (1.0 - DOLFIN_EPS)
 
49
 
 
50
# Define functionals for drag and lift
 
51
n = FacetNormal(mesh)
 
52
D = -p*n[0]*ds
 
53
L = p*n[1]*ds
 
54
 
 
55
# Assemble functionals over sub domain
 
56
fish = Fish()
 
57
drag = assemble(D, mesh=mesh, exterior_facet_domains=fish)
 
58
lift = assemble(L, mesh=mesh, exterior_facet_domains=fish)
 
59
 
 
60
print "Lift: %f" %lift
 
61
print "Drag: %f" %drag