In this code, I used the following code to build a local mesh
real[int, int] nodes = [[coordX0, coordY0], [coordX1,coordY1],
[coordX2, coordY2]];
//cout << nodes << "\n";
border f0(t=0,1){P.x=nodes(1,0)*t + nodes(0,0)*(1-t);P.y=nodes(1,1)*t + nodes(0,1)*(1-t); label=1;};
border f1(t=0,1){P.x=nodes(2,0)*t + nodes(1,0)*(1-t);P.y=nodes(2,1)*t + nodes(1,1)*(1-t); label=2;};
border f2(t=0,1){P.x=nodes(0,0)*t + nodes(2,0)*(1-t);P.y=nodes(0,1)*t + nodes(2,1)*(1-t); label=3;};
Th[K] = buildmesh(f0(localh) + f1(localh) + f2(localh), fixedborder=true);
//---------------------------------------------------Solving for etaf
varf ahKf(etafK, varphi)
= int2d(Th[K])(kappa*grad(etafK)'*grad(varphi));
varf fhKf(etafK, varphi)
= int2d(Th[K])(f*varphi);
varf bhKf(vf, varphi)
= int2d(Th[K])(vf*varphi);
// local matrices for etaf
matrix Akf = ahKf(VhK,VhK);
matrix Bkf = bhKf(V0,VhK);
real[int] Fkf(ndofVhK);
Fkf = fhKf(0, VhK);
// block matrix LHS
matrix AAkf = [ [ Akf, Bkf ],
[ Bkf', 0 ] ];
real[int] rhsetaf(ndofVhK+nDoFu0), soletafk(ndofVhK+nDoFu0);
real[int] blmf(nDoFu0), lmf(nDoFu0);
//second equation to hybrid form equal 0
blmf = 0;
// block vector RHS
rhsetaf = [ Fkf, blmf ];
set(AAkf, solver=sparsesolver);
soletafk = AAkf^-1*rhsetaf; //solve linear system
[etafK[], lmf] = soletafk; //set the value
cout << "etaF elem " << K << ": " << etafK[] << "\n";
I don’t know if its related, but for etaF I should have only one value, instead I’m getting a vector of values as you can see above. Could anyone help me with this?