How to plot 3d beam after deformation, not colormap

Uncategorized
1

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
2

you can try with:

# plot solution
from vedo.dolfin import plot
plot(u, mode="displacements with arrows", text='displaced mesh', interactive=False)
plot(SED, text='strain energy density', new=True, size=(900,400), pos=(1100,0), zoom=2)

Screenshot from 2021-04-05 12-50-34

3

Hi marcomusy, Thank you for the solution. It is beautiful and exactly what I wanted. Unfortunately I use Fenics Docker, which does not have vedo. But I sort of figured out in parameter, to use wrap by vector filter.

4

i’m not very familiar with docker, but you can also generate an offscreen image with

from vedo.dolfin import plot
cam = dict(pos=(2.04, -0.751, 0.336),
           focalPoint=(0.563, 0.169, -0.0380),
           viewup=(-0.162, 0.137, 0.977),
           distance=1.78,
           clippingRange=(0.676, 3.15))

plt = plot(u, mode="displacements with arrows", text='displaced mesh', camera=cam, offscreen=True)
plt.screenshot('screenshot.png')

…if it helps

5

Thank you so much marcomusy. Right now it seems the docker image does not have vedo package; and I am not skillful enough to build a new docker image. I have written down your solution in case I will run another FEniCS system on Anaconda.

6

You can simply install vedo in the docker container by calling pip3 install vedo

7

Thank you dokken. Would the vedo stay with the image, or I need to install it every time when a container is created?

8

This depends on your usage of docker.
If you always run docker run -it ... you are creating a new container each time. However, let’s say you run docker run -ti .... —name=dolfin_container ....
and install Vedo in this container, you can start this container (with vedo) later by running docker container start -i dolfin_container.

9

That makes sense. Thank you @dokken and @marcomusy .