TAPsModelDefBasedOnFEMHex.cpp

Go to the documentation of this file.

/******************************************************************************
TAPsModelDefBasedOnFEMHex.cpp
******************************************************************************/
/******************************************************************************
SUKITTI PUNAK   (02/26/2010)
UPDATE          (12/05/2010)
******************************************************************************/
#include "TAPsModelDefBasedOnFEMHex.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, typename DATA>
ModelDefBasedOnFEMHex<T,DATA>::ModelDefBasedOnFEMHex ( std::string const & modelFile )
    : ModelDefBasedOnFEM<T,DATA>( modelFile )
    , m_ElementsGridLocs( NULL )
{
    m_NumOfElementsPerGridSize[0] = 
    m_NumOfElementsPerGridSize[1] = 
    m_NumOfElementsPerGridSize[2] = 0;
}
//-----------------------------------------------------------------------------
//template <typename T, typename DATA>
//ModelDefBasedOnFEMHex<T,DATA>::ModelDefBasedOnFEMHex ( ModelDefBasedOnFEMHex<T,DATA> const &orig )
//  :
//{}
//-----------------------------------------------------------------------------
template <typename T, typename DATA>
ModelDefBasedOnFEMHex<T,DATA>::~ModelDefBasedOnFEMHex ()
{
    DeallocateElementsGridLocations();
}
//-----------------------------------------------------------------------------
template <typename T, typename DATA>
std::string ModelDefBasedOnFEMHex<T,DATA>::StrInfo () const
{
    std::ostringstream ss;
    ss << "ModelDefBasedOnFEMHex<" << typeid(T).name() << ">";
    ss << "\n";
    return ss.str();
}
//-----------------------------------------------------------------------------
template <typename T, typename DATA>
void ModelDefBasedOnFEMHex<T,DATA>::AllocateElementsGridLocations ( unsigned int xWidth, unsigned int yHeight, unsigned int zDepth )
{
    m_NumOfElementsPerGridSize[0] = xWidth;
    m_NumOfElementsPerGridSize[1] = yHeight;
    m_NumOfElementsPerGridSize[2] = zDepth;
    if ( m_ElementsGridLocs )   DeallocateElementsGridLocations();
    try {
        m_ElementsGridLocs = new int**[m_NumOfElementsPerGridSize[0]];
        for ( unsigned int i = 0; i < m_NumOfElementsPerGridSize[0]; ++i ) {
            m_ElementsGridLocs[i] = new int*[m_NumOfElementsPerGridSize[1]];
            for ( unsigned int j = 0; j < m_NumOfElementsPerGridSize[1]; ++j ) {
                m_ElementsGridLocs[i][j] = new int[m_NumOfElementsPerGridSize[2]];
                for ( unsigned int k = 0; k < m_NumOfElementsPerGridSize[2]; ++k ) {
                    m_ElementsGridLocs[i][j][k] = -1;
                }
            }
        }
    }
    catch ( Exception & e ) {
        throw e;
    }
}
//-----------------------------------------------------------------------------
template <typename T, typename DATA>
void ModelDefBasedOnFEMHex<T,DATA>::DeallocateElementsGridLocations ()
{
    if ( m_ElementsGridLocs ) {
        for ( unsigned int i = 0; i < m_NumOfElementsPerGridSize[0]; ++i ) {
            for ( unsigned int j = 0; j < m_NumOfElementsPerGridSize[1]; ++j ) {
                delete [] m_ElementsGridLocs[i][j];
            }
            delete [] m_ElementsGridLocs[i];
        }
    }
    delete [] m_ElementsGridLocs;
    m_ElementsGridLocs = NULL;
}
//-----------------------------------------------------------------------------
template <typename T, typename DATA>
int ModelDefBasedOnFEMHex<T,DATA>::GetElementIDFromGridLocation_wErrChk ( unsigned int xWidth, unsigned int yHeight, unsigned int zDepth ) const
{
    if ( m_ElementsGridLocs ) {
        if (    xWidth  < m_NumOfElementsPerGridSize[0]
            &&  yHeight < m_NumOfElementsPerGridSize[1]
            &&  zDepth  < m_NumOfElementsPerGridSize[2] )
        {
            return m_ElementsGridLocs[xWidth][yHeight][zDepth];
        }
        else {
            return -2;  // one of grid is out of range
        }
    }
    else {
        return -10; // m_ElementsGridLocs is not available
    }
}
//-----------------------------------------------------------------------------
template <typename T, typename DATA>
int ModelDefBasedOnFEMHex<T,DATA>::GetElementIDFromGridLocation ( unsigned int xWidth, unsigned int yHeight, unsigned int zDepth ) const
{
    return m_ElementsGridLocs[xWidth][yHeight][zDepth];
}
//-----------------------------------------------------------------------------
template <typename T, typename DATA>
void ModelDefBasedOnFEMHex<T,DATA>::ListAllElementIDsInGridLocationFormat_zYx () const
{
    std::cout << "START: ListAllElementIDsInGridLocationFormat_zYx "
        << "[0-"<<m_NumOfElementsPerGridSize[2]<<"] -- back to front (z-direction)"
        << "["<<m_NumOfElementsPerGridSize[1]<<"-0] -- bottom to top (y-direction)"
        << "[0-"<<m_NumOfElementsPerGridSize[0]<<"] -- left to right (x-direction)"
        << "\n";
    //int limit = 0;
    for ( int z = 0; z < static_cast<int>( m_NumOfElementsPerGridSize[2] ); ++z ) {
        std::cout << "Depth: " << z << "\n";
        for ( int y = static_cast<int>( m_NumOfElementsPerGridSize[1] )-1; y >= 0; --y ) {
            for ( int x = 0; x < static_cast<int>( m_NumOfElementsPerGridSize[0] ); ++x ) {
                std::cout << std::setw(6) << m_ElementsGridLocs[x][y][z];
            }
            std::cout << "\n";

            //++limit;
            //if ( limit > 1000 ) return;
        }
    }
    std::cout << "END: ListAllElementIDsInGridLocationFormat_zYx "
        << "[0-"<<m_NumOfElementsPerGridSize[2]<<"] -- back to front (z-direction)"
        << "["<<m_NumOfElementsPerGridSize[1]<<"-0] -- bottom to top (y-direction)"
        << "[0-"<<m_NumOfElementsPerGridSize[0]<<"] -- left to right (x-direction)"
        << "\n";
}
//-----------------------------------------------------------------------------
//=============================================================================
// Assignment Operator
//-----------------------------------------------------------------------------
//template <typename T, typename DATA>
//inline ModelDefBasedOnFEMHex<T,DATA> & ModelDefBasedOnFEMHex<T,DATA>::operator= ( ModelDefBasedOnFEMHex<T,DATA> const &orig )
//{ 
//    return *this;
//}
//-----------------------------------------------------------------------------
//=============================================================================









//-----------------------------------------------------------------------------
template <typename T, typename DATA>
bool ModelDefBasedOnFEMHex<T,DATA>::GenerateGridsWithDimension (
    T iGridWidth,       
    T iGridHeight,      
    T iGridDepth,       
    T   mass,           
    T   YoungsModulus,  
    T   PoissonsRatio,  
    bool bGenElementsInsideVolume,      
    bool bRecordsElementsGridLocations  
)
{
    std::cout << __FILE__ << " (LINE#" << __LINE__ << ") NOT IMPLEMENTED YET!\n";
    return false;
}


//-----------------------------------------------------------------------------
template <typename T, typename DATA>
bool ModelDefBasedOnFEMHex<T,DATA>::GenerateGridsWithSize (
    int iNumOfXSlices,  
    int iNumOfYSlices,  
    int iNumOfZSlices,  
    T   mass,           
    T   YoungsModulus,  
    T   PoissonsRatio,  
    bool bGenElementsInsideVolume,      
    bool bRecordsElementsGridLocations  
)
{
    if ( !m_pHEModel ) {
        return false;
    }

    // Generate grids
    //GridGenerator<T,DATA> GridGen;
    m_GridGen.SetModelUseToGenGrid( m_pHEModel );
    if ( m_GridGen.GenerateGridsWithSize( iNumOfXSlices, iNumOfYSlices, iNumOfZSlices ) ) {
        // Generate FEM mesh from the generated grids
        if ( GenerateMeshFromGridGenerator( m_GridGen, mass, YoungsModulus, PoissonsRatio, bGenElementsInsideVolume ) ) {
            // Clear the generated grids
            m_GridGen.DeleteGrids();

            // Set transformation support
            m_FEMMesh.SetTransformationSupportPointer( &(m_pHEModel->GetTransform()) );

            return true;
        }
        else {
            std::cout << "ModelDefBasedOnFEMHex<T,DATA>::GenerateGridsWithSize(...) FAILED to generate FEM mesh!\n";
            // Clear the generated grids
            m_GridGen.DeleteGrids();
            return false;
        }
    }
    else {
        std::cout << "ModelDefBasedOnFEMHex<T,DATA>::GenerateGridsWithSize(...) FAILED to generate grids!\n";
        return false;
    }
}


//-----------------------------------------------------------------------------
template <typename T, typename DATA>
bool ModelDefBasedOnFEMHex<T,DATA>::GenerateMeshFromGridGenerator (
    GridGenerator<T,DATA> & GridGen,    
    T   mass,                           
    T   YoungsModulus,                  
    T   PoissonsRatio,                  
    bool bGenElementsInsideVolume,      
    bool bRecordsElementsGridLocations  
)
{
    // FOR HALF-EDGE MODEL ONLY
    if ( !m_pHEModel ) {
        return false;
    }
    if ( !GridGen.IsGridGenerated() ) {
        return false;
    }

    GridGen.RecordHexBarycentricCoordinates();

    // Clear data
    // other data are clear by MeshHexahedraGen function
    m_BaryCoords.clear();

    // Use the static function of Meshexahedra Genetator to generate the hexahedral mesh from the grid generator
    if ( bRecordsElementsGridLocations ) {
        //std::cout << "ModelDefBasedOnFEMHex<T,DATA>::GenerateMeshFromGridGenerator() called with bRecordsElementsGridLocations is true.\n";
        AllocateElementsGridLocations( GridGen.GetGridSizeX()-1, GridGen.GetGridSizeY()-1, GridGen.GetGridSizeZ()-1 );
        if ( !FEM::MeshHexahedraGen<T,DATA>::MeshFromGridGenerator(
                m_FEMMesh, m_BaryCoords, GridGen, mass, YoungsModulus, PoissonsRatio, 
                bGenElementsInsideVolume, m_ElementsGridLocs ) )
        {
            DeallocateElementsGridLocations();
            return false;
        }
        ListAllElementIDsInGridLocationFormat_zYx();    // CHECKING
        //std::cout << "m_NumOfElementsPerGridSize: " << Vector3<unsigned int>( m_NumOfElementsPerGridSize ) << "\n";
    }
    else {
        //std::cout << "ModelDefBasedOnFEMHex<T,DATA>::GenerateMeshFromGridGenerator() called with bRecordsElementsGridLocations is false.\n";
        if ( !FEM::MeshHexahedraGen<T,DATA>::MeshFromGridGenerator(
                m_FEMMesh, m_BaryCoords, GridGen, mass, YoungsModulus, PoissonsRatio, 
                bGenElementsInsideVolume, NULL ) )
        {
            return false;
        }
    }

    // Set the link between HEVertex and FEM element
    // HEVertex'ID <--> FEM element location in the element list
    SetHEVertexIDToElementLocation();

    return true;
}

//-----------------------------------------------------------------------------
template <typename T, typename DATA>
void ModelDefBasedOnFEMHex<T,DATA>::AdvanceSimulationStatic ()
{
    CDUnit().MoveAllGrabbedParts();
    m_FEMMesh.AdvanceSimulationStatic();
    UpdateSurfaceMesh();
}

//-----------------------------------------------------------------------------
template <typename T, typename DATA>
void ModelDefBasedOnFEMHex<T,DATA>::AdvanceSimulationStatic_woFixedDisplacementNodes ()
{
    CDUnit().MoveAllGrabbedParts();
    m_FEMMesh.AdvanceSimulationStatic_woFixedDisplacementNodes();
    UpdateSurfaceMesh();
}

//-----------------------------------------------------------------------------
#ifdef  TAPs_SIM_DYNAMIC_SYSTEM
template <typename T, typename DATA>
void ModelDefBasedOnFEMHex<T,DATA>::AdvanceSimulationDynamic ()
{
    CDUnit().MoveAllGrabbedParts();
    m_FEMMesh.AdvanceSimulationDynamic();
    UpdateSurfaceMesh();
}
#endif//TAPs_SIM_DYNAMIC_SYSTEM

//-----------------------------------------------------------------------------
#ifdef  TAPs_SIM_DYNAMIC_SYSTEM
template <typename T, typename DATA>
void ModelDefBasedOnFEMHex<T,DATA>::AdvanceSimulationDynamic_woFixedDisplacementNodes ()
{
    CDUnit().MoveAllGrabbedParts();
    m_FEMMesh.AdvanceSimulationDynamic_woFixedDisplacementNodes();
    UpdateSurfaceMesh();
}
#endif//TAPs_SIM_DYNAMIC_SYSTEM

//-----------------------------------------------------------------------------
template <typename T, typename DATA>
void ModelDefBasedOnFEMHex<T,DATA>::UpdateSurfaceMesh ()
{
    // Iterate all vertices of the surface mesh
    //std::vector< HexBarycentricCoordsForVertexRef<T,DATA> >::iterator it = m_BaryCoords.begin();
    //std::vector< FEM::Hexahedron<T> > & E = m_FEMMesh.RetListOfElements();
    //std::vector< Vector3<T> > &   Nd = m_FEMMesh.RetListOfDeformedMeshNodes();
    //std::vector< Vector3<T> > &   Nu = m_FEMMesh.RetListOfUndeformedMeshNodes();
    //unsigned int e;
    //int C_1 = m_GridGen.GetGridSizeX()-1;
    //int R_1 = m_GridGen.GetGridSizeY()-1;
    for ( unsigned int v = 0; v < m_BaryCoords.size(); ++v ) {
        //e = m_ElementID[v];
        //e = m_BaryCoords[v].GetElementID();
        m_BaryCoords[v].ComputeNewVertexLocation();
        /*
        m_BaryCoords[v].ComputeNewVertexLocation(
            E[e].DeformedNode(0)->GetX(), E[e].DeformedNode(0)->GetY(), E[e].DeformedNode(0)->GetZ(),
            E[e].DeformedNode(1)->GetX(), E[e].DeformedNode(1)->GetY(), E[e].DeformedNode(1)->GetZ(),
            E[e].DeformedNode(2)->GetX(), E[e].DeformedNode(2)->GetY(), E[e].DeformedNode(2)->GetZ(),
            E[e].DeformedNode(3)->GetX(), E[e].DeformedNode(3)->GetY(), E[e].DeformedNode(3)->GetZ(),
            E[e].DeformedNode(4)->GetX(), E[e].DeformedNode(4)->GetY(), E[e].DeformedNode(4)->GetZ(),
            E[e].DeformedNode(5)->GetX(), E[e].DeformedNode(5)->GetY(), E[e].DeformedNode(5)->GetZ(),
            E[e].DeformedNode(6)->GetX(), E[e].DeformedNode(6)->GetY(), E[e].DeformedNode(6)->GetZ(),
            E[e].DeformedNode(7)->GetX(), E[e].DeformedNode(7)->GetY(), E[e].DeformedNode(7)->GetZ()
        );
        //*/
    }

    // Update all extra vertices
    HexBarycentricCoordsForVertexRef<T,DATA> * pHexBC;
    for ( unsigned int v = 0; v < m_SetOfExtraVertexCoords.size(); ++v ) {
        pHexBC = dynamic_cast< HexBarycentricCoordsForVertexRef<T,DATA>* >( m_SetOfExtraVertexCoords[v] );
        pHexBC->ComputeNewVertexLocation();
    }

    // Update BVH tree
    // The update of normals is not done by this function.
    // The update of normals must be done after this but before rendering.
    ModelDefBasedOnFEM<T,DATA>::UpdateSurfaceMesh();

    // Adjust the surface mesh, i.e. the vertices that is collided with another object.
    //CDUnit().AdjustSurfaceMeshByCDObjs();
    CDUnit().AdjustSurfaceMeshDueToGrabbedParts();

    if ( m_CD.GetCollisionStatus() ) {
        m_CD.SetCollisionStatus( false );
    }
}

//-----------------------------------------------------------------------------
template <typename T, typename DATA>
void ModelDefBasedOnFEMHex<T,DATA>::SetHEVertexIDToElementLocation ()
{
    if ( !m_pHEModel )  return;

    std::list< HEVertex<T> * > pVList;
    // Store all vertices for linking
    HEVertex<T> * V = m_pHEModel->GetVertexList()->Head();
    while ( V ) {
        pVList.push_back( &(*V) );
        V = V->Next();
    }

    //std::cout << __FILE__ << "\n";
    
    // Do the link
    int c = 0;
    for ( unsigned int i = 0; i < m_BaryCoords.size(); ++i ) {
        std::list< HEVertex<T> * >::iterator it = pVList.begin();
        while ( it != pVList.end() ) {
            if ( (void *)(m_BaryCoords[i].RetVertexPtr()) == (void *)((**it).RefToPosition()) ) {
                //(**it).SetID( m_ElementID[i] );
                (**it).SetID( m_BaryCoords[i].GetElementID() );
                pVList.erase( it );

                //std::cout << ++c << ": " << m_ElementID[i] << "\n";

                break;
            }
            ++it;
        }
    }
}

//-----------------------------------------------------------------------------
template <typename T, typename DATA>
bool ModelDefBasedOnFEMHex<T,DATA>::InitCD ()
{
    if ( ModelDefBasedOnFEM<T,DATA>::InitCD() ) {
        m_CD.Init( m_pHEModel, &m_FEMMesh );
        return true;
    }
    return false;
}

//-----------------------------------------------------------------------------
template <typename T, typename DATA>
HEVertex<T> * ModelDefBasedOnFEMHex<T,DATA>::AddExtraVertexInLocalSpace (
    Vector3<T> location,                    
    int elementID,                          
    FEM::Element<T> * elementPtr,           
    Vector3<unsigned int> gridLocations,    
    Vector3<T> coordinates                  
)
{
    FEM::Hexahedron<T> * hexPtr = dynamic_cast< FEM::Hexahedron<T> * >( elementPtr );
    if ( !hexPtr ) {
        std::cout << "Error: elementPtr cannot be converted to Hexahedron<T> *" << std::endl;
        return NULL;
    }

    HEVertex<T> * pV = new HEVertex<T>( location );
    HexBarycentricCoordsForVertexRef<T,DATA> * pBC = 
        new HexBarycentricCoordsForVertexRef<T,DATA>( 
                pV->RefToPosition(), 
                gridLocations[0], gridLocations[1], gridLocations[2],   // grid locations (col, row, & depth) (i.e. the hexahedral element's relative to the FEM mesh
                coordinates[0], coordinates[1], coordinates[2]  // parametric coordinates
            );
    pBC->SetElementID( elementID );
    pBC->SetPtrToHexahedron( hexPtr );
    pV->SetBarycentricPtr( pBC );
    pV->SetID( elementID ); // The id is for identifying the hexahedral element (i.e. the element id)

    m_SetOfExtraVertices.push_back( pV );
    m_SetOfExtraVertexCoords.push_back( pBC );

    return pV;
}

//-----------------------------------------------------------------------------
template <typename T, typename DATA>
bool ModelDefBasedOnFEMHex<T,DATA>::FindElementThatContainsPointInLocalSpace (
    Vector3<T> const & point,               
    unsigned int & elementID,               
    Vector3<unsigned int> & gridLocations,  
    Vector3<T> & coordinates                
)
{
    //std::cout << "START: FindElementThatContainsPointInLocalSpace" << std::endl;

    //std::cout << "GetNumberOfElementsPerGridSize_X: " << GetNumberOfElementsPerGridSize_X() << std::endl;
    //std::cout << "GetNumberOfElementsPerGridSize_Y: " << GetNumberOfElementsPerGridSize_Y() << std::endl;
    //std::cout << "GetNumberOfElementsPerGridSize_Z: " << GetNumberOfElementsPerGridSize_Z() << std::endl;

    bool bFound = false;
    for ( unsigned int x = 0; !bFound && x < GetNumberOfElementsPerGridSize_X(); ++x ) {
        //std::cout << "x: " << x << std::endl;
        for ( unsigned int y = 0; !bFound && y < GetNumberOfElementsPerGridSize_Y(); ++y ) {
            //std::cout << "\ty: " << y << std::endl;
            for ( unsigned int z = 0; !bFound && z < GetNumberOfElementsPerGridSize_Z(); ++z ) {
                //std::cout << "\t\tz: " << z << std::endl;
                int ID = GetElementIDFromGridLocation( x, y, z );

                //std::cout << "elementID: " << ID << std::endl;

                //*
                if ( ID >= 0 ) {
                    elementID = ID;
                    bFound = RefToFEMMesh().RefToElement( elementID ).EstimateParametricCoords( 
                        point[0], point[1], point[2], 
                        coordinates[0], coordinates[1], coordinates[2] ) 
                        > 0 ? true : false;
                    if ( bFound ) {
                        gridLocations[0] = x;
                        gridLocations[1] = y;
                        gridLocations[2] = z;
                        return true;
                    }
                }
                //*/
            }
        }
    }

    //std::cout << "END: FindElementThatContainsPointInLocalSpace" << std::endl;

    return bFound;
}

//-----------------------------------------------------------------------------
template <typename T, typename DATA>
void ModelDefBasedOnFEMHex<T,DATA>::FixOrUnfixNodesByRowColumnDepthRange (
    unsigned int x0,    
    unsigned int x1,    
    unsigned int y0,    
    unsigned int y1,    
    unsigned int z0,    
    unsigned int z1,    
    bool    doFix       
)
{
    unsigned int x, y, z;
    std::vector< HexBarycentricCoordsForVertexRef<T,DATA> >::iterator it = m_BaryCoords.begin();
    while ( it != m_BaryCoords.end() ) {
        it->GetGridLocations( x, y, z );
        if ( x0 <= x && x <= x1 ) {
            if ( y0 <= y && y <= y1 ) {
                if ( z0 <= z && z <= z1 ) {
                    Vector3<T> const * const pNode = it->RetHexahedron()->UndeformedNode( 0 );
                    std::vector< Vector3<T> >::iterator NodeList = m_FEMMesh.RetListOfUndeformedMeshNodes().begin();
                    unsigned int i = 0;
                    while ( NodeList != m_FEMMesh.RetListOfUndeformedMeshNodes().end() ) {
                        if ( pNode == &(*NodeList) ) {
                            m_FEMMesh.FixNode( i, doFix );
                            break;
                        }
                        ++NodeList;
                        ++i;
                    }
                }
            }
        }
        ++it;
    }
}


//=============================================================================
#if defined(__gl_h_) || defined(__GL_H__)
//-----------------------------------------------------------------------------
    template <typename T, typename DATA>
    void ModelDefBasedOnFEMHex<T,DATA>::DrawUndeformedFEMMesh ()
    {
        m_FEMMesh.DrawUndeformedMesh();
    }

    template <typename T, typename DATA>
    void ModelDefBasedOnFEMHex<T,DATA>::DrawDeformedFEMMesh ()
    {
        m_FEMMesh.DrawDeformedMesh();
    }

    template <typename T, typename DATA>
    void ModelDefBasedOnFEMHex<T,DATA>::DrawCollidedMeshParts ()
    {
        glPushMatrix();
        m_pHEModel->GetTransform().TransformByOpenGLForDrawing();
        for ( unsigned int i = 0; i < m_CD.RefToCollidedMeshParts().size(); ++i ) {
            m_CD.RefToCollidedMeshParts()[i].pHEVertex->Draw();
        }
        glPopMatrix();
    }
//-----------------------------------------------------------------------------
#endif // #if defined(__gl_h_) || defined(__GL_H__)
//=============================================================================

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