load "splitmesh6"
mesh Th=square(5,5);
mesh Th6=splitmesh6(Th);
mesh Th6b=splitmesh6PowellSabin(Th);// Add for François Bouchut  feb 2024 
plot(Th6,wait=1);
plot(Th6b,wait=1);

fespace Vh(Th,P1);
fespace Nh(Th,P0);

fespace Vh6(Th6,P1);
fespace Nh6(Th6,P0);
fespace RT6(Th6,RT0); // Raviart Thomas ordre 0

//  
varf vM6(u,v) = int2d(Th6,qforder=1)(u*v*3/area);
matrix  M610= vM6(Nh6,Vh6);
matrix  I61 = interpolate(Vh,Vh6);
matrix  S61 =  I61*M610;

Nh6 eta6=1;
Vh  eta;
eta[]= S61*eta6[];
plot(eta,wait=1);

RT6  [u6,v6];
RT6  [uu6,vv6];

Nh6 x6=x,y6=y;

//solve  PP([u6,v6],[uu6,vv6]) = intalledges(Th6)(  u6*uu6*N.x + v6*vv6*N.y)
//- intalledges(Th6)(  mean(x6)*uu6*N.x + mean(y6)*vv6*N.y);

  [u6,v6] = [mean(x6),mean(y6)]; // OK version 2.19 
 [uu6,vv6] = [x6,y6];
 plot( [u6,v6],  [uu6,vv6], wait=1);
 
 
