#include "TAPsCUDA_VertexListPLHMModelElasticRod.cu"

Include dependency graph for TAPsCUDA_VertexListPLHMModelElasticRod_Def.cu:

This graph shows which files directly or indirectly include this file:

Functions
BEGIN_NAMESPACE_TAPs__CUDA bool	InitailizeDataForPLHMElasticRodModel (unsigned int &cudaID, unsigned int numOfNodes, float (&posList), float (&prevPosList), float (&oriList), float (&prevOriList), float (&intForceList), float (&extForceList))
	Initialize CUDA Data for an Elastic Rod model (ModelElasticRod)
void	ClearDataForPLHMElasticRodModel (unsigned int &cudaID)
__global__ void	Global__ModelElasticRod_AdvSimPLHM_CU (unsigned int numOfNodes, unsigned int numOfThreads, float currentTime, float timeStep, float radius, float length, float length_ori, float mass, float k_mdr, float Kconstraint_3rdDirAlignTangent, float Kvdamping, float Ps, float3 Pb, float centerline, float orientation, float force_1, float force_2, float out_next_centerline, float out_next_orientation)
	CUDA Kernel -- Calculate total force of the node at vertexNo.
void	Global__PLHMModelElasticRod_AdvSim (unsigned int cudaID, unsigned int numOfNodes, unsigned int numOfThreads, float currentTime, float timeStep, int numOfSubSteps, float radius, float length, float length_ori, float mass, float material_density, float Kconstraint_3rdDirAlignTangent, float Kvdamping, float Ps, float Pb_x, float Pb_y, float Pb_z, float (&host_pos_data), float (&host_prev_pos_data), float (&host_ori_data), float (&host_prev_ori_data), float (&host_int_force_data), float (&host_ext_force_data))

Function Documentation

void ClearDataForPLHMElasticRodModel ( unsigned int & cudaID )

Clear CUDA Data for an Elastic Rod Model (ModelElasticRod) Clear CUDA Data for an Elastic Rod Model (ModelElasticRod)

Parameters:

cudaID CUDA ID for the object

Definition at line 83 of file TAPsCUDA_VertexListPLHMModelElasticRod_Def.cu.

References DATA_GlobalPool, and DATA_Pool::DataForVertexListPLHM.

Referenced by ElasticRod_CompByCUDA< T >::CleanPLHM().

{
    if ( cudaID != 0 ) {
        delete DATA_GlobalPool.DataForVertexListPLHM[cudaID];
        cudaID = 0;
    }
}

Here is the caller graph for this function:

__global__ void Global__ModelElasticRod_AdvSimPLHM_CU	(	unsigned int	numOfNodes,
		unsigned int	numOfThreads,
		float	currentTime,
		float	timeStep,
		float	radius,
		float	length,
		float	length_ori,
		float	mass,
		float	k_mdr,
		float	Kconstraint_3rdDirAlignTangent,
		float	Kvdamping,
		float	Ps,
		float3	Pb,
		float *	centerline,
		float *	orientation,
		float *	force_1,
		float *	force_2,
		float *	out_next_centerline,
		float *	out_next_orientation
	)

CUDA Kernel -- Calculate total force of the node at vertexNo.

To advance the simulation of the elastic rod model these things occur:

Device__CalForce_ModelElasticRod is called to compute the total forces of centerline and orientation.
Device__CalStretchForce_ModelElasticRod is called by Device__CalForce_ModelElasticRod to compute the stretch force of centerline.
Device__CalBendForce_ModelElasticRod is called by Device__CalForce_ModelElasticRod to compute the bend force of orientation.
Device__CalTheConstraintForce_Centerline_ModelElasticRod and Device__CalTheConstraintForce_Orientation_ModelElasticRod are called by Device__CalForce_ModelElasticRod to compute the constraint forces of centerline and orientation for aligning centerline's tangent with orientation's 3rd direction.
Device__CalDamplingVelocity_ModelElasticRod is called by Device__CalForce_ModelElasticRod to compute the damping velocity of centerline.
Device__CalAdvCenterline_ModelElasticRod is called by Device__CalForce_ModelElasticRod to advance the centerline by a semi-implicit Euler integration.
Device__CalAdvOrientation_ModelElasticRod is called by Device__CalForce_ModelElasticRod to advance the orientation by a semi-implicit Euler integration. CUDA Kernel -- For ModelElasticRod AdvSim Function

Definition at line 131 of file TAPsCUDA_VertexListPLHMModelElasticRod_Def.cu.

References Device__CalForce_ModelElasticRod(), Device__CalTorque_ModelElasticRod(), Device__EulerInt_NewOri_ModelElasticRod(), Device__EulerInt_NewPos_ModelElasticRod(), Device__EulerInt_NewVel_ModelElasticRod(), QuaternionUnit(), and XYZ().

{
    // Block index
    int bx = blockIdx.x;
    // Thread index
    int tx = threadIdx.x;
    // Vertex number
    int vertexNo = (numOfThreads * bx + tx);
    
    // DEBUG -- vertex#0 is set somewhere else or fixed
    if ( vertexNo == 0 ) return;
    
    if ( vertexNo >= numOfNodes )   return;
    // Index for output
    int idx0 = vertexNo * 4;
    int idx1 = idx0+1;
    int idx2 = idx1+1;
    int idx3 = idx2+1;

    // 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 );
    // 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 );
    
    // Get data from device's memory
    float4 curr_pos = make_float4( centerline[idx0], centerline[idx1], centerline[idx2], centerline[idx3] );
    float4 prev_pos, next_pos;
    float4 curr_ori = make_float4( orientation[idx0], orientation[idx1], orientation[idx2], orientation[idx3] );
    float4 prev_ori, next_ori;
    if ( vertexNo > 0 ) {
        prev_pos = make_float4( centerline[idx0-4], centerline[idx1-4], centerline[idx2-4], centerline[idx3-4] );
        prev_ori = make_float4( orientation[idx0-4], orientation[idx1-4], orientation[idx2-4], orientation[idx3-4] );
    }
    if ( vertexNo < numOfNodes-1 ) {
        next_pos = make_float4( centerline[idx0+4], centerline[idx1+4], centerline[idx2+4], centerline[idx3+4] );
        next_ori = make_float4( orientation[idx0+4], orientation[idx1+4], orientation[idx2+4], orientation[idx3+4] );
    }
    float4 int_force = make_float4( force_1[idx0], force_1[idx1], force_1[idx2], force_1[idx3] );
    float4 ext_force = make_float4( force_2[idx0], force_2[idx1], force_2[idx2], force_2[idx3] );

    float3 force  = make_float3( 0.0f, 0.0f, 0.0f );
    float4 torque = make_float4( 0.0f, 0.0f, 0.0f, 0.0f );

    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
        );

        // 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_next_centerline[idx0] = new_pos.x;
        out_next_centerline[idx1] = new_pos.y;
        out_next_centerline[idx2] = new_pos.z;
        out_next_centerline[idx3] = curr_pos.w;

        out_next_orientation[idx0] = new_ori.x;
        out_next_orientation[idx1] = new_ori.y;
        out_next_orientation[idx2] = new_ori.z;
        out_next_orientation[idx3] = new_ori.w;
    }
    else {
        out_next_centerline[idx0] = curr_pos.x;
        out_next_centerline[idx1] = curr_pos.y;
        out_next_centerline[idx2] = curr_pos.z;
        out_next_centerline[idx3] = curr_pos.w;

        out_next_orientation[idx0] = curr_ori.x;
        out_next_orientation[idx1] = curr_ori.y;
        out_next_orientation[idx2] = curr_ori.z;
        out_next_orientation[idx3] = curr_ori.w;
    }
}

Here is the call graph for this function:

void Global__PLHMModelElasticRod_AdvSim	(	unsigned int	cudaID,
		unsigned int	numOfNodes,
		unsigned int	numOfThreads,
		float	currentTime,
		float	timeStep,
		int	numOfSubSteps,
		float	radius,
		float	length,
		float	length_ori,
		float	mass,
		float	material_density,
		float	Kconstraint_3rdDirAlignTangent,
		float	Kvdamping,
		float	Ps,
		float	Pb_x,
		float	Pb_y,
		float	Pb_z,
		float *&	host_pos_data,
		float *&	host_prev_pos_data,
		float *&	host_ori_data,
		float *&	host_prev_ori_data,
		float *&	host_int_force_data,
		float *&	host_ext_force_data
	)

CUDA Wrapper Function -- For ModelElasticRod AdvanceSimulation Function Utilizing CUDA's Page-Locked Host Memory

Parameters:

cudaID	CUDA ID
numOfNodes	number of nodes
numOfThreads	number of threads per CUDA thread block
currentTime	current time
timeStep	time step
numOfSubSteps	number of sub-steps
radius	radius
length	rest length
length_ori	orientation's rest length
mass	point mass
material_density	material density
Kconstraint_3rdDirAlignTangent	Kconstraint for aligning centerline's tangent with orientation's 3rd direction
Kvdamping	centerline's velocity damper
Ps	potential stretch constant
Pb_x	potential bend constant -- x
Pb_y	potential bend constant -- y
Pb_z	potential bend constant -- z
host_pos_data	host's position data
host_prev_pos_data	host's previous position data
host_ori_data	host's orientation data
host_prev_ori_data	host's previous orientation data
host_int_force_data	host's (internal) force data -- xyzw
host_ext_force_data	host's (external) force data -- xyzw

Definition at line 310 of file TAPsCUDA_VertexListPLHMModelElasticRod_Def.cu.

Referenced by ElasticRod_CompByCUDA< T >::AdvSimPLHM().

{
    //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 );
    
    //printf( "090921 -- PLHM\n" );

    // Allocate device memories for out data
    unsigned int size = sizeof(float) * numOfNodes * 4;

    // // Bind Textures
    // DATA_GlobalPool.DataForVertexListPLHM[cudaID]->BindVertexList( numOfNodes );
    // DATA_GlobalPool.DataForVertexListPLHM[cudaID]->BindPrevVertexList( numOfNodes );
    // //DATA_GlobalPool.DataForVertexListPLHM[cudaID]->BindHomeVertexList( numOfNodes );
    // DATA_GlobalPool.DataForVertexListPLHM[cudaID]->BindOrientationList( numOfNodes );
    // DATA_GlobalPool.DataForVertexListPLHM[cudaID]->BindPrevOrientationList( numOfNodes );
    // DATA_GlobalPool.DataForVertexListPLHM[cudaID]->BindForceList_1( numOfNodes );
    // DATA_GlobalPool.DataForVertexListPLHM[cudaID]->BindForceList_2( numOfNodes );
    
    // 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;
    
    #if ER_TIMING_ADV_SIM_PLHM != 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 (w/ PLHM): numOfSubSteps %i\n", i );

            #if ER_TIMING_ADV_SIM_PLHM != 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_next_orientation
            );
            */

            /*
            // 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_next_centerline
            );
            */
            
            // Call the CUDA kernel (Choice 1) -- combine centerline and orientation
            Global__ModelElasticRod_AdvSimPLHM_CU<<< Dg, Db, Ns, S >>>( 
                numOfNodes, numOfThreads, currentTime, subTimeStep, 
                radius, length, length_ori, mass, k_mdr,
                Kconstraint_3rdDirAlignTangent, Kvdamping,
                Ps, Pb,
                DATA_GlobalPool.DataForVertexListPLHM[cudaID]->GetDevVertexList(),
                DATA_GlobalPool.DataForVertexListPLHM[cudaID]->GetDevOrientationList(),
                DATA_GlobalPool.DataForVertexListPLHM[cudaID]->GetDevForceList_1(),
                DATA_GlobalPool.DataForVertexListPLHM[cudaID]->GetDevForceList_2(),
                DATA_GlobalPool.DataForVertexListPLHM[cudaID]->GetDevPrevVertexList(),
                DATA_GlobalPool.DataForVertexListPLHM[cudaID]->GetDevPrevOrientationList()
            );
            
            // Swap pointers for Vertex List and Previous Vertex List
            // So that the Previous Vertex List is remain unchanged in the CUDA memory.
            DATA_GlobalPool.DataForVertexListPLHM[cudaID]->SwapVertexList();

            currentTime += subTimeStep;

            #if ER_TIMING_ADV_SIM_PLHM != 0
                cudaEventRecord( stop, 0 );
                cudaEventSynchronize( stop ); 
                cudaEventElapsedTime( &time, start, stop );
                if ( !isFirstRun ) {
                    totalTime += time;
                    ++counts;
                }
            #endif
        }
    }
    
    #if ER_TIMING_ADV_SIM_PLHM != 0
        if ( isFirstRun )   isFirstRun = false;
        cudaEventDestroy( start );
        cudaEventDestroy( stop );
        printf( "CUDA Kernal (w/ PLHM) TotalTime %f; Counts %i\n", totalTime, counts );
    #endif

    // Need this to wait for all threads to finish
    //cudaThreadSynchronize();
    cudaStreamSynchronize(S);
    
    //cudaSetDeviceFlags( cudaDeviceBlockingSync );
    
    // // Unbind Textures
    // DATA_GlobalPool.DataForVertexListPLHM[cudaID]->UnbindVertexList();
    // DATA_GlobalPool.DataForVertexListPLHM[cudaID]->UnbindPrevVertexList();
    // //DATA_GlobalPool.DataForVertexListPLHM[cudaID]->UnbindHomeVertexList();
    // DATA_GlobalPool.DataForVertexListPLHM[cudaID]->UnbindOrientationList();
    // DATA_GlobalPool.DataForVertexListPLHM[cudaID]->UnbindPrevOrientationList();
    // DATA_GlobalPool.DataForVertexListPLHM[cudaID]->UnbindForceList_1();
    // DATA_GlobalPool.DataForVertexListPLHM[cudaID]->UnbindForceList_2();

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

    // Reset the pointers to data
    host_pos_data           = DATA_GlobalPool.DataForVertexListPLHM[cudaID]->GetVertexList();
    host_prev_pos_data      = DATA_GlobalPool.DataForVertexListPLHM[cudaID]->GetPrevVertexList();
    host_ori_data           = DATA_GlobalPool.DataForVertexListPLHM[cudaID]->GetOrientationList();
    host_prev_ori_data      = DATA_GlobalPool.DataForVertexListPLHM[cudaID]->GetPrevOrientationList();
    //host_int_force_data       = DATA_GlobalPool.DataForVertexListPLHM[cudaID]->GetForceList_1();
    //host_ext_force_data       = DATA_GlobalPool.DataForVertexListPLHM[cudaID]->GetForceList_2();
}

Here is the call graph for this function:

Here is the caller graph for this function:

BEGIN_NAMESPACE_TAPs__CUDA bool InitailizeDataForPLHMElasticRodModel	(	unsigned int &	cudaID,
		unsigned int	numOfNodes,
		float *&	posList,
		float *&	prevPosList,
		float *&	oriList,
		float *&	prevOriList,
		float *&	intForceList,
		float *&	extForceList
	)

Initialize CUDA Data for an Elastic Rod model (ModelElasticRod)

Initialize CUDA Data for an Elastic Rod Model (ModelElasticRod)

Parameters:

cudaID	CUDA ID for the object
numOfNodes	number of nodes
posList	list of centerlines' position -- xyzw
prevPosList	list of centerlines' previous position -- xyzw
oriList	list of orientations -- xyzw
prevOriList	list of previous orientations -- xyzw
intForceList	list of (internal) forces -- xyzw
extForceList	list of (external) forces -- xyzw

Definition at line 17 of file TAPsCUDA_VertexListPLHMModelElasticRod_Def.cu.

References AddToGlobal_Pool_Of_DATA_Vertex_ListPLHM(), DATA_GlobalPool, DATA_Vertex_ListPLHM::GetForceList_1(), DATA_Vertex_ListPLHM::GetForceList_2(), DATA_Vertex_ListPLHM::GetOrientationList(), DATA_Vertex_ListPLHM::GetPrevOrientationList(), DATA_Vertex_ListPLHM::GetPrevVertexList(), DATA_Vertex_ListPLHM::GetVertexList(), and DATA_Pool::SizeOfGlobal_PoolPLHM.

Referenced by ElasticRod_CompByCUDA< T >::InitPLHM().

{
    int size = numOfNodes * 4 * sizeof(float);

    int version;
    cudaRuntimeGetVersion( &version );
    printf( "CUDA runtime version: %i\n", version );

    /*
    // Size of memory for (xyzw) vertices
    float * h_p, d_p;
    cudaError_t err;
    // Pinned memory allocated on the CPU
    // Report and exit the program if error
    // Get the device pointers for the pinned CPU memory mapped into the GPU memory space
    // Report and exit the program if error
    err = cudaHostAlloc( (void **)& h_p, size, cudaHostAllocMapped );
    CheckError( err, "cudaHostAlloc for posList", true );
    err = cudaHostGetDevicePointer( (void **)& d_p, (void *) h_p, 0 );
    CheckError( err, "get the device pointer for posList", true );
    //*/

    // 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_PoolPLHM;

    // Allocate CUDA data for the object.
    DATA_Vertex_ListPLHM * newData = new DATA_Vertex_ListPLHM( 
        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_ListPLHM( newData );
    
    // Set the data pointers
    posList         = newData->GetVertexList();
    prevPosList     = newData->GetPrevVertexList();
    oriList         = newData->GetOrientationList();
    prevOriList     = newData->GetPrevOrientationList();
    intForceList    = newData->GetForceList_1();
    extForceList    = newData->GetForceList_2();
    
    return true;
}

Here is the call graph for this function:

Here is the caller graph for this function:

Functions

Function Documentation