Hi, I just began learning FEA/FEniCS. I followed the tutorial of linear elasticity of clamped 3d beam bending by its gravity. The codes run OK except that plot(u, mode=‘displacement’ gives an error of"NotImplementedError: It is not currently possible to manually set the aspect on 3D axes". Google shows it is because the latest matplotlib has problems of setting equal aspect for 3d. So I commented that line, and exported u to a .pvd file. In paraview it shows the original beam with some colormap. How can I plot the actual bended beam like what is shown in the tutorial, instead of a colormaping? I am not sure if this is a FEniCS or paraview question. Thank you!

The link to the tutorial:

https://fenicsproject.org/pub/tutorial/html/._ftut1008.html

and my codes are in the following (almost identical to the tutorial).

```
from fenics import *
from ufl import nabla_div, nabla_grad
#import matplotlib.pyplot as plt
# 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
#Boxmesh(x0,y0,z0, x,y,z,nx,ny,nz) create tetrahegonal mesh between two points with number of cells
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
# DirichletBC(FunctionSpace, value, boundary)
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)
def sigma(u):
return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)
# define vairational problem
u = TrialFunction(V)
d = u.geometric_dimension() #space dimension
v = TestFunction(V)
f = Constant((0, 0, -rho*g)) # gravity
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') #Implemtn name error
# plot stress
s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d) # dviatoric 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 magnitue of displacement
u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, V)
print('min/max u: ', u_magnitude.vector().get_local().min(), u_magnitude.vector().get_local().max())
# compute strain energy density
strain_energy_density = 0.5*inner(sigma(u), epsilon(u))
SED= project(strain_energy_density, V)
File('ft06_elasticity_u.pvd') << u
File('ft06_elasticity_vonMises.pvd') << von_Mises
File('ft06_elasticity_magnitude.pvd') << u_magnitude
```