# Internal Forces at Elasticity Problem

**URL:** <https://community.freefem.org/t/internal-forces-at-elasticity-problem/1811>\
**Category:** General Discussion\
**Created:** [June 15, 2022, 2:52pm UTC](https://community.freefem.org/t/internal-forces-at-elasticity-problem/1811 "2022-06-15T14:52:53Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![Stefanosv](https://avatars.discourse-cdn.com/v4/letter/s/ee7513/32.png) [@Stefanosv](https://community.freefem.org/u/Stefanosv)\
**Post date:** [June 15, 2022, 2:52pm UTC](https://community.freefem.org/t/internal-forces-at-elasticity-problem/1811/1 "2022-06-15T14:52:53Z")

</div>

Hello,  
I have the following code written with the many help of Frédéric Hecht. Other than calculating the deformation, I need to calculate the internal forces that are applied to each node of the mesh.  
How would you suggest to calculate the internal forces?

> // Parameters  
> real E1 = 21.5;  
> real sigma = 0.29;  
> real gravity = -0.05;
> 
> // 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  
> 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;
> 
> // Plot  
> plot([uu, vv], wait=1);  
> //plot([uu, vv], wait=1, bb=[[-0.5, 2.5], [2.5, -0.5]]);
> 
> // Move mesh  
> mesh th1 = movemesh(th, [x+uu, y+vv]);  
> plot(th1, wait=1);
> 
> include “ffmatlib-Edited.idp”  
> savemesh(th,“export\_mesh.msh”);  
> ffSaveVh(th,Ph,“export\_vh.txt”);  
> ffSaveData(E, “efile.txt”)

---

<div class="post-metadata">

**Author:** ![Stefanosv](https://avatars.discourse-cdn.com/v4/letter/s/ee7513/32.png) [@Stefanosv](https://community.freefem.org/u/Stefanosv)\
**Post date:** [June 19, 2022, 11:17am UTC](https://community.freefem.org/t/internal-forces-at-elasticity-problem/1811/2 "2022-06-19T11:17:12Z")

</div>

I was thinking about calculating the stress of each node and multiplying it with the area of each node.