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!
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). 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 = Esigma/((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)( gravityvy );
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
3 Likes
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
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 )
1 Like
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
Dear fivanci,
I ran you code and it worked for I needed. Thank you so much, it was really helpfull.
Jorge Morvan
1 Like
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
2 Likes
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!
1 Like
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).