Hi everyone,

I am solving an axisymmetric problem in FEniCs and I am sure that the correct solution is achieved. From this solution, I need to calculate some quantities, and the formula for this calculation includes dividing by the radial coordinate. Because of that, in the symmetry axis, where the radial coordinate becomes zero, the obtained values do not make sense. I attach a file showing this behaviour.

I cannot post my code, so I will post a minimal example where this problem occurs, despite it is not an axisymmetric problem. In this example, the Laplace equation is solved in a unit square, being the solution u = x^2. This solution is imposed as a Dirichlet condition on the whole boundary.

from fenics import *

from dolfin import *

import matplotlib.pyplot as plt

import numpy as np

# Read mesh and subdomains

mesh = UnitSquareMesh(32, 32)

# Define finite element space

V = FunctionSpace(mesh, “CG”, 1)

# Define Dirichlet boundary

def boundary(x):

return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS or x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS

# Define variational problem

u = TrialFunction(V)

v = TestFunction(V)

f = Constant(-2.0)

a = -dot(grad(u), grad(v))*dx
L = - f*v*dx

# Boundary conditions

bc = DirichletBC(V, Expression(‘pow(x[0],2)’, degree=2), boundary)

# Compute solution

u = Function(V)

solve(a == L, u, bc)

# Save quantity J to file in XDMF format

varx = Expression(‘x[0]’, degree = 1)

W = VectorFunctionSpace(mesh, “CG”, 1, dim=2)

J = project(as_vector(grad(u)/varx),W)

with XDMFFile(mesh.mpi_comm(), ‘gradu.xdmf’) as XFile2:

XFile2.write_checkpoint(J, ‘gradu_varx’, 0)

After the solution is obtained, the quantity grad(u)/x was calculated. As u=x^2, this quantity should be grad(u)/x = (2, 0) at every point of the domain. I attach a picture of the x-component of this quantity (should be 2), as it can be visualized in ParaView:

Clearly, there is a problem in the axis (x=0) due to the division by x. Is there a way to avoid this problem in FEniCS and obtain a good visualization? I know this problem may be attached by increasing the finite element order to 2, but this causes other issues in my real problem, so it would be preferable to avoid this procedure.

Many thanks in advance.