Hello,
I am using dolfinx on python to generate displacement fields on a mesh of a cube generated with the gmsh GUI.
I am following the Hyperelasticity tutorial as well as this tutorial with some explanations but deprecated functions.
I successfully implemented the neo-Hookean model with a body load corresponding to the gravity and a traction force applied on the upper side of the cube. I applied Dirichlet boundary conditions on two lateral sides. (cf downloaded picture).
As you can see, the plot corresponds to the deformation field of the initial mesh (red glyphs). My question is: how can I recover the coordinates of these glyphs (vectors on each node) and their value/direction in order to extract this displacement field ?
I provide you a bit of code below :
from dolfinx.plot import create_vtk_mesh
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, assemble_scalar,
form, locate_dofs_geometrical, locate_dofs_topological)
import pyvista
from petsc4py import PETSc
import ufl
import numpy as np
bottom_cells, side_1, side_2, side_cells2, top_cells = ft.find(9), ft.find(13), ft.find(14), ft.find(12), ft.find(8)
V = dolfinx.fem.VectorFunctionSpace(mesh, ("CG", 2))
# side_dofs1 = locate_dofs_topological(V, ft.dim, side_cells1)
bot_dofs = locate_dofs_topological(V, ft.dim, bottom_cells)
top_dofs = locate_dofs_topological(V, ft.dim, top_cells)
u_bc = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
bcs = [dirichletbc(u_bc, bot_dofs, V), dirichletbc(u_bc, top_dofs, V)]
u = Function(V)
v = ufl.TestFunction(V)
# Identity tensor
d = len(u) # Spatial dimension
I = ufl.variable(ufl.Identity(d))
# Deformation gradient
F = ufl.variable(ufl.grad(u) + I)
# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)
# Invariants of deformation tensors
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))
# Elasticity parameters
# E = Young modulus
# nu is to be adjusted as exercise to observe locking.
E, nu = 5000, 0.495
mu = Constant(mesh, E/(2.0*(1.0 + nu)))
lmbda = Constant(mesh, E*nu/((1.0 + nu)*(1 - 2.0*nu)))
# Body load in undeformed configuration
g, rho = -9.81, 1
# b = Constant(mesh, [0.0, 0.0, rho*g])
B = Constant(mesh, PETSc.ScalarType((0, 0, rho*g)))
T = Constant(mesh, PETSc.ScalarType((100, 0, 0)))
# Stored strain energy density (compressible neo-Hookean model)
psi = ufl.variable((mu / 2.0) * (Ic - 2) - mu * ufl.ln(J) + (lmbda / 2.0) * (ufl.ln(J))**2)
# dx = ufl.Measure("dx")
# Pi = psi*dx - ufl.inner(b, u)*dx
# First Piola-Kirchhoff Stress
P = ufl.diff(psi, F)
metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=mesh, subdomain_data=ft, metadata=metadata)
dx = ufl.Measure("dx", domain=mesh, subdomain_data = ct, metadata=metadata)
# Define form F (we want to find u such that F(u) = 0)
F = ufl.inner(ufl.grad(v), P)*dx(10) - ufl.inner(v, T)*ds(14) - ufl.inner(v, B)*dx(10)
Pi = psi*dx(10) - ufl.inner(u, T)*ds(14)
# - ufl.inner(u, B)*dx(10)
F_res = ufl.derivative(Pi, u, v)
Jac = ufl.derivative(F_res, u)
problem1 = dolfinx.fem.petsc.NonlinearProblem(F, u, bcs)
from dolfinx import nls
solver1 = nls.petsc.NewtonSolver(mesh.comm, problem1)
# Set Newton solver options
solver1.atol = 1e-8
solver1.rtol = 1e-8
solver1.convergence_criterion = "incremental"
num_its, converged = solver1.solve(u)```
from xml import dom
import pyvista
from dolfinx import plot
pyvista.set_jupyter_backend("pythreejs")
domain = mesh
grid = pyvista.UnstructuredGrid(*plot.create_vtk_mesh(domain, domain.topology.dim))
p = pyvista.Plotter()
# Create grid defined by the function space for visualization of the function
topology, cells, geometry = plot.create_vtk_mesh(u.function_space)
function_grid = pyvista.UnstructuredGrid(topology, cells, geometry)
var_name = f"u"
values = np.zeros((geometry.shape[0], 3))
values[:, :len(u)] = u.x.array.reshape(geometry.shape[0], len(u))
function_grid[var_name] = values
# function_grid.set_active_vectors(var_name)
# Warp mesh by deformation
warped = function_grid.warp_by_vector(var_name, factor=5)
glyphs = function_grid.glyph(var_name, factor=5)
# warped = grid.warp_by_vector(var_name, factor=1.5)
actor_1 = p.add_mesh(warped)
actor_0 = p.add_mesh(function_grid, style="wireframe", color="k")
p.add_mesh(glyphs, color="r")
# Add mesh to plotter and visualize
# actor = p.add_mesh(warped, show_edges=True)
pyvista.OFF_SCREEN = False
if not pyvista.OFF_SCREEN:
p.show()
else:
pyvista.start_xvfb()
figure_as_array = p.screenshot(f"diffusion_{t:.2f}.png")
# Clear plotter for next plot
p.remove_actor(actor_1)
