kill(all);

load("eigen");

andamento spostamenti su tria3

define(
        area(x_i,y_i,x_j,y_j,x_k,y_k),
        determinant(matrix([1,1,1],[x_i,x_j,x_k],[y_i,y_j,y_k]))/2
)$

define(
        N_i(x,y,x_i,y_i,x_j,y_j,x_k,y_k),
        area( x , y ,x_j,y_j,x_k,y_k)/area(x_i,y_i,x_j,y_j,x_k,y_k)
)$

define(
        N_j(x,y,x_i,y_i,x_j,y_j,x_k,y_k),
        area(x_i,y_i, x , y ,x_k,y_k)/area(x_i,y_i,x_j,y_j,x_k,y_k)
)$

define(
        N_k(x,y,x_i,y_i,x_j,y_j,x_k,y_k),
        area(x_i,y_i,x_j,y_j, x , y )/area(x_i,y_i,x_j,y_j,x_k,y_k)
)$

funzioni spostamento sul triangolo 3 nodi

u(x,y):=
        fullratsimp(
         u_i*N_i(x,y,x_i,y_i,x_j,y_j,x_k,y_k)
        +u_j*N_j(x,y,x_i,y_i,x_j,y_j,x_k,y_k)
        +u_k*N_k(x,y,x_i,y_i,x_j,y_j,x_k,y_k)
        );

v(x,y) :=
        fullratsimp(
         v_i*N_i(x,y,x_i,y_i,x_j,y_j,x_k,y_k)
        +v_j*N_j(x,y,x_i,y_i,x_j,y_j,x_k,y_k)
        +v_k*N_k(x,y,x_i,y_i,x_j,y_j,x_k,y_k)
        );

deformazioni (costanti..)

eps_x : diff(u(x,y),x);

eps_y : diff(v(x,y),y);

gamma_xy : diff(u(x,y),y) + diff(v(x,y),x);

eps:[eps_x,eps_y,gamma_xy];

delta:[u_i,v_i,u_j,v_j,u_k,v_k];

ricavo matrice B

B:coefmatrix(eps,delta);

matrice D, formulazione plane strain

D:matrix([1-nu,nu,0],[nu,1-nu,0],[0,0,(1-2*nu)/2])*E/(1+nu)/(1-2*nu);

costruisco matrice di rigidezza

define(
K(x_i,y_i,x_j,y_j,x_k,y_k,E,nu) , factor( 
             area(x_i,y_i,x_j,y_j,x_k,y_k)
            *transpose(B).D.B
)
);

my_x_i:0;
my_y_i:0;
my_x_j:1;
my_y_j:0;
my_x_k:0;
my_y_k:1;
my_E  :1000;
my_nu :3/10;

myK : K(
my_x_i,
my_y_i,
my_x_j,
my_y_j,
my_x_k,
my_y_k,
my_E  ,
my_nu
);

[vals,vecs] : eigenvectors(myK);

plotto

load(draw);

colorlist:[red         ,yellow ,       green    ,       blue    ,          cyan ,  magenta  ,  turquoise     ,   pink ,  coral  ];

plottami:[];

iscorri:0;
for i : 1 thru length(vals[2]) do(
for j : 1 thru vals[2][i] do ( /* scorro sulle molteplicità degli autovalori*/
tmp_norm:sqrt(lsum(x^2,x,vecs[i][j])),
disp(ev(vals[1][i],numer)),
tmp_u_i:vecs[i][j][1]/tmp_norm,
tmp_v_i:vecs[i][j][2]/tmp_norm,
tmp_u_j:vecs[i][j][3]/tmp_norm,
tmp_v_j:vecs[i][j][4]/tmp_norm,
tmp_u_k:vecs[i][j][5]/tmp_norm,
tmp_v_k:vecs[i][j][6]/tmp_norm,
iscorri:iscorri+1,
dx:(iscorri-1)*2,
plottami:append(plottami,
             [
             color=black,
             polygon([[my_x_i+dx,my_y_i],[my_x_j+dx,my_y_j],[my_x_k+dx,my_y_k]]),
             color=colorlist[iscorri],
             head_length = 0.1,
             vector([my_x_i+dx,my_y_i],[tmp_u_i, tmp_v_i ]),
             vector([my_x_j+dx,my_y_j],[tmp_u_j, tmp_v_j ]),
             vector([my_x_k+dx,my_y_k],[tmp_u_k, tmp_v_k ])]
)

))
$

wxdraw2d(      xrange=[-2,12],yrange=[-2,4],
             fill_color=white,
             plottami
);