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.

Improve this question
$\endgroup$
8
  • 1
    $\begingroup$ Welcome to Scicomp.SE! You haven't actually written what's going wrong... $\endgroup$
    – Christian Clason
    Apr 20, 2015 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, 2015 at 10:51
  • $\begingroup$ Sorry, that should have been $\partial_n u = 0$. $\endgroup$
    – Christian Clason
    Apr 20, 2015 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, 2015 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, 2015 at 10:55

1 Answer 1

Reset to default
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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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