Hi all. I am trying to implement primal hybrid formulation. I have the following formulation:
Find a pair (u, \lambda) \in X_{h}\times M_{h} such that
a(u,v)+b(v,\lambda) = f(v)\quad \forall v\in X_{h}
b(u,\mu) = 0 \quad \mu \in M_{h}
where a(u,v)= \sum_{K\in \mathcal{T}_{h}} \int_{K} \nabla u\cdotp \nabla v dx
and b(v,\mu) = -\sum_{K\in \mathcal{T}_{h}} \int_{\partial K} v\mu ds.
X_{h} = \{v \in L^{2}(\Omega) \,\mbox{such that}\, v\mid_{K}\in P_{1}(K) \forall K\in \mathcal{T}_{h} \}
and M_{h} = \{ w\in L^{2}(\partial K) \,\mbox{such that} \, w\mid_{\partial K_{1}} + w\mid_{\partial K_{2}} = 0 \,\mbox{for any adjacent}\, K_{1}\, \mbox{and}\, K_{2} \in \mathcal{T}_{h} \mbox{and}\, w\mid_{\partial K}\in P_{0}(\partial K)\}.
Here is my code:
from fenics import*
import numpy as np
import matplotlib.pyplot as plt
u_ex = Expression('x[0]*(1-x[0])*sin(pi*x[1])', degree=3)
f = Expression('-sin(pi*x[1])*(2+pi*pi*x[0]*(1-x[0]))', degree=3)
mesh = UnitSquareMesh(16,16)
Mh = FiniteElement('DG',mesh.ufl_cell(),1)
Qh = FiniteElement('Discontinuous Lagrange Trace', mesh.ufl_cell(), 0)
element = MixedElement([Mh,Qh])
V = FunctionSpace(mesh,element)
bc = DirichletBC(V.sub(0), u_ex, "on_boundary")
u, p = TrialFunctions(V)
v, q = TestFunctions(V)
a0 = inner(grad(u),grad(v))*dx - avg(p)*(v('+')-v('-'))*dS - (p('+')-p('-'))*avg(v)*dS - p*v*ds
a1 = - avg(q)*(u('+')-u('-'))*dS - avg(u)*(q('+')-q('-'))*dS - q*u*ds
L = f*v*dx
file_u = File('primalhybrid/u.pvd')
A0 = assemble(a0)
A1 = assemble(a1)
A = A0 + A1
b = assemble(L)
bc.apply(A,b)
U = Function(V)
solve(A, U.vector(), b)
(u, p) = U.split(deepcopy=True)
file_u << u
V1 = FunctionSpace(mesh,'DG',1)
u_i = interpolate(u_ex, V1)
error_L2 = errornorm(u_i,u,'L2')
print('error_L2=', error_L2)
plot(u)
plot(u_i)
I think u is not taking the boundary condition correctly, because of that I am not getting convergence.
If someone wants to check these spaces or formulation then one can look at the following paper:
https://www.jstor.org/stable/2006423