This is somehow hidden in FEniCS. It will depend on how you defined your integration measures, in the bilinear form for instance. Gauss points are not attached to a mesh or an interpolation space but to an integration measure such as dx. Without specifying anything, FEniCS will automatically deduce the required quadrature degree for the form you have defined.
If you have a “standard” bilinear form such as a = dot(grad(u), grad(v))*dx and u, v are of degree 2 let’s say, then it will be enough to compute the bilinear form with a degree 2 quadrature (i.e. 3 Gauss points on a triangle).
You can pass yourself the required degree d by using dx(metadata={"quadrature_degree": d}).
Why do you want to display the Gauss points exactly ?
If you want to see where are the Gauss points of degree d=2 for instance on a 2D mesh you can use a Quadrature function space with the corresponding degree such as:
from dolfin import *
import matplotlib.pyplot as plt
mesh = UnitSquareMesh(1,1, "crossed")
Ve = FiniteElement("Quadrature", mesh.ufl_cell(), degree=2, quad_scheme='default')
V = FunctionSpace(mesh, Ve)
plot(mesh)
dofs = V.tabulate_dof_coordinates()
for (i, dof) in enumerate(dofs):
plt.plot(*dof, "xk")
plt.annotate(str(i), xy=dof, xytext=dof+0.02)