Home › Forum › Community Help › Programming with SOFA › Forcefield’s addDForce
 This topic has 3 replies, 2 voices, and was last updated 5 years, 1 month ago by Hugo.

AuthorPosts

5 August 2015 at 12 h 02 min #3444
Dear all,
Currently, the amount of force in a
VaccumShereForceField
is linearly proportional to the distance from the center of the sphere. I want to change it to a gravitylike force field (inversely proportional to the square of the distance) by changing the relevant addForce(…) and addDForce(…) functions. I have already changed the addForce(…) function successfully but do not know how to change addDForce(…) function. I would be grateful if you describe the underlying math for computing addDForce(…).Regards,
Sharif5 August 2015 at 16 h 20 min #3448Dear Sharif,
I will try to explain it in a short way. The use of
addForce
andaddDForce
depends on the integration scheme you are using (implicit or explicit). In the explicit case, you solve the following system: M a(t+dt) = f(t) from the Newton’s law. If we solve the equation regarding the velocity v, we have: M ∆v = dt * f(t). In SOFA, the right hand side f(t) is computed from the
addForce
function.  In the implicit case, you solve the following system: M a(t+dt) = f(t+dt) and this changes the computation since f(t+dt) is unknown. Following the Taylor series, you have:
f(t+dt) = f(t) + (df(t)/dx) * ∆x
f(t+dt) = f(t) + (df(t)/dx) * dt*(v+∆v)
f(t+dt) = f(t) + K * dt*(v+∆v).
In the implicit case, we then have:
(M – dt*dt * K) ∆v = dt* f(t) + dt*dt* K*v
A ∆v = b
In SOFA, f(t) is computed from theaddForce
function and theaddDForce
computes the term dt*K, used for both right and lefthand side.
Therefore:
 in case of an explicit scheme, only the
addForce
is called;  in case of an implicit scheme, both
addForce
andaddDForce
are called.
Now, how to code the
addDForce
from youraddForce
?
Basically, theaddDForce
term is the derivative of theaddForce
term, regarding the DOF.
Let’s assume that the addForce is : f(t) = x² where x is you DOF (position).
Then, the addDForce is: df(t) = 2 * x * dx * dt.
This should be coded:for (i = 0 ; i lowerthan nbPoints; i++) df[i] += 2 * x[i] * dx[i] * mparams.kFactor;
Hope this helps and do not hesitate to ask for more.
Cheers,Hugo
5 August 2015 at 17 h 30 min #3459Dear Talbot,
Thank you for your comprehensive response. Can you please explain more on the differentiation procedure? What should we do if our force field depends on location as a vector? In “VaccumSphereForceField”, force is computed according to the following (bold variables are vectors):
f = (dp/dp)*(dpr)*stiffnessv*(dpr)*damping
f
is the force vector,dp
the position vector of a particle relative to the center of the field,r
is the radius of the sphere where the forcefield applies,stiffness
is a userinput parameter specifying the intensity of the field,v
is the velocity of the particle, anddamping
is another userinput parameter which determines the viscosity of the field. As it is seen, our force field is best described by the following:f = f(dp, v)
According to what you said, we should compute
df/dt * dx(t+dt)
foraddDForce
. How is it when we deal with vectors? Dogradians
apply in this case?As another question, what is
mparams.kFactor
and how is it set in SOFA?As a suggestion, it would be perfect if this forum supports mathematical notations. It is hard to be precise with mathematics using normal notations!
Best regards,
Sharif20 August 2015 at 13 h 01 min #3485Dear Sharif,
No according to what I wrote, you should compute
df(t)/dx * dx(t+dt)
forAddDForce
. It means you need to derive the expression f(dp,v) regarding the position of the particles dp.kFactor is set by the time integration scheme, e.g. the EulerImplicit solver. this kFactor corresponds to dt the time step.
About the forum, I know it is complicated to write math equations but unfortunately we did not find a better solution for the moment. We will try to improve this.
Best regards,
Hugo
 In the explicit case, you solve the following system: M a(t+dt) = f(t) from the Newton’s law. If we solve the equation regarding the velocity v, we have: M ∆v = dt * f(t). In SOFA, the right hand side f(t) is computed from the

AuthorPosts
 You must be logged in to reply to this topic.