Basically as the title says Im learning FreeFem for a university project and I want to better understand what the program is actually doing. I have the following program modified from https://modules.freefem.org/modules/elasticity/:
load "msh3"
load "gmsh"
//* mesh
int Clamped = 27; //Beam fixed
mesh3 Th = readmesh3("rails");
plot(Th, wait=true);
//* Some parameters defined here...
real E = 210e9; //* Young's Modulus
real Nu = 0.27; //* poisson's ratio
real Rho = 7800; //* density
real Gravity = 9.81; //* gravitational acceleration
//* fespace
func Pk = P1;
fespace Uh(Th, [Pk, Pk, Pk]);
Uh [ux, uy, uz];
Uh [vx, vy, vz];
Uh [uxp, uyp, uzp];
Uh [uxpp, uypp, uzpp];
//* Macro
real sqrt2 = sqrt(2.);
macro Epsilon(ux, uy, uz) [dx(ux), dy(uy), dz(uz),
(dz(uy)+dy(uz))/sqrt2,
(dz(ux)+dx(uz))/sqrt2,
(dy(ux)+dx(uy))/sqrt2] //
macro Divergence(ux, uy, uz) (dx(ux) + dy(uy) + dz(uz)) //
//* Problem
real Mu = E/(2.*(1.+Nu));
real Lambda = E*Nu/((1.+Nu)*(1.2.*Nu));
varf vElasticity ([ux, uy, uz], [vx, vy, vz])
= int3d(Th)(
Lambda * Divergence(vx, vy, vz) * Divergence(ux, uy, uz)
+ 2. * Mu * (
Epsilon(vx, vy, vz)' * Epsilon(ux, uy, uz)
)
)
+ int3d(Th)(
Rho * Gravity * vz
)
 int2d(Th, 30)(
30000*vy
)
+ on(Clamped, ux=0, uy=0, uz=0)
+ on(26, ux=0, uy=0, uz=0)
;
matrix Elasticity = vElasticity(Uh, Uh, solver=sparsesolver);
real[int] ElasticityBoundary = vElasticity(0, Uh);
ux[] = Elasticity^1 * ElasticityBoundary;
//*Movemesh
real dmax = ux[].max;
cout << "max displacement = " << dmax << endl;
real coef = 50;
//Th = movemesh(Th, [x+ux, y+uy, z+uz]);
//[ux, uy, uz] = [ux, uy, uz];
mesh3 Thm = movemesh3(Th, transfo=[x+ux*coef, y+uy*coef, z+uz*coef]);
Thm = change(Thm);
//*Plot
//plot([ux, uy, uz], value=true);
plot(Th, Thm, value=true, wait=true, cmm="coef amplification = "+coef);
A few of the points Im confused about are:

under fespace, what does
Uh [uxp, uyp, uzp]
andUh [uxpp, uypp, uzpp]
do/represent? These arenâ€™t used elsewhere so does not sure what it does. 
when solving the matrix vector equation Ax=b, why do the results get written to
ux[]
? I have applied the forces in x,y and z direction but it doesnâ€™t seem to matter? so why the x inux[]
? 
When applying external forces I am using
int2d(Th, No.)
. What doesint2d
represent? I also foundint3d
,int1d
andint0d
. Wouldint0d
for example be a point force? And the extra number I have to pass, does that correspond to a line, point or surface on my mesh? (Im preparing in gmsh fyi). 
The program I have can display the displacements and not the stresses. I am actually mostly interested in these stresses. What modifications would I need to calculate the stresses in the structure?
Sorry if my code is a bit bad, Im very new to this software and am trying to understand it better. Any help is greatly appreciated!
p.s. This is my result so far: