I need to interpolate the plastic strain and vonmises stress at specific coordinate using shape function because my element size is not consistent while having different volume fraction of inclusion in the matrix.

If I understand you correctly, you could try a stress recovery algorithm to obtain the VM stress and plastic strain at the elements nodal locations. Essentially a simple form of these algorithms could be the inverse isoparametric mapping. Then you can use the local interpolation functions to compute the VM stress and plastic strain at any point in the element. There are a few things to worry about in trying to do this, namely getting unphysical (negative) values for the plastic strain, but you could make it work.

Thank  you

If you can import your data (plastic strain, VM stress, nodal coordinates ) into matlab (or a similar program). Then you can map the response (i.e. plastic strain or VM stress) as a function of x and y coordinates. Once you do that you can easily obtain an interpolated value at any location you want.