Composite spaces, on() and tgv=-1

General Discussion
1

It seems that there is a bug on composite spaces when using them with Dirichlet BC set with on() and tgv=-1.
In the example of “Stokes with P2-iso-P1 elements” given in the documentation, setting tgv=-1 gives a wrong result.
Writing the matrix by blocks
M=\pmatrix{A & B \\ B^T & C},
the lines of A corresponding to an index i involved in the boundary condition are set correctly with 1 on the diagonal and zero otherwise.
But the line i of B should be set to zero and it is not the case.

2

That is indeed a bug, thank you for reporting it. I will let you know when it is fixed.

3

Thank you.
Here is indeed a simplified test code to see the problem

int nn = 1; // number of edge in each direction
mesh Th=square(nn,nn);
//savemesh(Th,"mesh.msh");
//plot(Th,wait=true);

fespace Vh(Th,P1);
fespace Xh=Vh*Vh;

tgv=-1;
varf test(<[u],[p]>,<[v],[q]>)
=
//-int2d(Th)(u*v)
//int2d(Th)(p*q)
//int2d(Th)(-u*q)
int2d(Th)(-p*v)
+int2d(Th)(v-q)
+on(1,u=1)
;

matrix M = test(Xh,Xh);
real[int] b = test(0,Xh);
//set(M,solver=UMFPACK);

cout << "M= " << M << endl;
cout << "b= " << b << endl;
4

Hi professor,

Much thanks for this discussion that help me find the bugs when I do an ALE problem, I first use tgv=-1 and forget to enforce the boundary condition to the off-diagonal block of the Jacobian matrix, the convergency of Newton iteration is very slow.

My method is to set tgv=1e30, not exact Dirichlet boundary condition but the results are acceptable and easy to do parallel computation, not needed to update varf or Mat. Another way I find is to zero-out the corresponding rows in the off-diagonal block by hand, as in this post https://community.freefem.org/t/preferred-way-to-construct-block-matrices-with-exact-dirichlet-bc-tgv-1/222/2.