3d elasticity problem

General Discussion
1

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
load "msh3"
load "medit"
load "gmsh"

// 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);
2

My correction trivial,

zoumpou.edp (3.1 KB)

the output:


cd /Users/hecht/Downloads ; FreeFem++-CoCoa /Users/hecht/Downloads/zoumpou.edp  
brochet:~ hecht$ cd /Users/hecht/Downloads ; FreeFem++-CoCoa /Users/hecht/Downloads/zoumpou.edp  
do script " cd '/Users/hecht/Downloads'; /Applications/FreeFem++.app/Contents/ff-4.14-dev/bin/FreeFem++ 'zoumpou.edp' -glut '/Applications/FreeFem++.app/Contents/ff-4.14-dev/bin/ffglut' ;"
-- FreeFem++ v4.14 (Ven  7 jui 2024 11:35:37 CEST - git v4.14-60-g660dbd36)
   file : zoumpou.edp
 Load: lg_fem lg_mesh lg_mesh3 eigenvalue 
    1 : //// Load necessary FreeFem++ libraries
    2 : load "msh3"
    3 : load "medit"
    4 : load "gmsh"
    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)
     Copyright (c) LJLL, 1999-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
 Loading data file(s)
End of mesh
  Input seconds:     0.01

 medit1() 

 Building scene(s)
  Creating scene 1
   Loading default options
  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)
     Copyright (c) LJLL, 1999-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
 Loading data file(s)
End of mesh
  Input seconds:     0.01

 medit1() 

 Building scene(s)
  Creating scene 1
   Loading default options
  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)
     Copyright (c) LJLL, 1999-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
 Loading data file(s)
End of mesh
.sol: Dimension 3 (mesh)3 (lecture)1 
  Input seconds:     0.01

 medit1() 


 Building scene(s)
  Creating scene 1
   Loading default options
  Scene seconds:     0.02