hybridisation.py

mesh = firedrake.UnitSquareMesh(1, 1)
RTe = firedrake.FiniteElement("RT", firedrake.triangle, 1)
BrokenRT = firedrake.FunctionSpace(mesh, firedrake.BrokenElement(RTe))
DG = firedrake.FunctionSpace(mesh, "DG", 0)
TraceRT = firedrake.FunctionSpace(mesh, "CG", 1)

W = BrokenRT * DG
sigma, u = firedrake.TrialFunctions(W)
tau,v  = firedrake.TestFunctions(W)

# u = firedrake.TrialFunction(DG)
# v = firedrake.TestFunction(DG)

lambdar = firedrake.TrialFunction(TraceRT)
gammar = firedrake.TestFunction(TraceRT)

n = firedrake.FacetNormal(mesh)

mass = firedrake.dot(sigma, tau)*firedrake.dx

div = firedrake.div(sigma)*v*firedrake.dx
grad = firedrake.div(tau)*u*firedrake.dx

trace = firedrake.jump(tau, n=n)*lambdar('+')*firedrake.dS
trace_ext = firedrake.dot(tau, n)*lambdar*firedrake.ds

positive_trace = firedrake.dot(tau, n)('+')*lambdar('+')*firedrake.dS

cell_hdiv = mass + div + grad
Cell_hdiv = firedrake.assemble(cell_hdiv, nest=False).M.values

Trace = Matrix(trace_ext)

X = Trace.T * Matrix(cell_hdiv).inv * Trace
trace = firedrake.assemble(trace_ext, nest=False).M.values
print 'Global Schur complement'
print np.dot(trace.T, np.dot(np.linalg.inv(Cell_hdiv), trace))
print 'Cellwise schur complement'
print assemble(X).values

Global Schur complement
[[  1.0000e+00   5.0000e-01   5.5511e-17   5.0000e-01]
 [  5.0000e-01   1.0000e+00   5.0000e-01   0.0000e+00]
 [  2.7756e-16   5.0000e-01   1.0000e+00   5.0000e-01]
 [  5.0000e-01   0.0000e+00   5.0000e-01   1.0000e+00]]
Cellwise schur complement
[[  1.0000e+00   5.0000e-01   1.3878e-16   5.0000e-01]
 [  5.0000e-01   1.0000e+00   5.0000e-01   0.0000e+00]
 [  1.3878e-16   5.0000e-01   1.0000e+00   5.0000e-01]
 [  5.0000e-01   0.0000e+00   5.0000e-01   1.0000e+00]]