1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
|
"""This demo solves the equations of static linear elasticity for a
pulley subjected to centripetal accelerations. The solver uses
smoothed aggregation algerbaric multigrid."""
# Copyright (C) 2014 Garth N. Wells
#
# This file is part of DOLFIN.
#
# DOLFIN is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# DOLFIN is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
from dolfin import *
import matplotlib.pyplot as plt
# Test for PETSc
if not has_linear_algebra_backend("PETSc"):
print("DOLFIN has not been configured with PETSc. Exiting.")
exit()
# Set backend to PETSC
parameters["linear_algebra_backend"] = "PETSc"
def build_nullspace(V, x):
"""Function to build null space for 3D elasticity"""
# Create list of vectors for null space
nullspace_basis = [x.copy() for i in range(6)]
# Build translational null space basis
V.sub(0).dofmap().set(nullspace_basis[0], 1.0);
V.sub(1).dofmap().set(nullspace_basis[1], 1.0);
V.sub(2).dofmap().set(nullspace_basis[2], 1.0);
# Build rotational null space basis
V.sub(0).set_x(nullspace_basis[3], -1.0, 1);
V.sub(1).set_x(nullspace_basis[3], 1.0, 0);
V.sub(0).set_x(nullspace_basis[4], 1.0, 2);
V.sub(2).set_x(nullspace_basis[4], -1.0, 0);
V.sub(2).set_x(nullspace_basis[5], 1.0, 1);
V.sub(1).set_x(nullspace_basis[5], -1.0, 2);
for x in nullspace_basis:
x.apply("insert")
# Create vector space basis and orthogonalize
basis = VectorSpaceBasis(nullspace_basis)
basis.orthonormalize()
return basis
# Load mesh from file
mesh = Mesh()
XDMFFile(MPI.comm_world, "../pulley.xdmf").read(mesh)
# Function to mark inner surface of pulley
def inner_surface(x, on_boundary):
r = 3.75 - x[2]*0.17
return (x[0]*x[0] + x[1]*x[1]) < r*r and on_boundary
# Rotation rate and mass density
omega = 300.0
rho = 10.0
# Loading due to centripetal acceleration (rho*omega^2*x_i)
f = Expression(("rho*omega*omega*x[0]", "rho*omega*omega*x[1]", "0.0"),
omega=omega, rho=rho, degree=2)
# Elasticity parameters
E = 1.0e9
nu = 0.3
mu = E/(2.0*(1.0 + nu))
lmbda = E*nu/((1.0 + nu)*(1.0 - 2.0*nu))
# Stress computation
def sigma(v):
return 2.0*mu*sym(grad(v)) + lmbda*tr(sym(grad(v)))*Identity(len(v))
# Create function space
V = VectorFunctionSpace(mesh, "Lagrange", 1)
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
a = inner(sigma(u), grad(v))*dx
L = inner(f, v)*dx
# Set up boundary condition on inner surface
c = Constant((0.0, 0.0, 0.0))
bc = DirichletBC(V, c, inner_surface)
# Assemble system, applying boundary conditions and preserving
# symmetry)
A, b = assemble_system(a, L, bc)
# Create solution function
u = Function(V)
# Create near null space basis (required for smoothed aggregation
# AMG). The solution vector is passed so that it can be copied to
# generate compatible vectors for the nullspace.
null_space = build_nullspace(V, u.vector())
# Attach near nullspace to matrix
as_backend_type(A).set_near_nullspace(null_space)
# Create PETSC smoothed aggregation AMG preconditioner and attach near
# null space
pc = PETScPreconditioner("petsc_amg")
# Use Chebyshev smoothing for multigrid
PETScOptions.set("mg_levels_ksp_type", "chebyshev")
PETScOptions.set("mg_levels_pc_type", "jacobi")
# Improve estimate of eigenvalues for Chebyshev smoothing
PETScOptions.set("mg_levels_esteig_ksp_type", "cg")
PETScOptions.set("mg_levels_ksp_chebyshev_esteig_steps", 50)
# Create CG Krylov solver and turn convergence monitoring on
solver = PETScKrylovSolver("cg", pc)
solver.parameters["monitor_convergence"] = True
# Set matrix operator
solver.set_operator(A);
# Compute solution
solver.solve(u.vector(), b);
# Save solution to VTK format
File("elasticity.pvd", "compressed") << u
# Save colored mesh partitions in VTK format if running in parallel
if MPI.size(mesh.mpi_comm()) > 1:
File("partitions.pvd") << MeshFunction("size_t", mesh, mesh.topology().dim(), \
MPI.rank(mesh.mpi_comm()))
# Project and write stress field to post-processing file
W = TensorFunctionSpace(mesh, "Discontinuous Lagrange", 0)
stress = project(sigma(u), V=W)
File("stress.pvd") << stress
# Plot solution
plot(u)
plt.show()
|
