\Delta \mathbf{u} = 0, \;\;\Delta c = 0, with the boundary condition \mathbf{u} = \nabla c and \mathbf{n} \cdot \nabla c =0,

where \mathbf{n} is the normal vector of the boundary. The vector \mathbf{u} has two components u and v.

I want to solve \mathbf{u} and c in a coupling way.

This boundary condition is analogous to the classical Robin boundary condition. I know how to implement the Robin boundary condition in FreeFEM++, but I have no idea how to implement this boundary condition.

I would appreciate it greatly if anyone can give me some suggestions or hints.

Since your problem is vectorial, you should write the Robin boundary condition in a vector-matrix form. As it is now, it is incomplete (as far as I understand, of course). Can you write another equation for \mathbf{n} \cdot \nabla u depending on (u,v)?

Thanks for your kind reply. You are right, the minimum problem I raised before is incomplete. I just update the question (see the above-updated question). What really confuses me is the boundary condition \mathbf{u} = \nabla c. Thank you in advance.

Sorry for the late reply,
Actually your variable is the vector (u,c) and the boundary term involves both \textbf{n} \cdot \nabla u and \textbf{n} \cdot \nabla c. Right now you are not specifying anything about \textbf{n} \cdot \nabla u.
The boundary condition \nabla c = \textbf{u} doesn’t make much sense in this context. Are you sure of it?

Sorry, the minimum problem I simplified is incomplete, so I will give the original problem:

I want to perform a linear stability analysis. The following is the linearized perturbation equation. I try to use the shift-invert method to find eigenvalues.

I learned several example codes solving eigenvalues and found that we can only construct the matrix \mathbf{A} (\mathbf{A}u=\lambda\mathbf{B}u) in a coupling way.

So, when I implemented this problem in FreeFEM++, the boundary condition marked in blue stopped me.

I would greatly appreciate it if you could give me some suggestions on how to deal with this boundary condition. Can the FreeFEM++ solves eigenvalues in a decoupled way?

You should try to write down a weak form and precise a few things: 2D or 3D problem (though I guess 2D)? What is p in line 1 (it looks like the Stokes problem)? Divergence of \mathbf{u} is zero: how do you enforce that property, by a particular type of finite element?
Your equation in blue, if correct, suggests \mathbf{u} is tangential on C1.

We write down the weak form of the linearized governing equation (As you said, the first equation is the classical Stokes equation) and try to implement them in FreeFEM. But I fail to obtain eigenvalues because I do not know how to implement the coupling equation marked in blue.

Do you have any suggestions for the implementation of this tangential boundary condition?

We should try to write a small FreeFem++ script from your equations. For that:
You have variables (p,\mathbf{u},c). What finite elements did you choose? Do you add a stabilization term? (see FreeFEM - Stokes) It also seems a boundary term is missing in the first equation.
The boundary conditions are incomplete: no equation for p and the two equations for C1 are self-contradicting (or we can simplify the blue one on the boundary).
Finally: what is the eigenvalue in these equations?

(1) \sigma in the above equation is the eigenvalue value, c = c_b + \hat{c}e^{\sigma t} .

(2) I choose P2 and P1 for (\mathbf{u}, c) and p, repectively. Due to the absence of reference pressure, I add a stabilization term just as you suggested.

(3) The velocity boundary condition on C1 can be further simplified as \hat{\mathbf{u}} |_{C1} = \nabla c, using the concentration boundary condition \mathbf{n}\cdot\nabla c = 0.

(4) I try to implement the velocity boundary condition in FreeFEM++

As I wrote before, the first line of your variational formulation misses a boundary term: \int_{\partial \Omega} (\frac{\partial \mathbf{u}}{\partial \mathbf{n}} - p \mathbf{n}) \cdot \mathbf{v}.
The boundary condition for \mathbf{u} cannot be simply replaced in this boundary term, it seems.
You must specify a boundary condition for pressure p as well.
The eigenvalue equation you obtain has a rather singular right-hand-side (matrix B would have non-zero elements only for the dofs corresponding to concentration c).
My suggestion would be to try an iterative solution where you solve in sequence the Stokes problem (with the velocity BC working because now the concentration is fixed at every iteration step) then the concentration equation (with fixed velocity at every step). But I don’t know if it would converge or not.

(1) We drop \int \left ( \frac{\partial \bar{\mathbf{u}}}{\partial \mathbf{n}} -\bar{ p}\mathbf{n} \right) since \mathbf{v} = 0 on a Dirchlet boundary condition for \bar{\mathbf{u}}.

(2) Do you know any FreeFEM example that solves the eigenvalue iteratively?

(1) Then you cannot at the same time impose \mathbf{u}=0 and the BC in blue, can you?
(2) No, I don’t no a specific example in FreeFem++. As a remark, the power iteration method (Power iteration - Wikipedia) or the Rayleigh quotient iteration (Rayleigh quotient iteration - Wikipedia) are classical examples.