Based on the tutorial Hyperelasticity — FEniCSx tutorial (jorgensd.github.io), I am trying to reproduce the results with the compressible Ogden model instead of neo–Hookean model. However, principle stretches \lambda_1, \lambda_2, and \lambda_3 (eigenvalues of the right Cauchy–Green tensor C=F^TF) are used in the Ogden model.
(1) So my first question is do we have any functions to directly represent the eigenvalues of a tensor in the variational form in the latest version of FEniCS-X?
(2) In addition, I am trying to transfer the function given in the post Copy of ufl.algebra.product to the FEniCSX version to evaluate the eigenvalues indirectly. And I would like to share my code for future readers:
import numpy as np
import dolfinx
from dolfinx.cpp.mesh import CellType
from mpi4py import MPI
from petsc4py import PETSc
import ufl
mesh = dolfinx.BoxMesh(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([10, 5, 5])], \
[2, 1, 1], CellType.hexahedron)
V = dolfinx.VectorFunctionSpace(mesh, ("CG", 2))
M = 1 # number of terms
mu = [1e3, 1e3, 1e3]
alpha = [1, 2, 3]
beta = [1, 2, 3]
def left(x):
return np.isclose(x[0], 0)
def right(x):
return np.isclose(x[0], 10)
left_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1, left)
right_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1, right)
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full(len(left_facets), 1, dtype=np.int32), np.full(len(right_facets), 2, dtype=np.int32)])
sorted_facets = np.argsort(marked_facets)
facet_tag = dolfinx.MeshTags(mesh, mesh.topology.dim-1, marked_facets[sorted_facets], marked_values[sorted_facets])
u_bc = dolfinx.Function(V)
with u_bc.vector.localForm() as loc:
loc.set(0)
left_dofs = dolfinx.fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.indices[facet_tag.values==1])
bcs = [dolfinx.DirichletBC(u_bc, left_dofs)]
B = dolfinx.Constant(mesh, (0, 0, 0)) # body force
T = dolfinx.Constant(mesh, (0, 0, 0)) # traction
u = dolfinx.Function(V) # test function (displacment field)
v = ufl.TestFunction(V) # trial function
d = len(u) # spatial dimension
I = ufl.variable(ufl.Identity(d)) # identity tensor
F = ufl.variable(I + ufl.grad(u)) # deformation gradient
C = ufl.variable(F.T * F) # right Cauchy-Green tensor
Ic = ufl.variable(ufl.tr(C)) # first invariant
J = ufl.variable(ufl.det(F)) # third invariant
def EigenValues3D(T, tol=1e-12):
# determine perturbation from tolerance
eps = 2*tol
# invariants
I1 = ufl.tr(T) # first invariant
I2 = (I1**2 - ufl.inner(T,T)) / 2 # second invariant
I3 = ufl.det(T) # third invariant
p = I1**2 - 3*I2 # preliminary value for p
p = ufl.conditional(ufl.lt(p,tol), abs(p)+eps, p)
# add numerical perturbation to p to ensure the positiveness if close to zero
q = 27/2*I3 + I1**3 - 9/2*I1*I2 # preliminary value for q
q = ufl.conditional(ufl.lt(abs(q),tol), q+ufl.sign(q)*eps, q)
# add numerical perturbation with sign to q if close to zero
# determine angle phi for calculation of roots
phiNom2 = 27 * ( 1/4*I2**2*(p-I2) + I3*(27/4*I3-q) )
# preliminary value for squared nominator of expression for angle phi
phiNom2 = ufl.conditional(ufl.lt(phiNom2,tol), abs(phiNom2)+eps, phiNom2)
# add numerical perturbation to ensure non-zero nominator expression for angle phi
phi = 1/3*ufl.atan_2(ufl.sqrt(phiNom2),q) # calculate angle phi
# calculate polynomial roots
lambda1 = 1/3 * ( ufl.sqrt(p)*2*ufl.cos(phi) + I1 )
lambda2 = 1/3 * (-ufl.sqrt(p)*( ufl.cos(phi) + ufl.sqrt(3)*ufl.sin(phi) ) + I1)
lambda3 = 1/3 * (-ufl.sqrt(p)*( ufl.cos(phi) - ufl.sqrt(3)*ufl.sin(phi) ) + I1)
# return polynomial roots (eigenvalues)
return lambda1, lambda2, lambda3
lambda1, lambda2, lambda3 = EigenValues3D(C) # eigenvalues
# strain energy density of compressible Ogden model
psi = sum([mu[a]/alpha[a]*(lambda1**alpha[a] + lambda2**alpha[a] + lambda3**alpha[a] - 3) for a in range(M)]) \
+ sum([mu[a]/alpha[a]/beta[a]*( J**(-alpha[a]*beta[a]) - 1 ) for a in range(M)])
P = ufl.diff(psi, F) # first Piola-Kirchhoff stress
metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", metadata=metadata)
R = ufl.inner(ufl.grad(v), P)*dx - ufl.inner(v, B)*dx - ufl.inner(v, T)*ds(2)
problem = dolfinx.fem.NonlinearProblem(R, u, bcs)
solver = dolfinx.NewtonSolver(MPI.COMM_WORLD, problem)
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental" # 'residual'
tval0 = -1.5
for n in range(1, 5):
T.value[2] = n * tval0
num_its, converged = solver.solve(u)
assert(converged)
u.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
print(f"Time step {n}, Number of iterations {num_its}, Load {T.value}")
print(u.compute_point_values().round(3))
This code can give us the results without any bugs when I run it on the work station.
Time step 1, Number of iterations 6, Load [ 0. 0. -1.5]
Time step 2, Number of iterations 4, Load [ 0. 0. -3.]
Time step 3, Number of iterations 4, Load [ 0. 0. -4.5]
Time step 4, Number of iterations 4, Load [ 0. 0. -6.]
[[-0.793 0.31 0.249]
[-1.443 0.322 0.123]
[-0.793 -0.31 0.249]
[-1.443 -0.322 0.123]
[-0.698 0.317 -0.379]
[-1.312 0.323 -0.523]
[-0.698 -0.317 -0.379]
[-1.312 -0.323 -0.523]
[ 0. -0. 0. ]
[ 0. 0. -0. ]
[ 0. -0. 0. ]
[ 0. 0. 0. ]]
However, when I run it on my laptop, there are some timeout errors.
Traceback (most recent call last):
File "/usr/local/lib/python3.8/dist-packages/ffcx/codegeneration/jit.py", line 62, in get_cached_module
open(c_filename, "x")
FileExistsError: [Errno 17] File exists: '/root/.cache/fenics/libffcx_forms_7afeadb1fbb16ced3e564cb0b1649784e87bd431.c'
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "/shared/FEM_3D_Large_Deformation/Test.py", line 91, in <module>
problem = dolfinx.fem.NonlinearProblem(R, u, bcs)
File "/usr/local/dolfinx-real/lib/python3.8/dist-packages/dolfinx/fem/problem.py", line 178, in __init__
self._a = fem.form.Form(J, form_compiler_parameters=form_compiler_parameters,
File "/usr/local/dolfinx-real/lib/python3.8/dist-packages/dolfinx/fem/form.py", line 42, in __init__
self._ufc_form, module, self._code = jit.ffcx_jit(
File "/usr/local/dolfinx-real/lib/python3.8/dist-packages/dolfinx/jit.py", line 61, in mpi_jit
return local_jit(*args, **kwargs)
File "/usr/local/dolfinx-real/lib/python3.8/dist-packages/dolfinx/jit.py", line 215, in ffcx_jit
r = ffcx.codegeneration.jit.compile_forms([ufl_object], parameters=p_ffcx, **p_jit)
File "/usr/local/lib/python3.8/dist-packages/ffcx/codegeneration/jit.py", line 152, in compile_forms
obj, mod = get_cached_module(module_name, form_names, cache_dir, timeout)
File "/usr/local/lib/python3.8/dist-packages/ffcx/codegeneration/jit.py", line 84, in get_cached_module
raise TimeoutError(f"""JIT compilation timed out, probably due to a failed previous compile.
TimeoutError: JIT compilation timed out, probably due to a failed previous compile.
Try cleaning cache (e.g. remove /root/.cache/fenics/libffcx_forms_7afeadb1fbb16ced3e564cb0b1649784e87bd431.c) or increase timeout parameter.
So my second question is how can we increase the timeout parameter in JIT compilation?