clear; xdel(0); xdel(1); xdel(2);
//  on charge les programmes dont on a besoin
//
exec('Q1.sci'); exec('Q2.sci'); exec('gauss.sci'); 
exec('maillage_2D_Q2.sci');
exec('trace_maillage_2D.sci'); exec('trace_u_2D.sci');
exec('cal_indice.sci');  exec('interpol.sci'); 

//
//  nbre de points pour l'integration numerique
Ng = 3;
[XG,PG] = gauss(Ng); 
//
nx = 10;  ny = 10;  // nbre d'element dans les directions x et y
lx = 1.; ly = 1.; // taille du carre

[CONNEC,COOR,Nelt,Nno,Npt] =  maillage_2D_Q2(1,1,nx,ny);
trace_maillage_2D(CONNEC,COOR,0);

//
//  creation de la matrice  NUMER et ADRESS
//
//  Conditions aux limites 
//  u = 0 sur tout le bord
NUMER=zeros(Npt,1);
NI = 0;
nc = 0;
for i=1:Npt
   x = COOR(i,1); y = COOR(i,2);
   if  y <> 0 & y <> 1 then            
     NI = NI+1;
     NUMER(i) = NI;
   else
     NUMER(i) = Npt-nc;
     nc = nc + 1;   
   end;
end;      
//
//  creation de la matrice ADRESS
//
ADRESS = zeros(Nelt,Nno);
for i =1:Nelt
  for j = 1:Nno
    ADRESS(i,j) = NUMER(CONNEC(i,j));
  end;
end;
//
//  creation du vecteur UC
//  
UC = zeros(Npt-NI,1);  

//
//    Fin initialisation des tableaux

A_glob = spzeros(Npt,Npt); F_glob = zeros(Npt,1);

for iel = 1:Nelt  // boucle sur les elements
//  
// ---------------------------------------------------------------------------- 
//  PARTIE   A ECRIRE
//-------------------------------------------------------------------
// calcul du syteme elementaire sur chaque element
//
   A_el = zeros(Nno,Nno);
   F_el = zeros(Nno,1);

   x  = COOR(CONNEC(iel,:),1); y = COOR(CONNEC(iel,:),2);
 

   for ig=1:Ng    // double boucle pour l'integration numerique
    xg = XG(ig);
    for jg=1:Ng
       yg = XG(jg);
//
// calcul de la matrice jacobienne
//
          [F,DxiFP,DetaFP] = Q1(xg,yg);
          j11 = x(1:4)'*DxiFP; 
          j12 = x(1:4)'*DetaFP;   
          j21 = y(1:4)'*DxiFP;  
          j22 = y(1:4)'*DetaFP;
//
          DT  = [j11 j12; j21 j22];
          jac = det(DT);
          B = inv(DT)';
//
       [FP,DxiFP,DetaFP] = Q2(xg,yg);
   
       for i=1:Nno
          DFPi = B*[DxiFP(i); DetaFP(i)];
          F_el(i) = F_el(i) + PG(ig)*PG(jg)*FP(i)*jac;
          for j=1:Nno
            DFPj = B*[DxiFP(j); DetaFP(j)];
            A_el(i,j) = A_el(i,j) + PG(ig)*PG(jg)*( DFPj' * DFPi) *jac;
          end;
       end;
    end;
    end;
//
// fin calcul systeme elementaire pour l'element iel
//
//--------------------------------------------------------------
//  FIN PARTIE A ECRIRE
//--------------------------------------------------------------
// on doit assembler

   for i=1:Nno
      i_glob = ADRESS(iel,i);
      F_glob(i_glob) = F_glob(i_glob) + F_el(i);
    end;

   A_el=sparse(A_el);
   [ij,v,mn] = spget(A_el);
   for i = 1:length(v)
     ij(i,1) = ADRESS(iel,ij(i,1));
     ij(i,2) = ADRESS(iel,ij(i,2));
   end;
   A_el = sparse(ij,v,[Npt Npt]); 
   A_glob = A_glob + A_el;

end;  // fin boucle sur les elements


A = A_glob(1:NI,1:NI);  F = F_glob(1:NI);
M12 = A_glob(1:NI,NI+1:$);
F = F - M12*UC;

[h,rk] = lufact(A);
U = lusolve(h,F);
ludel(h);

u_cal = [U; UC];

[u,x,y]=trace_u_2D(u_cal,COOR,NUMER,1,20,20);

//
// comparaison avec la solution y(1-y)/2
//
Nplot = 2*ny+1;
y = COOR(1:Nplot:Npt,2);
u_plot = u_cal(NUMER(1:Nplot:Npt));
u_e = y .*(1 - y)/2;
err = norm(u_plot - u_e)/norm(u_e)

xset("window",2);
xbasc();
plot2d(y, u_plot,style = -1);

y = linspace(0.,1.,50);
u_e = y .*(1 - y)/2;
plot2d(y, u_e,style = 1);
xtitle('Ecoulement de Poiseuille','y', 'u');
