# Equation with Finite Element Spaces to calculate Internal Forces

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

``````