Dear all,

I used the following code to extract the stiffness matrix of the grid as shown in the figure, with the left boundary fixed.

The structure is relatively simple. The stiffness matrix can be pushed manually. It should be an 8 * 8 matrix, but only 56 elements can be calculated by FreeFEM. See the figure below.

I guess it’s because some 0 elements are omitted, isn’t it?

Here is my code.

//Parameters

real Rho = 8000.; //Density

real E =210.; //Young modulus

real Nu = 0.27; //Poisson ratio

real Gravity = -9.81; //Gravity

mesh Th = square(1,1);

//Fespace

func Pk = P1;

fespace Uh(Th, [Pk, Pk]);

Uh [ux, uy];

//Macro

real sqrt2 = sqrt(2.);

macro Epsilon(ux, uy) [dx(ux), dy(uy), (dy(ux)+dx(uy))/sqrt2] //

macro Divergence(ux, uy) (dx(ux) + dy(uy)) //

//Problem

real Mu = E/(2.*(1.+Nu));
real Lambda = E*Nu/((1.+ Nu)*(1.-2.*Nu));

varf vElasticity ([ux,uy], [vx, vy])

= int2d(Th)(

Lambda * Divergence(vx, vy) * Divergence(ux, uy)

+ 2. * Mu * (

Epsilon(vx, vy)’ * Epsilon(ux, uy)

)

)

+ int2d(Th)(

Rho * Gravity * vy

)

+ on(4, uy=0,ux=0)

;

matrix Elasticity = vElasticity(Uh, Uh, solver=sparsesolver);

//Movemesh

Th = movemesh(Th, [x+ux, y+uy]);

[ux, uy] = [ux, uy];

//Plot

plot(Th,wait = 1);

cout <<" u = "<< Elasticity << endl;