// example39.edp from LaplaceRT.edp
mesh Th=square(100,100);
fespace Vh(Th,RT0);
fespace Ph(Th,P0);

func gd = 1.;
func gn = 1.;
func f = 1.;

Vh [u1,u2],[v1,v2];
Ph p,q; 

problem laplaceMixed([u1, u2, p], [v1, v2, q], solver=UMFPACK) =
   int2d(Th)( p*q*0e-10 + u1*v1 + u2*v2 + p*(dx(v1) + dy(v2)) + 
             (dx(u1) + dy(u2))*q )
 + int2d(Th) ( f*q )
 - int1d(Th,1,2,3)( gd*(v1*N.x + v2*N.y) )   //  int on gamma 
 + on(4, u1 = gn*N.x, u2 = gn*N.y);

 laplaceMixed;

 plot([u1,u2], coef=0.1, wait=true, value=true); 
 plot(p, fill=true, wait=true, value=true);