3
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ Welcome to Scicomp.SE! You haven't actually written what's going wrong... $\endgroup$ – Christian Clason Apr 20 '15 at 10:26
  • 2
    $\begingroup$ What's the domain you are solving on? From what you describe, I suspect it's a problem with the boundary conditions. Note that the biharmonic equation is a fourth-order equation and hence you need two boundary conditions. The natural boundary conditions are $\Delta u=0$ on $\Gamma_D$, so that's what your weak formulation imposes (and $u=\sin(4\pi x)\sin(4\pi y)$ presumably satisfies). Depending on the domain, this isn't valid for your other $u$ and so you get the wrong solution. In this case, you need to add the corresponding nonhomogeneous Neumann condition to your weak form. $\endgroup$ – Christian Clason Apr 20 '15 at 10:51
  • $\begingroup$ Sorry, that should have been $\partial_n u = 0$. $\endgroup$ – Christian Clason Apr 20 '15 at 11:17
  • $\begingroup$ I don't understand what you mean by "I am trying to solve" when you pre-define the solution, u. Is u a so-called "manufactured solution" where you substitute u into your PDE and Dirichlet BC to compute f and g and then use that f and g to solve the FE equations? $\endgroup$ – Bill Greene Apr 20 '15 at 17:45
  • $\begingroup$ @ChristianClason: My domain is $[0 \times 1]^2$. ahh may be that's where I am going wrong. However, I thought that since I am solving dirichlet problem, I won't need to add the neumann boundary integrand? $\endgroup$ – user21674 Apr 21 '15 at 10:55
2
$\begingroup$

I agree with Christian Clason that your problem is likely with the boundary conditions. $\partial_n u$ is another essential boundary condition (like $u$) so, in general, you will also need to use your penalty approach to satisfy it. Apparently, the prescribed $\partial_n u$ for your $u=\sin(4\pi x)\sin(4\pi y)$ case is consistent with the zero natural boundary conditions that you are getting by default from the weak form.

$\endgroup$

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.