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
return x < DOLFIN_EPS or x > 1.0 - DOLFIN_EPS or x < DOLFIN_EPS or x > 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 = - fv*dx
bc = DirichletBC(V, Expression(‘pow(x,2)’, degree=2), boundary)
u = Function(V)
solve(a == L, u, bc)
Save quantity J to file in XDMF format
varx = Expression(‘x’, 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.