# Problem when visualizing near the symmetry axis

**URL:** <https://fenicsproject.discourse.group/t/problem-when-visualizing-near-the-symmetry-axis/5336>\
**Category:** mesh\
**Created:** [March 17, 2021, 4:20pm UTC](https://fenicsproject.discourse.group/t/problem-when-visualizing-near-the-symmetry-axis/5336 "2021-03-17T16:20:19Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![ims](https://avatars.discourse-cdn.com/v4/letter/i/eb9ed0/32.png) [@ims](https://fenicsproject.discourse.group/u/ims)\
**Post date:** [March 17, 2021, 4:20pm UTC](https://fenicsproject.discourse.group/t/problem-when-visualizing-near-the-symmetry-axis/5336/1 "2021-03-17T16:20:19Z")

</div>

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:

 ![axisym](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/d/d6b5c27b806567cc63382b16a4ecd4e7179a8b5d.png)

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.

---

<div class="post-metadata">

**Author:** ![hernan](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/hernan/32/67_2.png) [@hernan](https://fenicsproject.discourse.group/u/hernan)\
**Post date:** [March 17, 2021, 6:27pm UTC](https://fenicsproject.discourse.group/t/problem-when-visualizing-near-the-symmetry-axis/5336/2 "2021-03-17T18:27:55Z")

</div>

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

```python
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:

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

```

---

<div class="post-metadata">

**Author:** ![ims](https://avatars.discourse-cdn.com/v4/letter/i/eb9ed0/32.png) [@ims](https://fenicsproject.discourse.group/u/ims)\
**Post date:** [March 18, 2021, 8:28am UTC](https://fenicsproject.discourse.group/t/problem-when-visualizing-near-the-symmetry-axis/5336/3 "2021-03-18T08:28:51Z")

</div>

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!”.

---

<div class="post-metadata">

**Author:** ![ims](https://avatars.discourse-cdn.com/v4/letter/i/eb9ed0/32.png) [@ims](https://fenicsproject.discourse.group/u/ims)\
**Post date:** [March 22, 2021, 3:20pm UTC](https://fenicsproject.discourse.group/t/problem-when-visualizing-near-the-symmetry-axis/5336/4 "2021-03-22T15:20:45Z")

</div>

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