# Elasticity problem

**URL:** <https://community.freefem.org/t/elasticity-problem/2689>\
**Category:** General Discussion\
**Created:** [September 7, 2023, 9:34pm UTC](https://community.freefem.org/t/elasticity-problem/2689 "2023-09-07T21:34:30Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![moufide](https://avatars.discourse-cdn.com/v4/letter/m/48db29/32.png) [@moufide](https://community.freefem.org/u/moufide)\
**Post date:** [September 7, 2023, 9:34pm UTC](https://community.freefem.org/t/elasticity-problem/2689/1 "2023-09-07T21:34:30Z")

</div>

Bonjour cher Professeur

Ma compilation affiche je n’ai pas compris

– Square mesh : nb vertices =22 , nb triangles = 20 , nb boundary edges 22  
– FESpace: Nb of Nodes 44 Nb of DoF 132  
– Solve :  
min -1.3237e-12 max 0.345492  
current line = 68  
Exec error : Try to get unset x,y, …  
– number :1  
Exec error : Try to get unset x,y, …  
– number :1  
err code 8 , mpirank 0  
try getConsole C:\Users\ADMIN\Desktop\simulations\codee.edp  
save log in : ‘C:\Users\ADMIN\Desktop\simulations\codee.log’  
wait enter ?

voici mon code:

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],[u1oldd, u2oldd, u3oldd]=[sin(2_pi_x)_sin(2_pi_y)sin(2pi_z), sin(2_pi_x)_sin(2_pi_y)sin(2pi_z), 0]  
,[u1old, u2old, u3old]=[sin(2_pi_x)_sin(2_pi_y)sin(2pi_z), sin(2_pi_x)_sin(2_pi_y)sin(2pi_z), 0], [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)) //

int kk=0;

// Problem  
real mu = E/(2\*(1+sigma));  
real lambda = E_sigma/((1+sigma)_(1-2\*sigma));  
real dt=0.00000001, beta=0,T=0.000001, alpha=0;

real[int] instT(floor(T/dt) + 1);  
real[int] NORML2(floor(T/dt) + 1);

real[int] NORML(floor(T/dt) + 1);

for(real t=0; t \<T; t+=dt )  
{

solve Lame ([u1, u2, u3], [v1, v2, v3])  
= int3d(Th)((rho/dt^2)_[u1, u2, u3]'_[v1, v2, v3]  
+ 0.5_lambda_div(u1, u2, u3)_div(v1, v2, v3)  
+0.5_2._mu_( epsilon(u1, u2, u3)'\*epsilon(v1, v2, v3) )

```
              )
        +int2d(Th,1,2,3,4,6)((beta/dt)*[u1old, u2old, u3old]'*[v1, v2, v3]) 
    
        -int3d(Th)( (alpha/dt)*([u1old, u2old, u3old]'*[v1, v2, v3])
                    + 2*(rho/dt^2)*[u1old, u2old, u3old]'*[v1, v2, v3]
                    -(rho/dt^2)*[u1oldd, u2oldd, u3oldd]'*[v1, v2, v3]
                    - 0.5* lambda*div(u1oldd, u2oldd, u3oldd)*div(v1, v2, v3)
                    - 0.5*2.*mu*( epsilon(u1oldd, u2oldd, u3oldd)'*epsilon(v1, v2, v3) )
                     +(alpha/dt)*[u1old, u2old, u3old]'*[v1, v2, v3]
                 )

        -int2d(Th,1,2,3,4,6)( (beta/dt)*[u1oldd, u2oldd, u3oldd]'*[v1, v2, v3])
        + on(5, u1=0, u2=0, u3=0) ;

```

cout \<\< " " \<\< u1old\<\< u2old\<\< u3old\<\< u1\<\< u2 \<\<endl;

[up1,up2,up3]=[u1old-u1oldd, u2old-u2oldd, u3old-u3oldd];

instT[kk]=t;  
NORML2[kk]=int3d(Th)((rho/dt^2)_([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) ) )  
;

NORML[kk]= int2d(Th,1,2,3,4,6)((beta/dt)_[u1, u2, u3]'_[u1, u2, u3]) ;

[u1oldd, u2oldd, u3oldd]=[u1old, u2old, u3old];  
[u1old, u2old, u3old]=[u1, u2, u3];

cout \<\< " " \<\< NORML2[kk]\<\< " " \<\< instT[kk] \<\<endl;

kk++;

```
}

```

cout \<\< " " \<\< NORML2 \<\<endl;

J’aimerais envoyer le fichier.loq mais ça dit que je suis nouveau je ne peux pas envoyer un fichier.  
S’il ya lien pour faire l’envoie je serai ravi.

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [September 7, 2023, 10:06pm UTC](https://community.freefem.org/t/elasticity-problem/2689/2 "2023-09-07T22:06:06Z")

</div>

did you try making ny and nz bigger like the 10 I suggested?

---

<div class="post-metadata">

**Author:** ![moufide](https://avatars.discourse-cdn.com/v4/letter/m/48db29/32.png) [@moufide](https://community.freefem.org/u/moufide)\
**Post date:** [September 7, 2023, 10:19pm UTC](https://community.freefem.org/t/elasticity-problem/2689/3 "2023-09-07T22:19:44Z")

</div>

No i will try it now to see.

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [September 7, 2023, 10:43pm UTC](https://community.freefem.org/t/elasticity-problem/2689/4 "2023-09-07T22:43:09Z")

</div>

that alone won’t fix your problem see if this helps,

[code\_freefem.edp.edp](https://community.freefem.org/uploads/short-url/wEPI0HIrnxnLLsqh4AF4jZ8S8nW.edp) (2.3 KB)

---

<div class="post-metadata">

**Author:** ![moufide](https://avatars.discourse-cdn.com/v4/letter/m/48db29/32.png) [@moufide](https://community.freefem.org/u/moufide)\
**Post date:** [September 7, 2023, 10:44pm UTC](https://community.freefem.org/t/elasticity-problem/2689/5 "2023-09-07T22:44:29Z")

</div>

I increased but I still observe the same problem where the numerical results do not reflect the theoretical results.

---

<div class="post-metadata">

**Author:** ![moufide](https://avatars.discourse-cdn.com/v4/letter/m/48db29/32.png) [@moufide](https://community.freefem.org/u/moufide)\
**Post date:** [September 7, 2023, 11:22pm UTC](https://community.freefem.org/t/elasticity-problem/2689/6 "2023-09-07T23:22:54Z")

</div>

no the problem still persists

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [September 7, 2023, 11:54pm UTC](https://community.freefem.org/t/elasticity-problem/2689/7 "2023-09-07T23:54:46Z")

</div>

If you ran the last code I posted do you get a result that quickly  
decreases towards zero? Now its just a matter of getting the numbers  
right? Just looking at the results they could be exponentially decreasing  
which may be what you expect except for the time constant being wrong.  
It may help to post the theory too.

---

<div class="post-metadata">

**Author:** ![moufide](https://avatars.discourse-cdn.com/v4/letter/m/48db29/32.png) [@moufide](https://community.freefem.org/u/moufide)\
**Post date:** [September 8, 2023, 4:56am UTC](https://community.freefem.org/t/elasticity-problem/2689/8 "2023-09-08T04:56:13Z")

</div>

\Omega est un ouvert de \mathbb{R}^3  
\begin{equation}  
\begin{array}{ll}  
\rho \frac{\partial^2 u}{\partial t^2}- div (Ce(u))=0\\  
u=0 ~~\mbox{sur}~~ \Gamma\_D\\  
Ce(u).n= -\beta u ~~\mbox{sur}~~ \Gamma\_N  
\end{array}  
\end{equation}

\texbf{Formulation variationnelle}

\begin{equation}  
\displaystyle \int\_\Omega \rho \frac{\partial^2 u}{\partial t^2}\cdot v+ \displaystyle \int\_\Omega Ce(u): Ce(v)= -\beta\displaystyle \int\_{\Gamma\_N} \frac{\partial u}{\partial t} \cdot v  
\end{equation}

\textbf{\estimation d’Energie}

On pose  
Energie = \displaystyle \int\_\Omega \rho| \frac{\partial u}{\partial t}|^2 + \displaystyle \int\_\Omega Ce(u): Ce(u)  
\begin{equation}  
\displaystyle \frac{\partial }{\partial t} \left(  
\displaystyle \int\_\Omega \rho| \frac{\partial u}{\partial t}|^2 + \displaystyle \int\_\Omega Ce(u): Ce(u)  
\right)= -\beta\displaystyle \int\_{\Gamma\_N}| \frac{\partial u}{\partial t}|^2  
\end{equation}

Donc ici on a une décroissance simple et non une décroissance exponentielle.

Deplus lorsque \beta=0 la dérivée de l’énergie est nulle par conséquent l’énergie doit être une constante.  
C’est ce que dit la theory

---

<div class="post-metadata">

**Author:** ![moufide](https://avatars.discourse-cdn.com/v4/letter/m/48db29/32.png) [@moufide](https://community.freefem.org/u/moufide)\
**Post date:** [September 8, 2023, 11:18am UTC](https://community.freefem.org/t/elasticity-problem/2689/9 "2023-09-08T11:18:38Z")

</div>

Oui je suis censée obtenu que l’énergie est constant.

Par besoin d’une décroissance.

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [September 8, 2023, 1:31pm UTC](https://community.freefem.org/t/elasticity-problem/2689/10 "2023-09-08T13:31:50Z")

</div>

Can you describe the error or post the working code? I was curious now.  
Thanks.

---

<div class="post-metadata">

**Author:** ![moufide](https://avatars.discourse-cdn.com/v4/letter/m/48db29/32.png) [@moufide](https://community.freefem.org/u/moufide)\
**Post date:** [September 8, 2023, 1:46pm UTC](https://community.freefem.org/t/elasticity-problem/2689/11 "2023-09-08T13:46:52Z")

</div>

I haven’t found a good code yet… I myself cannot detect an error in the code that you published… I don’t know why but the result is not correct… when the parameters beta =0 the energy is theoretically constant.but this code prove that the energy are decrease.

---

<div class="post-metadata">

**Author:** ![marchywka](https://avatars.discourse-cdn.com/v4/letter/m/ee59a6/32.png) [@marchywka](https://community.freefem.org/u/marchywka)\
**Post date:** [September 8, 2023, 9:13pm UTC](https://community.freefem.org/t/elasticity-problem/2689/12 "2023-09-08T21:13:42Z")

</div>

Your kinetic term doesn’t seem to have a dt in it. If I reduce the time  
step to .0001 or so it still loses some energy. I guess you could check the  
time and length scales and see if the mesh and time step are fine enough  
or see how the energy loss varies with those parameters.  
Apparently your initial conditions are at zero derivative or maximum deflection  
with no instaneous motion.