I am trying to solve simple scalar biharmonic equation using bubnov-galerkin finite element method. I am using $H^2$ conforming basis functions. I was wondering that if anyone can give me some pointers on how can I further debug my algorithm.
I proceeded by converting the strong form of equation to weak form: \begin{align} \Delta^2 u &= f \\ u &= g \;\;\; \text{on} \;\;\; \Gamma_D \\ \Gamma_N &= \emptyset \end{align}
As an example problem, I am trying to solve: \begin{align} u &= \cos(4 \pi x) \cos(4 \pi y) \end{align} Boundary conditions are implemented using penalty method. The weak form is follows:
\begin{align} a(u,v) &= L(v) \\ a(u,v) &= \sum_{\Omega_K \in \mathcal T} \int_{\Omega_K} \Delta \psi_{i} \Delta \psi_{j} + \sum_{E \in \mathcal E} \gamma \int_{E} \psi_{i} \psi_{h,j} \\ L(v) &= \sum_{\Omega_K \in \mathcal T} \int_{\Omega_K} \psi_{i} f + \sum_{E \in \mathcal E} \gamma \int_{E} \psi_{i} g\\ \end{align}
Problem: When I solve this equation on a square domain, my L2 error barely converges (convergence rate ~ 0.3). However, if I solve $u = \sin(4 \pi x) \sin(4 \pi y)$, I get correct convergence rates. I have tried the following things to debug my code:
I solves the poisson equation, so replaced the first integrand with the stiffness integrand. I get correct convergence rates. I concluded from this that my penalty method implementation is right
I am using Cartesian grid, so Jacobian lines up with expectations
Since, my shape functions are defined in parametric ($\xi,\eta$) space, I had to work out the transformation for the laplacian. I tested this transformation on polar coordinates.
Many thanks.
Apologies. Editted.