// 
//  Reponse ex2 si le maillage est fait par des Triangles
//
clear
xbasc(0); xbasc(1);
//
//  on charge les programmes dont on a besoin
//
exec('P1.sci'); 
exec('hammer.sci'); 
exec('maillage_2D_P1.sci');
// Pour tracer
exec('trace_maillage_2D.sci'); exec('trace_u_2D.sci');
exec('cal_indice.sci');  exec('interpol.sci'); 
//
nx = 8;  ny = 8;  // nbre d'element dans les directions x et y
lx = 1.; ly = 1.; // taille du carre

[CONNEC,COOR,Nelt,Nno,Npt] =  maillage_2D_P1(1,1,nx,ny);
trace_maillage_2D(CONNEC,COOR,0);
// 
//  nbre de points pour l'integration numerique

Ng = 4; [XG,PG] = hammer(Ng);

//
//  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 x <> 0 & x <> 1 & 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);  
for i=1:Npt
  x = COOR(i,1); y = COOR(i,2);
  if y == 0 | y == 1 | x == 0 | x == 1 then
     UC(NUMER(i)-NI) = exp(2*x) * sin(y);
  end;
end;

//
//    Fin initialisation des tableaux
//
A_glob = spzeros(Npt,Npt); F_glob = zeros(Npt,1);

for iel = 1:Nelt  // boucle sur les elements
//  
// calcul du systeme 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); yg = XG(ig);
//
// calcul de la matrice jacobienne
//
      [FP,DxiFP,DetaFP] = P1(xg,yg);
      j11 = x(1:3)'*DxiFP; 
      j12 = x(1:3)'*DetaFP;   
      j21 = y(1:3)'*DxiFP;  
      j22 = y(1:3)'*DetaFP;
//
      DT  = [j11 j12; j21 j22];
      jac = det(DT);
      B = inv(DT)';
   
      for i=1:Nno
         DFPi = B*[DxiFP(i); DetaFP(i)];
         for j=1:Nno
            DFPj = B*[DxiFP(j); DetaFP(j)];
            A_el(i,j) = A_el(i,j) + ...
                        PG(ig)*( DFPj(1)*DFPi(1) + DFPj(2)*DFPi(2) + 3*FP(j)*FP(i) ) *jac;
         end;
      end;
  end;
//
// fin calcul systeme elementaire pour l'element iel
//
// 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,21,21);

// test de comparaison
//
U_c = u_cal(NUMER); x = COOR(:,1); y = COOR(:,2);
U_th = exp(2*x) .* sin(y);

err_r = norm(U_c - U_th)/norm(U_th)


