# Request help with a 3D model
**URL:** https://community.freefem.org/t/request-help-with-a-3d-model/3935
**Category:** General Discussion
**Created:** [May 24, 2025, 2:30pm UTC](https://community.freefem.org/t/request-help-with-a-3d-model/3935 "2025-05-24T14:30:45Z")
**Posts on this page:** 5
**Page:** 1
### Author:  [@EChase](https://community.freefem.org/u/EChase)
#### Post date: [May 24, 2025, 2:30pm UTC](https://community.freefem.org/t/request-help-with-a-3d-model/3935/1 "2025-05-24T14:30:45Z")
Good morning,
Please be patient with this request. I am trying to learn Freefem, and decided to try to convert the 2D example in the documentation (under the system of elasticity) into a 3d problem. I am stuck at the point where the boundary conditions are applied.
Here is my attempt at the model:
// Three dimensional plate model
// a trial model
//
// metric model kg,m,s unit system
load “msh3”
// r3 is the surface label vector
int[int] r3=[11,12,13,14,15,16];
// face y=0 is 11
// face x=1 is 12
// face y=1 is 13
// face x=0 is 14
// face z=0 is 15
// face z=1 is 16
int r5 = 5;
// r5 is the volume label
// Material Properties
real Yngs = 210e9; // Pa
real nu = 0.3;
real rho = 7750; // kg/m^3
real mu = Yngs/(2\*(1+nu)); // lame coefficient mu
real lambda = Yngs_nu/((1+nu)_(1-2\*nu)); // lame coefficient lambda
// Geometry
real Wdth=0.42; // width
real Lngth=0.42; // length
real Th=0.01; // Thickness
int NW = 20; // Grid resolution width
int NL = 20; // Grid resolution length
int NT = 5; // grid resolution thickness
// Meshing
mesh3 Sp = cube(NW,NL,NT, [Wdth_x, Lngth_y, Th\*z], label=r3, flags=3, region=r5);
plot(Sp);
// gravity type load
real f = -1;
// Define Fespace
fespace Vh(Sp, P2);
Vh u, v, w;
Vh uu, vv, ww;
// define the differential operators
real sqrt2=(2.)^0.5;
macro epsilon(u1,u2,u3) [dx(u1), dy(u2), (dy(u1)+dx(u2))/sqrt2, dz(u3), (dz(u2)+dy(u3))/sqrt2, (dz(u1)+dx(u2))/sqrt2] // diff op
macro div(u,v,w) (dx(u)+dy(v)+dz(w))
// Problem
solve lame([u, v, w], [uu, vv, ww])
= int3d(Sp)(
lambda\*(div(u, v, w) \* div(uu, vv, ww)
+ 2._mu \* ( epsilon(u,v,w)’ \* epsilon(uu, vv, ww) )
)
- int3d(Sp)(
f_ww
)
+ on(14, u=0)
+ on(14, v=0)
+ on(14, w=0)
;
plot([u, v, w],wait 1, coef=coef);
// move mesh
mesh SP2 = movemesh(Sp, [x+u_coef, y+v_coef, z+w\*coef]);
Freefem crashes at the on command, and I am stuck.
I would appreciate some assistance in getting this example model to work, and would greatly appreciate any advice in the model in general.
Thank you in advance for your time and assistance.
Ed
---
### Author:  [@fb77](https://community.freefem.org/u/fb77)
#### Post date: [May 24, 2025, 6:38pm UTC](https://community.freefem.org/t/request-help-with-a-3d-model/3935/2 "2025-05-24T18:38:33Z")
You just miss a closing parenthese for the term with lambda in factor.
---
### Author:  [@Monirul25](https://community.freefem.org/u/Monirul25)
#### Post date: [May 25, 2025, 3:15am UTC](https://community.freefem.org/t/request-help-with-a-3d-model/3935/3 "2025-05-25T03:15:46Z")
Here, is your Code with several corrections.
[FEM\_elasticity-3D.edp](https://community.freefem.org/uploads/short-url/e5eRenfpFozSKyGPW4sMN8x5Z5q.edp) (2.0 KB)
( **Now, error is coming higher as i have not included the source term in the problem. If you give that i hope error will come nice as exact solution and approximate solution are matching** ).
If you give me your equations or provide your source, i will check order of convergence.
Thanks in advance!.
---
### Author:  [@EChase](https://community.freefem.org/u/EChase)
#### Post date: [May 30, 2025, 12:11pm UTC](https://community.freefem.org/t/request-help-with-a-3d-model/3935/4 "2025-05-30T12:11:00Z")
Thank you both very much!
My equations are all written down on paper, and I need to formalize them. I will post the math and ask you to look at it. Thank you for making the offer.
Thank you
---
### Author:  [@EChase](https://community.freefem.org/u/EChase)
#### Post date: [June 24, 2025, 10:44am UTC](https://community.freefem.org/t/request-help-with-a-3d-model/3935/6 "2025-06-24T10:44:12Z")
Math Notes:
Grad u= ∇u= ∂u/∂x i+ ∂u/∂y j+ ∂u/∂z k
div u= ∇∙u= ∂u/∂x+ ∂u/∂y+ ∂u/∂z
curl u= ∇×u=((∂u\_3)/∂y- (∂u\_2)/∂z,(∂u\_3)/∂z- (∂u\_1)/∂x,(∂u\_2)/∂x- (∂u\_1)/∂y)
curl u= ∇×u= (∂/∂x,∂/∂y,∂/∂z)×(u\_1,u\_2,u\_3 )
∇×u= ((∂u\_3)/∂y- (∂u\_2)/∂z)i-((∂u\_3)/∂z- (∂u\_1)/∂x)j+((∂u\_2)/∂x- (∂u\_1)/∂y)k
Kronecker Delta=δ\_ij= ■(1 if i=j@0 otherwise)
Relationship between stress tensor and strain tensor
stress tensor= σ\_ij= λδ\_ij ∇∙u+2μϵ\_ij (u)
strain tensor= ϵ\_ij= 1/2 (〖∂u〗\_i/〖∂x〗\_j +〖∂u〗\_j/〖∂x〗\_i )
Lame constants
λ= Eν/(1+ν)(1-2ν)
μ=E/2(1+ν)
Where:
E = Young’s Modulus
= poisson’s ratio
Displacement formula : Navier-Lamé equations
(λ+μ) u\_(k,ki)+μu\_(i,kk)+F\_i=ρ (∂^2 u\_i)/〖∂t〗^2
(λ+μ)∇∇∙u+μ∇^2×u+F\_i=ρ (∂^2 u\_i)/〖∂t〗^2
(λ+μ)grad div u+μcurl curl u+F\_i=ρ (∂^2 u\_i)/〖∂t〗^2
For static conditions
ρ (∂^2 u\_i)/〖∂t〗^2 =0
FreeFEM does not use the Navier-Lame equation for a structural problem, because the associated variational form does not give the right boundary conditions. Instead:
-div(σ)=f in Ω
As FreeFEM solves the div value to begin with, the matrix definition of is required.
σ\_ij= λδ\_ij (∂u/∂x+ ∂u/∂y+ ∂u/∂z)+2/2 μ(〖∂u〗\_i/〖∂x〗\_j +〖∂u〗\_j/〖∂x〗\_i )
σ\_11= λ(〖∂u〗\_x/∂x)+2μ(〖∂u〗\_x/∂x)
σ\_12= μ(〖∂u〗\_x/∂y+〖∂u〗\_y/∂x)
σ\_13= μ(〖∂u〗\_x/∂z+〖∂u〗\_z/∂x)
σ\_22= λ( 〖∂u〗\_y/∂y)+2μ(〖∂u〗\_y/∂y)
σ\_23= μ(〖∂u〗\_y/∂z+〖∂u〗\_z/∂y)
σ\_33= λ( 〖∂u〗\_z/∂z)+2μ(〖∂u〗\_z/∂z)
These values are used in the div and epsilon matrices which are fed to the div and epsilon matrices which are fed to the solve function
epsilon=[■(〖∂u〗\_x/∂x&〖∂u〗\_x/∂y+〖∂u〗\_y/∂x&〖∂u〗\_x/∂z+〖∂u〗\_z/∂x@〖∂u〗\_x/∂y+〖∂u〗\_y/∂x&〖∂u〗\_y/∂y&〖∂u〗\_y/∂z+〖∂u〗\_z/∂y@〖∂u〗\_x/∂z+〖∂u〗\_z/∂x&〖∂u〗\_y/∂z+〖∂u〗\_z/∂y& 〖∂u〗\_z/∂z)]
div=[■(〖∂u〗\_x/∂x&0&0@0&〖∂u〗\_y/∂y&0@0&0& 〖∂u〗\_z/∂z)]
Note: The coefficients ,2 and are handles in the solve sequence. These matrices just separate the terms.
real sqrt2=(2.)^0.5;
macro epsilon(u1,u2,u3) [dx(u1), dy(u2), (dy(u1)+dx(u2))/sqrt2, dz(u3), (dz(u2)+dy(u3))/sqrt2, (dz(u1)+dx(u2))/sqrt2]
// epsilon(u1,u2)’\*epsilon(v1,v2) = epsilon(u): epsilon(v)
macro div(u,v,w) (dx(u)+dy(v)+dz(w))
and the solve function looks like this:
// Problem
solve lame([u, v, w], [uu, vv, ww])
= int3d(Sp)(lambda\*(div(u, v, w) \* div(uu, vv, ww)))
+ int3d(Sp)(2.\*mu \* ( epsilon(u,v,w)’ \* epsilon(uu, vv, ww)))
…
Minus forces
…
On boundary conditions
…
;
lame; // solve
These are the notes that go with the script. It is a simple cantilevered beam structural statics model. I was wondering if the definitions are correct, especially that for epsilon. Does the order of the terms matter?