However, there, the applied stress field is constant over the domain.

What is the best way to initialize â€śmy_constant_tensorâ€ť in each element cell / material point? The stress data is a numpy array of dimensions (#material points, stress dimension).

The interpolate function seems to do the job. But I am wondering about the dimension x.shape[1] in your example. Printing it gives x.shape[1]=600, with 200 triangle cells. How do I obtain/set the number of integration points in this case?

You are creating a CG-1 quantity, which has one degree of freedom per vertex.
If you want something that is constant per cell, use a DG 0 function space.