Applying Multiple Point Forces

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

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);
// 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);

1 Like

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.