How to change a function from cartesian to polar?

Hi! In the below MWE, nh is a function in the cartesian frame on the boundary of a circle. I need to plot nh as a function of \theta. I tried to find some posts in the community page but didn’t come across any. I tried a transformation matrix but it didn’t work. Can anyone please help me with this?

MWE:

from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt

from dolfin import *
import meshio
from dolfin import Mesh, XDMFFile, File, MeshValueCollection, cpp

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}

import numpy as np
import meshio
from mshr import Circle, generate_mesh
from dolfin import Mesh, File, MeshFunction, Point, BoundaryMesh, SubDomain, plot, File
import matplotlib.pyplot as plt
from dolfin import *
from mshr import *

class boundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary

class disk(SubDomain):
    def inside(self, x, on_boundary):
        return True


mesh = generate_mesh(Circle(Point(0, 0), 3), 50)
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
boundary().mark(boundary_markers, 2)
disk().mark(surface_markers, 1)

ds = Measure('ds', subdomain_data=boundary_markers)
dx = Measure('dx', subdomain_data=surface_markers)
n = FacetNormal(mesh)

W1 = FunctionSpace(mesh, "Lagrange", 1)
n = FacetNormal(mesh)
G , mu = 1, 0.1
u_D=Constant(0.0)

bc = DirichletBC(W1, u_D, boundary_markers, 2)

# Define variational problem
u = TrialFunction(W1)
v = TestFunction(W1)
f = Constant(-G/mu) # f=-G/mu
a = dot(grad(u), grad(v))*dx
L = -f*v*dx

# Compute solution
u = Function(W1)
solve(a == L, u, bc)

# Calculating shear stress S
Ux=0.0
Uy=0.0
Uz=u
U = as_vector((Uy, u))
N = as_matrix([[0.0,0.5],
                 [0.0,0.0]])
Sn = mu/2 * dot(grad(U) + grad(U).T, (N.T)*n)
print(assemble(Sn[1]*ds))
S=Sn[1]
u = TrialFunction(W1)
v = TestFunction(W1)
a = inner(u,v)*ds
l = inner(S, v)*ds
A = assemble(a, keep_diagonal=True)
L = assemble(l)
A.ident_zeros()
nh = Function(W1)
solve(A, nh.vector(), L)
File("c2ps.pvd") << nh  # This (nh) is the function I need to convert to a polar frame.

x = mesh.coordinates()[:,0]  # Defining x as the x-coordinates of the mesh
y = mesh.coordinates()[:,1]  # Defining y as the y-coordinates of the mesh
def r(x, y):                 # Defining r as a function of x and y
  return [sqrt(x*x + y*y)]   # Returning the formula for r
def theta(x, y):             # Defining theta as a function of x and y
  return [atan(y,x)]         # Returning the formula for theta

How to change the x and y in nh to nh(\theta)?

Any suggestion would be great.

One approach would be to just iterate over a range of \theta values, compute the corresponding (x,y) points, and evaluate the function nh at those points:

# Allow extrapolation of nh in case a point falls slightly
# outside of the mesh.
nh.set_allow_extrapolation(True)

# Iterate over a bunch of points around the circular
# boundary and evaluate nh at each one.
import math
from numpy import array, zeros
N_theta = 100
dtheta = 2*math.pi/N_theta
r = 3.0
nh_theta = zeros(N_theta)
theta = zeros(N_theta)
for i in range(0,N_theta):
    theta[i] = i*dtheta
    x = np.array([r*math.cos(theta[i]), r*math.sin(theta[i])])
    nh_theta[i] = nh(x)

# Write two-column file that can be plotted easily with Gnuplot.
nh_theta_file = open("nh_theta","w")
nh_theta_str = ""
for i in range(0,N_theta):
    nh_theta_str += str(theta[i])+" "+str(nh_theta[i])+"\n"
nh_theta_file.write(nh_theta_str)
nh_theta_file.close()

The resulting file named nh_theta could be plotted with Gnuplot using the command
plot "nh_theta" with lines.

Hi kamensky, thanks for your reply. I tried the suggested code but it returned an error. As I am not sure what “Gnuplot” means, I tried the conventional plotting command. Below is the full script followed by the error:

Script:

from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt
import numpy as np
import sympy as sym
from sympy import symbols
from sympy import atan2,Abs,cos,sin

from dolfin import *
import meshio
from dolfin import Mesh, XDMFFile, File, MeshValueCollection, cpp



# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}


import numpy as np
import meshio
from mshr import Circle, generate_mesh
from dolfin import Mesh, File, MeshFunction, Point, BoundaryMesh, SubDomain, plot, File
import matplotlib.pyplot as plt
from dolfin import *
from mshr import *


class boundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary


class disk(SubDomain):
    def inside(self, x, on_boundary):
        return True


mesh = generate_mesh(Circle(Point(0, 0), 3), 50)
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
boundary().mark(boundary_markers, 2)
disk().mark(surface_markers, 1)


ds = Measure('ds', subdomain_data=boundary_markers)
dx = Measure('dx', subdomain_data=surface_markers)
n = FacetNormal(mesh)

W1 = FunctionSpace(mesh, "Lagrange", 1)
n = FacetNormal(mesh)
G , mu = 1, 0.1
u_D=Constant(0.0)

bc = DirichletBC(W1, u_D, boundary_markers, 2)

# Define variational problem
u = TrialFunction(W1)
v = TestFunction(W1)
f = Constant(-G/mu) # f=-G/mu
a = dot(grad(u), grad(v))*dx
L = -f*v*dx

# Compute solution
u = Function(W1)
solve(a == L, u, bc)

file = File("fsp1.pvd");
file << u;

# Calculating shear stress S
Ux=0.0
Uy=0.0
Uz=u
U = as_vector((Uy, u))
N = as_matrix([[0.0,0.5],
                 [0.0,0.0]])
Sn = mu/2 * dot(grad(U) + grad(U).T, (N.T)*-n)
print(assemble(Sn[0]*ds), assemble(Sn[1]*ds))
S=Sn[1]
u = TrialFunction(W1)
v = TestFunction(W1)
a = inner(u,v)*ds
l = inner(S, v)*ds
A = assemble(a, keep_diagonal=True)
L = assemble(l)

A.ident_zeros()
nh = Function(W1)

solve(A, nh.vector(), L)

File("fsp2.pvd") << nh

# Allow extrapolation of nh in case a point falls slightly
# outside of the mesh.
nh.set_allow_extrapolation(True)

# Iterate over a bunch of points around the circular
# boundary and evaluate nh at each one.
import math
from numpy import array, zeros
N_theta = 100
dtheta = 2*math.pi/N_theta
r = 3.0
nh_theta = zeros(N_theta)
theta = zeros(N_theta)
for i in range(0,N_theta):
    theta[i] = i*dtheta
    x = np.array([r*math.cos(theta[i]), r*math.sin(theta[i])])
    nh_theta[i] = nh(x)

# Write two-column file that can be plotted easily with Gnuplot.
nh_theta_file = open("nh_theta","w")
nh_theta_str = ""
for i in range(0,N_theta):
    nh_theta_str += str(theta[i])+" "+str(nh_theta[i])+"\n"
nh_theta_file.write(nh_theta_str)
nh_theta_file.close()
#The resulting file named nh_theta could be plotted with Gnuplot using the command
#plot "nh_theta" with lines.
plot(theta)
plot(nh_theta)
plt.show()

Error message (this is a partial error message, the actual message was quite large):

RuntimeError: Don't know how to plot given object:
  [0.         0.06283185 0.12566371 0.18849556 0.25132741 0.31415927
 0.37699112 0.43982297 0.50265482 0.56548668 0.62831853 0.69115038
 0.75398224 0.81681409 0.87964594 0.9424778  1.00530965 1.0681415
 1.13097336 1.19380521 1.25663706 1.31946891 1.38230077 1.44513262
 1.50796447 1.57079633 1.63362818 1.69646003 1.75929189 1.82212374
 1.88495559 1.94778745 2.0106193  2.07345115 2.136283   2.19911486
 2.26194671 2.32477856 2.38761042 2.45044227 2.51327412 2.57610598
 2.63893783 2.70176968 2.76460154 2.82743339 2.89026524 2.95309709
 3.01592895 3.0787608  3.14159265 3.20442451 3.26725636 3.33008821
 3.39292007 3.45575192 3.51858377 3.58141563 3.64424748 3.70707933
 3.76991118 3.83274304 3.89557489 3.95840674 4.0212386  4.08407045
 4.1469023  4.20973416 4.27256601 4.33539786 4.39822972 4.46106157
 4.52389342 4.58672527 4.64955713 4.71238898 4.77522083 4.83805269
 4.90088454 4.96371639 5.02654825 5.0893801  5.15221195 5.2150438
 5.27787566 5.34070751 5.40353936 5.46637122 5.52920307 5.59203492
 5.65486678 5.71769863 5.78053048 5.84336234 5.90619419 5.96902604
 6.03185789 6.09468975 6.1575216  6.22035345]
and projection failed:
  'numpy.ndarray' object has no attribute 'ufl_domain'

The arrays theta and nh_theta cannot be plotted using the plot function from DOLFIN. However, you can use them as lists of x and y coordinates for any 2D plotting functionality of your choice. I happen to be most familiar with Gnuplot, which is a separate program you would run afterward, using the outputted file as input. You could also plot an nh_theta vs. theta graph directly from the Python script using the matplotlib.pyplot.plot function:

plt.plot(theta,nh_theta)

Note that this is different from DOLFIN’s plot function.

Thanks, the program executed perfectly. But there is an issue.

My actual domain is not a circle. I created the circle for creating an MWE for the community. My actual domain is a arbitrary closed contour. How to modify the code to deal with any arbitrary closed domain? For the case of an arbitrary closed domain, r is not 3. So the part of the code r = 3.0 will not work.

For an arbitrary closed domain, I tried to change r = 3.0 to r = ds/dtheta, which produced the error

r = ds/dtheta
TypeError: unsupported operand type(s) for /: 'Measure' and 'float'

How to modify my r in such cases?

Again thanks for the help @kamensky . Much appreciated.

Hi @kamensky , I came up with the following code for calculating r for an arbitrary close domain:

import math
from numpy import array, zeros
N_theta = 100
dtheta = 2*math.pi/N_theta
#r = 3.0
nh_theta = zeros(N_theta)
theta = zeros(N_theta)
r = zeros(N_theta)
for i in range(0,N_theta):
    theta[i] = i*dtheta
    r[i] = ds/theta[i]
    x = np.array([r[i]*math.cos(theta[i]), r[i]*math.sin(theta[i])])
    nh_theta[i] = nh(x)

Now the error reads:

 r[i] = ds/theta[i]
TypeError: unsupported operand type(s) for /: 'Measure' and 'float'

How to fix this? Any suggestion will be extremely helpful.