Error in mixed hybrid Poisson equation with discontinuous Lagrange trace

Hi.

I tried to implement hybridized mixed method for the Poisson equation using discontinuous Lagrange trace space. However, in contrast to standard mixed method, the method does not give a correct answer. I used Dolfin 1.6.0. Please take a look.

Added: Standard mixed method and the hybridized one in this code have to give same answer for sigma and u but errors of the numerical solutions with this code are "nan".

Regards, Jeonghun J. Lee

# hybridized mixed method for Poisson equation

from dolfin import *
import numpy


#exact solution
u_ex = Expression('x[0]*(1-x[0])*sin(pi*x[1])', degree = 3)
sigma_ex = Expression(('(1 - 2*x[0])*sin(pi*x[1])', 'pi*x[0]*(1-x[0])*cos(pi*x[1])'), degree = 3)
f = Expression('-sin(pi*x[1])*(2+pi*pi*x[0]*(1-x[0]))', degree = 3)


def compute(nx):
    mesh = UnitSquareMesh(nx,nx)

    # define function spaces
    Mh = VectorElement("DG", mesh.ufl_cell(), 1)
    Vh = FiniteElement("DG", mesh.ufl_cell(), 0)
    Qh = FiniteElement("Discontinuous Lagrange Trace", mesh.ufl_cell(), 1)
    element = MixedElement([Mh,Vh,Qh])
    V = FunctionSpace(mesh, element)

    #def uex_bdy(x, on_boundary):
    #    return on_boundary

    bc = DirichletBC(V.sub(2), u_ex, "on_boundary")

    # test and trial functions
    sigma, u, p = TrialFunctions(V)
    tau, v, q = TestFunctions(V)

    # bilinear form
    n = FacetNormal(mesh)

    a0 = inner(sigma, tau)*dx + u*div(tau)*dx + div(sigma)*v*dx
    a1 = avg(p)*jump(tau,n)*dS + jump(sigma,n)*avg(q)*dS + p*dot(tau, n)*ds + dot(sigma,n)*q*ds

### To verify convergence of standard mixed method
#    Mh = FiniteElement('BDM', mesh.ufl_cell(), 1)
#    Vh = FiniteElement('DG', mesh.ufl_cell(), 0)
#    element = MixedElement([Mh, Vh])
#    V = FunctionSpace(mesh, element)
#  
#    (tau, v) = TestFunctions(V)
#    (sigma, u) = TrialFunctions(V)
#  
#    a = inner(sigma, tau)*dx + inner(u, div(tau))*dx + inner(div(sigma), v)*dx
###

    L = f*v*dx

    #A = assemble(a)
    A0 = assemble(a0)
    A1 = assemble(a1)

    A = A0 + A1
    b = assemble(L)
    bc.apply(A, b)
    U = Function(V)
    solve(A, U.vector(), b)
    (sigma, u, p) = U.split(deepcopy = True)

    # Compute errors
    V1 = FunctionSpace(mesh, 'BDM', 2) 
    sigma_i = interpolate(sigma_ex, V1)
    V2 = FunctionSpace(mesh, 'DG', 1)
    u_i = interpolate(u_ex, V2)
    sigma_err = sqrt(abs(assemble(inner(sigma_i - sigma, sigma_i - sigma)*dx)))
    u_err = sqrt(assemble((u_i - u)*(u_i - u)*dx))
    errors = {'u_L2_err': u_err, 'sigma_L2_err': sigma_err}
    print 'u_err', u_err
    print 'sigma_err', sigma_err

    return errors

h = []  # element sizes
E = []  # errors


for nx in [4, 8, 16, 32, 64]:
  h.append(1.0/nx)
  E.append(compute(nx))  # list of dicts

# Compute convergence rates
from math import log as ln          # log is a dolfin name too
error_types = E[0].keys()
for error_type in sorted(error_types):
  l = len(E)
  print '\nError norm based on', error_type
  for i in range(0, l):
    if i==0:
      print 'Mesh size=%.10f Error=%.10f Convergence rate= -- ' % (h[i], E[i][error_type])
    else:
      r = ln(E[i][error_type]/E[i-1][error_type])/ln(0.5)
      print 'Mesh size=%.10f Error=%.10f Convergence rate=%.2f' % (h[i], E[i][error_type], r)

Comments (15)