TAPsCUDA_VertexListModelElasticRod_Def.cu

Go to the documentation of this file.

/******************************************************************************
TAPsCUDA_VertexListModelElasticRod_Def.cu

CUDA Definition for ModelElasticRod (i.e., cpp file).

SUKITTI PUNAK   (09/03/2009)
UPDATE          (09/24/2009)
******************************************************************************/

#include "TAPsCUDA_VertexListModelElasticRod.cu"

BEGIN_NAMESPACE_TAPs__CUDA
//=============================================================================
// CUDA Initialization and Clear Functions
//-----------------------------------------------------------------------------
bool InitailizeDataForElasticRodModel ( 
    unsigned int & cudaID,      // CUDA ID for the object
    unsigned int numOfNodes,    // number of nodes
    float * posList,            // list of centerlines' position -- xyzw
    float * prevPosList,        // list of centerlines' previous position -- xyzw
    float * oriList,            // list of orientations -- xyzw
    float * prevOriList,        // list of previous orientations -- xyzw
    float * intForceList,       // list of (internal) forces -- xyzw
    float * extForceList        // list of (external) forces -- xyzw
)
{
    // Here the Elastic Rod Model uses Vertex List Data
    // However, it has implicit connection (each vertex connect to the previous and next vertices).
    // So the connection list is set to zero.

    // Home vertex is unused for now.  Plan is to use home vertex for sticking suture to a surface.

    // Assign CUDA ID to the object
    // The object has to use cudaID to communicate with the CUDA.
    cudaID = DATA_GlobalPool.SizeOfGlobal_Pool;

    // Allocate CUDA data for the object.
    DATA_Vertex_List * newData = new DATA_Vertex_List( 
        numOfNodes,     // total number of nodes
        false,          // include home positions into the data
        0,              // maximum vertex connection
        false,          // include simulation flags constraint into the data
        true,           // include orientations
        true,           // include previous orientations
        true,           // include forces (set 1) as for internal force
        true            // include forces (set 1) as for external force
    );
    AddToGlobal_Pool_Of_DATA_Vertex_List( newData );

    // Size of memory for (xyzw) vertices
    int size = numOfNodes * 4 * sizeof(float);

    // Copy Vertex List from host data to device data
    CUDA_SAFE_CALL( cudaMemcpy( newData->GetVertexList(), posList, size, cudaMemcpyHostToDevice ) );
    // Copy Previous Vertex List from host data to device data
    CUDA_SAFE_CALL( cudaMemcpy( newData->GetPrevVertexList(), prevPosList, size, cudaMemcpyHostToDevice ) );
    // Copy Home Vertex List from host data to device data
    //CUDA_SAFE_CALL( cudaMemcpy( newData->GetHomeVertexList(), homeVertexList, size, cudaMemcpyHostToDevice ) );
    // Copy Vertex Connection List from host data to device data
    CUDA_SAFE_CALL( cudaMemcpy( newData->GetOrientationList(), oriList, size, cudaMemcpyHostToDevice ) );
    CUDA_SAFE_CALL( cudaMemcpy( newData->GetPrevOrientationList(), prevOriList, size, cudaMemcpyHostToDevice ) );
    CUDA_SAFE_CALL( cudaMemcpy( newData->GetForceList_1(), intForceList, size, cudaMemcpyHostToDevice ) );
    CUDA_SAFE_CALL( cudaMemcpy( newData->GetForceList_2(), extForceList, size, cudaMemcpyHostToDevice ) );

    return true;
}


// Clear CUDA Data for an Elastic Rod model (ModelElasticRod)
void ClearDataForElasticRodModel ( unsigned int & cudaID )
{
    if ( cudaID != 0 ) {
        delete DATA_GlobalPool.DataForVertexList[cudaID];
        cudaID = 0;
    }
}
//-----------------------------------------------------------------------------
// CUDA Initialization and Clear Functions
//=============================================================================


/*
__global__
void Global__ModelElasticRod_AdvSim_CU_Centerline ( 
    unsigned int numOfNodes,    // number of nodes
    unsigned int numOfThreads,  // number of threads
    float currentTime,          // current time
    float timeStep,             // time step
    float radius,               // radius
    float length,               // centerline's rest length
    float length_ori,           // orientation's rest length
    float mass,                 // point mass
    float k_mdr,                // A constant value in proportion to material density and radius
    //float Kt,                 // kinetic translational constant
    //float3 Kr,                // kinetic rotational constant -- xyz
    //float Dt,                 // translational constant
    //float3 Dr,                // rotational dissipation constant -- xyz
    float Kconstraint_3rdDirAlignTangent,   // Kconstraint for aligning centerline's tangent with orientation's 3rd direction
    float Kvdamping,            // centerline's velocity damper
    float Ps,                   // potential stretch constant
    float3 Pb,                  // potential bend constant -- xyz
    float * out_data_1      // output data for centerlines
)
{
    // Block index
    int bx = blockIdx.x;
    // Thread index
    int tx = threadIdx.x;
    // Vertex number
    int vertexNo = (numOfThreads * bx + tx);
    if ( vertexNo >= numOfNodes )   return;
    // Index for output
    int idx = (numOfThreads * bx + tx) * 4;

    // Fetch data from texture linear memory
    float4 current_pos  = tex1Dfetch( CudaTexVertexList, vertexNo );
    float3 force = make_float3( 0.0f, 0.0f, 0.0f );
    
    if ( current_pos.w == 1 ) {
        force = Device__CalForce_ModelElasticRod( 
            vertexNo,       // vertex number
            numOfNodes,     // number of vertices
            radius,         // radius
            length,         // rest length
            mass,           // point mass
            //Kt,           // kinetic translational constant
            //Kr,           // kinetic rotational constant -- xyz
            //Dt,           // translational constant
            //Dr,           // rotational dissipation constant -- xyz
            Kconstraint_3rdDirAlignTangent, // Kconstraint for aligning centerline's tangent with orientation's 3rd direction
            Kvdamping,      // centerline's velocity damper
            Ps,             // potential stretch constant
            Pb              // potential bend constant -- xyz
        );

        // Use a semi Euler integration to find the new velocity and position
        float3 new_vel = Device__EulerInt_NewVel_ModelElasticRod (
            force,              // force
            make_float3(0.0f, 0.0f, 0.0f),  // velocity
            mass,               // mass
            timeStep            // time step
        );
        float3 new_pos = Device__EulerInt_NewPos_ModelElasticRod (
            new_vel,            // velocity
            XYZ(current_pos),   // position
            timeStep            // time step
        );
        
        //  WARNING: out data have to be set at the end, otherwise it won't work!!! ???
        out_data_1[idx  ] = new_pos.x;
        out_data_1[idx+1] = new_pos.y;
        out_data_1[idx+2] = new_pos.z;
        out_data_1[idx+3] = current_pos.w;
    }
    else {
        out_data_1[idx  ] = current_pos.x;
        out_data_1[idx+1] = current_pos.y;
        out_data_1[idx+2] = current_pos.z;
        out_data_1[idx+3] = current_pos.w;
    }
    // Synchronize to make sure the data are loaded
    //__syncthreads();
    
    //out_data_1[idx  ] = current_pos.x;
    //out_data_1[idx+1] = current_pos.y;
    //out_data_1[idx+2] = current_pos.z;
    //out_data_1[idx+3] = current_pos.w;
}
*/


/*
__global__
void Global__ModelElasticRod_AdvSim_CU_Orientation ( 
    unsigned int numOfNodes,    // number of nodes
    unsigned int numOfThreads,  // number of threads
    float currentTime,          // current time
    float timeStep,             // time step
    float radius,               // radius
    float length,               // centerline's rest length
    float length_ori,           // orientation's rest length
    float mass,                 // point mass
    float k_mdr,                // A constant value in proportion to material density and radius
    //float Kt,                 // kinetic translational constant
    //float3 Kr,                // kinetic rotational constant -- xyz
    //float Dt,                 // translational constant
    //float3 Dr,                // rotational dissipation constant -- xyz
    float Kconstraint_3rdDirAlignTangent,   // Kconstraint for aligning centerline's tangent with orientation's 3rd direction
    float Kvdamping,            // centerline's velocity damper
    float Ps,                   // potential stretch constant
    float3 Pb,                  // potential bend constant -- xyz
    float * out_data_2      // output data for centerlines
)
{
    // Block index
    int bx = blockIdx.x;
    // Thread index
    int tx = threadIdx.x;
    // Vertex number
    int vertexNo = (numOfThreads * bx + tx);
    if ( vertexNo >= numOfNodes - 1 )   return;
    // Index for output
    int idx = (numOfThreads * bx + tx) * 4;

    // Fetch data from texture linear memory
    float4 current_pos  = tex1Dfetch( CudaTexVertexList, vertexNo );
    float4 current_ori = tex1Dfetch( CudaTexOrientationList, vertexNo );
    float3 force = make_float3( 0.0f, 0.0f, 0.0f );
    float4 torque = make_float4( 0.0f, 0.0f, 0.0f, 0.0f );
    
    if ( current_pos.w == 1 ) {
        torque = Device__CalTorque_ModelElasticRod( 
            vertexNo,       // vertex number
            numOfNodes,     // number of vertices
            radius,         // radius
            length,         // rest length
            mass,           // point mass
            //Kt,           // kinetic translational constant
            //Kr,           // kinetic rotational constant -- xyz
            //Dt,           // translational constant
            //Dr,           // rotational dissipation constant -- xyz
            Kconstraint_3rdDirAlignTangent, // Kconstraint for aligning centerline's tangent with orientation's 3rd direction
            Kvdamping,      // centerline's velocity damper
            Ps,             // potential stretch constant
            Pb              // potential bend constant -- xyz
        );
        
        //torque = make_float4( 0.0f, 0.0f, 0.0f, 0.0f );
        
        // Use a semi Euler integration to find the new orientation
        //torque = make_float4( 0.0f, 0.0f, 0.0f, 0.0f );   // DEBUG
        float4 new_ori = Device__EulerInt_NewOri_ModelElasticRod (
            k_mdr,          // an constant value in proportion to material density and radius
            length_ori,     // orientation's rest length
            current_ori,    // orientation
            torque,         // torque in 4-dimension
            timeStep        // time step
        );
        
        // Normalize the new orientation
        //new_ori = current_ori;    // DEBUG
        new_ori = QuaternionUnit( new_ori );
        
        //  WARNING: out data have to be set at the end, otherwise it won't work!!! ???
        out_data_2[idx  ] = new_ori.x;
        out_data_2[idx+1] = new_ori.y;
        out_data_2[idx+2] = new_ori.z;
        out_data_2[idx+3] = new_ori.w;
    }
    else {
        out_data_2[idx  ] = current_ori.x;
        out_data_2[idx+1] = current_ori.y;
        out_data_2[idx+2] = current_ori.z;
        out_data_2[idx+3] = current_ori.w;
    }
    // Synchronize to make sure the data are loaded
    //__syncthreads();
}
*/

__global__
void Global__ModelElasticRod_AdvSim_CU (
    unsigned int numOfNodes,    // number of nodes
    unsigned int numOfThreads,  // number of threads
    float currentTime,          // current time
    float timeStep,             // time step
    float radius,               // radius
    float length,               // centerline's rest length
    float length_ori,           // orientation's rest length
    float mass,                 // point mass
    float k_mdr,                // A constant value in proportion to material density and radius
    //float Kt,                 // kinetic translational constant
    //float3 Kr,                // kinetic rotational constant -- xyz
    //float Dt,                 // translational constant
    //float3 Dr,                // rotational dissipation constant -- xyz
    float Kconstraint_3rdDirAlignTangent,   // Kconstraint for aligning centerline's tangent with orientation's 3rd direction
    float Kvdamping,            // centerline's velocity damper
    float Ps,                   // potential stretch constant
    float3 Pb,                  // potential bend constant -- xyz
    float * out_data_1,     // output data for centerlines
    float * out_data_2      // output data for orientations
)
{
    // Block index
    int bx = blockIdx.x;
    // Thread index
    int tx = threadIdx.x;
    // Vertex number
    int vertexNo = (numOfThreads * bx + tx);
    
    #ifdef  __DEVICE_EMULATION__
    printf( "(vertexNo: %i\n", vertexNo );
    #endif//__DEVICE_EMULATION__
    
    // DEBUG -- vertex#0 is set somewhere else or fixed
    //if ( vertexNo == 0 ) return;

    if ( vertexNo >= numOfNodes )   return;
    // Index for output
    int idx = vertexNo * 4;

    // Fetch data from texture linear memory
    float4 curr_pos  = tex1Dfetch( CudaTexVertexList, vertexNo );
    float4 curr_ori  = tex1Dfetch( CudaTexOrientationList, vertexNo );
    float4 int_force = tex1Dfetch( CudaTexForceList_1, vertexNo );
    float4 ext_force = tex1Dfetch( CudaTexForceList_2, vertexNo );
    float3 force  = make_float3( 0.0f, 0.0f, 0.0f );
    float4 torque = make_float4( 0.0f, 0.0f, 0.0f, 0.0f );
    
    float4 prev_pos = tex1Dfetch( CudaTexVertexList, vertexNo-1 );
    float4 next_pos = tex1Dfetch( CudaTexVertexList, vertexNo+1 );
    float4 prev_ori = tex1Dfetch( CudaTexOrientationList, vertexNo-1 );
    float4 next_ori = tex1Dfetch( CudaTexOrientationList, vertexNo+1 );

    #ifdef  __DEVICE_EMULATION__
    printf( "prev_pos: %f %f %f\n", prev_pos.x, prev_pos.y, prev_pos.z );
    printf( "prev_ori: %f %f %f %f\n", prev_ori.x, prev_ori.y, prev_ori.z, prev_ori.w );
    printf( "curr_pos: %f %f %f\n", curr_pos.x, curr_pos.y, curr_pos.z );
    printf( "curr_ori: %f %f %f %f\n", curr_ori.x, curr_ori.y, curr_ori.z, curr_ori.w );
    printf( "next_pos: %f %f %f\n", next_pos.x, next_pos.y, next_pos.z );
    printf( "next_ori: %f %f %f %f\n", next_ori.x, next_ori.y, next_ori.z, next_ori.w );
    printf( "int_force: %f %f %f %f\n", int_force.x, int_force.y, int_force.z, int_force.w );
    printf( "ext_force: %f %f %f %f\n", ext_force.x, ext_force.y, ext_force.z, ext_force.w );
    #endif
    
    if ( curr_pos.w == 1 ) {
        force = Device__CalForce_ModelElasticRod( 
            vertexNo,       // vertex number
            numOfNodes,     // number of vertices
            prev_pos,       // previous position
            curr_pos,       // current position
            next_pos,       // next position
            prev_ori,       // previous orientation
            curr_ori,       // current orientation
            radius,         // radius
            length,         // rest length
            mass,           // point mass
            //Kt,           // kinetic translational constant
            //Kr,           // kinetic rotational constant -- xyz
            //Dt,           // translational constant
            //Dr,           // rotational dissipation constant -- xyz
            Kconstraint_3rdDirAlignTangent, // Kconstraint for aligning centerline's tangent with orientation's 3rd direction
            Kvdamping,      // centerline's velocity damper
            Ps,             // potential stretch constant
            Pb              // potential bend constant -- xyz
        );
        
        torque = Device__CalTorque_ModelElasticRod(
            vertexNo,       // vertex number
            numOfNodes,     // number of vertices
            curr_pos,       // current position
            next_pos,       // next position
            prev_ori,       // previous orientation
            curr_ori,       // current orientation
            next_ori,       // next orientation
            radius,         // radius
            length,         // rest length
            mass,           // point mass
            //Kt,           // kinetic translational constant
            //Kr,           // kinetic rotational constant -- xyz
            //Dt,           // translational constant
            //Dr,           // rotational dissipation constant -- xyz
            Kconstraint_3rdDirAlignTangent, // Kconstraint for aligning centerline's tangent with orientation's 3rd direction
            Kvdamping,      // centerline's velocity damper
            Ps,             // potential stretch constant
            Pb              // potential bend constant -- xyz
        );
        
        force.x += ext_force.x;
        force.y += ext_force.y;
        force.z += ext_force.z;

        force.x += int_force.x;
        force.y += int_force.y;
        force.z += int_force.z;

        // Use a semi Euler integration to find the new velocity and position
        float3 new_vel = Device__EulerInt_NewVel_ModelElasticRod (
            force,          // force
            make_float3(0.0f, 0.0f, 0.0f),  // velocity
            mass,           // mass
            timeStep        // time step
        );
        float3 new_pos = Device__EulerInt_NewPos_ModelElasticRod (
            new_vel,        // velocity
            XYZ(curr_pos),  // position
            timeStep        // time step
        );
        
        #ifdef  __DEVICE_EMULATION__
        printf( "mass: %f\n", mass );
        printf( "timeStep: %f\n", timeStep );
        printf( "Force: %f %f %f\n", force.x, force.y, force.z );
        printf( "New Vel: %f %f %f\n", new_vel.x, new_vel.y, new_vel.z );
        printf( "New Pos: %f %f %f\n", new_pos.x, new_pos.y, new_pos.z );
        #endif

        // Use a semi Euler integration to find the new orientation
        float4 new_ori = Device__EulerInt_NewOri_ModelElasticRod (
            k_mdr,          // an constant value in proportion to material density and radius
            length_ori,     // orientation's rest length
            curr_ori,       // orientation
            torque,         // torque in 4-dimension
            timeStep        // time step
        );

        // Normalize the new orientation
        new_ori = QuaternionUnit( new_ori );
        
        //  WARNING: out data have to be set at the end, otherwise it won't work!!!
        out_data_1[idx  ] = new_pos.x;
        out_data_1[idx+1] = new_pos.y;
        out_data_1[idx+2] = new_pos.z;
        out_data_1[idx+3] = curr_pos.w;

        out_data_2[idx  ] = new_ori.x;
        out_data_2[idx+1] = new_ori.y;
        out_data_2[idx+2] = new_ori.z;
        out_data_2[idx+3] = new_ori.w;
    }
    else {
        out_data_1[idx  ] = curr_pos.x;
        out_data_1[idx+1] = curr_pos.y;
        out_data_1[idx+2] = curr_pos.z;
        out_data_1[idx+3] = curr_pos.w;

        out_data_2[idx  ] = curr_ori.x;
        out_data_2[idx+1] = curr_ori.y;
        out_data_2[idx+2] = curr_ori.z;
        out_data_2[idx+3] = curr_ori.w;
    }
}
//-----------------------------------------------------------------------------
// CUDA Device Functions
//=============================================================================


//=============================================================================
// CUDA Wrapper Functions
//-----------------------------------------------------------------------------
void Global__ModelElasticRod_AdvSim ( 
    unsigned int cudaID,        // CUDA ID
    unsigned int numOfNodes,    // number of nodes
    unsigned int numOfThreads,  // number of threads per CUDA thread block
    float currentTime,          // current time
    float timeStep,             // time step
    int numOfSubSteps,          // number of sub-steps
    float radius,               // radius
    float length,               // rest length
    float length_ori,           // orientation's rest length    
    float mass,                 // point mass
    float material_density,     // material density
    //float Kt,                 // kinetic translational constant
    //float Kr_x,               // kinetic rotational constant -- x
    //float Kr_y,               // kinetic rotational constant -- y
    //float Kr_z,               // kinetic rotational constant -- z
    //float Dt,                 // translational constant
    //float Dr_x,               // rotational dissipation constant -- x
    //float Dr_y,               // rotational dissipation constant -- y
    //float Dr_z,               // rotational dissipation constant -- z
    float Kconstraint_3rdDirAlignTangent,   // Kconstraint for aligning centerline's tangent with orientation's 3rd direction
    float Kvdamping,            // centerline's velocity damper
    float Ps,                   // potential stretch constant
    float Pb_x,                 // potential bend constant -- x
    float Pb_y,                 // potential bend constant -- y
    float Pb_z,                 // potential bend constant -- z
    float * host_pos_data,      // host's position data
    float * host_prev_pos_data, // host's previous position data
    float * host_ori_data,      // host's orientation data
    float * host_prev_ori_data, // host's previous orientation data
    float * host_int_force_data,    // host's (internal) force data -- xyzw
    float * host_ext_force_data     // host's (external) force data -- xyzw
)
{
    //printf( "FILE:%s LINE:%d Global__ModelElasticRod_AdvSim\n", __FILE__, __LINE__ );
    //printf( "currentTime: %g, timeStep: %g,\n", currentTime, timeStep );
    //printf( "ptMass: %g, Ks: %g, Kd: %g,\n", ptMass, Ks, Kd );
    //printf( "Number of vertices: %d\n", numOfVertices );
    //printf( "Number of threads: %d\n", numOfThreads );
    //printf( "Number of grids: %d\n", ( numOfVertices + numOfThreads - 1 ) / numOfThreads );

    // Allocate device memories for out data
    unsigned int size = sizeof(float) * numOfNodes * 4;
    float * out_data_1;     // for centerlines' position
    CUDA_SAFE_CALL( cudaMalloc( (void **)&out_data_1, size ) );
    float * out_data_2;     // for orientations
    CUDA_SAFE_CALL( cudaMalloc( (void **)&out_data_2, size ) );

    // Copy Vertex List (and Previous Vertex List) from host data to device data, 
    // since vertices are dynamically changed.
    // While home vertices and connection list are unchanged.
    CUDA_SAFE_CALL( cudaMemcpy( DATA_GlobalPool.DataForVertexList[cudaID]->GetVertexList(), host_pos_data, size, cudaMemcpyHostToDevice ) );
    // The Previous Vertex List can be saved in CUDA memory.
    // Hence, the Previous Vertex List does not have to be updated.
    //CUDA_SAFE_CALL( cudaMemcpy( DATA_GlobalPool.DataForVertexList[cudaID]->GetPrevVertexList(), host_prev_pos_data, size, cudaMemcpyHostToDevice ) );
    CUDA_SAFE_CALL( cudaMemcpy( DATA_GlobalPool.DataForVertexList[cudaID]->GetOrientationList(), host_ori_data, size, cudaMemcpyHostToDevice ) );
    CUDA_SAFE_CALL( cudaMemcpy( DATA_GlobalPool.DataForVertexList[cudaID]->GetForceList_1(), host_int_force_data, size, cudaMemcpyHostToDevice ) );
    CUDA_SAFE_CALL( cudaMemcpy( DATA_GlobalPool.DataForVertexList[cudaID]->GetForceList_2(), host_ext_force_data, size, cudaMemcpyHostToDevice ) );

    // Kernel invocation
    unsigned int numOfThreadBlocks = ( numOfNodes + numOfThreads - 1 ) / numOfThreads;
    dim3 Dg( numOfThreadBlocks, 1, 1 );
    dim3 Db( numOfThreads, 1, 1 );
    size_t Ns = 0;
    cudaStream_t S = 0;

    // Bind Textures
    DATA_GlobalPool.DataForVertexList[cudaID]->BindVertexList( numOfNodes );
    DATA_GlobalPool.DataForVertexList[cudaID]->BindPrevVertexList( numOfNodes );
    //DATA_GlobalPool.DataForVertexList[cudaID]->BindHomeVertexList( numOfNodes );
    DATA_GlobalPool.DataForVertexList[cudaID]->BindOrientationList( numOfNodes );
    DATA_GlobalPool.DataForVertexList[cudaID]->BindPrevOrientationList( numOfNodes );
    DATA_GlobalPool.DataForVertexList[cudaID]->BindForceList_1( numOfNodes );
    DATA_GlobalPool.DataForVertexList[cudaID]->BindForceList_2( numOfNodes );

    #if ER_TIMING_ADV_SIM != 0
        static bool isFirstRun = true;
        cudaEvent_t start, stop;
        float time;
        static unsigned int counts = 0;
        static float totalTime = 0.0f;
        cudaEventCreate( &start );
        cudaEventCreate( &stop );
    #endif
    
    {
        //float k_mdr = material_density * radius * radius * K_PI;
        
        float k_mdr = material_density * radius * radius * 3.1415926535897932384626433832795f;
        
        //printf( "%f = %f * %f * %f * %f\n", k_mdr, material_density, radius, radius, K_PI );

        float subTimeStep = timeStep / numOfSubSteps;
        float3 Pb = make_float3( Pb_x, Pb_y, Pb_z );
        for ( int i = 1; i <= numOfSubSteps; ++i ) {
            //printf( "CUDA: numOfSubSteps %i\n", i );
        
            #if ER_TIMING_ADV_SIM != 0
                cudaEventRecord( start, 0 );
            #endif
        
            /*
            // Call the CUDA kernel
            Global__ModelElasticRod_AdvSim_CU_Orientation<<< Dg, Db, Ns, S >>>( 
                numOfNodes, numOfThreads, currentTime, subTimeStep, 
                radius, length, length_ori, mass, k_mdr,
                Kconstraint_3rdDirAlignTangent, Kvdamping,
                Ps, Pb,
                out_data_2
            );
            */

            /*
            // Call the CUDA kernel
            Global__ModelElasticRod_AdvSim_CU_Centerline<<< Dg, Db, Ns, S >>>( 
                numOfNodes, numOfThreads, currentTime, subTimeStep, 
                radius, length, length_ori, mass, k_mdr,
                Kconstraint_3rdDirAlignTangent, Kvdamping,
                Ps, Pb,
                out_data_1
            );
            */
            
            // Call the CUDA kernel (Choice 1) -- combine centerline and orientation
            Global__ModelElasticRod_AdvSim_CU<<< Dg, Db, Ns, S >>>( 
                numOfNodes, numOfThreads, currentTime, subTimeStep, 
                radius, length, length_ori, mass, k_mdr,
                Kconstraint_3rdDirAlignTangent, Kvdamping,
                Ps, Pb,
                out_data_1, out_data_2
            );
            
            //cudaThreadSynchronize();
            
            // Swap pointers for Vertex List and Previous Vertex List
            // So that the Previous Vertex List is remain unchanged in the CUDA memory.
            DATA_GlobalPool.DataForVertexList[cudaID]->SwapVertexList();

            // Copy data to the device's memory for (current) Vertex List
            CUDA_SAFE_CALL( cudaMemcpy( DATA_GlobalPool.DataForVertexList[cudaID]->GetVertexList(), out_data_1, size, cudaMemcpyDeviceToDevice ) );
            CUDA_SAFE_CALL( cudaMemcpy( DATA_GlobalPool.DataForVertexList[cudaID]->GetOrientationList(), out_data_2, size, cudaMemcpyDeviceToDevice ) );

            currentTime += subTimeStep;
            
            #if ER_TIMING_ADV_SIM != 0
                cudaEventRecord( stop, 0 );
                cudaEventSynchronize( stop ); 
                cudaEventElapsedTime( &time, start, stop );
                if ( !isFirstRun ) {
                    totalTime += time;
                    ++counts;
                }
            #endif
        }
    }

    #if ER_TIMING_ADV_SIM != 0
        if ( isFirstRun )   isFirstRun = false;
        cudaEventDestroy( start );
        cudaEventDestroy( stop );
        printf( "CUDA Kernal TotalTime %f; Counts %i\n", totalTime, counts );
    #endif

    // Unbind Textures
    DATA_GlobalPool.DataForVertexList[cudaID]->UnbindVertexList();
    DATA_GlobalPool.DataForVertexList[cudaID]->UnbindPrevVertexList();
    //DATA_GlobalPool.DataForVertexList[cudaID]->UnbindHomeVertexList();
    DATA_GlobalPool.DataForVertexList[cudaID]->UnbindOrientationList();
    DATA_GlobalPool.DataForVertexList[cudaID]->UnbindPrevOrientationList();
    DATA_GlobalPool.DataForVertexList[cudaID]->UnbindForceList_1();
    DATA_GlobalPool.DataForVertexList[cudaID]->UnbindForceList_2();

    // Check if kernel execution generated and error
    CUT_CHECK_ERROR( "Kernel execution failed!" );

    // Copy output data to host data Previous Vertex List
    // The host will have to swap its pointers (for current and previous vertex list)
    CUDA_SAFE_CALL( cudaMemcpy( host_prev_pos_data, out_data_1, size, cudaMemcpyDeviceToHost ) );
    CUDA_SAFE_CALL( cudaMemcpy( host_prev_ori_data, out_data_2, size, cudaMemcpyDeviceToHost ) );

    //DATA_GlobalPool.DataForVertexList[cudaID]->PrintDebug( numOfVertices, host_vertex_data );

    // Free device memories
    CUDA_SAFE_CALL( cudaFree( out_data_1 ) );
    CUDA_SAFE_CALL( cudaFree( out_data_2 ) );
}
//-----------------------------------------------------------------------------
// CUDA Wrapper Functions
//=============================================================================
END_NAMESPACE_TAPs__CUDA
//-----------------------------------------------------------------------------
//34567890123456789012345678901234567890123456789012345678901234567890123456789
//--+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----