IPC0 Biharmonic problem
Hi,
I hope everyone is having a good start to the new year.
I have the following relatively minimal code for IPC0 biharmonic problem:
from ufl import *
from dune.fem.space import lagrange
from dune.fem import integrate
from dune.fem.scheme import galerkin
from dune.alugrid import aluSimplexGrid
import math
import numpy as np
from dune.ufl import DirichletBC, Constant
from scipy.sparse.linalg import spsolve
laplace = lambda v : div(grad(v))
vertices = np.array( [(-1,-1), (1,-1), (1,1), (-1,1),(0,0)] )
triangles = np.array( [(0,1,4), (1,2,4), (2,4,3), (3,4,0)] )
gridView = aluSimplexGrid({"vertices": vertices, "simplices": triangles})
init_refines = 2
gridView.hierarchicalGrid.globalRefine(init_refines)
order = 2
order_constant = Constant(order)
iterations = 4
for i in range(iterations):
gridView.hierarchicalGrid.globalRefine()
space = lagrange(gridView, order=order)
u = TrialFunction(space); v = TestFunction(space)
x = SpatialCoordinate(space)
normal = FacetNormal(space)
h = FacetArea(space)#problem setup
inverseConstant = Constant(10) # 10 should be fine...
beta = inverseConstant * order_constant**2 / (h)
exact = (sin(math.pi * x[0]) *sin(math.pi*x[1]) )**2 ## for example
f = laplace(laplace(exact))
a = laplace(u)*laplace(v)*dx - f*v *dx
#'clamped' conditions penalised
a+= beta*inner(grad(u-exact),normal)*inner(grad(v),normal)*ds
#interior penalty for jumps of gradients
a+= beta * jump (grad(u),normal) * jump (grad(v),normal)*dS
## non-conforming terms
a+=-avg(laplace(u)) * jump(grad(v),normal) * dS - avg(laplace(v)) * jump(grad(u), normal) * dS
a+=- (laplace(u)) * inner(grad(v),normal) * ds - (laplace(v)) * inner(grad(u), normal) * ds
a-=-laplace(v) * inner(grad(exact),normal) * ds
solScheme = galerkin([a==0,DirichletBC(space,exact)])
solution = space.function(name = "solution")
interpolatedExact = space.interpolate(exact,name = "interpolated")
solScheme.solve(target = solution)
L2 = integrate((solution-exact)**2, gridView, order = 2 * order )
H1 = integrate(inner(grad(solution-exact),grad(solution-exact)), gridView, order = 2 * (order-1) )
H2 = integrate(laplace(solution-exact)**2, gridView, order = 2 * (order-2) )
print("order =", order, "\tL2Error =", math.sqrt(L2),\
"\tH1Error =", math.sqrt(H1),\
"\tH2Error =", math.sqrt(H2))
which does not give the expected convergence. Of course I am not 'correctly' measuring the error.
The ufl form worked out fine when I tried it elsewhere, but maybe that has a slightly different syntax which makes a problem here.
Hopefully this ends up being some embarrassing syntax error on my part, rather than there being some challenging underlying issue.
Best,
Philip