splitmesh6_edp.edp

load "splitmesh6" 
mesh Th=square(5,5); 
mesh Th6=splitmesh6(Th); 
plot(Th6,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);