Hello everyone, I have studied this linear elasticity tutorialhttps://jsdokken.com/dolfinx-tutorial/chapter2/linearelasticity_code.html, but it only calculates the von Mises stress. I would like to know how to solve for the entire stress field and save the results. Thank you!
For example, to solve for the stress field in the tutorial.
You mean outputting the stress tensor sigma? You can do that in a similar fashion to von Mises stresses. Just create a tensor function space, and use an expression with sigma as input to interpolate.
In fact, I have made attempts, but the TensorFunctionSpace is discontinued from the dolfinx.fem module. And I used the code below, but an error occurred. I don’t know what the cppV parameter is, which prevents me from solving for the stress field.
gdim = mesh.geometry.dim
Function_space_for_sigma = FunctionSpace(mesh, ("DG", 1, (gdim, gdim)))
expr= Expression(sigma(u_n,p_n), Function_space_for_sigma.element.interpolation_points())
sigma_values = Function(Function_space_for_sigma)
sigma_values.interpolate(expr)
Function_space_for_sigma = FunctionSpace(mesh, ("DG", 1, (gdim, gdim)))
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
TypeError: FunctionSpace.__init__() missing 1 required positional argument: 'cppV'
Could you provide me with an example based on the tutorial in the link above? Thank you very much!
Thank you for your reply. I followed your advice. In fact, I am solving the stress field for a thermomechanical problem, so I created the following minimal example based on the linear elasticity tutorial you provided and solved for the stress field：
import ufl
import numpy as np
from mpi4py import MPI
from dolfinx.fem.petsc import LinearProblem
from dolfinx import mesh, fem, io, default_scalar_type
from dolfinx.fem import (functionspace, Function, Expression)
import basix
L = 1
W = 0.2
E = 2e5
nu = 0.3
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
alpha = 1.2e-5
Delta_T = 100
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, W])],[20, 6, 6], cell_type=mesh.CellType.hexahedron)
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))
def clamped_boundary(x):
return np.isclose(x[0], 0)
fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)
u_D = np.array([0, 0, 0], dtype=default_scalar_type)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)
def eps(v):
return 0.5 * (ufl.grad(v) + ufl.grad(v).T)
def sigma(v, dT):
return (lmbda * ufl.tr(eps(v)) - alpha * (3 * lmbda + 2 * mu) * dT) * ufl.Identity(3) + 2.0 * mu * eps(v)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
Wint = ufl.inner(sigma(u, Delta_T), eps(v)) * ufl.dx
a = ufl.lhs(Wint)
L = ufl.rhs(Wint)
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
# stress field
cell = domain.basix_cell()
tensor_el = basix.ufl.element("Lagrange", cell, 2, shape=(domain.geometry.dim, domain.geometry.dim))
W = functionspace(domain, tensor_el)
expr = Expression(sigma(uh, Delta_T), W.element.interpolation_points())
sigma_values = Function(W)
sigma_values.interpolate(expr)
with io.XDMFFile(domain.comm, "results/stress.xdmf", "w") as file:
file.write_mesh(domain)
file.write_function(sigma_values)
However, an error occurred. I apologize for my lack of expertise and for being unable to derive the solution from your suggestions. Could you please provide an example of solving the stress field for a thermomechanical problem based on my code? Thank you.
Traceback (most recent call last):
File "/home/dyfluid/Desktop/mises3D/f_solver.py", line 50, in <module>
file.write_function(sigma_values)
File "/home/dyfluid/anaconda3/envs/fenicsx-env/lib/python3.11/site-packages/dolfinx/io/utils.py", line 252, in write_function
super().write_function(getattr(u, "_cpp_object", u), t, mesh_xpath)
RuntimeError: Degree of output Function must be same as mesh degree. Maybe the Function needs to be interpolated?
I think this is simply an error regarding the functional space where the stress lives. Given a piecewise linear displacement on each element, the stress field which in your case is proportional to the deformation (i.e., the derivative of the displacement field) is naturally discontinuous by element, so that you should write:
tensor_el = basix.ufl.element("DG", cell, 0, shape=(domain.geometry.dim, domain.geometry.dim))
This is because you are using XDMFFile, which is limited in what function spaces it can represent (as shown by the error message).
Secondly, as @BouteillerP states, the stresses are naturally discontinuous, and the stresses should in your case be piecewise constant: