// allocation de memoire
//stacksize(4000000);
//taille de la discretisation
Nx=44;Ny=44;
//taille physique
Lx=5;Ly=1;
// pas de discretisation
dx=Lx/(Nx-1);
dy=Ly/(Ny-1);
//
// Exemple:Construction d'une matrice creuse diagonale
 Id=diag(sparse(ones(Nx*Ny,1)),0);
//


//Matrice locale de coefficients du systeme (contribution sur un element)
// la matrice locale depend des coordonnees des 4 sommets
// ici chaque element a les memes dimensions=> matrice constante
rig=diag((1/dx^2+1/dy^2)*ones(4,1),0);
rig(1,2)=-1/dx^2;
rig(1,4)=-1/dy^2;
rig(2,1)=-1/dx^2;
rig(2,3)=-1/dy^2;
rig(3,2)=-1/dy^2;
rig(3,4)=-1/dx^2;
rig(4,3)=-1/dx^2;
rig(4,1)=-1/dy^2;


// initialisationd de la matrice globale
matg=diag(sparse(zeros(Nx*Ny,1)),0);
//matg=Id;// variante d'initialisation
// boucle sur les elements (au nombre de (Nx-1)*(Ny-1)
for i=1:Nx-1,
  for j=1:Ny-1
    // liste des indices de sommets de l'element considere (peut etre remplace par un tableau donne par un mailleur
      s=zeros(4);
      s(1)=(j-1)*Nx+i; s(2)=(j-1)*Nx+i+1;s(3)=(j)*Nx+i+1;s(4)=(j)*Nx+i;
      // on pourrait appeler ici la construction de rig si celle-ci depend des sommets
      
      // assemblage de la matrice globale
      for iloc=1:4,
        for jloc=1:4
           matg(s(iloc),s(jloc))=matg(s(iloc),s(jloc))+rig(iloc,jloc);
        end  
      end    
  
  end
end
// calcul du terme source
smm=-ones(Nx,Ny);
smm(Nx/4+1:3*(Nx/4),Ny/4+1:3*(Ny/4))=1;// ici on utilise une affectation litigieuse acceptable seulement car sm a ete initialise avant
sm=matrix(smm,Nx*Ny,1);//reshape en matlab
sm=sm-ones(1,Nx*Ny)*sm/(Nx*Ny);// on a retranche sa moyenne
// initialisation element fini: sm a partir du calcul des contributions sur chaque element
// .......
f=ones(Nx*Ny,1);
f(Nx*Ny/2+1:Nx*Ny)=-1;
sm=zeros(Nx*Ny,1);
for i=1:Nx-1,
  for j=1:Ny-1
    // liste des indices de sommets de l'element considere (peut etre remplace par un tableau donne par un mailleur
      s=zeros(4);
      s(1)=(j-1)*Nx+i; s(2)=(j-1)*Nx+i+1;s(3)=(j)*Nx+i+1;s(4)=(j)*Nx+i;
      // assemblage du terme source
      for iloc=1:4,
           sm(s(iloc))=sm(s(iloc))+f(s(iloc))*0.25; 
      end    
  
  end
end
sm=sm-ones(1,Nx*Ny)*sm/(Nx*Ny);// on a retranche sa moyenne
// solution de Matg u = sm
u=matg\sm;
// affichage de la solution
utab=matrix(u,Nx,Ny);
xi=dx*[0:Nx-1]';
yj=dy*[0:Ny-1]';
plot3d(xi,yj,utab);

