%*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
%*~*										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.nodes = [0 0 ; 0.1 0; 0.2 0; 0.3 0; 0.4 0; 0.5 0; 0.6 0; 0.7 0; 0.8 0; 0.9 0; 1 0];
% MESH.elem = [1 2; 2 3; 3 4 ; 4 5; 5 6; 6 7; 7 8; 8 9; 9 10];
% MESH.dof = [3 3 3 3 3 3 3 3 3 3];
% MESH.mat = [1 1 1 1 1 1 1 1 1];

%definizione della mesh per trave una campata
nEF = 10;
xcoordMesh = linspace(0, 2, nEF + 1);
MESH.nodes = [xcoordMesh.', zeros(length(xcoordMesh), 1)];
MESH.elem = [];
for k = 1 : nEF
    MESH.elem = [MESH.elem; k, k + 1];
end
MESH.dof = 3 * ones(1, size(MESH.nodes, 1));
MESH.mat = 3 * ones(1, size(MESH.elem, 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)
axis equal
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.55;
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];
%A = FEM.defineInternalConstraints(MESH, matrixInternalConstraints);
A = [];

%vincoli esterni (vincoli a terra)
fixed = [1, 2, 3, 32];
valueFixed = [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 = 5;
loadGlobalConc = [0; -10; 0];
%fGlobal = FEM.addConcentratedLoad(MESH, fGlobal, nodeLoaded, loadGlobalConc);

%aggiungi carico distribuito
fFun = fLin;
elemLoaded = [1 : 1 : size(MESH.elem, 1)];
loadParameters = [0, 0, 2, 2];
fGlobal = FEM.addDistributedLoad(MESH, matrixGlobalDofs, fGlobal, fFun, Tfun, elemLoaded, loadParameters);
distributedLoadInfo.elemLoaded = [distributedLoadInfo.elemLoaded, elemLoaded];
distributedLoadInfo.loadParameters = [distributedLoadInfo.loadParameters, loadParameters];
distributedLoadInfo.elemLoaded = [distributedLoadInfo.elemLoaded, elemLoaded];
loadParametersInfo = [];
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 = [0.01, 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 = [0.01, 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
[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.003;
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 = 100;
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 = 3e9;
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, 2.55);
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];
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 = 5e0;
    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.25 2.25 -0.5 0.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