clear all
close all

% ***********************   NB   ***********************
%Il programma costruisce la geometria per l'analisi con FEM tramite ANSYS
%del pannello sandwich date le dimensioni. Nella dimensione dello spessore
%del pannello si sviluppa una sola cella.
%Una prima parte è dedicata alla vera e propria costruzione della
%geometria, una seconda parte è dedicata a impostare mesh, vincoli e
%carichi in modo che non sia necessario ad ogni cambio di dimensioni andare
%a variarli a mano.
%Sono previste due condizioni di carico: carico lineare al mid-span e
%momento puro ai bordi.
% ***********************   NB   ***********************
% Il programma definisce i k-points partendo dall'origine del sistema
% di riferimento, prosegue su x (nel senso della lunghezza). Arrivato
% all'estremità del pannello si sposta su z (senso della larghezza) e
% ripete come sopra. L'ultima dimensione su cui si sposta è y (spessore)
% ******************************************************

%Indicare l'indirizzo della cartella in cui salvare il file txt della
%geometria
%ggpath='C:\Users\silvano\Documents\Paolo\Universidad\Masterarbeit zwei\ANSYS_models\Scripts\Comparison_models\';
ggpath='C:\Users\silvano\Documents\Paolo\Universidad\Masterarbeit zwei\ANSYS_models\Scripts\Sw.panels_scripts\';
%Definire numero di bordi vincolati (constrained edges). Se due soltanto
%(caso di sandwich beams) i bordi sono agli estremi lungo la x
%beams: c_ed=2; panels: c_ed=4
c_ed=4;
%Definire quale cella si usa col parametro cell=1;2;3;4;5
%dimensioni della cella: x,z base, y altezza
cell=4;
if cell==1
    bx=22;
    bz=bx;
    h=11;
elseif cell==2
    bx=22;
    bz=bx;
    h=15;
elseif cell==3
    bx=16.1;
    bz=bx;
    h=11;
elseif cell==4
    bx=30;
    bz=bx;
    h=15;
elseif cell==5
    bx=26.1;
    bz=bx;
    h=15;
end
%Definire lunghezza x e larghezza z della beam o panel. Lo spessore è pari
%all'altezza della cella
%Ricorda che L=176mm; poi si passa a 2L,3L,[...],8L...
%Multipli della cella altrimenti il programma non parte!
%NB: se cell==3 inserire lx e lz come "lx=(n_celle)*bx"
%lx=8*176;
lx=1500;
%lz=180;
lz=lx;
lg=300;
n=round(lx/bx);
m=round(lz/bz);
t=lg/bx;
n1=lx/bx;
m1=lz/bz;
if floor(n)==n && floor(m)==m && (c_ed==2 | c_ed==4)%
    %Di seguito matrici dei nomi dei k-points e delle loro coordinate.
    %Nella medesima posizione (i,j,s) delle matrici trovo nome, coordinate
    %x e z del singolo k-point. La matrice dei nomi ha 3 dimensioni per la
    %geometria di un volume. X e Z solo due perché i due pannelli sono
    %perfettamente sovrapposti. y è un vettore poiché i pannelli sono piani
    %e paralleli.
    K_NAME=zeros(2*m+1,2*n+1,2);
    X=zeros(2*m+1,2*n+1);%,2);
    y=[-h/2,h/2];
    Z=zeros(2*m+1,2*n+1);%,2);
    s=0;
    for s=1:2
        j=0;
        for j=1:2*n+1
            i=0;
            for i=1:2*m+1
                K_NAME(i,j,s)=5+5*(j-1)+5*(2*n+1)*(i-1)+5*(2*n+1)*(2*m+1)*(s-1);
                if s==1
                    X(i,j,s)=bx/2*(j-1);
                    Z(i,j,s)=bz/2*(i-1);
                end
            end
        end
    end
    %Definisco i k-points su cui si poggeranno i vincoli e che si trovano
    %alle estremità del pannello (senso della lunghezza) a y=0.
    %Definisco i NOMI dei k-points seguendo la logica muovo x poi z
    %C,c="constraint"
    C_K_NAME=zeros(2*m+1,2*n+1);
    C_K_NAME(1,1)=5+K_NAME(2*m+1,2*n+1,2);
    i=0;
    for i=1:2*m+1
        j=0;
        if i==1 && c_ed==2
            C_K_NAME(1,2*n+1)=5+C_K_NAME(1,1);
        elseif i==1 && c_ed==4
            for j=2:2*n+1
                C_K_NAME(i,j)=5+C_K_NAME(i,j-1);
            end
        elseif i==2*m+1 && c_ed==4
            C_K_NAME(i,1)=5+C_K_NAME(i-1,2*n+1);
            for j=2:2*n+1
                C_K_NAME(i,j)=5+C_K_NAME(i,j-1);
            end
        else
            C_K_NAME(i,1)=5+C_K_NAME(i-1,2*n+1);
            C_K_NAME(i,2*n+1)=5+C_K_NAME(i,1);
        end
    end
    %Con sprintf definisco il nome del file txt su cui vado a riportare la
    %geometria. Nel nome vengono specificate dimensioni e numero di celle
    %fileID=fopen(sprintf('Panel_l×b×h--%.1f×%.1f×%d_mm3--%d×%d×%d_cells',lx,lz,h,n,m,1),'w');
    fileID=fopen([ggpath sprintf('Panel__%d__l×b×h--%.1f×%.1f×%d_mm3--%d×%d×%d_cells.txt',cell,lx,lz,h,n,m,1)],'w');
    %Cicli per stampare i k-points nel file
    fprintf(fileID,'! K-points dei pannelli\r\n');
    s=0;
    for s=1:2
        i=0;
        for i=1:2*m+1
            j=0;
            for j=1:2*n+1
                fprintf(fileID,'k,%d,%.1f,%.1f,%.1f\r\n',K_NAME(i,j,s),X(i,j),y(s),Z(i,j));
            end
        end
    end
    %Cicli per stampare i k-points dei vincoli
    fprintf(fileID,'\r\n! K-points per i vincoli\r\n');
    i=0;
    for i=1:2*m+1
        j=0;
        for j=1:2*n+1
            if C_K_NAME(i,j)>0
                fprintf(fileID,'k,%d,%.1f,%d,%.1f\r\n',C_K_NAME(i,j),X(i,j),0,Z(i,j));
            end
        end
    end
    %Cicli per stampare le aree fra k-points dei pannelli
    fprintf(fileID,'\r\n! Aree dei pannelli\r\n');
    s=0;
    for s=1:2
        i=0;
        for i=1:2*m
            j=0;
            for j=1:2*n
                fprintf(fileID,'a,%d,%d,%d,%d\r\n',K_NAME(i,j,s),K_NAME(i,j+1,s),K_NAME(i+1,j+1,s),K_NAME(i+1,j,s));
            end
        end
    end
    %Cicli per stampare le aree di bordo del pannello
    %Le aree di bordo alle estremità secondo il senso x sono sempre
    %presenti; le altre compaiono solo con c_ed=4
    fprintf(fileID,'\r\n! Aree di bordo\r\n');
    j=0;
    for j=1:2*n:2*n+1
        i=0;
        for i=1:2*m
            s=0;
            for s=1:2
                fprintf(fileID,'a,%d,%d,%d,%d\r\n',C_K_NAME(i,j),K_NAME(i,j,s),K_NAME(i+1,j,s),C_K_NAME(i+1,j));
            end
        end
    end
    if c_ed==4
        fprintf(fileID,'\r\n');
        i=0;
        for i=1:2*m:2*m+1
            j=0;
            for j=1:2*n
                for s=1:2
                    fprintf(fileID,'a,%d,%d,%d,%d\r\n',C_K_NAME(i,j),K_NAME(i,j,s),K_NAME(i,j+1,s),C_K_NAME(i,j+1));
                end
            end
        end
    end
    %Cicli per stampare le aree del graticcio
    fprintf(fileID,'\r\n! Aree del grid\r\n');
    s=1;
    for i=1+2*t:2*t:2*m+1-2*t%in una direzione
        for j=1:2*n
            fprintf(fileID,'a,%d,%d,%d,%d\r\n',K_NAME(i,j,s),K_NAME(i,j+1,s),K_NAME(i,j+1,s+1),K_NAME(i,j,s+1));
        end
    end
    for j=1+2*t:2*t:2*n+1-2*t%nell'altra direzione
        for i=1:2*m
            fprintf(fileID,'a,%d,%d,%d,%d\r\n',K_NAME(i,j,s),K_NAME(i+1,j,s),K_NAME(i+1,j,s+1),K_NAME(i,j,s+1));
        end
    end
    %
    %2a PARTE
    %Mesh
    fprintf(fileID,'\r\n! MESHING\r\n');
    fprintf(fileID,'\r\n!/prep7\r\n');
    fprintf(fileID,'\r\n! Definisco con un parametro il numero di divisioni della mesh delle aree\r\nNDIV=30\r\n');
    %Mesh truss core
    fprintf(fileID,'\r\n!Mesh grid core\r\n');
    fprintf(fileID,'asel,s,loc,y,%d,%d\r\n',-floor((h-1)/2),floor((h-1)/2));
    fprintf(fileID,'asel,u,loc,x\r\nasel,u,loc,x,%.1f\r\n',lx);
    fprintf(fileID,'asel,u,loc,z\r\nasel,u,loc,z,%.1f\r\n',lz);
    fprintf(fileID,'aesize,all,NDIV\r\nmat,1\r\nsecnum,1\r\ntype,1\r\namesh,all\r\n');
    %Mesh pannelli
    fprintf(fileID,'\r\n!Mesh glass panes\r\n');
    fprintf(fileID,'asel,s,loc,y,%.1f\r\nasel,a,loc,y,%.1f\r\n',-h/2,h/2);
    fprintf(fileID,'asel,a,loc,x\r\nasel,a,loc,x,%.1f\r\n',lx);
    fprintf(fileID,'asel,a,loc,z\r\nasel,a,loc,z,%.1f\r\n',lz);
    fprintf(fileID,'aesize,all,NDIV\r\nmat,2\r\nsecnum,2\r\ntype,2\r\namesh,all\r\n');
    %Vincoli
    fprintf(fileID,'\r\n! VINCOLI\r\n');
    %Appoggio semplice
    fprintf(fileID,'\r\n! Vincoli sui due bordi laterali\r\n! Appoggio semplice\r\n');
    fprintf(fileID,'nsel,s,loc,y\r\nd,all,uy\r\nnsel,r,loc,x\r\nnplot\r\nd,all,ux\r\nd,all,uz\r\n');
    fprintf(fileID,'\r\nnsel,all\r\nnplot\r\n');
    %Carichi
    fprintf(fileID,'\r\n! CARICHI\r\n');
    %Wind load distribuito sul pannello inferiore
    fprintf(fileID,'\r\n! Wind load distribuito sul pannello inferiore, 1kN/m^2\r\n');
    fprintf(fileID,'asel,s,loc,y,%.1f\r\nsfa,all,,pres,-1/1000\r\n',-h/2);
    fprintf(fileID,'\r\nnsel,all\r\nnplot\r\n');
    %Carico al mid-span
    fprintf(fileID,'\r\n! Carico lineare al mid-span\r\n! F=100N\r\n!APPLOAD=100\r\n');
    fprintf(fileID,'!nsel,s,loc,x,%.2f\r\n!nsel,r,loc,y,%.1f\r\n',lx/2,-h/2);
    fprintf(fileID,'!*get,NNODES,node,0,count\r\n');
    fprintf(fileID,'!nsel,u,loc,z\r\n!nsel,u,loc,z,%.1f\r\n',lz);
    fprintf(fileID,'!f,all,fy,APPLOAD/(NNODES-1)\r\n');
    fprintf(fileID,'!nsel,s,loc,x,%.2f\r\n!nsel,r,loc,y,%.1f\r\n',lx/2,-h/2);
    fprintf(fileID,'!nsel,u,loc,z,1,%.1f\r\n',lz-1);
    fprintf(fileID,'!f,all,fy,APPLOAD/(NNODES-1)/2\r\n');
    fprintf(fileID,'\r\n!nsel,all\r\n!nplot\r\n');
    %Momento puro ai bordi
    fprintf(fileID,'\r\n! Momento applicato al bordo\r\n! F=1600N\r\n!APPLOAD=1600\r\n');
    fprintf(fileID,'!nsel,s,loc,x,%.1f\r\n!nsel,r,loc,y,%.1f\r\n',lx,-h/2);
    fprintf(fileID,'!*get,NNODES,node,0,count\r\n!f,all,fx,-APPLOAD/NNODES\r\n');
    fprintf(fileID,'!nsel,s,loc,x\r\n!nsel,r,loc,y,%.1f\r\n',-h/2);
    fprintf(fileID,'!*get,NNODES,node,0,count\r\n!f,all,fx,-APPLOAD/NNODES\r\n\r\n');
    fprintf(fileID,'!nsel,s,loc,x,%.1f\r\n!nsel,r,loc,y,%.1f\r\n',lx,h/2);
    fprintf(fileID,'!*get,NNODES,node,0,count\r\n!f,all,fx,-APPLOAD/NNODES\r\n');
    fprintf(fileID,'!nsel,s,loc,x\r\n!nsel,r,loc,y,%.1f\r\n',h/2);
    fprintf(fileID,'!*get,NNODES,node,0,count\r\n!f,all,fx,-APPLOAD/NNODES\r\n');
    fprintf(fileID,'\r\n!nsel,all\r\n!nplot\r\n');
    %Soluzione e postprocessing
    fprintf(fileID,'\r\n! SOLUTION AND POSTPROC\r\n');
    fprintf(fileID,'\r\n/solu\r\nsolve\r\n\r\n/post1\r\n');
    fprintf(fileID,'\r\n! Controllo reazione vincolare\r\n! Check value |R_y|=(wind.load)x(sup.)=%.3fN\r\n!prrfor,fy\r\n',1/1000*lx*lz);
    %Studio convergenza
    fprintf(fileID,'\r\n! STUDIO CONVERGENZA\r\n');
    fprintf(fileID,'\r\n/GRAPHICS,FULL\r\n');
    fprintf(fileID,'! Calcolo del PRERR el. shell (vetro)\r\n');
    fprintf(fileID,'!esel,s,mat,,2\r\n!prerr\r\n!esel,all\r\n');
    fprintf(fileID,'\r\n! Spostamento nodo centrale, linea mid-span, pannello inferiore\r\n');
    fprintf(fileID,'!nsel,s,loc,y,%.1f\r\n!nsel,r,loc,x,%.2f\r\n!nsel,r,loc,z,%.2f\r\n!prnsol,u,y\r\n',-h/2,lx/2,lz/2);
    fprintf(fileID,'! Considera elementi shell (tensione nel pannello)\r\n! Tensione eqv nodo su linea mid-span, pannello inferiore\r\n! Il nodo è scostato di mezza cella da quello centrale, non tocca il truss core\r\n');
    fprintf(fileID,'!nsel,s,loc,y,%.1f\r\n!nsel,r,loc,x,%.2f\r\n!nsel,r,loc,z,%.2f\r\n!prnsol,s,prin\r\n',-h/2,lx/2,lz/2-bz/2);
    fprintf(fileID,'\r\n!nsel,all\r\n');
    %Risultati
    fprintf(fileID,'\r\n! SPOSTAMENTI, DEFORMATA, TENSIONI\r\n');
    fprintf(fileID,'\r\n!pldisp,2\r\n');
    %fprintf(fileID,'\r\n! Spostamento dei punti mid-span di entrambi i pannelli\r\n');
    %fprintf(fileID,'! superiore\r\nnsel,s,loc,y,%.1f\r\nnsel,r,loc,x,%.2f\r\nprdisp\r\nnsel,all\r\n',h/2,lx/2);
    %fprintf(fileID,'! inferiore\r\nnsel,s,loc,y,%.1f\r\nnsel,r,loc,x,%.2f\r\nprdisp\r\nnsel,all\r\n',-h/2,lx/2);
    fprintf(fileID,'\r\n! Spostamento del punto centrale di entrambi i pannelli\r\n');
    fprintf(fileID,'! superiore\r\nnsel,s,loc,y,%.1f\r\nnsel,r,loc,x,%.2f\r\nnsel,r,loc,z,%.2f\r\nprdisp\r\nnsel,all\r\n',h/2,lx/2,lz/2);
    fprintf(fileID,'! inferiore\r\nnsel,s,loc,y,%.1f\r\nnsel,r,loc,x,%.2f\r\nnsel,r,loc,z,%.2f\r\nprdisp\r\nnsel,all\r\n',-h/2,lx/2,lz/2);
    fprintf(fileID,'\r\n!plnsol,s,eqv\r\n!plesol,s,eqv');
    %Fine
    fclose(fileID);
else
    fprintf('\n    ******** ATTENZIONE! ********\n Inserire dimensioni multiple dei valori di cella.\n Scegliere un numero di bordi vincolati fra c_ed=2 e c_ed=4.\n');
end

%Per visualizzare nella command window il file txt
%type (sprintf('Panel_l×b×h--%d×%d×%d_mm3--%d×%d×%d_cells',lx,lz,h,n,m,1))