# Applying Multiple Point Forces
**URL:** https://community.freefem.org/t/applying-multiple-point-forces/1792
**Category:** General Discussion
**Created:** [May 31, 2022, 11:29am UTC](https://community.freefem.org/t/applying-multiple-point-forces/1792 "2022-05-31T11:29:29Z")
**Posts on this page:** 9
**Page:** 1
### Author:  [@Stefanosv](https://community.freefem.org/u/Stefanosv)
#### Post date: [May 31, 2022, 11:29am UTC](https://community.freefem.org/t/applying-multiple-point-forces/1792/1 "2022-05-31T11:29:29Z")
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](https://doc.freefem.org/models/elasticity.html) 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](https://community.freefem.org/t/apply-force-at-a-single-node/909) 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.
---
### Author:  [@frederichecht](https://community.freefem.org/u/frederichecht)
#### Post date: [June 2, 2022, 9:10am UTC](https://community.freefem.org/t/applying-multiple-point-forces/1792/2 "2022-06-02T09:10:47Z")
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
```auto
// 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);
```
---
### Author:  [@Stefanosv](https://community.freefem.org/u/Stefanosv)
#### Post date: [June 2, 2022, 6:58pm UTC](https://community.freefem.org/t/applying-multiple-point-forces/1792/3 "2022-06-02T18:58:45Z")
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.
---
### Author:  [@Stefanosv](https://community.freefem.org/u/Stefanosv)
#### Post date: [July 13, 2022, 8:29am UTC](https://community.freefem.org/t/applying-multiple-point-forces/1792/4 "2022-07-13T08:29:22Z")
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).
```auto
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);
```
---
### Author:  [@frederichecht](https://community.freefem.org/u/frederichecht)
#### Post date: [July 15, 2022, 12:22pm UTC](https://community.freefem.org/t/applying-multiple-point-forces/1792/5 "2022-07-15T12:22:08Z")
I make a bug in the script,
Me I put a vectorial finite element
```auto
fespace Vh(th, [P1, P1]);
```
to 2 scalars finite element
```auto
fespace Vh(th, P1]);
```
so this code works:
```auto
// 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
```auto
// 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.
```
---
### Author:  [@Stefanosv](https://community.freefem.org/u/Stefanosv)
#### Post date: [July 15, 2022, 11:52pm UTC](https://community.freefem.org/t/applying-multiple-point-forces/1792/6 "2022-07-15T23:52:19Z")
The issue I have with using
```auto
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?
```auto
for(int i=0; i < K.n / 2; i++){
ux(i) = uu[](2*i);
uy(i) = vv[](2*i+1);
}
```
---
### Author:  [@frederichecht](https://community.freefem.org/u/frederichecht)
#### Post date: [July 17, 2022, 10:33am UTC](https://community.freefem.org/t/applying-multiple-point-forces/1792/7 "2022-07-17T10:33:16Z")
Yes but this is strongly dépendant of internat data structure, use:
```auto
ux=uu;
uy=vv;
```
---
### Author:  [@Stefanosv](https://community.freefem.org/u/Stefanosv)
#### Post date: [July 18, 2022, 10:41am UTC](https://community.freefem.org/t/applying-multiple-point-forces/1792/8 "2022-07-18T10:41:40Z")
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
---
### Author:  [@frederichecht](https://community.freefem.org/u/frederichecht)
#### Post date: [July 18, 2022, 12:49pm UTC](https://community.freefem.org/t/applying-multiple-point-forces/1792/9 "2022-07-18T12:49:00Z")
Yes because ux,uy must be finite element function,
```auto
fespace Wh(th, P2);
Wh ux,uy;
...
```
then you have traduct a vectorial Finite element function in scalar finite element function.