Hello,
I am new to FEniCSx and now studying with demo code.
Among those, I met some error by just following the code on the website :
https://jsdokken.com/dolfinx-tutorial/chapter2/hyperelasticity.html
import numpy as np
import ufl
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh, plot
L = 20.0
domain = mesh.create_box(MPI.COMM_WORLD,[[0.0,0.0,0.0], [L, 1, 1]], [20, 5, 5], mesh.CellType.hexahedron)
V = fem.VectorFunctionSpace(domain, ("CG", 2))
#We create two python functions for determining the facets to apply boundary conditions to
def left(x):
return np.isclose(x[0], 0)
def right(x):
return np.isclose(x[0], L)
fdim = domain.topology.dim -1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)
#Next, we create a marker based on these two functions
# Concatenate and sort the arrays based on facet indices. Left facets marked with 1, right facets with two
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])
#We then create a function for supplying the boundary condition on the left side, which is fixed.
u_bc = np.array((0,) * domain.geometry.dim, dtype=PETSc.ScalarType)
#To apply the boundary condition, we identity the dofs located on the facets marked by the MeshTag.
left_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(1))
bcs = [fem.dirichletbc(u_bc, left_dofs, V)]
#Next, we define the body force on the reference configuration (B), and nominal (first Piola-Kirchhoff) traction (T).
B = fem.Constant(domain, PETSc.ScalarType((0, 0, 0)))
T = fem.Constant(domain, PETSc.ScalarType((0, 0, 0)))
#Define the test and solution functions on the space
v = ufl.TestFunction(V)
u = fem.Function(V)
#Define kinematic quantities used in the problem
# Spatial dimension
d = len(u)
# Identity tensor
I = ufl.variable(ufl.Identity(d))
# Deformation gradient
F = ufl.variable(I + ufl.grad(u))
# 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))
#Define the elasticity model via a stored strain energy density function
#, and create the expression for the first Piola-Kirchhoff stress:
# Elasticity parameters
E = PETSc.ScalarType(1.0e4)
nu = PETSc.ScalarType(0.3)
mu = fem.Constant(domain, E/(2*(1 + nu)))
lmbda = fem.Constant(domain, E*nu/((1 + nu)*(1 - 2*nu)))
# Stored strain energy density (compressible neo-Hookean model)
psi = (mu / 2) * (Ic - 3) - mu * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J))**2
# Stress
# Hyper-elasticity
P = ufl.diff(psi, F)
#Comparison to linear elasticity
#To illustrate the difference between linear and hyperelasticity, the following lines can be uncommented to solve the linear elasticity problem.
# P = 2.0 * mu * ufl.sym(ufl.grad(u)) + lmbda * ufl.tr(ufl.sym(ufl.grad(u))) * I
#Define the variational form with traction integral over all facets with value 2. We set the quadrature degree for the integrals to 4.
metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
# Define form F (we want to find u such that F(u) = 0)
F = ufl.inner(ufl.grad(v), P)*dx - ufl.inner(v, B)*dx - ufl.inner(v, T)*ds(2)
#As the varitional form is non-linear and written on residual form, we use the non-linear problem class from DOLFINx to set up required structures to use a Newton solver.
problem = fem.petsc.NonlinearProblem(F, u, bcs)
The error is as follows:
ArityMismatch Traceback (most recent call last)
Cell In[43], line 87
84 F = ufl.inner(ufl.grad(v), P)*dx - ufl.inner(v, B)*dx - ufl.inner(v, T)*ds(2)
85 #As the varitional form is non-linear and written on residual form, we use the non-linear problem class from DOLFINx to set up required structures to use a Newton solver.
---> 87 problem = fem.petsc.NonlinearProblem(F, u, bcs)
File /usr/local/lib/python3.10/dist-packages/ufl/algorithms/check_arities.py:177, in check_form_arity(form, arguments, complex_mode)
175 def check_form_arity(form, arguments, complex_mode=False):
176 for itg in form.integrals():
--> 177 check_integrand_arity(itg.integrand(), arguments, complex_mode)
...
...
File /usr/local/lib/python3.10/dist-packages/ufl/algorithms/check_arities.py:170, in check_integrand_arity(expr, arguments, complex_mode)
168 for arg, conj in arg_tuples:
169 if arg.number() == 0 and not conj:
--> 170 raise ArityMismatch("Failure to conjugate test function in complex Form")
171 elif arg.number() > 0 and conj:
172 raise ArityMismatch("Argument {0} is spuriously conjugated in complex Form".format(arg))
ArityMismatch: Failure to conjugate test function in complex Form
Yes, I know the same topic is already discussed in the website; I saw this page
to fix the same issue when I was studying with Heat equation demo code, by using ufl.inner. However, it seems that there is no such point in hyperelasticity code, so I am stuck.
Could you help me to solve this issue, and how to cope with such type of error in general?
Thanks a lot.
