This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Integral error norms

0 votes

Hi,

I am trying to verify that my FEniCS code is correct with respect to the errors. I have attached a piece of code which solves the Laplacian on a unit square domain with Dirichlet boundary conditions. The boundary conditions are prescribed by u(x,y) = (1,1) so the forcing term is F(x,y) = (0,0). I would expect the errors to be roughly 1e-12 or smaller for both the L2 and H1 integral norms. Where am I going wrong in my code?

Thanks!

from dolfin import *
import numpy as np
import pandas as pd
n = 6
Dim = np.zeros((n,1))
ErrorL2 = np.zeros((n,1))
ErrorH1 = np.zeros((n,1))
OrderL2 = np.zeros((n,1))
OrderH1 = np.zeros((n,1))

for x in range(1,n+1):

    mesh = UnitSquareMesh(2**x,2**x)
    V = VectorFunctionSpace(mesh, "CG", 2)

    class u_in(Expression):
        def __init__(self):
            self.p = 1
        def eval_cell(self, values, x, ufc_cell):
            values[0] = 1.0
            values[1] = 2.0
        def value_shape(self):
            return (2,)
    class F_in(Expression):
        def __init__(self):
            self.p = 1
        def eval_cell(self, values, x, ufc_cell):
            values[0] = 0.0
            values[1] = 0.0
        def value_shape(self):
            return (2,)

    u0 = u_in()
    F = F_in()

    u = TrialFunction(V)
    v = TestFunction(V)

    a = inner(grad(u), grad(v))*dx
    L = inner(F, v)*dx

    def boundary(x, on_boundary):
        return on_boundary

    bc = DirichletBC(V, u0, boundary)
    u = Function(V)

    solve(a == L, u, bcs=bc,
      solver_parameters={"linear_solver": "lu"},
      form_compiler_parameters={"optimize": True})

    Vexact = VectorFunctionSpace(mesh, "CG", 4)
    ue = interpolate(u0, Vexact)

    e = ue - u
    Dim[x-1] = V.dim()
    ErrorL2[x-1] = sqrt(abs(assemble(inner(e,e)*dx)))
    ErrorH1[x-1] = sqrt(abs(assemble(inner(grad(e),grad(e))*dx)))

    if (x > 1):
        OrderL2[x-1] = np.log2(ErrorL2[x-1]/ErrorL2[x-2])
        OrderH1[x-1] = np.log2(ErrorH1[x-1]/ErrorH1[x-2])


TableTitles = ["DoF","L2-erro","L2-order","H1-error","H1-order"]
TableValues = np.concatenate((Dim,ErrorL2,OrderL2,ErrorH1,OrderH1),axis=1)
Table = pd.DataFrame(TableValues, columns = TableTitles)
pd.set_option('precision',3)
print Table
asked May 7, 2016 by mwathen FEniCS Novice (710 points)
...