Equation with Finite Element Spaces to calculate Internal Forces

General Discussion
1

I’ve been trying to calculate the internal forces each node gets after solving the mesh with the elasticity code provided at the tutorials. I have calculated the stiffness matrix and want to use the equation F= K*U to calculate the forces. To my understanding U and F are finite element spaces made by uu,vv and f1, f2. After solving everything and extracting f1, f2 from F it turns out they are filled with 0 and I can’t find what I’ve done wrong. My code follows below

real E1 = 21.5;
real sigma = 0.29;
real gravity = 0;

// Mesh
int bottombeam = 2;
border a(t=2, 0){x=0; y=t ;label=1;}	// left beam
border b(t=0, 10){x=t; y=0 ;label=bottombeam;}	// bottom of beam
border c(t=0, 2){x=10; y=t ;label=1;}	// rigth beam
border d(t=0, 10){x=10-t; y=2; label=3;}	// top beam
real xp = 5,yp=1;
real[int,int] P=[[xp,yp]]; // add point (xp,yp) in a mesh ..
mesh th = buildmesh(b(20) + c(5) + d(20) + a(5),points=P,fixeborder=1);
plot(th,wait=1);

fespace Ph(th,P0); // constant par triangle
Ph E=E1;

//fespace U(th, [P1, P1]);
//U [ux,uy];

fespace f(th, [P1, P1]);
f [f1,f2];

// Fespace
fespace Vh(th, [P1, P1]);
Vh [uu, vv], [w, s];


//  find midlle (xp,yp) dof 
w[]=0:w[].n-1; //  put dof number in array 
int dof1=w(xp,yp)+0.5;
int dof2=s(xp,yp)+0.5;
assert( abs(dof1-w(xp,yp)) + abs(dof2-s(xp,yp)) < 1e-6); 
cout << " dof "<< dof1 << " "<< dof2 << endl; 
real[int] brhs(Vh.ndof); // the add rhs ..
brhs = 0;
brhs[dof2]=  10;

// Macro
real sqrt2 = sqrt(2.);
macro epsilon(u1, u2) [dx(u1), dy(u2), (dy(u1) + dx(u2))/sqrt2] // EOM
macro div(u, v) (dx(u) + dy(v)) // EOM

// Problem (with solve)
func mu = E/(2*(1 + sigma));
func lambda = E*sigma/((1 + sigma)*(1 - 2*sigma));
//Not used cout << "Lambda = " << lamda << endl;
//Not used cout << "Mu = " << mu << endl;
cout << "Gravity = " << gravity << endl;
solve Elasticity ([uu, vv], [w, s], solver=sparsesolver)
	= int2d(th)(
		  lambda*div(w, s)*div(uu, vv)
		+ 2.*mu*(epsilon(w, s)' * epsilon(uu, vv))
	)
	- int2d(th)(gravity*s)
	+ brhs
	+ on(1, uu=0, vv=0)
	;
cout << "Max displacement = " << uu[].linfty << endl;


//Stiffness Matrix
varf bin(uu,vv) = int2d(th)(dx(uu) * dx(vv) + dy(uu) * dy(vv));
matrix K = bin(Vh, Vh);

real[int] U(K.n);
for(int i=0; i < K.n / 2; i++){
    U(2*i) = uu[](i);
    U(2*i+1) = vv[](i);
}

real[int] F(K.n);
F = (K * U);

for(int i=0; i < K.n / 2; i++){
    f1[](i) = F(2*i);
    f2[](i) = F(2*i+1);
}