TAPsSpringHalfEdgeModel.cpp

Go to the documentation of this file.

/******************************************************************************
TAPsSpringHalfEdgeModel.cpp

A Polygonal Model based on half-edge data structure.
(See TAPsSpringHalfEdgeModel.hpp for details)

SUKITTI PUNAK   (02/27/2006)
UPDATE          (04/10/2006)
******************************************************************************/
#include "TAPsSpringHalfEdgeModel.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__OpenGL
//=============================================================================
//-----------------------------------------------------------------------------
// default constructor
template <typename T>
SpringHalfEdgeModel<T>::SpringHalfEdgeModel ()
    : HalfEdgeModel<T>(), 
      m_pviListVertexRing1( NULL ), 
      m_pviListConnectedVertex( NULL ), 
      m_pviListSpring( NULL ), 
      m_iRigidity( 1 ), 
      m_tPtMass( 1 ), 
      m_tKStiffness( 1 ), 
      m_tKDamper( 1 ) 
{
    //----------------------------------------------------------------
    //std::cout << "Number of Vertex: " << GetNoVertices();
    //std::cout << "Number of Vertex: " << GetNoVertices();
    #ifdef  TAPs_ENABLE_DEBUG
    std::cout << "SpringHalfEdgeModel<" << typeid(T).name() << "> Constructor\n";
    #endif//TAPs_ENABLE_DEBUG
}
//-----------------------------------------------------------------------------
// destructor
template <typename T>
SpringHalfEdgeModel<T>::~SpringHalfEdgeModel ()
{
    DeleteAllData();
    //----------------------------------------------------------------
    // For simulation
    delete m_prODESolver;
    //----------------------------------------------------------------
    #ifdef  TAPs_ENABLE_DEBUG
    std::cout << "SpringHalfEdgeModel<" << typeid(T).name() << "> Destructor\n";
    #endif//TAPs_ENABLE_DEBUG
}
//-----------------------------------------------------------------------------
// Initialize
// a virtual fn from Model class
template <typename T>
void SpringHalfEdgeModel<T>::Initialize ()
{
    HalfEdgeModel<T>::Initialize();
    //--------------------------------------------------------------------
    #ifdef  TAPs_ENABLE_DEBUG
    std::cout << "++++++++++++++Initialize HALF-EDGE SPRING MODEL++++++++++++++\n";
    #endif//TAPs_ENABLE_DEBUG
}
//-----------------------------------------------------------------------------
// Create Mesh Springs
template <typename T>
void SpringHalfEdgeModel<T>::CreateMeshSprings ( 
            int rigidity, T mass, T Ks, T Kd )
{
    PrepareData();
    //--------------------------------------------------------------------
    // For Simulation
    if ( ( m_prODESolver = new Simulation::ODESolverEuler<T>() ) == NULL ) {
    //if ( ( m_prODESolver = new Simulation::ODESolverMidpoint<T>() ) == NULL ) {
    //if ( ( m_prODESolver = new Simulation::ODESolverRungeKutta4<T>() ) == NULL ) {
        #ifdef TAPs_ENABLE_DEBUG
        std::cerr   << "ERROR => SpringHalfEdgeModel Constructor:" 
                    << " Cannot allocate memory for an ODE solver!" 
                    << std::endl;
        #endif
        delete this;
        return;
    }
    //--------------------------------------------------------------------
    m_iStateSize = static_cast<int>( m_vpParticleRef.size() );
    m_prODESolver->SetSize( m_iStateSize * 6 );
    x0   = &xData_1;
    xEnd = &xData_2;
    x0->resize(   m_iStateSize * 6 );
    xEnd->resize( m_iStateSize * 6 );
    StateToArray( &(*xEnd)[0] );
    //--------------------------------------------------------------------
    #ifdef  TAPs_ENABLE_DEBUG
    std::cout << "++++++++++++++HALF-EDGE SPRING MODEL++++++++++++++\n";
    std::cout << "++++++++++++++ --> m_iStateSize = " << m_iStateSize << "\n";
    std::cout << "++++++++++++++HALF-EDGE SPRING MODEL++++++++++++++\n";
    #endif//TAPs_ENABLE_DEBUG
    //--------------------------------------------------------------------
    SetListOf1RingVertices();
    SetMassOfParticles( mass );
    CreateSprings( Ks, Kd, rigidity );  // A recursive method for creating springs

    #ifdef  TAPs_ENABLE_DEBUG
    std::cout << "Spring Mesh Rigidity: " << m_iRigidity << "\n";
    std::cout << "The Total Number of Springs is " << static_cast<int>(m_vpSpringRef.size()) << "\n";
    #endif//TAPs_ENABLE_DEBUG

    /*
    //--------------------------------------------------------------------
    // DEBUG PRINTOUT
    for ( int i = 0; i < static_cast<int>( m_vpParticleRef.size() ); ++i ) {
        std::cout << *m_vpParticleRef[i] << "\n";
    }
    for ( int i = 0; i < static_cast<int>( m_vpSpringRef.size() ); ++i ) {
        std::cout << *m_vpSpringRef[i] << "\n";
    }
    //*/

    //--------------------------------------------------------------------
    /*
    for ( int i = 0; i < GetNoVertices(); ++i ) {
        std::cout << "Vertex#" << i << ":";
        for ( int j = 0; j < static_cast<int>(m_pviListConnectedVertex[i].size()); ++j ) {
            std::cout << " " << m_pviListConnectedVertex[i][j];
        }
        std::cout << "\n";
    }
    //*/
}
//-----------------------------------------------------------------------------
//******************************************************************************

//******************************************************************************
// Helper Fn(s)
//-----------------------------------------------------------------------------
template <typename T>
void SpringHalfEdgeModel<T>::ChangeConstantStiffnessOfAllSprings ( T val )
{
    std::vector< SpringRef<T> * >::iterator pos;
    for ( pos = m_vpSpringRef.begin(); pos < m_vpSpringRef.end(); ++pos ) {
        (*pos)->SetConstantK( val );
    }
}
//-----------------------------------------------------------------------------
template <typename T>
void SpringHalfEdgeModel<T>::ChangeConstantDampingOfAllSprings ( T val )
{
    std::vector< SpringRef<T> * >::iterator pos;
    for ( pos = m_vpSpringRef.begin(); pos < m_vpSpringRef.end(); ++pos ) {
        (*pos)->SetDampingD( val );
    }
}
//-----------------------------------------------------------------------------
template <typename T>
void SpringHalfEdgeModel<T>::ChangeMassOfAllPoints ( T val )
{
    std::vector< ParticleRef<T> * >::iterator pos;
    for ( pos = m_vpParticleRef.begin(); pos < m_vpParticleRef.end(); ++pos ) {
        (*pos)->SetMass( val );
    }
}
//-----------------------------------------------------------------------------
//******************************************************************************

//******************************************************************************
// Helper Fn(s)
//-----------------------------------------------------------------------------
template <typename T>
void SpringHalfEdgeModel<T>::DeleteAllData ()
{
    //m_ptrHEVertex.clear();
    delete [] m_pviListVertexRing1;
    delete [] m_pviListConnectedVertex;
    delete [] m_pviListSpring;
    //m_listParticleRef.clear();
    //m_listSpringRef.clear();
    //m_ptrSpringRef.clear();
    for ( int i = 0; i < static_cast<int>( m_vpParticleRef.size() ); ++i ) {
        delete m_vpParticleRef[i];
    }
    for ( int i = 0; i < static_cast<int>( m_vpSpringRef.size() ); ++i ) {
        delete m_vpSpringRef[i];
    }
}
//-----------------------------------------------------------------------------
template <typename T>
void SpringHalfEdgeModel<T>::SetListOf1RingVertices ()
{
    HEHalfEdge<T> * firstHalfEdge, * currentHalfEdge;
    HEVertex<T> *   heVertex = GetVertexList()->Head();
    HEVertex<T> *   currentVertex;

    #ifdef  TAPs_ENABLE_DEBUG
    std::cout << "Progress: ";
    #endif//TAPs_ENABLE_DEBUG

    #ifdef  TAPs_ENABLE_DEBUG
    int count = 0;
    #endif//TAPs_ENABLE_DEBUG

    int size = static_cast<int>( m_vpHEVertex.size() );
    for ( int v = 0; v < size; ++v ) {

        #ifdef  TAPs_ENABLE_DEBUG
        if ( ++count % 100 == 0 )   std::cout << ".";
        #endif//TAPs_ENABLE_DEBUG

        currentHalfEdge = firstHalfEdge = heVertex->IncidentHalfEdge();
        do {
            currentVertex = currentHalfEdge->Pair()->Vertex();
            for ( int i = 0; i < size; ++i ) {
                if ( currentVertex == m_vpHEVertex[i] ) {
                    m_pviListVertexRing1[v].push_back( i );
                    break;
                }
            }
            currentHalfEdge = currentHalfEdge->Pair()->Next();
        } while ( firstHalfEdge != currentHalfEdge );
        heVertex = heVertex->Next();

        //std::cout << "One Ring Size " 
        //  << static_cast<int>( m_pviListVertexRing1[v].size() ) << ":";
        //for ( int i = 0; i < static_cast<int>( m_pviListVertexRing1[v].size() ); ++i ) {
        //  std::cout << " " << m_pviListVertexRing1[v][i];
        //}
        //std::cout << "\n";
    }

    #ifdef  TAPs_ENABLE_DEBUG
    std::cout << " Done\n";
    #endif//TAPs_ENABLE_DEBUG
}
//-----------------------------------------------------------------------------
template <typename T>
void SpringHalfEdgeModel<T>::PrepareData ()
{
    int size = GetVertexList()->Size();
    if ( m_pviListVertexRing1 )     delete [] m_pviListVertexRing1;
    if ( m_pviListConnectedVertex ) delete [] m_pviListConnectedVertex;
    if ( m_pviListSpring )          delete [] m_pviListSpring;
    m_pviListVertexRing1     = new std::vector<int>[ size ];
    m_pviListConnectedVertex = new std::vector<int>[ size ];
    m_pviListSpring          = new std::vector<int>[ size ];
    //-----------------------------------------------------------------------
    // Setup pointers to HEVertex
    HEVertex<T> * heVertex = GetVertexList()->Head();
    while ( heVertex ) {
        m_vpHEVertex.push_back( heVertex );
        m_vpParticleRef.push_back( new ParticleRef<T>( heVertex->GetProtectedPosition() ) );
        //m_vpParticleRef.push_back( new ParticleRef<T>( const_cast< Vector3<T> & >( heVertex->GetPosition() ) ) );
        heVertex = heVertex->Next();
    }
}
//-----------------------------------------------------------------------------
template <typename T>
void SpringHalfEdgeModel<T>::SetMassOfParticles ( T mass )
{
    //std::list< ParticleRef<T> >::iterator pos;
    //for ( pos = m_listParticleRef.begin(); pos != m_listParticleRef.end(); ++pos ) {
    //  pos->SetMass( mass );
    //}
    for ( int i = 0; i < static_cast<int>( m_vpParticleRef.size() ); ++i ) {
        m_vpParticleRef[i]->SetMass( mass );
    }
}
//-----------------------------------------------------------------------------
template <typename T>
bool SpringHalfEdgeModel<T>::CheckASpringExistance ( int p1, int p2 )
{
    int size = static_cast<int>( m_pviListConnectedVertex[p1].size() );
    for ( int i = 0; i < size; ++i ) {
        if ( m_pviListConnectedVertex[p1][i] == p2 )    return true;
    }
    return false;
}
//-----------------------------------------------------------------------------
template <typename T>
void SpringHalfEdgeModel<T>::CreateASpring ( 
                T Ks, T Kd, int p1, int p2, int hops )
{
    int pos = static_cast<int>( m_vpSpringRef.size() );
    m_vpSpringRef.push_back( 
        new SpringRef<T>( 
                *(m_vpParticleRef[p1]), 
                *(m_vpParticleRef[p2]), 
                Ks, Kd, 
                hops ) );
    m_pviListConnectedVertex[p1].push_back( p2 );
    m_pviListConnectedVertex[p2].push_back( p1 );
    m_pviListSpring[p1].push_back( pos );
    m_pviListSpring[p2].push_back( pos );
    //std::cout << "# of Hops: " << hops << "\n";
    //std::cout << "Create Spring with " << hops << " hop(s)\n";
}
//-----------------------------------------------------------------------------
template <typename T>
void SpringHalfEdgeModel<T>::CreateSprings ( T Ks, T Kd, int rigidity )
{
    m_iRigidity = rigidity;
    int size = static_cast<int>( m_vpParticleRef.size() );
    bool * isMark = new bool[size];
    for ( int i = 0; i < size; ++i ) {
        isMark[i] = false;
    }
    //--------------------------------------------------------------------
    for ( int i = 0; i < size; ++i ) {
        std::vector<int> particleSet;       // start as empty set
        //std::cout << static_cast<int>(particleSet.size()) << "\n";
        BFSCreateSprings(   Ks, Kd, 0, 
                            i, i,           // particle#1, particle#2
                            particleSet, 
                            isMark );
        for ( int j = 0; j < static_cast<int>( particleSet.size() ); ++j ) {
            isMark[ particleSet[j] ] = false;
        }
    }
}
//-----------------------------------------------------------------------------
template <typename T>
void SpringHalfEdgeModel<T>::BFSCreateSprings ( 
                T Ks, T Kd, int rigidity,       // Ks, Kd, & hop distance
                int p1, int p2,                 // particle#1 & #2
                std::vector<int> & particleSet, // set of connected particles
                bool * isMark )                 // marked particles in particleSet
{
    if ( rigidity >= m_iRigidity ) return;
    //--------------------------------------------------------------------
    // Keep the particle reference in particleSet
    // And mark it as visited
    particleSet.push_back( p2 );
    isMark[ p2 ] = true;
    //--------------------------------------------------------------------
    // Traverse each edge of p1
    ++rigidity;
    int p3;
    for ( int i = 0; i < static_cast<int>( m_pviListVertexRing1[p2].size() ); ++i ) {
        p3 = m_pviListVertexRing1[p2][i];
        if ( !isMark[ p3 ] ) {
            if ( !CheckASpringExistance( p1, p3 ) ) {
                CreateASpring( Ks, Kd, p1, p3, rigidity );
            }
        }
    }
    // This makes it become Breadth First Search
    for ( int i = 0; i < static_cast<int>( m_pviListVertexRing1[p2].size() ); ++i ) {
        p3 = m_pviListVertexRing1[p2][i];
        BFSCreateSprings( Ks, Kd, rigidity, p1, p3, particleSet, isMark );
    }
    //--------------------------------------------------------------------
}
/*
//-----------------------------------------------------------------------------
template <typename T>
void SpringHalfEdgeModel<T>::DFSCreateSprings ( 
                T Ks, T Kd, int rigidity,       // Ks, Kd, & hop distance
                int p1, int p2,                 // particle#1 & #2
                std::vector<int> & particleSet, // set of connected particles
                bool * isMark )                 // marked particles in particleSet
{
    if ( rigidity >= m_iRigidity ) return;
    //--------------------------------------------------------------------
    // Keep the particle reference in particleSet
    // And mark it as visited
    particleSet.push_back( p2 );
    isMark[ p2 ] = true;
    //--------------------------------------------------------------------
    // Traverse each edge of p1
    ++rigidity;
    int p3;
    for ( int i = 0; i < static_cast<int>( m_pviListVertexRing1[p2].size() ); ++i ) {
        p3 = m_pviListVertexRing1[p2][i];
        if ( !isMark[ p3 ] ) {
            if ( !CheckASpringExistance( p1, p3 ) ) {
                CreateASpring( Ks, Kd, p1, p3, rigidity );
            }
            DFSCreateSprings( Ks, Kd, rigidity, p1, p3, particleSet, isMark );
        }
    }
    //--------------------------------------------------------------------
}
//*/
//-----------------------------------------------------------------------------
//******************************************************************************


//=============================================================================
// For Simulation
//-----------------------------------------------------------------------------
// Advance Simulation
template <typename T>
void SpringHalfEdgeModel<T>::AdvanceSimulation ( T tCurrent, T tNext )
{
    //-----------------------------------------------------
    // Swap pointers
    Simulation::VectorSet<T> * tmpPtr = x0;
    x0   = xEnd;
    xEnd = tmpPtr;
    //------------------------------------------------------
    // ODE Solver
    m_prODESolver->Solve( *x0, *xEnd, tCurrent, tNext, DxDt, this );
    //------------------------------------------------------
    // Copy d/dt X(tNext) into state variables xEnd
    ArrayToState( &(*xEnd)[0] );
    //------------------------------------------------------
    //SimGravity( tCurrent, tNext );
    //------------------------------------------------------
    GetBVHTree()->Update();

    CalAndSetNormals();
}
//-----------------------------------------------------------------------------
// Advance Simulation
template <typename T>
void SpringHalfEdgeModel<T>::AdvanceSimulation ( Simulation::SimClock<T> & simClock )
{
    std::cout << "SpringHalfEdgeModel<T>::AdvanceSimulation ( Simulation::SimClock<T> & simClock ) -- NOT IMPLEMENTED YET!!!\n";
}
//-----------------------------------------------------------------------------
// DxDt
template <typename T>
bool SpringHalfEdgeModel<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
    SpringHalfEdgeModel<T> *pThis = static_cast<SpringHalfEdgeModel<T> *>( userData );
    assert( pThis );
    pThis->ClearAllForces();
    //------------------------------------------------------
    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 SpringHalfEdgeModel<T>::StateToArray ( T *dst )
{
    for ( int i = 0; i < m_iStateSize; ++i ) {
        *(dst++) = m_vpParticleRef[i]->GetPosition()[0];
        *(dst++) = m_vpParticleRef[i]->GetPosition()[1];
        *(dst++) = m_vpParticleRef[i]->GetPosition()[2];
        *(dst++) = m_vpParticleRef[i]->GetVelocity()[0];
        *(dst++) = m_vpParticleRef[i]->GetVelocity()[1];
        *(dst++) = m_vpParticleRef[i]->GetVelocity()[2];
    }
}
//-----------------------------------------------------------------------------
// ArrayToState
template <typename T>
void SpringHalfEdgeModel<T>::ArrayToState ( T *src )
{
    for ( int i = 0; i < m_iStateSize; ++i ) {
        m_vpParticleRef[i]->SetPosition( *src, *(src+1), *(src+2) );
        src += 3;
        m_vpParticleRef[i]->SetVelocity( *src, *(src+1), *(src+2) );
        src += 3;
    }
}
//-----------------------------------------------------------------------------
// DdtStateToArray
template <typename T>
void SpringHalfEdgeModel<T>::DdtStateToArray ( T *xdot, T dt )
{
    for ( int i = 0; i < m_iStateSize; ++i ) {
        if ( m_vpParticleRef[i]->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_vpParticleRef[i]->GetVelocity()[0];
            *(xdot++) = m_vpParticleRef[i]->GetVelocity()[1];
            *(xdot++) = m_vpParticleRef[i]->GetVelocity()[2];
            *(xdot++) = m_vpParticleRef[i]->GetForce()[0];
            *(xdot++) = m_vpParticleRef[i]->GetForce()[1];
            *(xdot++) = m_vpParticleRef[i]->GetForce()[2];
        }
    }
}
//-----------------------------------------------------------------------------
// ClearForces
template <typename T>
void SpringHalfEdgeModel<T>::ClearAllForces ()
{
    for ( int i = 0; i < static_cast<int>( m_vpParticleRef.size() ); ++i ) {
        m_vpParticleRef[i]->SetForce( 0, 0, 0 );
    }
}
//-----------------------------------------------------------------------------
//=============================================================================


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