% coordinates of the 3 nodes in the initial state
% node 1
xk1_deb = x1(TRI(k,1)); yk1_deb = y1(TRI(k,1)); zk1_deb = z1(TRI(k,1));
% node 2
xk2_deb = x1(TRI(k,2)); yk2_deb = y1(TRI(k,2)); zk2_deb = z1(TRI(k,2));
% node 3
xk3_deb = x1(TRI(k,3)); yk3_deb = y1(TRI(k,3)); zk3_deb = z1(TRI(k,3));

% coordinates of the 3 nodes in the current state
% node 1
xk1_fin = x2(TRI(k,1)); yk1_fin = y2(TRI(k,1)); zk1_fin = z2(TRI(k,1));
% node 2
xk2_fin = x2(TRI(k,2)); yk2_fin = y2(TRI(k,2)); zk2_fin = z2(TRI(k,2));
% node 3
xk3_fin = x2(TRI(k,3)); yk3_fin = y2(TRI(k,3)); zk3_fin = z2(TRI(k,3));

% coordonnées du centre de gravité état initial
xGdeb(k) = (xk1_deb+xk2_deb+xk3_deb)/3 ;
yGdeb(k) = (yk1_deb+yk2_deb+yk3_deb)/3 ;
zGdeb(k) = (zk1_deb+zk2_deb+zk3_deb)/3 ;

% coordonnées du centre de gravité état courant
xGfin(k) = (xk1_fin+xk2_fin+xk3_fin)/3 ;
yGfin(k) = (yk1_fin+yk2_fin+yk3_fin)/3 ;
zGfin(k) = (zk1_fin+zk2_fin+zk3_fin)/3 ;

%%%% partie utilisée pour calculer l'aire
% coordonnées des noeuds dans repère centré au centre de gravité
%node 1 initial state
xk1_deb = x1(TRI(k,1))-xGdeb(k); yk1_deb = y1(TRI(k,1))-yGdeb(k); zk1_deb = z1(TRI(k,1))-zGdeb(k);
%node 2 initial state
xk2_deb = x1(TRI(k,2))-xGdeb(k); yk2_deb = y1(TRI(k,2))-yGdeb(k); zk2_deb = z1(TRI(k,2))-zGdeb(k);
%node 3 initial state
xk3_deb = x1(TRI(k,3))-xGdeb(k); yk3_deb = y1(TRI(k,3))-yGdeb(k); zk3_deb = z1(TRI(k,3))-zGdeb(k);

%node 1 current state
xk1_fin = x2(TRI(k,1))-xGfin(k); yk1_fin = y2(TRI(k,1))-yGfin(k); zk1_fin = z2(TRI(k,1))-zGfin(k);
%node 2 current state
xk2_fin = x2(TRI(k,2))-xGfin(k); yk2_fin = y2(TRI(k,2))-yGfin(k); zk2_fin = z2(TRI(k,2))-zGfin(k);
%node 3 current state
xk3_fin = x2(TRI(k,3))-xGfin(k); yk3_fin = y2(TRI(k,3))-yGfin(k); zk3_fin = z2(TRI(k,3))-zGfin(k);

% coordonnées des noeuds dans la base barycentrique
%état initial
Vect1(k,:)=(vect1(TRI(k,1),:)+vect1(TRI(k,2),:)+vect1(TRI(k,3),:))/3;
Vect2(k,:)=(vect2(TRI(k,1),:)+vect2(TRI(k,2),:)+vect2(TRI(k,3),:))/3;
Vect3(k,:)=(vect3(TRI(k,1),:)+vect3(TRI(k,2),:)+vect3(TRI(k,3),:))/3;

%état courant
Vect1b(k,:)=(vect1b(TRI(k,1),:)+vect1b(TRI(k,2),:)+vect1b(TRI(k,3),:))/3;
Vect2b(k,:)=(vect2b(TRI(k,1),:)+vect2b(TRI(k,2),:)+vect2b(TRI(k,3),:))/3;
Vect3b(k,:)=(vect3b(TRI(k,1),:)+vect3b(TRI(k,2),:)+vect3b(TRI(k,3),:))/3;

%node 1
xxk1_deb = xk1_deb*Vect1(k,1)+yk1_deb*Vect1(k,2)+zk1_deb*Vect1(k,3);
xxk1_fin = xk1_fin*Vect1b(k,1)+yk1_fin*Vect1b(k,2)+zk1_fin*Vect1b(k,3);

yyk1_deb = xk1_deb*Vect2(k,1)+yk1_deb*Vect2(k,2)+zk1_deb*Vect2(k,3);
yyk1_fin = xk1_fin*Vect2b(k,1)+yk1_fin*Vect2b(k,2)+zk1_fin*Vect2b(k,3);

zzk1_deb = xk1_deb*Vect3(k,1)+yk1_deb*Vect3(k,2)+zk1_deb*Vect3(k,3);
zzk1_fin = xk1_fin*Vect3b(k,1)+yk1_fin*Vect3b(k,2)+zk1_fin*Vect3b(k,3);


%node 2
xxk2_deb = xk2_deb*Vect1(k,1)+yk2_deb*Vect1(k,2)+zk2_deb*Vect1(k,3);
xxk2_fin = xk2_fin*Vect1b(k,1)+yk2_fin*Vect1b(k,2)+zk2_fin*Vect1b(k,3);

yyk2_deb = xk2_deb*Vect2(k,1)+yk2_deb*Vect2(k,2)+zk2_deb*Vect2(k,3);
yyk2_fin = xk2_fin*Vect2b(k,1)+yk2_fin*Vect2b(k,2)+zk2_fin*Vect2b(k,3);

zzk2_deb = xk2_deb*Vect3(k,1)+yk2_deb*Vect3(k,2)+zk2_deb*Vect3(k,3);
zzk2_fin = xk2_fin*Vect3b(k,1)+yk2_fin*Vect3b(k,2)+zk2_fin*Vect3b(k,3);


%node 3
xxk3_deb = xk3_deb*Vect1(k,1)+yk3_deb*Vect1(k,2)+zk3_deb*Vect1(k,3);
xxk3_fin = xk3_fin*Vect1b(k,1)+yk3_fin*Vect1b(k,2)+zk3_fin*Vect1b(k,3);

yyk3_deb = xk3_deb*Vect2(k,1)+yk3_deb*Vect2(k,2)+zk3_deb*Vect2(k,3);
yyk3_fin = xk3_fin*Vect2b(k,1)+yk3_fin*Vect2b(k,2)+zk3_fin*Vect2b(k,3);

zzk3_deb = xk3_deb*Vect3(k,1)+yk3_deb*Vect3(k,2)+zk3_deb*Vect3(k,3);
zzk3_fin = xk3_fin*Vect3b(k,1)+yk3_fin*Vect3b(k,2)+zk3_fin*Vect3b(k,3);

% area and volume of each element
surf_fin(k,step) = abs(det([1 xxk1_fin yyk1_fin; 1 xxk2_fin yyk2_fin; 1 xxk3_fin yyk3_fin]) / 2) ;
surf_ini(k) = abs(det([1 xxk1_deb yyk1_deb; 1 xxk2_deb yyk2_deb; 1 xxk3_deb yyk3_deb]) / 2) ;

%%%% fin partie utilisée pour calculer l'aire

% Curvilinear Coordinates
% Partie utilisé retrouver G_1, G_2, G^1, G^2, g_1, g_2, g^1, g^2 pour
% chaque element triangulaire pour retrouver F et finalement E

Delta = xk2_deb.*yk3_deb-xk3_deb.*yk2_deb-(xk1_deb.*yk3_deb-xk3_deb.*yk1_deb)+xk1_deb.*yk2_deb-xk2_deb.*yk1_deb;

% calcul de G_1
G_1(k,1)=1;
G_1(k,2)=0;
G_1(k,3)=zk1_deb*(yk2_deb-yk3_deb)/Delta+zk2_deb*(yk3_deb-yk1_deb)/Delta+zk3_deb*(yk1_deb-yk2_deb)/Delta;

% calcul de G_2
G_2(k,1)=0;
G_2(k,2)=1;
G_2(k,3)=zk1_deb*(xk2_deb-xk3_deb)/Delta+zk2_deb*(xk3_deb-xk1_deb)/Delta+zk3_deb*(xk1_deb-xk2_deb)/Delta;

% calcul de g_1
g_1(k,1)=xk1_fin*(yk2_deb-yk3_deb)/Delta+xk2_fin*(yk3_deb-yk1_deb)/Delta+xk3_fin*(yk1_deb-yk2_deb)/Delta;
g_1(k,2)=yk1_fin*(yk2_deb-yk3_deb)/Delta+yk2_fin*(yk3_deb-yk1_deb)/Delta+yk3_fin*(yk1_deb-yk2_deb)/Delta;
g_1(k,3)=zk1_fin*(yk2_deb-yk3_deb)/Delta+zk2_fin*(yk3_deb-yk1_deb)/Delta+zk3_fin*(yk1_deb-yk2_deb)/Delta;

% calcul de g_2
g_2(k,1)=xk1_fin*(xk2_deb-xk3_deb)/Delta+xk2_fin*(xk3_deb-xk1_deb)/Delta+xk3_fin*(xk1_deb-xk2_deb)/Delta;
g_2(k,2)=yk1_fin*(xk2_deb-xk3_deb)/Delta+yk2_fin*(xk3_deb-xk1_deb)/Delta+yk3_fin*(xk1_deb-xk2_deb)/Delta;
g_2(k,3)=zk1_fin*(xk2_deb-xk3_deb)/Delta+zk2_fin*(xk3_deb-xk1_deb)/Delta+zk3_fin*(xk1_deb-xk2_deb)/Delta;

% calcul de G^1
% il faut le vecteur normal initial N = vect3
N=Vect3(k,:); %N=(0,0,1);
intermG1(k,1)=G_2(k,2)*N(3)-G_2(k,3)*N(2);
intermG1(k,2)=G_2(k,3)*N(1)-G_2(k,1)*N(3);
intermG1(k,3)=G_2(k,1)*N(2)-G_2(k,2)*N(1);
normintermG1(k) = G_1(k,1)*intermG1(k,1)+G_1(k,2)*intermG1(k,2)+G_1(k,3)*intermG1(k,3);
G1(k,:)=intermG1(k,:)/normintermG1(k);

% calcul de G^2
% il faut le vecteur normal initial N = vect3
N=Vect3(k,:); %N=(0,0,1);
intermG2(k,1)=G_1(k,2)*N(3)-G_1(k,3)*N(2);
intermG2(k,2)=G_1(k,3)*N(1)-G_1(k,1)*N(3);
intermG2(k,3)=G_1(k,1)*N(2)-G_1(k,2)*N(1);
normintermG2(k) = G_2(k,1)*intermG2(k,1)+G_2(k,2)*intermG2(k,2)+G_2(k,3)*intermG2(k,3);
G2(k,:)=-intermG2(k,:)/normintermG2(k);

% calcul de g^1
% il faut le vecteur normal courant n = vect3b
n=Vect3b(k,:); %calculé avec la base barycentrique
intermg1(k,1)=g_2(k,2)*n(3)-g_2(k,3)*n(2);
intermg1(k,2)=g_2(k,3)*n(1)-g_2(k,1)*n(3);
intermg1(k,3)=g_2(k,1)*n(2)-g_2(k,2)*n(1);
normintermg1(k) = g_1(k,1)*intermg1(k,1)+g_1(k,2)*intermg1(k,2)+g_1(k,3)*intermg1(k,3);
g1(k,:)=intermg1(k,:)/normintermg1(k);

% calcul de g^2
% il faut le vecteur normal courant n = vect3b
n=Vect3b(k,:); %calculé avec la base barycentrique
intermg2(k,1)=g_1(k,2)*n(3)-g_1(k,3)*n(2);
intermg2(k,2)=g_1(k,3)*n(1)-g_1(k,1)*n(3);
intermg2(k,3)=g_1(k,1)*n(2)-g_1(k,2)*n(1);
normintermg2(k) = g_2(k,1)*intermg2(k,1)+g_2(k,2)*intermg2(k,2)+g_2(k,3)*intermg2(k,3);
g2(k,:)=-intermg2(k,:)/normintermg2(k);

%Composantes de F
F11 = g_1(k,1)*G1(k,1)+g_2(k,1)*G2(k,1) ;
F12 = g_1(k,1)*G1(k,2)+g_2(k,1)*G2(k,2) ;
F21 = g_1(k,2)*G1(k,1)+g_2(k,2)*G2(k,1) ;
F22 = g_1(k,2)*G1(k,2)+g_2(k,2)*G2(k,2) ;
F13 = g_1(k,1)*G1(k,3)+g_2(k,1)*G2(k,3) ;
F23 = g_1(k,2)*G1(k,3)+g_2(k,2)*G2(k,3) ;
F31 = g_1(k,3)*G1(k,1)+g_2(k,3)*G2(k,1) ;
F32 = g_1(k,3)*G1(k,2)+g_2(k,3)*G2(k,2) ;
F33 = g_1(k,3)*G1(k,3)+g_2(k,3)*G2(k,3) ;

%Construction du tenseur F
F_tensor = [ F11 F12 F13 ; F21 F22 F23 ; F31 F32 F33 ] ;

F11_tensor(k,:)=F_tensor(1,1);
F22_tensor(k,:)=F_tensor(2,2);
F12_tensor(k,:)=F_tensor(1,2);
F21_tensor(k,:)=F_tensor(2,1);


DEF_GRADIENT(k,step)=F11_tensor(k);
DEF_GRADIENT(Nbel+k,step)=F22_tensor(k);
DEF_GRADIENT(2*Nbel+k,step)=F12_tensor(k);
DEF_GRADIENT(3*Nbel+k,step)=F21_tensor(k);