Boundary conditions in a finite element eigenvalue problem

Question

I've been reading multiple papers and related posts for a while now, but I can't seem to find a specific answer to the issues I'm having so I hope someone can clarify things here. I'll provide some background first before diving into the issues at hand.

Background

Basically, I'm developing a 1D finite element code which calculates eigenvalues and eigenvectors after performing a Fourier analysis on a set of equations. There are quite a few equations/variables, but for simplicity say there are 3 of each and name the variables $u_1, u_2$ and $u_3$. I'm using either quadratic or cubic basis functions depending on the variable (see the figure below, left is quadratic, right is cubic)

Hence it basically comes down to this for a single variable: due to the nature of the basis functions we have a 4x4 matrix, subdivided into four 2x2 blocks that are combining the basis functions, where for example $h3$ stands for $C3$ or $Q3$ depending on the variable.

Here $\alpha$ is a combination of the finite element integrals evaluated using Gaussian quadrature etc. Hence, for 3 variables, this yields a 12x12 block for every point in the grid. Matrix assembly is quite straightforward using the above scheme, and eventually this results in a matrix eigenvalue problem of the form $ A\mathbf{x} = \omega B \mathbf{x}$, where $A$ and $B$ are block-tridiagonal matrices. $B$ is real and symmetric, $A$ is complex and is, in general, not symmetric. Say for simplicity that the domain is $x \in [0, 1]$.

Boundary conditions

Partial integration of the finite element integrals introduces natural boundary conditions, which are already added to their respective elements in the $A$ matrix during assembly. Applying for example Dirichlet boundary conditions on both sides of the domain is easy for 'standard' finite element problems. However, here we have a general eigenvalue problem. This post contains a good overview of the possible approaches. I should note that I use method 2 described there since my problem can vary, which implies that reducing the system would mean handling the different cases separately and I would like to keep it general. Also, the matrix $B$ is not allowed to be singular, since it boils down to solving the problem $B^{-1}A\mathbf{x} = \omega \mathbf{x}$, as for example described here.

Now, assume that $u_1$ and $u_2$ are quadratic, and $u_3$ is cubic. The Dirichlet conditions could state that $u_2(0) = u_2(1) = u_3(0) = u_3(1) = 0$, and I'm handling that by omitting the basis functions that are not zero on the boundaries. For the cubic ones those are $C1$ and $C2$, for the quadratic ones $Q3$ and $Q4$. Looking at the grid above this means setting the odd rows/columns to 0 for the cubic variable $u_3$ and the even rows/columns to 0 for the quadratic variable $u_2$. Nothing is done for $u_1$. Afterwards a 1 is introduced on $B$'s diagonal and a large number on $A$'s diagonal to yield $u_i = (\omega / s)u_i$ where $s$ is the large number, resulting in $u_1 = 0$ for $\omega \neq s$. This works out quite nicely, since the eigenfunctions of $u_2$ and $u_3$ are zero on the boundaries as expected and look fine.

The question

Say I want to replace one Dirichlet condition on the right side by a Neumann condition, for example $\partial_x u_2(1) = 0$. This is again not hard for a 'standard' problem, but I don't see how to do this for an eigenvalue problem. How would one approach implementing this? I'm having two thoughts about this:

1. Zero out rows/columns containing the derivative of the basis functions for that specific variable.
2. Perhaps these enter through the natural boundary conditions, which would have to be handled explicitly in that case instead of just adding them to the A-matrix. See also this post.

Is one of these the correct way to handle Neumann conditions?

---

I'm using LAPACK routines and Fortran to assemble and solve the system, by the way. I'm also linking these posts, which are related (and quite informative) but didn't help:

1. Generalised eigenvalue problem
2. Quadratic eigenvalue problem
3. Boundary conditions and Galerkin method

---

Edit: As per @nicoguaro's suggestion I split the OP and kept the most relevant question here, the others will be linked later on.

Answer 1

The boundary conditions don't depend on the choice of your basis but on the formulation you have for your problem. If you have a "standard" finite element formulation, you don't need to do anything to apply (homogeneous) Neumann boundary conditions, they are already satisfied by your system.

In the most common formulation, Neumann boundary conditions are natural boundary conditions. That means that they appear in your variational formulation and are satisfied implicitly by your problem. For the Helmholtz equation, the weak form is something like the following

$$\int\nabla u\nabla w\, \mathrm{d}V + \int w\frac{\partial u}{\partial \hat{\mathbf{n}}}\,\mathrm{d}S = \omega^2 \int u w\, \mathrm{d}V\, ,$$

but in the case of homogeneous boundary conditions, that you should have in an eigenvalue problem, the second term in the left-hand-side is zero because $\partial u/\partial \hat{\mathrm{n}}=0$.