Vectorial PDE, conditions, stiffness matrix

General Discussion
1

Hello FreeFem community,
I am looking for how to impose conditions in vectorial PDE problems.

More preciesely, I am working on a NS problem. Assuming my variables are [u,v,p], in a domain, say [0,2]x[0,1], I look for how to impose conditions like:

  • u=A in [1,2]x[0,1],
  • v=B in [0,1]x[0,1],
  • p=C on the segment [0,0]-[2,1],

where A, B, C are constants (or functions; these conditions are just for the purpose of presenting the question and we let aside any related mathematical issue).

For scalar problems. say wrt u, I can easily impose my condition - by identifying the index of nodes involved in the condition for u and modifying accordingly the stiffness matrix M and the Rhs vector.

For vectorial problems I can proceed similarly, but I need to know the precise structure of M. For example, to set

  • u(X(i))=A,
  • v(X(j))=B,
  • p(X(k))=C.

Namely, is the matrix organised as N blocks of 3, or 3 blocs of N? So should I do:

  • M(3i,3i) = M(3j+1,3j+1) = M(3k+2,3k+2) = tgv,
  • M(i,i) = M(N+j,N+j) = M(2N+k,2N+k) = tgv,
  • or ?

Here, N=number of nodes (so M is a (3N)x(3N) matrix).

Thank you for your help.

ps: For convenience, here is a piece of code I use
// domain, mesh, FE spaces
//
mesh Th = square(20,10, [2*x, y]); // rectangle [0,2]x[0,1]
fespace Vh(Th,P1b);
fespace Qh(Th,P1);
fespace Zh(Thf,[P1b,P1b,P1]);
//
varf a([u,v,p],[uh,vh,ph])
= int2d(Th)(…) + on(1, u=Cu) + on(2,v=Cb)+…;
//
matrix A(Zh,Zh);
//////
A(i,i)=tgv; // ?
// … ?

2

Degrees of freedom are interleaved, so it’s like N blocks of 3, but not exactly, because you are using different fespaces for the velocities and pressure.

3

Also, don’t forget to do a search on the forum before posting, this is a recurring question with many available answers, e.g., Order of variables in matrix generated with varf with a [P1, P1, P0] fespace.

4

Hello,
Here you can do simply as adding in your varf
+int2d(Th,labreg)(tgv*u*uh)
and a similar term for the rhs, where labreg is the label of your region.

If you don’t like to define region labels you can do (for the left region x<1)

fespace Vh0(Th,P0);
Vh0 regl;
regl=(x<1. ? 1. : 0.);

and use
+int2d(Th)(tgv*u*uh*regl)

It is not easy to handle what happens close to the interface between the regions, since usually u is in a space imposing continuity.

5

Thank you for suggesting this elegant solution.

6

Thank you for the answer and the references.