How to resolve ArityMismatch error in the nonlinear problem

General
1

Following is my code for a nonlinear spatial timoshenko beam, I defined a 1d beam embedded in 3d space. I am getting ArityMismatch error from this code, which I do not not how to resolve.

from dolfinx import mesh, io, log, default_scalar_type, fem, plot
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
import numpy as np
import meshio
from mpi4py import MPI
import ufl
from basix.ufl import element, mixed_element
from ufl import (Measure, SpatialCoordinate, TestFunctions, TrialFunctions, div, exp, inner, grad, as_matrix, as_vector, dot, Dx, dx)



# Mesh parameters
length = 1.0
N = 5
points = np.zeros((N+1, 3))
points[:, 0] = np.linspace(0, length, N+1)
cells = np.array([[i, i+1] for i in range(N)])

# Create and save the mesh in meshio format
meshio.write("mesh.xdmf", meshio.Mesh(points, {"line": cells}))

# Load mesh in dolfinx
with io.XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as infile:
    beam_mesh = infile.read_mesh(name="Grid")

# Print mesh information to verify
print(f"Number of vertices: {beam_mesh.topology.index_map(0).size_local}")
print(f"Number of edges: {beam_mesh.topology.index_map(1).size_local}")

zeta = element("CG", beam_mesh.basix_cell(), degree=2, shape=(3,))
theta = element("CG", beam_mesh.basix_cell(), degree=2, shape=(3,))
kappa = element("CG", beam_mesh.basix_cell(), degree=2, shape=(3,))
phi = element("CG", beam_mesh.basix_cell(), degree=2, shape=(3,))
rotation = element("CG", beam_mesh.basix_cell(), degree=2, shape=(3,3))
displacement = element("CG", beam_mesh.basix_cell(), degree=2, shape=(3,))

element = mixed_element([zeta, theta])
V = fem.functionspace(beam_mesh, element)
#V_r = fem.functionspace(beam_mesh, rotation)
#V_kappa = fem.functionspace(beam_mesh, kappa)
#V_phi = fem.functionspace(beam_mesh, phi)
#V_u = fem.functionspace(beam_mesh, displacement)

u = fem.Function(V)
z, t = ufl.split(u)
#print(len(z))

def rotation_matrix(theta):
    """Constructs a rotation matrix using the vector theta."""
    theta_x, theta_y, theta_z = theta[0], theta[1], theta[2]
    # Approximate rotation matrix for small rotations
    return as_matrix([[1, -theta_z, theta_y],
                      [theta_z, 1, -theta_x],
                      [-theta_y, theta_x, 1]])

# Compute the rotation matrix R from theta
#R = fem.Function(V_r)
#print(R.ufl_shape)
R = rotation_matrix(t)

def curvature_vector(R):
    dr = Dx(R,0)
    k = dot(R.T, dr)
    #print(phi.ufl_shape)
    return as_vector([k[2,1], k[0,2], k[1,0]])

#k_vec = fem.Function(V_kappa)
k_vec = curvature_vector(R)
#print(phi_vec.ufl_shape)

def distortion_vector(R, zeta):
   # Constructs a vector to define distortion vector
    dzeta_1 = Dx(zeta[0], 0)
    dzeta_2 = Dx(zeta[1], 0)
    dzeta_3 = Dx(zeta[2], 0)
    dzeta = as_vector([1+dzeta_1, dzeta_2, dzeta_3])
    
    return (dot(R.T,dzeta))

phi_vec = distortion_vector(R,z)
print(len(phi_vec))

def d_zeta(zeta):
   # Constructs a vector to define distortion vector
    dzeta_1 = Dx(zeta[0], 0)
    dzeta_2 = Dx(zeta[1], 0)
    dzeta_3 = Dx(zeta[2], 0)
    
    return as_vector([dzeta_1, dzeta_2, dzeta_3])

dz = d_zeta(z)

def d_zeta_bar(zeta):
   # Constructs a vector to define distortion vector
    dzeta_1 = Dx(zeta[0], 0)
    dzeta_2 = Dx(zeta[1], 0)
    dzeta_3 = Dx(zeta[2], 0)
    
    return as_vector([1+dzeta_1, dzeta_2, dzeta_3])

x= SpatialCoordinate(beam_mesh)

x_s = as_vector([0, x[1], x[2]])

def skew_matrix(theta):
    """Constructs a rotation matrix using the vector theta."""
    theta_x, theta_y, theta_z = theta[0], theta[1], theta[2]
    # Approximate rotation matrix for small rotations
    return as_matrix([[0, -theta_z, theta_y],
                      [theta_z, 0, -theta_x],
                      [-theta_y, theta_x, 0]])

#S = fem.Function(V_r)
S = skew_matrix(t)

#ud = fem.Function(V_u)
ud = z + dot(S,x_s)
#print(len(ud))

#testfunctions
(z_,t_)  = TestFunctions(V)

#dR = fem.Function(V_r)
#print(R.ufl_shape)
dR = rotation_matrix(t_)
print(dR.ufl_shape)

dz_ = d_zeta(z_)
dzb_ = d_zeta_bar(z_)
dt_ = d_zeta(t_)

#dk_ = TestFunctions(V_kappa)
dk_ = dot(R.T, dt_)
dS = skew_matrix(t_)
dSzb = skew_matrix(dzb_)
#dphi_ = TestFunctions(V_phi)
dphi_ = dot(R.T, dz_ + dot(dSzb,t_))
print(len(dphi_))

#dud = TestFunctions(V_u)
dud = z_ + dot(dS,x_s)






strain_n = as_vector([dz[0] + 0.5 * (dz[0] **2 + dz[1] **2 + dz[2] **2),
                      phi_vec[0]*k_vec[2] - k_vec[0]*phi_vec[2],
                      phi_vec[0]*k_vec[1] - k_vec[0]*phi_vec[1],
                      0.5 * (k_vec[0] **2 + k_vec[2] **2), 
                      0.5 * (k_vec[0] **2 + k_vec[1] **2), 
                      k_vec[1]*k_vec[2]]) 
strain_s = as_vector([phi_vec[1], phi_vec[2], k_vec[0]])
stress_n = dot(E_n_tensor, strain_n)
stress_s = dot(E_s_tensor, strain_s)

dstrain_n = as_vector([dz_[0] + dz_[0]*dz[0] + dz_[1]*dz[1] + dz_[2]*dz[2],
                       dk_[2]*phi_vec[0] + dphi_[0]*k_vec[2] - dk_[0]*phi_vec[2] - dphi_[2]*k_vec[0],
                       dk_[1]*phi_vec[0] + dphi_[0]*k_vec[1] - dk_[0]*phi_vec[1] - dphi_[1]*k_vec[0],
                       dk_[0]*k_vec[0] + dk_[2]*k_vec[2],
                       dk_[0]*k_vec[0] + dk_[1]*k_vec[1],
                       dk_[1]*k_vec[2] + dk_[2]*k_vec[1]]
                       )
dstrain_s = as_vector([dphi_[1], dphi_[2], dk_[0]])

#boundary terms

#Next, we define the body force on the reference configuration (B), and nominal (first Piola-Kirchhoff) traction (T).

B = fem.Constant(beam_mesh, default_scalar_type((0, 0, 0)))
f = fem.Constant(beam_mesh, default_scalar_type((0, 0, 0)))

#defining two python functions to determine the facets to apply boundary conditions

def left(x):
    return np.isclose(x[0], 0)


def right(x):
    return np.isclose(x[0], length)


fdim = beam_mesh.topology.dim - 1
facets_left = mesh.locate_entities_boundary(beam_mesh, fdim, left_end)
V_z, _ = V.sub(0).collapse()  # Subspace for displacement
V_theta, _ = V.sub(1).collapse()  # Subspace for rotation

dofs_z_left = fem.locate_dofs_topological((V.sub(0), V_z), fdim, facets_left)
dofs_theta_left = fem.locate_dofs_topological((V.sub(1), V_theta), fdim, facets_left)

# Define zero-valued functions for the boundary conditions
zero_z = fem.Function(V_z)  # Function for displacement boundary condition (w = 0)
zero_theta = fem.Function(V_theta)  # Function for rotation boundary condition (phi = 0)
zero_z.interpolate(lambda x: np.zeros((3, x.shape[0])))  
zero_theta.interpolate(lambda x: np.zeros((3, x.shape[0])))  

bc_z_left = fem.dirichletbc(zero_z, dofs_z_left, V.sub(0))  # BC for displacement
bc_theta_left = fem.dirichletbc(zero_theta, dofs_theta_left, V.sub(1))  # BC for rotation

bcs = [bc_z_left, bc_theta_left]

facets_right = mesh.locate_entities_boundary(beam_mesh, fdim, right_end)

metadata = {"quadrature_degree": 4}
#ds = Measure("ds", domain=mesh, subdomain_data=facets, metadata=metadata)
dx = Measure("dx", domain=beam_mesh, metadata=metadata)

# Check if any facets were found at the right end
if len(facets_right) == 0:
    raise ValueError("No boundary facets found at the right end of the beam.")

# Create an array of markers (use 1 for the right end)
values_right = np.full(len(facets_right), 1, dtype=np.int32)

# Create the meshtags object for the boundary
boundary_markers = mesh.meshtags(beam_mesh, fdim, facets_right, values_right)

# Define boundary measure with subdomain markers
ds = Measure("ds", domain=beam_mesh, subdomain_data=boundary_markers, metadata=metadata)

# Define the body force term in the weak form
body_force_term = inner(dud, B) * dx 
traction_term = inner(dud, f) * ds(1)

F = inner(dstrain_n, stress_n) * dx
F+= inner(dstrain_s, stress_s) * dx
F-= body_force_term
F-= traction_term

problem = NonlinearProblem(F, ud, bcs=bcs, jit_options = {"timeout": 300})

Following is the error message

ArityMismatch                             Traceback (most recent call last)
Cell In[95], line 1
----> 1 problem = NonlinearProblem(F, u, bcs=bcs, jit_options = {"timeout": 300})

File /opt/tljh/user/lib/python3.10/site-packages/dolfinx/fem/petsc.py:897, in NonlinearProblem.__init__(self, F, u, bcs, J, form_compiler_options, jit_options)
    868 def __init__(
    869     self,
    870     F: ufl.form.Form,
   (...)
    875     jit_options: typing.Optional[dict] = None,
    876 ):
    877     """Initialize solver for solving a non-linear problem using Newton's method`.
    878 
    879     Args:
   (...)
    895 
    896     """
--> 897     self._L = _create_form(
    898         F, form_compiler_options=form_compiler_options, jit_options=jit_options
    899     )
    901     # Create the Jacobian matrix, dF/du
    902     if J is None:

File /opt/tljh/user/lib/python3.10/site-packages/dolfinx/fem/forms.py:249, in form(form, dtype, form_compiler_options, jit_options, entity_maps)
    246         return list(map(lambda sub_form: _create_form(sub_form), form))
    247     return form
--> 249 return _create_form(form)

File /opt/tljh/user/lib/python3.10/site-packages/dolfinx/fem/forms.py:244, in form.<locals>._create_form(form)
    241 """Recursively convert ufl.Forms to dolfinx.fem.Form, otherwise
    242 return form argument"""
    243 if isinstance(form, ufl.Form):
--> 244     return _form(form)
    245 elif isinstance(form, collections.abc.Iterable):
    246     return list(map(lambda sub_form: _create_form(sub_form), form))

File /opt/tljh/user/lib/python3.10/site-packages/dolfinx/fem/forms.py:186, in form.<locals>._form(form)
    184 if mesh is None:
    185     raise RuntimeError("Expecting to find a Mesh in the form.")
--> 186 ufcx_form, module, code = jit.ffcx_jit(
    187     mesh.comm, form, form_compiler_options=form_compiler_options, jit_options=jit_options
    188 )
    190 # For each argument in form extract its function space
    191 V = [arg.ufl_function_space()._cpp_object for arg in form.arguments()]

File /opt/tljh/user/lib/python3.10/site-packages/dolfinx/jit.py:51, in mpi_jit_decorator.<locals>.mpi_jit(comm, *args, **kwargs)
     47 @functools.wraps(local_jit)
     48 def mpi_jit(comm, *args, **kwargs):
     49     # Just call JIT compiler when running in serial
     50     if comm.size == 1:
---> 51         return local_jit(*args, **kwargs)
     53     # Default status (0 == ok, 1 == fail)
     54     status = 0

File /opt/tljh/user/lib/python3.10/site-packages/dolfinx/jit.py:201, in ffcx_jit(ufl_object, form_compiler_options, jit_options)
    199 # Switch on type and compile, returning cffi object
    200 if isinstance(ufl_object, ufl.Form):
--> 201     r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
    202 elif isinstance(ufl_object, ufl.AbstractFiniteElement):
    203     r = ffcx.codegeneration.jit.compile_elements([ufl_object], options=p_ffcx, **p_jit)

File /opt/tljh/user/lib/python3.10/site-packages/ffcx/codegeneration/jit.py:256, in compile_forms(forms, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries, visualise)
    253     for name in form_names:
    254         decl += form_template.format(name=name)
--> 256     impl = _compile_objects(
    257         decl,
    258         forms,
    259         form_names,
    260         module_name,
    261         p,
    262         cache_dir,
    263         cffi_extra_compile_args,
    264         cffi_verbose,
    265         cffi_debug,
    266         cffi_libraries,
    267         visualise=visualise,
    268     )
    269 except Exception as e:
    270     try:
    271         # remove c file so that it will not timeout next time

File /opt/tljh/user/lib/python3.10/site-packages/ffcx/codegeneration/jit.py:383, in _compile_objects(decl, ufl_objects, object_names, module_name, options, cache_dir, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries, visualise)
    379 libraries = _libraries + cffi_libraries if cffi_libraries is not None else _libraries
    381 # JIT uses module_name as prefix, which is needed to make names of all struct/function
    382 # unique across modules
--> 383 _, code_body = ffcx.compiler.compile_ufl_objects(
    384     ufl_objects, prefix=module_name, options=options, visualise=visualise
    385 )
    387 ffibuilder = cffi.FFI()
    389 ffibuilder.set_source(
    390     module_name,
    391     code_body,
   (...)
    394     libraries=libraries,
    395 )

File /opt/tljh/user/lib/python3.10/site-packages/ffcx/compiler.py:108, in compile_ufl_objects(ufl_objects, options, object_names, prefix, visualise)
    106 # Stage 1: analysis
    107 cpu_time = time()
--> 108 analysis = analyze_ufl_objects(ufl_objects, options["scalar_type"])  # type: ignore
    109 _print_timing(1, time() - cpu_time)
    111 # Stage 2: intermediate representation

File /opt/tljh/user/lib/python3.10/site-packages/ffcx/analysis.py:94, in analyze_ufl_objects(ufl_objects, scalar_type)
     91     else:
     92         raise TypeError("UFL objects not recognised.")
---> 94 form_data = tuple(_analyze_form(form, scalar_type) for form in forms)
     95 for data in form_data:
     96     elements += data.unique_sub_elements

File /opt/tljh/user/lib/python3.10/site-packages/ffcx/analysis.py:94, in <genexpr>(.0)
     91     else:
     92         raise TypeError("UFL objects not recognised.")
---> 94 form_data = tuple(_analyze_form(form, scalar_type) for form in forms)
     95 for data in form_data:
     96     elements += data.unique_sub_elements

File /opt/tljh/user/lib/python3.10/site-packages/ffcx/analysis.py:180, in _analyze_form(form, scalar_type)
    177 complex_mode = np.issubdtype(scalar_type, np.complexfloating)
    179 # Compute form metadata
--> 180 form_data = ufl.algorithms.compute_form_data(
    181     form,
    182     do_apply_function_pullbacks=True,
    183     do_apply_integral_scaling=True,
    184     do_apply_geometry_lowering=True,
    185     preserve_geometry_types=(ufl.classes.Jacobian,),
    186     do_apply_restrictions=True,
    187     do_append_everywhere_integrals=False,  # do not add dx integrals to dx(i) in UFL
    188     complex_mode=complex_mode,
    189 )
    191 # Determine unique quadrature degree and quadrature scheme
    192 # per each integral data
    193 for id, integral_data in enumerate(form_data.integral_data):
    194     # Iterate through groups of integral data. There is one integral
    195     # data for all integrals with same domain, itype, subdomain_id
   (...)
    199     # all integrals in this integral data group, i.e. must be the
    200     # same for for the same (domain, itype, subdomain_id)

File /opt/tljh/user/lib/python3.10/site-packages/ufl/algorithms/compute_form_data.py:427, in compute_form_data(form, do_apply_function_pullbacks, do_apply_integral_scaling, do_apply_geometry_lowering, preserve_geometry_types, do_apply_default_restrictions, do_apply_restrictions, do_estimate_degrees, do_append_everywhere_integrals, complex_mode)
    424 preprocessed_form = reconstruct_form_from_integral_data(self.integral_data)
    426 # TODO: Test how fast this is
--> 427 check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode)
    429 # TODO: This member is used by unit tests, change the tests to
    430 # remove this!
    431 self.preprocessed_form = preprocessed_form

File /opt/tljh/user/lib/python3.10/site-packages/ufl/algorithms/check_arities.py:211, in check_form_arity(form, arguments, complex_mode)
    209 """Check the arity of a form."""
    210 for itg in form.integrals():
--> 211     check_integrand_arity(itg.integrand(), arguments, complex_mode)

File /opt/tljh/user/lib/python3.10/site-packages/ufl/algorithms/check_arities.py:192, in check_integrand_arity(expr, arguments, complex_mode)
    190 arguments = tuple(sorted(set(arguments), key=lambda x: (x.number(), x.part())))
    191 rules = ArityChecker(arguments)
--> 192 arg_tuples = map_expr_dag(rules, expr, compress=False)
    193 args = tuple(a[0] for a in arg_tuples)
    194 if args != arguments:

File /opt/tljh/user/lib/python3.10/site-packages/ufl/corealg/map_dag.py:35, in map_expr_dag(function, expression, compress, vcache, rcache)
     15 def map_expr_dag(function, expression, compress=True, vcache=None, rcache=None):
     16     """Apply a function to each subexpression node in an expression DAG.
     17 
     18     If the same function is called multiple times in a transformation
   (...)
     33         The result of the final function call
     34     """
---> 35     (result,) = map_expr_dags(
     36         function, [expression], compress=compress, vcache=vcache, rcache=rcache
     37     )
     38     return result

File /opt/tljh/user/lib/python3.10/site-packages/ufl/corealg/map_dag.py:103, in map_expr_dags(function, expressions, compress, vcache, rcache)
    101     r = handlers[v._ufl_typecode_](v)
    102 else:
--> 103     r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
    105 # Optionally check if r is in rcache, a memory optimization
    106 # to be able to keep representation of result compact
    107 if compress:

File /opt/tljh/user/lib/python3.10/site-packages/ufl/algorithms/check_arities.py:57, in ArityChecker.sum(self, o, a, b)
     55 """Apply to sum."""
     56 if a != b:
---> 57     raise ArityMismatch(
     58         f"Adding expressions with non-matching form arguments {_afmt(a)} vs {_afmt(b)}."
     59     )
     60 return a

ArityMismatch: Adding expressions with non-matching form arguments () vs ('v_0',).
2

anyone have any solution or guidelines to this?

3

Unfortunately your example is quite long which makes it difficult for most advanced users to justify spending time debugging all your code.

I’d suggest a process of elimination. Reduce your code (particularly the finite element formulation) as much as you can until you locate the precise source of the problem. The arity mismatch error usually implies an incorrect combination of test and trial functions in your formulation.

1 Like
4

okay, I figured out where the problem is happening, I am attaching that particular portion below

# Mesh parameters
length = 1.0
N = 5

L = 1.0
points = np.array([[0.0, 0.0, 0.0], [L, 0.0, 0.0]])
# Create a 1D mesh for the beam
beam_mesh = mesh.create_interval(MPI.COMM_WORLD, 5, points[:,0])

z = element("CG", beam_mesh.basix_cell(), degree=1, shape=(3,))
theta = element("CG", beam_mesh.basix_cell(), degree=1, shape=(3,))

element = mixed_element([z, theta])

V = fem.functionspace(beam_mesh, element)

# Create the function space and functions from the mixed element
u = fem.Function(V)
z,t = ufl.split(u)
def rotation_matrix(theta):
    """Constructs a rotation matrix using the vector theta."""
    theta_x, theta_y, theta_z = theta[0], theta[1], theta[2]
    # Approximate rotation matrix for small rotations
    return as_matrix([[1, -theta_z, theta_y],
                      [theta_z, 1, -theta_x],
                      [-theta_y, theta_x, 1]])

# Compute the rotation matrix R from theta

R = rotation_matrix(t)
#print(R.ufl_shape)

def curvature_vector(R):
    dr = Dx(R,0)
    k = dot(R.T, dr)
    #print(phi.ufl_shape)
    return as_vector([k[2,1], k[0,2], k[1,0]])


k_vec = curvature_vector(R)

def distortion_vector(R, zeta):
   # Constructs a vector to define distortion vector
    dzeta_1 = Dx(zeta[0], 0)
    dzeta_2 = Dx(zeta[1], 0)
    dzeta_3 = Dx(zeta[2], 0)
    dzeta = as_vector([1+dzeta_1, dzeta_2, dzeta_3])
    
    return (dot(R.T,dzeta))

phi_vec = distortion_vector(R,z)
print(len(phi_vec))

def d_zeta(zeta):
   # Constructs a vector to define distortion vector
    dzeta_1 = Dx(zeta[0], 0)
    dzeta_2 = Dx(zeta[1], 0)
    dzeta_3 = Dx(zeta[2], 0)
    
    return as_vector([dzeta_1, dzeta_2, dzeta_3])

dz = ufl.variable(d_zeta(z))


#testfunctions
v = ufl.TestFunction(V)
z_, t_ = ufl.split(v)
print(t_.ufl_shape)
dR = rotation_matrix(t_)
print(dR.ufl_shape)

dz_ = d_zeta(z_)
dzb_ = d_zeta_bar(z_)
dt_ = d_zeta(t_)

#dk_ = TestFunctions(V_kappa)
#dk_ = dot(dR.T, dt_)
dk_ = curvature_vector(dR)
dS = skew_matrix(t_)
dSzb = skew_matrix(dzb_)
#dphi_ = TestFunctions(V_phi)
dphi_ = dot(dR.T, dz_ + dot(dSzb,t_))

here I am deriving the vectors, dphi_ and dk_ based on the test functions z_ and t_, even thoigh all the other vectors are working fine, these two vectors are not. I am getting value error stating too many values to unpack (expected 2). While the same functions I derived from function form z,t are working fine.

5

Please add all your import statements, so that others can run the code and reproduce the error.

6

Here you go

from dolfinx import mesh, io, log, default_scalar_type, fem, plot
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem
from dolfinx.nls.petsc import NewtonSolver
import numpy as np
import meshio
from mpi4py import MPI
import ufl
from basix.ufl import element, mixed_element
from ufl import (Measure, SpatialCoordinate, TestFunctions, TrialFunctions, div, exp, inner, grad, as_matrix, as_vector, dot, Dx, dx)

# Mesh parameters
length = 1.0
N = 5

L = 1.0
points = np.array([[0.0, 0.0, 0.0], [L, 0.0, 0.0]])
# Create a 1D mesh for the beam
beam_mesh = mesh.create_interval(MPI.COMM_WORLD, 5, points[:,0])

z = element("CG", beam_mesh.basix_cell(), degree=1, shape=(3,))
theta = element("CG", beam_mesh.basix_cell(), degree=1, shape=(3,))

element = mixed_element([z, theta])

V = fem.functionspace(beam_mesh, element)

# Create the function space and functions from the mixed element
u = fem.Function(V)
z,t = ufl.split(u)
def rotation_matrix(theta):
    """Constructs a rotation matrix using the vector theta."""
    theta_x, theta_y, theta_z = theta[0], theta[1], theta[2]
    # Approximate rotation matrix for small rotations
    return as_matrix([[1, -theta_z, theta_y],
                      [theta_z, 1, -theta_x],
                      [-theta_y, theta_x, 1]])

# Compute the rotation matrix R from theta

R = rotation_matrix(t)
#print(R.ufl_shape)

def curvature_vector(R):
    dr = Dx(R,0)
    k = dot(R.T, dr)
    #print(phi.ufl_shape)
    return as_vector([k[2,1], k[0,2], k[1,0]])


k_vec = curvature_vector(R)

def distortion_vector(R, zeta):
   # Constructs a vector to define distortion vector
    dzeta_1 = Dx(zeta[0], 0)
    dzeta_2 = Dx(zeta[1], 0)
    dzeta_3 = Dx(zeta[2], 0)
    dzeta = as_vector([1+dzeta_1, dzeta_2, dzeta_3])
    
    return (dot(R.T,dzeta))

phi_vec = distortion_vector(R,z)
print(len(phi_vec))

def d_zeta(zeta):
   # Constructs a vector to define distortion vector
    dzeta_1 = Dx(zeta[0], 0)
    dzeta_2 = Dx(zeta[1], 0)
    dzeta_3 = Dx(zeta[2], 0)
    
    return as_vector([dzeta_1, dzeta_2, dzeta_3])

dz = ufl.variable(d_zeta(z))


#testfunctions
v = ufl.TestFunction(V)
z_, t_ = ufl.split(v)
print(t_.ufl_shape)
dR = rotation_matrix(t_)
print(dR.ufl_shape)

dz_ = d_zeta(z_)
dzb_ = d_zeta_bar(z_)
dt_ = d_zeta(t_)

#dk_ = TestFunctions(V_kappa)
#dk_ = dot(dR.T, dt_)
dk_ = curvature_vector(dR)
dS = skew_matrix(t_)
dSzb = skew_matrix(dzb_)
#dphi_ = TestFunctions(V_phi)
dphi_ = dot(dR.T, dz_ + dot(dSzb,t_))
7

Unfortunately this is not a minimal working example. The functions d_zeta_bar and skew_matrix are not provided. If I remove the calls to these functions I cannot reproduce your error.

8

sorry for the incomplete one, here I updated with the necessary functions. Hope this will work.

9

I checked the shape of dk_ and dphi_ here, which inclines with k_vec and phi_vec. yet while using dk_ and dphi_ I am getting following value_error.

ValueError                                Traceback (most recent call last)
Cell In[5], line 111
    108 u = fem.Function(V)
    110 #problem = NonlinearProblem(F, u,bcs=bcs)
--> 111 problem = NonlinearProblem(F_test_dk, u, bcs=bcs)
    112 solver = NewtonSolver(beam_mesh.comm, problem)
    113 # Set Newton solver options
    114 #solver.atol = 1e-8
    115 #solver.rtol = 1e-8
    116 #solver.convergence_criterion = "incremental"

File /opt/tljh/user/lib/python3.10/site-packages/dolfinx/fem/petsc.py:897, in NonlinearProblem.__init__(self, F, u, bcs, J, form_compiler_options, jit_options)
    868 def __init__(
    869     self,
    870     F: ufl.form.Form,
   (...)
    875     jit_options: typing.Optional[dict] = None,
    876 ):
    877     """Initialize solver for solving a non-linear problem using Newton's method`.
    878 
    879     Args:
   (...)
    895 
    896     """
--> 897     self._L = _create_form(
    898         F, form_compiler_options=form_compiler_options, jit_options=jit_options
    899     )
    901     # Create the Jacobian matrix, dF/du
    902     if J is None:

File /opt/tljh/user/lib/python3.10/site-packages/dolfinx/fem/forms.py:249, in form(form, dtype, form_compiler_options, jit_options, entity_maps)
    246         return list(map(lambda sub_form: _create_form(sub_form), form))
    247     return form
--> 249 return _create_form(form)

File /opt/tljh/user/lib/python3.10/site-packages/dolfinx/fem/forms.py:244, in form.<locals>._create_form(form)
    241 """Recursively convert ufl.Forms to dolfinx.fem.Form, otherwise
    242 return form argument"""
    243 if isinstance(form, ufl.Form):
--> 244     return _form(form)
    245 elif isinstance(form, collections.abc.Iterable):
    246     return list(map(lambda sub_form: _create_form(sub_form), form))

File /opt/tljh/user/lib/python3.10/site-packages/dolfinx/fem/forms.py:186, in form.<locals>._form(form)
    184 if mesh is None:
    185     raise RuntimeError("Expecting to find a Mesh in the form.")
--> 186 ufcx_form, module, code = jit.ffcx_jit(
    187     mesh.comm, form, form_compiler_options=form_compiler_options, jit_options=jit_options
    188 )
    190 # For each argument in form extract its function space
    191 V = [arg.ufl_function_space()._cpp_object for arg in form.arguments()]

File /opt/tljh/user/lib/python3.10/site-packages/dolfinx/jit.py:51, in mpi_jit_decorator.<locals>.mpi_jit(comm, *args, **kwargs)
     47 @functools.wraps(local_jit)
     48 def mpi_jit(comm, *args, **kwargs):
     49     # Just call JIT compiler when running in serial
     50     if comm.size == 1:
---> 51         return local_jit(*args, **kwargs)
     53     # Default status (0 == ok, 1 == fail)
     54     status = 0

File /opt/tljh/user/lib/python3.10/site-packages/dolfinx/jit.py:201, in ffcx_jit(ufl_object, form_compiler_options, jit_options)
    199 # Switch on type and compile, returning cffi object
    200 if isinstance(ufl_object, ufl.Form):
--> 201     r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
    202 elif isinstance(ufl_object, ufl.AbstractFiniteElement):
    203     r = ffcx.codegeneration.jit.compile_elements([ufl_object], options=p_ffcx, **p_jit)

File /opt/tljh/user/lib/python3.10/site-packages/ffcx/codegeneration/jit.py:276, in compile_forms(forms, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries, visualise)
    274     except Exception:
    275         pass
--> 276     raise e
    278 obj, module = _load_objects(cache_dir, module_name, form_names)
    279 return obj, module, (decl, impl)

File /opt/tljh/user/lib/python3.10/site-packages/ffcx/codegeneration/jit.py:256, in compile_forms(forms, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries, visualise)
    253     for name in form_names:
    254         decl += form_template.format(name=name)
--> 256     impl = _compile_objects(
    257         decl,
    258         forms,
    259         form_names,
    260         module_name,
    261         p,
    262         cache_dir,
    263         cffi_extra_compile_args,
    264         cffi_verbose,
    265         cffi_debug,
    266         cffi_libraries,
    267         visualise=visualise,
    268     )
    269 except Exception as e:
    270     try:
    271         # remove c file so that it will not timeout next time

File /opt/tljh/user/lib/python3.10/site-packages/ffcx/codegeneration/jit.py:383, in _compile_objects(decl, ufl_objects, object_names, module_name, options, cache_dir, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries, visualise)
    379 libraries = _libraries + cffi_libraries if cffi_libraries is not None else _libraries
    381 # JIT uses module_name as prefix, which is needed to make names of all struct/function
    382 # unique across modules
--> 383 _, code_body = ffcx.compiler.compile_ufl_objects(
    384     ufl_objects, prefix=module_name, options=options, visualise=visualise
    385 )
    387 ffibuilder = cffi.FFI()
    389 ffibuilder.set_source(
    390     module_name,
    391     code_body,
   (...)
    394     libraries=libraries,
    395 )

File /opt/tljh/user/lib/python3.10/site-packages/ffcx/compiler.py:108, in compile_ufl_objects(ufl_objects, options, object_names, prefix, visualise)
    106 # Stage 1: analysis
    107 cpu_time = time()
--> 108 analysis = analyze_ufl_objects(ufl_objects, options["scalar_type"])  # type: ignore
    109 _print_timing(1, time() - cpu_time)
    111 # Stage 2: intermediate representation

File /opt/tljh/user/lib/python3.10/site-packages/ffcx/analysis.py:94, in analyze_ufl_objects(ufl_objects, scalar_type)
     91     else:
     92         raise TypeError("UFL objects not recognised.")
---> 94 form_data = tuple(_analyze_form(form, scalar_type) for form in forms)
     95 for data in form_data:
     96     elements += data.unique_sub_elements

File /opt/tljh/user/lib/python3.10/site-packages/ffcx/analysis.py:94, in <genexpr>(.0)
     91     else:
     92         raise TypeError("UFL objects not recognised.")
---> 94 form_data = tuple(_analyze_form(form, scalar_type) for form in forms)
     95 for data in form_data:
     96     elements += data.unique_sub_elements

File /opt/tljh/user/lib/python3.10/site-packages/ffcx/analysis.py:180, in _analyze_form(form, scalar_type)
    177 complex_mode = np.issubdtype(scalar_type, np.complexfloating)
    179 # Compute form metadata
--> 180 form_data = ufl.algorithms.compute_form_data(
    181     form,
    182     do_apply_function_pullbacks=True,
    183     do_apply_integral_scaling=True,
    184     do_apply_geometry_lowering=True,
    185     preserve_geometry_types=(ufl.classes.Jacobian,),
    186     do_apply_restrictions=True,
    187     do_append_everywhere_integrals=False,  # do not add dx integrals to dx(i) in UFL
    188     complex_mode=complex_mode,
    189 )
    191 # Determine unique quadrature degree and quadrature scheme
    192 # per each integral data
    193 for id, integral_data in enumerate(form_data.integral_data):
    194     # Iterate through groups of integral data. There is one integral
    195     # data for all integrals with same domain, itype, subdomain_id
   (...)
    199     # all integrals in this integral data group, i.e. must be the
    200     # same for for the same (domain, itype, subdomain_id)

File /opt/tljh/user/lib/python3.10/site-packages/ufl/algorithms/compute_form_data.py:427, in compute_form_data(form, do_apply_function_pullbacks, do_apply_integral_scaling, do_apply_geometry_lowering, preserve_geometry_types, do_apply_default_restrictions, do_apply_restrictions, do_estimate_degrees, do_append_everywhere_integrals, complex_mode)
    424 preprocessed_form = reconstruct_form_from_integral_data(self.integral_data)
    426 # TODO: Test how fast this is
--> 427 check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode)
    429 # TODO: This member is used by unit tests, change the tests to
    430 # remove this!
    431 self.preprocessed_form = preprocessed_form

File /opt/tljh/user/lib/python3.10/site-packages/ufl/algorithms/check_arities.py:211, in check_form_arity(form, arguments, complex_mode)
    209 """Check the arity of a form."""
    210 for itg in form.integrals():
--> 211     check_integrand_arity(itg.integrand(), arguments, complex_mode)

File /opt/tljh/user/lib/python3.10/site-packages/ufl/algorithms/check_arities.py:192, in check_integrand_arity(expr, arguments, complex_mode)
    190 arguments = tuple(sorted(set(arguments), key=lambda x: (x.number(), x.part())))
    191 rules = ArityChecker(arguments)
--> 192 arg_tuples = map_expr_dag(rules, expr, compress=False)
    193 args = tuple(a[0] for a in arg_tuples)
    194 if args != arguments:

File /opt/tljh/user/lib/python3.10/site-packages/ufl/corealg/map_dag.py:35, in map_expr_dag(function, expression, compress, vcache, rcache)
     15 def map_expr_dag(function, expression, compress=True, vcache=None, rcache=None):
     16     """Apply a function to each subexpression node in an expression DAG.
     17 
     18     If the same function is called multiple times in a transformation
   (...)
     33         The result of the final function call
     34     """
---> 35     (result,) = map_expr_dags(
     36         function, [expression], compress=compress, vcache=vcache, rcache=rcache
     37     )
     38     return result

File /opt/tljh/user/lib/python3.10/site-packages/ufl/corealg/map_dag.py:103, in map_expr_dags(function, expressions, compress, vcache, rcache)
    101     r = handlers[v._ufl_typecode_](v)
    102 else:
--> 103     r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
    105 # Optionally check if r is in rcache, a memory optimization
    106 # to be able to keep representation of result compact
    107 if compress:

File /opt/tljh/user/lib/python3.10/site-packages/ufl/algorithms/check_arities.py:78, in ArityChecker.product(self, o, a, b)
     74 for x in b:
     75     if x[0].number() in anumbers:
     76         raise ArityMismatch(
     77             "Multiplying expressions with overlapping form argument number "
---> 78             f"{x[0].number()}, argument is {_afmt(x)}."
     79         )
     80 # Combine argument lists
     81 c = tuple(sorted(set(a + b), key=lambda x: (x[0].number(), x[0].part())))

File /opt/tljh/user/lib/python3.10/site-packages/ufl/algorithms/check_arities.py:19, in _afmt(atuple)
     17 def _afmt(atuple):
     18     """Return a string representation of an arity tuple."""
---> 19     return tuple(f"conj({arg})" if conj else str(arg) for arg, conj in atuple)

File /opt/tljh/user/lib/python3.10/site-packages/ufl/algorithms/check_arities.py:19, in <genexpr>(.0)
     17 def _afmt(atuple):
     18     """Return a string representation of an arity tuple."""
---> 19     return tuple(f"conj({arg})" if conj else str(arg) for arg, conj in atuple)

ValueError: too many values to unpack (expected 2)
10

Again, please note that the cod you are quoting above is only a partial snippet of the complete code.
What stands out to me is:

You cannot solve for ud, as ud is not a dolfinx.fem.Function.
Please consider tutorials on non-linear problems, such as
https://docs.fenicsproject.org/dolfinx/v0.9.0/python/demos/demo_cahn-hilliard.html
https://jsdokken.com/dolfinx-tutorial/chapter2/hyperelasticity.html

11

that u_d was a typo there, I am providing you the code portion that I am using for test cases are given below:

from dolfinx import mesh, io, log, default_scalar_type, fem, plot
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem
from dolfinx.nls.petsc import NewtonSolver
import numpy as np
import meshio
from mpi4py import MPI
import ufl
from basix.ufl import element, mixed_element
from ufl import (Measure, SpatialCoordinate, TestFunctions, TrialFunctions, div, exp, inner, grad, as_matrix, as_vector, dot, Dx, dx)

# Mesh parameters
length = 1.0
N = 5

L = 1.0
points = np.array([[0.0, 0.0, 0.0], [L, 0.0, 0.0]])
# Create a 1D mesh for the beam
beam_mesh = mesh.create_interval(MPI.COMM_WORLD, 5, points[:,0])

z = element("CG", beam_mesh.basix_cell(), degree=1, shape=(3,))
theta = element("CG", beam_mesh.basix_cell(), degree=1, shape=(3,))

element = mixed_element([z, theta])

V = fem.functionspace(beam_mesh, element)

# Create the function space and functions from the mixed element
u = fem.Function(V)
z,t = ufl.split(u)
def rotation_matrix(theta):
    """Constructs a rotation matrix using the vector theta."""
    theta_x, theta_y, theta_z = theta[0], theta[1], theta[2]
    # Approximate rotation matrix for small rotations
    return as_matrix([[1, -theta_z, theta_y],
                      [theta_z, 1, -theta_x],
                      [-theta_y, theta_x, 1]])

# Compute the rotation matrix R from theta

R = rotation_matrix(t)
#print(R.ufl_shape)

def curvature_vector(R):
    dr = Dx(R,0)
    k = dot(R.T, dr)
    #print(phi.ufl_shape)
    return as_vector([k[2,1], k[0,2], k[1,0]])


k_vec = curvature_vector(R)

def distortion_vector(R, zeta):
   # Constructs a vector to define distortion vector
    dzeta_1 = Dx(zeta[0], 0)
    dzeta_2 = Dx(zeta[1], 0)
    dzeta_3 = Dx(zeta[2], 0)
    dzeta = as_vector([1+dzeta_1, dzeta_2, dzeta_3])
    
    return (dot(R.T,dzeta))

phi_vec = distortion_vector(R,z)
print(len(phi_vec))

def d_zeta(zeta):
   # Constructs a vector to define distortion vector
    dzeta_1 = Dx(zeta[0], 0)
    dzeta_2 = Dx(zeta[1], 0)
    dzeta_3 = Dx(zeta[2], 0)
    
    return as_vector([dzeta_1, dzeta_2, dzeta_3])

dz = ufl.variable(d_zeta(z))

def d_zeta_bar(zeta):
   # Constructs a vector to define distortion vector
    dzeta_1 = Dx(zeta[0], 0)
    dzeta_2 = Dx(zeta[1], 0)
    dzeta_3 = Dx(zeta[2], 0)
    
    return as_vector([1+dzeta_1, dzeta_2, dzeta_3])

def skew_matrix(theta):
    
    theta_x, theta_y, theta_z = theta[0], theta[1], theta[2]
    # Approximate rotation matrix for small rotations
    return as_matrix([[0, -theta_z, theta_y],
                      [theta_z, 0, -theta_x],
                      [-theta_y, theta_x, 0]])

#S = fem.Function(V_r)
S = skew_matrix(t)


#testfunctions
v = ufl.TestFunction(V)
z_, t_ = ufl.split(v)
print(t_.ufl_shape)
dR = rotation_matrix(t_)
print(dR.ufl_shape)

dz_ = d_zeta(z_)
dzb_ = d_zeta_bar(z_)
dt_ = d_zeta(t_)

#dk_ = TestFunctions(V_kappa)
#dk_ = dot(dR.T, dt_)
dk_ = curvature_vector(dR)
dS = skew_matrix(t_)
dSzb = skew_matrix(dzb_)
#dphi_ = TestFunctions(V_phi)
dphi_ = dot(dR.T, dz_ + dot(dSzb,t_))

#defining two python functions to determine the facets to apply boundary conditions

def left(x):
    return np.isclose(x[0], 1e-13)


def right(x):
    return np.isclose(x[0], length)


fdim = beam_mesh.topology.dim - 1
print(beam_mesh.topology.dim)
facets_left = mesh.locate_entities_boundary(beam_mesh, fdim, left)
V_z, _= V.sub(0).collapse()  # Subspace for displacement
V_theta, _= V.sub(1).collapse()  # Subspace for rotation


dofs_z_left = fem.locate_dofs_geometrical((V.sub(0),V_z), left)
dofs_theta_left = fem.locate_dofs_geometrical((V.sub(1),V_theta), left)

# Define zero-valued functions for the boundary conditions
zero_z = fem.Function(V_z)  # Function for displacement boundary condition (w = 0)
zero_theta = fem.Function(V_theta)  # Function for rotation boundary condition (phi = 0)
with zero_z.vector.localForm() as z_local:
    z_local.set(0.0)

with zero_theta.vector.localForm() as theta_local:
    theta_local.set(0.0)
    
#zero_z.interpolate(lambda x: np.zeros((3, x.shape[0])))  
#zero_theta.interpolate(lambda x: np.zeros((3, x.shape[0]))

bc_z_left = fem.dirichletbc(zero_z, dofs_z_left, V.sub(0))  # BC for displacement
bc_theta_left = fem.dirichletbc(zero_theta, dofs_theta_left, V.sub(1))  # BC for rotation

bcs = [bc_z_left, bc_theta_left]

facets_right = mesh.locate_entities_boundary(beam_mesh, fdim, right)

metadata = {"quadrature_degree": 4}
#ds = Measure("ds", domain=mesh, subdomain_data=facets, metadata=metadata)
dx = Measure("dx", domain=beam_mesh, metadata=metadata)

# Check if any facets were found at the right end
if len(facets_right) == 0:
    raise ValueError("No boundary facets found at the right end of the beam.")

# Create an array of markers (use 1 for the right end)
values_right = np.full(len(facets_right), 1, dtype=np.int32)

# Create the meshtags object for the boundary
boundary_markers = mesh.meshtags(beam_mesh, fdim, facets_right, values_right)

# Define boundary measure with subdomain markers
ds = Measure("ds", domain=beam_mesh, subdomain_data=boundary_markers, metadata=metadata)

#F_test = G * dstrain_n[0] * stress_n[0] * dx
#F_simple = G * inner(k_vec[2], dz_[2]) * dx  # Simplest possible form with a single component 
#F_test = inner(dphi_, dphi_) * dx
F_test_dk = inner(dz[2]**2, dk_[1]) * dx
u = fem.Function(V)

#problem = NonlinearProblem(F, u,bcs=bcs)
problem = NonlinearProblem(F_test_dk, u, bcs=bcs)
solver = NewtonSolver(beam_mesh.comm, problem)
# Set Newton solver options
#solver.atol = 1e-8
#solver.rtol = 1e-8
#solver.convergence_criterion = "incremental"
solver.solve(u)
12

u should be part of your variational form, as F_test_dk is differentiated by u to form the jacobian for solving the nonlinear problem with a newton type method. See for instance

https://jsdokken.com/dolfinx-tutorial/chapter4/newton-solver.html#definition-of-residual-and-jacobian

or
http://jsdokken.com/FEniCS-workshop/src/unified_form_language/ufl_forms.html#non-linear-functional

13

here, F_simple works fine though. What I am observing is, when I am deriving dphi_, dk_ in terms of test funtions z_ and t_ they are showing value error. but every other functions are working fine. k_vec and phi_vec are derived from function u in similar way to dk_ and dphi_ but they are working fine. dz_, which is derivative of test function z_ is wroking fine as well.

14

So, I see the issue is with test functions appearing here in nonlinear combinations, which the compiler is struggling to interpret. Is there any way to express the test functions in nonlinearcombinations? for example, when I print d_phi I get the following


dphi_0: (grad(v_0[0]))[0]
dphi_1: (grad(v_0[0]))[0] * v_0[5]
dphi_2: (grad(v_0[0]))[0] * -1 * v_0[4]

here dphi_0 works, but dphi_1 and dphi_2 do not work as they are nonlinear combinations of test functions v_0[0], v_0[1], v_0[2], etc/

15

You can’t have non-linear relationships between test functions.

Note that a weak form is either based on:

  1. Minimization of an energy J(u)
  2. Based of a strong form of a PDE

For method one, we want to find the minimum, thus we differentiate J wrt. u in some direction v, see;
https://jsdokken.com/FEniCS-workshop/src/unified_form_language/ufl_forms.html
Thus you have a linear relationship with respect to v.
Similarly, this happens through the transformation of the strong to weak/variational form,as you multiply by v and integrate over the domain of interest. The relationship is also linear