/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 13.04.2 ] */

/* [wxMaxima: input   start ] */
kill(all);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
L dovendo andare sotto radice si assume positivo cosi' da evitare di lavorare con i moduli.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
assume( L>0 );
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Definizione delle coordinate nodali dell'elemento e dell'inclinazione alpha rispetto alla verticale.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
alpha:%pi/6;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Definizione dei vettori colonna contenenti le coordinate dei nodi i j k l;
in questo modo chiamando x(1) individuo univocamente la coordinata x del nodo i e analogamente per gli altri nodi.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
x : transpose(matrix(  [1 , 5 ,5 - L*sin(alpha),  1]  ));
y : transpose(matrix(  [1 , 1 ,1 + L*cos(alpha),  5]  ));
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Definizione di un fattore di amplificazione che rende visibili le azioni esterne sul grafico che plotteremo.
Senza questo comando il grafico che rappresenta il carico sull'elemento risulterebbe non fruibile.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
ampfac : 1;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Definizione dei carichi distribuiti sulle superfici;
ipotesi(come da Richiami teorici):
1_carico distribuito P sul solo lato jk 
2_retta d'azione di P in direzione x 
3_variazione di carico dP a farfalla lungo y.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
La formula che definisce il carico sul lato jk in direzione x 'sx_jk(x,y)' 
è composta da una costante che rappresenta il carico P e da una componente dP che vale 1 sul nodo k e -1 su j.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
sx_ij(x,y) := 0;
sy_ij(x,y) := 0;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
sx_jk(x,y) := p+dp*( y - 1 - L*cos(alpha)/2 )/(L*cos(alpha)/2) ;
sy_jk(x,y) := 0;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
sx_kl(x,y) := 0;
sy_kl(x,y) := 0;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
sx_li(x,y) := 0;
sy_li(x,y) := 0;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Visualizzazione dell'elemento e dei carichi applicati.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
-voglio disegnare delle linee spezzate, una funzione per l'elemento e una per il carico.
-ho definito numericamente i parametri L, alfa, p, dp per la graficazione.
-i due comandi 'discrete' servono a definire per punti il comportamento delle funzioni.

input di wxplot2d: 

-ev() serve per potere inserie i parametri numerici di cui sopra.

input 1: lista delle curve che voglio plottare, gli elementi devono essere separati da virgole e contenuti in parentesi quadre.

input 2: definisce il range entro il quale costruire il grafico;  i termini "'x" e "'y" con gli apici servono per distinguere 
le variabili x e y precedentemente definite come vettori colonna, dagli assi del grafico.

input 3: label associati alle curve plottate, scritti in ordine di inserimento.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
wxplot2d(
ev(
[
 [
    discrete,
    [x[1,1],x[2,1],x[3,1],x[4,1],x[1,1]],
    [y[1,1],y[2,1],y[3,1],y[4,1],y[1,1]]
 ] ,
 [
    discrete,
    [
        x[2,1],  
        x[2,1]+sx_jk(x[2,1],y[2,1])*ampfac,  
        x[3,1]+sx_jk(x[3,1],y[3,1])*ampfac,
        x[3,1]
    ],
    [
        y[2,1],
        y[2,1]+sy_jk(x[2,1],y[2,1])*ampfac,
        y[3,1]+sy_jk(x[3,1],y[3,1])*ampfac,
        y[3,1]
    ]
 ]
],L=4,p=1,dp=0.5
),
['x,0,10],
['y,0,10],
[legend,"elemento","sx,sy"]
);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Definizione delle funzioni di forma nelle coordinate naturali xi e eta.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
Ni(xi,eta):=(1-xi)*(1-eta)/4;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
Nj(xi,eta):=(1+xi)*(1-eta)/4;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
Nk(xi,eta):=(1+xi)*(1+eta)/4;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
Nl(xi,eta):=(1-xi)*(1+eta)/4;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Matrice N dei moti elementari di spostamento.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
Nmat(xi,eta) := matrix(
[Ni(xi,eta),         0,Nj(xi,eta),         0,Nk(xi,eta),         0,Nl(xi,eta),         0],
[         0,Ni(xi,eta),         0,Nj(xi,eta),         0,Nk(xi,eta),         0,Nl(xi,eta)]
);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
In maxima si possono definire le funzioni in 2 modi:

1_con ":" assegna al membro di sinistra ciò che sta a destra senza valutarlo

2_con "define" prima valuta l'espressione di destra e poi la assegna al membro di sinistra.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Definizione del vettore delle funzioni di forma nodali contenente tutte le funzioni di forma, utile per le derivazioni successive.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
define (
 Nvec(xi,eta) ,
 matrix(
  [Ni(xi,eta),Nj(xi,eta),Nk(xi,eta),Nl(xi,eta)]
 )
);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Per effettuare il cambio di variabile definisco le coordinate fisiche (x,y) in funzione delle coordinate naturali (xi,eta).

(il "." è il prodotto tra matrici)
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
x(xi,eta):=Nvec(xi,eta) . x;
y(xi,eta):=Nvec(xi,eta) . y;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Il codice seguente effettua le derivate parziali di x e y in xi e eta.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
define ( dx_dxi (xi,eta) , diff(x(xi,eta),xi ) );
define ( dx_deta(xi,eta) , diff(x(xi,eta),eta) );
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
define ( dy_dxi (xi,eta) , diff(y(xi,eta),xi ) );
define ( dy_deta(xi,eta) , diff(y(xi,eta),eta) );
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Integrazione:
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Integrazione sul perimetro: considero il solo contributo non nullo associato al lato jk.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
L'integrando e' da integrarsi in eta da -1 a 1.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Sul lato jk xi=1, eta varia da -1 a 1

t e' lo spessore dell'elemento.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
int_jk è l'integrando definito a partire dalla formula vista in teoria.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
int_jk : 

t*
matrix( 

  [  
        sx_jk(  x(1,eta)  ,   y(1,eta)  )  ,  
        sy_jk(  x(1,eta)  ,   y(1,eta)  )  
  ]  

)
. Nmat(1,eta) 
* sqrt ( dx_deta(1,eta)^2 + dy_deta(1,eta)^2 )  ;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
per semplificarlo un po'...
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
intjk:fullratsimp(%);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
"integrate" svolge l'integrale: il risultato sarà in forma esatta. Qualora non vi sia una primitiva il codice darà errore.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
integrate(int_jk,eta,-1,1);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
Fs_INT :fullratsimp(%);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Vettore composto dalle azioni nodali in direzione x e y sui nodi i j k l.
[x.i y.i x.j y.j x.k y.k x.l y.l]
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Per aggirare l'eventuale irrisolubilita' dell'integrale i codici fem usano l'approssimazione mediante la quadratura di Gauss perchè 
conduce sempre ad una soluzione vicina a quella esatta.
   [wxMaxima: comment end   ] */

/* [wxMaxima: comment start ]
Si riporta la soluzione alternativa con Gauss di secondo ordine(a 2 punti).
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
 1*  ev(int_jk,  eta= 1/sqrt(3))
+1*  ev(int_jk,  eta=-1/sqrt(3))   ;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
Fs_QUAD: fullratsimp(%);
/* [wxMaxima: input   end   ] */

/* Maxima can't load/batch files which end with a comment! */
"Created with wxMaxima"$