# 3d elasticity problem

In my code when i use the type:solve i get results for the displacement but when i try to use varf in my problem the displecements are 0. Also i get a different result in the plotting of the movemesh although in the second case i get 0.

``````// Load necessary FreeFem++ libraries

// Parameters
int bottombeam = 2; //label of bottombeam
real E = 210;
real sigma = 0.29;
real gravity = -9.81;
real coef = 0.2;
real uMax= 10;
real D=3;
real rho = 1;

// Time parameters
real T = 1.0; // Total time
real dt = 0.01; // Time step
int nsteps = T/dt; // Number of time steps
real freq = 1.0; // Frequency of the wave
real amp = 1.0; // Amplitude of the wave

// Mesh
border a(t=0, 1) {x=t; y=0; label=1;};
border b(t=0, 2) {x=1; y=t; label=2;};
border c(t=1, 3) {x=t; y=2; label=3;};
border d(t=0, 2) {x=3; y=2-t; label=4;};
border e(t=0, 1) {x=3+t; y=0; label=5;};
border f(t=0, 1) {x=4; y=3*t; label=6;};
border g(t=1, 0) {x=4*t; y=3; label=7;};
border h(t=1, 0) {x=0; y=3*t; label=8;};

mesh th1 = buildmesh(a(10) + b(15) + c(20) + d(15) + e(10) + f(20) + g(30) + h(20));
plot(th1, wait=true, dim=2, fill=true);

// Construction 3d mesh
func zmin = 0;
func zmax = 2.;
int MaxLayer = 10;

mesh3 th2 = buildlayers(th1, MaxLayer, zbound=[zmin,zmax]);
medit("th1", th2);

// Fespace (Solid)
fespace Uh(th2,[P1,P1,P1]);
Uh [uux,uuy,uuz];
Uh [vx,vy,vz];
Uh [uuxp, uuyp, uuzp];
Uh [uuxpp, uuypp, uuzpp];

// Macro
real sqrt2=sqrt(2.);
macro epsilon(uux,uuy,uuz)  [dx(uux),dy(uuy),dz(uuz),(dz(uuy)+dy(uuz))/sqrt2,(dz(uux)+dx(uuz))/sqrt2,(dy(uuz)+dx(uuy))/sqrt2]   // EOM
macro Divergence(uux,uuy,uuz) ( dx(uux)+dy(uuy)+dz(uuz) )   // EOM

real mu = E/(2*(1+sigma));
real lambda = E*sigma/((1+sigma)*(1-2*sigma));

cout <<"Lambda = "<< gravity << endl;
cout << "Mu = " << mu << endl;

solve vElasticity([uux,uuy,uuz],[vx,vy,vz])=
int3d(th2)(
lambda*Divergence(vx,vy,vz)*Divergence(uux,uuy,uuz)
+2.*mu*(epsilon(uux,uuy,uuz)'*epsilon(vx,vy,vz) )
)
- int3d(th2) (rho*gravity*vy)
+ on(1,2,3,4,5, uux=0, uuy=0 , uuz=0)
;

cout << "Max displacement = " << uux[].linfty << endl;
cout << "Max displacement = " << uuy[].linfty << endl;
cout << "Max displacement = " << uuz[].linfty << endl;

//2nd case
varf vElasticity([uux,uuy,uuz],[vx,vy,vz])=
int3d(th2)(
lambda*Divergence(vx,vy,vz)*Divergence(uux,uuy,uuz)
+2.*mu*(epsilon(uux,uuy,uuz)'*epsilon(vx,vy,vz) )
)
- int3d(th2) (rho*gravity*vy)
+ on(1,2,3,4,5, uux=0, uuy=0 , uuz=0)
;

cout << "Max displacement = " << uux[].linfty << endl;
cout << "Max displacement = " << uuy[].linfty << endl;
cout << "Max displacement = " << uuz[].linfty << endl;

matrix Elasticity = vElasticity(Uh, Uh, solver=sparsesolver);
real[int] ElasticityBoundary = vElasticity(0, Uh);
uux[] = Elasticity^-1 * ElasticityBoundary;

th2 = movemesh(th2,[x+uux, y+uuy, z+uuz]);
[uux, uuy, uuz] = [uux, uuy, uuz];
mesh3 th4 = movemesh3(th2, transfo=[x+uux, y+uuy, z+uuz]);

plot([uux, uuy, uuz], value = true, cmm="u");
plot(th2, wait=true, fill=true, value=true, dim=2);
medit("displacement", th4);
plot(th4, value = true, wait = 1);
medit("utotal",th2,[uux,uuy,uuz],order=1,wait=1);
``````

My correction trivial,

zoumpou.edp (3.1 KB)

the output:

``````
-- FreeFem++ v4.14 (Ven  7 jui 2024 11:35:37 CEST - git v4.14-60-g660dbd36)
file : zoumpou.edp
1 : //// Load necessary FreeFem++ libraries
5 :
6 : //// Parameters
7 : int bottombeam = 2; ////label of bottombeam
8 : real E = 210;
9 : real sigma = 0.29;
10 : real gravity = -9.81;
11 : real coef = 0.2;
12 : real uMax= 10;
13 : real D=3;
14 : real rho = 1;
15 :
16 : //// Time parameters
17 : real T = 1.0; //// Total time
18 : real dt = 0.01; //// Time step
19 : int nsteps = T/dt; //// Number of time steps
20 : real freq = 1.0; //// Frequency of the wave
21 : real amp = 1.0; //// Amplitude of the wave
22 :
23 : //// Mesh
24 : border a(t=0, 1) {x=t; y=0; label=1;};
25 : border b(t=0, 2) {x=1; y=t; label=2;};
26 : border c(t=1, 3) {x=t; y=2; label=3;};
27 : border d(t=0, 2) {x=3; y=2-t; label=4;};
28 : border e(t=0, 1) {x=3+t; y=0; label=5;};
29 : border f(t=0, 1) {x=4; y=3*t; label=6;};
30 : border g(t=1, 0) {x=4*t; y=3; label=7;};
31 : border h(t=1, 0) {x=0; y=3*t; label=8;};
32 :
33 : mesh th1 = buildmesh(a(10) + b(15) + c(20) + d(15) + e(10) + f(20) + g(30) + h(20));
34 : plot(th1, wait=true, dim=2, fill=true);
35 :
36 : //// Construction 3d mesh
37 : func zmin = 0;
38 : func zmax = 2.;
39 : int MaxLayer = 10;
40 :
41 : mesh3 th2 = buildlayers(th1, MaxLayer, zbound=[zmin,zmax]);
42 : medit("th1", th2);
43 :
44 :
45 : //// Fespace (Solid)
46 : fespace Uh(th2,[P1,P1,P1]);
47 : fespace Vh(th2,[P1]);
48 : Uh [uux,uuy,uuz];
49 : Uh [vx,vy,vz];
50 : Uh [uuxp, uuyp, uuzp];
51 : Uh [uuxpp, uuypp, uuzpp];
52 :
53 : //// Macro
54 : real sqrt2=sqrt(2.);
55 : macro epsilon(uux,uuy,uuz)  [dx(uux),dy(uuy),dz(uuz),(dz(uuy)+dy(uuz))/sqrt2,(dz(uux)+dx(uuz))/sqrt2,(dy(uuz)+dx(uuy))/sqrt2]    )  //// EOM
56 : macro Divergence(uux,uuy,uuz) ( dx(uux)+dy(uuy)+dz(uuz) )    )  //// EOM
57 :
58 : real mu = E/(2*(1+sigma));
59 : real lambda = E*sigma/((1+sigma)*(1-2*sigma));
60 :
61 : cout <<"Lambda = "<< gravity << endl;
62 : cout << "Mu = " << mu << endl;
63 :
64 : solve vElasticity([uux,uuy,uuz],[vx,vy,vz])=
65 :   int3d(th2)(
66 :             lambda*Divergence(vx,vy,vz)       ( dx(vx)+dy(vy)+dz(vz) )   *Divergence(uux,uuy,uuz)       ( dx(uux)+dy(uuy)+dz(uuz) )
67 :             +2.*mu*(epsilon(uux,uuy,uuz)       [dx(uux),dy(uuy),dz(uuz),(dz(uuy)+dy(uuz))/sqrt2,(dz(uux)+dx(uuz))/sqrt2,(dy(uuz)+dx(uuy))/sqrt2]   '*epsilon(vx,vy,vz)       [dx(vx),dy(vy),dz(vz),(dz(vy)+dy(vz))/sqrt2,(dz(vx)+dx(vz))/sqrt2,(dy(vz)+dx(vy))/sqrt2]    )
68 :               )
69 :   - int3d(th2) (rho*gravity*vy)
70 :   + on(1,2,3,4,5, uux=0, uuy=0 , uuz=0)
71 :   ;
72 : Vh u1=uux,u2=uuy,u3=uuz;
73 : //// to see signe problem:
74 : cout << "Min/Max displacement x = " << u1[].min << " " << u1[].max << endl;
75 : cout << "Min/Max displacement y = " << u2[].min << " " << u2[].max<< endl;
76 : cout << "Min/Max displacement z = " << u3[].min << " " << u3[].max<< endl;
77 : ////2nd case
78 : varf vElasticity1([uux,uuy,uuz],[vx,vy,vz])=
79 :   int3d(th2)(
80 :             lambda*Divergence(vx,vy,vz)       ( dx(vx)+dy(vy)+dz(vz) )   *Divergence(uux,uuy,uuz)       ( dx(uux)+dy(uuy)+dz(uuz) )
81 :             +2.*mu*(epsilon(uux,uuy,uuz)       [dx(uux),dy(uuy),dz(uuz),(dz(uuy)+dy(uuz))/sqrt2,(dz(uux)+dx(uuz))/sqrt2,(dy(uuz)+dx(uuy))/sqrt2]   '*epsilon(vx,vy,vz)       [dx(vx),dy(vy),dz(vz),(dz(vy)+dy(vz))/sqrt2,(dz(vx)+dx(vz))/sqrt2,(dy(vz)+dx(vy))/sqrt2]    )
82 :               )
83 :   + int3d(th2) (rho*gravity*vy)
84 :   + on(1,2,3,4,5, uux=0, uuy=0 , uuz=0)
85 :   ;
86 :
87 :  u1=uux;u2=uuy;u3=uuz;
88 : cout << "Min/Max displacement x = " << u1[].min << " " << u1[].max << endl;
89 : cout << "Min/Max displacement y = " << u2[].min << " " << u2[].max<< endl;
90 : cout << "Min/Max displacement z = " << u3[].min << " " << u3[].max<< endl;
91 : matrix Elasticity = vElasticity1(Uh, Uh, solver=sparsesolver);
92 : real[int] ElasticityBoundary = vElasticity1(0, Uh);
93 : uux[] = Elasticity^-1 * ElasticityBoundary;
94 :
95 :
96 : th2 = movemesh(th2,[x+uux, y+uuy, z+uuz]);
97 : [uux, uuy, uuz] = [uux, uuy, uuz];
98 : mesh3 th4 = movemesh3(th2, transfo=[x+uux, y+uuy, z+uuz]);
99 :
100 : plot([uux, uuy, uuz], value = true, cmm="u");
101 : plot(th2, wait=true, fill=true, value=true, dim=2);
102 : medit("displacement", th4);
103 : plot(th4, value = true, wait = 1);
104 : medit("utotal",th2,[uux,uuy,uuz],order=1,wait=1); sizestack + 1024 =7008  ( 5984 )

--  mesh:  Nb of Triangles =   1072, Nb of Vertices 607
version de medit ffmedit -popen -filebin 1  th1
-- Medit,  Release 3.0a (Nov. 30, 2007)
compiled:  Ven 7 jui 2024 11:36:37 CEST (with ff++ 4.14).

medit with binary version of popen : mesh(es)
mesh_name= th1
End of mesh
Input seconds:     0.01

medit1()

Building scene(s)
Creating scene 1
Scene seconds:     0.02

Rendering scene(s)
FALLBACK (log once): Fallback to SW fragment processing (outputInfo.polygonModeMismatch)
FALLBACK (log once): Fallback to SW fragment processing, m_disable_code: 800000
FALLBACK (log once): Fallback to SW fragment processing in drawCore, m_disable_code: 800000

Total running seconds:  0.11
Thank you for using Medit.
-- FESpace: Nb of Nodes 6677 Nb of DoF 20031
-- FESpace: Nb of Nodes 6677 Nb of DoF 6677
Lambda = -9.81
Mu = 81.3953
-- Solve :
min -0.128858  max 0.0740267
Min/Max displacement x = -0.0303418 0.0513578
Min/Max displacement y = -0.128858 -6.10332e-61
Min/Max displacement z = -0.0747728 0.0740267
Min/Max displacement x = -0.0303418 0.0513578
Min/Max displacement y = -0.128858 -6.10332e-61
Min/Max displacement z = -0.0747728 0.0740267
-- FESpace: Nb of Nodes 6677 Nb of DoF 20031
Plot::  Sorry no ps version for this type of plot 7
FALLBACK (log once): Fallback to SW vertex for line stipple
FALLBACK (log once): Fallback to SW vertex processing, m_disable_code: 2000
FALLBACK (log once): Fallback to SW vertex processing in drawCore, m_disable_code: 2000
Plot::  Sorry no ps version for this type of plot 5
version de medit ffmedit -popen -filebin 1  displacement
-- Medit,  Release 3.0a (Nov. 30, 2007)
compiled:  Ven 7 jui 2024 11:36:37 CEST (with ff++ 4.14).

medit with binary version of popen : mesh(es)
mesh_name= displacement
End of mesh
Input seconds:     0.01

medit1()

Building scene(s)
Creating scene 1
Scene seconds:     0.02

Rendering scene(s)

Total running seconds:  0.17
Thank you for using Medit.
Plot::  Sorry no ps version for this type of plot 5
version de medit ffmedit -popen -filebin -addsol 1  utotal
-- Medit,  Release 3.0a (Nov. 30, 2007)
compiled:  Ven 7 jui 2024 11:36:37 CEST (with ff++ 4.14).

medit with binary version of popen : mesh(es) and solution(s)
mesh_name= utotal