# Simulation of a beam

**URL:** <https://community.freefem.org/t/simulation-of-a-beam/3261>\
**Category:** General Discussion\
**Created:** [May 22, 2024, 5:31am UTC](https://community.freefem.org/t/simulation-of-a-beam/3261 "2024-05-22T05:31:46Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Oramastein1](https://avatars.discourse-cdn.com/v4/letter/o/e5b9ba/32.png) [@Oramastein1](https://community.freefem.org/u/Oramastein1)\
**Post date:** [May 22, 2024, 5:31am UTC](https://community.freefem.org/t/simulation-of-a-beam/3261/1 "2024-05-22T05:31:46Z")

</div>

Hello, good morning everyone, the truth is I’m new to this, I was supported to make a code for the simulation of a beam simply supported with a point load, but to run the code and display it does not contemplate the deformation of the beam, what I do not know if it is a coding issue because when changing the parameters does not change the shape of the beam, I do not know if you could support me. Thanks.

// Parámetros de la viga  
real L = 10; // Longitud de la viga  
real h = 1; // Altura de la viga  
real E = 70e9; // Módulo de Young  
real I = (h^4)/12; // Momento de inercia  
real P = -1000; // Carga puntual (negativa indica hacia abajo)  
real a = L/2; // Posición de la carga puntual desde el extremo izquierdo

// Malla  
border b1(t=0, L){x=t; y=0;};  
border b2(t=0, h){x=L; y=t;};  
border b3(t=L, 0){x=t; y=h;};  
border b4(t=h, 0){x=0; y=t;};  
mesh Th = buildmesh(b1(10) + b2(1) + b3(10) + b4(1));

// Espacio de funciones  
fespace Vh(Th, P2);  
Vh u, v;

// Definición de las condiciones de frontera  
varf vBoundary(u, v) = on(1,2, u=0);

// Definición de la carga puntual como una delta de Dirac en el punto ‘a’  
func delta = (abs(x-a)\<1e-2);

// Formulación variacional  
varf vPDE(u, v) =  
int2d(Th)(  
E_I_(dx(u)_dx(v))  
)  
+ on(1, 2, u=0)  
- int1d(Th, 1)(  
delta_P\*v  
);

// Matriz y vector  
matrix A = vPDE(Vh, Vh);  
real[int] b = vBoundary(0, Vh);

// Resolución del sistema  
set(A, solver=sparsesolver);  
real[int] sol = A^-1\*b;

// Proyección de la solución en el espacio de funciones  
u = sol;

// Visualización de la deformación de la viga  
plot(u, fill=true, value=true);

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [May 22, 2024, 9:12am UTC](https://community.freefem.org/t/simulation-of-a-beam/3261/2 "2024-05-22T09:12:53Z")

</div>

Hello,  
I think that there are the following issues in your code:

If you want the boundary conditions to be set on left and right, you must put `+on(2,4,u=0)` instead of `+on(1,2,u=0)`

You probably miss `+dy(u)*dy(v)` in addition to `dx(u)*dx(v)`

Your delta is too narrow, you should put `(abs(x-a)<L/10)`, knowing that you have 10 intervals on the boundary.

Moreover, to simplify you should use `solve vPDE(u, v)` instead of `varf vPDE(u, v)`, so that you can remove the lines `matrix A =...` until `u[] = sol;`

---

<div class="post-metadata">

**Author:** ![Oramastein1](https://avatars.discourse-cdn.com/v4/letter/o/e5b9ba/32.png) [@Oramastein1](https://community.freefem.org/u/Oramastein1)\
**Post date:** [May 22, 2024, 6:01pm UTC](https://community.freefem.org/t/simulation-of-a-beam/3261/3 "2024-05-22T18:01:48Z")

</div>

thank you very much, the tips helped me a lot, in fact I have almost the simulation, I just can’t get the top part to deflect, how can I do it?

 ![Captura de pantalla 2024-05-22 115431](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/e/ea74747c2a325db5cb9c563e8b0c6b0ecff43796.png)

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [May 23, 2024, 9:35am UTC](https://community.freefem.org/t/simulation-of-a-beam/3261/4 "2024-05-23T09:35:31Z")

</div>

I do not understand well the equations you want to solve. You should compare precisely the equations and the lines of code to see where the problem comes from.

In your code the Dirac mass is on boundary 1, which is the lower boundary. Thus it is at x=L/2, y=0. This may be the reason why the deformation is only on the bottom. If you want also a Dirac on the upper boundary you have to put int1d(Th,1,3). Maybe you would like a Dirac on the whole line x=L/2, y arbitrary, then you have to use int2d() instead of int1d().