# 3D elasticity problem with different Young's Modulus in the field

**URL:** <https://community.freefem.org/t/3d-elasticity-problem-with-different-youngs-modulus-in-the-field/2753>\
**Category:** General Discussion\
**Created:** [October 18, 2023, 12:39pm UTC](https://community.freefem.org/t/3d-elasticity-problem-with-different-youngs-modulus-in-the-field/2753 "2023-10-18T12:39:15Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![Gooyoon](https://avatars.discourse-cdn.com/v4/letter/g/ba9def/32.png) [@Gooyoon](https://community.freefem.org/u/Gooyoon)\
**Post date:** [October 18, 2023, 12:39pm UTC](https://community.freefem.org/t/3d-elasticity-problem-with-different-youngs-modulus-in-the-field/2753/1 "2023-10-18T12:39:15Z")

</div>

Dear FreeFem++ developers

I’m really new to FreeFem++

I want to proceed with simulation in which young’s modulus transforms other objects by dividing them into three parts as shown in the picture below, how should I perform?

I tried to apply modulus differently depending on the label by dividing the label, but I failed.  
I hope the block on the left has a high modulus so that it doesn’t deform even if it’s strong and stays in its original shape.  
I want that middle block can only deform.

 ![1-1](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/c/c8d26e809b33516cd4a2f8501cb1c584640d2ebd.png)  
 ![1-2](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/3/3c64fa7c0b3d67c2e08919c086dcfbdc070b456b.png)

Below is the code I wrote.

```auto
load "msh3"

// Mesh Generation

real S1 = 1;
real H1 = 1;

int move = 1;
int fixed = 10;

border a1(t=0., 1.){x=0; y=(H1-(2*H1)*t)*1.e-3; label = 1;}
border a2(t=0., 1.){x=(S1*t)*1.e-3; y=(-H1)*1.e-3; label = 1;}
border a3(t=0., 1.){x=(S1)*1.e-3; y=(-H1-(-2*H1)*t)*1.e-3; label = 1;}
border a4(t=0., 1.){x=(S1+(-S1)*t)*1.e-3; y=(H1)*1.e-3; label = 1;}

border b1(t=0., 1.){x=S1*1.e-3; y=(H1-(2*H1)*t)*1.e-3; label = 1;}
border b2(t=0., 1.){x=(S1+(S1*t))*1.e-3; y=(-H1)*1.e-3; label = 1;}
border b3(t=0., 1.){x=(2*S1)*1.e-3; y=(-H1-(-2*H1)*t)*1.e-3; label = 1;}
border b4(t=0., 1.){x=(2*S1+(-S1)*t)*1.e-3; y=(H1)*1.e-3; label = 1;}

border c1(t=0., 1.){x=2*S1*1.e-3; y=(H1-(2*H1)*t)*1.e-3; label = 1;}
border c2(t=0., 1.){x=(2*S1+(S1*t))*1.e-3; y=(-H1)*1.e-3; label = 1;}
border c3(t=0., 1.){x=(3*S1)*1.e-3; y=(-H1-(-2*H1)*t)*1.e-3; label = 1;}
border c4(t=0., 1.){x=(3*S1+(-S1)*t)*1.e-3; y=(H1)*1.e-3; label = 1;}

mesh ThS1 = buildmesh (a1(10) + a2(5) + a3(10) + a4(5));
mesh ThS2 = buildmesh (b1(10) + b2(5) + b3(10) + b4(5));
mesh ThS3 = buildmesh (c1(10) + c2(5) + c3(10) + c4(5));

func zmin = -0.1*1.e-3;
func zmax = 0.1*1.e-3;
int MaxLayer = 1;

int[int] rmid = [1, 1];
int[int] rmid2 = [1, 20];
int[int] rmidfix = [1, 10];

int[int] rup = [0, 1];
int[int] rup2 = [0, 20];
int[int] rupfix = [0, 10];

int[int] rdown = [0, 1];
int[int] rdown2 = [0, 20];
int[int] rdownfix = [0, 10];

mesh3 Th1 = buildlayers(ThS1, MaxLayer, zbound=[zmin, zmax], labelmid=rmid, labelup=rupfix, labeldown=rdown);
mesh3 Th2 = buildlayers(ThS2, MaxLayer, zbound=[zmin, zmax], labelmid=rmid, labelup=rup, labeldown=rdown);
mesh3 Th3 = buildlayers(ThS3, MaxLayer, zbound=[zmin, zmax], labelmid=rmid2, labelup=rup2, labeldown=rdown2);

mesh3 Th = Th1 + Th2 + Th3;

plot (Th, wait=1);

// Define Values

real Fx;
real Fy;
real Fz = 200.;

real ym1 = 10.; // Young's Modulus 1
real pr1 = 0.; // Poisson ratio 1
real ym2 = 1.; // Young's Modulus 2
real pr2 = 0.; // Poisson ratio 2 

real Mu1 = ym1/2*(1+pr1);
real Lambda1 = ym1*pr1/((1+pr1)*(1-2*pr1));

real Mu2 = ym2/2*(1+pr2);
real Lambda2 = ym2*pr2/((1+pr2)*(1-2*pr2));

// Fespace

fespace Vh(Th,P1);
Vh u1,u2,u3,u4;
Vh v1,v2,v3;
Vh lambda = (label = 1) ? Lambda1 : Lambda2; // If label = 1 : Low modulus
Vh mu = (label = 1) ? Mu1 : Mu2; // If label = 1 : Low modulus

// 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] // EOM
macro div(u1,u2,u3) ( dx(u1)+dy(u2)+dz(u3) ) // EOM

// Problem

problem Lame([u1,u2,u3],[v1,v2,v3])=
  int3d(Th)(  
	    lambda*div(u1,u2,u3)*div(v1,v2,v3)	
	    +2.*mu*( epsilon(u1,u2,u3)'*epsilon(v1,v2,v3) ) //')
	      )
  - int3d(Th) (
            [Fx, Fy, Fz]' * [v1, v2, v3]
        )
    + on(fixed,u1=0,u2=0,u3=0)
  ;

// Solve

Lame;

u4 = u4 + u3;

Th=movemesh(Th,[x+u1,y+u2,z+u3]);
u1=u1; u2=u2; u3=u3;

plot(Th, prev=1, wait=1 );

```

Thank you for your cooperation.

Best regards,  
Gooyoon

---

<div class="post-metadata">

**Author:** ![vlaude](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/vlaude/32/939_2.png) [@vlaude](https://community.freefem.org/u/vlaude)\
**Post date:** [October 19, 2023, 8:10am UTC](https://community.freefem.org/t/3d-elasticity-problem-with-different-youngs-modulus-in-the-field/2753/2 "2023-10-19T08:10:16Z")

</div>

Dear Gooyoon,  
Maybe you can try to change the logical test in the definition of lambda and mu (the Lamé constants), to use a condition on coordinates instead of label. I am not sure if I understand correctly your code, but S1 and H1 define the three different regions in the mesh? Then one of the logical conditions could be (x\<S1).  
Furthermore, you may use functions (func) to define the Lamé constants. Now you define them as finite element functions (of type Vh, with P1 FE functions). This choice has implications for the computation of integrals with int3d(Th)(). You can also define a P0 finite element space especially for them, if they are piecewise constants.  
Vincent

---

<div class="post-metadata">

**Author:** ![frederichecht](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/frederichecht/32/15_2.png) [@frederichecht](https://community.freefem.org/u/frederichecht)\
**Post date:** [October 19, 2023, 9:34am UTC](https://community.freefem.org/t/3d-elasticity-problem-with-different-youngs-modulus-in-the-field/2753/3 "2023-10-19T09:34:47Z")

</div>

Remark, in FreeFEM the label are only for boundary condition, you can use region name

and use P0 finite element of define lambda and mu

the correction of you script:

[Gooyoon.edp](https://community.freefem.org/uploads/short-url/kLDjC3WAKQn0OMdrXC0AbFIDdfV.edp) (3.2 KB)

---

<div class="post-metadata">

**Author:** ![Gooyoon](https://avatars.discourse-cdn.com/v4/letter/g/ba9def/32.png) [@Gooyoon](https://community.freefem.org/u/Gooyoon)\
**Post date:** [October 20, 2023, 6:10am UTC](https://community.freefem.org/t/3d-elasticity-problem-with-different-youngs-modulus-in-the-field/2753/4 "2023-10-20T06:10:08Z")

</div>

Dear Vincent

Thank you for your answer.

According to the professor’s answer below, I should define a P0 finite element space especially for the constants (the Lamé constants).  
And I think I should use the regions instead of the labels to designate the section.  
So I want to designate (region = 1) as low modulus and (region = 10) as high modulus.

Here is the part of the code I defined different constants depend on the regions.

```auto
fespace Ph(Th,P0);

Ph lambda = (region = 10) ? Lambda1 : Lambda2; // If region = 10 : High modulus
Ph mu = (region = 10) ? Mu1 : Mu2; // If region = 10 : High modulus

```

However, the Lamé constants specified in (region = 10) does not apply.  
I don’t know why…

Here is the code I ran.

[Gooyoon.edp](https://community.freefem.org/uploads/short-url/1cvSlL5VD5djRRJ1qnGOPqy1CdF.edp) (3.0 KB)

Gooyoon

---

<div class="post-metadata">

**Author:** ![Gooyoon](https://avatars.discourse-cdn.com/v4/letter/g/ba9def/32.png) [@Gooyoon](https://community.freefem.org/u/Gooyoon)\
**Post date:** [October 20, 2023, 6:36am UTC](https://community.freefem.org/t/3d-elasticity-problem-with-different-youngs-modulus-in-the-field/2753/5 "2023-10-20T06:36:12Z")

</div>

Thank you for your answer, professor!  
I’m sorry for the late reply. I’ve tried it in many ways.

I understand I should define different constant values by **regions** instead of labels.

So I modified the code based on your advice and I changed the situation a little bit.  
But modulus setting according to region was not applied.  
How do I declare variables (the Lamé constants) according to region?

The result is at below.

 ![1-1](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/6/655f47b2b4536dd6bbb94331bb8ee7304d6ce430.png)  
 ![1-2](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/a/ae17a785a0992f35af1afda58ff171e39154349e.png)

**This is what I expect.**

 ![1-3](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/e/e650865db5e8e0ea6762ff137c4f8d85dd6b75b1.png)

```auto
load "msh3"

// Mesh Generation

real S1 = 1;
real H1 = 1;

int low = 1;
int high = 10;
int fixed = 15;

border a1(t=0., 1.){x=0; y=(H1-(2*H1)*t)*1.e-3; label = fixed;}
border a2(t=0., 1.){x=(S1*t)*1.e-3; y=(-H1)*1.e-3; label = high;}
border a3(t=0., 1.){x=(S1)*1.e-3; y=(-H1-(-2*H1)*t)*1.e-3; label = high;}
border a4(t=0., 1.){x=(S1+(-S1)*t)*1.e-3; y=(H1)*1.e-3; label = high;}

border b1(t=0., 1.){x=S1*1.e-3; y=(H1-(2*H1)*t)*1.e-3; label = high;}
border b2(t=0., 1.){x=(S1+(S1*t))*1.e-3; y=(-H1)*1.e-3; label = low;}
border b3(t=0., 1.){x=(2*S1)*1.e-3; y=(-H1-(-2*H1)*t)*1.e-3; label = high;}
border b4(t=0., 1.){x=(2*S1+(-S1)*t)*1.e-3; y=(H1)*1.e-3; label = low;}

border c1(t=0., 1.){x=2*S1*1.e-3; y=(H1-(2*H1)*t)*1.e-3; label = high;}
border c2(t=0., 1.){x=(2*S1+(S1*t))*1.e-3; y=(-H1)*1.e-3; label = high;}
border c3(t=0., 1.){x=(3*S1)*1.e-3; y=(-H1-(-2*H1)*t)*1.e-3; label = high;}
border c4(t=0., 1.){x=(3*S1+(-S1)*t)*1.e-3; y=(H1)*1.e-3; label = high;}

mesh ThS1 = buildmesh (a1(10) + a2(5) + a3(10) + a4(5),fixedborder=1); // to be sure that the point on boundary do not move 
mesh ThS2 = buildmesh (b1(10) + b2(5) + b3(10) + b4(5),fixedborder=1);
mesh ThS3 = buildmesh (c1(10) + c2(5) + c3(10) + c4(5),fixedborder=1);

func zmin = -0.1*1.e-3;
func zmax = 0.1*1.e-3;
int MaxLayer = 1;

int[int] rghigh = [0, 10];
int[int] rglow = [0, 1];

ThS1 = change(ThS1, region=rghigh);
ThS2 = change(ThS2, region=rglow);
ThS3 = change(ThS3, region=rghigh);

mesh3 Th1 = buildlayers(ThS1, MaxLayer, zbound=[zmin, zmax], region=rghigh);
mesh3 Th2 = buildlayers(ThS2, MaxLayer, zbound=[zmin, zmax], region=rglow);
mesh3 Th3 = buildlayers(ThS3, MaxLayer, zbound=[zmin, zmax], region=rghigh);

mesh3 Th = Th1 + Th2 + Th3;

cout << " regions ="<< regions(Th) << endl;
plot (Th, wait=1);

// Define Values

real Fx;
real Fy;
real Fz = 50.;

real ym1 = 10.; // Young's Modulus High
real pr1 = 0.; // Poisson ratio 1
real ym2 = 1.; // Young's Modulus Low
real pr2 = 0.; // Poisson ratio 2 

real Mu1 = ym1/2*(1+pr1);
real Lambda1 = ym1*pr1/((1+pr1)*(1-2*pr1));
real Mu2 = ym2/2*(1+pr2);
real Lambda2 = ym2*pr2/((1+pr2)*(1-2*pr2));

// Fespace

fespace Vh(Th,P1);
fespace Ph(Th,P0);
Vh u1,u2,u3,u4;
Vh v1,v2,v3;
Ph lambda = (region = 10) ? Lambda1 : Lambda2; // If region = 10 : High modulus
Ph mu = (region = 10) ? Mu1 : Mu2; // If region = 10 : High modulus

// 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] // EOM
macro div(u1,u2,u3) ( dx(u1)+dy(u2)+dz(u3) ) // EOM

// Problem

problem Lame([u1,u2,u3],[v1,v2,v3])=
  int3d(Th)(  
	    lambda*div(u1,u2,u3)*div(v1,v2,v3)	
	    +2.*mu*( epsilon(u1,u2,u3)'*epsilon(v1,v2,v3) ) //')
	      )
  - int3d(Th) (
            [Fx, Fy, Fz]' * [v1, v2, v3]
        )
    + on(fixed,u1=0,u2=0,u3=0)
  ;

// Solve

Lame;

Th=movemesh(Th,[x+u1,y+u2,z+u3]);
u1=u1; u2=u2; u3=u3;

plot(Th, prev=1, wait=1);

```

Here is the code I modified.

---

<div class="post-metadata">

**Author:** ![vlaude](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/vlaude/32/939_2.png) [@vlaude](https://community.freefem.org/u/vlaude)\
**Post date:** [October 20, 2023, 7:45am UTC](https://community.freefem.org/t/3d-elasticity-problem-with-different-youngs-modulus-in-the-field/2753/6 "2023-10-20T07:45:58Z")

</div>

Dear Gooyon,  
Just a quick answer: maybe you should correct the logical comparison test, the operator is “==” whereas “=” is for assignment.  
Please try:

> \*\*Ph lambda = (\*\ ***region ==** **10) ? Lambda1 : Lambda2; // If region = 10 : High modulus**  
> \*\*Ph mu = (\*\ ***region ==** **10) ? Mu1 : Mu2; // If region = 10 : High modulus**

Vincent

---

<div class="post-metadata">

**Author:** ![frederichecht](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/frederichecht/32/15_2.png) [@frederichecht](https://community.freefem.org/u/frederichecht)\
**Post date:** [October 20, 2023, 8:23am UTC](https://community.freefem.org/t/3d-elasticity-problem-with-different-youngs-modulus-in-the-field/2753/7 "2023-10-20T08:23:00Z")

</div>

The result with the correct derand a smaller Fz = 50., ici ok

Remember this is linear elasticity so the figure is wrong.

the code :

[Gooyoon.edp](https://community.freefem.org/uploads/short-url/tKy3sgJOD6MUl1uMIpJyn4JufHC.edp) (3.3 KB)

---

<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:** [October 20, 2023, 8:48am UTC](https://community.freefem.org/t/3d-elasticity-problem-with-different-youngs-modulus-in-the-field/2753/8 "2023-10-20T08:48:16Z")

</div>

Did you post the right code? The regions don’t appear to reflect the materials  
properties until you change as per earlier poster and Fz is still 50 creating  
extreme deformation.  
It may be easier to see by just dumping to cout,

```auto
Ph lambda = (region == 2) ? Lambda1 : Lambda2; // If label = 1 : Low modulus
Ph mu = (region == 2) ? Mu1 : Mu2; // If label = 1 : Low modulus
Vh l1=lambda;
plot(l1,fill=1,wait=1);
// Macro
cout<<" mu "<<mu[]<<endl;

```

---

<div class="post-metadata">

**Author:** ![frederichecht](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/frederichecht/32/15_2.png) [@frederichecht](https://community.freefem.org/u/frederichecht)\
**Post date:** [October 21, 2023, 1:34pm UTC](https://community.freefem.org/t/3d-elasticity-problem-with-different-youngs-modulus-in-the-field/2753/9 "2023-10-21T13:34:16Z")

</div>

I think so , but I can male mistake.

---

<div class="post-metadata">

**Author:** ![Gooyoon](https://avatars.discourse-cdn.com/v4/letter/g/ba9def/32.png) [@Gooyoon](https://community.freefem.org/u/Gooyoon)\
**Post date:** [November 3, 2023, 5:59am UTC](https://community.freefem.org/t/3d-elasticity-problem-with-different-youngs-modulus-in-the-field/2753/10 "2023-11-03T05:59:31Z")

</div>

All the problems have been resolved! It really worked  
Thank you for all your help. I really appreciated.

Best regards,

Gooyoon