Oscillations in stress field due to plastic strain (CG)

I want to solve a coupled problem: I solve one equation (equilibrium, $div(\sigma)=0$) and use the solution to calculate a ($\sigma$-dependent) velocity field. This velocity field is used to solve a transport equation. The solution of the transport equation is then used to update the plastic strain and with it the stress field - and then so on.

I've come across a problem. It seems, that the choice of function spaces is crucial. And i may be doing domething wrong:

For the equilibrium

$$\sigma = C : (\epsilon - \epsilon_p)$$

is solved. With $\epsilon$ beeing

$$\epsilon_{ij} = u_{i,j} + u_{j,i}$$

and the displacement $u$ being in CG2, $\epsilon$ has to be in DG1. $\epsilon_p$ is coming from the transport equation.

The transport equation

$$ \frac{\partial\rho}{\partial t} + \frac{\partial(v \rho)}{\partial x} = -sinkterm $$

is solved with FCT (flux corrected transport). The solution $\rho$ has to be in CG1. Since $\rho$ is used to calculate $\epsilon_{p}$, $\epsilon_{p}$ also is in CG1.

It seems that the plastic strain in CG1 is causing oscillations in the stress field. Since my actual code is too long, i created two minimal examples, where in one example the plastic strain is in CG1 and in the other one it is in DG1. You can see the results below:

comparison

The solution seems to be wrong for plastic strain in CG1. In areas with high gradients, the solution is looking bad. I am a little bit lost here.

Is it possible to somehow use the plastic strain solution in CG1 and get a smooth stress field? Maybe transform from CG1 to DG1? (project($\epsilon_p$,DG1) doesn't work)
Is it possible to combine proper element types to get the desired results?
Will using a mixed FE formulation help us? Or are there any other ways?

Posting both examples ($\epsilon_p$ in CG1 and DG1) is too long, and a mesh is needed. I post the code of the CG1 example and here are download links for both examples with mesh:

CG1:
https://www.dropbox.com/s/78k9ev0j89idvbd/CG1.zip?dl=0

DG1:
https://www.dropbox.com/s/qm6byracq9cwh5w/DG1.zip?dl=0

Working example for case $\epsilon_p$ in CG1 (with oscillations)

# Import all packages from fenics into cache
from fenics import *
import os


abspath = os.path.abspath(__file__)
dname = os.path.dirname(abspath)
os.chdir(dname)


mesh = Mesh("2D_Mesh.xml")
celllist  = MeshFunction('size_t',mesh,'2D_Mesh_physical_region.xml');
#facetlist = MeshFunction('size_t',mesh,'2D_Mesh_facet_region.xml');


domains = CellFunction('size_t', mesh)
domains.set_all(0)
domains.array()[celllist.array()==31] = 1


dA = Measure("ds", domain=mesh, subdomain_data=facets)
dV = Measure("dx", domain=mesh, subdomain_data=domains)

N = FacetNormal(mesh)

#----------------------------------------------------------------------------#

ScalarSpaceCG1 = FunctionSpace(mesh,'CG',1)
VectorSpaceCG2 = VectorFunctionSpace(mesh,'CG',2)
TensorSpaceDG1 = TensorFunctionSpace(mesh,'DG',1)
TensorSpaceCG1 = TensorFunctionSpace(mesh,'CG',1)

#----------------------------------------------------------------------------#
nu = 0.3
E = 21.0
G = E/(2.0*(1.0+nu))
lamda = 2.0*G*nu / (1.0-2.0*nu)
mu = G
#----------------------------------------------------------------------------#

gamma = Function(ScalarSpaceCG1)
dofmap = ScalarSpaceCG1.dofmap()
dofs = []

for cell in cells(mesh): # compute dofs in the domains
        if domains[cell] == 1:
                dofs.extend(dofmap.cell_dofs(cell.index()))
        else:
                pass

dofs = list(set(dofs))


gamma_Expression = Expression("0.05", element = ScalarSpaceCG1.ufl_element())
gamma_Expression = interpolate(gamma_Expression, ScalarSpaceCG1)

gamma.vector()[dofs] = gamma_Expression.vector()[dofs]

#----------------------------------------------------------------------------#

n_1 = Constant((1.0,0.0))
m_1 = Constant((0.0,-1.0))

epsilon_plastic = -gamma*(  0.5*(outer(n_1,m_1) + outer(m_1,n_1))  )
epsilon_plastic = project( epsilon_plastic, TensorSpaceCG1 )

#----------------------------------------------------------------------------#

def epsilon(u):
    return as_tensor(0.5*(u[i].dx(j)+u[j].dx(i)), (i,j))

def epsilon_elastic(epsilon, epsilon_plastic):
    return as_tensor( epsilon[i,j] - epsilon_plastic[i,j], (i,j))

def sigma(epsilon_elastic):
    return as_tensor(lamda*epsilon_elastic[k,k]*delta[i,j] + 2.0*mu*epsilon_elastic[i,j], (i,j))

#----------------------------------------------------------------------------#

delta = Identity(2)

u = TrialFunction(VectorSpaceCG2)
del_u = TestFunction(VectorSpaceCG2)

a = sigma(epsilon_elastic(epsilon(u), epsilon_plastic))[j,i]*del_u[i].dx(j)*dV

Sigma_Global = as_matrix([[0.0, 0.0], [0.0, 0.0]]) # No Load
L = Sigma_Global[i,j]*N[j]*del_u[i]*dA

eq = a - L

a = lhs(eq)
L = rhs(eq)

Displacement = Function(VectorSpaceCG2)

solve (a==L, Displacement, bcs=[])

#----------------------------------------------------------------------------#

file_stress = File('stress.pvd')
file_stress << project( sigma(epsilon_elastic(epsilon(Displacement), epsilon_plastic)), TensorSpaceDG1 )

file_epsilon_plastic = File('epsilon_plastic.pvd')
file_epsilon_plastic << epsilon_plastic

Thanks a lot.