Problem when visualizing near the symmetry axis

mesh
1

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 = - fv*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:

axisym
axisym797×492 22 KB

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.

2

You could define varx through the ufl function SpatialCoordinate as follows:

x = SpatialCoordinate(mesh)
varx = x[0]

Using this approach, varx is only evaluated at the quarature points of each element, which avoids the division by zero. Additionally, you can pass a finite element to your Expression to avoid the evaluation at nodal dofs:

Q=FiniteElement('Quadrature',triangle,degree=quad_degree,quad_scheme='default')
varx = Expression('x[0]', element=Q)
3

Thanks for your answer. But it actually doesn’t solve my problem. If x is set to SpatialCoordinate(mesh) and then varx = x[0], the visualization problem in the axis remains the same. In addition, if pass the finite element Q to the Expression (setting quad_degree = 3, for instance), I get an “AssertionError: Mismatch of quadrature points!”.

4

Does anyone know how to solve this problem? Any help is appreciated.