TAPsSimPropForMultiPartMeshModel_HalfEdge.cpp

Go to the documentation of this file.

/******************************************************************************
TAPsSimPropForMultiPartMeshModel_HalfEdge.cpp

SimPropForMultiPartMeshModel_HalfEdge class (cpp file).

SUKITTI PUNAK   (07/25/2008)
UPDATE          (10/07/2009)
******************************************************************************/
#include "TAPsSimPropForMultiPartMeshModel_HalfEdge.hpp"
// Using Inclusion Model (i.e. definitions are included in declarations)
//                       (this name.cpp is included in name.hpp)
// Each friend is defined directly inside its declaration.

BEGIN_NAMESPACE_TAPs
//=============================================================================
//-----------------------------------------------------------------------------

//-----------------------------------------------------------------------------
template <typename T>
SimPropForMultiPartMeshModel_HalfEdge<T>::SimPropForMultiPartMeshModel_HalfEdge ( 
        T massOfPoint,                      
        T springStiffness,                  
        T springDamper,                     
        TAPs::Enum::ModelType modelType,    
        T predefinedTimeStep,               
        int numSimSubSteps,                 
        bool enableColDet,                  
        bool enableSimulation               
)
    :   SimPropForMultiPartMeshModel<T>( 
            springStiffness,
            springDamper,
            massOfPoint,
            modelType,
            predefinedTimeStep,
            numSimSubSteps,
            enableColDet,
            enableSimulation
        )
        , m_prODESolver( NULL )
    #ifdef  TAPs_USE_CUDA
        , m_cudaVertexList( NULL )
    #endif  TAPs_USE_CUDA
{}

/*
//-----------------------------------------------------------------------------
template <typename T>
SimPropForMultiPartMeshModel_HalfEdge<T>::SimPropForMultiPartMeshModel_HalfEdge ( 
    SimPropForMultiPartMeshModel_HalfEdge<T> const &obj )
{
    // Inherited Variables
    m_tKStiffness   = obj.m_tKStiffness;
    m_tKDamper      = obj.m_tKDamper;
    m_tPtMass       = obj.m_tPtMass;
    m_bColDetOn     = obj.m_bColDetOn;
    m_bSimOn        = obj.m_bSimOn;

    // Local Variables
    m_vpSpringRef       = obj.m_vpSpringRef;
    m_HEVertexPtrList   = obj.m_HEVertexPtrList;

    // SKIP ALL OF THESE VARIABLES
    //===============================================================
    // For simulation by ODE solver
    //---------------------------------------------------------------
    //TAPs::Simulation::ODESolver<T> *  m_prODESolver;  //!< ODE solver
    //Simulation::VectorSet<T> *    x0;             //!< pointer to x0 Data
    //Simulation::VectorSet<T> *    xEnd;           //!< pointer to xEnd Data
    //Simulation::VectorSet<T>      xData_1;        //!< temp data with size of m_iStateSize * 6
    //Simulation::VectorSet<T>      xData_2;        //!< temp data with size of m_iStateSize * 6
    //int   m_iStateSize;       //!< #Mass Points (#mesh vertices) where 6 = 3 positions + 3 velocities (for x,y,z)
    //int   m_iLimitStepCount;  //!< Limit step count

    m_tPredefinedTimeStep   = obj.m_tPredefinedTimeStep;
    m_iNumSimSubSteps       = obj.m_iNumSimSubSteps;

    TEST();
#ifdef  TAPs_USE_CUDA
    CUDA_Initialize();
#endif//TAPs_USE_CUDA
}
//*/

//-----------------------------------------------------------------------------
template <typename T>
SimPropForMultiPartMeshModel_HalfEdge<T>::~SimPropForMultiPartMeshModel_HalfEdge ()
{
    if ( m_prODESolver ) delete m_prODESolver;

    #ifdef  TAPs_USE_CUDA
    CUDA_Cleanup();
    #endif//TAPs_USE_CUDA
}

//-----------------------------------------------------------------------------
template <typename T>
std::string SimPropForMultiPartMeshModel_HalfEdge<T>::StrInfo () const
{
    std::ostringstream ss;
    ss  << "SimPropForMultiPartMeshModel_HalfEdge<" << typeid(T).name() << "> ==>"
        << " Spring Stiffness(" << m_tKStiffness
        << ") Damper(" << m_tKDamper
        << ") Mass of Point (" << m_tPtMass
        << ") Collision Detection (";
        if ( m_bColDetOn )  ss << "On";
        else                ss << "Off";
    ss  << ") Simulation (";
        if ( m_bSimOn )     ss << "On";
        else                ss << "Off";
    ss  << ".\n";
    return ss.str();
}

/*
//-----------------------------------------------------------------------------
template <typename T>
SimPropForMultiPartMeshModel_HalfEdge<T> & SimPropForMultiPartMeshModel_HalfEdge<T>::operator= ( 
    SimPropForMultiPartMeshModel_HalfEdge<T> const &obj )
{
    m_tKStiffness   = obj.m_tKStiffness;
    m_tKDamper      = obj.m_tKDamper;
    m_tPtMass       = obj.m_tPtMass;
    m_bColDetOn     = obj.m_bColDetOn;
    m_bSimOn        = obj.m_bSimOn;

    return *this;
}
//*/

//-----------------------------------------------------------------------------
// DxDt
template <typename T>
bool SimPropForMultiPartMeshModel_HalfEdge<T>::DxDt ( 
    T                       dt,         // I/P: time step
    Simulation::VectorSet<T> &  x,      // I/p: array  x = {pos, vel}
    Simulation::VectorSet<T> &  xdot,   // O/P: array xdot = {vel, accel}
    void *                  userData    // O/P: array of user data
)
{
    //-----------------------------------------------------
    // Convert userData ptr to ModelStrand ptr
    SimPropForMultiPartMeshModel_HalfEdge<T> *pThis = static_cast<SimPropForMultiPartMeshModel_HalfEdge<T> *>( userData );
    assert( pThis );
    pThis->ClearAllForces();
    //-----------------------------------------------------
    // Calculate and Set Spring Forces
    for ( int i = 0; i < static_cast<int>( pThis->m_vpSpringRef.size() ); ++i ) {
        pThis->m_vpSpringRef[i]->CalAndSetForce();
    }
    //-----------------------------------------------------
    pThis->DdtStateToArray( &xdot[0], dt );
    //-----------------------------------------------------
    return true;
}

//-----------------------------------------------------------------------------
// StateToArray
template <typename T>
void SimPropForMultiPartMeshModel_HalfEdge<T>::StateToArray ( T *dst )
{
    for ( int i = 0; i < m_iStateSize; ++i ) {
        *(dst++) = m_HEVertexPtrList[i]->GetParticleRef()->GetPosition()[0];
        *(dst++) = m_HEVertexPtrList[i]->GetParticleRef()->GetPosition()[1];
        *(dst++) = m_HEVertexPtrList[i]->GetParticleRef()->GetPosition()[2];
        *(dst++) = m_HEVertexPtrList[i]->GetParticleRef()->GetVelocity()[0];
        *(dst++) = m_HEVertexPtrList[i]->GetParticleRef()->GetVelocity()[1];
        *(dst++) = m_HEVertexPtrList[i]->GetParticleRef()->GetVelocity()[2];
    }
}

//-----------------------------------------------------------------------------
// ArrayToState
template <typename T>
void SimPropForMultiPartMeshModel_HalfEdge<T>::ArrayToState ( T *src )
{
    for ( int i = 0; i < m_iStateSize; ++i ) {
        m_HEVertexPtrList[i]->GetParticleRef()->SetPosition( *src, *(src+1), *(src+2) );
        src += 3;
        m_HEVertexPtrList[i]->GetParticleRef()->SetVelocity( *src, *(src+1), *(src+2) );
        src += 3;
    }
}

//-----------------------------------------------------------------------------
// DdtStateToArray
template <typename T>
void SimPropForMultiPartMeshModel_HalfEdge<T>::DdtStateToArray ( T *xdot, T dt )
{
    for ( int i = 0; i < m_iStateSize; ++i ) {
        if ( m_HEVertexPtrList[i]->GetParticleRef()->GetFixStatus() ) {
            *(xdot++) = 
            *(xdot++) = 
            *(xdot++) = 
            *(xdot++) = 
            *(xdot++) = 
            *(xdot++) = 0;
        }
        else {
            //m_prForce[i] = m_StructForce[i].gravity;
            //m_prVelocity[i][1] = m_prForce[i][1] * m_prVertex[i][1] * dt;
            //----------------------------------------
            *(xdot++) = m_HEVertexPtrList[i]->GetParticleRef()->GetVelocity()[0];
            *(xdot++) = m_HEVertexPtrList[i]->GetParticleRef()->GetVelocity()[1];
            *(xdot++) = m_HEVertexPtrList[i]->GetParticleRef()->GetVelocity()[2];
            *(xdot++) = m_HEVertexPtrList[i]->GetParticleRef()->GetForce()[0];
            *(xdot++) = m_HEVertexPtrList[i]->GetParticleRef()->GetForce()[1];
            *(xdot++) = m_HEVertexPtrList[i]->GetParticleRef()->GetForce()[2];
        }
    }
}

//-----------------------------------------------------------------------------
// ClearForces
template <typename T>
void SimPropForMultiPartMeshModel_HalfEdge<T>::ClearAllForces ()
{
    for ( int i = 0; i < static_cast<int>( m_HEVertexPtrList.size() ); ++i ) {
        m_HEVertexPtrList[i]->GetParticleRef()->SetForce( 0, 0, 0 );
        m_HEVertexPtrList[i]->GetHomeParticle().SetForce( 0, 0, 0 );
    }
}

//-----------------------------------------------------------------------------
// SetODESolver
template <typename T>
bool SimPropForMultiPartMeshModel_HalfEdge<T>::SetODESolver ( 
    TAPs::Enum::ODESolver solverMethod )
{
    // Delete the current ODE solver
    if ( m_prODESolver )    delete m_prODESolver;

    bool result = true;
    //---------------------------------------------------------------
    // For Simulation
    switch ( solverMethod ) {
        case TAPs::Enum::EULER:
            if ( ( m_prODESolver = new Simulation::ODESolverEuler<T>() ) == NULL ) {
                result = false;
            }
            break;
        case TAPs::Enum::MIDPOINT:
            if ( ( m_prODESolver = new Simulation::ODESolverMidpoint<T>() ) == NULL ) {
                result = false;
            }
            break;
        case TAPs::Enum::RUNGE_KUTTA_4:
            if ( ( m_prODESolver = new Simulation::ODESolverRungeKutta4<T>() ) == NULL ) {
                result = false;
            }
            break;
        default:
            result = false;
            break;
    }

    if ( result ) {
        m_iStateSize = static_cast<int>( m_HEVertexPtrList.size() );
        int size = m_iStateSize * 6;
        m_prODESolver->SetSize( size );
        x0   = &xData_1;
        xEnd = &xData_2;
        x0->resize( size );
        xEnd->resize( size );
        StateToArray( &(*xEnd)[0] );
    }
    else {
        #ifdef TAPs_ENABLE_DEBUG
        std::cerr   << "ERROR => HETriMeshOneModelMultiParts Constructor:" 
                    << " Cannot allocate memory for an ODE solver!" 
                    << std::endl;
        #endif
        delete this;
    }

    return result;
}

//-----------------------------------------------------------------------------
// AdvSim
template <typename T>
void SimPropForMultiPartMeshModel_HalfEdge<T>::AdvSim ( T tCurrent, T tNext )
{
    if ( IsSimOn() == false )   return;

    switch ( GetModelType() ) {
        case TAPs::Enum::MODEL_DEFORM:
            //std::cout << "\tMODEL_DEFORM\n";
            //break;

        #ifdef  TAPs_USE_CUDA
            CUDA_CopyVertexToMem();
            TAPs::CUDA::Global__SimPropForMultiPartMeshModel_HalfEdge_AdvSim( 
                m_HEVertexPtrList.size(),       
                64,                             
                tCurrent, tNext, 
                m_cudaVertexList 
            );
            CUDA_CopyMemToVertex();
        #else //TAPs_USE_CUDA
            {
                T tSubStepTime = ( tNext - tCurrent ) / m_iNumSimSubSteps;
                tNext = tCurrent + tSubStepTime;
                for ( int i = 0; i < m_iNumSimSubSteps; ++i )
                {
                    //-----------------------------------------------------
                    // Swap pointers
                    Simulation::VectorSet<T> * tmpPtr = x0;
                    x0   = xEnd;
                    xEnd = tmpPtr;
                    //-----------------------------------------------------
                    // ODE Solver
                    m_prODESolver->Solve( *x0, *xEnd, tCurrent, tNext, DxDt, this );
                    //-----------------------------------------------------
                    // Copy xEnd to object state (position and velocity of particles)
                    ArrayToState( &(*xEnd)[0] );
                    //-----------------------------------------------------
                    //SimGravity( tCurrent, tNext );
                    //-----------------------------------------------------
                    //GetBVHTree()->Update();
                    //CalAndSetNormals();
                    tCurrent = tNext;
                    tNext += tSubStepTime;
                }
            }
        #endif//TAPs_USE_CUDA
            break;

        case TAPs::Enum::MODEL_DEFORM_ELASTIC:
            //std::cout << "\tMODEL_DEFORM_ELASTIC\n";
            //break;
            {
                std::vector< HEVertexPtrWithExtraInfo<T> * >::iterator  pos = m_HEVertexPtrList.begin();
                while ( pos != m_HEVertexPtrList.end() ) {
                    //(*pos)->SetHomePosition( (*pos)->GetHomePosition() + Vector3<T>(0,0,-0.001) );
                    (*pos)->GetParticleRef()->SetPosition( (*pos)->GetHomePosition() );
                    ++pos;
                }
            }
            break;
        case TAPs::Enum::MODEL_DEFORM_PLASTIC:
            std::cout << "\tMODEL_DEFORM_PLASTIC -- NOT IMPLEMENTED YET\n";
            break;
        case TAPs::Enum::MODEL_RIGID:
            std::cout << "\tMODEL_RIGID -- NOT IMPLEMENTED YET\n";
            break;
        case TAPs::Enum::MODEL_FIXED:
            std::cout << "\tMODEL_FIXED -- NOT IMPLEMENTED YET\n";
            break;
    }
}

//-----------------------------------------------------------------------------


#ifdef  TAPs_USE_CUDA
//-----------------------------------------------------------------------------

template <typename T>
bool SimPropForMultiPartMeshModel_HalfEdge<T>::CUDA_Initialize ()
{
    CUDA_Cleanup();

    //m_cudaVertexList = new float[ m_HEVertexPtrList.size() * 3 ];
    unsigned int size = sizeof(float) * m_HEVertexPtrList.size() * 3;
    m_cudaVertexList = (float *)malloc( size );

    std::cout << "SimPropForMultiPartMeshModel_HalfEdge -- m_cudaVertexList size: " << m_HEVertexPtrList.size() * 3 << "\n";

    if ( !m_cudaVertexList ) {
        std::cout << "ERROR: Could not allocate memory for CUDA!" << std::endl;
        exit( -1 );
    }
    return true;    
}

template <typename T>
void SimPropForMultiPartMeshModel_HalfEdge<T>::CUDA_Cleanup ()
{
    if ( m_cudaVertexList ) {
        //delete [] m_cudaVertexList;
        free( m_cudaVertexList );
        m_cudaVertexList = NULL;
    }
}

template <typename T>
void SimPropForMultiPartMeshModel_HalfEdge<T>::CUDA_CopyVertexToMem ()
{
    if ( !m_cudaVertexList )    return;
    for ( unsigned int i = 0, p = 0; i < static_cast<unsigned int>( m_HEVertexPtrList.size() ); ++i ) {
        m_cudaVertexList[p++] = m_HEVertexPtrList[i]->GetHEVertexPtr()->GetPosition()[0];
        m_cudaVertexList[p++] = m_HEVertexPtrList[i]->GetHEVertexPtr()->GetPosition()[1];
        m_cudaVertexList[p++] = m_HEVertexPtrList[i]->GetHEVertexPtr()->GetPosition()[2];
    }
}

template <typename T>
void SimPropForMultiPartMeshModel_HalfEdge<T>::CUDA_CopyMemToVertex ()
{
    if ( !m_cudaVertexList )    return;
    for ( unsigned int i = 0, p = 0; i < static_cast<unsigned int>( m_HEVertexPtrList.size() ); ++i ) {
        m_HEVertexPtrList[i]->GetHEVertexPtr()->SetPosition( 
            m_cudaVertexList[p],
            m_cudaVertexList[p+1],
            m_cudaVertexList[p+2]
        );
        p += 3;
    }
}

//-----------------------------------------------------------------------------
#endif//TAPs_USE_CUDA


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