TAPsVolPresTriModel.cpp

Go to the documentation of this file.

/******************************************************************************
TAPsVolPresTriModel.cpp

A Triangulated Model that is deformable with volume preservation.

SUKITTI PUNAK   (11/02/2004)
UPDATE          (11/11/2004)
********************************************************************************
+====================+
| DESCRIPTION / NOTE |
+====================+

The model is a triangulated mesh of V (number of vertices) masses, which each 
mass linked to its neighbors by massless springs of natural length greater than 
zeroes.

Since the model is triangulated mesh, the linkage between neighbors comprised 
of two kinds of linking springs:
  -- Structural springs linking vertices of a face f along its edges
     (i.e. metric deformation)
     e.g.    2-----1      2
              \   /      /|
               \ /      / |
                0      0--1
  -- Flexion springs linking diagonally between two faces
     (i.e. bending deformation)
     e.g.     2-----1    this example shows the link of vertex 2 and 3
               \ + / \
                \ / + \
                 0-----3
  (-- No Shear springs (i.e. shearing deformation))

                    (2)-----------(2)-----------(2)
                    / \           /|\           / \
                   /   \         / | \         /   \
                  /     \       /  |  \       /     \
                 /       \     /   |   \     /       \
                /         \   /    |    \   /         \
               /           \ /     |     \ /           \
             (2)-----------(1)-----|-----(1)-----------(2)
             / \  .        / \     |     / \        .  / \
            /   \   .     /   \    |    /   \     .   /   \
           /     \    .  /     \   |   /     \  .    /     \
          /       \     / .     \  |  /     . \     /       \
         /         \   /    .    \ | /    .    \   /         \
        /           \ /        .  \|/  .        \ /           \
      (2)-----------(1)-----------(0)-----------(1)----------(2)
        \           / \        .  /|\  .        / \           /
         \         /   \    .    / | \   .     /   \         /
          \       /     \ .     /  |  \    .  /     \       /
           \     /     . \     /   |   \     / .     \     /
            \   /   .     \   /    |    \   /    .    \   /
             \ / .         \ /     |     \ /        .  \ /
             (2)-----------(1)-----|-----(1)-----------(2)
               \           / \     |     / \           /
                \         /   \    |    /   \         /
                 \       /     \   |   /     \       /
                  \     /       \  |  /       \     /
                   \   /         \ | /         \   /
                    \ /           \|/           \ /
                    (2)------------(2)----------(2)
  
  Data Structure:
  ===============
  Need edge list
  Each edge incidents to two faces, find cross list



    In case of Quadrateral Mesh
  (0,0)---------(0,1)---------(0,2)
    |  \       /  |  \       /  |
    |    \   /    |    \   /    |
    |      +      |      +      |
    |    /   \    |    /   \    |
    |  /       \  |  /       \  |
  (1,0)---------(1,1)---------(1,2)
    |  \       /  |  \       /  |
    |    \   /    |    \   /    |
    |      +      |      +      |
    |    /   \    |    /   \    |
    |  /       \  |  /       \  |
  (2,0)---------(2,1)---------(2,2)

******************************************************************************/
#include "TAPsVolPresTriModel.hpp"
// Using Inclusion Model (i.e. definitions are included in declarations)
//                       (this name.cpp is included in name.hpp)
// Each friend is defined directly inside its declaration.

BEGIN_NAMESPACE_TAPs__OpenGL
//=============================================================================
//-----------------------------------------------------------------------------
// default constructor
template <typename T>
VolPresTriModel<T>::VolPresTriModel ()
    : XTrigonalModel<T>(), m_tVolume(), m_vtCentroid(),
      m_iNoEdges( 0 ), m_prEdge( NULL ),
      // Spring
      m_pprParticleVertex( NULL ),
      m_pprSpringVertex( NULL )
{
    #ifdef  TAPs_ENABLE_DEBUG
    std::cout << "VolPresTriModel<" << typeid(T).name() << "> Constructor\n";
    #endif//TAPs_ENABLE_DEBUG
}
//-----------------------------------------------------------------------------
// destructor
template <typename T>
VolPresTriModel<T>::~VolPresTriModel ()
{
    delete [] m_prEdge;
    m_prEdge = NULL;

    #ifdef  TAPs_ENABLE_DEBUG
    std::cout << "VolPresTriModel<" << typeid(T).name() << "> Destructor\n";
    #endif//TAPs_ENABLE_DEBUG

    // Delete All ParticleVertex
    for ( int i = 0; i < m_iNoVertices; ++i ) {
        delete m_pprParticleVertex[i];
    }
    delete [] m_pprParticleVertex;
    // Delete All SpringVertex
    for ( int i = 0; i < m_iNoEdges; ++i ) {
        delete m_pprSpringVertex[i];
    }
    delete [] m_pprSpringVertex;
}
//-----------------------------------------------------------------------------
// Initialize
template <typename T>
void VolPresTriModel<T>::Initialize ()
{
    // Model Initialization
    DetermineAndSortRings();
    CalAndSetNormals();
    // DetermineAndSortRings();
    ApplyMaterial();
    DetermineEdgeList();
}
//-----------------------------------------------------------------------------
// Find Edge List
template <typename T>
void VolPresTriModel<T>::DetermineEdgeList ()
{
    std::vector<int> *edgeRecord;
    edgeRecord = new std::vector<int>[ m_iNoVertices ];
    int iOrig, iDest;   // original vertex# and destination vertex# of edge
    int iV;     // temp vertex number
    int numberOfEdges = 0;

    // Get edges in a temp data
    for ( int f = 0; f < m_iNoFaces; ++f ) {
        iV = m_prFace[f].GetVertexNo( 0 );
        for ( int v = 1; v < m_prFace[f].GetNoVertices(); ++v ) {
            iOrig = iV;
            iDest = m_prFace[f].GetVertexNo( v );
            
            // Make the original vertex always smaller than the destination vertex
            if ( iOrig > iDest ) {
                int tmp = iOrig;
                iOrig = iDest;
                iDest = tmp;
            }
            iV = iDest;     // keep the current vertex

            // Record the edge
            bool notStored = true;
            for ( int e = 0; e < static_cast<int>( edgeRecord[iOrig].size() ); ++e ) {
                if ( edgeRecord[iOrig][e] == iDest ) {
                    notStored = false;
                    break;
                }
            }
            if ( notStored ) {
                edgeRecord[iOrig].push_back( iDest );
                ++numberOfEdges;
                //std::cout << "Edge: from " << iOrig << " to " << iDest << std::endl;
            }
        }
    }

    // Set edges
    int edgeNo = 0;
    m_iNoEdges = numberOfEdges;
    m_prEdge = new Edge[numberOfEdges];
    for ( int v = 0; v < m_iNoVertices; ++v ) {
        for ( int e = 0; e < static_cast<int>( edgeRecord[v].size() ); ++e ) {
            m_prEdge[edgeNo++].SetOrigAndDest( v, edgeRecord[v][e] );
            //std::cout << "Edge: from " << v << " to " << edgeRecord[v][e] << std::endl;
        }
    }

    CreateParticleVerticesFromVertices();
    CreateSpringVerticesFromEdges();
    CalVolume2();

    /*
    for ( int i = 0; i < m_iNoEdges; ++i ) {
        std::cout << "Spring#" << i << ":" << *m_pprSpring[i] << "\n";
    }
    */
}
//-----------------------------------------------------------------------------
// Create Particles From Vertices
template <typename T>
void VolPresTriModel<T>::CreateParticleVerticesFromVertices ()
{
    m_pprParticleVertex = new ParticleVertex<T>*[m_iNoVertices];
    for ( int i = 0; i < m_iNoVertices; ++i ) {
        m_pprParticleVertex[i] = new ParticleVertex<T> ( 
            2.0,                // Mass
            m_prXVertex[i],     // Vertex (including Position and Normal Vectors)
            Vector3<T>(),       // Velocity
            Vector3<T>(),       // Acceleration
            Vector3<T>()        // Force
        );
    }
}
//-----------------------------------------------------------------------------
// Create Springs From Edges
template <typename T>
void VolPresTriModel<T>::CreateSpringVerticesFromEdges ()
{
    m_pprSpringVertex = new SpringVertex<T>*[m_iNoEdges];
    for ( int i = 0; i < m_iNoEdges; ++i ) {
        m_pprSpringVertex[i] = new SpringVertex<T> ( 
            *m_pprParticleVertex[m_prEdge[i].GetOriginal()],    // Vertex#1
            *m_pprParticleVertex[m_prEdge[i].GetDestination()], // Vertex#2
            10.0,   // K (spring constant)
            0.25    // D (damping constant)
        );
    }
}
//-----------------------------------------------------------------------------
// Calculate Volume
template <typename T>
T VolPresTriModel<T>::CalVolume ()
{
    Vector3<T> va, vb, vc;
    Vector3<T> vd;
    Vector3<T> vmom;
    T tv;
    T tVolume = 0.0;
    for ( int i = 0; i < m_iNoFaces; ++i ) {
        va.SetXYZ ( 
            m_prXVertex[m_prFace[i].GetVertexNo(0)][0],
            m_prXVertex[m_prFace[i].GetVertexNo(0)][1],
            m_prXVertex[m_prFace[i].GetVertexNo(0)][2]
        );
        vb.SetXYZ ( 
            m_prXVertex[m_prFace[i].GetVertexNo(1)][0],
            m_prXVertex[m_prFace[i].GetVertexNo(1)][1],
            m_prXVertex[m_prFace[i].GetVertexNo(1)][2]
        );
        vc.SetXYZ ( 
            m_prXVertex[m_prFace[i].GetVertexNo(2)][0],
            m_prXVertex[m_prFace[i].GetVertexNo(2)][1],
            m_prXVertex[m_prFace[i].GetVertexNo(2)][2]
        );

        tVolume += tv = ( va * ( vb ^ vc ) ) / 6.0;
        vd = ( va + vb + vc ) / 4.0;
        vmom += vd * tv;
    }
    //m_tVolume = tVolume;
    m_vtCentroid = vmom / tVolume;
    return tVolume;
}
//-----------------------------------------------------------------------------
// Calculate Volume Too
template <typename T>
T VolPresTriModel<T>::CalVolume2 ()
{
    T wgt[4] = { 0, 0, 0, 0 };
    T v[3][3];
    T dv[3][3];

    dv[0][0] = dv[0][1] = dv[0][2] = 0;
    dv[1][0] = dv[1][1] = dv[1][2] = 0;
    dv[2][0] = dv[2][1] = dv[2][2] = 0;
    //*/

    for ( int i = 0; i < m_iNoFaces; ++i ) {
        v[0][0] = m_prXVertex[ m_prFace[i].GetVertexNo(0) ][0];
        v[0][1] = m_prXVertex[ m_prFace[i].GetVertexNo(0) ][1];
        v[0][2] = m_prXVertex[ m_prFace[i].GetVertexNo(0) ][2];
        v[1][0] = m_prXVertex[ m_prFace[i].GetVertexNo(1) ][0];
        v[1][1] = m_prXVertex[ m_prFace[i].GetVertexNo(1) ][1];
        v[1][2] = m_prXVertex[ m_prFace[i].GetVertexNo(1) ][2];
        v[2][0] = m_prXVertex[ m_prFace[i].GetVertexNo(2) ][0];
        v[2][1] = m_prXVertex[ m_prFace[i].GetVertexNo(2) ][1];
        v[2][2] = m_prXVertex[ m_prFace[i].GetVertexNo(2) ][2];

        d3mom0 ( v, dv, wgt );  // call subroutine
    }

    m_tVolume = wgt[0];
    return m_tVolume;
}
//-----------------------------------------------------------------------------
// Preserve Volume
template <typename T>
void VolPresTriModel<T>::PreserveVolume ( T wgt[4] )
{
    wgt[0] = wgt[1] = wgt[2] = wgt[3] = 0;
    T v[3][3];
    T dv[3][3];

    /*
    dv[0][0] = dv[0][1] = dv[0][2] = 0;
    dv[1][0] = dv[1][1] = dv[1][2] = 0;
    dv[2][0] = dv[2][1] = dv[2][2] = 0;
    */

    for ( int i = 0; i < m_iNoFaces; ++i ) {
        v[0][0] = m_prXVertex[ m_prFace[i].GetVertexNo(0) ][0];
        v[0][1] = m_prXVertex[ m_prFace[i].GetVertexNo(0) ][1];
        v[0][2] = m_prXVertex[ m_prFace[i].GetVertexNo(0) ][2];
        v[1][0] = m_prXVertex[ m_prFace[i].GetVertexNo(1) ][0];
        v[1][1] = m_prXVertex[ m_prFace[i].GetVertexNo(1) ][1];
        v[1][2] = m_prXVertex[ m_prFace[i].GetVertexNo(1) ][2];
        v[2][0] = m_prXVertex[ m_prFace[i].GetVertexNo(2) ][0];
        v[2][1] = m_prXVertex[ m_prFace[i].GetVertexNo(2) ][1];
        v[2][2] = m_prXVertex[ m_prFace[i].GetVertexNo(2) ][2];

        dv[0][0] = m_prXVertex[ m_prFace[i].GetVertexNo(0) ].GetNormal().GetX();
        dv[0][1] = m_prXVertex[ m_prFace[i].GetVertexNo(0) ].GetNormal().GetY();
        dv[0][2] = m_prXVertex[ m_prFace[i].GetVertexNo(0) ].GetNormal().GetZ();
        dv[1][0] = m_prXVertex[ m_prFace[i].GetVertexNo(1) ].GetNormal().GetX();
        dv[1][1] = m_prXVertex[ m_prFace[i].GetVertexNo(1) ].GetNormal().GetY();
        dv[1][2] = m_prXVertex[ m_prFace[i].GetVertexNo(1) ].GetNormal().GetZ();
        dv[2][0] = m_prXVertex[ m_prFace[i].GetVertexNo(2) ].GetNormal().GetX();
        dv[2][1] = m_prXVertex[ m_prFace[i].GetVertexNo(2) ].GetNormal().GetY();
        dv[2][2] = m_prXVertex[ m_prFace[i].GetVertexNo(2) ].GetNormal().GetZ();
        //*/
        
        d3mom0 ( v, dv, wgt );  // call subroutine
    }

    T diffVolume = wgt[0] - m_tVolume;
    if ( diffVolume != 0 ) {
    //if ( fabs( diffVolume ) > Math<T>::ROUND_ERROR ) {
        T sigma = SP_Maths::rootOfCubicPolynomial( wgt[3], wgt[2], wgt[1], diffVolume, 1 );

        #ifdef  TAPs_ENABLE_DEBUG
        std::cout << "\nsigma: " << sigma << "\n";
        #endif//TAPs_ENABLE_DEBUG

        for ( int i = 0; i < m_iNoVertices; ++i ) {
            m_prXVertex[i].SetPosition( m_prXVertex[i].GetPosition() + sigma * m_prXVertex[i].GetNormal() );
            //m_prXVertex[i].SetPosition( m_prXVertex[i].GetPosition() * sigma );
        }
    }
    
    // Recalculate Normals
    CalAndSetFaceNormalsNotNormalized();
    CalAndSetVertexNormals();
    NormalizeFaceNormals();

    // For Checking
    CalVolume2();
}
//-----------------------------------------------------------------------------
/* IN : coefficients  v  of degree  1  patch
 *   and coefficients dv  of the pertubation
 * OUT: weights   wgt  of moment 0: 
 *     vol = wgt[0]+wgt[1]*s+wgt[2]*s^2+wgt[3]*s^3 
 * generated by Jorg Peters and Maple */
//#define PTS 3 
//#define XYZ 3
//double d3mom0 (double v[PTS][XYZ], double
//               dv[PTS][XYZ], double wgt[4] )
template <typename T>
inline void VolPresTriModel<T>::d3mom0 ( T v[][3], T dv[][3], T wgt[4] )
//void VolPresTriModel<T>::d3mom0 ()
{
    T t[184];
    //-------------------------------------------------------------------------
    t[1]  = v[1][2];
    t[2]  = v[2][0];
    t[3]  = t[1]*t[2];
    t[4]  = v[1][1];
    t[6]  = v[0][2];
    t[7]  = t[6]*t[2];
    t[9]  = v[1][0];
    t[10] = t[6]*t[9];
    t[11] = v[2][1];
    t[13] = v[0][1];
    t[15] = v[0][0];
    t[16] = t[6]*t[15];
    t[20] = t[1]*t[9];
    t[23] = t[1]*t[15];
    t[27] = v[2][2];
    t[28] = t[27]*t[2];
    t[30] = t[27]*t[9];
    t[33] = t[27]*t[15];
    t[37] = -t[3]*t[4]-t[7]*t[4]+t[10]*t[11]-t[10]*t[13]-t[16]*t[11]+t[7]*t[13]+t[16]*t[4]+t[20]*t[11]-t[20]*t[13]-
        t[23]*t[11]+t[3]*t[13]+t[23]*t[4]-t[28]*t[4]+t[30]*t[11]-t[30]*t[13]-t[33]*t[11]+t[28]*t[13]+t[33]*t[4];
    wgt[0] += t[37]/6.0;
    //-------------------------------------------------------------------------
    t[38] = dv[1][0];
    t[40] = dv[2][0];
    t[43] = dv[1][1];
    t[45] = dv[2][1];
    t[47] = dv[0][1];
    t[51] = dv[0][0];
    t[56] = t[38]*t[11]-t[40]*t[4]+t[40]*t[13]+t[15]*t[43]+t[9]*t[45]-t[9]*t[47]-t[38]*t[13]-t[15]*t[45]-t[51]*
        t[11]+t[2]*t[47]-t[2]*t[43]+t[51]*t[4];
    t[58] = dv[0][2];
    t[65] = -t[2]*t[4]+t[9]*t[11]-t[9]*t[13]-t[15]*t[11]+t[2]*t[13]+t[15]*t[4];
    t[68] = dv[1][2];
    t[71] = dv[2][2];
    wgt[1] += t[6]*t[56]/6.0 + t[58]*t[65]/6.0 + t[1]*t[56]/6.0 + t[68]*t[65]/6.0 + t[27]*t[56]/6.0 + 
        t[71]*t[65]/6.0;
    //-------------------------------------------------------------------------
    t[74] = t[71]*t[15];
    t[76] = t[68]*t[15];
    t[78] = t[58]*t[40];
    t[81] = t[71]*t[9];
    t[83] = t[68]*t[38];
    t[87] = t[71]*t[2];
    t[89] = t[58]*t[15];
    t[91] = t[68]*t[2];
    t[93] = t[58]*t[2];
    t[96] = t[74]*t[43]-t[76]*t[45]+t[78]*t[13]+t[76]*t[43]-t[81]*t[47]-t[83]*t[13]-t[78]*t[4]+t[83]*t[11]+t[87]*
        t[47]-t[89]*t[45]-t[91]*t[43]+t[93]*t[47]+t[81]*t[45];
    t[98] = t[58]*t[38];
    t[100] = t[71]*t[51];
    t[102] = t[6]*t[38];
    t[104] = t[27]*t[51];
    t[106] = t[27]*t[38];
    t[108] = t[27]*t[40];
    t[111] = t[6]*t[51];
    t[113] = t[1]*t[38];
    t[115] = t[71]*t[40];
    t[121] = t[98]*t[11]-t[100]*t[11]-t[102]*t[47]-t[104]*t[45]+t[106]*t[45]-t[108]*t[43]+t[108]*t[47]+t[111]
        *t[43]-t[113]*t[47]-t[115]*t[4]+t[91]*t[47]-t[87]*t[43]+t[89]*t[43]+t[113]*t[45];
    t[124] = t[68]*t[40];
    t[126] = t[68]*t[51];
    t[130] = t[1]*t[51];
    t[133] = t[6]*t[40];
    t[136] = t[71]*t[38];
    t[142] = -t[124]*t[4]+t[126]*t[4]+t[100]*t[4]-t[74]*t[45]-t[130]*t[45]+t[102]*t[45]+t[133]*t[47]-t[93]*
        t[43]-t[136]*t[13]+t[130]*t[43]-t[98]*t[13]-t[111]*t[45]-t[133]*t[43];
    t[144] = t[58]*t[9];
    t[147] = t[58]*t[51];
    t[150] = t[1]*t[40];
    t[154] = t[68]*t[9];
    t[162] = t[144]*t[45]-t[144]*t[47]-t[147]*t[11]+t[147]*t[4]+t[150]*t[47]-t[150]*t[43]+t[124]*t[13]+t[154]
        *t[45]-t[154]*t[47]-t[126]*t[11]+t[104]*t[43]-t[106]*t[47]+t[136]*t[11]+t[115]*t[13];
    wgt[2] += t[162]/6.0 + t[142]/6.0 + t[121]/6.0 + t[96]/6.0;
    //-------------------------------------------------------------------------
    t[183] = t[78]*t[47]+t[147]*t[43]+t[98]*t[45]-t[98]*t[47]-t[147]*t[45]-t[78]*t[43]+t[124]*t[47]+t[126]*
        t[43]+t[83]*t[45]-t[83]*t[47]-t[126]*t[45]-t[124]*t[43]+t[115]*t[47]+t[100]*t[43]+t[136]*t[45]-t[136]*t[47]-t[100]*
    t[45]-t[115]*t[43];
    wgt[3] += t[183]/6.0;
    //-------------------------------------------------------------------------
}
//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
// StepSimulation
template <typename T>
void VolPresTriModel<T>::StepSimulation( T dt ) {
    T divider = 10.0;
    
    //CalForces( Vector3<T>( (rand()%1000)/divider, (rand()%10)/divider, (rand()%10)/divider ) );
    CalForces( Vector3<T>() );
    
    //CalForces();
    // Euler's Method for Integration for ODE Solver
    for ( int n = 0; n < m_iNoVertices; ++n ) {
        m_pprParticleVertex[n]->SetAcceleration( m_pprParticleVertex[n]->GetForce() / m_pprParticleVertex[n]->GetMass() );
        m_pprParticleVertex[n]->SetVelocity( m_pprParticleVertex[n]->GetVelocity() 
                                            + ( m_pprParticleVertex[n]->GetAcceleration() * dt ) );
        m_pprParticleVertex[n]->SetPosition( m_pprParticleVertex[n]->GetPosition() 
                                            + ( m_pprParticleVertex[n]->GetVelocity() * dt ) );
    }
}
//-----------------------------------------------------------------------------
// Calculate All Forces Acting on Particles
template <typename T>
void VolPresTriModel<T>::CalForces( Vector3<T> vWind, T g )
{
    //CalParticleNormals( 0 );  // calculate particle normals use for finding force due to wind
    static int simCounter = 0;

    T zero = 0, y = 0;
    for ( int n = 0; n < m_iNoVertices; ++n ) {
        if ( m_pprParticleVertex[n]->GetFixStatus() ) {     // If the particle is fixed, skip it!
            m_pprParticleVertex[n]->SetForce( 0, 0, 0 );
            continue;
        }
        // Clear and Calculate Force due to Acceleration of Gravity --------
        //     Fij = mass * g; where g is a vector
        y = m_pprParticleVertex[n]->GetMass() * (-g);

        /*
        m_pprParticleVertex[n]->SetForce( 0, y, 0 );
        */
        if ( simCounter++ < 10 && ( n == 0 || n == 4 ) ) {
            m_pprParticleVertex[n]->SetPosition( m_pprParticleVertex[n]->GetPosition() * 0.9999999999 );
        }
        //*/

        // Add Force due to Viscous Damping --------------------------------
        //     Fij = -Kvd * Vij
        m_pprParticleVertex[n]->SetForce( m_pprParticleVertex[n]->GetForce() 
                + 0.1/*m_tKvd (viscous damping)*/ * m_pprParticleVertex[n]->GetVelocity() );

        // Add Force due to viscous fluid moving (wind) --------------------
        //     Fij = Kvfm * [Nij * (vWind - Vij)] * Nij
        //         where Nij is the unit normal on the surface at particle Pij
        m_pprParticleVertex[n]->SetForce( m_pprParticleVertex[n]->GetForce()
                + 0.1/*m_tKvfm (viscous fluid moving damping)*/ 
                    * ( m_pprParticleVertex[n]->GetNormal() * ( vWind - m_pprParticleVertex[n]->GetVelocity() ) )
                    * m_pprParticleVertex[n]->GetNormal() );

        // Dump
        //cout << "Mass: " << m_pcParticle[m][n]->m_tMass << "    Force: " << m_pcParticle[m][n]->m_vtForce << "\n";
    }
    
    AddSpringForces();  // calculate and add spring forces to particle forces
}
//-----------------------------------------------------------------------------
// Calculate and Add the Force of a Spring
template <typename T>
inline void VolPresTriModel<T>::AddSpringForces()
{
    for ( int i = 0; i < m_iNoEdges; ++i ) {

        //------------------------------------------------------------
        // Local Variables
        // spring current length vector
        Vector3<T> cL = ( m_pprSpringVertex[i]->GetParticleVertexOne().GetPosition()
                        - m_pprSpringVertex[i]->GetParticleVertexTwo().GetPosition() );
        // spring current length scalar
        T scl = cL.Length();
        // velocity difference of the two particles linked by the spring
        Vector3<T> V = ( m_pprSpringVertex[i]->GetParticleVertexOne().GetVelocity()
                       - m_pprSpringVertex[i]->GetParticleVertexTwo().GetVelocity() );

        // F1 = -{Ks(l - r) + Kd[(V1-V2)*L]/l}*L/l
        // F2 = -F1
        //   where Ks = spring stiffness, Kd = spring damping,
        //     V1 and V2 are velocity of particles linked by the spring
        //     r = spring rest length
        //     L = vector of difference of positions of particle#1 and #2
        //     l = magnitude of L
        Vector3<T> f  = -( 
                  ( m_pprSpringVertex[i]->GetConstantK() * ( scl - m_pprSpringVertex[i]->GetRestLengthL() ) ) 
                + ( m_pprSpringVertex[i]->GetDampingD() * ((V*cL) / scl) ) ) 
                * (cL / scl);
        //*/

        //std::cout << "Spring Force: " << f << std::endl;
        /*
        if ( !m_pprSpringVertex[i]->GetParticleVertexOne().GetFixStatus() )
            m_pprSpringVertex[i]->GetParticleVertexOne().SetForce( f );
        if ( !m_pprSpringVertex[i]->GetParticleVertexTwo().GetFixStatus() )
            m_pprSpringVertex[i]->GetParticleVertexTwo().SetForce( -f );
        */
        if ( !m_pprSpringVertex[i]->GetParticleVertexOne().GetFixStatus() )
            m_pprSpringVertex[i]->GetParticleVertexOne().SetForce( 
                      m_pprSpringVertex[i]->GetParticleVertexOne().GetForce()
                    + m_pprSpringVertex[i]->GetForce() 
                    + f
            );
        if ( !m_pprSpringVertex[i]->GetParticleVertexTwo().GetFixStatus() )
            m_pprSpringVertex[i]->GetParticleVertexTwo().SetForce( 
                      m_pprSpringVertex[i]->GetParticleVertexTwo().GetForce()
                    - m_pprSpringVertex[i]->GetForce() 
                    - f 
            );
        //*/
    }
}
//-----------------------------------------------------------------------------
//=============================================================================
END_NAMESPACE_TAPs__OpenGL
//-----------------------------------------------------------------------------
//345678901234567890123456789012345678901234567890123456789012345678901234567890
//--+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8