Evaluating operator based on pairs of nodes

General
1

I have to calculate the eigenvalues/-functions of an operator that is based on the L2-distance between two points. My idea was to evaluate the operator at each pair of nodes to obtain discrete eigenvectors and map these to the function space afterwards. However, my implementation does not yield a symmetric matrix, which should be the case as the elements are only based on the distance between two points, so the ordering of these points should not matter.


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

mesh = RectangleMesh(Point(0,0), Point(1, 1), 10, 10)
V = FunctionSpace(mesh, 'Lagrange', 1)
x = SpatialCoordinate(mesh)
v = TestFunction(V)
x_i = Constant(0)
y_i = Constant(0)
k = sqrt((x[0]-x_i)**2+(x[1]-y_i)**2)
k_on_mesh = k*v*dx

C = np.zeros((mesh.num_vertices(), mesh.num_vertices()))
coords = mesh.coordinates()
dof2v = dof_to_vertex_map(V)
for i in range(len(coords)):
    mp = Vertex(mesh, dof2v[i]).midpoint()
    x_i.assign(mp[0])
    y_i.assign(mp[1])
    C[i, :] = assemble(k_on_mesh)

vals, vecs = np.linalg.eig(C)

plt.imshow(C - np.transpose(C))
plt.colorbar()
plt.show()

Specifically the expression k * v * dx bothers me, as, when I plot the rows of C, these do not correspond with values of k taking one point fixed and varying the other across the mesh nodes (which they should). How can I fix this?