I installed Fenics using Conda on Ubuntu 18.04, and receive the following error while running their ft06_elasticity.py example.
I have tried to find a solution or workaround to this in the documentation, but I cannot even find a nabla_div() function description anywhere.
The Fenics documentation states the following:
nabla_grad
The gradient and divergence operators now have a prefix nabla_. This is strictly not necessary in the present problem, but recommended in general for vector PDEs arising from continuum mechanics, if you interpret ∇ as a vector in the PDE notation; see the box about nabla_grad in the section Variational formulation.
"""
FEniCS tutorial demo program: Linear elastic problem.
-div(sigma(u)) = f
The model is used to simulate an elastic beam clamped at
its left end and deformed under its own weight.
"""
from __future__ import print_function
from fenics import *
# Scaled variables
L = 1; W = 0.2
mu = 1
rho = 1
delta = W/L
gamma = 0.4*delta**2
beta = 1.25
lambda_ = beta
g = gamma
# Create mesh and define function space
mesh = BoxMesh(Point(0, 0, 0), Point(L, W, W), 10, 3, 3)
V = VectorFunctionSpace(mesh, 'P', 1)
# Define boundary condition
tol = 1E-14
def clamped_boundary(x, on_boundary):
return on_boundary and x[0] < tol
bc = DirichletBC(V, Constant((0, 0, 0)), clamped_boundary)
# Define strain and stress
def epsilon(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
#return sym(nabla_grad(u))
def sigma(u):
return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)
# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension() # space dimension
v = TestFunction(V)
f = Constant((0, 0, -rho*g))
T = Constant((0, 0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds
# Compute solution
u = Function(V)
solve(a == L, u, bc)
# Plot solution
plot(u, title='Displacement', mode='displacement')
# Plot stress
s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d) # deviatoric stress
von_Mises = sqrt(3./2*inner(s, s))
V = FunctionSpace(mesh, 'P', 1)
von_Mises = project(von_Mises, V)
plot(von_Mises, title='Stress intensity')
# Compute magnitude of displacement
u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, V)
plot(u_magnitude, 'Displacement magnitude')
print('min/max u:',
u_magnitude.vector().array().min(),
u_magnitude.vector().array().max())
# Save solution to file in VTK format
File('elasticity/displacement.pvd') << u
File('elasticity/von_mises.pvd') << von_Mises
File('elasticity/magnitude.pvd') << u_magnitude
# Hold plot
interactive()
Traceback (most recent call last):
File "fenics_ft06_elasticity.py", line 48, in <module>
a = inner(sigma(u), epsilon(v))*dx
File "fenics_ft06_elasticity.py", line 40, in sigma
return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)
NameError: name 'nabla_div' is not defined
I found that replacing 'nabla_div(u)' with just 'div(u)' solved that error. However, it did lead straight to the next error:
Traceback (most recent call last):
File "fenics_ft06_elasticity.py", line 56, in <module>
plot(u, title='Displacement', mode='displacement')
File "/home/ron/miniconda3/envs/fenicsproject/lib/python3.7/site-packages/dolfin/common/plotting.py", line 438, in plot
return _plot_matplotlib(object, mesh, kwargs)
File "/home/ron/miniconda3/envs/fenicsproject/lib/python3.7/site-packages/dolfin/common/plotting.py", line 282, in _plot_matplotlib
ax.set_aspect('equal')
File "/home/ron/miniconda3/envs/fenicsproject/lib/python3.7/site-packages/matplotlib/axes/_base.py", line 1281, in set_aspect
'It is not currently possible to manually set the aspect '
NotImplementedError: It is not currently possible to manually set the aspect on 3D axes