# How to compute body forces
**URL:** https://community.freefem.org/t/how-to-compute-body-forces/1473
**Category:** General Discussion
**Created:** [February 5, 2022, 9:45am UTC](https://community.freefem.org/t/how-to-compute-body-forces/1473 "2022-02-05T09:45:34Z")
**Posts on this page:** 3
**Page:** 1
### Author:  [@distractor](https://community.freefem.org/u/distractor)
#### Post date: [February 5, 2022, 9:45am UTC](https://community.freefem.org/t/how-to-compute-body-forces/1473/1 "2022-02-05T09:45:34Z")
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.
```auto
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)](https://www.sciencedirect.com/book/9781856176330/the-finite-element-method-its-basis-and-fundamentals), 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:
```auto
// 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?
---
### Author:  [@vlaude](https://community.freefem.org/u/vlaude)
#### Post date: [February 6, 2022, 1:32pm UTC](https://community.freefem.org/t/how-to-compute-body-forces/1473/2 "2022-02-06T13:32:52Z")
Your first equation is Navier equation, Hooke’s law is the linear constitutive relation between stress and strain tensors. The body forces are prescribed, normally you don’t ‘compute’ them. Now I understand you want to estimate to what degree the equilibrium equation (2) is satisfied by the FEM solution. Note that you use a weak displacement formulation, so you can’t expect it to hold exactly.
That the average residual value you compute is of the order of 10^7 does not mean that there is necessarily a problem. You are using E=2.23 GPa=223 10^7 Pa to set the Lamé constants lambda and mu, so you should compare stresses to a similar scale.
Maybe you can observe the values of the different terms of the stress divergence to check if they actually are larger than 10^7 and cancel each other when you sum them up?
Now, you can also simplify your experiment by imposing displacement along only one axis, similar to a compression-expansion experiment:
```auto
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)
+ on(2,u1=0.025)
;
```
The result has very small average residual value. Clearly, this way you avoid the fast deformation around the clamped side resulting from `on(1,u1=0,u2=0,u3=0)`.
---
### Author:  [@distractor](https://community.freefem.org/u/distractor)
#### Post date: [February 6, 2022, 2:01pm UTC](https://community.freefem.org/t/how-to-compute-body-forces/1473/3 "2022-02-06T14:01:36Z")
Yes, normally you do not compute body forces. But this is only a small part of an iterative process, where computation of body forces is the main thing.
> [@vlaude](#):
>
> Maybe you can observe the values of the different terms of the stress divergence to check if they actually are larger than 10^7 and cancel each other when you sum them up?
I am not sure I understand what you mean by this. Could you be more specific?
> [@vlaude](#):
>
> Now, you can also simplify your experiment by imposing displacement along only one axis, similar to a compression-expansion experiment.
Yes, sure, but that’s a different problem. I may have better idea but not enough freefem knowledge. I had an idea to rewrite the body force term from the Navier equation using some calculus and divergence theorem, see [MFEM - Finite Element Discretization Library](https://mfem.org/fem_weak_form/) chapter Weak divergence. This mathematically seems correct, yet I am lacking some freefem knowledge to implement that.