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,[8x,2y]);
//Parameters
real Rho = 8000.; //Density
real E1 = 210.e9; //Solid material Young modulus
real E0 = E11e-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 = ENu/((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.5t; label=Free;};
border b3(t=0., L){x=L-0.5t; y=0.5H; label=Free;};
border b4(t=0., H){x=0.5L; y=0.5H+0.5t; label=Free;};
border b5(t=0., L){x=0.5L-0.5t; 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 = E11e-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 = ENu/((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;