I have a axisymmetric problem with 1/r^2 singularities at the centerline. I have perused some of this forum’s thread’s on the subject :

- This old post clearly states there should be a factor of r in addition to the canonical UFL dx
- Dealing with removable singularities in cylindrical coordinates - #10 by Xuefeng_LIU seems to state that there shouldn’t be a problem since forms are evaluated on quadrature points not vertices
- Stokes equations from Cartesian to Cylindrical, results mismatch - #4 by kamilova-maths has a sliding wall at r=0 which allows for Dirichlet conditions

I have resolved to do a transformation \phi\rightarrow r\phi for my test functions in order to take care of the singularity (this seems to amount to multiplying the variational formulation by r). However, this does not produce the expected improvement ; I don’t understand how UFL simplifies some expressions but not others. Here’s a first example :

```
import dolfinx, ufl
from mpi4py.MPI import COMM_WORLD
import numpy
mesh=dolfinx.UnitSquareMesh(COMM_WORLD,10,10)
FE=ufl.FiniteElement("Lagrange",mesh.ufl_cell(),2)
V=dolfinx.FunctionSpace(mesh,FE)
u=ufl.TrialFunction(V)
v=ufl.TestFunction(V)
uD = dolfinx.Function(V)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
uD.x.scatter_forward()
fdim = mesh.topology.dim - 1
# Create facet to cell connectivity required to determine boundary facets
mesh.topology.create_connectivity(fdim, mesh.topology.dim)
boundary_facets = numpy.where(numpy.array(dolfinx.cpp.mesh.compute_boundary_facets(mesh.topology)) == 1)[0]
boundary_dofs = dolfinx.fem.locate_dofs_topological(V, fdim, boundary_facets)
bc = dolfinx.DirichletBC(uD, boundary_dofs)
x = ufl.SpatialCoordinate(mesh)
f=dolfinx.Constant(mesh,-6)
a1=ufl.inner(ufl.grad(u)/x[0],ufl.grad(v))*x[0]*ufl.dx
a2=ufl.inner(ufl.grad(u), ufl.grad(v))*ufl.dx
L=ufl.inner(f,v)*ufl.dx
problem1 = dolfinx.fem.LinearProblem(a1, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh1 = problem1.solve()
problem2 = dolfinx.fem.LinearProblem(a2, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh2 = problem2.solve()
L2_error = ufl.inner(uh1- uh2, uh1- uh2) * ufl.dx
print(numpy.sqrt(dolfinx.fem.assemble_scalar(L2_error)))
```

In Docker dolfinx with real mode this gives an error of about 10^{-15}, so machine precision. The same code with :

```
a1=ufl.inner(ufl.grad(u)/x[0]**2,ufl.grad(v))*x[0]*ufl.dx
a2=ufl.inner(ufl.grad(u)/x[0], ufl.grad(v))*ufl.dx
```

Leads to much higher error ; about 10^{-4}. In my application, I think this imperfect simplification is my dominant source of error. Anyone has an idea on how to “force” UFL to simplify a form ?