Plotting along the radial and circumferential directions

mesh
1

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:

annulus
annulus1920Ă—920 77 KB

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.

2

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
3

@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.

4

Legacy FEniCS has two similar-looking but different ways to “split” Functions 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
5

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