%%
clear all
close all
clc

elm_nodes=4;
bdf='C:\Users\Nicco\Documents\tesi\Modello_beam\largo_ifem_beams_thk_update_alu.bdf';

str_coord='ijk';

[nEle,nNode,Cne,index_quad,index_tria,t,Pnn,gcoord_index,thk,Off,bc,Cbar,...
                    Off_beam,Pbarl,Pbarl_index,mat_e,gcoord,mcoord]=iFem_pre_v2(bdf,str_coord);
% Questa funzione importa il bdf, Cne è connectivity, Pnn è coordinate nodi
% mcoord importante, sono 60 matrici, una per ogni elemento, con versore di
% asse locale per ogni elemento, per righe. Serve poi in Rot a trasformare coordinate
% nodi da locale a globale.
q = zeros(nNode*6,1);% vettore dei gradi di libertà nodali, riferimento globale, 6 DOFs per ogni nodo (84 nodi)
A = load('C:\Users\Nicco\Documents\tesi\Codici\U_beam_red.txt');
B = load('C:\Users\Nicco\Documents\tesi\Codici\UR_beam_red.txt');
q(1:6:end)=A(:,2); % In q ho 6 DOFs nodo 1, 6 DOFs nodo 2, ecc
q(2:6:end)=A(:,3);
q(3:6:end)=A(:,4);
q(4:6:end)=B(:,2);
q(5:6:end)=B(:,3);
q(6:6:end)=B(:,4);
straing = zeros(nEle*6*9,1);
strainconzeri = zeros(nEle*6*5,1);
matrice = zeros(nEle*6*9,nNode*6);
Rotespansa = zeros(24,24);
%qtot=zeros(6*4*nEle,1);
%qvero=zeros(nNode*6,1);
 for e=1:nEle
    %DEFINING THE GLOBAL COORDINATES OF THE NODES AND THE ROTATION MATIX
    % FOR EACH ELEMENT
    glocoord = Pnn(gcoord_index(Cne(e,2:elm_nodes+1)),:);
    Rot = mcoord(:,:,e);
    % Rot per trasformare coordinate
    
    % DEFINING THE TOP (Z2) AND BOTTOM (Z1) SURFACES OF THE ELEMENT
    z2=t(e)+Off(e);
    z1=-t(e)+Off(e);

    % TRANSFORMING THE GLOBAL COORDINATES TO THE LOCAL ONES
    locCoord = glocoord * Rot';

    X = locCoord(:,1);
    Y = locCoord(:,2);

    [Bm(:,:,e),Bk(:,:,e),Bg(:,:,e)]=compute_B_quad(X,Y,2,3,0,0);
    % Derivate funzioni di forma nei punti di Gauss (membranali, flessione,
    % taglio che non mi servono)
    % X e Y coordinate dei nodi nel sistema locale, 2 è flag per dire come
    % lavorare, ovvero calcola nei punti di Gauss, 3 punti di Gauss, che in
    % realtà è 3*3, 0 e 0 per flag = 1.
    for i=1:8
    Rotespansa((i-1)*3+1:(i-1)*3+3,(i-1)*3+1:(i-1)*3+3) = Rot;
    % E' una matrice 24*24 con 8 blocchi 3*3 in diagonale che hanno il
    % vettore di orientamento del SR
    end
    nodie = Cne(e,2:5); % Per ogni elemento salvo i 4 nodi della connectivity, corrispondenti alle colonne della matrice. 
    matricettasu = Bm(:,:,e)+z2*Bk(:,:,e); %Per ogni elemento salvo la matrice che deve moltiplicare i DOFs
    matricettagiu = Bm(:,:,e)+z1*Bk(:,:,e);%stessa cosa ma a -spessore
    matricettasu  = matricettasu*Rotespansa;
    matricettagiu  = matricettagiu*Rotespansa;
    k = 1;
    while k < 5 % Ovvero 4 nodi per ogni elemento
    % matrice salva per ogni 6 colonne i gradi di libertà di un nodo,
    % Per esempio, il nodo 7 avrà 6 DOFs dalle colonne 37 a 42 di
    % matrice, corrispondenti alle righe da 37 a 42 di q. 
    matrice(1+(e-1)*54:1+(e-1)*54+26,1+(nodie(k)-1)*6:1+(nodie(k)-1)*6+5) = matricettasu(:,1+(k-1)*6:1+(k-1)*6+5);
    matrice(1+(e-1)*54+27:1+(e-1)*54+53,1+(nodie(k)-1)*6:1+(nodie(k)-1)*6+5) = matricettagiu(:,1+(k-1)*6:1+(k-1)*6+5);
    k=k+1;
    end

% In questo momento le epsilon sono (9*3=27) per ogni punto di Gauss in un
% elemento e *2 (top e bottom), quindi 54 per ogni elemento. Il vettore
% qtot ha (6*4=24) ovvero 6DOFs per nodo per ogni elemento (*nEle). Quindi
% la matrice sarà una (54*24) espansa per nEle. 

end
straing = matrice*q;
% straing ha primi 9 elementi= epsxx(1 elemento) nei 9 punti di Gauss,
% secondi 9 elementi = epsyy(1 elemento) ecc... 
% Sono 60(nEle) pacchetti da (9*3*2)=54, ogni pacchetto ha prima 27 top poi
% 27 bottom per ogni elemento

% Ora devo trovare gli strain nei nodi dai punti di Gauss
 %Scrivo le funzioni di forma per l'elemento interno che ha come nodi i
 %punti di Gauss
% coordinate elemento nel SR locale
 %zie = [-1 1 0 1 -1];
 %etae = [-1 -1 0 1 1];
 % trasformazione coordinate elemento di Gauss
 %zig = sqrt(3)*zie;
 %etag = sqrt(3)*etae;
 % Voglio gli indici dei nodi di Gauss secondo la numerazione che corrispon
 % dono ai nodi dell'elemento, in senso antiorario

% Gli strain nei centroidi sono il quinto elemento ogni 9 di straing
% Così ottengo straincentoidi(1:6) = 3 epsxx,yy,xy top e 3 bottom per il
% primo elemento, e così via
straincentroidi = zeros(6*nEle,1);
for e=1:6*nEle
straincentroidi(e) = straing(9*(e-1)+5);
end
% Voglio sotto forma di matrice gli strain nei centroidi, voglio 3 colonne
% xx,yy,xy e tutti gli elementi da 1:nEle ripetuti due volte, prima top e
% poi bottom
straincentroidimat = zeros(nEle*2,3);
for e=1:nEle
straincentroidimat(e,:) = straincentroidi(6*(e-1)+1:6*(e-1)+3);
straincentroidimat(nEle+e,:) = straincentroidi(6*(e-1)+4:6*(e-1)+6);
end

% errori tra strains nei centroidi
% carico il file con gli strain da laboratorio
strain_lab = load('C:\Users\Nicco\Documents\tesi\Modello_beam\Cantele_shell_beams_shell.rpt');
% Quali elementi sono sensorizzati?
index_1 = [3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47];%c'è
index_2 = [50:72];
index_3 = [121:144];
index_4 = [194:216];
index_5 = [265:288];%c'è
index_6 = [337:360];%c'è
index_7 = [385:408];
index_8 = [457:480];
index_9 = [529:552];
index_10 = [554 556 558 560 562 564 566 568 570 572 574 576 578 580 582 584 586 588 590 592 594 596 598 600];
index_11 = [602:624];%c'è
%index_12 = [673:696];
%index_13 = [745:768];
%index_14 = [817:840];
%index_15 = [889:912];
index_16 = [242:264];
% Fibra 1
errori_xx1 = zeros(length(index_1),1);
%errori_yy = zeros(nEle*2,1);
%errori_xy = zeros(nEle*2,1);
% Calcolo l'errore degli epsxx solo negli elementi sensorizzati
straincentroidimatsensor1 = zeros(length(index_1),1);
strain_lab1 = zeros(length(index_1),1);
for i=1:length(index_1)
e = index_1(i);
straincentroidimatsensor1(i) = straincentroidimat(e,1);
strain_lab1(i) = strain_lab(e,2);
errori_xx1(i) = abs(straincentroidimat(e,1)-strain_lab(e,2))/abs(strain_lab(e,2))*100;
end
figure(1)
plot(index_1,straincentroidimatsensor1*1e6,'-or')
hold on
plot(index_1, strain_lab1*1e6,'-ok')
xlabel('Elements along the Fiber')
ylabel('Longitudinal micro-Strain')
legend('Calculated micro-Strain','Experimental micro-Strain')
title('Fiber 1')
figure(2)
plot(index_1, errori_xx1, '-ob')
xlabel('Elements along the Fiber')
ylabel('Error [%]')
title('Fiber 1')
% Fibra 2
errori_xx2 = zeros(length(index_2),1);
%errori_yy = zeros(nEle*2,1);
%errori_xy = zeros(nEle*2,1);
% Calcolo l'errore degli epsxx solo negli elementi sensorizzati
straincentroidimatsensor2 = zeros(length(index_2),1);
strain_lab2 = zeros(length(index_2),1);
for i=1:length(index_2)
e = index_2(i);
straincentroidimatsensor2(i) = straincentroidimat(e,1);
strain_lab2(i) = strain_lab(e,2);
errori_xx2(i) = abs(straincentroidimat(e,1)-strain_lab(e,2))/abs(strain_lab(e,2))*100;
end
figure(3)
plot(index_2,straincentroidimatsensor2*1e6,'-or')
hold on
plot(index_2, strain_lab2*1e6,'-ok')
xlabel('Elements along the Fiber')
ylabel('Longitudinal micro-Strain')
legend('Calculated micro-Strain','Experimental micro-Strain')
title('Fiber 2')
figure(4)
plot(index_2, errori_xx2, '-ob')
xlabel('Elements along the Fiber')
ylabel('Error [%]')
title('Fiber 2')
% Fibra 3
errori_xx3 = zeros(length(index_3),1);
%errori_yy = zeros(nEle*2,1);
%errori_xy = zeros(nEle*2,1);
% Calcolo l'errore degli epsxx solo negli elementi sensorizzati
straincentroidimatsensor3 = zeros(length(index_3),1);
strain_lab3 = zeros(length(index_3),1);
for i=1:length(index_3)
e = index_3(i);
straincentroidimatsensor3(i) = straincentroidimat(e,1);
strain_lab3(i) = strain_lab(e,2);
errori_xx3(i) = abs(straincentroidimat(e,1)-strain_lab(e,2))/abs(strain_lab(e,2))*100;
end
figure(5)
plot(index_3,straincentroidimatsensor3*1e6,'-or')
hold on
plot(index_3, strain_lab3*1e6,'-ok')
xlabel('Elements along the Fiber')
ylabel('Longitudinal micro-Strain')
legend('Calculated micro-Strain','Experimental micro-Strain')
title('Fiber 3')
figure(6)
plot(index_3, errori_xx3, '-ob')
xlabel('Elements along the Fiber')
ylabel('Error [%]')
title('Fiber 3')
% Fibra 4
errori_xx4 = zeros(length(index_4),1);
%errori_yy = zeros(nEle*2,1);
%errori_xy = zeros(nEle*2,1);
% Calcolo l'errore degli epsxx solo negli elementi sensorizzati
straincentroidimatsensor4 = zeros(length(index_4),1);
strain_lab4 = zeros(length(index_4),1);
for i=1:length(index_4)
e = index_4(i);
straincentroidimatsensor4(i) = straincentroidimat(e,1);
strain_lab4(i) = strain_lab(e,2);
errori_xx4(i) = abs(straincentroidimat(e,1)-strain_lab(e,2))/abs(strain_lab(e,2))*100;
end
figure(7)
plot(index_4,straincentroidimatsensor4*1e6,'-or')
hold on
plot(index_4, strain_lab4*1e6,'-ok')
xlabel('Elements along the Fiber')
ylabel('Longitudinal micro-Strain')
legend('Calculated micro-Strain','Experimental micro-Strain')
title('Fiber 4')
figure(8)
plot(index_4, errori_xx4, '-ob')
xlabel('Elements along the Fiber')
ylabel('Error [%]')
title('Fiber 4')
% Fibra 5
errori_xx5 = zeros(length(index_5),1);
%errori_yy = zeros(nEle*2,1);
%errori_xy = zeros(nEle*2,1);
% Calcolo l'errore degli epsxx solo negli elementi sensorizzati
straincentroidimatsensor5 = zeros(length(index_5),1);
strain_lab5 = zeros(length(index_5),1);
for i=1:length(index_5)
e = index_5(i);
straincentroidimatsensor5(i) = straincentroidimat(e,1);
strain_lab5(i) = strain_lab(e,2);
errori_xx5(i) = abs(straincentroidimat(e,1)-strain_lab(e,2))/abs(strain_lab(e,2))*100;
end
figure(9)
plot(index_5,straincentroidimatsensor5*1e6,'-or')
hold on
plot(index_5, strain_lab5*1e6,'-ok')
xlabel('Elements along the Fiber')
ylabel('Longitudinal micro-Strain')
legend('Calculated micro-Strain','Experimental micro-Strain')
title('Fiber 5')
figure(10)
plot(index_5, errori_xx5, '-ob')
xlabel('Elements along the Fiber')
ylabel('Error [%]')
title('Fiber 5')
% Fibra 6
errori_xx6 = zeros(length(index_6),1);
%errori_yy = zeros(nEle*2,1);
%errori_xy = zeros(nEle*2,1);
% Calcolo l'errore degli epsxx solo negli elementi sensorizzati
straincentroidimatsensor6 = zeros(length(index_6),1);
strain_lab6 = zeros(length(index_6),1);
for i=1:length(index_6)
e = index_6(i);
straincentroidimatsensor6(i) = straincentroidimat(e,1);
strain_lab6(i) = strain_lab(e,2);
errori_xx6(i) = abs(straincentroidimat(e,1)-strain_lab(e,2))/abs(strain_lab(e,2))*100;
end
figure(11)
plot(index_6,straincentroidimatsensor6*1e6,'-or')
hold on
plot(index_6, strain_lab6*1e6,'-ok')
xlabel('Elements along the Fiber')
ylabel('Longitudinal micro-Strain')
legend('Calculated micro-Strain','Experimental micro-Strain')
title('Fiber 6')
figure(12)
plot(index_6, errori_xx6, '-ob')
xlabel('Elements along the Fiber')
ylabel('Error [%]')
title('Fiber 6')
% Fibra 7
errori_xx7 = zeros(length(index_7),1);
%errori_yy = zeros(nEle*2,1);
%errori_xy = zeros(nEle*2,1);
% Calcolo l'errore degli epsxx solo negli elementi sensorizzati
straincentroidimatsensor7 = zeros(length(index_7),1);
strain_lab7 = zeros(length(index_7),1);
for i=1:length(index_7)
e = index_7(i);
straincentroidimatsensor7(i) = straincentroidimat(e,1);
strain_lab7(i) = strain_lab(e,2);
errori_xx7(i) = abs(straincentroidimat(e,1)-strain_lab(e,2))/abs(strain_lab(e,2))*100;
end
figure(13)
plot(index_7,straincentroidimatsensor7*1e6,'-or')
hold on
plot(index_7, strain_lab7*1e6,'-ok')
xlabel('Elements along the Fiber')
ylabel('Longitudinal micro-Strain')
legend('Calculated micro-Strain','Experimental micro-Strain')
title('Fiber 7')
figure(14)
plot(index_7, errori_xx7, '-ob')
xlabel('Elements along the Fiber')
ylabel('Error [%]')
title('Fiber 7')
% Fibra 8
errori_xx8 = zeros(length(index_8),1);
%errori_yy = zeros(nEle*2,1);
%errori_xy = zeros(nEle*2,1);
% Calcolo l'errore degli epsxx solo negli elementi sensorizzati
straincentroidimatsensor8 = zeros(length(index_8),1);
strain_lab8 = zeros(length(index_8),1);
for i=1:length(index_8)
e = index_8(i);
straincentroidimatsensor8(i) = straincentroidimat(e,1);
strain_lab8(i) = strain_lab(e,2);
errori_xx8(i) = abs(straincentroidimat(e,1)-strain_lab(e,2))/abs(strain_lab(e,2))*100;
end
figure(15)
plot(index_8,straincentroidimatsensor8*1e6,'-or')
hold on
plot(index_8, strain_lab8*1e6,'-ok')
xlabel('Elements along the Fiber')
ylabel('Longitudinal micro-Strain')
legend('Calculated micro-Strain','Experimental micro-Strain')
title('Fiber 8')
figure(16)
plot(index_8, errori_xx8, '-ob')
xlabel('Elements along the Fiber')
ylabel('Error [%]')
title('Fiber 8')
% Fibra 9
errori_xx9 = zeros(length(index_9),1);
%errori_yy = zeros(nEle*2,1);
%errori_xy = zeros(nEle*2,1);
% Calcolo l'errore degli epsxx solo negli elementi sensorizzati
straincentroidimatsensor9 = zeros(length(index_9),1);
strain_lab9 = zeros(length(index_9),1);
for i=1:length(index_9)
e = index_9(i);
straincentroidimatsensor9(i) = straincentroidimat(e,1);
strain_lab9(i) = strain_lab(e,2);
errori_xx9(i) = abs(straincentroidimat(e,1)-strain_lab(e,2))/abs(strain_lab(e,2))*100;
end
figure(17)
plot(index_9,straincentroidimatsensor9*1e6,'-or')
hold on
plot(index_9, strain_lab9*1e6,'-ok')
xlabel('Elements along the Fiber')
ylabel('Longitudinal micro-Strain')
legend('Calculated micro-Strain','Experimental micro-Strain')
title('Fiber 9')
figure(18)
plot(index_9, errori_xx9, '-ob')
xlabel('Elements along the Fiber')
ylabel('Error [%]')
title('Fiber 9')
% Fibra 10
errori_xx10 = zeros(length(index_10),1);
%errori_yy = zeros(nEle*2,1);
%errori_xy = zeros(nEle*2,1);
% Calcolo l'errore degli epsxx solo negli elementi sensorizzati
straincentroidimatsensor10 = zeros(length(index_10),1);
strain_lab10 = zeros(length(index_10),1);
for i=1:length(index_10)
e = index_10(i);
straincentroidimatsensor10(i) = straincentroidimat(e,1);
strain_lab10(i) = strain_lab(e,2);
errori_xx10(i) = abs(straincentroidimat(e,1)-strain_lab(e,2))/abs(strain_lab(e,2))*100;
end
figure(19)
plot(index_10,straincentroidimatsensor10*1e6,'-or')
hold on
plot(index_10, strain_lab10*1e6,'-ok')
xlabel('Elements along the Fiber')
ylabel('Longitudinal micro-Strain')
legend('Calculated micro-Strain','Experimental micro-Strain')
title('Fiber 10')
figure(20)
plot(index_10, errori_xx10, '-ob')
xlabel('Elements along the Fiber')
ylabel('Error [%]')
title('Fiber 10')
% Fibra 11
errori_xx11 = zeros(length(index_11),1);
%errori_yy = zeros(nEle*2,1);
%errori_xy = zeros(nEle*2,1);
% Calcolo l'errore degli epsxx solo negli elementi sensorizzati
straincentroidimatsensor11 = zeros(length(index_11),1);
strain_lab11 = zeros(length(index_11),1);
for i=1:length(index_11)
e = index_11(i);
straincentroidimatsensor11(i) = straincentroidimat(e,1);
strain_lab11(i) = strain_lab(e,2);
errori_xx11(i) = abs(straincentroidimat(e,1)-strain_lab(e,2))/abs(strain_lab(e,2))*100;
end
figure(21)
plot(index_11,straincentroidimatsensor11*1e6,'-or')
hold on
plot(index_11, strain_lab11*1e6,'-ok')
xlabel('Elements along the Fiber')
ylabel('Longitudinal micro-Strain')
legend('Calculated micro-Strain','Experimental micro-Strain')
title('Fiber 11')
figure(22)
plot(index_11, errori_xx11, '-ob')
xlabel('Elements along the Fiber')
ylabel('Error [%]')
title('Fiber 11')
% Fibra 16
errori_xx16 = zeros(length(index_16),1);
%errori_yy = zeros(nEle*2,1);
%errori_xy = zeros(nEle*2,1);
% Calcolo l'errore degli epsxx solo negli elementi sensorizzati
straincentroidimatsensor16 = zeros(length(index_16),1);
strain_lab16 = zeros(length(index_16),1);
for i=1:length(index_16)
e = index_16(i);
straincentroidimatsensor16(i) = straincentroidimat(e,1);
strain_lab16(i) = strain_lab(e,2);
errori_xx16(i) = abs(straincentroidimat(e,1)-strain_lab(e,2))/abs(strain_lab(e,2))*100;
end
figure(31)
plot(index_16,straincentroidimatsensor16*1e6,'-or')
hold on
plot(index_16, strain_lab16*1e6,'-ok')
xlabel('Elements along the Fiber')
ylabel('Longitudinal micro-Strain')
legend('Calculated micro-Strain','Experimental micro-Strain')
title('Fiber 16')
figure(32)
plot(index_16, errori_xx16, '-ob')
xlabel('Elements along the Fiber')
ylabel('Error [%]')
title('Fiber 16')
%%
% Fibra 12
errori_xx12 = zeros(length(index_12),1);
%errori_yy = zeros(nEle*2,1);
%errori_xy = zeros(nEle*2,1);
% Calcolo l'errore degli epsxx solo negli elementi sensorizzati
straincentroidimatsensor12 = zeros(length(index_12),1);
strain_lab12 = zeros(length(index_12),1);
for i=1:length(index_12)
e = index_12(i);
straincentroidimatsensor12(i) = straincentroidimat(e,1);
strain_lab12(i) = strain_lab(e,2);
errori_xx12(i) = abs(straincentroidimat(e,1)-strain_lab(e,2))/abs(strain_lab(e,2))*100;
end
figure(23)
plot(index_12,straincentroidimatsensor12*1e6,'-or')
hold on
plot(index_12, strain_lab12*1e6,'-ok')
xlabel('Elements along the Fiber')
ylabel('Longitudinal micro-Strain')
legend('Calculated micro-Strain','Experimental micro-Strain')
title('Fiber 12')
figure(24)
plot(index_12, errori_xx12, '-ob')
xlabel('Elements along the Fiber')
ylabel('Error [%]')
title('Fiber 12')
% Fibra 13
errori_xx13 = zeros(length(index_13),1);
%errori_yy = zeros(nEle*2,1);
%errori_xy = zeros(nEle*2,1);
% Calcolo l'errore degli epsxx solo negli elementi sensorizzati
straincentroidimatsensor13 = zeros(length(index_13),1);
strain_lab13 = zeros(length(index_13),1);
for i=1:length(index_13)
e = index_13(i);
straincentroidimatsensor13(i) = straincentroidimat(e,1);
strain_lab13(i) = strain_lab(e,2);
errori_xx13(i) = abs(straincentroidimat(e,1)-strain_lab(e,2))/abs(strain_lab(e,2))*100;
end
figure(25)
plot(index_13,straincentroidimatsensor13*1e6,'-or')
hold on
plot(index_13, strain_lab13*1e6,'-ok')
xlabel('Elements along the Fiber')
ylabel('Longitudinal micro-Strain')
legend('Calculated micro-Strain','Experimental micro-Strain')
title('Fiber 13')
figure(26)
plot(index_13, errori_xx13, '-ob')
xlabel('Elements along the Fiber')
ylabel('Error [%]')
title('Fiber 13')
%% Fibra 14
errori_xx14 = zeros(length(index_14),1);
%errori_yy = zeros(nEle*2,1);
%errori_xy = zeros(nEle*2,1);
% Calcolo l'errore degli epsxx solo negli elementi sensorizzati
straincentroidimatsensor14 = zeros(length(index_14),1);
strain_lab14 = zeros(length(index_14),1);
for i=1:length(index_14)
e = index_14(i);
straincentroidimatsensor14(i) = straincentroidimat(e,1);
strain_lab14(i) = strain_lab(e,2);
errori_xx14(i) = abs(straincentroidimat(e,1)-strain_lab(e,2))/abs(strain_lab(e,2))*100;
end
figure(27)
plot(index_14,straincentroidimatsensor14*1e6,'-or')
hold on
plot(index_14, strain_lab14*1e6,'-ok')
xlabel('Elements along the Fiber')
ylabel('Longitudinal micro-Strain')
legend('Calculated micro-Strain','Experimental micro-Strain')
title('Fiber 14')
figure(28)
plot(index_14, errori_xx14, '-ob')
xlabel('Elements along the Fiber')
ylabel('Error [%]')
title('Fiber 14')
% Fibra 15
errori_xx15 = zeros(length(index_15),1);
%errori_yy = zeros(nEle*2,1);
%errori_xy = zeros(nEle*2,1);
% Calcolo l'errore degli epsxx solo negli elementi sensorizzati
straincentroidimatsensor15 = zeros(length(index_15),1);
strain_lab15 = zeros(length(index_15),1);
for i=1:length(index_15)
e = index_15(i);
straincentroidimatsensor15(i) = straincentroidimat(e,1);
strain_lab15(i) = strain_lab(e,2);
errori_xx15(i) = abs(straincentroidimat(e,1)-strain_lab(e,2))/abs(strain_lab(e,2))*100;
end
figure(29)
plot(index_15,straincentroidimatsensor15*1e6,'-or')
hold on
plot(index_15, strain_lab15*1e6,'-ok')
xlabel('Elements along the Fiber')
ylabel('Longitudinal micro-Strain')
legend('Calculated micro-Strain','Experimental micro-Strain')
title('Fiber 15')
figure(30)
plot(index_15, errori_xx15, '-ob')
xlabel('Elements along the Fiber')
ylabel('Error [%]')
title('Fiber 15')
% Fibra 16
errori_xx16 = zeros(length(index_16),1);
%errori_yy = zeros(nEle*2,1);
%errori_xy = zeros(nEle*2,1);
% Calcolo l'errore degli epsxx solo negli elementi sensorizzati
straincentroidimatsensor16 = zeros(length(index_16),1);
strain_lab16 = zeros(length(index_16),1);
for i=1:length(index_16)
e = index_16(i);
straincentroidimatsensor16(i) = straincentroidimat(e,1);
strain_lab16(i) = strain_lab(e,2);
errori_xx16(i) = abs(straincentroidimat(e,1)-strain_lab(e,2))/abs(strain_lab(e,2))*100;
end
figure(31)
plot(index_16,straincentroidimatsensor16*1e6,'-or')
hold on
plot(index_16, strain_lab16*1e6,'-ok')
xlabel('Elements along the Fiber')
ylabel('Longitudinal micro-Strain')
legend('Calculated micro-Strain','Experimental micro-Strain')
figure(32)
plot(index_16, errori_xx16, '-ob')
xlabel('Elements along the Fiber')
ylabel('Error [%]')
%%
% Matrice N è una 4*9, 4 nodi ottenuti da 9 Gauss
% Le coordinate dei nodi nel SR dei punti di Gauss sono sqrt(5/3), più grandi
% di 1 poiché 1 sono i punti di Gauss. 
chin = [-sqrt(5/3) sqrt(5/3) sqrt(5/3) -sqrt(5/3)];
etan = [-sqrt(5/3) -sqrt(5/3) sqrt(5/3) sqrt(5/3)];
% Ci sono 9 funzioni di forma, N1 deve moltiplicare gli strain del primo
% punto di Gauss e così via, mettendo nelle 4 righe le coordinate dei 4
% nodi nel SR di Gauss
for i = 1:4
 N(i,1) = 1/4*(chin(i)-1)*(etan(i)-1)*chin(i)*etan(i);
 N(i,7) = 1/4*(1+chin(i))*(etan(i)-1)*chin(i)*etan(i);
 N(i,9) = 1/4*(1+chin(i))*(1+etan(i))*chin(i)*etan(i);
 N(i,3) = 1/4*(chin(i)-1)*(1+etan(i))*chin(i)*etan(i);
 N(i,4) = 1/2*(1-chin(i)^2)*(etan(i)-1)*etan(i);
 N(i,6) = 1/2*(1-chin(i)^2)*(etan(i)+1)*etan(i);
 N(i,2) = 1/2*(chin(i)-1)*(1-etan(i)^2)*chin(i);
 N(i,8) = 1/2*(chin(i)+1)*(1-etan(i)^2)*chin(i);
 N(i,5) = (1-chin(i)^2)*(1-etan(i)^2);
end

%straingverticiecen = zeros(nEle*5*3*2,1); %creo un vettore con i 3 strain sopra e sotto ai vertici del quad e al centro (5*3*2)
%for e=1:nEle %per ogni elemento
    %for j=1:6 %per ogni strain 
    %straingverticiecen((e-1)*30+5*(j-1)+1:(e-1)*30+5*j) = straing((e-1)*54+9*(j-1)+1:2:(e-1)*54+9*j);
    %end
%end
%Prendo ogni 2 posizioni per ogni 9 punti di Gauss, così prendo il
%i 4 nodi(1,3,7,9) e il centroide (5). 
%Questo vettore straingaussverticiecentroide è fatto 30 posizioni per un elemento, 3*2 strain per ogni
%nodo e per il centroide (3*2*5), ho primi 6 per primo elemento e via così
% Voglio togliere gli zeri dei centroidi
%straingvertici = zeros(nEle*6*4,1);
%straingvertici(1:2) = straingverticiecen(1:2);
%straingvertici(end-1:end) = straingverticiecen(end-1:end);
%for k=0:6*nEle-2
    %straingvertici(3+4*k:6+4*k) = straingverticiecen(4+5*k:7+5*k);
%end

%Ora voglio gli strain nei nodi
strainnodalixxtop = zeros(nEle*4,1);
strainnodaliyytop = zeros(nEle*4,1);
strainnodalixytop = zeros(nEle*4,1);
strainnodalixxbottom = zeros(nEle*4,1);
strainnodaliyybottom = zeros(nEle*4,1);
strainnodalixybottom = zeros(nEle*4,1);
for e=1:nEle
    %for j=1:6
 %strainnodali((e-1)*24+4*(j-1)+1:(e-1)*24+4*j) = N*straingvertici((e-1)*24+4*(j-1)+1:(e-1)*24+4*j); 
    strainnodalixxtop((e-1)*4+1:(e-1)*4+4) = N*straing((e-1)*9*6+1:(e-1)*9*6+9); 
    strainnodaliyytop((e-1)*4+1:(e-1)*4+4) = N*straing((e-1)*9*6+10:(e-1)*9*6+10+8);
    strainnodalixytop((e-1)*4+1:(e-1)*4+4) = N*straing((e-1)*9*6+19:(e-1)*9*6+19+8);
    strainnodalixxbottom((e-1)*4+1:(e-1)*4+4) = N*straing((e-1)*9*6+28:(e-1)*9*6+28+8);
    strainnodaliyybottom((e-1)*4+1:(e-1)*4+4) = N*straing((e-1)*9*6+37:(e-1)*9*6+37+8);
    strainnodalixybottom((e-1)*4+1:(e-1)*4+4) = N*straing((e-1)*9*6+46:(e-1)*9*6+46+8);
    %end
end
% Sfrutto il fatto che straing ha 54 strain per elemento, 9 xx top, 9 yy
% top e così via. Ogni pacchetto di 9 strain di Gauss ottengo 4 strain
% nodali. Avrò strainnodalixxtop per esempio che avrà 4 epsxxtop nei 4
% nodi del 1 elemento, poi 4 epxxtop nei 4 nodi del secondo elemento 
% e così via. 4nodi*60 elementi, avrò 6 vettori, uno per ogni eps.
strainnodali = zeros(nEle*4*6,1);
for e=1:nEle
    strainnodali((e-1)*24+1:(e-1)*24+4)=strainnodalixxtop((e-1)*4+1:(e-1)*4+4);
    strainnodali((e-1)*24+5:(e-1)*24+8)=strainnodaliyytop((e-1)*4+1:(e-1)*4+4);
    strainnodali((e-1)*24+9:(e-1)*24+12)=strainnodalixytop((e-1)*4+1:(e-1)*4+4);
    strainnodali((e-1)*24+13:(e-1)*24+16)=strainnodalixxbottom((e-1)*4+1:(e-1)*4+4);
    strainnodali((e-1)*24+17:(e-1)*24+20)=strainnodaliyybottom((e-1)*4+1:(e-1)*4+4);
    strainnodali((e-1)*24+21:(e-1)*24+24)=strainnodalixybottom((e-1)*4+1:(e-1)*4+4);
end
% Qui metto i 6 vettori in fila
% ! Qui sono fatti tali che ho 4 epsxx per i 4 nodi del primo elemento, poi 4
% epsyy e così via, 3 per ogni nodo sopra e sotto, 4*6=24 per ogni elemento.
% Quindi i primi 24 componenti sono xxtop, yytop, xytop, xxbottom ecc per
% il primo elemento, e via così per il secondo
strainnodalimat = zeros(nEle*4,6);
  for e=1:nEle
          for k=1:6
strainnodalimat((e-1)*4+1:(e-1)*4+4,k) = strainnodali((e-1)*24+(k-1)*4+1:(e-1)*24+(k-1)*4+4);
    end
  end
  %  Ho creato qui una matrice che ha sulle righe 4 nodi di ogni elemento
  %  (4*nEle) e sulle colonne le 3 eps top e le 3 eps bottom (3*2)
  % L'ordine dei 4 nodi di ogni pacchetto dovrebbe essere quello antiorario
  % della connectivity
l=1;
i = zeros(nNode,1);
  strainmediati = zeros(nNode,6);
for e=1:nEle
    for k=2:5
strainmediati(Cne(e,k),:)=strainmediati(Cne(e,k),:)+strainnodalimat(l,:); % se non è antiorario per Gauss è un casino 
i(Cne(e,k))=i(Cne(e,k))+1; %contatore aumenta ogni volta che sommo in un nodo
l = l+1;
    end
end
% Io ho considerato che procedendo per ogni elemento si vada in senso
% antiorario sia per la connectivity sia per i punti di Gauss, ma non so se
% sia giusto. Così facendo, io avanzo nella matrice con i nodi e so già che troverò i nodi della
% connectivity 
for l=1:nNode
strainmediati(l,:)=strainmediati(l,:)/i(l); %media per ogni nodo (ogni riga)
end



 
% Bm ha 60 matrici 27*24, poiché 3 epsilon per 9 punti di Gauss (3*9=27).
% Prime 9 righe sono punti di Gauss epsilonxx, seconde 9 epsyy, ultimi 9
% epsxy. 24 perché 4 nodi*6 DOFs. 60 sono gli elementi.
% Devo ottenere grandezze=FDF*DOFs. 
% Prima devo trasformare DOFs globali in locali, poi li posso moltiplicare
% per FDF derivate. 
% Occhio a z perché membranale + z*flessionale, c'è t che è spessore. 
% Report Disp sono spostamenti e Rot sono rotazioni. 
% Nella matrice (?) si hanno 6 DOFs del primo nodo, poi 6 per secondo e così
% via. 
% Mi ricavo DOFs nodali e posso confrontare con Patran, dato che ho il bdf.
% Sono due carichi concentrati negli estremi. Lui utilizza ijk come sistema
% di riferimento, x da nodo 1 a nodo 2. 
% Attenzione che con queste FDF trovo strain nei punti di Gauss, uno sarà
% il centroide, io faccio confronto con centroide. Centroide dovrebbe
% essere il quinto. In quella matrice il 5, 5+9, ecc sarà il centroide. Poi
% così posso confrontare con Patran.
% Perché ho due modelli identici? Perché ijk in uno è a 90° e nell'altro
% inclinato di 90°. x diventa y e y diventa -x. Si vede su mcoord. Devo
% ottenere quelle che erano y uguali a -x e resto non cambia. 
% Prossimo step è di calcolare deformazioni nei nodi a partire da Gauss. 
% Useremo interpolazione biquadratica, usiamo FDF come se fossero un
% elemento finito. Siccome codice è iFEM, queste FDF sono quelle che già
% calcoliamo per calcolare gli spostamenti. 
% Ora lavoro per elementi, calcolare DOFs nodali dell'elemento e calcolare
% gli strain, una cosa ottima sarebbe creare vettore DOFs q chge sarà (84*6)*1,
% mentre matrice (6*9*60)*(84*6), 6 poiché 3 top e 3 bottom = 6. 
