Mesh.coordinates() vs. u_nodal_values.array() and using a numerical source function

Question

Hello,

Suppose I have calculated a solution of Poisson's equation on a certain mesh. To find u-values in an array, and mesh coordinates, I use:

u_array=u_nodal_values.array()
mesh_points=mesh.coordinates()

However, when I try to make a scatter plot of the solution u, like:

x0=mesh_points[:,0]
x1=mesh_points[:,1]
fig=pyplot.figure
ax=fig.add_subplot(1,1,1,projection='3D')
surf=ax.scatter3D(x0,x1,u_array)

I get all points randomly distributed over space, while for example

u_array2=np.array([u(Point(x,y)) for x,y in mesh_points])

works fine in a scatter plot. So why do the u_nodal_values.array() not correspond to the associated mesh.coordinates()?

This question is somewhat related to a second one. Suppose I have a numerical expression for each of my N nodal points, called d. This d is a N times 1 vector, corresponding correctly to the mesh.coordinates(), so I can plot this vector like

surf=ax.scatter3D(x0,x1,d)

Now, I want this vector d as the source function in my FEM system. So instead of an analytical source f=Expression('function of x[0] and x[1]'), I want this vector as a source for all my N nodal points. Is it possible to use this numerical expression as a source function in FEniCS?

Thanks for any help,

Adriaan

Answer 1

Hi, there is a difference between how coordinates of vertices in mesh are ordered and the order of degrees of freedom of functions from functions spaces over that mesh. In this call

surf=ax.scatter3D(x0,x1,u_array)

you are assuming that the order is the same, but as I said this does not hold. To get the order correct you should check out DofMap class. The following code can be used to scatter plot functions. Note that it is plotting values at all dofs not just those at vertices.

from dolfin import *                                                             
import numpy as np                                                               
from mpl_toolkits.mplot3d import Axes3D                                          
import matplotlib.pyplot as plt                                                  

mesh = CircleMesh(Point(0., 0.), 1, 0.1)                                         
V = FunctionSpace(mesh, "CG", 2)                                                 

n = V.dim()                                                                      
d = mesh.geometry().dim()                                                        

dof_coordinates = V.dofmap().tabulate_all_coordinates(mesh)                      
dof_coordinates.resize((n, d))                                                   
dof_x = dof_coordinates[:, 0]                                                    
dof_y = dof_coordinates[:, 1]                                                    

# use FEniCS to get some data to plot                                            
out = Expression("sin(pi*x[0])*cos(pi*x[1])")                                    
u = interpolate(out, V)                                                          

fig = plt.figure()                                                               
ax = fig.add_subplot(111, projection='3d')                                       
ax.scatter(dof_x, dof_y, u.vector().array(), c='b', marker='.')                  
plt.show()                                                                       

# now compute the data to be used by FEniCS                                      
z = np.exp(-(dof_x**2 + dof_y**2))                                               
u.vector()[:] = z                                                                

plot(u, interactive=True)      

Comments

Thanks a lot! This works fine to plot the u.vector().array().
Can you also explain how to use the dofs to use the numerical source function (the second part of my question)? As far as I know, a source function f should be either a constant or an expression of x[0] and x[1]. Instead, I have no analytic expression, but only numerical 'source values' at all mesh vertices, stored in a vector. Do you know how to 'tell' FEniCS the source values at each node?

---

Hi, to use your vector of values in variational formulation you should transform it to a Function object. Here's an example

from dolfin import *
import numpy as np

# solve -laplace(u) = f in [0, 1]^2 with homog. Dirichlet boundary conditions
# and funny 'numerical' rhs

mesh = RectangleMesh(-1, -1, 1, 1, 10, 10, 'right')

# source term calculation
F = FunctionSpace(mesh, 'CG', 1)  # so that dofs are only in mesh vertices
n = F.dim()
d = mesh.geometry().dim()

f = Function(F)

# dofmap approach
F_dof_coordinates = F.dofmap().tabulate_all_coordinates(mesh)
F_dof_coordinates.resize((n, d))
F_dof_x = F_dof_coordinates[:, 0]
F_dof_y = F_dof_coordinates[:, 1]

f_values = np.zeros(n)  
f_values = abs(F_dof_x + F_dof_y)
f.vector()[:] = f_values
plot(f, title='Source term')

# vertex to dofmap approach
vertex_values = np.zeros(mesh.num_vertices())
for vertex in vertices(mesh):
  x = vertex.x(0)
  y = vertex.x(1)
  vertex_values[vertex.index()] = abs(x + y)

f.vector()[:] = vertex_values[dof_to_vertex_map(F)]
plot(f, title='Source term')

# variational problem definition
V = FunctionSpace(mesh, 'CG', 2)

u = TrialFunction(V)
v = TestFunction(V)

a = inner(grad(u), grad(v))*dx
L = f*v*dx
bc = DirichletBC(V, Constant(0.), DomainBoundary())

u = Function(V)
solve(a == L, u, bc)

plot(u, title='Solution')
interactive() 

---

This does the job, now I can use a vector as a source function. Thanks for the example!

---

Regarding the example above, I had the case of V being a VectorFunctionSpace.
Using

mesh = UnitSquareMesh (10,10)
V = VectorFunctionSpace(mesh, "Lagrange", 1)      
n = V.dim()                                #  is 11*11*2 = 242                                                           
d = mesh.geometry().dim()            #  is  2                                             
dof_coordinates = V.dofmap().tabulate_all_coordinates(mesh)   

dof_coordinates is an array in which the vertices of the mesh are repeated two times, something like:

[ 0  0    0   0       0   0.1   0   0.1    0  0.2    0   0.2  etc    ]

To use

dof_coordinates.resize((n, d))                                                   
dof_x = dof_coordinates[:, 0]                                                   
dof_y = dof_coordinates[:, 1]    
v_array = v.vector().array()  # v is in the same space V
for i in range(n):
print 'v(%8g,%8g) = %g ' % (dof_x[i] , dof_y[i] , v_array[i]  ) 

works well.
In other words, it seems that each vertex coordinate is repeated two times, one for the x- and one for the y-component of v_array.
Is this always true, i.e., would the short code here give always a consistent outcome?

dof_coordinates = V.dofmap().tabulate_all_coordinates(mesh) 
dof_coordinates.resize((n/2, d*2))    
dof_x = dof_coordinates[:, 0]        # equivalent to [:,2]                                           
dof_y = dof_coordinates[:, 1]        # equivalent to [:,3]                                               
v_array.resize(n/2,2)
for i in range(n/2):
      print 'v(%8g,%8g) = %g , %g' % (dof_x[i] , dof_y[i] , v_array[i][0] , v_array[i][1]   )  

R.

---

Hi, I think it is safer not to assume anything about the ordering and use dofmap of components of V

from dolfin import *

mesh = UnitSquareMesh(10, 10)

V = VectorFunctionSpace(mesh, 'CG', 1)

u = interpolate(Expression(('1 + x[0]', '2 + x[1]')), V)
u_v = u.vector()

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

x_dofs = V.sub(0).dofmap().dofs()
y_dofs = V.sub(1).dofmap().dofs()

for x_dof, y_dof in zip(x_dofs, y_dofs):
  print dof_coordinates[x_dof], dof_coordinates[y_dof], u_v[x_dof], u_v[y_dof] 

---

Ok.
I think you are right about not assuming anything.
Thanks for your valuable advice.
R.

---

As of FEniCS 2016.1.0, this line stopped working for me:

dof_coordinates = V.dofmap().tabulate_all_coordinates(mesh)

I replaced it with the code below, which seems to work:

dof_coordinates = V.tabulate_dof_coordinates()