Hi, I’m solving a problem in elasticity and I want to plot the displacement, however when I plot it in the undeformed configuration I obtain a bad plot.
I’m just using the command plot:
plot(u)
plt.show()
Is it possible to obtain a good contour plot for the undeformed configuration?
If not how can I plot my solution using the contour plot in matplotlib? (I was thinking about using griddata to interpolate, but I don’t know how to do).
Thanks
Without a reproducible example, it is unlikely that you will get any help.
Could you first save it to file and visualize it in Paraview to ensure that it is just the plotting that is bad, and not the solution itself?
Hi, thank you for your reply.
I’m just reproducing this example
2D Hyper-elasticity — FEniCS solid tutorial 2019.0.2 documentation,
so the answer must be correct. Indeed the plot in the deformed configuration is equal and it is good.
I have removed ‘displacement’ from fe.plot amd I get this plot not good. I would like to have one plot as the one in the deformed configuration but maintaining the initial undeformed one. How can I do?
As I don’t have the mesh, I cannot reproduce it.
Could you try visualizing it in Paraview, as it is a better visualization engine than matplotlib.
Note that since you are using the plot command with matplotlib, you can send in various commands to the plot function to tweak the plot.
Thank you, I have created the mesh just by using mshr. If you want I can write here my code.
Do you know what command do I have to use in order to have the right colours like this
You havent provided the code, and you have tried my suggestion check by opening it in Paraview. Please try it as a sanity check.
I am not able to use Paraview, I have always used only the plot tools provided by Fenics or Matplotlib. If you want I can provide my code
Please provide the code then. Bote that mshr has been deprecated for a long time.
from dolfin import *
import fenics as fe
import matplotlib.pyplot as plt
import numpy as np
from mshr import *
#import ufl
fe.set_log_level(13)
# --------------------
# Functions and classes
# --------------------
def bottom(x, on_boundary):
return (on_boundary and fe.near(x[1], 0.0))
# Strain function
def epsilon(u):
return 0.5*(fe.grad(u) + fe.grad(u).T)
# Stress function
def sigma(u):
return lmbda*fe.div(u)*fe.Identity(2) + 2*mu*epsilon(u)
# --------------------
# Parameters
# --------------------
# Lame's constants
#lmbda = 1.25
#mu = 20
E = 20.0e6
mu = 0.5*E
rho_0 = 200.0
l_x, l_y = 5.0, 5.0 # Domain dimensions
n_x, n_y = 500, 500 # Number of elements
# Load
g_int = 1.0e5
b_int = -10.0
# --------------------
# Geometry
# --------------------
#mesh_quad = fe.UnitSquareMesh.create(n_x, n_y, fe.CellType.Type.quadrilateral)
domain = Rectangle(Point(0., 0.), Point(5, 5)) - Circle(Point(2.5, 2.5), 1.8, 100)
mesh = generate_mesh(domain, 500)
fe.plot(mesh)
plt.show()
# Definition of Neumann condition domain
boundaries = fe.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
top = fe.AutoSubDomain(lambda x: fe.near(x[1], 5.0))
top.mark(boundaries, 1)
ds = fe.ds(subdomain_data=boundaries)
# --------------------
# Function spaces
# --------------------
V = fe.VectorFunctionSpace(mesh, "CG", 1)
u_tr = fe.TrialFunction(V)
u_test = fe.TestFunction(V)
u = fe.Function(V)
g = fe.Constant((0.0, g_int))
b = fe.Constant((0.0, b_int))
N = fe.Constant((0.0, 1.0))
aa, bb, cc, dd, ee = 0.5*mu, 0.0, 0.0, mu, -1.5*mu
# --------------------
# Boundary conditions
# --------------------
bc = fe.DirichletBC(V, fe.Constant((0.0, 0.0)), bottom)
# --------------------
# Weak form
# --------------------
I = fe.Identity(2)
F = I + fe.grad(u) # Deformation gradient
C = F.T*F # Right Cauchy-Green tensor
J = fe.det(F) # Determinant of deformation fradient
#psi = (aa*[fe.tr](http://fe.tr/)(C) + bb*[fe.tr](http://fe.tr/)(ufl.cofac(C)) + cc*J**2 - dd*fe.ln(J))*fe.dx - fe.dot(b, u)*fe.dx + fe.inner(f, u)*ds(1)
#n = fe.dot(ufl.cofac(F), N)
#surface_def = fe.sqrt(fe.inner(n, n))
psi = (aa* fe.inner(F, F) + ee - dd* fe.ln(J))* fe.dx - rho_0* J* fe.dot(b, u)* fe.dx +fe.inner(g,u)* ds(1) #+ surface_def* fe.inner(g, u)* ds(1)
# --------------------
# Solver
# --------------------
Form = fe.derivative(psi, u, u_test)
Jac = fe.derivative(Form, u, u_tr)
problem = fe.NonlinearVariationalProblem(Form, u, bc, Jac)
solver = fe.NonlinearVariationalSolver(problem)
prm = solver.parameters
#prm["newton_solver"]["error_on_convergence"] = False
#fe.solve(Form == 0, u, bc, J=Jac, solver_parameters={"error_on_convergence": False})
solver.solve()
print(np.amax(u.vector()[:]))
# --------------------
# Post-process
# --------------------
fe.plot(u,mode='displacement')
plt.show()
fe.plot(u)
plt.show()
That’s my code, if you run it you can see that the first plot provides the deformed configuration correctly, instead in the second the undeformed is not good.
Do you know what can I change in the second plot command?
Or you know if it is possible to plot u using matplotlib?
I have also avoided to use ufl.cofac because it doesn’t work and I don’t know why
Your code isn’t properly formatted. I tried to fix it, but it is missing several indentation. Please Edit it to ensure that the code can be copy pasted by someone and executed with no modification.
Ok, I modify it immediately
#add code here
from dolfin import *
import fenics as fe
import matplotlib.pyplot as plt
import numpy as np
from mshr import *
#import ufl
fe.set_log_level(13)
# --------------------
# Functions and classes
# --------------------
def bottom(x, on_boundary):
return (on_boundary and fe.near(x[1], 0.0))
# Strain function
def epsilon(u):
return 0.5*(fe.grad(u) + fe.grad(u).T)
# Stress function
def sigma(u):
return lmbda*fe.div(u)*fe.Identity(2) + 2*mu*epsilon(u)
# --------------------
# Parameters
# --------------------
# Lame's constants
#lmbda = 1.25
#mu = 20
E = 20.0e6
mu = 0.5*E
rho_0 = 200.0
l_x, l_y = 5.0, 5.0 # Domain dimensions
n_x, n_y = 500, 500 # Number of elements
# Load
g_int = 1.0e5
b_int = -10.0
# --------------------
# Geometry
# --------------------
#mesh_quad = fe.UnitSquareMesh.create(n_x, n_y, fe.CellType.Type.quadrilateral)
domain = Rectangle(Point(0., 0.), Point(5, 5)) - Circle(Point(2.5, 2.5), 1.8, 100)
mesh = generate_mesh(domain, 500)
fe.plot(mesh)
plt.show()
# Definition of Neumann condition domain
boundaries = fe.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
top = fe.AutoSubDomain(lambda x: fe.near(x[1], 5.0))
top.mark(boundaries, 1)
ds = fe.ds(subdomain_data=boundaries)
# --------------------
# Function spaces
# --------------------
V = fe.VectorFunctionSpace(mesh, "CG", 1)
u_tr = fe.TrialFunction(V)
u_test = fe.TestFunction(V)
u = fe.Function(V)
g = fe.Constant((0.0, g_int))
b = fe.Constant((0.0, b_int))
N = fe.Constant((0.0, 1.0))
aa, bb, cc, dd, ee = 0.5*mu, 0.0, 0.0, mu, -1.5*mu
# --------------------
# Boundary conditions
# --------------------
bc = fe.DirichletBC(V, fe.Constant((0.0, 0.0)), bottom)
# --------------------
# Weak form
# --------------------
I = fe.Identity(2)
F = I + fe.grad(u) # Deformation gradient
C = F.T*F # Right Cauchy-Green tensor
J = fe.det(F) # Determinant of deformation fradient
#psi = (aa*fe.tr(C) + bb*fe.tr(ufl.cofac(C)) + cc*J**2 - dd*fe.ln(J))*fe.dx - fe.dot(b, u)*fe.dx + fe.inner(f, u)*ds(1)
#n = fe.dot(ufl.cofac(F), N)
#surface_def = fe.sqrt(fe.inner(n, n))
psi = (aa*fe.inner(F, F) + ee - dd*fe.ln(J))*fe.dx - rho_0*J*fe.dot(b, u)*fe.dx +fe.inner(g,u)*ds(1) #+ surface_def*fe.inner(g, u)*ds(1)
# --------------------
# Solver
# --------------------
Form = fe.derivative(psi, u, u_test)
Jac = fe.derivative(Form, u, u_tr)
problem = fe.NonlinearVariationalProblem(Form, u, bc, Jac)
solver = fe.NonlinearVariationalSolver(problem)
prm = solver.parameters
#prm["newton_solver"]["error_on_convergence"] = False
#fe.solve(Form == 0, u, bc, J=Jac, solver_parameters={"error_on_convergence": False})
solver.solve()
print(np.amax(u.vector()[:]))
# --------------------
# Post-process
# --------------------
fe.plot(u,mode='displacement')
plt.show()
Thank you Dokken, this should be correct.
I am referring to this
n = fe.dot(ufl.cofac(F), N)
surface_def = fe.sqrt(fe.inner(n, n))
Now it should work correctly
Do you know also why ufl gives an error in the part with # ?
Again, you need to use 3x` syntax for encapsulation, as you can see above the indentations are all wrong.
example
```python
#add code here
def f(x):
return x
```
Secondly it is unclear what part you are talking about with
As multiple things are commented out
Now it should be correct
This seems like a matplotlib issue, not a DOLFIN issue,as here is the same field shown in Paraview, using:
with fe.XDMFFile("test.xdmf") as xdmf:
xdmf.write(u)
plt.figure()
fe.plot(u) # , mode="displacement")
plt.savefig("disp.png")
yields
Thank you, I have solved by deciding to plot the two components of the displacement separately.