//  programme avec des elements P2   et dU/dx(1) = 0 et U(0) = 0
//
//  on charge les programmes dont on a besoin
//
exec('P2.sci'); exec('gauss.sci');
//
// nbre points pour l'integration numerique
Ng = 3;
[XG,PG] = gauss(Ng); 
//
//
Nelt = 10; // nbre d'elements

Nno = 3; // nbre de noeuds par element

Npt = 2*Nelt+1; // nbre de points

Dx = 1./(Npt-1);


COOR = matrix((0:Dx:1),Npt,1);  //  Coordonnees des points

CONNEC = zeros(Nelt,Nno);  // matrice de connectivite      
for i=1:Nelt
   CONNEC(i,1) = 1 + 2*(i-1); 
   CONNEC(i,2) = 3 + 2*(i-1);  
   CONNEC(i,3) = 2 + 2*(i-1); 
end;
//
//    Fin initialisation des tableaux
//

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

for iel = 1:Nelt    
   // calcul du syteme elementaire sur chaque element
   A_el = zeros(Nno,Nno);
   F_el = zeros(Nno,1);

   x1 = COOR(CONNEC(iel,1)); x2 = COOR(CONNEC(iel,2));
   jac = (x2 - x1)/2.;   // le jacobien est constant sur l'element

   for ig=1:Ng    // boucle pour l'integration numerique
      [FP,DxFP] = P2(XG(ig));
      
      DxFP = DxFP/jac;

       for i=1:Nno
          F_el(i) = F_el(i) + PG(ig)*FP(i)*jac;
          for j=1:Nno
             A_el(i,j) = A_el(i,j) + PG(ig)*(DxFP(j)*DxFP(i) + FP(i)*FP(j))*jac;
          end;
       end;
    end;
   // fin calcul systeme elementaire pour l'element iel

// on doit assembler

   for i=1:Nno
      i_glob = CONNEC(iel,i);
      F_glob(i_glob) = F_glob(i_glob) + F_el(i);
      for j=1:Nno
         j_glob = CONNEC(iel,j);
         A_glob(i_glob,j_glob) = A_glob(i_glob, j_glob) + A_el(i,j);
      end;
    end;
end;  // fin boucle sur les elements

//
// on impose les deux conditions aux limites de Diriclet
//   U(0) = 0

A = [zeros(1,Npt); A_glob(2:Npt,:)]; // <======
A(1,1) = 1;                          // <======
F = [0; F_glob(2:Npt)];              // <======

u_cal = A\F;

xbasc();
x = COOR;
plot2d(x, u_cal);

// comparaison avec la solution exacte
//
x = COOR;
u_exact = 1- (exp(x) + exp(2-x))/(1+exp(2)); 

erreur = norm(u_cal-u_exact)/norm(u_exact)

xbasc();
plot2d(x,u_cal,style=-1);
x = linspace(0,1,100);
u_exact = 1- (exp(x) + exp(2-x))/(1+exp(2));
plot2d(x,u_exact,style=1);
xtitle('Solution P2; Nelt ='+string(Nelt)+', Npt ='+string(Npt),'x','u');


