# Apply force at a single node

**URL:** https://community.freefem.org/t/apply-force-at-a-single-node/909
**Category:** General Discussion
**Created:** [April 13, 2021, 12:58am UTC](https://community.freefem.org/t/apply-force-at-a-single-node/909 "2021-04-13T00:58:30Z")
**Posts on this page:** 9
**Page:** 1

<div class="post-metadata">

### Author: ![jmorvan](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/jmorvan/32/2928_2.png) [@jmorvan](https://community.freefem.org/u/jmorvan)
#### Post date: [April 13, 2021, 12:58am UTC](https://community.freefem.org/t/apply-force-at-a-single-node/909/1 "2021-04-13T00:58:30Z")

</div>

Hello everyone!

I’m a new freefem user. I would like to know if is it possible to apply a force at only one node instead of a whole border of the mesh in a 2D elasticity problem.

Thanks!

---

<div class="post-metadata">

### Author: ![fivanci](https://avatars.discourse-cdn.com/v4/letter/f/e95f7d/32.png) [@fivanci](https://community.freefem.org/u/fivanci)
#### Post date: [April 13, 2021, 5:31am UTC](https://community.freefem.org/t/apply-force-at-a-single-node/909/2 "2021-04-13T05:31:38Z")

</div>

Hi !  
Yes, it is possible, but, as far as I am aware, you need to do this ‘by hand’ - construct matrix and rhs explicitly with “varf” rather than using “solve” or “problem”. If you provide a MWE and point to what is your problem, it would be easier to check out what exactly do you need.

Below, I slightly modified MWE from the FreeFEM++ documentation (Mathematical models → Elasticity → Beam under gravity [Elasticity](https://doc.freefem.org/models/elasticity.html)). Elastic beam is attached to the wall, and it falls down under gravity. For illustration, I modified the forcing term in just one dof which corresponds to bottom right corner of the beam and add force in the gravity-opposite direction.

real E = 21.5, sigma = 0.29, gravity = -0.1;

mesh Th = square(20,4,[5\*x,y]); plot(Th,wait=1,cmm=“initial elastic beam mesh”);  
fespace Vh(Th,[P1,P1]);  
Vh [ux,uy], [vx,vy]; // [ux,uy] is the displacement field

int forceNode; // here, dof corresponding to uy at bottom right corner on mesh Th; subroutine to  
// extract the dof corresponding to y-component of displacement in the bottom right corner  
{  
varf bottomID([ux,uy],[vx,vy]) = on(1,ux=0,uy=1);  
varf rightID([ux,uy],[vx,vy]) = on(2,ux=0,uy=1);  
ux[] = bottomID(0,Vh,tgv=1); vx[] = rightID(0,Vh,tgv=1);  
for(int k=0; k\<Vh.ndof; ++k) if(ux[][k]\>1e-5 && vx[][k]\>1e-5){  
forceNode=k;  
break; }  
}

macro epsilon(u1, u2) [dx(u1), dy(u2), (dy(u1)+dx(u2))/sqrt(2.)] //  
macro div(u,v) (dx(u) + dy(v)) //

// Problem  
real mu = E/(2\*(1+sigma));  
real lambda = E_sigma/((1+sigma)_(1-2\*sigma));

varf elasticity([ux,uy],[vx,vy])  
= int2d(Th)(  
lambda\*div(vx,vy)\*div(ux,uy)  
+ 2._mu_( epsilon(vx,vy)'_epsilon(ux,uy) )  
)  
+ on(4,ux=1,uy=1);  
varf rhs([ux,uy],[vx,vy]) = int2d(Th)( gravity_vy );  
real[int] bcID = elasticity(0,Vh);

matrix A = elasticity(Vh,Vh,solver=UMFPACK); // system matrix  
real[int] b = rhs(0,Vh); // right hand side  
b = bcID ? 0 : b; // Dirichlet b.c.  
b[forceNode] += 0.2; // modifying forcing term in bottom right dof; try comment this line

ux[] = A^-1\*b;  
mesh ThNew = movemesh(Th, [x+ux,y+uy]); plot(Th,ThNew,wait=true,cmm=“deformed beam mesh”);

Try commenting the line ‘b[forceNode] += 0.2;’ in above code.  
Hope this helps,  
if

---

<div class="post-metadata">

### Author: ![jmorvan](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/jmorvan/32/2928_2.png) [@jmorvan](https://community.freefem.org/u/jmorvan)
#### Post date: [April 13, 2021, 2:04pm UTC](https://community.freefem.org/t/apply-force-at-a-single-node/909/3 "2021-04-13T14:04:55Z")

</div>

Hi fivanci!

Thank you so much for your reply. I’ve been trying to do this the last few days but I just couldn’t find a way. I will run the code you’ve sent in order to understand how it works. It was very helpfull. Thanks again!

Jorge Morvan

---

<div class="post-metadata">

### Author: ![frederichecht](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/frederichecht/32/15_2.png) [@frederichecht](https://community.freefem.org/u/frederichecht)
#### Post date: [April 14, 2021, 7:19am UTC](https://community.freefem.org/t/apply-force-at-a-single-node/909/4 "2021-04-14T07:19:13Z")

</div>

I fact take force on one point is same has take a Dirac force,

Add you have example of that with Poisson equation (Laplacian) in  
tutorial/Laplace-RHS-Dirac.edp (cf. [FreeFem-sources/Laplace-RHS-Dirac.edp at master · FreeFem/FreeFem-sources · GitHub](https://github.com/FreeFem/FreeFem-sources/blob/master/examples/tutorial/Laplace-RHS-Dirac.edp) )

---

<div class="post-metadata">

### Author: ![jmorvan](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/jmorvan/32/2928_2.png) [@jmorvan](https://community.freefem.org/u/jmorvan)
#### Post date: [April 14, 2021, 9:07pm UTC](https://community.freefem.org/t/apply-force-at-a-single-node/909/5 "2021-04-14T21:07:27Z")

</div>

Dear Frédéric,

Thanks for your reply. I tried to adapt the code you’ve sent to the elasticity problem, however I don’t know how to choose the force components. In this same example, I think you provide the position where the force is applied with xdelta and ydelta and its value with cdelta. It’s not clear to me how to chosse the direction of this force.

Thanks again,  
Jorge Morvan

---

<div class="post-metadata">

### Author: ![jmorvan](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/jmorvan/32/2928_2.png) [@jmorvan](https://community.freefem.org/u/jmorvan)
#### Post date: [April 14, 2021, 9:10pm UTC](https://community.freefem.org/t/apply-force-at-a-single-node/909/6 "2021-04-14T21:10:21Z")

</div>

Dear fivanci,

I ran you code and it worked for I needed. Thank you so much, it was really helpfull.

Jorge Morvan

---

<div class="post-metadata">

### Author: ![fivanci](https://avatars.discourse-cdn.com/v4/letter/f/e95f7d/32.png) [@fivanci](https://community.freefem.org/u/fivanci)
#### Post date: [April 15, 2021, 1:09am UTC](https://community.freefem.org/t/apply-force-at-a-single-node/909/7 "2021-04-15T01:09:28Z")

</div>

Glad it worked for you. Keep in mind though, this works only when the point you apply the extra force to is actual degree of freedom. If it is some random point in the mesh (not an actual dof), you should consider the approach with Dirichlet delta function as suggested by Prof. Hecht (with interpolation matrix). Of course, if the point matches with dof, then these two should be the same.

if

---

<div class="post-metadata">

### Author: ![RaMattoso](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/ramattoso/32/499_2.png) [@RaMattoso](https://community.freefem.org/u/RaMattoso)
#### Post date: [January 18, 2022, 7:38pm UTC](https://community.freefem.org/t/apply-force-at-a-single-node/909/8 "2022-01-18T19:38:19Z")

</div>

Dear Frédéric,

I’m using the code you shared as example.  
I have Dirac’s masses as source in Helmholtz problem.  
The thing is: during my optimization process I need to get (from the others involved matrices) only values related to the Dirac’s masses’ position.

Can I obtain this information from the interpolate matrix D or should I try solve this in a completely different form?

Thanks in advance!

---

<div class="post-metadata">

### Author: ![frederichecht](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/frederichecht/32/15_2.png) [@frederichecht](https://community.freefem.org/u/frederichecht)
#### Post date: [January 25, 2022, 3:50pm UTC](https://community.freefem.org/t/apply-force-at-a-single-node/909/9 "2022-01-25T15:50:39Z")

</div>

In the matrix D the position are coded in the matrix coef and in the row and column index.

If dirac mass is one a node \ell then the `matrix D = interpolate(Vh,xdelta,ydelta); //` the interpolation matrix of the position, the matrix D have only one column and D \_{i0}= \delta\_{i\ell} where \delta\_{ij} is the Kronecker symbol.  
Otherwise just understand what is interpolation matrix (value at some point).
