//  Programme qui resoud   div(sigma )  = 0 dans un rectangle 
//   Element P1-P1
//
// detruit les anciennes variables et efface les fenetres graphiques
clear 
for i=1:6
   xbasc(i);
end; 
//
//  function scilab  que l'on va utiliser 
//
   getf('maillage_2D_P1.sci'); 
   getf('trace_maillage_2D.sci');  
   getf('hammer.sci');
   getf('cond_limit_2D.sci');
   getf('P1.sci'); getf('bulle.sci'); 
   getf('assemblage.sci');   getf('cal_sol_lin.sci');
//
// pour tracer
//   
   getf('cal_indice.sci');  getf('interpol.sci'); 
   getf('trace_u_2D.sci');  getf('trace_vecteur.sci'); getf('trace_p_2D.sci');
//
//   calcul du systeme elementaire sur l'element 
//
    getf('cal_Melm_P1pP1.sci'); 
//
//  on utilise des variables globales pour 
//
   global CONNEC COOR Nelt Nno Npt;
   global NUMER ADRESS Nno_u Nno_p Npt_u Npt_p;
   global UC NI;
   global XG PG Ng;
//
//  nbre de points pour l'integration numerique
   Ng = 12;
   [XG, PG] = hammer(Ng); 
//
   nx = 10;  ny = 10;  // nbre d'element dans les directions x et y
//
// maillage
//
   [CONNEC,COOR,Nelt,Nno,Npt] = maillage_2D_P1(1,1,nx,ny);
   trace_maillage_2D(CONNEC,COOR,0);

   CONNEC_u = CONNEC; COOR_u = COOR; Nno_u = Nno; Npt_u = Npt;

   CONNEC_p = CONNEC; COOR_p = COOR; Nno_p = Nno; Npt_p = Npt;

   [NUMER, ADRESS, UC, NI] = cond_limit_2D(COOR_u,CONNEC_u,COOR_p,CONNEC_p,Npt_p);
   Ndof = length(NUMER);

//
//     Debut Cacul
//
// valeur initiale
//

   U = zeros(Ndof,1); Fext = zeros(Ndof,1);
//
//  calcul de la solution pour gradP fixe
//
   U = cal_sol_lin(cal_Melm_P1P1,Fext,U);

             trace_vecteur(U,COOR_u,NUMER,3,21,21);
   [u,x,y] = trace_u_2D(U,COOR_u,NUMER,1,4,21,21); 
   [v,x,y] = trace_u_2D(U,COOR_u,NUMER,2,5,21,21); 
   [p,x,y] = trace_p_2D(U,COOR_p,NUMER,Npt_u,6,21,21);

//
// teste si solution analytique
//
   Ux_c = U(NUMER(1:Npt_u)); Uy_c = U(NUMER(Npt_u+(1:Npt_u))); 
   P_c  = U(NUMER(2*Npt_u+(1:Npt_p)));

   x = COOR_u(:,1); y = COOR_u(:,2);
   Ux_th = exp(x) .* (y .* sin(y) - cos(y));
   Uy_th = exp(x) .* y .* cos(y);
  
   x = COOR_p(:,1); y = COOR_p(:,2);
   x1 = x(1); y1 = y(1); 
   P_th = 2 * exp(x) .* cos(y);
  
   test_p = norm(P_c -P_th)/norm(P_th)
   test_x = norm(Ux_c -Ux_th)/norm(Ux_th)
   test_y = norm(Uy_c -Uy_th)/norm(Uy_th)
   
   if test_p > 1 then 
      Val_p = [COOR_p P_c-P_th];
   end; 



