Return UFL object from function

UFL
,
1

If I want to package my often repeated UFL objects in functions I notice the system becomes singular, which I guess is because UFL isn’t getting the objects together properly. Does anyone know a ‘correct way’ to do that, or is it ill-advised?

As an MWE I can take the hyperelasticity demo:

from dolfinx import log, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
#import pyvista
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.functionspace(domain, ("Lagrange", 2, (domain.geometry.dim, )))

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)

# 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])

u_bc = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)
left_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(1))
bcs = [fem.dirichletbc(u_bc, left_dofs, V)]
B = fem.Constant(domain, default_scalar_type((0, 0, 0)))
T = fem.Constant(domain, default_scalar_type((0, 0, 0)))
v = ufl.TestFunction(V)
u = fem.Function(V)
# Elasticity parameters
E = default_scalar_type(1.0e4)
nu = default_scalar_type(0.3)
mu = fem.Constant(domain, E / (2 * (1 + nu)))
lmbda = fem.Constant(domain, E * nu / ((1 + nu) * (1 - 2 * nu))) 
metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)

this represents the construction of the problem parameters.

Now if the following is left as-is it solves the original problem. But if I change False to True it solves the problem with variables generated by functions:

metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
if False:
    def F(u):
        # Spatial dimension
        d = len(u)
    
        # Identity tensor
        I = ufl.variable(ufl.Identity(d))
    
        # Deformation gradient
        return ufl.variable(I + ufl.grad(u))

    def psi(u): 
        dudx = F(u)
        # Invariants of deformation tensors
        Ic = ufl.variable(ufl.tr(dudx.T*dudx))
        J = ufl.variable(ufl.det(dudx))
        return (mu / 2) * (Ic - 3) - mu * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J))**2
        
        
    def PK1(u):
        return ufl.diff(psi(u), F(u))
    # Define form F (we want to find u such that F(u) = 0)
    Form = ufl.inner(ufl.grad(v), PK1(u)) * dx - ufl.inner(v, B) * dx - ufl.inner(v, T) * ds(2)
else:
    # 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))
    # 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)

    Form = ufl.inner(ufl.grad(v), P) * dx - ufl.inner(v, B) * dx - ufl.inner(v, T) * ds(2)

and the resulting matrix is singular.


problem = NonlinearProblem(Form, u, bcs,)
solver = NewtonSolver(domain.comm, problem)

ksp = solver.krylov_solver
opts = PETSc.Options()  # type: ignore
option_prefix = ksp.getOptionsPrefix()
sys = PETSc.Sys()  # type: ignore

opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu" 
if sys.hasExternalPackage("mumps"):
    opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
# Set Newton solver options
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"

log.set_log_level(log.LogLevel.INFO)
tval0 = -1.5
for n in range(1, 10):
    T.value[2] = n * tval0
    num_its, converged = solver.solve(u)
    assert (converged)
    u.x.scatter_forward()
    

[2025-09-29 14:09:35.688] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-09-29 14:09:35.688] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-09-29 14:09:35.688] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-09-29 14:09:35.688] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-09-29 14:09:35.692] [info] Requesting connectivity (2, 0) - (3, 0)
[2025-09-29 14:09:35.692] [info] Requesting connectivity (3, 0) - (2, 0)
[2025-09-29 14:09:35.692] [info] Column ghost size increased from 0 to 0
[2025-09-29 14:09:35.694] [info] PETSc Krylov solver starting to solve system.
[2025-09-29 14:09:35.694] [info] PETSc Krylov solver starting to solve system.
[2025-09-29 14:09:35.695] [info] Newton iteration 2: r (abs) = inf (tol = 1e-08), r (rel) = -nan (tol = 1e-08)
3

This is because whenever you call ufl.variable, you create a new variable, thus the dudx in psi is has a different variable than the one from F(u) in ufl.diff(psi(u), F(u)).

In theory (unfortunately, not in pratice, ref: Expose `ufl.Label` by jorgensd · Pull Request #423 · FEniCS/ufl · GitHub) one could make a unique label to pass in to each instance of ufl.operators.Variable(F(u), label), that would render them the same objects under equality and hopefully under differentiation (to be checked).

4

OK Thanks. A shame because in higher order time-stepping it would make the code much more compact without sacrificing readability, and reduce typos. Perhaps I could implement a list structure that would point to the objects correctly.

5

@mg-tub once Expose `ufl.Label` by jorgensd · Pull Request #423 · FEniCS/ufl · GitHub is merged, you can create variables with a specific label, making them identical, i.e.

    space = FunctionSpace(domain, element)
    u = TrialFunction(space)
    F = sin(u)
    F_var = variable(F)
    C = F_var**2

    # Check that standard variable gives correct result
    ref_expr = 2.0 * F_var
    dCdF = apply_derivatives(diff(C, F_var))
    assert dCdF == ref_expr

    # Check that variable with same label gives same result
    F_var_2 = Variable(F, label=F_var.label())
    dCdF_2 = apply_derivatives(diff(C, F_var_2))
    assert dCdF_2 == ref_expr

    # Check that new variable gives 0 result
    F_var_3 = variable(F)
    dCdF_3 = apply_derivatives(diff(C, F_var_3))
    assert dCdF_3 == 0

    # Check that variable with different label gives 0 result
    F_var_4 = Variable(F, label=Label(888))
    dCdF_4 = apply_derivatives(diff(C, F_var_4))
    assert dCdF_4 == 0
6

@mg-tub I just released a new version 2025.2.0 of UFL that exposes the ufl.Label, so that you can redefine a variable, with the same level, and achieve equivalence when differentiating with respect to different Python objects .

7

Oh very nice. Thanks a lot. Any idea how long that takes to go into conda-forge and the various containers e.g. spack / docker ?

8

Docker will happen tomorrow morning.

Conda hopefully tomorrow (Need to merge: fenics-ufl v2025.2.0 by regro-cf-autotick-bot · Pull Request #13 · conda-forge/fenics-ufl-feedstock · GitHub) and then it takes a few hours.

Spack usually takes a bit longer to get the tagged version, but you can install it as spack add py-fenics-ufl@main and it should work.

1 Like