Hello,
I’m a new Free Fem user and this is my first dive to programming after using MatLAB for my university.
I am studying about topology optimization and I want to use the Elasticity module to extract values for Displacement that will be used on my optimization algorithm.
The issue I’m having is that I will need to apply concentrated forces on single nodes. I’ve tried to use the answer from this post but I can’t seem to be able to change the value of the load or the node where it is applied. Most of my problems come with lack of experience of using FreeFem and also the lack of understanding of the commands.
Could someone guide me with explaining how the code linked above “works” and give me some documentation about the commands being used?
Thanks in advance.
in 2d, If you want to add Dirac force on a node it is more simple, but you need to find the
2 dof associated to the node( call dof1,dof2).
This is the tricky part in this example I add a point in the mesh to be sure
after you have just to add the the right hand side a vector b such that
b=0;
b[dof1]=f1;
b[dof2]=f2;
bee careful with the signe of force ( in problem put a minus) because
we solve " A u + b = 0" not A u = b .
where [f1, f2] is the forces
// Beam under it own weight
// Parameters
real E = 21.5;
real sigma = 0.29;
real gravity = -0.05;
// Mesh
int bottombeam = 2;
border a(t=2, 0){x=0; y=t ;label=1;} // left beam
border b(t=0, 10){x=t; y=0 ;label=bottombeam;} // bottom of beam
border c(t=0, 2){x=10; y=t ;label=1;} // rigth beam
border d(t=0, 10){x=10-t; y=2; label=3;} // top beam
real xp = 5,yp=1;
real[int,int] P=[[xp,yp]]; // add point (xp,yp) in a mesh ..
mesh th = buildmesh(b(20) + c(5) + d(20) + a(5),points=P,fixeborder=1);
plot(th,wait=1);
// Fespace
fespace Vh(th, [P1, P1]);
Vh [uu, vv], [w, s];
// find midlle (xp,yp) dof
w[]=0:w[].n-1; // put dof number in array
int dof1=w(xp,yp)+0.5;
int dof2=s(xp,yp)+0.5;
assert( abs(dof1-w(xp,yp)) + abs(dof2-s(xp,yp)) < 1e-6); // the point (xp,yp ) is a node point
cout << " dof "<< dof1 << " "<< dof2 << endl;
real[int] brhs(Vh.ndof); // the add rhs ..
brhs = 0;
brhs[dof2]= - 100000;
// Macro
real sqrt2 = sqrt(2.);
macro epsilon(u1, u2) [dx(u1), dy(u2), (dy(u1) + dx(u2))/sqrt2] // EOM
macro div(u, v) (dx(u) + dy(v)) // EOM
// Problem (with solve)
real mu = E/(2*(1 + sigma));
real lambda = E*sigma/((1 + sigma)*(1 - 2*sigma));
cout << "Lambda = " << lambda << endl;
cout << "Mu = " << mu << endl;
cout << "Gravity = " << gravity << endl;
solve Elasticity ([uu, vv], [w, s], solver=sparsesolver)
= int2d(th)(
lambda*div(w, s)*div(uu, vv)
+ 2.*mu*(epsilon(w, s)' * epsilon(uu, vv))
)
//- int2d(th)(gravity*s)
+ brhs
+ on(1, uu=0, vv=0)
;
cout << "Max displacement = " << uu[].linfty << endl;
// Plot
plot([uu, vv], wait=1);
plot([uu, vv], wait=1, bb=[[-0.5, 2.5], [2.5, -0.5]]);
// Move mesh
mesh th1 = movemesh(th, [x+uu, y+vv]);
plot(th1, wait=1);
I can’t thank you enough for your post. Your help was invaluable and your notes next to the code were really helpful for understanding the procedure.
Sorry for re-opening the topic but I’ve noticed an issue. By using P1 to create the fespace Vh the values of displacement uu and vv (displacement on x and y) are equal which is false. Now if I use P2 to create the fespace the displacement come out correctly but the direct force is not applied in the correct place (middle of the beam).
real E = 21.5;
real sigma = 0.29;
real gravity = 0;
// Mesh
int bottombeam = 2;
border a(t=2, 0){x=0; y=t ;label=1;} // left beam
border b(t=0, 10){x=t; y=0 ;label=bottombeam;} // bottom of beam
border c(t=0, 2){x=10; y=t ;label=1;} // rigth beam
border d(t=0, 10){x=10-t; y=2; label=3;} // top beam
real xp = 5,yp=1;
real[int,int] P=[[xp,yp]]; // add point (xp,yp) in a mesh ..
mesh th = buildmesh(b(20) + c(5) + d(20) + a(5),points=P,fixeborder=1);
plot(th,wait=1);
// Fespace
fespace Vh(th, P2);
Vh uu, vv;
Vh w, s;
// find midlle (xp,yp) dof
w[]=0:w[].n-1; // put dof number in array
int dof1=w(xp,yp)+0.5;
int dof2=s(xp,yp)+0.5;
assert( abs(dof1-w(xp,yp)) + abs(dof2-s(xp,yp)) < 1e-6);
cout << " dof "<< dof1 << " "<< dof2 << endl;
real[int] brhs(Vh.ndof); // the add rhs ..
brhs = 0;
brhs[dof2]= 10;
//brhs[dof1]= 5;
// Macro
real sqrt2 = sqrt(2.);
macro epsilon(u1, u2) [dx(u1), dy(u2), (dy(u1) + dx(u2))/sqrt2] // EOM
macro div(u, v) (dx(u) + dy(v)) // EOM
// Problem (with solve)
real mu = E/(2*(1 + sigma));
real lambda = E*sigma/((1 + sigma)*(1 - 2*sigma));
cout << "Gravity = " << gravity << endl;
solve Elasticity ([uu, vv], [w, s], solver=sparsesolver)
= int2d(th)(
lambda*div(w, s)*div(uu, vv)
+ 2.*mu*(epsilon(w, s)' * epsilon(uu, vv))
)
- int2d(th)(gravity*s)
+ brhs
+ on(1, uu=0, vv=0)
;
// Output
real dxmin = uu[].min;
real dymin = vv[].min;
cout << " - dep. max x = "<< dxmin << " y=" << dymin << endl;
//Stiffness Matrix
varf bin(uu,vv) = int2d(th)(dx(uu) * dx(vv) + dy(uu) * dy(vv));
matrix K = bin(Vh, Vh);
real[int] U(K.n);
for(int i=0; i < K.n / 2; i++){
U(2*i) = uu[](i);
U(2*i+1) = vv[](i);
}
real[int] F(K.n);
F = K * U;
real[int] f1(K.n / 2);
real[int] f2(K.n / 2);
for(int i=0; i < K.n / 2; i++){
f1(i) = F(2*i);
f2(i) = F(2*i+1);
}
{ofstream fout("U.txt");
fout << U << endl;
}
{ofstream fout("f1.txt");
fout << f1 << endl;
}
{ofstream fout("f2.txt");
fout << f2 << endl;
}
cout << "Max f = " << f1.linfty << endl;
{ofstream fout("Stiffness_matrix.txt");
fout << K << endl;
}
//Plot
plot([uu, vv], wait=1);
//plot([uu, vv], wait=1, bb=[[-0.5, 2.5], [2.5, -0.5]]);
// Move mesh
mesh th1 = movemesh(th, [x+uu, y+vv]);
plot(th1, wait=1);
I make a bug in the script,
Me I put a vectorial finite element
fespace Vh(th, [P1, P1]);
to 2 scalars finite element
fespace Vh(th, P1]);
so this code works:
// find midlle (xp,yp) dof
w[]=0:w[].n-1; // put dof number in array
int dof1=w(xp,yp)+0.5;
int dof2=s(xp,yp)+0.5;
assert( abs(dof1-w(xp,yp)) + abs(dof2-s(xp,yp)) < 1e-6);
cout << " dof "<< dof1 << " "<< dof2 << endl;
otherwise dof1 and dof2 are the same number and I put no explanation how the coupling for inite element space is build in problem.
Please do
// Fespace
fespace Vh(th, [P2,P2]);
Vh [uu, vv];
Vh [w, s];
// find midlle (xp,yp) dof
w[]=0:w[].n-1; // put dof number in array
int dof1=w(xp,yp)+0.5;
int dof2=s(xp,yp)+0.5;
assert( abs(dof1-w(xp,yp)) + abs(dof2-s(xp,yp)) < 1e-6);
cout << " dof "<< dof1 << " "<< dof2 << endl;
real[int] brhs(Vh.ndof); // the add rhs ..
brhs = 0;
brhs[dof2]= 10; // vertical force
//brhs[dof1]= 5;// horizontal force here 0.
The issue I have with using
fespace Vh(th, [P1, P1]);
or
fespace Vh(th, [P2, P2]);
is that when I export the values of displacement on the axis x and y (the values uu and vv) they are equal which is not possible. The bizarre thing is that ploting with [x+uu, y+vv] is correct which makes no sense if uu and vv the same.
Do I have to run something like this to seperate the values of ux and uy from uu and vv?
for(int i=0; i < K.n / 2; i++){
ux(i) = uu[](2*i);
uy(i) = vv[](2*i+1);
}
Yes but this is strongly dépendant of internat data structure, use:
ux=uu;
uy=vv;
By using the ux=uu; and uy=vv; I get the error.
Exec error : Try to get unset x,y, …
– number :1
And I’m not really sure how that will help the issue cause seeing the values of uu and vv you realize that uu=vv
Yes because ux,uy must be finite element function,
fespace Wh(th, P2);
Wh ux,uy;
...
then you have traduct a vectorial Finite element function in scalar finite element function.