clear all
close all
clc

load data_TP

Nstep=size(x,2);
Nbel=size(TRI,1);

hmesh=0.2;
x1=x_0;y1=y_0;z1=z_0;

initial_normals ; % this function derives the normal vector vect_n at each node in the initial state

% local orthonormal basis in the initial configuration
% vect1 and vect2 are the vector basis in the plane and vect3 is the normal vector
% vect2 est aligné avec la direction y du plan du capteur CCD
vect10 = [ones(length(x1),1),zeros(length(x1),2)] ;
vect20 = [zeros(length(x1),1),ones(length(x1),1),zeros(length(x1),1)] ;
vect3 = vect_n ;

%produit scalaire de vect20 et vect3
dot23=vect20.*vect3;
dot23=sum(dot23,2)*ones(1,3);

vect22 = vect20-dot23.*vect3;
norm_vect_2 = sqrt(vect22(:,1).^2+vect22(:,2).^2+vect22(:,3).^2) ;
vect2 = vect22./(norm_vect_2*ones(1,3)) ;
vect1 = [vect2(:,2).*vect3(:,3)-vect2(:,3).*vect3(:,2) , vect2(:,3).*vect3(:,1)-vect2(:,1).*vect3(:,3) , ...
    vect2(:,1).*vect3(:,2)-vect2(:,2).*vect3(:,1) ] ;

for step=1:Nstep, %1:Nstep
    
    x2=x(:,step);y2=y(:,step);z2=z(:,step);
    
    final_normals ; % this function derives the normal vector vect_n at each node in the current state
    
    % local basis in the final configuration
    % vect1b and vect2b are the vector basis in the plane and vect3b is the normalvector
    vect10b = vect1 ;
    vect20b = vect2 ;
    vect3b = vect_n ;
    
    %produit scalaire de vect20 et vect3
    dot23b=vect20b.*vect3b;
    dot23b=sum(dot23b,2)*ones(1,3);
    
    vect22b = vect20b-dot23b.*vect3b;
    norm_vect_2b = sqrt(vect22b(:,1).^2+vect22b(:,2).^2+vect22b(:,3).^2) ;
    vect2b = vect22b./(norm_vect_2b*ones(1,3)) ;
    vect1b = [vect2b(:,2).*vect3b(:,3)-vect2b(:,3).*vect3b(:,2) , vect2b(:,3).*vect3b(:,1)-vect2b(:,1).*vect3b(:,3) , ...
        vect2b(:,1).*vect3b(:,2)-vect2b(:,2).*vect3b(:,1) ] ;
    
    for k =1:Nbel, %1:Nbel,
        
        Mydifferentiate
        
        %Tenseur E
        E_tensor = (F_tensor'*F_tensor-eye(3))/2 ;
        
        %Composantes de E
        E11(k,step)=E_tensor(1,1);
        E22(k,step)=E_tensor(2,2);
        E12(k,step)=E_tensor(1,2);
        
        EPSILON(k,step)=E11(k,step);
        EPSILON(Nbel+k,step)=E22(k,step);
        EPSILON(2*Nbel+k,step)=E12(k,step);
        
        
        
    end
    
    Myplot(xGfin,yGfin,zGfin,E11(:,step))
    
end