# Topology optimization based on variable density method
**URL:** https://community.freefem.org/t/topology-optimization-based-on-variable-density-method/1969
**Category:** General Discussion
**Created:** [August 29, 2022, 7:18am UTC](https://community.freefem.org/t/topology-optimization-based-on-variable-density-method/1969 "2022-08-29T07:18:09Z")
**Posts on this page:** 3
**Page:** 1
### Author:  [@Yangshen](https://community.freefem.org/u/Yangshen)
#### Post date: [August 29, 2022, 7:18am UTC](https://community.freefem.org/t/topology-optimization-based-on-variable-density-method/1969/1 "2022-08-29T07:18:09Z")
Dear all,
I want to apply an initial density field with freefem, which works when my design domain is rectangular. And the resulting density matrix is a matrix related to the number of nodes. The code is as follows.
mesh Sh=square(4,5,[8_x,2_y]);
//Parameters
real Rho = 8000.; //Density
real E1 = 210.e9; //Solid material Young modulus
real E0 = E1_1e-9; //empty material Young modulus
real Nu = 0.27; //Poisson ratio
real Gravity = -9.81; //Gravity
fespace Vh(Sh,[P1,P1]);
Vh[ux,uy],[vx,vy];
fespace Vh1(Sh,P1);
Vh1 E,theta;
theta = 1; //Define initial density
E = theta_(E1-E0)+E0;
//Macro
real sqrt2 = sqrt(2.);
macro Epsilon(ux, uy) [dx(ux), dy(uy), (dy(ux)+dx(uy))/sqrt2] //
macro Divergence(ux, uy) (dx(ux) + dy(uy)) //
//Problem
real Mu = E/(2._(1.+Nu));
real Lambda = E_Nu/((1.+ Nu)\*(1.-2.\*Nu));
varf vElasticity ([ux,uy], [vx, vy])
= int2d(Sh)(
Lambda \* Divergence(vx, vy) \* Divergence(ux, uy)
+ 2. \* Mu \* (
Epsilon(vx, vy)’ \* Epsilon(ux, uy)
)
)
+ int2d(Sh)(
Rho \* Gravity \* vy
)
+ on(4, ux=0, uy=0)
;
plot(theta,fill=1);
cout \<\<“u=”\<\< theta \<\< endl;
Node density matrix
It seems that freefem is obtained by applying density to nodes, and cell density is obtained by interpolation of node density. If I change the design domain to an L-beam, it doesn’t seem to work. I don’t know why.
real L = 20.; //Beam length
real H = 10.; //Beam height
int Fixed = 1; //Beam fixed label
int Free = 2; //Beam free label
border b1(t=0., L){x=t; y=0.; label=Free;};
border b2(t=0., H){x=L; y=0.5_t; label=Free;};
border b3(t=0., L){x=L-0.5_t; y=0.5_H; label=Free;};
border b4(t=0., H){x=0.5_L; y=0.5_H+0.5_t; label=Free;};
border b5(t=0., L){x=0.5_L-0.5_t; y=H; label=Free;};
border b6(t=0., H){x=0.; y=H-t; label=Fixed;};
mesh Sh = buildmesh(b1(10) + b2(10) + b3(10) + b4(10) + b5(10) + b6(10));
//Parameters
real Rho = 8000.; //Density
real E1 = 210.e9; //Young modulus
real E0 = E1_1e-9;
real Nu = 0.27; //Poisson ratio
real Gravity = -9.81; //Gravity
fespace Vh(Sh,[P1,P1]);
Vh[ux,uy],[vx,vy];
fespace Vh1(Sh,P1);
Vh1 E,theta;
theta = 1;
E = theta_(E1-E0)+E0;
//Macro
real sqrt2 = sqrt(2.);
macro Epsilon(ux, uy) [dx(ux), dy(uy), (dy(ux)+dx(uy))/sqrt2] //
macro Divergence(ux, uy) (dx(ux) + dy(uy)) //
//Problem
real Mu = E/(2._(1.+Nu));
real Lambda = E_Nu/((1.+ Nu)\*(1.-2.\*Nu));
varf vElasticity ([ux,uy], [vx, vy])
= int2d(Sh)(
Lambda \* Divergence(vx, vy) \* Divergence(ux, uy)
+ 2. \* Mu \* (
Epsilon(vx, vy)’ \* Epsilon(ux, uy)
)
)
+ int2d(Sh)(
Rho \* Gravity \* vy
)
+ on(Fixed, ux=0, uy=0)
;
plot(Sh,wait=1);
cout \<\<“u=”\<\< theta \<\< endl;
---
### Author:  [@frederichecht](https://community.freefem.org/u/frederichecht)
#### Post date: [September 9, 2022, 8:50am UTC](https://community.freefem.org/t/topology-optimization-based-on-variable-density-method/1969/2 "2022-09-09T08:50:25Z")
theta is a fe function so R^2\mapsto R not a function form R \mapsto R
> [@Yangshen](#):
>
> _E = theta_(E1-E0)+E0;
FreeFEM say no error, but this have no sens, sorry.
---
### Author:  [@Yangshen](https://community.freefem.org/u/Yangshen)
#### Post date: [September 9, 2022, 10:18am UTC](https://community.freefem.org/t/topology-optimization-based-on-variable-density-method/1969/3 "2022-09-09T10:18:25Z")
Thank you for your reply.