Assembly error involving dgjumps for MatrixValued fes

I am getting the error message when trying to Assemble the following BilinearForm involving dgjumps of MatrixValued fes:

NgException: SparseMatrixTM::AddElementMatrix: illegal dnumsin Assemble BilinearForm ‘biform_from_py’

from netgen.geom2d import SplineGeometry #To create 2D geometries
from netgen.meshing import MeshingParameters #To generate meshes
from ngsolve import *

def hesse(v):
    return v.Operator("hesse") #Broken Hessian

# Geometry
geo = SplineGeometry()
geo.AddCircle(c=(0,0), r=1, bc="circ")
ngmesh = geo.GenerateMesh(maxh=1, quad_dominated=False)
mesh = Mesh (ngmesh)
n = specialcf.normal(mesh.dim)

Xh = H1(mesh, order = 2, dgjumps = True) #Space for displacements
Xl = MatrixValued(H1(mesh, order=1, dgjumps=True),2)

y1= Xh.TrialFunction()
x1 = Xl.TestFunction()

# DG jump and mean
jump_grady1 = grad(y1)-grad(y1).Other()
meanx1_n = 0.5 * (x1+x1.Other()) * n
dS = dx(skeleton=True)

# Bilinear Form
B2 = BilinearForm(trialspace=Xh,testspace=Xl)
B2 += -( InnerProduct(hesse(y1),x1) )* dx
B2 += (meanx1_n*jump_grady1)*dS
B2.Assemble()

It seems that the NGSolve function MatrixValued is not properly passing the dgjump=true flag. Is there any way to fix this?

If one tries computing the Transpose of the Matrix B2 by inverting the role of Trial and Test functions, I am simply getting a

Segmentation fault (core dumped)

from netgen.geom2d import SplineGeometry #To create 2D geometries
from netgen.meshing import MeshingParameters #To generate meshes
from ngsolve import *

def hesse(v):
    return v.Operator("hesse") #Broken Hessian

# Geometry
geo = SplineGeometry()
geo.AddCircle(c=(0,0), r=1, bc="circ")
ngmesh = geo.GenerateMesh(maxh=1, quad_dominated=False)
mesh = Mesh (ngmesh)
n = specialcf.normal(mesh.dim)

Xh = H1(mesh, order = 2, dgjumps = True) #Space for displacements
Xl = MatrixValued(H1(mesh, order=1, dgjumps=True),2)

y1= Xh.TestFunction()
x1 = Xl.TrialFunction()

# DG jump and mean
jump_grady1 = grad(y1)-grad(y1).Other()
meanx1_n = 0.5 * (x1+x1.Other()) * n
dS = dx(skeleton=True)

# Bilinear Form
B2 = BilinearForm(trialspace=Xl,testspace=Xh)
B2 += -( InnerProduct(hesse(y1),x1) )* dx
B2 += (meanx1_n*jump_grady1)*dS
B2.Assemble()