TAPsElasticRod_CompByCUDA.cpp

Go to the documentation of this file.

/******************************************************************************
TAPsElasticRod_CompByCUDA.cpp
******************************************************************************/
/******************************************************************************
SUKITTI PUNAK   (09/02/2009)
******************************************************************************/
#include "TAPsElasticRod_CompByCUDA.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
//=============================================================================
//-----------------------------------------------------------------------------

// Constructor
template <typename T>
ElasticRod_CompByCUDA<T>::ElasticRod_CompByCUDA () : 
    m_pIdxCurr( 0 ),
    m_pIdxPrev( 0 ),
    m_pIdxNext( 0 ),
    m_pParameters( NULL ),
    m_pNodeList( NULL ),
    m_cudaPosList( NULL ),
    m_cudaPrevPosList( NULL ),
    m_cudaOriList( NULL ),
    m_cudaPrevOriList( NULL ),
    m_cudaIntForceList( NULL ),
    m_cudaExtForceList( NULL ),
    m_cudaID( 0 ),
    m_bUsePLHM( false )
{}


// Destructor
template <typename T>
ElasticRod_CompByCUDA<T>::~ElasticRod_CompByCUDA ()
{
    DeleteCUDA();
}


// Create CUDA
template <typename T>
void ElasticRod_CompByCUDA<T>::CreateCUDA (
    int *                       pIdxCurr,       
    int *                       pIdxPrev,       
    int *                       pIdxNext,       
    ElasticRodParameters<T> *   pParameters,    
    SetOfElasticRodNodes *      pNodeList,      
    bool                        bUsePLHM        
)
{
    m_pIdxCurr    = pIdxCurr;
    m_pIdxPrev    = pIdxPrev;
    m_pIdxNext    = pIdxNext;
    m_pParameters = pParameters;
    m_pNodeList   = pNodeList;

    // Setting utilizing CUDA's page-locked host memory
    m_bUsePLHM = bUsePLHM;

//#if CUDART_VERSION >= 2020
    if ( m_bUsePLHM ) {
        int device;
        TAPs::CUDA::GetDevice( device );
        m_bUsePLHM = TAPs::CUDA::CanDeviceMapHostMemory( device );
        std::cout << "m_bUsePLHM: " << m_bUsePLHM << " printed from by File: " << __FILE__ << "; Line: " << __LINE__ << "\n";
    }
//#else
//  m_bUsePLHM = false;
//#endif

    try {
        if ( m_bUsePLHM )   InitPLHM();
        else                Init();
    }
    catch ( char * e ) {
        throw e;
    }
}


// Delete CUDA
template <typename T>
void ElasticRod_CompByCUDA<T>::DeleteCUDA ()
{
    if ( m_bUsePLHM )   CleanPLHM();
    else                Clean();
}


// Advance simulation by time step
template <typename T>
void ElasticRod_CompByCUDA<T>::AdvanceSimulation ()
{
    if ( m_bUsePLHM )   AdvSimPLHM();
    else                AdvSim();
}


// Advance simulation by time step
template <typename T>
void ElasticRod_CompByCUDA<T>::AdvSim ()
{
    // Copy the (current) vertex positions of the object to memory for transfering to CUDA.
    CopyPosToMem();
    CopyIntForceToMem();
    CopyExtForceToMem();
    CopyOriToMem();
    #ifdef  TAPs_ADVANCED_SIMULATION
        //CopySimFlagsToMem();
        //CopyPosConstraintToMem();
    #endif//TAPs_ADVANCED_SIMULATION

    // Call the wrapper function
    TAPs::CUDA::Global__ModelElasticRod_AdvSim( 
        m_cudaID,
        m_pParameters->NumOfNodes,  
        64,                         
        0,                          
        m_pParameters->TimeStep,    
        1,                          
        m_pParameters->Radius,      
        m_pParameters->LinkLength,  
        m_pParameters->LinkLength,  
        m_pParameters->MassOfPoint, 
        m_pParameters->MaterialDensity, 
        //m_pParameters->Kt,        //!< kinetic translational constant
        //m_pParameters->Kr[0],     //!< kinetic rotational constant -- x
        //m_pParameters->Kr[1],     //!< kinetic rotational constant -- y
        //m_pParameters->Kr[2],     //!< kinetic rotational constant -- z
        //m_pParameters->Dt,        //!< translational constant
        //m_pParameters->Dr[0],     //!< rotational dissipation constant -- x
        //m_pParameters->Dr[1],     //!< rotational dissipation constant -- y
        //m_pParameters->Dr[2],     //!< rotational dissipation constant -- z
        m_pParameters->Kconstraint_3rdDirAlignTangent,  
        m_pParameters->Kvdamping,   
        m_pParameters->Ps,          
        m_pParameters->Pb[0],       
        m_pParameters->Pb[1],       
        m_pParameters->Pb[2],       
        m_cudaPosList,          
        m_cudaPrevPosList,      
        m_cudaOriList,          
        m_cudaPrevOriList,      
        m_cudaIntForceList,     
        m_cudaExtForceList      
    );


    //*
    // The new vertex positions are returned by CUDA to memory for previous vertex positions.
    // So the new vertex positions must be copied to the vertex positions of the object.
    CopyMemToPrevPos();
    CopyMemToPrevOri();
    // Then the previous and current vertex pointers (that point to memory for transfering to CUDA) are swapped.
    // So that now the previous becomes current and the current becomes previous.
    SwapBuffers();
    //*/

    /*
    // DEBUG -- PRINTOUT
    {
        std::cout << "CUDA\n";
        SetOfElasticRodNodes::iterator node = m_pNodeList->begin();
        int i = 0;
        int p = 0;
        int P = *m_pIdxNext;
        while ( node != m_pNodeList->end() ) {
            std::cout << "C" << i << ": " << m_cudaPrevPosList[p+0] << " " << m_cudaPrevPosList[p+1] << " " << m_cudaPrevPosList[p+2] << "\n";
            p += 4;
            ++i;
            ++node;
        }
    }
    {
        SetOfElasticRodNodes::iterator node = m_pNodeList->begin();
        int i = 0;
        int p = 0;
        int P = *m_pIdxNext;
        while ( node != m_pNodeList->end() ) {
            std::cout << "O" << i << ": " << m_cudaPrevOriList[p+0] << " " << m_cudaPrevOriList[p+1] << " " << m_cudaPrevOriList[p+2] << " " << m_cudaPrevOriList[p+3] << "\n";
            p += 4;
            ++i;
            ++node;
        }
    }
    //*/


    #ifdef  TAPs_ADVANCED_SIMULATION
    // The new strand's simulation flags list is returned by CUDA.
    // It must be copied to the simulation flags list of the strand.
        //CUDA_CopyMemToSimFlags();
        //CUDA_CopyMemToPosConstraint();

        // NEED TO CHANGE THE VALUES STORED IN m_pAdvSimData

    #endif//TAPs_ADVANCED_SIMULATION
}


// Advance simulation by time step
template <typename T>
void ElasticRod_CompByCUDA<T>::AdvSimPLHM ()
{
    // Copy the (current) vertex positions of the object to memory for transfering to CUDA.
    CopyPosToMem();
    CopyIntForceToMem();
    CopyExtForceToMem();
    CopyOriToMem();
    #ifdef  TAPs_ADVANCED_SIMULATION
        //CopySimFlagsToMem();
        //CopyPosConstraintToMem();
    #endif//TAPs_ADVANCED_SIMULATION

    // DEBUG
    //printf( "BEFORE SIM -- HOST PTRs %x %x; %x %x \n", m_cudaPosList, m_cudaPrevPosList, m_cudaOriList, m_cudaPrevOriList );

    // Call the wrapper function
    TAPs::CUDA::Global__PLHMModelElasticRod_AdvSim( 
        m_cudaID,
        m_pParameters->NumOfNodes,  
        64,                         
        0,                          
        m_pParameters->TimeStep,    
        1,                          
        m_pParameters->Radius,      
        m_pParameters->LinkLength,  
        m_pParameters->LinkLength,  
        m_pParameters->MassOfPoint, 
        m_pParameters->MaterialDensity, 
        //m_pParameters->Kt,        //!< kinetic translational constant
        //m_pParameters->Kr[0],     //!< kinetic rotational constant -- x
        //m_pParameters->Kr[1],     //!< kinetic rotational constant -- y
        //m_pParameters->Kr[2],     //!< kinetic rotational constant -- z
        //m_pParameters->Dt,        //!< translational constant
        //m_pParameters->Dr[0],     //!< rotational dissipation constant -- x
        //m_pParameters->Dr[1],     //!< rotational dissipation constant -- y
        //m_pParameters->Dr[2],     //!< rotational dissipation constant -- z
        m_pParameters->Kconstraint_3rdDirAlignTangent,  
        m_pParameters->Kvdamping,   
        m_pParameters->Ps,          
        m_pParameters->Pb[0],       
        m_pParameters->Pb[1],       
        m_pParameters->Pb[2],       
        m_cudaPosList,          
        m_cudaPrevPosList,      
        m_cudaOriList,          
        m_cudaPrevOriList,      
        m_cudaIntForceList,     
        m_cudaExtForceList      
    );

    // DEBUG
    //printf( "AFTER SIM -- HOST PTRs %x %x; %x %x \n\n", m_cudaPosList, m_cudaPrevPosList, m_cudaOriList, m_cudaPrevOriList );

    // The new vertex positions are returned by CUDA to memory for previous vertex positions.
    // So the new vertex positions must be copied to the vertex positions of the object.
    //CopyMemToPrevPos();
    //CopyMemToPrevOri();
    // The swapping is done in Global__PLHMModelElasticRod_AdvSim fn
    //SwapBuffers();

    CopyMemPLHMToPos();
    CopyMemPLHMToOri();
    //SwapBuffers();

    #ifdef  TAPs_ADVANCED_SIMULATION
    // The new strand's simulation flags list is returned by CUDA.
    // It must be copied to the simulation flags list of the strand.
        //CUDA_CopyMemToSimFlags();
        //CUDA_CopyMemToPosConstraint();

        // NEED TO CHANGE THE VALUES STORED IN m_pAdvSimData

    #endif//TAPs_ADVANCED_SIMULATION
}


#if defined(__gl_h_) || defined(__GL_H__)
// Generate generalized cylinder vertex data for drawing
template <typename T>
void ElasticRod_CompByCUDA<T>::GenCylinderVertexDataForDrawing ( GLuint vbo_GL, unsigned int numOfCrossSectionVertices )
{
    if ( m_bUsePLHM ) {
        TAPs::CUDA::GL__GenCylinderForElasticRodModel_PLHM( 
            m_cudaID,
            m_pParameters->NumOfNodes,  // number of vertices
            64,                         // number of threads per block
            vbo_GL,                     // OpenGL vertex buffer object
            numOfCrossSectionVertices,  // number of cross section's vertices
            m_pParameters->Radius
        );
    }
    else {
        TAPs::CUDA::GL__GenCylinderForElasticRodModel( 
            m_cudaID,
            m_pParameters->NumOfNodes,  // number of vertices
            64,                         // number of threads per block
            vbo_GL,                     // OpenGL vertex buffer object
            numOfCrossSectionVertices,  // number of cross section's vertices
            m_pParameters->Radius
        );
    }
}
#endif//#if defined(__gl_h_) || defined(__GL_H__)

    
template <typename T>
void ElasticRod_CompByCUDA<T>::Init () throw (...)
{
    Clean();

    unsigned int size = sizeof(float) * (m_pParameters->NumOfNodes) * 4;

    // Allocate memory for the centerline's position list
    m_cudaPosList = (float *)malloc( size );
    if ( !m_cudaPosList )   throw "Couldn't allocate memory for ModelElasticRod's m_cudaPosList!";

    // Allocate memory for the centerline's previous position list
    m_cudaPrevPosList = (float *)malloc( size );
    if ( !m_cudaPrevPosList )   throw "Couldn't allocate memory for ModelElasticRod's m_cudaPrevPosList!";

    // Allocate memory for the orientation list
    m_cudaOriList = (float *)malloc( size );
    if ( !m_cudaOriList )   throw "Couldn't allocate memory for ModelElasticRod's m_cudaOriList!";

    // Allocate memory for the previous orientation list
    m_cudaPrevOriList = (float *)malloc( size );
    if ( !m_cudaPrevOriList )   throw "Couldn't allocate memory for ModelElasticRod's m_cudaPreOriList!";

    // Allocate memory for the internal force list
    m_cudaIntForceList = (float *)malloc( size );
    if ( !m_cudaIntForceList )  throw "Couldn't allocate memory for ModelElasticRod's m_cudaIntForceList!";

    // Allocate memory for the external force list
    m_cudaExtForceList = (float *)malloc( size );
    if ( !m_cudaExtForceList )  throw "Couldn't allocate memory for ModelElasticRod's m_cudaExtForceList!";

    SetOfElasticRodNodes::iterator node = m_pNodeList->begin();
    int p = 0;
    int C = *m_pIdxCurr;
    int P = *m_pIdxNext;
    while ( node != m_pNodeList->end() ) {
        // x
        m_cudaPosList[p]        = node->Centerline[C].GetPosition()[0];
        m_cudaPrevPosList[p]    = node->Centerline[P].GetPosition()[0];
        m_cudaIntForceList[p]   = node->Centerline[C].GetForce()[0];
        m_cudaExtForceList[p]   = node->ExternalForce[0];
        m_cudaOriList[p]        = node->Orientation[C].GetR();
        m_cudaPrevOriList[p++]  = node->Orientation[P].GetR();
        // y
        m_cudaPosList[p]        = node->Centerline[C].GetPosition()[1];
        m_cudaPrevPosList[p]    = node->Centerline[P].GetPosition()[1];
        m_cudaIntForceList[p]   = node->Centerline[C].GetForce()[1];
        m_cudaExtForceList[p]   = node->ExternalForce[1];
        m_cudaOriList[p]        = node->Orientation[C].GetI();
        m_cudaPrevOriList[p++]  = node->Orientation[P].GetI();
        // z
        m_cudaPosList[p]        = node->Centerline[C].GetPosition()[2];
        m_cudaPrevPosList[p]    = node->Centerline[P].GetPosition()[2];
        m_cudaIntForceList[p]   = node->Centerline[C].GetForce()[2];
        m_cudaExtForceList[p]   = node->ExternalForce[2];
        m_cudaOriList[p]        = node->Orientation[C].GetJ();
        m_cudaPrevOriList[p++]  = node->Orientation[P].GetJ();
        // w
        m_cudaPosList[p]        = 1;
        m_cudaPrevPosList[p]    = 1;
        m_cudaIntForceList[p]   = 1;
        m_cudaExtForceList[p]   = 1;
        m_cudaOriList[p]        = node->Orientation[C].GetK();
        m_cudaPrevOriList[p++]  = node->Orientation[P].GetK();
        ++node;
    }

    TAPs::CUDA::InitailizeDataForElasticRodModel(
        m_cudaID,                   
        m_pParameters->NumOfNodes,  
        m_cudaPosList,              
        m_cudaPrevPosList,          
        m_cudaOriList,              
        m_cudaPrevOriList,          
        m_cudaIntForceList,         
        m_cudaExtForceList          
    );
}


template <typename T>
void ElasticRod_CompByCUDA<T>::Clean ()
{
    if ( m_cudaPosList ) {
        free( m_cudaPosList );
        m_cudaPosList = NULL;
    }
    if ( m_cudaPrevPosList ) {
        free( m_cudaPrevPosList );
        m_cudaPrevPosList = NULL;
    }
    if ( m_cudaOriList ) {
        free( m_cudaOriList );
        m_cudaOriList = NULL;
    }
    if ( m_cudaPrevOriList ) {
        free( m_cudaPrevOriList );
        m_cudaPrevOriList = NULL;
    }
    if ( m_cudaIntForceList ) {
        free( m_cudaIntForceList );
        m_cudaIntForceList = NULL;
    }
    if ( m_cudaExtForceList ) {
        free( m_cudaExtForceList );
        m_cudaExtForceList = NULL;
    }

    if ( m_cudaID != 0 ) {
        TAPs::CUDA::ClearDataForElasticRodModel( m_cudaID );
    }
}


template <typename T>
void ElasticRod_CompByCUDA<T>::CopyPosToMem ()
{
    if ( m_cudaPosList ) {
        SetOfElasticRodNodes::iterator node = m_pNodeList->begin();
        int p = 0;
        int C = *m_pIdxCurr;
        while ( node != m_pNodeList->end() ) {
            m_cudaPosList[p++] = node->Centerline[C].GetPosition()[0];
            m_cudaPosList[p++] = node->Centerline[C].GetPosition()[1];
            m_cudaPosList[p++] = node->Centerline[C].GetPosition()[2];
            //if ( p != 0 ) {
                m_cudaPosList[p++] = 1;
            //}
            //else {
            //  m_cudaPosList[p++] = 0;
            //}
            ++node;
        }
    }
}


template <typename T>
void ElasticRod_CompByCUDA<T>::CopyPrevPosToMem ()
{
    if ( m_cudaPrevPosList ) {
        SetOfElasticRodNodes::iterator node = m_pNodeList->begin();
        int p = 0;
        int P = *m_pIdxNext;
        while ( node != m_pNodeList->end() ) {
            m_cudaPrevPosList[p++] = node->Centerline[P].GetPosition()[0];
            m_cudaPrevPosList[p++] = node->Centerline[P].GetPosition()[1];
            m_cudaPrevPosList[p++] = node->Centerline[P].GetPosition()[2];
            //if ( p != 0 ) {
                m_cudaPrevPosList[p++] = 1;
            //}
            //else {
            //  m_cudaPrevPosList[p++] = 0;
            //}
            ++node;
        }
    }
}


template <typename T>
void ElasticRod_CompByCUDA<T>::CopyOriToMem ()
{
    if ( m_cudaOriList ) {
        SetOfElasticRodNodes::iterator node = m_pNodeList->begin();
        int p = 0;
        int C = *m_pIdxCurr;
        while ( node != m_pNodeList->end() ) {
            m_cudaOriList[p++] = node->Orientation[C].GetR();
            m_cudaOriList[p++] = node->Orientation[C].GetI();
            m_cudaOriList[p++] = node->Orientation[C].GetJ();
            m_cudaOriList[p++] = node->Orientation[C].GetK();
            ++node;
        }
    }
}


template <typename T>
void ElasticRod_CompByCUDA<T>::CopyPrevOriToMem ()
{
    if ( m_cudaOriList ) {
        SetOfElasticRodNodes::iterator node = m_pNodeList->begin();
        int p = 0;
        int P = *m_pIdxNext;
        while ( node != m_pNodeList->end() ) {
            m_cudaPrevOriList[p++] = node->Orientation[P].GetR();
            m_cudaPrevOriList[p++] = node->Orientation[P].GetI();
            m_cudaPrevOriList[p++] = node->Orientation[P].GetJ();
            m_cudaPrevOriList[p++] = node->Orientation[P].GetK();
            ++node;
        }
    }
}


template <typename T>
void ElasticRod_CompByCUDA<T>::CopyIntForceToMem ()
{
    if ( m_cudaIntForceList ) {
        SetOfElasticRodNodes::iterator node = m_pNodeList->begin();
        int p = 0;
        int C = *m_pIdxCurr;
        while ( node != m_pNodeList->end() ) {
            m_cudaIntForceList[p++] = node->Centerline[C].GetForce()[0];
            m_cudaIntForceList[p++] = node->Centerline[C].GetForce()[1];
            m_cudaIntForceList[p++] = node->Centerline[C].GetForce()[2];
            m_cudaIntForceList[p++] = 1;
            ++node;
        }
    }
}


template <typename T>
void ElasticRod_CompByCUDA<T>::CopyExtForceToMem ()
{
    if ( m_cudaExtForceList ) {
        SetOfElasticRodNodes::iterator node = m_pNodeList->begin();
        int p = 0;
        int C = *m_pIdxCurr;
        while ( node != m_pNodeList->end() ) {
            m_cudaExtForceList[p++] = node->ExternalForce[0];
            m_cudaExtForceList[p++] = node->ExternalForce[1];
            m_cudaExtForceList[p++] = node->ExternalForce[2];
            m_cudaExtForceList[p++] = 1;
            ++node;
        }
    }
}


template <typename T>
void ElasticRod_CompByCUDA<T>::CopyMemToPos ()
{
    if ( m_cudaPosList ) {
        SetOfElasticRodNodes::iterator node = m_pNodeList->begin();
        int p = 0;
        int C = *m_pIdxCurr;
        while ( node != m_pNodeList->end() ) {
            node->Centerline[C].SetPosition( m_cudaPosList[p], m_cudaPosList[p+1], m_cudaPosList[p+2] );
            p += 4;
            ++node;
        }
    }
}


template <typename T>
void ElasticRod_CompByCUDA<T>::CopyMemToPrevPos ()
{
    if ( m_cudaPrevPosList ) {
        SetOfElasticRodNodes::iterator node = m_pNodeList->begin();
        int p = 0;
        int P = *m_pIdxNext;
        while ( node != m_pNodeList->end() ) {
            node->Centerline[P].SetPosition( m_cudaPrevPosList[p], m_cudaPrevPosList[p+1], m_cudaPrevPosList[p+2] );
            p += 4;
            ++node;
        }
    }
}


template <typename T>
void ElasticRod_CompByCUDA<T>::CopyMemToOri ()
{
    if ( m_cudaOriList ) {
        SetOfElasticRodNodes::iterator node = m_pNodeList->begin();
        int p = 0;
        int C = *m_pIdxCurr;
        while ( node != m_pNodeList->end() ) {
            node->Orientation[C].SetR( m_cudaOriList[p++] );
            node->Orientation[C].SetI( m_cudaOriList[p++] );
            node->Orientation[C].SetJ( m_cudaOriList[p++] );
            node->Orientation[C].SetK( m_cudaOriList[p++] );
            ++node;
        }
    }
}


template <typename T>
void ElasticRod_CompByCUDA<T>::CopyMemToPrevOri ()
{
    if ( m_cudaOriList ) {
        SetOfElasticRodNodes::iterator node = m_pNodeList->begin();
        int p = 0;
        int P = *m_pIdxNext;
        while ( node != m_pNodeList->end() ) {
            node->Orientation[P].SetR( m_cudaPrevOriList[p++] );
            node->Orientation[P].SetI( m_cudaPrevOriList[p++] );
            node->Orientation[P].SetJ( m_cudaPrevOriList[p++] );
            node->Orientation[P].SetK( m_cudaPrevOriList[p++] );
            ++node;
        }
    }
}


template <typename T>
void ElasticRod_CompByCUDA<T>::CopyMemToExtForce ()
{
    if ( m_cudaExtForceList ) {
        SetOfElasticRodNodes::iterator node = m_pNodeList->begin();
        int p = 0;
        int C = *m_pIdxCurr;
        while ( node != m_pNodeList->end() ) {
            node->ExternalForce[0] = m_cudaExtForceList[p++];
            node->ExternalForce[1] = m_cudaExtForceList[p++];
            node->ExternalForce[2] = m_cudaExtForceList[p++];
            ++p;
            ++node;
        }
    }
}


template <typename T>
void ElasticRod_CompByCUDA<T>::CopyMemToIntForce ()
{
    if ( m_cudaExtForceList ) {
        SetOfElasticRodNodes::iterator node = m_pNodeList->begin();
        int p = 0;
        int C = *m_pIdxCurr;
        while ( node != m_pNodeList->end() ) {
            node->Centerline[C].GetForce[0] = m_cudaIntForceList[p++];
            node->Centerline[C].GetForce[1] = m_cudaIntForceList[p++];
            node->Centerline[C].GetForce[2] = m_cudaIntForceList[p++];
            ++p;
            ++node;
        }
    }
}


// Initialization function for utilizing CUDA's Page-Locked Host Memory
template <typename T>
void ElasticRod_CompByCUDA<T>::InitPLHM () throw (...)
{
    CleanPLHM();

    // All of these pointers will be set to point to the allocated memory on host by cudaHostAlloc function calls
    // called by TAPs::CUDA::InitailizeDataForPLHMElasticRodModel function
    //  m_cudaPosList,              //!< list of centerlines' position
    //  m_cudaPrevPosList,          //!< list of centerlines' previous position
    //  m_cudaOriList,              //!< list of orientaions
    //  m_cudaPrevOriList,          //!< list of previous orientations
    //  m_cudaIntForceList,         //!< list of internal forces
    //  m_cudaExtForceList          //!< list of external forces
    TAPs::CUDA::InitailizeDataForPLHMElasticRodModel(
        m_cudaID,                   
        m_pParameters->NumOfNodes,  
        m_cudaPosList,              
        m_cudaPrevPosList,          
        m_cudaOriList,              
        m_cudaPrevOriList,          
        m_cudaIntForceList,         
        m_cudaExtForceList          
    );

    SetOfElasticRodNodes::iterator node = m_pNodeList->begin();
    int p = 0;
    int C = *m_pIdxCurr;
    int P = *m_pIdxNext;
    while ( node != m_pNodeList->end() ) {
        // x
        m_cudaPosList[p]        = node->Centerline[C].GetPosition()[0];
        m_cudaPrevPosList[p]    = node->Centerline[P].GetPosition()[0];
        m_cudaIntForceList[p]   = node->Centerline[C].GetForce()[0];
        m_cudaExtForceList[p]   = node->ExternalForce[0];
        m_cudaOriList[p]        = node->Orientation[C].GetR();
        m_cudaPrevOriList[p++]  = node->Orientation[P].GetR();
        // y
        m_cudaPosList[p]        = node->Centerline[C].GetPosition()[1];
        m_cudaPrevPosList[p]    = node->Centerline[P].GetPosition()[1];
        m_cudaIntForceList[p]   = node->Centerline[C].GetForce()[1];
        m_cudaExtForceList[p]   = node->ExternalForce[1];
        m_cudaOriList[p]        = node->Orientation[C].GetI();
        m_cudaPrevOriList[p++]  = node->Orientation[P].GetI();
        // z
        m_cudaPosList[p]        = node->Centerline[C].GetPosition()[2];
        m_cudaPrevPosList[p]    = node->Centerline[P].GetPosition()[2];
        m_cudaIntForceList[p]   = node->Centerline[C].GetForce()[2];
        m_cudaExtForceList[p]   = node->ExternalForce[2];
        m_cudaOriList[p]        = node->Orientation[C].GetJ();
        m_cudaPrevOriList[p++]  = node->Orientation[P].GetJ();
        // w
        m_cudaPosList[p]        = 1;
        m_cudaPrevPosList[p]    = 1;
        m_cudaIntForceList[p]   = 1;
        m_cudaExtForceList[p]   = 1;
        m_cudaOriList[p]        = node->Orientation[C].GetK();
        m_cudaPrevOriList[p++]  = node->Orientation[P].GetK();
        ++node;
    }
}


// Clean function for utilizing CUDA's Page-Locked Host Memory
template <typename T>
void ElasticRod_CompByCUDA<T>::CleanPLHM ()
{
    if ( m_cudaID != 0 ) {
        TAPs::CUDA::ClearDataForPLHMElasticRodModel( m_cudaID );
    }
}


template <typename T>
void ElasticRod_CompByCUDA<T>::CopyMemPLHMToPos ()
{
    if ( m_cudaPrevPosList ) {
        SetOfElasticRodNodes::iterator node = m_pNodeList->begin();
        int p = 0;
        int P = *m_pIdxNext;
        while ( node != m_pNodeList->end() ) {
            node->Centerline[P].SetPosition( m_cudaPosList[p], m_cudaPosList[p+1], m_cudaPosList[p+2] );
            p += 4;
            ++node;
        }
    }
}


template <typename T>
void ElasticRod_CompByCUDA<T>::CopyMemPLHMToOri ()
{
    if ( m_cudaOriList ) {
        SetOfElasticRodNodes::iterator node = m_pNodeList->begin();
        int p = 0;
        int P = *m_pIdxNext;
        while ( node != m_pNodeList->end() ) {
            node->Orientation[P].SetR( m_cudaOriList[p++] );
            node->Orientation[P].SetI( m_cudaOriList[p++] );
            node->Orientation[P].SetJ( m_cudaOriList[p++] );
            node->Orientation[P].SetK( m_cudaOriList[p++] );
            ++node;
        }
    }
}


template <typename T>
void ElasticRod_CompByCUDA<T>::SwapBuffers ()
{
    float * tmpPtr = m_cudaPosList;
    m_cudaPosList = m_cudaPrevPosList;
    m_cudaPrevPosList = tmpPtr;

    tmpPtr = m_cudaOriList;
    m_cudaOriList = m_cudaPrevOriList;
    m_cudaPrevOriList = tmpPtr;
}

//=============================================================================
END_NAMESPACE_TAPs
//-----------------------------------------------------------------------------
//34567890123456789012345678901234567890123456789012345678901234567890123456789
//--+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----