%% Analisi di una struttura reticolare iperstatica
%
% DATI:
% F = 10 kN;         L = 1 m;                   E = 210 GPa;
% \Delta t = 20°C;        A = 150 mm^2;       \alpha = 11.6 * 10^(-6)
% \Delta = 2 mm;
%
%% Ordinamento delle variabili (solo promemoria)
% Le variabili in gioco saranno ordinate come segue:
% u = [uAx; uAy; uBx; uBy; uCx; uCy; uDx; uDy; uGx; uGy]
%
% uf = [uBy; uCx; uCy; uDx; uDy; uGy]
% Forze = [FBy; FCx; FCy; FDx; FDy; FGy] -->  Forze = [0; 0; -F; 0; 0; -F]
%
% uv = [uAx; uAy; uBx; uGx] --> DeltaEst = [-Delta; 0; 0; 0]
% R = [RAx; RAy; RBx; RGx]
%
% N  = [N1;  N2;  N3;  N4;  N5;  N6;  N7;  N8];
% deltaaste  = [d1;  d2;  d3;  d4;  d5;  d6;  d7;  d8]; <-- allungamenti nelle aste
% DAnelast  = [da1;  da2;  da3;  da4;  da5;  da6;  da7;  da8]; -->
%                  = [   0;      0;  +alfa Dt L3;  0;  -alfa Dt L5;  0;  0;  0];
%

% funzione per il calcolo del versore B-A
    vers = @(B,A) (B-A) / sqrt((B-A)'*(B-A));
% funzione per il calcolo del modulo B-A
    modulo = @(B,A) sqrt((B-A)'*(B-A));

%% 1. Definizione delle azioni/matrici di congruenza/matrice di rigidezza
%% Geometria

    syms L El Ar % El = modulo elastico, Ar = area
    assume([0 0 0]<=[L El Ar])

% coordinate dei punti (componenti per colonna)
    A = [0; 5*L];   B = [0; 0]; C = [3*L; 5*L];
    D = [3*L; 3*L]; G = [6*L; 5*L];

%% Azioni

syms F Delta % si indica genericamente con F l'intenzità della forza
       % il valore di riferiemnto del carico applicato
Forze  = [0; 0; -F; 0; 0; -F];

syms alpha Dt Delta % coeff. di dilatazione termica e Delta temperatura
DAnelast = alpha*Dt* [0; 0; modulo(D,A);  0;  -modulo(D,B);  0;  0;  0]; % allungamento anelastico delle aste (dist. termica)
DeltaEst  = [-Delta; 0; 0; 0];   % cedimenti vincolari esterni

%
%% Matrici di congruenza

% Matrice C, Cf e Cv
    vAB = vers(A,B);
    vCA = vers(C,A);
    vDA = vers(D,A);
    vCB = vers(C,B);
    vDB = vers(D,B);
    vDC = vers(D,C);
    vGC = vers(G,C);
    vGD = vers(G,D);
    zero = [0; 0];
    
    Ctot = [ vAB', -vAB', zero', zero', zero';...
            -vCA', zero',  vCA', zero', zero';...
            -vDA', zero', zero',  vDA', zero';...
            zero', -vCB',  vCB', zero', zero';...
            zero', -vDB', zero',  vDB', zero';...
            zero', zero', -vDC',  vDC', zero';...
            zero', zero', -vGC', zero',  vGC';...
            zero', zero', zero', -vGD',  vGD']
    IndLib  = [4, 5, 6, 7, 8, 10];
    IndVinc = [1, 2, 3, 9];
    
    Cf = Ctot(:,[IndLib])
    Cv = Ctot(:,[IndVinc])

%
%% Matrice di rigidezza assiale

    mAB = modulo(A,B);
    mCA = modulo(C,A);
    mDA = modulo(D,A);
    mCB = modulo(C,B);
    mDB = modulo(D,B);
    mDC = modulo(D,C);
    mGC = modulo(G,C);
    mGD = modulo(G,D);
    
    kel = El*Ar*diag([1/mAB, 1/mCA, 1/mDA, ...
                      1/mCB, 1/mDB, 1/mDC, ...
                      1/mGC, 1/mGD])
% clear mBA mDA mDB mCB mCA mCD
% clear vBA vDA vDB vCB vCA vCD vAx vAy vBx vCy vDy
% clear Ctot IndVinc IndLib

%
%% Calcolo attraverso le equazioni del metodo degli spostamenti
%% 2. Calcoli preliminari (metodo degli spostamenti)

Forze
Fhat = Forze - Cf'*kel*Cv*DeltaEst + Cf'*kel*DAnelast

%

Krig = Cf'*kel*Cf

%
%% 3. Calcolo degli spostamenti nodali (metodo degli spostamenti)

uf = inv(Krig)*Fhat;

%% 4. Calcolo degli allungamenti assiali (metodo degli spostamenti)

dassiali = Cf*uf + Cv*DeltaEst;

%% 5. Calcolo degli sforzi normali (metodo degli spostamenti)

SfNorm = kel*(dassiali-DAnelast);

%% 6. Calcolo delle reazioni vincolari (metodo degli spostamenti)

Reaz = Cv'*SfNorm;

%
%% Cacolo

% sostituzione per lunghezze generiche
DaSost = [El,     Ar,  L,    alpha,        Dt, F,     Delta];
Valori = [210000, 150, 1000, 11.6*10^(-6), 20, 10000, 2];

table({'N_1=';'N_2=';'N_3=';'N_4=';'N_5=';'N_6=';'N_7=';'N_8='},...
          double(subs(SfNorm,DaSost,Valori)),...
          {'N';'N';'N';'N';'N';'N';'N';'N'})
table({'H_A=';'V_A=';'H_B=';'H_G='},...
       double(subs(Reaz,DaSost,Valori)),...
       {'N';'N';'N';'N'})
table({'delta_1=';'delta_2=';'delta_3=';'delta_4=';'delta_5=';'delta_6='; 'delta_7=';'delta_8='},...
    double(subs(dassiali,DaSost,Valori)),...
    {'mm';'mm';'mm';'mm';'mm';'mm';'mm';'mm'})
table({ 'uBy'; 'uCx'; 'uCy'; 'uDx'; 'uDy'; 'uGy'},...
    double(subs(uf,DaSost,Valori)),...
    {'mm';'mm';'mm';'mm';'mm';'mm'})