TAPsCUDA_DevFns_ModelElasticRod_Def.cu File Reference

#include "TAPsCUDA_DevFns_ModelElasticRod.cu"

Include dependency graph for TAPsCUDA_DevFns_ModelElasticRod_Def.cu:

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

Go to the source code of this file.

Functions
BEGIN_NAMESPACE_TAPs__CUDA __device__ float3	Device__EulerInt_NewVel_ModelElasticRod (float3 Force, float3 Velocity, float mass, float h)
	An explicit Euler integration for calculating new velocity.
__device__ float3	Device__EulerInt_NewPos_ModelElasticRod (float3 Velocity, float3 Position, float h)
	An explicit Euler integration for calculating new position.
__device__ float4	Device__EulerInt_NewOri_ModelElasticRod (float k_mdr, float l_rest, float4 Orientation, float4 Torque, float h)
	A semi explicit Euler integration for calculating new orientation.
__device__ float3	Device__CalStretchForce_ModelElasticRod (float Ps, float Lrest, float Lcurr, float3 Fdir)
	Calculate the stretch force of centerline.
__device__ float4	Device__CalBendTorque_ModelElasticRod (float4 O1, float4 O2, float3 Pb)
	Calculate the bend torque of orientation.
__device__ float3	Device__CalTheConstraintForce_ModelElasticRod (float3 Cen21, float l_Cen21, float inv_l_Cen21, float4 Ori, float l_rest, float Kc)
__device__ float4	Device__CalTheConstraintTorque_ModelElasticRod (float3 Cen21, float l_Cen21, float inv_l_Cen21, float4 Ori, float l_rest, float Kc)
__device__ float3	Device__CalForce_ModelElasticRod (int vertexNo, unsigned int numOfNodes, float4 Pos0, float4 Pos1, float4 Pos2, float4 Ori0, float4 Ori1, float radius, float length, float mass, float Kconstraint_3rdDirAlignTangent, float Kvdamping, float Ps, float3 Pb)
	Calculate the total force of the node at vertexNo.
__device__ float4	Device__CalTorque_ModelElasticRod (int vertexNo, unsigned int numOfNodes, float4 Pos1, float4 Pos2, float4 Ori0, float4 Ori1, float4 Ori2, float radius, float length, float mass, float Kconstraint_3rdDirAlignTangent, float Kvdamping, float Ps, float3 Pb)
	Calculate the total torque of the node at vertexNo.

---

Function Documentation

__device__ float4 Device__CalBendTorque_ModelElasticRod	(	float4	O1,
		float4	O2,
		float3	Pb
	)

Calculate the bend torque of orientation.

Parameters:

O1	orientation 1
O2	orientation 2
Pb	potential bend constant -- xyz

Definition at line 114 of file TAPsCUDA_DevFns_ModelElasticRod_Def.cu.

References Add(), InnerProduct(), Mul(), and Sub().

Referenced by Device__CalTorque_ModelElasticRod().

{
    //return make_float4( 0.0f, 0.0f, 0.0f, 0.0f );
    
    float4 OA = Add( O2, O1 );  // O2 + O1
    float4 OS = Sub( O2, O1 );  // O2 - O1
    
    // B_k(q_{j+1}+q{j})
    float4 B1 = make_float4( -OA.y,  OA.x,  OA.w, -OA.z );
    float4 B2 = make_float4( -OA.z, -OA.w,  OA.x,  OA.y );
    float4 B3 = make_float4( -OA.w,  OA.z, -OA.y,  OA.x );
    
    // 2 * K_k ( B_k(q_{j+1}+q{j}) inner_product_with (q_{j+1}-q{j}) )
    float k1 = 2.0f * Pb.x * InnerProduct( B1, OS );
    float k2 = 2.0f * Pb.y * InnerProduct( B2, OS );
    float k3 = 2.0f * Pb.z * InnerProduct( B3, OS );
    
    // Differentiation of ( B_k(q_{j+1}+q{j}) inner_product_with (q_{j+1}-q{j}) ) by q_{j}
    float4 dB1 = make_float4(  -O2.y,  O2.x,  O2.w, -O2.z );
    float4 dB2 = make_float4(  -O2.z, -O2.w,  O2.x,  O2.y );
    float4 dB3 = make_float4(  -O2.w,  O2.z, -O2.y,  O2.x );
    
    return Add( Mul( dB1, k1 ), Add( Mul( dB2, k2 ), Mul( dB3, k3 ) ) );
}

Here is the call graph for this function:

Here is the caller graph for this function:

__device__ float3 Device__CalForce_ModelElasticRod	(	int	vertexNo,
		unsigned int	numOfNodes,
		float4	Pos0,
		float4	Pos1,
		float4	Pos2,
		float4	Ori0,
		float4	Ori1,
		float	radius,
		float	length,
		float	mass,
		float	Kconstraint_3rdDirAlignTangent,
		float	Kvdamping,
		float	Ps,
		float3	Pb
	)

Calculate the total force of the node at vertexNo.

< potential stretch constant

< rest length

< current length

< difference of two position -- i.e. unpolished force

< centerline2 - centerline1

< length of (centerline2 - centerline1)

< the inverse length of (centerline2 - centerline1)

< orientation

< rest length

< Kconstraint for aligning centerline's tangent with orientation's 3rd direction

< potential stretch constant

< rest length

< current length

< difference of two position -- i.e. unpolished force

< centerline2 - centerline1

< length of (centerline2 - centerline1)

< the inverse length of (centerline2 - centerline1)

< orientation

< rest length

< Kconstraint for aligning centerline's tangent with orientation's 3rd direction

Parameters:

vertexNo	vertex number
numOfNodes	number of vertices
Pos0	previous position
Pos1	current position
Pos2	next position
Ori0	previous orientation
Ori1	current orientation
radius	radius
length	rest length
mass	point mass
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

Definition at line 217 of file TAPsCUDA_DevFns_ModelElasticRod_Def.cu.

References Device__CalStretchForce_ModelElasticRod(), Device__CalTheConstraintForce_ModelElasticRod(), Length(), Sub(), and XYZ().

Referenced by Global__ModelElasticRod_AdvSim_CU(), and Global__ModelElasticRod_AdvSimPLHM_CU().

{
    //float4 Pos1 = tex1Dfetch( CudaTexVertexList, vertexNo );      // position
    float3 P1 = XYZ( Pos1 );
    float3 total_force  = make_float3( 0.0f, 0.0f, 0.0f );
    
    if ( vertexNo > 0 ) {   // skip the first one
        // Cal the stretch force due to connection with the previous centerline
        //float4 Pos0 = tex1Dfetch( CudaTexVertexList, vertexNo-1 );    // position
        float3 P0 = XYZ( Pos0 );
        float3 P0_P1 = Sub( P0, P1 );
        
        #ifdef  __DEVICE_EMULATION__
        printf( "P1: %f %f %f\n", P1.x, P1.y, P1.z );
        printf( "P0: %f %f %f\n", P0.x, P0.y, P0.z );
        #endif
        
        float l_curr_01 = Length( P0_P1 );
        float3 force = Device__CalStretchForce_ModelElasticRod( 
            Ps,         
            length,     
            l_curr_01,  
            P0_P1       
        );
        total_force.x += force.x;
        total_force.y += force.y;
        total_force.z += force.z;
        
        #ifdef  __DEVICE_EMULATION__
        printf( "(#v>0) Fs[%i]: %f %f %f\n", vertexNo, force.x, force.y, force.z );
        #endif//__DEVICE_EMULATION__

        if ( vertexNo < numOfNodes - 2 ) {  // skip the last two
            // Cal the constraint force due to connection with the previous centerline and orientation
            //float4 Ori0 = tex1Dfetch( CudaTexOrientationList, vertexNo-1 );   // orientation
            float3 P10 = make_float3( -P0_P1.x, -P0_P1.y, -P0_P1.z );
            force = Device__CalTheConstraintForce_ModelElasticRod (
                P10,            
                l_curr_01,      
                1.0/l_curr_01,  
                Ori0,           
                length,         
                Kconstraint_3rdDirAlignTangent  
            );
            total_force.x -= force.x;
            total_force.y -= force.y;
            total_force.z -= force.z;

            #ifdef  __DEVICE_EMULATION__
            printf( "(#v>0) Fc[%i]: %f %f %f\n", vertexNo, force.x, force.y, force.z );
            #endif//__DEVICE_EMULATION__
        }
    }

    if ( vertexNo < numOfNodes - 1 ) {  // skip the last one
        // Cal the stretch force due to connection with the next centerline
        //float4 Pos2 = tex1Dfetch( CudaTexVertexList, vertexNo+1 );    // position
        float3 P2 = XYZ( Pos2 );
        float3 P2_P1 = Sub( P2, P1 );
        
        #ifdef  __DEVICE_EMULATION__
        printf( "P1: %f %f %f\n", P1.x, P1.y, P1.z );
        printf( "P2: %f %f %f\n", P2.x, P2.y, P2.z );
        #endif
        
        float l_curr_21 = Length( P2_P1 );
        float3 force = Device__CalStretchForce_ModelElasticRod( 
            Ps,         
            length,     
            l_curr_21,  
            P2_P1       
        );
        total_force.x += force.x;
        total_force.y += force.y;
        total_force.z += force.z;

        #ifdef  __DEVICE_EMULATION__
        printf( "(#v<#N-1) Fs[%i]: %f %f %f\n", vertexNo, force.x, force.y, force.z );
        #endif//__DEVICE_EMULATION__
        
        if ( vertexNo < numOfNodes - 2 ) {  // skip the last two
            // Cal the constraint force due to connection with the next centerline and orientation
            //float4 Ori1 = tex1Dfetch( CudaTexOrientationList, vertexNo ); // orientation
            force = Device__CalTheConstraintForce_ModelElasticRod (
                P2_P1,          
                l_curr_21,      
                1.0/l_curr_21,  
                Ori1,           
                length,         
                Kconstraint_3rdDirAlignTangent  
            );
            total_force.x += force.x;
            total_force.y += force.y;
            total_force.z += force.z;
            
            #ifdef  __DEVICE_EMULATION__
            printf( "(#v<#N-1) Fc[%i]: %f %f %f\n", vertexNo, force.x, force.y, force.z );
            #endif//__DEVICE_EMULATION__
        }
    }
    
    return total_force;
}

Here is the call graph for this function:

Here is the caller graph for this function:

__device__ float3 Device__CalStretchForce_ModelElasticRod	(	float	Ps,
		float	Lrest,
		float	Lcurr,
		float3	Fdir
	)

Calculate the stretch force of centerline.

Parameters:

Ps	potential stretch constant
Lrest	rest length
Lcurr	current length
Fdir	difference of two position -- i.e. unpolished force

Definition at line 85 of file TAPsCUDA_DevFns_ModelElasticRod_Def.cu.

Referenced by Device__CalForce_ModelElasticRod().

{

    #ifdef  __DEVICE_EMULATION__
    {
        float mag = (Lcurr/Lrest - 1.0f) * Ps;
        float3 Fs = make_float3( mag*Fdir.x/Lcurr, mag*Fdir.y/Lcurr, mag*Fdir.z/Lcurr );
        printf( "Ps: %f\n", Ps );
        printf( "Lrest: %f\n", Lrest );
        printf( "Lcurr: %f\n", Lcurr );
        printf( "Lcurr/Lrest: %f\n", Lcurr/Lrest );
        printf( "Lcurr/Lrest - 1.0f: %f\n", Lcurr/Lrest-1.0f );
        printf( "mag: %f\n", mag );
        printf( "Fs: %f %f %f\n", Fs.x, Fs.y, Fs.z );
    }
    #endif

    float mag = (Lcurr/Lrest - 1.0f) * Ps;
    return make_float3( mag*Fdir.x/Lcurr, mag*Fdir.y/Lcurr, mag*Fdir.z/Lcurr );
}

Here is the caller graph for this function:

__device__ float3 Device__CalTheConstraintForce_ModelElasticRod	(	float3	Cen21,
		float	l_Cen21,
		float	inv_l_Cen21,
		float4	Ori,
		float	l_rest,
		float	Kc
	)

Calculate the constraint force applied to centerline to align centerline's tangent with orientation's 3rd direction

Parameters:

Cen21	centerline2 - centerline1
l_Cen21	length of (centerline2 - centerline1)
inv_l_Cen21	the inverse length of (centerline2 - centerline1)
Ori	orientation
l_rest	rest length
Kc	Kconstraint for aligning centerline's tangent with orientation's 3rd direction

Definition at line 147 of file TAPsCUDA_DevFns_ModelElasticRod_Def.cu.

References InnerProduct(), Mul(), and Sub().

Referenced by Device__CalForce_ModelElasticRod().

{
    float l_Cen21_square = l_Cen21 * l_Cen21;
    
    // the third director of orientation
    float3 Dir3 = make_float3(
        2.0f*( Ori.y*Ori.w + Ori.x*Ori.z ),
        2.0f*( Ori.z*Ori.w - Ori.x*Ori.y ),
        Ori.x*Ori.x - Ori.y*Ori.y - Ori.z*Ori.z + Ori.w*Ori.w
    );
    
    // rest_len * Kappa * ( (r_{i+1}-r_{i})/||r_{i+1}-r_{i}|| - d3(q_{i}) )
    float3 Vc = Mul( Sub( Mul( Cen21, inv_l_Cen21 ), Dir3 ), l_rest * Kc );
    
    // Compute the constraint force for the centerline
    float3 Vcr = Mul( Vc, inv_l_Cen21 / l_Cen21_square );
    float xd = Cen21.x, yd = Cen21.y, zd = Cen21.z;
    float xd_yd = xd*yd;
    float xd_zd = xd*zd;
    float yd_zd = yd*zd;
    return make_float3(
        InnerProduct( Vcr, make_float3( l_Cen21_square - xd*xd, xd_yd, xd_zd ) ),
        InnerProduct( Vcr, make_float3( xd_yd, l_Cen21_square - yd*yd, yd_zd ) ),
        InnerProduct( Vcr, make_float3( xd_zd, yd_zd, l_Cen21_square - zd*zd ) )
    );
}

Here is the call graph for this function:

Here is the caller graph for this function:

__device__ float4 Device__CalTheConstraintTorque_ModelElasticRod	(	float3	Cen21,
		float	l_Cen21,
		float	inv_l_Cen21,
		float4	Ori,
		float	l_rest,
		float	Kc
	)

Calculate the constraint torque applied to orientation to align centerline's tangent with orientation's 3rd direction

Parameters:

Cen21	centerline2 - centerline1
l_Cen21	length of (centerline2 - centerline1)
inv_l_Cen21	the inverse length of (centerline2 - centerline1)
Ori	orientation
l_rest	rest length
Kc	Kconstraint for aligning centerline's tangent with orientation's 3rd direction

Definition at line 185 of file TAPsCUDA_DevFns_ModelElasticRod_Def.cu.

References InnerProduct(), Mul(), and Sub().

Referenced by Device__CalTorque_ModelElasticRod().

{
    // the third director of orientation
    float3 Dir3 = make_float3(
        2.0f*( Ori.y*Ori.w + Ori.x*Ori.z ),
        2.0f*( Ori.z*Ori.w - Ori.x*Ori.y ),
        Ori.x*Ori.x - Ori.y*Ori.y - Ori.z*Ori.z + Ori.w*Ori.w
    );
    
    // rest_len * Kappa * ( (r_{i+1}-r_{i})/||r_{i+1}-r_{i}|| - d3(q_{i}) )
    float3 Vc = Mul( Sub( Mul( Cen21, inv_l_Cen21 ), Dir3 ), l_rest * Kc );

    // Compute the constraint force (torque) for the orientation
    float3 Vcq = Mul( Vc, 2.0f );
    return make_float4(
        InnerProduct( Vcq, make_float3( Ori.z, -Ori.y,  Ori.x ) ),
        InnerProduct( Vcq, make_float3( Ori.w, -Ori.x, -Ori.y ) ),
        InnerProduct( Vcq, make_float3( Ori.x,  Ori.w, -Ori.z ) ),
        InnerProduct( Vcq, make_float3( Ori.y,  Ori.z,  Ori.w ) )
    );
}

Here is the call graph for this function:

Here is the caller graph for this function:

__device__ float4 Device__CalTorque_ModelElasticRod	(	int	vertexNo,
		unsigned int	numOfNodes,
		float4	Pos1,
		float4	Pos2,
		float4	Ori0,
		float4	Ori1,
		float4	Ori2,
		float	radius,
		float	length,
		float	mass,
		float	Kconstraint_3rdDirAlignTangent,
		float	Kvdamping,
		float	Ps,
		float3	Pb
	)

Calculate the total torque of the node at vertexNo.

< centerline2 - centerline1

< length of (centerline2 - centerline1)

< the inverse length of (centerline2 - centerline1)

< orientation

< rest length

< Kconstraint for aligning centerline's tangent with orientation's 3rd direction

< orientation 1

< orientation 2

< potential bend constant -- xyz

< orientation 1

< orientation 2

< potential bend constant -- xyz

Parameters:

vertexNo	vertex number
numOfNodes	number of vertices
Pos1	current position
Pos2	next position
Ori0	previous orientation
Ori1	current orientation
Ori2	next orientation
radius	radius
length	rest length
mass	point mass
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

Definition at line 343 of file TAPsCUDA_DevFns_ModelElasticRod_Def.cu.

References Device__CalBendTorque_ModelElasticRod(), Device__CalTheConstraintTorque_ModelElasticRod(), Length(), Sub(), and XYZ().

Referenced by Global__ModelElasticRod_AdvSim_CU(), and Global__ModelElasticRod_AdvSimPLHM_CU().

{
    //float4 Ori1 = tex1Dfetch( CudaTexOrientationList, vertexNo ); // orientation
    float4 total_torque = make_float4( 0.0f, 0.0f, 0.0f, 0.0f );

    //if ( vertexNo >= 0 )
    //  return total_torque;

    // Skip the last one
    if ( vertexNo < numOfNodes-1 ) 
    {
        {
            // Cal the constraint torque due to connection with the before and after centerlines
            // --- centerline1 ---- orientation1 --- centerline2 ---
            //float4 Pos1 = tex1Dfetch( CudaTexVertexList, vertexNo );  // position 
            float3 P1 = XYZ( Pos1 );
            //float4 Pos2 = tex1Dfetch( CudaTexVertexList, vertexNo+1 );    // position
            float3 P2 = XYZ( Pos2 );
            float3 P2_P1 = Sub( P2, P1 );
            float l_curr_21 = Length( P2_P1 );
            float4 torque = Device__CalTheConstraintTorque_ModelElasticRod (
                P2_P1,          
                l_curr_21,      
                1.0/l_curr_21,  
                Ori1,           
                length,         
                Kconstraint_3rdDirAlignTangent  
            );
            total_torque.x += torque.x;
            total_torque.y += torque.y;
            total_torque.z += torque.z;
            total_torque.w += torque.w;
            
            #ifdef  __DEVICE_EMULATION__
            printf( "(#v>0) Ctorque[%i]: %f %f %f %f\n", vertexNo, torque.x, torque.y, torque.z, torque.w );
            #endif//__DEVICE_EMULATION__
        }
        
        if ( vertexNo >= 0 )    // skip the first one
        {
            // Cal the bend torque due to connection with the previous orientation
            //float4 Ori0 = tex1Dfetch( CudaTexOrientationList, vertexNo-1 );   // orientation
            float4 torque = Device__CalBendTorque_ModelElasticRod(
                Ori1,   
                Ori0,   
                Pb      
            );
            total_torque.x += torque.x;
            total_torque.y += torque.y;
            total_torque.z += torque.z;
            total_torque.w += torque.w;
            
            #ifdef  __DEVICE_EMULATION__
            printf( "(#v>0) Btorque[%i]: %f %f %f %f\n", vertexNo, torque.x, torque.y, torque.z, torque.w );
            #endif//__DEVICE_EMULATION__
        }

        if ( vertexNo < numOfNodes - 2 )    // skip the last two
        {
            // Cal the bend torque due to connection with the next orientation
            //float4 Ori2 = tex1Dfetch( CudaTexOrientationList, vertexNo+1 );   // orientation
            float4 torque = Device__CalBendTorque_ModelElasticRod(
                Ori1,   
                Ori2,   
                Pb      
            );
            total_torque.z += torque.z;
            total_torque.x += torque.x;
            total_torque.w += torque.w;
            total_torque.y += torque.y;
            
            #ifdef  __DEVICE_EMULATION__
            printf( "(#v<#N-2) Btorque[%i]: %f %f %f %f\n", vertexNo, torque.x, torque.y, torque.z, torque.w );
            #endif//__DEVICE_EMULATION__
        }
    }
    return total_torque;
}

Here is the call graph for this function:

Here is the caller graph for this function:

__device__ float4 Device__EulerInt_NewOri_ModelElasticRod	(	float	k_mdr,
		float	l_rest,
		float4	Orientation,
		float4	Torque,
		float	h
	)

A semi explicit Euler integration for calculating new orientation.

Parameters:

k_mdr	an constant value in proportion to material density and radius
l_rest	orientation's rest length
Orientation	orientation
Torque	torque in 4-dimension
h	time step

Definition at line 47 of file TAPsCUDA_DevFns_ModelElasticRod_Def.cu.

References Add(), Mul(), QuaternionConjugate(), and QuaternionMul().

Referenced by Global__ModelElasticRod_AdvSim_CU(), and Global__ModelElasticRod_AdvSimPLHM_CU().

{
    // Compute the components of the inverse of the inertia tersor matrix
    // Only the diagonal components are non-zero, since the rod's material is assumed to be isotropic.
    float c = l_rest * k_mdr;
    float Jinv_0 = 4.0f / c;
    float Jinv_1 = Jinv_0;
    float Jinv_2 = 2.0f / k_mdr;

    
    // Compute the angular velocity from the 4d torque
    //   1/2 * Q_j^T * torque_in_4Dim
    //   (Q_j^T * torque_in_4Dim) equals "the multiplication of conjugation_of_q_j with torque_in_4Dim"
    float4 torque_transformed = Mul( QuaternionMul( QuaternionConjugate( Orientation ), Torque ), 0.5f );
    // Angular velocity := 1/2 * Jinv * torque_transformed_into_3Dim
    float4 angularVelocity = make_float4(
        0.0f,   // deliberately set to zero
        0.5f * Jinv_0 * torque_transformed.y, 
        0.5f * Jinv_1 * torque_transformed.z, 
        0.5f * Jinv_2 * torque_transformed.w
    );

    // Compute the velocity of the orientation
    float4 q_velocity = Mul( QuaternionMul( Orientation, angularVelocity ), 0.5f );

    // Compute and return the next orientation
    return Add( Orientation, Mul( q_velocity, h ) );
}

Here is the call graph for this function:

Here is the caller graph for this function:

__device__ float3 Device__EulerInt_NewPos_ModelElasticRod	(	float3	Velocity,
		float3	Position,
		float	h
	)

An explicit Euler integration for calculating new position.

Parameters:

Velocity	velocity
Position	position
h	time step

Definition at line 35 of file TAPsCUDA_DevFns_ModelElasticRod_Def.cu.

References Add(), and Mul().

Referenced by Global__ModelElasticRod_AdvSim_CU(), and Global__ModelElasticRod_AdvSimPLHM_CU().

{
    return Add( Position, Mul( Velocity, h ) );
}

Here is the call graph for this function:

Here is the caller graph for this function:

BEGIN_NAMESPACE_TAPs__CUDA __device__ float3 Device__EulerInt_NewVel_ModelElasticRod	(	float3	Force,
		float3	Velocity,
		float	mass,
		float	h
	)

An explicit Euler integration for calculating new velocity.

Parameters:

Force	force
Velocity	velocity
mass	mass
h	time step

Definition at line 19 of file TAPsCUDA_DevFns_ModelElasticRod_Def.cu.

References Add(), and Mul().

Referenced by Global__ModelElasticRod_AdvSim_CU(), and Global__ModelElasticRod_AdvSimPLHM_CU().

{
    // F = mA; A = F/m; V = A*h; P = V*h
    float scale = h / mass;
    float3 vel_chg = Mul( Force, scale );
    return Add( Velocity, vel_chg );
}

Here is the call graph for this function:

Here is the caller graph for this function: