%    testfem2.m  :   test de  int de 0 a L   de u'' v'' = F v(L)
%
clear;n=input('nombre de subdivisions '),
xi=symmul('L',symdiv((1:n),n)); xxi=symsub('x',xi);
%
%base de Lagrange:
%  le (2i-1)eme element est relatif a la VALEUR au ieme point,
%  le (2i)eme element est relatif a la DERIVEE  au ieme point.
%
%   Lg  = restrictions au sous-intervalle de gauche
%   Ld  = restrictions au sous-intervalle de droite
Lg=zeros(1,2*n-1);Lg=sym(Lg);        Ld=zeros(1,2*n-1);Ld=sym(Ld);
%   Lg(2i-1)=(1+n(x-xi)/L)^2(1-2n(x-xi)/L) ,
%                      Ld(2i-1)=(1-n(x-xi)/L)^2(1+2n(x-xi)/L) ,
%   Lg(2i)=(1+n(x-xi)/L)^2 (x-xi)  , Ld(2*i)=(1-n(x-xi)/L)^2 (x-xi) .
for j=1:n-1,
    aux=sym(xxi,1,j);
    aux1=symop(1,'+',n,'*',aux,'/','L');aux1=sympow(aux1,2);
    aux2=symmul(aux1,aux);   Lg=sym(Lg,1,2*j,aux2);
    aux2=symop(1,'-',2,'*',n,'*',aux,'/','L');aux1=symmul(aux1,aux2);
     Lg=sym(Lg,1,2*j-1,aux1);
    aux1=symop(1,'-',n,'*',aux,'/','L');aux1=sympow(aux1,2);
    aux2=symmul(aux1,aux);   Ld=sym(Ld,1,2*j,aux2);
    aux2=symop(1,'+',2,'*',n,'*',aux,'/','L');aux1=symmul(aux1,aux2);
     Ld=sym(Ld,1,2*j-1,aux1);
end;
% cas du dernier point
    aux=sym(xxi,1,n);
    aux1=symop(1,'+',n,'*',aux,'/','L');aux1=sympow(aux1,2);
                             Lg=sym(Lg,1,2*n-1,aux1);
% a droite de l'avant dernier point:
    aux1=sym(xxi,1,n-1);
    aux2=symop(1,'-',n*n,'*',aux1,'*',aux1,'/','L^2');
                             Ld=sym(Ld,1,2*n-3,aux2);
    aux2=symop(0,'-',n,'*',aux,'*',aux1,'/','L');
                             Ld=sym(Ld,1,2*n-2,aux2);
  Lg=simplify(Lg),       Ld=simplify(Ld),
%
%  matrice de rigidite
%
A=zeros(2*n-1,2*n-1);A=sym(A);
%
for k=1:2*n-1,
    ii=floor((k+1)/2);   % point associe a Lk
     xii0=sym(xi,1,ii);
     if ii>1,xim1=sym(xi,1,ii-1);else,xim1=0;end;
     if ii<n,xip1=sym(xi,1,ii+1);else,xip1=sym(xi,1,n);end;
    auxgk=sym(Lg,1,k);   auxdk=sym(Ld,1,k);
    auxgk2=diff(auxgk,'x',2); auxdk2=diff(auxdk,'x',2);
 for p=1:2*n-1;   % recherche des Lp appropries
   %  ne pas faire comme ca si n est grand!!!
 matel=0;
  % A(k,p)=int de xi(ii-1) a xi(ii+1)  de Lk''  Lp''
    pp=floor((p+1)/2);   % point associe a Lp
     xpp0=sym(xi,1,pp);
     if pp>1,xpm1=sym(xi,1,pp-1);else,xpm1=0;end;
     if pp<n,xpp1=sym(xi,1,pp+1);else,xpp1=sym(xi,1,n);end;
    auxgp=sym(Lg,1,p);   auxdp=sym(Ld,1,p);
    auxgp2=diff(auxgp,'x',2); auxdp2=diff(auxdp,'x',2);
     if xpp1==xii0,
      matel=symadd(matel,int( symmul(auxgk2,auxdp2),'x',xim1,xii0) );
     end;
     if xpp0==xii0,
      matel=symadd(matel,int( symmul(auxgk2,auxgp2),'x',xim1,xii0) );
      matel=symadd(matel,int( symmul(auxdk2,auxdp2),'x',xii0,xip1) );
     end;
     if xpm1==xii0,
      matel=symadd(matel,int( symmul(auxdk2,auxgp2),'x',xii0,xip1) );
     end;
   A=sym(A,k,p,matel);
 end;
end;
A,
% second membre
phi=zeros(1,2*n-1);phi=sym(phi);phi=sym(phi,1,2*n-1,'F/(EI)');
% solution
u=symdiv(phi,A),
J=symop('-1/2','*',sym(u,1,2*n-1)),
