kill(all);

assume(L>0);

coordinate nodali elemento, lato jk di lunghezza L a inclinazione alpha rispetto alla verticale

alpha:%pi/6;

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

fattore di scala visualizzazione azioni esterne

ampfac : 1;

definizione carichi distribuiti di superficie;

carico distribuito su solo lato jk con azione in direzione x variabile in y da p-dp a p+dp

sx_ij(x,y) := 0;
sy_ij(x,y) := 0;

sx_jk(x,y) := p+dp*( y - 1 - L*cos(alpha)/2 )/(L*cos(alpha)/2) ;
sy_jk(x,y) := 0;

sx_kl(x,y) := 0;
sy_kl(x,y) := 0;

sx_li(x,y) := 0;
sy_li(x,y) := 0;

visualizzazione elemento e carichi applicati

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]),x[3,1]+sx_jk(x[3,1],y[3,1]),x[3,1]]*ampfac,
           [y[2,1],y[2,1]+sy_jk(x[2,1],y[2,1]),y[3,1]+sy_jk(x[3,1],y[3,1]),y[3,1]]*ampfac
 ]
],L=4,p=1,dp=0.5
),
['x,0,10],
['y,0,10],
[legend,"elemento","sx,sy"]
);

funzioni di forma iso4

N1(xi,eta):=(1-xi)*(1-eta)/4;

N2(xi,eta):=(1+xi)*(1-eta)/4;

N3(xi,eta):=(1+xi)*(1+eta)/4;

N4(xi,eta):=(1-xi)*(1+eta)/4;

matrice N dei moti elementari di spostamento

Nmat(xi,eta) := matrix(
[N1(xi,eta),         0,N2(xi,eta),         0,N3(xi,eta),         0,N4(xi,eta),         0],
[         0,N1(xi,eta),         0,N2(xi,eta),         0,N3(xi,eta),         0,N4(xi,eta)]
);

vettore delle funzioni di forma nodali

define (
 Nvec(xi,eta) ,
 matrix(
  [N1(xi,eta),N2(xi,eta),N3(xi,eta),N4(xi,eta)]
 )
);

coordinate fisiche (x,y) in funzione delle coordinate naturali (xi,eta) 

x(xi,eta):=Nvec(xi,eta).x;
y(xi,eta):=Nvec(xi,eta).y;

derivate parziali della funzione di mappatura

define ( dx_dxi (xi,eta) , diff(x(xi,eta),xi ) );
define ( dx_deta(xi,eta) , diff(x(xi,eta),eta) );

define ( dy_dxi (xi,eta) , diff(y(xi,eta),xi ) );
define ( dy_deta(xi,eta) , diff(y(xi,eta),eta) );

integrazione sul perimetro: considero il solo contributo non nullo associato al lato jk

integrando, da integrarsi in eta da -1 a 1

int_jk : 

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 )  ;

integrale svolto in forma esatta

integrate(int_jk,eta,-1,1);

Fs_INT :fullratsimp(%);

integrale svolto per quadratura gaussiana a 2 punti

1*ev(int_jk,eta=1/sqrt(3))+1*ev(int_jk,eta=-1/sqrt(3))   ;

Fs_QUAD: fullratsimp(%);