Wrong order of plotted node values in matplotlib even after using dofmap().dofs()

1

Hi,
I am trying to save the results of a 1D problem with mixed elements. After using dofmap().dofs() the order is still wrong when plotting. The code is shown below:

import matplotlib.pyplot as plt
from fenics import *
import numpy as np

rho_p = 3000.0                
d_p   = 100.0e-6;              
m_p   = rho_p*pi*d_p**3/6.0   
l        = 30.0e-3            
numElems = 100                
abc = 2.0      
dfg = abc/60000.0
n   = dfg/m_p
r_0 = 1.0e-3
h_d      = 17.0e-3                           
U        = 0.017 
Dz_F     = 10.0e-3;
PI     = 3.1415926535
theta  = 10.95*PI/180.0

# Define function space for system
deg = 2
mesh = IntervalMesh(numElems, 0.0,l/2.0)
P1 = FiniteElement('Lagrange', interval, deg)

element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)

R = Expression("r_0 + (h_d-Dz_F/2.0)*tan(theta)", degree=deg, r_0=r_0, h_d=h_d, Dz_F=Dz_F, theta=theta)
u_ic='tan(theta)*U*x[0]/R'
c_ic='2.0*n/U/pi/pow(R,2)*exp( -2.0*pow(x[0],2)/pow(R,2) )'
w_0 = Expression((u_ic, c_ic), degree=deg, theta=theta, U=U, R=R, n=n, pi=PI)


#Write values in text file
w_p = project(w_0, V)

dof_coordinates = V.tabulate_dof_coordinates()
n = V.dim()
d = mesh.geometry().dim()
dof_coordinates.resize((n, d))

u_dofs = V.sub(0).dofmap().dofs()
c_dofs = V.sub(1).dofmap().dofs()

plt.plot(dof_coordinates[c_dofs], w_p.vector()[c_dofs])
plt.show()

Thank you for your time,
Alex

2

You can try

plot(w_p.sub(1))

and also

print(w_p.sub(1).compute_vertex_values())
3

Hi volkerk thank you for the reply.

w_p.sub(1).compute_vertex_values() only works when the degree of the elements is 1 and the values are displayed right to left.

Using w_p.sub(1) gives the following error:
NotImplementedError: Cannot take length of non-vector expression

4 
w_p.sub(1)compute_vertex_values()

works with arbitrary order elements. It calculates the values of the function at the mesh vertices.
https://fenicsproject.org/docs/dolfin/latest/python/_autogenerated/dolfin.cpp.function.html#dolfin.cpp.function.GenericFunction

I do not get the error for w_p.sub(1) when using your code:

import matplotlib.pyplot as plt
from fenics import *
import numpy as np

rho_p = 3000.0                
d_p   = 100.0e-6;              
m_p   = rho_p*pi*d_p**3/6.0   
l        = 30.0e-3            
numElems = 100                
abc = 2.0      
dfg = abc/60000.0
n   = dfg/m_p
r_0 = 1.0e-3
h_d      = 17.0e-3                           
U        = 0.017 
Dz_F     = 10.0e-3;
PI     = 3.1415926535
theta  = 10.95*PI/180.0

# Define function space for system
deg = 2
mesh = IntervalMesh(numElems, 0.0,l/2.0)
P1 = FiniteElement('Lagrange', interval, deg)

element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)

R = Expression("r_0 + (h_d-Dz_F/2.0)*tan(theta)", degree=deg, r_0=r_0, h_d=h_d, Dz_F=Dz_F, theta=theta)
u_ic='tan(theta)*U*x[0]/R'
c_ic='2.0*n/U/pi/pow(R,2)*exp( -2.0*pow(x[0],2)/pow(R,2) )'
w_0 = Expression((u_ic, c_ic), degree=deg, theta=theta, U=U, R=R, n=n, pi=PI)


#Write values in text file
w_p = project(w_0, V)

dof_coordinates = V.tabulate_dof_coordinates()
n = V.dim()
d = mesh.geometry().dim()
dof_coordinates.resize((n, d))

u_dofs = V.sub(0).dofmap().dofs()
c_dofs = V.sub(1).dofmap().dofs()

print(w_p.sub(1).compute_vertex_values())

#plt.plot(dof_coordinates[c_dofs], w_p.vector()[c_dofs])
plot(w_p.sub(1))
plt.show()
5

Hi thank you the prompt reply.

Its strange I’m getting the same error in my installation even when using plot(w_p.sub(1))

Also even if plot(w_p.sub(1)) were working I’m looking to write this into a numpy array. I’m doing the plotting first to check the ordering

I think the issue with w_p.sub(1).compute_vertex_values() is that it prints the vertex values of the elements and not the inside nodes, if you try
plt.plot(dof_coordinates[c_dofs], w_p.sub(1).compute_node_values())
With deg=2 it doesn’t work because they have a different size

6

Yes, as i wrote, the .compute_vertex_values() function computes the values at the mesh vertices, not the dof nodes. Please describe your desired result precisely.

If you just want to save the function, why not use

function_file = XDMFFile(MPI.comm_world, './solution.xdmf')
function_file.write_checkpoint(w_p.sub(1), 'w_p.sub(1)')
function_file.close()

You can then view the solution e.g. in paraview.

I don’t know why w_p.sub(1) does not work for you, but you could try

w_0, w_1 = w_p.split(True)

and then work with w_1.

7

Hi, using the w_p.split(True) method compiles without error but it still gives the same error in the ordering as when using the original code.

I would like to export the results to MATLAB that’s why I would like to save a numpy array of the results in a text file.

8

The dofs are not ordered in increasing direction. You need to sort them and the output. For instance, consider:

import numpy
asc_order = numpy.argsort(w_p.vector()[c_dofs])
permuted_dofs = [dof_coordinates[c_dofs][i] for i in asc_order]
permuted_solution = [w_p.vector()[c_dofs][i] for i in asc_order]
plt.plot(permuted_dofs, permuted_solution, "-go")
9

Hi dokken, thank you for the reply.

Since the answer might not be always decreasing I used the same method but with the dofs. I had to take one of the size 1 dimensions to get it working:

dofs      = np.squeeze(dof_coordinates[c_dofs])
asc_order = np.argsort(dofs)
solution  = w_p.vector()[c_dofs]

plt.plot(dofs[asc_order], solution[asc_order])
plt.show()

I’m still not 100% sure why the plot shows the wrong ordering, even if the dofs are ordered wrong the corresponding solution values are not ordered in the same way woudn’t that give you the right plot anyways?

10

I agree that you should use argsort on the x-coordinate (a bit too fast from my side). To elaborate,
if you plot it with the option plt.plot(dof_coordinates[c_dofs], w_p.vector()[c_dofs], "ro") you will observe that the dofs are at the correct places. It is just that matplotlib by default uses lines to Connect points linearly (equivalent to “-ro”), Which makes the plot look funky.