Lumping scheme using tri6 element with integration order=2

Dear all,

I am trying to use tri6 elements in a 2D wave propagation problem like wave.py in the example/usersguide folder.
When I set integration order to 2, the HRZ lumping doesn't work, giving an error message: "Lumped mass matrix has zero entry (try order 1 elements or HRZ lumping)". I tried tri3 and rec8 elements with 2nd order integration. Both worked properly. I would like to know why this is happening and what I should do to make the script work.

Here is a minimal script:

from esys.pycad import *
from esys.pycad.gmsh import Design
from esys.finley import MakeDomain
from esys.escript import *
from esys.escript.linearPDEs import LinearPDE,SolverOptions

lam=3.462e9; mu=3.462e9; rho=1154.; alpha=0.7; t0=3.
dim =2; lx = 1.2; ly = 1.2 # sample size, 1.2m by 1.2m
p0=Point(.0, .0, .0); p1=Point(.0, ly, .0); p2=Point(lx, ly, .0); p3=Point(lx, .0, .0)
l01=Line(p0, p1); l12=Line(p1, p2); l23=Line(p2, p3); l30=Line(p3, p0)
c=CurveLoop(l01, l12, l23, l30); rec = PlaneSurface(c)
d=Design(dim=dim, element_size=lx/2., order=2)
d.addItems(rec)
mydomain=MakeDomain(d,integrationOrder=2)
x=mydomain.getX()
mypde=LinearPDE(mydomain)
mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)
kron=kronecker(dim)
dunit=numpy.array([1.,0.])
mypde.setValue(D=kron*rho, q=whereZero(x-[lx,ly])*[1,0])
u=Vector(0.,Solution(mydomain))
g=grad(u)
stress=lam*trace(g)*kron+mu*(g+transpose(g))
mypde.setValue(X=-stress, r=dunit*1e-3)
a=mypde.getSolution()

Cheers,

Hongyang