[SOLVED] Forcefield’s addDForce

[sharif]

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 gravity-like 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,
Sharif

[Hugo]

Dear Sharif,

I will try to explain it in a short way. The use of addForce and addDForce 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 the addForce function and the addDForce computes the term dt*K, used for both right- and left-hand side.

Therefore:

• in case of an explicit scheme, only the addForce is called;
• in case of an implicit scheme, both addForce and addDForce are called.

Now, how to code the addDForce from your addForce ?
Basically, the addDForce term is the derivative of the addForce 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

[sharif]

Dear 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|)*(|dp|-r)*stiffness-v*(|dp|-r)*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 user-input parameter specifying the intensity of the field, v is the velocity of the particle, and damping is another user-input 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) for addDForce. How is it when we deal with vectors? Do gradians 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,
Sharif

[Hugo]

Dear Sharif,

No according to what I wrote, you should compute df(t)/dx * dx(t+dt) for AddDForce. 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

[Author]

Posts