5 August 2015 at 12 h 02 min #3444sharifParticipant
- PhD Student
Currently, the amount of force in a
VaccumShereForceFieldis 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(…).
Sharif5 August 2015 at 16 h 20 min #3448HugoKeymaster
- SOFA Consortium
I will try to explain it in a short way. The use of
addDForcedepends 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
- 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
addForcefunction and the
addDForcecomputes the term dt*K, used for both right- and left-hand side.
- in case of an explicit scheme, only the
- in case of an implicit scheme, both
Now, how to code the
addDForceterm is the derivative of the
addForceterm, 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.
Hugo5 August 2015 at 17 h 30 min #3459sharifParticipant
- PhD Student
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
fis the force vector,
dpthe position vector of a particle relative to the center of the field,
ris the radius of the sphere where the forcefield applies,
stiffnessis a user-input parameter specifying the intensity of the field,
vis the velocity of the particle, and
dampingis 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
gradiansapply in this case?
As another question, what is
mparams.kFactorand 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!
Sharif20 August 2015 at 13 h 01 min #3485HugoKeymaster
- SOFA Consortium
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.
- You must be logged in to reply to this topic.