Hook law in weak formulation is

where f are body forces (gravity for example). If there are no body forces, then the last integral simply equals zero.

So I made an example, where one side of the cube as displacement prescribed.

```
load "medit"
include "cube.idp"
int[int] numberOfNodes = [4, 4, 4]; // Nodes per side.
real [int, int] boundaryRange = [[0.0, 1.0], [0.0, 1.0], [0.0, 1.0]]; // Side lengths.
int [int, int] labels = [[1, 2], [3, 4], [5, 6]]; // Labels.
mesh3 Th = Cube(numberOfNodes, boundaryRange, labels);
real E = 2230e6;
real sigma = 0.4;
real mu = E/(2*(1+sigma));
real lambda = E*sigma/((1+sigma)*(1-2*sigma));
fespace Vh(Th,P2);
Vh u1,u2,u3, v1,v2,v3;
real sqrt2=sqrt(2.);
macro epsilon(u1,u2,u3) [dx(u1),dy(u2),dz(u3),(dz(u2)+dy(u3))/sqrt2,(dz(u1)+dx(u3))/sqrt2,(dy(u1)+dx(u2))/sqrt2] // EOM
macro div(u1,u2,u3) ( dx(u1)+dy(u2)+dz(u3) ) // EOM
solve Lame([u1,u2,u3],[v1,v2,v3])=
int3d(Th)(
lambda*div(u1,u2,u3)*div(v1,v2,v3)
+2.*mu*( epsilon(u1,u2,u3)'*epsilon(v1,v2,v3) ) //')
)
// - int3d(Th) (gravity*v3)
+ on(1,u1=0,u2=0,u3=0)
+ on(2,u1=0.025)
;
real dmax = u3[].max;
cout << " max deplacement = " << dmax << endl;
real coef= 0.1/dmax;
int[int] ref2=[1,0,2,0];
mesh3 Thm=movemesh3(Th,transfo=[x+u1*coef,y+u2*coef,z+u3*coef],label=ref2);
Thm=change(Thm,label=ref2);
plot(Th,Thm, wait=1,cmm="coef amplification = "+coef );
```

The results look fine:

# Problem

According to O.C. Zienkiewicz, … J.Z. Zhu, in The Finite Element Method: its Basis and Fundamentals (Seventh Edition), 2013

meaning the sum of stress tensor divergence div(sigma) and body forces b should be zero. Now in my case, there are NO body forces, because displacement is prescribed, meaning the divergence of the stress tensor should equal zero.

**But it does not**.

Here is how I computed the stress divergence:

```
// UPDATE
macro getTensorDivX(sxx, syy, szz, syz, sxz, sxy) (dx(sxx) + dy(sxy) + dz(sxz)) // EOM
macro getTensorDivY(sxx, syy, szz, syz, sxz, sxy) (dx(sxy) + dy(syy) + dz(syz)) // EOM
macro getTensorDivZ(sxx, syy, szz, syz, sxz, sxy) (dx(sxz) + dy(syz) + dz(szz)) // EOM
macro sigma(u1, u2, u3) [(lambda + 2 * mu) * dx(u1) + lambda * dy(u2) + lambda * dz(u3),
lambda * dx(u1) + (lambda + 2 * mu) * dy(u2) + lambda * dz(u3),
lambda + dx(u1) + lambda * dy(u2) + (lambda + 2 * mu) * dz(u3),
2 * mu * (dz(u2)+dy(u3)),
2 * mu * (dz(u1)+dx(u3)),
2 * mu * (dy(u1)+dx(u2))] // EOM
fespace Nh(Th, P1);
Nh sxx, syy, szz, syz, sxz, sxy;
Nh rx, ry, rz, r;
// Compute stresses.
sxx = sigma(u1, u2, u3)[0];
syy = sigma(u1, u2, u3)[1];
szz = sigma(u1, u2, u3)[2];
syz = sigma(u1, u2, u3)[3];
sxz = sigma(u1, u2, u3)[4];
sxy = sigma(u1, u2, u3)[5];
// Compute residuums.
rx = getTensorDivX(sxx, syy, szz, syz, sxz, sxy);
ry = getTensorDivY(sxx, syy, szz, syz, sxz, sxy);
rz = getTensorDivZ(sxx, syy, szz, syz, sxz, sxy);
for (int i = 0; i < r.n; i++) {
r[][i] = dist(rx[][i], ry[][i], rz[][i]);
}
cout << "Average body force (stress divergence) = " << r[].sum / r.n << endl;
```

The average residuum value is of order 10^7, which is not even close to zero. Does anybody know how to compute body forces in freefem?