### Author:  [@moufide](https://community.freefem.org/u/moufide)
#### Post date: [September 7, 2023, 11:45am UTC](https://community.freefem.org/t/numerical-results-error/2688/1 "2023-09-07T11:45:43Z")
Hello everyone
I modified the code of the elasticity equation in 3D to evaluate the evolution of energy as a function of time…but I don’t know why it does not reflect what I obtained theoretically.
Normally when I put the beta=0 parameter, the energy must be constant and if beta is strictly greater than 0 the energy decreases.
load “msh3”
load “medit”
int nx =10, ny =1, nz = 1;
int[int] rup = [0,6]; // Haut du cube étiqueté par zéro
int[int] rdown = [0,5]; // Bas du cube étiqueté par zéro
int[int] rmid = [1,2,3,4]; // Parois latérales étiquetées
mesh3 Th = buildlayers(square(nx, ny, [0.5_x, 0.4_y]), nz, zbound=[0., 0.1],
labelmid=rmid, labelup = rup, labeldown = rdown);
real E = 21.5e4,rho=1600;
real sigma = 0.29;
real gravity = -0.05;
// Fespace
fespace Vh(Th, [P1, P1, P1]);
Vh [u1, u2, u3]=[10_sin(2_pi_x)sin(2pi_y)_sin(2_pi_z),10_ sin(2_pi_x)_sin(2_pi_y)sin(2pi_z), 10]
,[u1old, u2old, u3old]=[10_sin(2_pi_x)sin(2pi_y)_sin(2_pi_z),10_ sin(2_pi_x)_sin(2_pi_y)sin(2pi_z), 10],[u1oldd, u2oldd, u3oldd], [v1, v2, v3],[up1,up2,up3];
// Macro
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] //
macro div(u1, u2, u3) (dx(u1) + dy(u2) + dz(u3)) //
// Problem
int kk=0;
real mu = E/(2\*(1+sigma));
real lambda = E_sigma/((1+sigma)_(1-2\*sigma));
real dt=0.1, beta=0,T=1;
real[int] instT(floor(T/dt) + 1);
real[int] NORML2(floor(T/dt) + 1);
for(real t=0; t \<=T; t+=dt )
{
[u1oldd, u2oldd, u3oldd]=[u1old, u2old, u3old];
[u1old, u2old, u3old]=[u1, u2, u3];
solve Lame ([u1, u2, u3], [v1, v2, v3])
= int3d(Th)((rho/dt^2)_[u1, u2, u3]'_[v1, v2, v3]
+ lambda_div(u1, u2, u3)div(v1, v2, v3)
+ 2.mu( epsilon(u1, u2, u3)'epsilon(v1, v2, v3) )
)
-int3d(Th)(2(rho/dt^2)[u1old, u2old, u3old]'_[v1, v2, v3]
-(rho/dt^2)_[u1oldd, u2oldd, u3oldd]'_[v1, v2, v3])
```
+int2d(Th,1,2,3,4,6)((beta/dt)*[u1, u2, u3]'*[v1, v2, v3])
-int2d(Th,1,2,3,4,6)( (beta/dt)*[u1old, u2old, u3old]'*[v1, v2, v3])
+ on(5, u1=0, u2=0, u3=0) ;
```
[up1,up2,up3]=[u1-u1old, u2-u2old, u3-u3old];
NORML2[kk]=int3d(Th)(rho\*([up1,up2,up3]‘_[up1, up2, up3]))
+ int3d(Th)( lambda_div(u1, u2, u3)\*div(u1, u2, u3)
+ 2._mu_(epsilon(u1, u2, u3)’\*epsilon(u1, u2, u3) ) )
;
cout \<\< " " \<\< NORML2[kk] \<\
### Author:  [@marchywka](https://community.freefem.org/u/marchywka)
#### Post date: [September 7, 2023, 12:01pm UTC](https://community.freefem.org/t/numerical-results-error/2688/2 "2023-09-07T12:01:30Z")
can you upload the file? copy/paste always has errors.
I guess generally I’d try to breakdown the contributions to total energy
and see what they are doing. It may be something simple like a sign error.
Also you could have hidden sources or sinks
of energy at boundaries etc.
---
### Author:  [@moufide](https://community.freefem.org/u/moufide)
#### Post date: [September 7, 2023, 12:39pm UTC](https://community.freefem.org/t/numerical-results-error/2688/3 "2023-09-07T12:39:21Z")
Ok this is the attached the source code
(Attachment code\_freefem1 is missing)
---
### Author:  [@marchywka](https://community.freefem.org/u/marchywka)
#### Post date: [September 7, 2023, 12:49pm UTC](https://community.freefem.org/t/numerical-results-error/2688/4 "2023-09-07T12:49:58Z")
the post says attachment is miessing. Can you put it on github or something?
---
### Author:  [@moufide](https://community.freefem.org/u/moufide)
#### Post date: [September 7, 2023, 1:37pm UTC](https://community.freefem.org/t/numerical-results-error/2688/5 "2023-09-07T13:37:09Z")
ok i resend
[code\_freefem.edp.edp](https://community.freefem.org/uploads/short-url/9y6akJdIgTgGJnpQHFmqPsfV0KK.edp) (2.23 KB)
---
### Author:  [@moufide](https://community.freefem.org/u/moufide)
#### Post date: [September 7, 2023, 1:39pm UTC](https://community.freefem.org/t/numerical-results-error/2688/6 "2023-09-07T13:39:49Z")
this is the tex
(Attachment code.txt is missing)
---
### Author:  [@moufide](https://community.freefem.org/u/moufide)
#### Post date: [September 7, 2023, 1:42pm UTC](https://community.freefem.org/t/numerical-results-error/2688/7 "2023-09-07T13:42:49Z")
[code\_freefem.edp.edp](https://community.freefem.org/uploads/short-url/9y6akJdIgTgGJnpQHFmqPsfV0KK.edp) (2.23 KB)
---
### Author:  [@moufide](https://community.freefem.org/u/moufide)
#### Post date: [September 7, 2023, 4:32pm UTC](https://community.freefem.org/t/numerical-results-error/2688/8 "2023-09-07T16:32:07Z")
Est ce que vous avez reçu le fichier.
Je n’arrive pas a envoyé sur github.
---
### Author:  [@marchywka](https://community.freefem.org/u/marchywka)
#### Post date: [September 7, 2023, 6:36pm UTC](https://community.freefem.org/t/numerical-results-error/2688/9 "2023-09-07T18:36:58Z")
I didn’t get to the math as your value of concern seemed to drop
once I increased the numbers ny and nz and moved the variable kk
init and increment,
[code\_freefem.edp.edp](https://community.freefem.org/uploads/short-url/wEPI0HIrnxnLLsqh4AF4jZ8S8nW.edp) (2.3 KB)