I have seen that others have tried to do this. So I thought that I could share it
# Copyright 2023 edgar (fenicsproject.discourse.group)
# License: GNU General Public License version 3
# https://www.gnu.org/licenses/gpl-3.0.en.html
#
# Partially based and adapted from https://jsdokken.com/dolfinx-tutorial/
# which holds a Creative Commons Attribution 4.0 International License
# http://creativecommons.org/licenses/by/4.0/
import numpy as np
import pyvista
from mpi4py import MPI
from dolfinx import mesh
from dolfinx import fem
import ufl
from dolfinx.fem import FunctionSpace
from petsc4py.PETSc import ScalarType
from dolfinx import plot
from dolfinx import io
import dolfinx.fem.petsc as petsc
def clamped_boundary(x):
return np.isclose(x[0], 0)
def epsilon(u):
# ufl.sym: Equivalent to
# 0.5 * (∇(u) + (∇(u))^T)
# 0.5 * (ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
return ufl.sym(ufl.grad(u))
def sigma(u, lambda_, mu):
return (lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u))
+ 2 * mu * epsilon(u))
# 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
box_vertex1 = np.array([0, 0, 0])
box_vertex2 = np.array([L, W, W])
ncells = [20, 6, 6]
domain = mesh.create_box(
MPI.COMM_WORLD, [box_vertex1, box_vertex2],
ncells, cell_type=mesh.CellType.hexahedron)
fe_type = "Lagrange"
fe_degree = 1
gdim = domain.geometry.dim
V = FunctionSpace(domain, (fe_type, fe_degree, (gdim,)))
fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(
domain, fdim, clamped_boundary)
u_D = np.array([0, 0, 0], dtype=ScalarType)
boundary_dof = fem.locate_dofs_topological(
V, fdim, boundary_facets)
bc = fem.dirichletbc(u_D, boundary_dof, V)
T = fem.Constant(domain, ScalarType((0, 0, 0)))
ds = ufl.Measure("ds", domain=domain)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, ScalarType((0, 0, -rho * g)))
a = ufl.inner(sigma(u, lambda_, mu), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds
problem = petsc.LinearProblem(
a, L, bcs=[bc],
petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
topology, cell_types, geometry = plot.vtk_mesh(V)
vtu = pyvista.UnstructuredGrid(
topology, cell_types, geometry)
vtu["u"] = uh.x.array.reshape((geometry.shape[0], 3))
warped = vtu.warp_by_vector("u", factor=1.5)
pvpl = pyvista.Plotter()
undef_actor = pvpl.add_mesh(vtu, style="wireframe", color="k")
warped_actor = pvpl.add_mesh(warped, show_edges=True, scalars="u")
# https://tutorial.pyvista.org/tutorial/04_filters/solutions/e_glyphs.html
arrows = vtu.glyph(
scale=True, orient=True, tolerance=0.05, factor=0.25)
pvpl.add_mesh(arrows, show_scalar_bar=False)
pvpl.show_axes()
pvpl.save_graphic("Media/dolfinx_tuto_linear_elasticity.svg")
pvpl.show()
sigma_dev = sigma(uh, lambda_, mu) - (
1 / 3 * ufl.tr(sigma(uh, lambda_, mu))
* ufl.Identity(len(uh)))
sigma_vm = ufl.sqrt(3 / 2 * ufl.inner(sigma_dev, sigma_dev))
V_vm = fem.FunctionSpace(
domain, ("Discontinuous Lagrange", 0))
sigma_vm_expr = fem.Expression(
sigma_vm, V_vm.element.interpolation_points())
sigma_vm_fun = fem.Function(V_vm)
sigma_vm_fun.interpolate(sigma_vm_expr)
warped.cell_data["sigma_vm"] = sigma_vm_fun.vector.array
warped.set_active_scalars("sigma_vm")
pvpl = pyvista.Plotter()
pvpl.add_actor(undef_actor)
pvpl.add_mesh(warped, scalars="sigma_vm")
pvpl.show_axes()
pvpl.save_graphic("Media/dolfinx_tuto_linear_elasticity_sigma_fun_vm.svg")
pvpl.show()
tens2_dim = (domain.geometry.dim,) * 2
interp_ord = 1
fe_type = "Lagrange"
CG_lin_2ord = fem.FunctionSpace(domain, (fe_type, interp_ord, tens2_dim))
sigma_expr = fem.Expression(
sigma(uh, lambda_, mu), CG_lin_2ord.element.interpolation_points())
sigma_fun = fem.Function(CG_lin_2ord)
sigma_fun.interpolate(sigma_expr)
vtu["sigma"] = sigma_fun.x.array.reshape((geometry.shape[0], 3, 3))
vtu.set_active_tensors("sigma")
# warped["sigma"] = vtu.warp_by_vector("u", factor=1.5)
warped["sigma"] = vtu["sigma"]
# https://docs.pyvista.org/version/stable/examples/02-plot/multi-window.html
pyvista.global_theme.multi_rendering_splitting_position = 0.60
pvpl = pyvista.Plotter(shape="1|2")
pvpl.subplot(0)
pvpl.add_actor(undef_actor)
pvpl.add_mesh(warped, scalars="sigma",
scalar_bar_args={"title": "stress norm"})
pvpl.show_axes()
pvpl.subplot(1)
pvpl.add_mesh(warped, scalars="sigma", component=0,
scalar_bar_args={"title": "stress xx"})
pvpl.subplot(2)
pvpl.add_mesh(warped, scalars="sigma", component=8,
scalar_bar_args={"title": "stress zz"})
pvpl.save_graphic("Media/dolfinx_tuto_linear_elasticity_sigma_fun.svg")
pvpl.show()
Thanks for FEniCS! Thanks to @dokken, @nate, @jackhale, @dparsons parsons for paying attention, sharing your time and helping everyone here.
-
FEniCSx software
dolfinx: 0.7.0.dev0_r27554.cfeffe0-1
basix: 0.7.0.dev0_r945.1117a8d-1
ufl: 2023.2.0.dev0_r3562.77ae57c-1
ffcx: 0.7.0.dev0_r7077.1d27238-1 -
Dependencies
python: 3.11.5-2
python-numpy: 1.26.0-1
petsc: 3.19.5.1171.g37df9106526-1
hdf5-openmpi: 1.14.2-1
boost: 1.83.0-2
adios2: 2.8.3-5
scotch: 7.0.4-1
pybind11: 2.11.1-1
python-build: 1.0.1-1
python-cffi: 1.15.1-4
python-cppimport: 22.08.02.r6.g0849d17-1
python-installer: 0.7.0-3
python-mpi4py: 3.1.4-3
python-pytest: 7.4.2-1
python-scikit-build: 0.17.6-1
python-setuptools: 1:68.0.0-1
python-wheel: 0.40.0-3
xtensor: 0.24.0-2
xtensor-blas: 0.20.0-1 -
Operating system
Arch Linux: 6.5.4-arch2-1
