You can use the following with P1 elements. If your 1D space can be mapped from the unit interval, the analytical solutions $(\lambda, \xi)$ of $\int_D C(x, x^\prime) \xi(x) \; \mathrm{d}x = \lambda \xi(x^\prime)$ are known, I believe.
from dolfin import *
import numpy as np
mesh = IntervalMesh(64, -1., 1.)
V = FunctionSpace(mesh, 'CG', 1)
v = TestFunction(V)
u = TrialFunction(V)
# Covariance kernel function
k = Expression('exp(-std::abs(x_i - x[0])/lam)', x_i=0.0, lam=1)
# Covariance discrete matrix
C = np.zeros((mesh.num_vertices(), mesh.num_vertices()))
# Evaluate \int C(x, x') \xi(x') dx'
coords = mesh.coordinates()
dof2v = dof_to_vertex_map(V)
for j in range(len(coords)):
k.x_i = Vertex(mesh, dof2v[j]).midpoint()[0]
C[j, :] = assemble(k*v*dx)
You may want to use this to compute the Mercer basis functions. This can be done by doing a similar Galerkin projection of $\int_D C(x, x^\prime) \xi(x) \; \mathrm{d}x$ onto the 1D interval and solving the corresponding eigenvalue problem. Keep in mind that the system will be dense.
