← Back to branch summary

~ubuntu-branches/ubuntu/wily/dolfin/wily-proposed

« back to all changes in this revision

Viewing changes to demo/undocumented/elasticity/python/demo_elasticity.py

  • Committer: Package Import Robot
  • Author(s): Johannes Ring
  • Date: 2015-03-17 07:57:11 UTC
  • mfrom: (1.1.18) (19.1.24 experimental)
  • Revision ID: package-import@ubuntu.com-20150317075711-1v207zbty9qmygow
Tags: 1.5.0-1
http://bugs.debian.org/780359
* New upstream release (closes: #780359).
* debian/control:
  - Bump Standards-Version to 3.9.6 (no changes needed).
  - Bump X-Python-Version to >= 2.7.
  - Update package names for new SONAME 1.5 (libdolfin1.4 ->
    libdolfin1.5, libdolfin1.4-dbg -> libdolfin1.5-dbg and
    libdolfin1.4-dev -> libdolfin1.5-dev).
  - Bump minimum required version for python-instant, python-ufl and
    python-ffc to 1.5.0.
  - Add python-sympy and python-six to Depends for binary package
    python-dolfin.
  - Add dh-python to Build-Depends.
  - Remove libcgal-dev from {Build-}Depends.
* Remove CSGCGALMeshGenerator3D-oom.patch since CGAL is no longer used
  by DOLFIN.
* Move debian/libdolfin1.4.install -> debian/libdolfin1.5.install.
* debian/rules: No longer any non DFSG-free stuff to remove, so update
  get-orig-source target (update debian/watch accordingly).
* Update debian/copyright and debian/copyright_hints.

Show diffs side-by-side

added added

removed removed

Lines of Context:
demo/undocumented/elasticity/python/demo_elasticity.py
1
 
"""This demo program solves the equations of static linear elasticity
2
 
for a gear clamped at two of its ends and twisted 30 degrees."""
 
1
"""This demo solves the equations of static linear elasticity for a
 
2
pulley subjected to centripetal accelerations. The solver uses
 
3
smoothed aggregation algerbaric multigrid."""
3
4
 
4
 
# Copyright (C) 2007 Kristian B. Oelgaard
 
5
# Copyright (C) 2014 Garth N. Wells
5
6
#
6
7
# This file is part of DOLFIN.
7
8
#
17
18
#
18
19
# You should have received a copy of the GNU Lesser General Public License
19
20
# along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
20
 
#
21
 
# Modified by Anders Logg 2008-2011
22
 
#
23
 
# First added:  2007-11-14
24
 
# Last changed: 2011-06-28
25
21
 
 
22
from __future__ import print_function
26
23
from dolfin import *
27
24
 
 
25
# Test for PETSc
 
26
if not has_linear_algebra_backend("PETSc"):
 
27
    print("DOLFIN has not been configured with PETSc. Exiting.")
 
28
    exit()
 
29
 
 
30
# Set backend to PETSC
 
31
parameters["linear_algebra_backend"] = "PETSc"
 
32
 
 
33
def build_nullspace(V, x):
 
34
    """Function to build null space for 3D elasticity"""
 
35
 
 
36
    # Create list of vectors for null space
 
37
    nullspace_basis = [x.copy() for i in range(6)]
 
38
 
 
39
    # Build translational null space basis
 
40
    V.sub(0).dofmap().set(nullspace_basis[0], 1.0);
 
41
    V.sub(1).dofmap().set(nullspace_basis[1], 1.0);
 
42
    V.sub(2).dofmap().set(nullspace_basis[2], 1.0);
 
43
 
 
44
    # Build rotational null space basis
 
45
    V.sub(0).dofmap().set_x(nullspace_basis[3], -1.0, 1, V.mesh());
 
46
    V.sub(1).dofmap().set_x(nullspace_basis[3],  1.0, 0, V.mesh());
 
47
    V.sub(0).dofmap().set_x(nullspace_basis[4],  1.0, 2, V.mesh());
 
48
    V.sub(2).dofmap().set_x(nullspace_basis[4], -1.0, 0, V.mesh());
 
49
    V.sub(2).dofmap().set_x(nullspace_basis[5],  1.0, 1, V.mesh());
 
50
    V.sub(1).dofmap().set_x(nullspace_basis[5], -1.0, 2, V.mesh());
 
51
 
 
52
    for x in nullspace_basis:
 
53
        x.apply("insert")
 
54
 
 
55
    return VectorSpaceBasis(nullspace_basis)
 
56
 
 
57
 
28
58
# Load mesh and define function space
29
 
mesh = Mesh("../gear.xml.gz")
30
 
mesh.order()
31
 
V = VectorFunctionSpace(mesh, "CG", 1)
32
 
 
33
 
# Sub domain for clamp at left end
34
 
def left(x, on_boundary):
35
 
    return x[0] < 0.5 and on_boundary
36
 
 
37
 
# Dirichlet boundary condition for rotation at right end
38
 
class Rotation(Expression):
39
 
 
40
 
    def eval(self, values, x):
41
 
 
42
 
        # Center of rotation
43
 
        y0 = 0.5
44
 
        z0 = 0.219
45
 
 
46
 
        # Angle of rotation (30 degrees)
47
 
        theta = 0.5236
48
 
 
49
 
        # New coordinates
50
 
        y = y0 + (x[1] - y0)*cos(theta) - (x[2] - z0)*sin(theta)
51
 
        z = z0 + (x[1] - y0)*sin(theta) + (x[2] - z0)*cos(theta)
52
 
 
53
 
        # Clamp at right end
54
 
        values[0] = 0.0
55
 
        values[1] = y - x[1]
56
 
        values[2] = z - x[2]
57
 
 
58
 
    def value_shape(self):
59
 
        return (3,)
60
 
 
61
 
# Sub domain for rotation at right end
62
 
def right(x, on_boundary):
63
 
    return x[0] > 0.9 and on_boundary
 
59
mesh = Mesh("../pulley.xml.gz")
 
60
 
 
61
# Function to mark inner surface of pulley
 
62
def inner_surface(x, on_boundary):
 
63
    r = 3.75 - x[2]*0.17
 
64
    return (x[0]*x[0] + x[1]*x[1]) < r*r and on_boundary
 
65
 
 
66
# Rotation rate and mass density
 
67
omega = 300.0
 
68
rho = 10.0
 
69
 
 
70
# Loading due to centripetal acceleration (rho*omega^2*x_i)
 
71
f = Expression(("rho*omega*omega*x[0]", \
 
72
                "rho*omega*omega*x[1]", \
 
73
                "0.0"), omega=omega, rho=rho)
 
74
 
 
75
# Elasticity parameters
 
76
E = 1.0e9
 
77
nu = 0.3
 
78
mu = E/(2.0*(1.0 + nu))
 
79
lmbda = E*nu/((1.0 + nu)*(1.0 - 2.0*nu))
 
80
 
 
81
# Stress computation
 
82
def sigma(v):
 
83
    gdim = v.geometric_dimension()
 
84
    return 2.0*mu*sym(grad(v)) + lmbda*tr(sym(grad(v)))*Identity(gdim)
 
85
 
 
86
# Create function space
 
87
V = VectorFunctionSpace(mesh, "Lagrange", 1)
64
88
 
65
89
# Define variational problem
66
90
u = TrialFunction(V)
67
91
v = TestFunction(V)
68
 
f = Constant((0.0, 0.0, 0.0))
69
 
 
70
 
E  = 10.0
71
 
nu = 0.3
72
 
 
73
 
mu    = E / (2.0*(1.0 + nu))
74
 
lmbda = E*nu / ((1.0 + nu)*(1.0 - 2.0*nu))
75
 
 
76
 
def sigma(v):
77
 
    return 2.0*mu*sym(grad(v)) + lmbda*tr(sym(grad(v)))*Identity(v.geometric_dimension())
78
 
 
79
92
a = inner(sigma(u), grad(v))*dx
80
93
L = inner(f, v)*dx
81
94
 
82
 
# Set up boundary condition at left end
 
95
# Set up boundary condition on inner surface
83
96
c = Constant((0.0, 0.0, 0.0))
84
 
bcl = DirichletBC(V, c, left)
85
 
 
86
 
# Set up boundary condition at right end
87
 
r = Rotation()
88
 
bcr = DirichletBC(V, r, right)
89
 
 
90
 
# Set up boundary conditions
91
 
bcs = [bcl, bcr]
 
97
bc = DirichletBC(V, c, inner_surface)
 
98
 
 
99
# Assemble system, applying boundary conditions and preserving
 
100
# symmetry)
 
101
A, b = assemble_system(a, L, bc)
 
102
 
 
103
# Create solution function
 
104
u = Function(V)
 
105
 
 
106
# Create near null space basis (required for smoothed aggregation
 
107
# AMG). The solution vector is passed so that it can be copied to
 
108
# generate compatible vectors for the nullspace.
 
109
null_space = build_nullspace(V, u.vector())
 
110
 
 
111
# Create PETSC smoothed aggregation AMG preconditioner and attach near
 
112
# null space
 
113
pc = PETScPreconditioner("petsc_amg")
 
114
pc.set_nullspace(null_space)
 
115
 
 
116
# Create CG Krylov solver and turn convergence monitoring on
 
117
solver = PETScKrylovSolver("cg", pc)
 
118
solver.parameters["monitor_convergence"] = True
 
119
 
 
120
# Set matrix operator
 
121
solver.set_operator(A);
92
122
 
93
123
# Compute solution
94
 
u = Function(V)
95
 
solve(a == L, u, bcs, solver_parameters={"symmetric": True})
 
124
solver.solve(u.vector(), b);
96
125
 
97
126
# Save solution to VTK format
98
127
File("elasticity.pvd", "compressed") << u
105
134
# Project and write stress field to post-processing file
106
135
W = TensorFunctionSpace(mesh, "Discontinuous Lagrange", 0)
107
136
stress = project(sigma(u), V=W)
108
 
File("streess.pvd") << stress
 
137
File("stress.pvd") << stress
109
138
 
110
139
# Plot solution
111
140
plot(u, interactive=True)