Error cannot reduce much with increase in mesh size

1

I am implementing biharmonic pde with neumann boundary condition . I am getting error 0.4961486 with mesh size 1/128 which is too high. I don,t understand where I am doing wrong . Kindly help me . Thank You.

Here is my code :

import matplotlib.pyplot as plt
import numpy as np
from dolfin import *

# Next, some parameters for the form compiler are set::

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True

# A mesh is created, and a quadratic finite element function space::

# Make mesh ghosted for evaluation of DG terms
parameters["ghost_mode"] = "shared_facet"

# Create mesh and define function space
mesh = UnitSquareMesh(8,8)
V = FunctionSpace(mesh, "CG", 2)


class Source(UserExpression):
    def eval(self, values, x):
        values[0] = 112*pi**4*(sin(pi*x[0]))**4*(sin(pi*x[1]))**4 + 24*pi**4*(cos(pi*x[0]))**4*(sin(pi*x[1]))**4 + 24*pi**4*(cos(pi*x[1]))**4*(sin(pi*x[0]))**4 - 288*pi**4*(cos(pi*x[0]))**2*(sin(pi*x[0]))**2*(sin(pi*x[1]))**4 -288*pi**4*(cos(pi*x[1]))**2*(sin(pi*x[0]))**4*(sin(pi*x[1]))**2 + 288*pi**4*(cos(pi*x[0]))**2*(cos(pi*x[1]))**2*(sin(pi*x[0]))**2*(sin(pi*x[1]))**2
   
# On the finite element space ``V``, trial and test functions are
# created::

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)

# Define normal component, mesh size and right-hand side
h = CellDiameter(mesh)
h_avg = (h('+') + h('-'))/2.0
n = FacetNormal(mesh)
f = Source(degree=2)


# Penalty parameter
alpha = Constant(6.0)

x = SpatialCoordinate(mesh)
g =(sin(pi*x[0]))*(sin(pi*x[0]))*(sin(pi*x[0]))*(sin(pi*x[0]))*(sin(pi*x[1]))*(sin(pi*x[1]))*(sin(pi*x[1]))*(sin(pi*x[1]))

z = as_vector((4*pi*cos(pi*x[0])*(sin(pi*x[0])**3)*(sin(pi*x[1])**4),4*pi*cos(pi*x[1])*(sin(pi*x[0])**4)*(sin(pi*x[1])**3)))

                           
b = inner(div(z), div(grad(v)))*dx \
  - inner(avg(div(z)), jump(grad(v), n))*dS \
  - inner(jump(z, n), avg(div(grad(v))))*dS \
  + alpha/h_avg*inner(jump(z,n), jump(grad(v),n))*dS

l = (cos(pi*x[0]))*(cos(pi*x[1]))
p = as_vector((((-1)*pi*(cos(pi*x[1]))*(sin(pi*x[0]))),((-1)*pi*(cos(pi*x[0]))*(sin(pi*x[1])))))

c = inner(div(grad(l)), div(grad(v)))*dx \
  - inner(avg(div(grad(l))), jump(grad(v), n))*dS \
  - inner(jump(grad(l), n), avg(div(grad(v))))*dS \
  + alpha/h_avg*inner(jump(grad(l),n), jump(grad(v),n))*dS

a = inner(div(grad(u)), div(grad(v)))*dx \
  - inner(avg(div(grad(u))), jump(grad(v), n))*dS \
  - inner(jump(grad(u), n), avg(div(grad(v))))*dS \
  + alpha/h_avg*dot(jump(grad(u),n), jump(grad(v),n))*dS



# Define linear form
L = b - f*v*dx + c

# A :py:class:`Function <dolfin.functions.function.Function>` is created
# to store the solution and the variational problem is solved::

# Solve variational problem
u = Function(V)
solve(a == L, u)

x = SpatialCoordinate(mesh)
u_ex = (cos(pi*x[0]))*(cos(pi*x[1]))
error_L2 = sqrt(assemble(((u-u_ex)**2)*dx))
print('error_L2  =', error_L2)

# The computed solution is written to a file in VTK format and plotted to
# the screen. ::

# Save solution to file
file = File("biharmonic3.pvd")
file << u

# Plot solution
plot(u)
plt.show()
2

The terms for the exterior boundary might be missing in your linear and bilinear form.

3

I am taking neumann boundary which is included in right hand side.

4

“dS” only applies to interior facets.