mini_vs_p2p1.py

from dolfin import *
import math


def stokes(n):
    # Load mesh and subdomains
    mesh = UnitSquareMesh(n, n)

    u_analytical = Expression(("-sin(pi*x[0])", "0.0"))

    # Define function spaces
    P2 = VectorFunctionSpace(mesh, "Lagrange", 2)
    P1 = VectorFunctionSpace(mesh, "Lagrange", 1)
    B  = VectorFunctionSpace(mesh, "Bubble", 3)
    Q  = FunctionSpace(mesh, "CG",  1)
    Mini = (P1 + B)*Q
    P2P1 = P2*Q

    ele = Mini
    #ele = P2P1

    bc = DirichletBC(ele.sub(0), u_analytical, "on_boundary")


    g = Expression("-pi*cos(pi*x[0])")

    # Define variational problem
    (u, p) = TrialFunctions(ele)
    (v, q) = TestFunctions(ele)
    f = Constant((0, 0))
    a = (inner(grad(u), grad(v)) - div(v)*p + q*div(u))*dx
    L = inner(f, v)*dx + g*q*dx

    # Compute solution
    w = Function(ele)
    solve(a == L, w, bc)

    # Split the mixed solution using deepcopy
    # (needed for further computation on coefficient vector)
    (u, p) = w.split(True)

    error = errornorm(u_analytical, u)
    print "Error in u", error

    # Plot solution
    (u, p) = w.split()
    #plot(u)
    #plot(p)
    #interactive()

    return error


ns = [5, 10, 20, 40, 80]
errors = []
for n in ns:
    errors.append(stokes(n))

    if len(errors) > 1:
        print "****** Order of convergence: %s." % (math.log(errors[-2], 2)-math.log(errors[-1], 2))