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 = ENu/((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;