# About how to use the lame formula with time domain for 3D plate

**URL:** <https://community.freefem.org/t/about-how-to-use-the-lame-formula-with-time-domain-for-3d-plate/905>\
**Category:** General Discussion\
**Created:** [April 12, 2021, 7:05am UTC](https://community.freefem.org/t/about-how-to-use-the-lame-formula-with-time-domain-for-3d-plate/905 "2021-04-12T07:05:02Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![HGH](https://avatars.discourse-cdn.com/v4/letter/h/3be4f8/32.png) [@HGH](https://community.freefem.org/u/HGH)\
**Post date:** [April 12, 2021, 7:05am UTC](https://community.freefem.org/t/about-how-to-use-the-lame-formula-with-time-domain-for-3d-plate/905/1 "2021-04-12T07:05:02Z")

</div>

Hello, guys

Regarding how to use the lame formula (elasiticity) to simulate a 1m x 0.5m x 0.002m plate, the condition of the plate is a simple support condition, and the internal grid of the plate is a tetrahedron, and external force is applied to the center of the plate, and then as time increases, record the midpoint of the plate displacement.

The following is the program I wrote, but now I am very troubled about how to write the simple support conditions of the plate and the relationship between vibration and time

// Parameters  
real E = 71.6e9, nu = 0.23;  
real f = -100;//gravity

// Mesh  
load “msh3”  
load “tetgen”  
load “medit”

real x0,x1,y0,y1,x2,x3,y2,y3;  
x0=0.; x1=1.; y0=0.; y1=0.5;  
int[int] lab1 = [20,21,22,23];  
mesh Thsq1 = square(100,100,[x0+(x1-x0)\*x,y0+(y1-y0)\*y]);

x2=0.; x3=1.; y2=0.; y3=0.5;  
int[int] lab2 = [24,25,26,27];  
mesh Th1 = square(100,100,[x2+(x3-x2)\*x,y2+(y3-y2)\*y]);

func ZZ1min = 0;  
func ZZ1max = 0.002;  
func XX1 = x;  
func YY1 = y;

int[int] ref31h = [0,11];  
int[int] ref31b = [0,12];

meshS Th31h = movemesh23(Thsq1,transfo=[XX1,YY1,ZZ1max],region=ref31h,orientation=1);  
savemesh(Th31h,“Th31h.mesh”);  
meshS Th31b = movemesh23(Th1,transfo=[XX1,YY1,ZZ1min],region=ref31b,orientation=-1);  
Th31h = change(Th31h, label=lab1);  
Th31b = change(Th31b, label=lab2);

/////////////////////////////////  
x0=0.; x1=1.; y0=0.; y1=0.002;  
int[int] lab3 = [26,28,22,29];  
mesh Thsq2 = square(100,1,[x0+(x1-x0)\*x,y0+(y1-y0)\*y]);

x2=0.; x3=1.; y2=0.; y3=0.002;  
int[int] lab4 = [24,31,20,30];  
mesh Thsq2q = square(100,1,[x2+(x3-x2)\*x,y2+(y3-y2)\*y]);

func ZZ2 = y;  
func XX2 = x;  
func YY2min = 0.;  
func YY2max = 0.5;

int[int] ref32h = [0,13];  
int[int] ref32b = [0,14];

meshS Th32h = movemesh23(Thsq2,transfo=[XX2,YY2max,ZZ2],region=ref32h,orientation=-1);  
savemesh(Th32h,“Th32h.mesh”);  
meshS Th32b = movemesh23(Thsq2q,transfo=[XX2,YY2min,ZZ2],region=ref32b,orientation=1);  
Th32h = change(Th32h, label=lab3);  
Th32b = change(Th32b, label=lab4);

/////////////////////////////////  
x0=0.; x1=0.5; y0=0.; y1=0.002;  
int[int] lab5 = [25,31,21,28];  
mesh Thsq3 = square(100,1,[x0+(x1-x0)\*x,y0+(y1-y0)\*y]);

x2=0.; x3=0.5; y2=0.; y3=0.002;  
int[int] lab6 = [27,30,23,29];  
mesh Thsq3q = square(100,1,[x2+(x3-x2)\*x,y2+(y3-y2)\*y]);

func XX3min = 0;  
func XX3max = 1.;  
func YY3 = x;  
func ZZ3 = y;

int[int] ref33h = [0,15];  
int[int] ref33b = [0,16];

meshS Th33h = movemesh23(Thsq3,transfo=[XX3max,YY3,ZZ3],region=ref33h,orientation=1);  
meshS Th33b = movemesh23(Thsq3q,transfo=[XX3min,YY3,ZZ3],region=ref33b,orientation=-1);

Th33h = change(Th33h, label=lab5);  
Th33b = change(Th33b, label=lab6);

////////////////////////////////  
meshS Th33 = Th31h+Th31b+Th32h+Th32b+Th33h+Th33b; // “gluing” surface meshs to obtain the surface of cube  
medit(“glumesh”,Th33);  
Th33=change(Th33, rmInternalEdges=true);  
savemesh(Th33,“Th33.mesh”);

// build a mesh of a axis parallel box with TetGen  
real[int] domain = [1,0.5,0.002,145,0.001];  
mesh3 Th44 = tetg(Th33,switch=“pqaAAYYQ”,nbofregions=1,regionlist=domain); // Tetrahelize the interior of the cube with tetgen

medit(“tetg”,Th44);  
savemesh(Th44,“Th44.mesh”);

// Problem  
real mu = E/(2\*(1 + nu));  
real lambda = E_nu/((1 + nu)_(1 - 2\*nu));

//Fespace  
fespace Vh(Th44,[P1,P1,P1]);  
Vh [u1,u2,u3], [v1,v2,v3] ;

varf m([u1,u2,u3],[v1,v2,v3]) = int3d(Th44) ();  
matrix mass=m(Vh,Vh);

cout \<\< “lambda,mu,f =”\<\<lambda\<\< " " \<\< mu \<\< " " \<\< f \<\< endl;

//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] // epsilon(u1,u2,u3)’\*epsilon(v1,v2,v3) = epsilon(u):epsilon(v)  
macro div(u1,u2,u3) ( dx(u1)+dy(u2)+dz(u3) ) // EOM

solve Lame([u1,u2,u3],[v1,v2,v3])=  
int3d(Th44)(lambda\*div(u1,u2,u3)\*div(v1,v2,v3)+2._mu_( epsilon(u1,u2,u3)'_epsilon(v1,v2,v3) ) )-int3d(Th44) (f_v3) + on(ref32h,u1=0,u2=0,u3=0) + on(ref32b,u1=0,u2=0,u3=0)+ on(ref33b,u1=0,u2=0,u3=0)+ on(ref33h,u1=0,u2=0,u3=0) ;

//Display  
real dmax= u1[].max;  
cout \<\< " max deplacement = " \<\< dmax \<\< endl;

// movemesh  
real coef= 0.1/dmax;  
int[int] ref2=[100,101,102,103];  
mesh3 Thm=movemesh3(Th44,transfo=[x+u1_coef,y+u2_coef,z+u3\*coef],label=ref2);  
Thm=change(Thm,label=ref2);

//Plot  
plot(Th44,Thm, wait=true,cmm="coef amplification = "+coef );