Hello,

I have to implement a deformable body framework and this requires discretisization of a solid object using tetrahedral cells. From there on, each vertex of the tetrahedral mesh is used to compute various types of energies: length preserving (2 adjacent vertexes), area preserving (3 vertexes of a face), volume (4 vertexes) and so on (dihedral, angular, curvature)..

The theory says that in order to derive forces from the energy potentials, one can use the variational derivative of the energy to compute the force. This is F= -dE/dp
where p is vertex of the discrete mesh and F should be the force that's acting on that vertex.

For example, let's take the length preserving energy: E=1/2K ((|p_i-p_j|-L)/L)^2

This energy is defined globally, but it is clear that it can be decomposed into local potentials.

This means that only local potentials influence the force that's acting on a point, therefor one should find a way to derive the force from the local potentials.

Intuitively, these forces should be generated such that they minimize the total energy of the solid..
How can one derive a scalar field with respect to a 3 D vector (point,vertex) p and obtain a force?

If my intuition is right, one can define a sphere with radius h, centered at the vertex p that we're examining and see which point u on this sphere minimizes the energy..
The natural direction for the force should be p u.. and its magnitude somehow related to the difference of potentials.. perhaps given by the difference (E(p)-E(u))/h..
Please give me a suggestion..
Thank you!

Best,
Teodor

Working in three dimensions, x, y and z.

E = Σ(over i)Σ(over j) (1/2)K_ij (( | p_i − p_j | − L_ij ) / L_ij)^2
⇓
E_ij = (1/2)K_ij (( | p_i − p_j | − L_ij ) / L_ij)^2
⇓
E_ij = (1/2)(K_ij/L_ij^2)( | p_i − p_j | − L_ij )^2
⇓
F_ij = −∂E_ij/∂p_i_x = (K_ij/L_ij^2)( | p_i − p_j | − L_ij )(p_j_x − p_i_x)/( | p_i − p_j | )
⇓
F_ij = −∂E_ij/∂p_i = (K_ij/L_ij^2)( | p_i − p_j | − L_ij )(p_j − p_i)/( | p_i − p_j | )
⇓
F_i = Σ(over j) (K_ij/L_ij^2)( | p_i − p_j | − L_ij )(p_j − p_i)/( | p_i − p_j | )

(j and i have swapped places in one of the terms)

For the undeformed piece | p_i − p_j | = L_ij so F_ij and therefore F_i is zero. Therefore the undeformed state is a state that creates an extreme state of E (this is the variational principle). (If p_i = p_j the expression blows up because p_j − p_i
is a vector without length or direction so there is no preferred direction for the force. This may be a problem that has to be checked for and possibly dealt with if you do a computer implementation.)

Thank you for your reply..
there is still a problem which I didn't understand in the first place:
If E(p_i,p_j,p_k)=f(p_i,p_j,p_k) where f is a scalar valued function (e.g f= <((p_i-p_j)xp_k,p_j>), how is the variational derivative computed? (<u,v> means the dot product of u and v)
A hunch is that you substitute the components of P_i=(p_i_x,p_i_y,p_i_z) with (x,y,z) - variables, derive the function f with respect to x,y and z respectively to get 3 scalar values which form a vector then.. But is this right?

P.S.

I am not trying to have a certain example resolved, I am more interested in computing generally the variational derivative of a scalar valued function f:R^3->R_+ (positive real numbers).
To have a better picture of where this variational calculus issue comes in handy, I suggest looking over how the authors define energies to preserve several properties of a tetrahedron (in this article HTTP:// cg.informatik.uni-freiburg.de/publications/deformationCGI2004.pdf )

Best,
Teodor

See if you can prove these from x,y,z components of vectors a,b,c and scalar s,t?
D[s a, a_x] = s
D[s a + t b, a_x] = s
D[<a,a>, a_x] = 2 a_x
D[<a,b>, a_x] = D[<b,a>, a_x] = <(unit_x), b_x> = b_x
D[a x b, a_x] = (unit_x)xb
D[<a x b, c>, a_x] = <(unit_x)xb, c>
D[<a x b, c>, b_x] = <ax(unit_x), c>
D[<a x b, c>, c_x] = <axb, (unit_x)> = (axb)_x

Sorry I have to abuse notation -- I don't have time to typeset these well (or even double-check them).

http://www.libraryofmath.com/derivative-rules.html
http://en.wikipedia.org/wiki/Euclidean_vec...vector_function

These rules seem enough to derive most of the forces from the structure preserving potentials.. as soon as I have results, I will let you know whether they are practically accurate..
Again, thank you for posting these valuable hints here!

Hello,

I've implemented a simulation, but the results are unstable and lead to a very bad outcome (chaos).

Please, could you check if my computations are correct? It is just this line:

D[<axb,axc>; a_x]= <(unit_x)xb,axc>+<axb,(unit_x)xc>
i.e. is the variational derivative operator prone to the same rule for deriving a product?
(D( a . b )=D( a ).b + a.D( b ) from the definition of the Algebra of Derivatives ).

Best regards,

Teodor

This one is a little tricky, so I did it twice for you

D[<a x b,a x c>; a_x]
= D[<a,a><c,b>-<a,b><a,c>; a_x]
= 2<c,b>a_x - b_x<a,c> - c_x<a,b>
= 2 a_x <b,c> - b_x<a,c> - c_x<a,b>
= 2 a_x (b_x c_x + b_y c_y + b_z c_z) - b_x (a_x c_x + a_y c_y + a_z c_z) - c_x(a_x b_x + a_y b_y + a_z b_z)
= 2 a_x (b_y c_y + b_z c_z) - b_x (a_y c_y + a_z c_z) - c_x(a_y b_y + a_z b_z)


D[<a x b,a x c>; a_x]
= <D[a x b; a_x],(a x c)> + <D[a x c; a_x],(axb)>
= <(unit_x x b ), (a x c)> + <(unit_x x c), (a x b )>
= <(b_y(unit_z) - b_z(unit_y)), (a x c)> + <(c_y(unit_z) - c_z(unit_y)), (a x b )>
= b_y<(unit_z), (a x c)> - b_z<(unit_y), (a x c)> + c_y<(unit_z), (a x b )> - c_z<(unit_y), (a x b )>
= b_y(a x c)_z - b_z(a x c)_y + c_y(a x b )_z - c_z(a x b )_y
= b_y(a_x c_y - a_y c_x) + b_z(a_x c_z - a_z c_x) + c_y(a_x b_y - a_y b_x) + c_z(a_x b_z - a_z b_x)
= 2 a_x (b_y c_y + b_z c_z) - b_x (a_y c_y + a_z c_z) - c_x(a_y b_y + a_z b_z)

Hello..

I managed to prove this simple "lemma" last night.. using the same idea of the triple cross product and derivative algebra properties. Since I spent quite a lot of time trying to figure this out, I would like to wrap up the key results in the following lines.

The variational derivative operator D[W(x1,x2,..,xn);x_i] "maps" a scalar valued function W:R^3->R+, at a certain point, to a vector value. If the scalar field W is an energy-like function, the variational derivative D[E(x1,..,xn);x_i] will describe the force F_i acting on the point x_i so that the potential E decreases..

The D operator has the properties of general derivation operators.
1. linearity: D[aU+bV;x]=aD[U;x]+bD[V;x] where a,b are scalars, U, V potential functions, x a 3D point.

2. product rule D[U*V;x]=D[U;v]*V+U*D[V;x]

3. composition rule: D[f(|U|);x]=f'(|U|)*D[U;x], where f:R->R, f' is the usual derivative, |U| is the norm of U.. Instead of |U| one can write sqrt(<U,U>). More generally, if g:R^3->R+ maps U onto a positive real scalar, g being a bilinear operator (e.g. the usual dot product g(U)=<U,V>, where V is fixed in this case), we can write D[f(g(U));x]=f'(g(U))*D[g(U);x]. computing D[g(U);x] is what computer graphics problems usually require.

As pointed out by rpenner, there are several ways one can compute D[g(U);x].
Assume a,b,c are R^3 vectors and s,t are scalars.
The result that we get when applying the variational derivative operator is a vector, but we can compute it only component-wise (each component at a time).

If R=D[g(U);w] (instead of the vector x, we will use w not to confuse it with the x,y and z component of a 3D point) then R_x:=D[g(U);w_x] and the same goes for R_y and R_z (this is a slight abuse of notation). R=[R_x,R_y,R_z] is now the desired result vector (e.g. the force if -g(U) is the energy).

The fundamental laws needed to derive a linear operator (g in our case) are easily provable if we consider the previously stated properties (1, 2, 3). Some notable results:
by convention D[a;a_x]=ex and thus D[sa+tb;a_x]=s ex (rule number 1.). This is just a way of quickly replacing a with ex in the bilinear form g.

-- D[<a,b>;a_x]=<ex,b> , where ex=(1,0,0) (the x axis versor)
-- D[<a,a>;a_x]=<ex,a>+<a,ex>=2<ex,a> (from the 2nd property of derivatives 2.)
-- D[<axb,c>;a_x]=<exxb,c>. Based on the previous two results, we can compute several other combinations.
-- D[<axb,axc>;a_x]=D[<exXb,axc>;a_x]+D[<axb,exXc] (2nd property again)

I suggest the reader to look through computer graphics/physically based animation papers to see how various energies are defined. To compute the corresponding forces, one needs to have a list of properties for the dot product <;>, the cross product "x" and also fully understand the simple rules a variational derivative operator obeys. I have to say that these observations are not rigorous from a mathematicians point of view, but they "work" for an applied computer scientist. I tested the energy potentials vs. the variational derivative way of computing forces and all of them seem to behave fine (e.g. the tend to preserve the exact properties of a tetrahedron when a vertex is displaced "out of the blue").

All the credits go to rpenner for providing a wonderful "minitutorial".