TAPsCUDA_GL_ModelElasticRod_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/18/2009)
******************************************************************************/

#include "TAPsCUDA_GL_ModelElasticRod.cu"

BEGIN_NAMESPACE_TAPs__CUDA
//=============================================================================
// CUDA for Drawing by OpenGL
//-----------------------------------------------------------------------------
//#if defined(__gl_h_) || defined(__GL_H__)

__global__
void GL__GenCylinderForElasticRodModel_CU ( 
    unsigned int numOfNodes,    
    unsigned int numOfThreads,  
    float4 *     pVertexData,   
    unsigned int numOfCrossSectionVertices, 
    float        radius         
)
{
    // 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 = vertexNo * (numOfCrossSectionVertices+1) * 4 * 3;
    
    float4 curr_ori = tex1Dfetch( CudaTexOrientationList, vertexNo );
    matrix4x4 rotMat0 = QuatenionToRotationMatrix4x4( curr_ori );
    float4 next_ori = tex1Dfetch( CudaTexOrientationList, vertexNo+1 );
    matrix4x4 rotMat1 = QuatenionToRotationMatrix4x4( next_ori );
    float4 curr_pos = tex1Dfetch( CudaTexVertexList, vertexNo );
    matrix4x4 rotMatAvg = (rotMat1 + rotMat0) / 2.0f;
    float4 next_pos = tex1Dfetch( CudaTexVertexList, vertexNo+1 );
    float4 middlePt0 = (next_pos + curr_pos) / 2.0f;
    float4 next2_pos = tex1Dfetch( CudaTexVertexList, vertexNo+2 );
    float4 middlePt1 = (next2_pos + next_pos) / 2.0f;
    
    //__shared__ float3 cs_vertices[numOfCrossSectionVertices];
    float4 vertices[3][20];
    float4 normals[3][20];
    float4 texCoords[3][20];
    
    // Compute the cross section vertices, normals, and texture coordinates for middle point from pos0 to pos1
    for ( int i = 0; i < numOfCrossSectionVertices; ++i ) {
        vertices[0][i] = rotMat0 * tex2D( CudaTexGLElasticRod, i, 0.0f );
        vertices[0][i].x *= radius;
        vertices[0][i].y *= radius;
        vertices[0][i].z *= radius;
        vertices[0][i].x += middlePt0.x;
        vertices[0][i].y += middlePt0.y;
        vertices[0][i].z += middlePt0.z;
        normals[0][i]  = rotMat0 * tex2D( CudaTexGLElasticRod, i, 1.0f );
        texCoords[0][i] = tex2D( CudaTexGLElasticRod, i, 2.0f );
    }
    // Compute the cross section vertices, normals, and texture coordinates for pos1
    for ( int i = 0; i < numOfCrossSectionVertices; ++i ) {
        vertices[1][i] = rotMatAvg * tex2D( CudaTexGLElasticRod, i, 0.0f );
        vertices[1][i].x *= radius;
        vertices[1][i].y *= radius;
        vertices[1][i].z *= radius;
        vertices[1][i].x += next_pos.x;
        vertices[1][i].y += next_pos.y;
        vertices[1][i].z += next_pos.z;
        normals[1][i]  = rotMatAvg * tex2D( CudaTexGLElasticRod, i, 1.0f );
        texCoords[1][i] = tex2D( CudaTexGLElasticRod, i, 2.0f );
        texCoords[1][i].x = 0.5f;
    }
    // Compute the cross section vertices, normals, and texture coordinates for middle point from pos1 to pos2
    for ( int i = 0; i < numOfCrossSectionVertices; ++i ) {
        vertices[2][i] = rotMat1 * tex2D( CudaTexGLElasticRod, i, 0.0f );
        vertices[2][i].x *= radius;
        vertices[2][i].y *= radius;
        vertices[2][i].z *= radius;
        vertices[2][i].x += middlePt1.x;
        vertices[2][i].y += middlePt1.y;
        vertices[2][i].z += middlePt1.z;
        normals[2][i]  = rotMat1 * tex2D( CudaTexGLElasticRod, i, 1.0f );
        texCoords[2][i] = tex2D( CudaTexGLElasticRod, i, 2.0f );
        texCoords[2][i].x = 1.0f;
    }

    {
        // Two triangle strips from middlePt0 to pos1 and from pos1 to middlePt1
        // Draw triangle strip from middlePt0 to pos1
        for ( int i = 0; i < numOfCrossSectionVertices; ++i ) {
            pVertexData[idx++] = vertices[1][i];
            pVertexData[idx++] = normals[1][i];
            pVertexData[idx++] = texCoords[1][i];

            pVertexData[idx++] = vertices[0][i];
            pVertexData[idx++] = normals[0][i];
            pVertexData[idx++] = texCoords[0][i];
        }
        pVertexData[idx++] = vertices[1][0];
        pVertexData[idx++] = normals[1][0];
        pVertexData[idx++] = texCoords[1][0];

        pVertexData[idx++] = vertices[0][0];
        pVertexData[idx++] = normals[0][0];
        pVertexData[idx++] = texCoords[0][0];

        // Draw triangle strip from pos1 to middlePt1
        for ( int i = 0; i < numOfCrossSectionVertices; ++i ) {
            pVertexData[idx++] = vertices[2][i];
            pVertexData[idx++] = normals[2][i];
            pVertexData[idx++] = texCoords[2][i];

            pVertexData[idx++] = vertices[1][i];
            pVertexData[idx++] = normals[1][i];
            pVertexData[idx++] = texCoords[1][i];

        }
        pVertexData[idx++] = vertices[2][0];
        pVertexData[idx++] = normals[2][0];
        pVertexData[idx++] = texCoords[2][0];

        pVertexData[idx++] = vertices[1][0];
        pVertexData[idx++] = normals[1][0];
        pVertexData[idx++] = texCoords[1][0];
    }
}


__global__
void GL__GenCylinderForElasticRodModel_PLHM_CU ( 
    unsigned int numOfNodes,    
    unsigned int numOfThreads,  
    float *      pCenterlines,  
    float *      pOrientations, 
    float4 *     pVertexData,   
    unsigned int numOfCrossSectionVertices, 
    float        radius         
)
{
    // 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 = vertexNo * (numOfCrossSectionVertices+1) * 4 * 3;
    
    int idx0 = vertexNo * 4;
    int idx1 = idx0+1;
    int idx2 = idx1+1;
    int idx3 = idx2+1;
    
    //float4 curr_ori = tex1Dfetch( CudaTexOrientationList, vertexNo );
    float4 curr_ori = make_float4( pOrientations[idx0], pOrientations[idx1], pOrientations[idx2], pOrientations[idx3] );
    matrix4x4 rotMat0 = QuatenionToRotationMatrix4x4( curr_ori );
    //float4 next_ori = tex1Dfetch( CudaTexOrientationList, vertexNo+1 );
    float4 next_ori = make_float4( pOrientations[idx0+4], pOrientations[idx1+4], pOrientations[idx2+4], pOrientations[idx3+4] );
    matrix4x4 rotMat1 = QuatenionToRotationMatrix4x4( next_ori );
    //float4 curr_pos = tex1Dfetch( CudaTexVertexList, vertexNo );
    float4 curr_pos = make_float4( pCenterlines[idx0], pCenterlines[idx1], pCenterlines[idx2], pCenterlines[idx3] );
    matrix4x4 rotMatAvg = (rotMat1 + rotMat0) / 2.0f;
    //float4 next_pos = tex1Dfetch( CudaTexVertexList, vertexNo+1 );
    float4 next_pos = make_float4( pCenterlines[idx0+4], pCenterlines[idx1+4], pCenterlines[idx2+4], pCenterlines[idx3+4] );
    float4 middlePt0 = (next_pos + curr_pos) / 2.0f;
    //float4 next2_pos = tex1Dfetch( CudaTexVertexList, vertexNo+2 );
    float4 next2_pos = make_float4( pCenterlines[idx0+8], pCenterlines[idx1+8], pCenterlines[idx2+8], pCenterlines[idx3+8] );
    float4 middlePt1 = (next2_pos + next_pos) / 2.0f;
    
    //__shared__ float3 cs_vertices[numOfCrossSectionVertices];
    float4 vertices[3][20];
    float4 normals[3][20];
    float4 texCoords[3][20];
    
    // Compute the cross section vertices, normals, and texture coordinates for middle point from pos0 to pos1
    for ( int i = 0; i < numOfCrossSectionVertices; ++i ) {
        vertices[0][i] = rotMat0 * tex2D( CudaTexGLElasticRod, i, 0.0f );
        vertices[0][i].x *= radius;
        vertices[0][i].y *= radius;
        vertices[0][i].z *= radius;
        vertices[0][i].x += middlePt0.x;
        vertices[0][i].y += middlePt0.y;
        vertices[0][i].z += middlePt0.z;
        normals[0][i]  = rotMat0 * tex2D( CudaTexGLElasticRod, i, 1.0f );
        texCoords[0][i] = tex2D( CudaTexGLElasticRod, i, 2.0f );
    }
    // Compute the cross section vertices, normals, and texture coordinates for pos1
    for ( int i = 0; i < numOfCrossSectionVertices; ++i ) {
        vertices[1][i] = rotMatAvg * tex2D( CudaTexGLElasticRod, i, 0.0f );
        vertices[1][i].x *= radius;
        vertices[1][i].y *= radius;
        vertices[1][i].z *= radius;
        vertices[1][i].x += next_pos.x;
        vertices[1][i].y += next_pos.y;
        vertices[1][i].z += next_pos.z;
        normals[1][i]  = rotMatAvg * tex2D( CudaTexGLElasticRod, i, 1.0f );
        texCoords[1][i] = tex2D( CudaTexGLElasticRod, i, 2.0f );
        texCoords[1][i].x = 0.5f;
    }
    // Compute the cross section vertices, normals, and texture coordinates for middle point from pos1 to pos2
    for ( int i = 0; i < numOfCrossSectionVertices; ++i ) {
        vertices[2][i] = rotMat1 * tex2D( CudaTexGLElasticRod, i, 0.0f );
        vertices[2][i].x *= radius;
        vertices[2][i].y *= radius;
        vertices[2][i].z *= radius;
        vertices[2][i].x += middlePt1.x;
        vertices[2][i].y += middlePt1.y;
        vertices[2][i].z += middlePt1.z;
        normals[2][i]  = rotMat1 * tex2D( CudaTexGLElasticRod, i, 1.0f );
        texCoords[2][i] = tex2D( CudaTexGLElasticRod, i, 2.0f );
        texCoords[2][i].x = 1.0f;
    }

    {
        // Two triangle strips from middlePt0 to pos1 and from pos1 to middlePt1
        // Draw triangle strip from middlePt0 to pos1
        for ( int i = 0; i < numOfCrossSectionVertices; ++i ) {
            pVertexData[idx++] = vertices[1][i];
            pVertexData[idx++] = normals[1][i];
            pVertexData[idx++] = texCoords[1][i];

            pVertexData[idx++] = vertices[0][i];
            pVertexData[idx++] = normals[0][i];
            pVertexData[idx++] = texCoords[0][i];
        }
        pVertexData[idx++] = vertices[1][0];
        pVertexData[idx++] = normals[1][0];
        pVertexData[idx++] = texCoords[1][0];

        pVertexData[idx++] = vertices[0][0];
        pVertexData[idx++] = normals[0][0];
        pVertexData[idx++] = texCoords[0][0];

        // Draw triangle strip from pos1 to middlePt1
        for ( int i = 0; i < numOfCrossSectionVertices; ++i ) {
            pVertexData[idx++] = vertices[2][i];
            pVertexData[idx++] = normals[2][i];
            pVertexData[idx++] = texCoords[2][i];

            pVertexData[idx++] = vertices[1][i];
            pVertexData[idx++] = normals[1][i];
            pVertexData[idx++] = texCoords[1][i];

        }
        pVertexData[idx++] = vertices[2][0];
        pVertexData[idx++] = normals[2][0];
        pVertexData[idx++] = texCoords[2][0];

        pVertexData[idx++] = vertices[1][0];
        pVertexData[idx++] = normals[1][0];
        pVertexData[idx++] = texCoords[1][0];
    }
}


// Is the texture for cross section's vertex data initialized?
bool GL__IsInitailizeDataForElasticRodModel_Drawing ()
{ return CUDA_GL_ElasticRod::IsInitialized(); }

// Initialize the CUDA/GL for elastic rod's drawing
bool GL__InitailizeDataForElasticRodModel_Drawing (
    unsigned int    numberOfCrossSectionVertices,   // number of cross section's vertices
    float *         crossSectionVertexData,         // cross section's vertex data
    GLuint          vbo_GL                          // the OpenGL buffer object for vertex data
)
{
    return Cuda_GL_ElasticRod.Create( numberOfCrossSectionVertices, crossSectionVertexData, vbo_GL );
}

// Clear the CUDA/GL for elastic rod's drawing
void GL__ClearDataForElasticRodModel_Drawing ()
{
    Cuda_GL_ElasticRod.Clear();
}

// Compute the generalized cylinder for the elastic rod
void GL__GenCylinderForElasticRodModel (
    unsigned int cudaID,        // CUDA ID
    unsigned int numOfNodes,    // number of nodes
    unsigned int numOfThreads,  // number of threads per CUDA thread block
    GLuint       vbo_GL,        // OpenGL vertex buffer object
    unsigned int numOfCrossSectionVertices, // number of cross section's vertices
    float        radius         // radius
)
{
    // 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]->BindOrientationList( numOfNodes );
    Cuda_GL_ElasticRod.BindCrossSectionVertices();
    // Map to GL
    float4 * pVertexData;
    CUDA_SAFE_CALL( cudaGLMapBufferObject( (void**)&pVertexData, vbo_GL ) );
    CUT_CHECK_ERROR( "Mapping an OpenGL buffer object FAILED!" );

    // Call the CUDA kernel
    GL__GenCylinderForElasticRodModel_CU<<< Dg, Db, Ns, S >>>( 
        numOfNodes, numOfThreads, pVertexData, numOfCrossSectionVertices, radius
    );
    CUT_CHECK_ERROR( "GenCylinderForElasticRodModel Kernel execution failed!" );

    // Unbind Textures
    DATA_GlobalPool.DataForVertexList[cudaID]->UnbindVertexList();
    DATA_GlobalPool.DataForVertexList[cudaID]->UnbindOrientationList();
    Cuda_GL_ElasticRod.UnbindCrossSectionVertices();
    // Unmap from GL
    CUDA_SAFE_CALL( cudaGLUnmapBufferObject( vbo_GL ) );
}


// Compute the generalized cylinder for the elastic rod (with CUDA's PLHM)
void GL__GenCylinderForElasticRodModel_PLHM (
    unsigned int cudaID,        // CUDA ID
    unsigned int numOfNodes,    // number of nodes
    unsigned int numOfThreads,  // number of threads per CUDA thread block
    GLuint       vbo_GL,        // OpenGL vertex buffer object
    unsigned int numOfCrossSectionVertices, // number of cross section's vertices
    float        radius         // radius
)
{
    // 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 Texture
    Cuda_GL_ElasticRod.BindCrossSectionVertices();
    // Map to GL
    float4 * pVertexData;
    CUDA_SAFE_CALL( cudaGLMapBufferObject( (void**)&pVertexData, vbo_GL ) );
    CUT_CHECK_ERROR( "Mapping an OpenGL buffer object FAILED!" );

    // Call the CUDA kernel
    GL__GenCylinderForElasticRodModel_PLHM_CU<<< Dg, Db, Ns, S >>>( 
        numOfNodes, numOfThreads,
        DATA_GlobalPool.DataForVertexListPLHM[cudaID]->GetDevVertexList(),
        DATA_GlobalPool.DataForVertexListPLHM[cudaID]->GetDevOrientationList(),
        pVertexData, numOfCrossSectionVertices, radius
    );
    CUT_CHECK_ERROR( "GenCylinderForElasticRodModel Kernel execution failed!" );

    // Unbind Texture
    Cuda_GL_ElasticRod.UnbindCrossSectionVertices();
    // Unmap from GL
    CUDA_SAFE_CALL( cudaGLUnmapBufferObject( vbo_GL ) );
}


// class GL_ElasticRod ----------------------------------------------
bool            CUDA_GL_ElasticRod::g_isInitialized = false;
cudaArray *     CUDA_GL_ElasticRod::g_Ptr_Mem_CrossSectionVertexList = NULL;
unsigned int    CUDA_GL_ElasticRod::g_numOfCrossSectionVertices = 0;
unsigned int    CUDA_GL_ElasticRod::g_sizeOfCrossSectionVertexData = 0;

// Constructor
CUDA_GL_ElasticRod::CUDA_GL_ElasticRod ()
{}

// Destructor
CUDA_GL_ElasticRod::~CUDA_GL_ElasticRod ()
{
    Clear();
}

// Create
bool CUDA_GL_ElasticRod::Create (
    unsigned int    numberOfCrossSectionVertices,   // number of cross section's vertices
    float *         crossSectionVertexData,         // cross section's vertex data
    GLuint          vbo_GL                          // the OpenGL buffer object for vertex data
)
{
    g_numOfCrossSectionVertices = numberOfCrossSectionVertices;
    if ( !g_isInitialized ) {
    
        cudaError_t error = cudaGLSetGLDevice( cutGetMaxGflopsDeviceId() );
        //cudaError_t error = cudaGLSetGLDevice( 1 );
        if ( error != cudaSuccess ) {
            const char * str = cudaGetErrorString(error);
            printf( "CUDA ERROR: %i (%s)\n", error, str );
        }
        
        printf( "cutGetMaxGflopsDeviceId(): %i\n", cutGetMaxGflopsDeviceId() );
        
        size_t dimension_size = g_numOfCrossSectionVertices * sizeof(float4);   // in bytes
        g_sizeOfCrossSectionVertexData = dimension_size * 3;    // in bytes
        cudaChannelFormatDesc desc = cudaCreateChannelDesc( 32, 32, 32, 32, cudaChannelFormatKindFloat );
        CUDA_SAFE_CALL( cudaMallocArray( &g_Ptr_Mem_CrossSectionVertexList, &desc, g_numOfCrossSectionVertices, 3 ) );
        CUT_CHECK_ERROR( "cudaMallocArray for Cross Section's Vertex Data Failed!" );
        CUDA_SAFE_CALL( cudaMemcpyToArray( g_Ptr_Mem_CrossSectionVertexList, 0, 0, crossSectionVertexData, g_sizeOfCrossSectionVertexData, cudaMemcpyHostToDevice ) );
        CUT_CHECK_ERROR( "cudaMemcpyToArray for Cross Section's Vertex Data Failed!" );
        g_isInitialized = true;
    }
    
    cudaError_t ret = cudaGLRegisterBufferObject( vbo_GL );
    if ( ret != 0 ) {
        const char * str = cudaGetErrorString(ret);
        printf( "CUDA ERROR: %i (%s)\n", ret, str );
    }
    // ret = cudaGLSetBufferObjectMapFlags( vbo_GL, cudaGLMapFlagsWriteDiscard );
    // if ( ret != 0 ) {
        // const char * str = cudaGetErrorString(ret);
        // printf( "CUDA ERROR: %i (%s)\n", ret, str );
    // }

    return true;
}

// Clear
void CUDA_GL_ElasticRod::Clear ()
{
    UnbindCrossSectionVertices();
    if ( g_isInitialized ) {
        CUDA_SAFE_CALL( cudaFreeArray( g_Ptr_Mem_CrossSectionVertexList ) );
        g_isInitialized = false;
    }
}
// class CUDA_GL_ElasticRod ----------------------------------------------

//#endif //#if defined(__gl_h_) || defined(__GL_H__)
//-----------------------------------------------------------------------------
// CUDA for Drawing by OpenGL
//=============================================================================


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