← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to demo/pde/elasticity/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/elasticity/python/demo.py
 
1
# This demo program solves the equations of static
 
2
# linear elasticity for a gear clamped at two of its
 
3
# ends and twisted 30 degrees.
 
4
#
 
5
# Original implementation: ../cpp/main.cpp by Johan Jansson and Anders Logg
 
6
#
 
7
__author__ = "Kristian B. Oelgaard (k.b.oelgaard@tudelft.nl)"
 
8
__date__ = "2007-11-14 -- 2007-12-03"
 
9
__copyright__ = "Copyright (C) 2007 Kristian B. Oelgaard"
 
10
__license__  = "GNU LGPL Version 2.1"
 
11
 
 
12
from dolfin import *
 
13
 
 
14
# Load mesh and create finite element
 
15
mesh = Mesh("../../../../../data/meshes/gear.xml.gz")
 
16
element = VectorElement("Lagrange", "tetrahedron", 1)
 
17
 
 
18
# Dirichlet boundary condition for clamp at left end
 
19
class Clamp(Function):
 
20
    def __init__(self, element, mesh):
 
21
        Function.__init__(self, element, mesh)
 
22
 
 
23
    def eval(self, values, x):
 
24
        values[0] = 0.0
 
25
        values[1] = 0.0
 
26
        values[2] = 0.0
 
27
 
 
28
    def rank(self):
 
29
        return 1
 
30
 
 
31
    def dim(self, i):
 
32
        return 3
 
33
 
 
34
# Sub domain for clamp at left end
 
35
class Left(SubDomain):
 
36
    def inside(self, x, on_boundary):
 
37
        return bool(x[0] < 0.5 and on_boundary)
 
38
 
 
39
# Dirichlet boundary condition for rotation at right end
 
40
class Rotation(Function):
 
41
    def __init__(self, element, mesh):
 
42
        Function.__init__(self, element, mesh)
 
43
 
 
44
    def eval(self, values, x):
 
45
        # Center of rotation
 
46
        y0 = 0.5
 
47
        z0 = 0.219
 
48
 
 
49
        # Angle of rotation (30 degrees)
 
50
        theta = 0.5236
 
51
 
 
52
        # New coordinates
 
53
        y = y0 + (x[1] - y0)*cos(theta) - (x[2] - z0)*sin(theta)
 
54
        z = z0 + (x[1] - y0)*sin(theta) + (x[2] - z0)*cos(theta)
 
55
 
 
56
        # Clamp at right end
 
57
        values[0] = 0.0
 
58
        values[1] = y - x[1]
 
59
        values[2] = z - x[2]
 
60
 
 
61
    def rank(self):
 
62
        return 1
 
63
 
 
64
    def dim(self, i):
 
65
        return 3
 
66
 
 
67
# Sub domain for rotation at right end
 
68
class Right(SubDomain):
 
69
    def inside(self, x, on_boundary):
 
70
      return bool(x[0] > 0.9 and on_boundary)
 
71
 
 
72
# Initialise source function
 
73
f = Function(element, mesh, 0.0)
 
74
 
 
75
# Define variational problem
 
76
# Test and trial functions
 
77
v = TestFunction(element)
 
78
u = TrialFunction(element)
 
79
 
 
80
E  = 10.0
 
81
nu = 0.3
 
82
 
 
83
mu    = E / (2*(1 + nu))
 
84
lmbda = E*nu / ((1 + nu)*(1 - 2*nu))
 
85
 
 
86
def epsilon(v):
 
87
    return 0.5*(grad(v) + transp(grad(v)))
 
88
 
 
89
def sigma(v):
 
90
    return 2.0*mu*epsilon(v) + lmbda*mult(trace(epsilon(v)), Identity(len(v)))
 
91
 
 
92
a = dot(grad(v), sigma(u))*dx
 
93
L = dot(v, f)*dx
 
94
 
 
95
# Set up boundary condition at left end
 
96
c = Clamp(element, mesh)
 
97
left = Left()
 
98
bcl = DirichletBC(c, mesh, left)
 
99
 
 
100
# Set up boundary condition at right end
 
101
r = Rotation(element, mesh)
 
102
right = Right()
 
103
bcr = DirichletBC(r, mesh, right)
 
104
 
 
105
# Set up boundary conditions
 
106
bcs = [bcl, bcr]
 
107
 
 
108
# Set up PDE and solve
 
109
pde = LinearPDE(a, L, mesh, bcs)
 
110
sol = pde.solve()
 
111
 
 
112
# Plot solution
 
113
plot(sol, mode="displacement")
 
114
 
 
115
# Save solution to VTK format
 
116
vtk_file = File("elasticity.pvd")
 
117
vtk_file << sol
 
118
 
 
119
# Save solution to XML format
 
120
xml_file = File("elasticity.xml")
 
121
xml_file << sol
 
122