%*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
%*~*										PRELIMINARY DEFINITIONS										   *~*
%*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*

%delete previous results!
close all
clear all
clc

%working directory
%workDir = 'C:\Users\marot\Desktop\lezioniTrieste\Codici';

%add the directory of the libraries to the Matlab Path, in order to use the libraries of functions
%addLibrariesToPath(workDir);

global toleranceZero
toleranceZero = 1e-11;

%*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
%*~*										    MESH DEFINITION 										   *~*
%*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*

%print message
fprintf('-define the mesh...\n')

%definizione della mesh (coordinate nel sistema globale)

%Mesh 1) telaio con 3 EF
MESH.nodes = [0 0 ; 0 2; 2 2; 2 2; 2 0];
MESH.elem = [1 2; 2 3; 4 5];
MESH.dof = [3 3 3 3 3];
MESH.mat = 3*[1 1 1];

%Mesh 2) telaio con 6 EF
MESH.nodes = [0 0 ; 0 2; 2 2; 2 2; 2 0; 0.5 2; 1 2; 1.5 2 ];
MESH.elem = [1 2; 2 6; 6 7; 7 8; 8 3; 4 5];
MESH.dof = [3 3 3 3 3 3 3 3];
MESH.mat = 3*[1 1 1 1 1 1];

%calcolo delle lunghezze degli elementi
MESH = FEM.lengthElementMesh(MESH);

%calcolo dell'inclinazione degli elementi nel sistema di riferimento globale
MESH = FEM.angleElementMesh(MESH);

%matrice che crea un link tra Elementi Finiti e gradi di libertà
matrixGlobalDofs = FEM.linkLocalToGlobalDofs(MESH);

%grafico della struttura indeformata
figure(1)
title('Struttura Indeformata e mesh')
graphicOptions.color = FEM.libColor('black1');
graphicOptions.lineWidth = 2;
graphicOptions.lineStile = '-';
graphicOptions.marker = 'o';
graphicOptions.markerSize = 3;
FEM.plotMesh(MESH, graphicOptions);

%stampa sistema locale
graphicOptions.colorX = 'r';
graphicOptions.colorY = 'b';
graphicOptions.maxHeadSize = 0.3;
scaleFraction = 0.15;
FEM.plotRefSysLocal(MESH, scaleFraction, graphicOptions);




%*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
%*~*							FINITE ELEMENT FORMULATION DEFINITIONS									   *~*
%*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*

%definizione degli elementi da utilizzare
formulation = 'EULER_BERNOULLI';
Kfun = FEM.defineStiffnessMatrix(formulation);
Mfun = FEM.defineMassMatrix(formulation);
Tfun = FEM.defineRotationMatrix(formulation);
[fLin, fTermAx, fTermCurv] = FEM.defineElementLoad(formulation);




%*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
%*~*										CONSTRAINT DEFINITIONS										   *~*
%*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*

%print message
fprintf('-define the constraints of the problem...\n')

%vincoli interni (cerniere interne, pattini interni...)
matrixInternalConstraints = [7, 10; 8, 11; 9, 12]; %valido per entambe le mesh
A = FEM.defineInternalConstraints(MESH, matrixInternalConstraints);
%A = [];

%vincoli esterni (vincoli a terra)
fixed = [1, 2, 3, 13, 14, 15]; %valido per entambe le mesh
valueFixed = [0, 0, 0, 0, 0, 0];
uFixed = FEM.defineUFixed(MESH, fixed, valueFixed);
C = FEM.defineExternalConstraints(MESH, fixed);




%*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
%*~*										    LOAD DEFINITIONS										   *~*
%*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*

%print message
fprintf('-define the loads on the structure...\n')

%definizione iniziale dei carichi
fGlobal = sparse(zeros(sum(MESH.dof), 1));
distributedLoadInfo.elemLoaded = [];
distributedLoadInfo.loadParameters = {};
thermalAxLoadInfo.elemLoaded = [];
thermalAxLoadInfo.loadParameters = {};
thermalCurvLoadInfo.elemLoaded = [];
thermalCurvLoadInfo.loadParameters = {};

%aggiungi carico concentrato
nodeLoaded = 2;
loadGlobalConc = [1; 1; 2];
loadGlobalConc = [0; -10; 0];
loadGlobalConc = [0; 0; -10];
%fGlobal = FEM.addConcentratedLoad(MESH, fGlobal, nodeLoaded, loadGlobalConc);

nodeLoaded = 3;
loadGlobalConc = [1; 1; 2];
loadGlobalConc = [0; -10; 0];
loadGlobalConc = [0; 0; -10];
%fGlobal = FEM.addConcentratedLoad(MESH, fGlobal, nodeLoaded, loadGlobalConc);

%aggiungi carico distribuito
fFun = fLin;
elemLoaded = [2]; %mesh 1
elemLoaded = [2, 3, 4, 5]; %mesh 2
loadParameters = [0, 0, 2, 2];
fGlobal = FEM.addDistributedLoad(MESH, matrixGlobalDofs, fGlobal, fFun, Tfun, elemLoaded, loadParameters);
distributedLoadInfo.elemLoaded = [distributedLoadInfo.elemLoaded, elemLoaded];
for k = 1 : length(elemLoaded)
    distributedLoadInfo.loadParameters = [distributedLoadInfo.loadParameters, loadParameters];
end


%aggiungi carico termico (dilatazione assiale)
% fFun = fTermAx;
% elemLoaded = [1 : 1 : size(MESH.elem, 1)];
% loadParameters = [1, 0.1, 20];
% fGlobal = FEM.addThermalLoad(MESH, matrixGlobalDofs, fGlobal, fFun, Tfun, elemLoaded, loadParameters);
% thermalAxLoadInfo.elemLoaded = [thermalAxLoadInfo.elemLoaded, elemLoaded];
% for k = 1 : length(elemLoaded)
%     thermalAxLoadInfo.loadParameters = [thermalAxLoadInfo.loadParameters, loadParameters];
% end

%aggiungi carico termico (curvatura termica)
% fFun = fTermCurv;
% elemLoaded = [1 : 1 : size(MESH.elem, 1)];
% loadParameters = [1, 0.1, 20];
% fGlobal = FEM.addThermalLoad(MESH, matrixGlobalDofs, fGlobal, fFun, Tfun, elemLoaded, loadParameters);
% thermalCurvLoadInfo.elemLoaded = [thermalCurvLoadInfo.elemLoaded, elemLoaded];
% for k = 1 : length(elemLoaded)
%     thermalCurvLoadInfo.loadParameters = [thermalCurvLoadInfo.loadParameters, loadParameters];
% end


%%
%*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
%*~*										PROCESSING and SOLUTION                                        *~*
%*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*

%print message
fprintf('-compute the stiffness matrix and solve the linear system...\n')

%definizione della matrice di rigidezza
[KGlobal, MGlobal, KEFGlobal] = FEM.finiteElementAnalysis(MESH, matrixGlobalDofs, Kfun, Mfun, Tfun);



%soluzione del problema agli elementi finiti (ANALISI STATICA)
fprintf('- analisi statica...\n');
[U, K, F, R, Rconst] = FEM.solveStatic(KGlobal, fGlobal, fixed, valueFixed, A, C);



%soluzione del problema agli elementi finiti (ANALISI MODALE)
fprintf('- analisi modale...\n');
sigmaEig = 'sm';
numEig = 10;
[psi, lambda, omega, Umodal, naturalFrequency] = FEM.solveModal(MGlobal, KGlobal, valueFixed, A, C, sigmaEig, numEig);



%soluzione del problema agli elementi finiti (INTEGRAZIONE DIRETTA)
fprintf('- analisi dinamica per integrazione diretta...\n');

%define the time step
t0 = 0;
tEnd = 0.03;
numStep = 100;
singleDeltaT = (tEnd - t0) / numStep;
deltaT = singleDeltaT * ones(1, numStep);

%define the vector of time
vecT = [0];
for k = 1 : length(deltaT)
    vecT = [vecT sum(deltaT(1 : k))];
end
vecT = t0 + vecT;

%define the load function of time
syms t
omegaLoad = 300;
symTimeFun = cos(omegaLoad * t);
[Ft] = FEM.defineLoadTimePath(symTimeFun, fGlobal, deltaT, t0);


%define the parameter for the generalized-alpha method
rhoInfty = 1;
[alphaM, alphaF, theta1, theta2] = FEM.defineParametersGenAlpha(rhoInfty);

%define initial condition
U0 = sparse(sum(MESH.dof), 1);
dUt0 = sparse(sum(MESH.dof), 1);

%define the damping matrix
CGlobal = sparse(size(KGlobal, 1), size(KGlobal, 2));

%define the damping matrix (Rayleigh method)
% betaR = 0.5;
% CGlobal = ((1-betaR) * KGlobal) + (betaR * MGlobal);

%find dispacement, velocity and acceleration
[Ud, dUtd, dUttd] = FEM.timeIntegratorGeneralizedAlpha(KGlobal, CGlobal, MGlobal, Ft, valueFixed, A, C,...
												deltaT, U0, dUt0, alphaM, alphaF, theta1, theta2);
                                           
                                            
                                            
                                            
                                            
%% %*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
%*~*										POST-PROCESSING                                                *~*
%*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*

%print message
fprintf('-post processing and solution...\n')

%calcolo delle forze nodali per ciascun elemento
TfunPostProc = FEM.defineRotationMatrixPostProcessing(formulation);
[nodalForcesGlobal, nodalForcesLocal] = FEM.computeNodalForces(MESH, U, matrixGlobalDofs, KEFGlobal, TfunPostProc);

%calcolo delle funzioni che rappresentano i vari campi incogniti (soluzione elastica)
[solU, solV, solN, solT, solM] = FEM.defineSolutionField(formulation);

%calcolo delle funzioni che rappresentano i vari campi incogniti (soluzione dilatazione assiale)
[solNthermAx, solTtermthermAx, solMthermAx] = FEM.defineSolutionFieldThermAxial(formulation);

%calcolo delle funzioni che rappresentano i vari campi incogniti (soluzione dilatazione assiale)
[solNthermCurv, solTthermCurv, solMthermCurv] = FEM.defineSolutionFieldThermCurv(formulation);

%calcolo dei valori delle soluzioni dei campi spostamento ed azioni interne in dati punti su ciascun elemento finito
nSegment = 30; 
[coordX, valuesU, valuesV, valuesNel, valuesTel, valuesMel] = FEM.calculateValueSolution(MESH, U, matrixGlobalDofs, ... 
                                                        nSegment, nodalForcesLocal, distributedLoadInfo, ...
                                                             solU, solV, solN, solT, solM);

%calcolo dei valori delle soluzioni dei campi spostamento ed azioni interne in dati punti su ciascun elemento finito
[valuesNta, valuesTta, valuesMta] = FEM.calculateValueSolutionTherm(MESH, nSegment, ...
                                                    thermalAxLoadInfo, solNthermAx, solTtermthermAx, solMthermAx);

%calcolo dei valori delle soluzioni dei campi spostamento ed azioni interne in dati punti su ciascun elemento finito
[valuesNtc, valuesTtc, valuesMtc] = FEM.calculateValueSolutionTherm(MESH, nSegment, ...
                                                    thermalCurvLoadInfo, solNthermCurv, solTthermCurv, solMthermCurv);

%composizione dei cell array delle soluzioni
[valuesN, valuesT, valuesM] = FEM.sumCellArraySolution(valuesNel, valuesTel, valuesMel, valuesNta, valuesTta, valuesMta, ... 
                                                                                   valuesNtc, valuesTtc, valuesMtc);



%*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
%*~*										        PLOTS										           *~*
%*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*

%print message
fprintf('-create the plot...\n')

%grafico della deformata della struttura
figure(2)
axis equal
title('Deformata della Struttura')
graphicOptions.color = FEM.libColor('grey4');
graphicOptions.lineWidth = 2;
graphicOptions.lineStile = '-';
graphicOptions.marker = 'o';
graphicOptions.markerSize = 2;
FEM.plotMesh(MESH, graphicOptions);

graphicOptions.color = FEM.libColor('orange3');
graphicOptions.lineWidth = 3;
graphicOptions.lineStile = '-';
graphicOptions.marker = 'none';
graphicOptions.markerSize = 2;
scale = 1e9;
FEM.plotDeformed(MESH, coordX, valuesU, valuesV, scale, graphicOptions);



%grafico delle azioni interne (sforzo assiale)
figure(3)
axis equal
title('Grafico Azione Assiale')
graphicOptions.color = FEM.libColor('black1');
graphicOptions.lineWidth = 2;
graphicOptions.lineStile = '-';
graphicOptions.marker = 'o';
graphicOptions.markerSize = 2;
FEM.plotMesh(MESH, graphicOptions);

valueFun = valuesN;
graphicOptions.faceColor = FEM.libColor('green4');
graphicOptions.edgeColor = FEM.libColor('green4');
graphicOptions.faceAlpha = 0.85;
graphicOptions.edgeAlpha = 1;
scale = FEM.automaticScaling(MESH, valueFun, 0.15);
FEM.plotDiagram(MESH, coordX, valueFun, scale, graphicOptions);



%grafico delle azioni interne (Taglio)
figure(4)
axis equal
title('Grafico del Taglio')
graphicOptions.color = FEM.libColor('black1');
graphicOptions.lineWidth = 2;
graphicOptions.lineStile = '-';
graphicOptions.marker = 'o';
graphicOptions.markerSize = 2;
FEM.plotMesh(MESH, graphicOptions);

valueFun = valuesT;
graphicOptions.faceColor = FEM.libColor('blue3');
graphicOptions.edgeColor = FEM.libColor('blue3');
graphicOptions.faceAlpha = 0.85;
graphicOptions.edgeAlpha = 1;
scale = FEM.automaticScaling(MESH, valueFun, 0.15);
FEM.plotDiagram(MESH, coordX, valueFun, scale, graphicOptions);



%grafico delle azioni interne (Momento)
figure(5)
axis equal
title('Grafico del Momento')
graphicOptions.color = FEM.libColor('black1');
graphicOptions.lineWidth = 2;
graphicOptions.lineStile = '-';
graphicOptions.marker = 'o';
graphicOptions.markerSize = 2;
FEM.plotMesh(MESH, graphicOptions);

valueFun = valuesM;
graphicOptions.faceColor = FEM.libColor('red2');
graphicOptions.edgeColor = FEM.libColor('red2');
graphicOptions.faceAlpha = 0.85;
graphicOptions.edgeAlpha = 1;
scale = FEM.automaticScaling(MESH, valueFun, 0.15);
FEM.plotDiagram(MESH, coordX, valueFun, scale, graphicOptions);



%% *~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
%*~*										        PLOTS - MODAL								           *~*
%*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
%print message
fprintf('-post processing and solution (modal)...\n')

numSolModal = [1 2 3 4 5 6 7 8 9 10 11 12 13];
for j = 1 : length(numSolModal)
    %calcolo dei valori delle soluzioni dei campi spostamento ed azioni interne in dati punti su ciascun elemento finito
    nSegment = 30; 
    [coordXmod, valuesUmod, valuesVmod, valuesNel, valuesTel, valuesMel] = FEM.calculateValueSolution(MESH, Umodal(:, numSolModal(j)), matrixGlobalDofs, ... 
                                                            nSegment, nodalForcesLocal, distributedLoadInfo, ...
                                                                 solU, solV, solN, solT, solM);


    %grafico della deformata della struttura
    figure(20+j)
    axis equal
    title(strcat('Modo di Vibrare \omega =', sprintf('%f', naturalFrequency(numSolModal(j)))))
    graphicOptions.color = FEM.libColor('grey4');
    graphicOptions.lineWidth = 2;
    graphicOptions.lineStile = '-';
    graphicOptions.marker = 'o';
    graphicOptions.markerSize = 2;
    FEM.plotMesh(MESH, graphicOptions);

    graphicOptions.color = FEM.libColor('orange3');
    graphicOptions.lineWidth = 3;
    graphicOptions.lineStile = '-';
    graphicOptions.marker = 'none';
    graphicOptions.markerSize = 2;
    scale = 1e1;
    FEM.plotDeformed(MESH, coordXmod, valuesUmod, valuesVmod, scale, graphicOptions);
end



%% *~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
%*~*										        PLOTS - GENERALIZED ALPHA					           *~*
%*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
%print message
fprintf('-post processing and solution (Generalised Alpha Method)...\n')


numSolModal = [1 20 30 40 50];
for j = 1 : length(numSolModal)
    %calcolo dei valori delle soluzioni dei campi spostamento ed azioni interne in dati punti su ciascun elemento finito
    nSegment = 30; 
    [coordXmod, valuesUmod, valuesVmod, valuesNel, valuesTel, valuesMel] = FEM.calculateValueSolution(MESH, Ud(:, numSolModal(j)), matrixGlobalDofs, ... 
                                                            nSegment, nodalForcesLocal, distributedLoadInfo, ...
                                                                 solU, solV, solN, solT, solM);


    %grafico della deformata della struttura
    figure(100+j)
    axis equal
    title(strcat('Soluzione Dinamica Integrazione Diretta, t =', sprintf('%f', vecT(numSolModal(j)))))
    graphicOptions.color = FEM.libColor('grey4');
    graphicOptions.lineWidth = 2;
    graphicOptions.lineStile = '-';
    graphicOptions.marker = 'o';
    graphicOptions.markerSize = 2;
    FEM.plotMesh(MESH, graphicOptions);

    graphicOptions.color = FEM.libColor('orange3');
    graphicOptions.lineWidth = 3;
    graphicOptions.lineStile = '-';
    graphicOptions.marker = 'none';
    graphicOptions.markerSize = 2;
    scale = 1e9;
    FEM.plotDeformed(MESH, coordXmod, valuesUmod, valuesVmod, scale, graphicOptions);
end


figure(1000)

    
numSolModal = 1: length(Ud(1, :));
for j = 1 : length(numSolModal)
    %calcolo dei valori delle soluzioni dei campi spostamento ed azioni interne in dati punti su ciascun elemento finito
    nSegment = 30; 
    [coordXmod, valuesUmod, valuesVmod, valuesNel, valuesTel, valuesMel] = FEM.calculateValueSolution(MESH, Ud(:, numSolModal(j)), matrixGlobalDofs, ... 
                                                            nSegment, nodalForcesLocal, distributedLoadInfo, ...
                                                                 solU, solV, solN, solT, solM);
    

    %grafico della deformata della struttura
    clf
    axis equal
    axis([-0.5 2.5 -0.5 2.5])
    title(strcat('Soluzione Dinamica Integrazione Diretta, t =', sprintf('%f', vecT(numSolModal(j)))))
    graphicOptions.color = FEM.libColor('grey4');
    graphicOptions.lineWidth = 2;
    graphicOptions.lineStile = '-';
    graphicOptions.marker = 'o';
    graphicOptions.markerSize = 2;
    FEM.plotMesh(MESH, graphicOptions);
    
    graphicOptions.color = FEM.libColor('orange3');
    graphicOptions.lineWidth = 3;
    graphicOptions.lineStile = '-';
    graphicOptions.marker = 'none';
    graphicOptions.markerSize = 2;
    scale = 1e9;
    FEM.plotDeformed(MESH, coordXmod, valuesUmod, valuesVmod, scale, graphicOptions);
    
    hold off                                    
 	drawnow
 	pause(.1)
end