Hi all,
I’m trying to plot Von Mises stress for a simply supported Reissner-Mindlin plate.
I’m following exemple code from fenics-shells:
https://fenics-shells.readthedocs.io/en/latest/demo/reissner-mindlin-simply-supported/demo_reissner-mindlin-simply-supported.py.html
Thanks in advance!
See for instance: How to restrict solutions or stresses to subdomains? - #5 by steve
which project the Von mises stresses into an appropriate function space, and then plots it with matplotlib.
For 3D geometries, I suggest saving the stresses to XDMF file with write_checkpoint and use Paraview to visualize them.
Thank you Dokken!
I tried to implement it in my code:
import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np
from dolfin import *
from fenics_shells import *
from ufl import nabla_div
mesh = RectangleMesh(Point(0.0, 0.0), Point(1000.0, 2000.0), 64, 64)
element = MixedElement([VectorElement("Lagrange", triangle, 2),
FiniteElement("Lagrange", triangle, 1),
FiniteElement("N1curl", triangle, 1),
FiniteElement("N1curl", triangle, 1)])
U = ProjectedFunctionSpace(mesh, element, num_projected_subspaces=2)
U_F = U.full_space
u_ = Function(U_F)
theta_, w_, R_gamma_, p_ = split(u_)
u = TrialFunction(U_F)
u_t = TestFunction(U_F)
E = Constant(70000.0)
nu = Constant(0.22)
kappa = Constant(5.0/6.0)
t = Constant(10)
k = sym(grad(theta_))
D = (E*t**3)/(24.0*(1.0 - nu**2))
psi_M = D*((1.0 - nu)*tr(k*k) + nu*(tr(k))**2)
psi_T = ((E*kappa*t)/(4.0*(1.0 + nu)))*inner(R_gamma_, R_gamma_)
f = Constant(0.000001)
W_ext = inner(f*t**3, w_)*dx
gamma = grad(w_) - theta_
# We instead use the :py:mod:`fenics_shells` provided :py:func:`inner_e` function::
L_R = inner_e(gamma - R_gamma_, p_)
L = psi_M*dx + psi_T*dx + L_R - W_ext
F = derivative(L, u_, u_t)
J = derivative(F, u_, u)
sigma = E/(1.0-nu**2)*(u_.dx(0).dx(0) + nu*u_.dx(1).dx(1))*t/2
xdmf_file = XDMFFile("sigma")
sigma_out = project(sigma, U_F)
sigma_out.rename("sigma", "stress")
xdmf_file.write(sigma_out, 0.0)
A, b = assemble(U, J, -F)
# and apply simply-supported boundary conditions::
def all_boundary(x, on_boundary):
return on_boundary
def left(x, on_boundary):
return on_boundary and near(x[0], 0.0)
def right(x, on_boundary):
return on_boundary and near(x[0], 1000.0)
def bottom(x, on_boundary):
return on_boundary and near(x[1], 0.0)
def top(x, on_boundary):
return on_boundary and near(x[1], 2000.0)
# Simply supported boundary conditions.
bcs = [DirichletBC(U, Constant((0.0,0.0,0.0)), all_boundary),
DirichletBC(U.sub(0).sub(0), Constant(0.0), top),
DirichletBC(U.sub(0).sub(0), Constant(0.0), bottom),
DirichletBC(U.sub(0).sub(1), Constant(0.0), left),
DirichletBC(U.sub(0).sub(1), Constant(0.0), right)]
for bc in bcs:
bc.apply(A, b)
# and solve the linear system of equations before writing out the results to
# files in ``output/``::
u_p_ = Function(U)
solver = PETScLUSolver("mumps")
solver.solve(A, u_p_.vector(), b)
reconstruct_full_space(u_, u_p_, J, -F)
save_dir = "output/"
theta, w, R_gamma, p = u_.split()
fields = {"theta": theta, "w": w, "R_gamma": R_gamma, "p": p}
for name, field in fields.items():
field.rename(name, name)
field_file = XDMFFile("%s/%s.xdmf" % (save_dir, name))
field_file.write(field)
plt.xlabel("Tizio%s"%name)
plt.ylabel("Caio")
c = plot(field)
cb = plt.colorbar(c)
plt.savefig("output/1-%s.png"%name)
cb.remove()
print("{customTag}%s{/customTag}" % name)
plt.figure()
c = plot(sigma)
cb = plt.colorbar(c)
plt.savefig("sigma.png")
cb.remove()
# We check the result against an analytical solution calculated using a
# series expansion::
from fenics_shells.analytical.simply_supported import Displacement
w_e = Displacement(degree=3)
w_e.t = t.values()
w_e.E = E.values()
w_e.p = f.values()*t.values()**3
w_e.nu = nu.values()
print("Numerical out-of-plane displacement at centre: %.4e" % w((500, 1000)))
print("Analytical out-of-plane displacement at centre: %.4e" % w_e((500, 1000)))
# Unit testing
# ============
#
# ::
def test_close():
import numpy as np
assert(np.isclose(w((500, 500)), w_e((500, 500)), atol=1E-3, rtol=1E-3))
But I have this error:
Traceback (most recent call last):
File "demo_reissner-mindlin-simply-supported.py", line 69, in <module>
sigma_out = project(sigma, U_F)
File "/home/fenics/shared/fenics_shells/fem/solving.py", line 26, in project
u_V_F = Function(V)
NameError: name 'Function' is not defined
Thank you for help
I cannot find this line in the code you have supplied above, so I cannot really help you.
Anyhow, I think ProjectedFunctionSpace is part of the fenics-shells code, and I don’t know much about that implementation. Maybe @jackhale can enlighten us.
I’m sorry, this is the error:
Traceback (most recent call last):
File "demo_reissner-mindlin-simply-supported.py", line 69, in <module>
sigma_out = project(sigma, U_F)
File "/home/fenics/shared/fenics_shells/fem/solving.py", line 26, in project
u_V_F = Function(V)
NameError: name 'Function' is not defined
I’m one of the developers of FEniCS-Shells.
This is a bug in a piece of code not touched under our integrated tests, I will fix it.
In the meantime, explicitly use import dolfin and dolfin.project.
Looking at bigger picture, there are a number of issues I can see in your code:
This code only solves the out-of-plane Reissner-Mindlin problem. The in-plane (membrane) and out-of-plane (bending) elasticity problems decouple for an isotropic Reissner-Mindlin plate. I do not know the good definition of the von-Mises stress for a Reissner-Mindlin plate but it may involve both membrane and bending stresses. Worth checking?
The von-Mises stress is a quantity you get from post-processing the solution, i.e. you can only compute it after solving for the unknown functions u_p (and u_).
The bug in FEniCS-Shells is now fixed (specifically, the method project has been removed).