Solving mixed strain-displacement formulation with linear-quadratic interpolation

I have problems with stress oscillations and found a stabilized mixed formulation (displacement and strain). It uses linear-linear interpolation and thus doesn't fulfil the inf-sup-condition. But the stabilization accounts for this and this code seems to run.

Then i wanted to run the code with linear(strain)-quadratic(displacement) interpolations. I get errors, depending on the used solver. For example:

*** Error: Unable to successfully call PETSc function 'KSPSolve'.
*** Reason: PETSc error code is: 76 (Error in external library).
*** Where: This error was encountered inside /build/dolfin-4SStI2/dolfin-2016.2.0/dolfin/la/PETScLUSolver.cpp.
*** Process: 0
*** DOLFIN version: 2016.2.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------

I don't understand it, since the only thing i change is increasing the degree of interpolation function. (I want to avoid the stabilization by fulfilling the inf-sup-condition.)

Source of the mixed formulation: http://www.sciencedirect.com/science/article/pii/S0266352X14001876

Thanks a lot.

I made a simplified version as running example:

from fenics import *
import os
import numpy as np

abspath = os.path.abspath(__file__)
dname = os.path.dirname(abspath)
os.chdir(dname)
#----------------------------------------------------------------------------#
mesh = UnitSquareMesh(32, 32)

domains = CellFunction('size_t', mesh)
domains.set_all(0)

facets = FacetFunction('size_t', mesh)
facets.set_all(0)

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

Normale = FacetNormal(mesh)
#----------------------------------------------------------------------------#
nu = 0.45
G = 33 * 10**(-3) # [N/um^2]
lamda = 2.0*G*nu / (1.0-2.0*nu)
mu = G
lamda = ((2*mu)/(lamda + 2*mu))*lamda
#----------------------------------------------------------------------------#

ScalarSpaceCG1 = FunctionSpace(mesh,'CG',1)
TensorSpaceCG1 = TensorFunctionSpace(mesh,'CG',1)


###### Crucial Part #######
#--------------------------
TCG = TensorElement("Lagrange", mesh.ufl_cell(), 1)
#### Running version
VCG = VectorElement("Lagrange", mesh.ufl_cell(), 1)
### NOT running version
# VCG = VectorElement("Lagrange", mesh.ufl_cell(), 2)
######################

W = TCG * VCG
W = FunctionSpace(mesh, W)

(epsilon_h, u_h) = TrialFunctions(W)
(gamma_h, v_h) = TestFunctions(W)

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

delta = Identity(2)  

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

def epsilon_elastisch(epsilon_gesamt):
       return as_tensor(epsilon_gesamt[i,j], (i,j))

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

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


#----------------------------------------------------------------------------#
# Boundary conditions

def boundary_xleft(x):
    return -DOLFIN_EPS < x[0] < DOLFIN_EPS
bc1 = DirichletBC(W.sub(1).sub(0), Constant((0.0)), boundary_xleft)

def boundary_xright(x):
    return -DOLFIN_EPS +1 < x[0] < DOLFIN_EPS +1
bc2 = DirichletBC(W.sub(1).sub(0), Constant((0.1)), boundary_xright)

def boundary_xleft_corner(x):
    return ( (-DOLFIN_EPS < x[0] < DOLFIN_EPS) and (-DOLFIN_EPS < x[1] < DOLFIN_EPS))
bc3 = DirichletBC(W.sub(1), Constant((0.0,0.0)), boundary_xleft_corner,method="pointwise")


Sigma_Global_Feld = as_matrix([[0.0, 0.0], [0.0, 0.0]])  # [N/um^2]
#----------------------------------------------------------------------------#

Eq1 = -gamma_h[i,j]*C_mal_T2Stufe(epsilon_h-sym_grad(u_h))[i,j]*dV 

Eq2 = sym_grad(v_h)[i,j]*C_mal_T2Stufe( epsilon_h )[i,j]*dV    - Sigma_Global_Feld[i,j]*Normale[j]*v_h[i]*dA

Eq = Eq1 + Eq2

a = lhs(Eq)
L = rhs(Eq)

w = Function(W)

solve (a==L, w, bcs=[bc1,bc2,bc3], solver_parameters={"linear_solver": "mumps"}, form_compiler_parameters={"optimize": True,"cpp_optimize":True})


# Post process
#-------------
(epsilon_h, u_h) = w.split()
#----------------------------------------------------------------------------#
epsilon_gesamt = epsilon_h
epsilon_elastisch = epsilon_elastisch(epsilon_gesamt)
sigma = as_tensor(lamda*epsilon_elastisch[k,k]*delta[i,j] + 2.0*mu*epsilon_elastisch[i,j] , (i,j))
#----------------------------------------------------------------------------#
epsilon_elastisch =  project( epsilon_elastisch, TensorSpaceCG1, solver_type="cg",form_compiler_parameters={"cpp_optimize":True} )

sigma =  project( sigma, TensorSpaceCG1, solver_type="cg",form_compiler_parameters={"cpp_optimize":True} )
#----------------------------------------------------------------------------#
file_verschiebung = File('verschiebung.pvd')
file_verschiebung << u_h

file_stress = File('stress.pvd')
file_stress << sigma