Hello everyone,

I want to obtain plots for a function (**phi_exact**) along the radial and circumferential directions (i.e. along the radius and along the angle theta), in Fenics 2019.1.0.

Following is my code:

```
import matplotlib.pyplot as plt
from dolfin import *
import logging
import numpy as np
mesh = Mesh("../meshes/annulus_refined.xml")
s = FiniteElement("P", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh,s)
cs = interpolate(Expression("cos(atan2(x[1],x[0]))",degree=3),W)
r = interpolate(Expression("sqrt(x[0]*x[0] + x[1]*x[1])",degree = 3),W)
# Function to be plotted along the radial and circumferential directions
phi_exact = project(-(r - (1/r))*cs,W)
```

Following is how my mesh looks like:

Is there any function in Fenics that could provide the radial and circumferential points/ vertices over which the function **phi_exact** could be plotted? Any leads would be greatly appreciated.

The most straightforward approach is probably just to directly evaluate `phi_exact`

at the points you want, e.g., something like the following example:

```
import matplotlib.pyplot as plt
from dolfin import *
import logging
import numpy as np
# Replaced custom mesh with unit square for testing
#mesh = Mesh("../meshes/annulus_refined.xml")
N = 32
mesh = UnitSquareMesh(N,N)
s = FiniteElement("P", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh,s)
cs = interpolate(Expression("cos(atan2(x[1],x[0]))",degree=3),W)
r = interpolate(Expression("sqrt(x[0]*x[0] + x[1]*x[1])",degree = 3),W)
# Function to be plotted along the radial and circumferential directions
phi_exact = project(-(r - (1/r))*cs,W)
# Plotting "by hand", by sampling `phi_exact` at custom points:
import math
max_theta = 0.5*math.pi
Nstep = 128
theta = np.zeros(Nstep+1)
phi = np.zeros(Nstep+1)
rad = 0.5
for i in range(0,Nstep+1):
theta[i] = i*max_theta/Nstep
x0 = rad*math.cos(theta[i])
x1 = rad*math.sin(theta[i])
# The `__call__` method of `Function` is overloaded to evaluate at points.
phi[i] = phi_exact(np.array([x0,x1]))
# Could also write out arrays `theta` and `phi` in .csv format, to
# plot with other tools.
plt.plot(theta,phi)
plt.show()
```

There may also be a way to do something similar by writing `phi_exact`

to a file and applying filters in ParaView.

2 Likes

@kamensky , I am thankful to you for the appropriate response and I apologize for the late update.

I am trying to have a line plot of a variable, which is being read through the HDF5 module. The file read by the HDF5 module is the last timestamp of my numerical solution.

```
from dolfin import *
import matplotlib.pyplot as plt
mesh = Mesh('../../res/mesh_generation/annulus_refined.xml')
V = VectorElement("Lagrange", mesh.ufl_cell(), 2)
Q = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh,MixedElement([V,Q,Q,Q,Q]))
hdf5 = HDF5File(mesh.mpi_comm(), './state.h5', 'r')
qTilde = Function(W)
hdf5.read(qTilde, 'var')
u,p,c,q,phi = split(qTilde)
```

I wish to plot **phi** as a line plot (which has been obtained by **splitting** **qTilde**), to compare it with **phi_exact** in radial and circumferential directions. As far as I know **phi** is a ufl variable and **phi_exact** is doflin expression, so does the ufl variable **phi** needs to be converted to dolfin format?

Any leads would be greatly appreciated.

Legacy FEniCS has two similar-looking but different ways to â€śsplitâ€ť `Function`

s in mixed spaces. See the following example for a way to get `phi`

as a `Function`

:

```
from dolfin import *
import matplotlib.pyplot as plt
#mesh = Mesh('../../res/mesh_generation/annulus_refined.xml')
mesh = UnitSquareMesh(8,8)
V = VectorElement("Lagrange", mesh.ufl_cell(), 2)
Q = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh,MixedElement([V,Q,Q,Q,Q]))
#hdf5 = HDF5File(mesh.mpi_comm(), './state.h5', 'r')
qTilde = Function(W)
#hdf5.read(qTilde, 'var')
qTilde.assign(Constant((0,1,2,3,4,5)))
# Long-standing ambiguity in the Legacy FEniCS API:
#
# https://bitbucket.org/fenics-project/dolfin/issues/194/split-u-versus-usplit
#
# Difference between "split" as free function and method:
u,p,c,q,phi = split(qTilde)
print(type(phi))
_,_,_,_,phi_func = qTilde.split()
print(type(phi_func))
# Can evaluate `phi_func` at a spatial point:
from numpy import array
print(phi_func(array([0.3,0.3])))
```

1 Like

I am extremely thankful to you for pointing that out!!!