Error with fenics-ufl 2025.3.0.dev0

[marazzaf]

Hello,

I was having an issue with a mixed formulation using DG 0 elements that should have been fixed by this pull request.

However, when I install the corresponding version of fenics-ufl (2025.3.0.dev0) and run it I get an error when importing a mesh:

File "/home/marazzato/Research/Coding/Harmonic_maps/stereographic.py", line 11, in <module>
  base = Mesh(ngmesh)
File "petsc4py/PETSc/Log.pyx", line 250, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
File "petsc4py/PETSc/Log.pyx", line 251, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
File "/home/marazzato/venv-firedrake/lib/python3.10/site-packages/firedrake/mesh.py", line 3106, in Mesh
  topology = MeshTopology(plex, name=plex.getName(), reorder=reorder,
File "petsc4py/PETSc/Log.pyx", line 250, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
File "petsc4py/PETSc/Log.pyx", line 251, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
File "/home/marazzato/venv-firedrake/lib/python3.10/site-packages/firedrake/mesh.py", line 1107, in __init__
  super().__init__(plex, name, reorder, sfXB, perm_is, distribution_name, permutation_name, comm, submesh_parent=submesh_parent)
File "/home/marazzato/venv-firedrake/lib/python3.10/site-packages/firedrake/mesh.py", line 567, in __init__
  self._cell_numbering, _ = self.create_section(entity_dofs)
File "/home/marazzato/venv-firedrake/lib/python3.10/site-packages/firedrake/mesh.py", line 746, in create_section
  return dmcommon.create_section(self, nodes_per_entity, on_base=real_tensorproduct, block_size=block_size, boundary_set=boundary_set)
File "firedrake/cython/dmcommon.pyx", line 1248, in firedrake.cython.dmcommon.create_section
File "/usr/lib/python3.10/functools.py", line 981, in __get__
  val = self.func(instance)
File "/home/marazzato/venv-firedrake/lib/python3.10/site-packages/firedrake/mesh.py", line 1370, in cell_set
  size = list(self._entity_classes[self.cell_dimension(), :])
File "/home/marazzato/venv-firedrake/lib/python3.10/site-packages/firedrake/mesh.py", line 811, in cell_dimension
  return self.ufl_cell().topological_dimension()
TypeError: 'int' object is not callable

Any idea what's happening?

Thank you for your help.

[pbrubeck]

The fix should be backported to ufl release. This hasn't happened yet.

[pbrubeck]

Many people are on holiday, so this might take a while. If you would like a quick fix, I suggest you checkout the release branch and cherry-pick the commit

[marazzaf]

Hi Pablo,
Thank you for your answer.
I believe I did that.
I cloned the fenics/ufl repo on my computer and am on the main branch where your pull request was accepted.
Then I ran pip install in that folder.
When I run firedrake, I get that error now.
Any idea what's up?

[pbrubeck]

UFL main is only expected to work with Firedrake main, you get this error because you have Firedrake release and UFL main. You will not get this bug on a developer installation of Firedrake.

[pbrubeck]

My suggestion is to checkout the UFL release branch and cherry-pick the commit with the fix (assuming you pip installed UFL in editable mode)

[pbrubeck]

Try the branch from this PR FEniCS/ufl#459

[marazzaf]

There you go

from firedrake import *
from firedrake.petsc import PETSc
from firedrake.output import VTKFile
import sys

#Square domain
N = 10
H = 4
mesh = RectangleMesh(N, N, H, H, diagonal='crossed', name='mesh')

#Mesh size
h_func = Function(FunctionSpace(mesh, "DG", 0)).interpolate(CellDiameter(mesh))
hmax = h_func.dat.data.max()

#FEM space
V = VectorFunctionSpace(mesh, 'DG', 2, dim=3)
W = TensorFunctionSpace(mesh, 'DG', 0, symmetry=True)
Z = V * W

#Dirichlet boundary conditions
x = SpatialCoordinate(mesh)
g = as_vector((x[0], x[1], 0))

#Solution
sol = Function(Z, name='sol')
sol.sub(0).interpolate(g) #Storing initial guess
phi,Lag = split(sol)

#Bilinear form
y,lmbda = TrialFunctions(Z)
w,mu = TestFunctions(Z)
a = inner(grad(grad(y)), grad(grad(w))) * dx

#Gradient flow for isometry constraint
tau = hmax #time step
h = CellDiameter(mesh)
a += inner(grad(grad(y)), grad(grad(w))) / tau * dx #for H^2_h norm
a += inner(jump(grad(y)), jump(grad(w))) / avg(h) / tau * dS #for H^2_h norm
a += inner(jump(y), jump(w))/ avg(h)**3 / tau * dS #for H^2_h norm

#Penalty term for gradients
pen_0 = 5e3
pen_1 = 1100
a += pen_1 / avg(h) * inner(jump(grad(y)), jump(grad(w))) * dS
a += pen_0 / avg(h)**3 * inner(jump(y), jump(w)) * dS #for DG

#Adding the isometry constraint.
a += inner(lmbda, dot(grad(phi).T, grad(w)) + dot(grad(w).T, grad(phi))) * dx
a += inner(mu, dot(grad(phi).T, grad(y)) + dot(grad(y).T, grad(phi))) * dx

#Symmetry and consistency terms for nonconformity
n = FacetNormal(mesh)
a += -inner(avg(dot(grad(grad(y)), n)), jump(grad(w))) * dS
a += -inner(avg(dot(grad(grad(w)), n)), jump(grad(y))) * dS #symmetry term

#Imposing weakly the Homogeneous Dirichlet BC
a += pen_0 / h**3 * inner(y, w) * (ds(1) + ds(3))

#Imposing weakly the normal gradient
a += pen_1 / h * inner(dot(grad(y), n), dot(grad(w),n)) * (ds(1) + ds(3))

#RHS
#Volume load
f = Constant((0, 0, 2.5e-2))
L = inner(f, w) * dx

#Other terms
L -= inner(grad(grad(phi)), grad(grad(w))) * dx
L += inner(avg(dot(grad(grad(phi)), n)), jump(grad(w))) * dS
L += inner(avg(dot(grad(grad(w)), n)), jump(grad(phi))) * dS #symmetry term
L -= pen_1 / avg(h) * inner(jump(grad(phi)), jump(grad(w))) * dS
L -=  pen_0 / avg(h)**3 * inner(jump(phi), jump(w)) * dS

#Function to update the results
delta_sol = Function(Z)

solve(a == L, delta_sol)

[pbrubeck]

Of course. symmetry.

[marazzaf]

You think symmetry=True is the issue?

[pbrubeck]

Yes. I fixed it here 813f5de4966a07a563de911499895fa56fa50f9d

[marazzaf]

Thank you so much @pbrubeck! It works great now!