TAPsFEMMeshHexahedra.cpp

Go to the documentation of this file.

/******************************************************************************
TAPsFEMMeshHexahedra.cpp
******************************************************************************/
/******************************************************************************
SUKITTI PUNAK   (02/11/2010)
UPDATE          (09/03/2010)
******************************************************************************/
#include "TAPsFEMMeshHexahedra.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__FEM
//=============================================================================
// Constructors
//-----------------------------------------------------------------------------
template <typename T>
MeshHexahedra<T>::MeshHexahedra ()
{
    //std::cout << "START: Default Ctor of MeshHexahedra<T> ------------------\n";
    Init();
    //std::cout << "END: Default Ctor of MeshHexahedra<T> --------------------\n";
}
//-----------------------------------------------------------------------------
//template <typename T>
//MeshHexahedra<T>::MeshHexahedra ( MeshHexahedra<T> const &orig )
//{}
//-----------------------------------------------------------------------------
template <typename T>
MeshHexahedra<T>::~MeshHexahedra ()
{}  
//-----------------------------------------------------------------------------
template <typename T>
std::string MeshHexahedra<T>::StrInfo () const
{
    std::stringstream ss;
    std::vector< Hexahedron<T> >::const_iterator it = m_vElements.begin();
    ss  << "FEM::Hexahedra Mesh<" << typeid(T).name() << ">:\n";
    while ( it != m_vElements.end() ) {
        ss << *it;
        ++it;
    }
    return ss.str();
}
//-----------------------------------------------------------------------------

//-----------------------------------------------------------------------------
template <typename T>
void MeshHexahedra<T>::Init ()
{
    FEM::Hexahedron<T>::UpdateShapeFnsPartialDerivatives_1PtGuassRule();
    FEM::Hexahedron<T>::UpdateShapeFnsPartialDerivatives_2PtGuassRule();

    m_vDeformedMeshNodes.push_back( Vector3<T>(0,0,0) );
    m_vDeformedMeshNodes.push_back( Vector3<T>(0,0,0) );
    m_vDeformedMeshNodes.push_back( Vector3<T>(0,0,0) );
    m_vDeformedMeshNodes.push_back( Vector3<T>(0,0,0) );
    m_vDeformedMeshNodes.push_back( Vector3<T>(0,0,0) );
    m_vDeformedMeshNodes.push_back( Vector3<T>(0,0,0) );
    m_vDeformedMeshNodes.push_back( Vector3<T>(0,0,0) );
    m_vDeformedMeshNodes.push_back( Vector3<T>(0,0,0) );

    m_vUndeformedMeshNodes.push_back( Vector3<T>(0,0,0) );
    m_vUndeformedMeshNodes.push_back( Vector3<T>(1,0,0) );
    m_vUndeformedMeshNodes.push_back( Vector3<T>(1,1,0) );
    m_vUndeformedMeshNodes.push_back( Vector3<T>(0,1,0) );
    m_vUndeformedMeshNodes.push_back( Vector3<T>(0,0,1) );
    m_vUndeformedMeshNodes.push_back( Vector3<T>(1,0,1) );
    m_vUndeformedMeshNodes.push_back( Vector3<T>(1,1,1) );
    m_vUndeformedMeshNodes.push_back( Vector3<T>(0,1,1) );

    m_vElements.push_back( Hexahedron<T>( 
        &(RetListOfDeformedMeshNodes()),
        &(RetListOfUndeformedMeshNodes()),
        1.0,    // Young's modulus
        0.01,   // Poisson's ratio
        0, 1, 2, 3, 4, 5, 6, 7 )
    );
}
//-----------------------------------------------------------------------------

//=============================================================================
// Assignment Operator
//-----------------------------------------------------------------------------
//template <typename T>
//MeshHexahedra<T> & MeshHexahedra<T>::operator= ( MeshHexahedra<T> const &orig )
//{}
//-----------------------------------------------------------------------------
//=============================================================================


//=============================================================================
// Operations
//-----------------------------------------------------------------------------
//template <typename T>
//void MeshHexahedra<T>::ResetToDeformedState ()
//{
//  std::vector< Hexahedron<T> >::iterator it = m_vElements.begin();
//  while ( it != m_vElements.end() ) {
//      it->ResetToDeformedState();
//      ++it;
//  }
//} 
//-----------------------------------------------------------------------------
template <typename T>
void MeshHexahedra<T>::ResetToUndeformedState ()
{
    std::vector< Hexahedron<T> >::iterator it = m_vElements.begin();
    while ( it != m_vElements.end() ) {
        it->ResetToUndeformedState();
        ++it;
    }

    #ifdef  TAPs_SIM_DYNAMIC_SYSTEM
        ResetCurrentVelocities();
    #endif//TAPs_SIM_DYNAMIC_SYSTEM
}
//-----------------------------------------------------------------------------
template <typename T>
void MeshHexahedra<T>::AssembleStiffnessMatrix ()
{
    #ifdef  TAPs_USE_SYMMETRIC_MATRIX
        Mstiff.SetDirectAccessors();
        std::vector< Hexahedron<T> >::iterator it = m_vElements.begin();
        while( it != m_vElements.end() ) {
            for ( int n = 0; n < 8; ++n ) {
                for ( int m = 0; m < 8; ++m ) {
                    if ( it->ids[n] <= it->ids[m] ) {
                        Mstiff.DirectAccess( it->ids[n] )->UpdateAdd( it->ids[m], it->GetStiffnessMatrix( n, m ) );
                    }
                }
            }
            ++it;
        }
    #else //TAPs_USE_SYMMETRIC_MATRIX
        Mstiff.SetDirectAccessors();
        std::vector< Hexahedron<T> >::iterator it = m_vElements.begin();
        while( it != m_vElements.end() ) {
            for ( int n = 0; n < 8; ++n ) {
                for ( int m = 0; m < 8; ++m ) {
                    Mstiff.DirectAccess( it->ids[n] )->UpdateAdd( it->ids[m], it->GetStiffnessMatrix( n, m ) );
                }
            }
            ++it;
        }
    #endif//TAPs_USE_SYMMETRIC_MATRIX
    //*/

    // DEBUG
    //std::cout << "HexFEM Mstiff: " << Mstiff << "\n";
}
//-----------------------------------------------------------------------------
template <typename T>
void MeshHexahedra<T>::UpdateStrainDisplacementMatrix ()
{
    std::vector< Hexahedron<T> >::iterator it = m_vElements.begin();
    while ( it != m_vElements.end() ) {
        it->UpdateStrainDisplacementMatrixElements();
        ++it;
    }
}
//-----------------------------------------------------------------------------
//=============================================================================




//=============================================================================
// Get/Set Fns
//-----------------------------------------------------------------------------
template <typename T>
unsigned int MeshHexahedra<T>::GetNumOfElements () const
{
    return m_vElements.size();
}

template <typename T>
Hexahedron<T> const & MeshHexahedra<T>::RefToElement ( unsigned int i ) const
{
    assert( i < GetNumOfElements() );
    return m_vElements[i];
}

template <typename T>
Hexahedron<T> & MeshHexahedra<T>::RefToElement ( unsigned int i )
{
    assert( i < GetNumOfElements() );
    return m_vElements[i];
}
//-----------------------------------------------------------------------------
//=============================================================================




//=============================================================================
// Simulation
//-----------------------------------------------------------------------------
template <typename T>
void MeshHexahedra<T>::AdvanceSimulationStatic ()
{
    Mesh<T>::AdvanceSimulationStatic();
}

#ifdef  TAPs_SIM_DYNAMIC_SYSTEM
template <typename T>
void MeshHexahedra<T>::AdvanceSimulationDynamic ()
{
    Mesh<T>::AdvanceSimulationDynamic();
}
#endif//TAPs_SIM_DYNAMIC_SYSTEM

//-------------------------------------------------------------------
template <typename T>
void MeshHexahedra<T>::AllocateDataForSimulation ( T mass )
{
    Mesh<T>::AllocateDataForSimulation( mass );
}

//-----------------------------------------------------------------------------
//=============================================================================


//=============================================================================
// Checking
//-----------------------------------------------------------------------------
template <typename T>
std::string MeshHexahedra<T>::CheckVolumeOfEachElement ()
{
    unsigned int count = 0;
    std::stringstream ss;
    std::vector< Hexahedron<T> >::iterator it = m_vElements.begin();
    while ( it != m_vElements.end() ) {
        ss //<< "Element# " << count++ << " Deformed - Undeformed Volume: " 
            << "#" << count++ << ": " << it->CalDeformedVolume() << " - " << it->CalUndeformedVolume() << "; ";
        ++it;
    }
    return ss.str();
}   
//-----------------------------------------------------------------------------
//=============================================================================



//=============================================================================
// For interactions
#ifdef  TAPs_ADD_INTERACTIONS
//-----------------------------------------------------------------------------
template <typename T>
void MeshHexahedra<T>::ApplyPointLoad ( PointForce<T> const * const ptrPtForce )
{
    static T weights[8];
    Hexahedron<T>::InterpolateShapeFunctions( ptrPtForce->RefToPosition()[0], ptrPtForce->RefToPosition()[1], ptrPtForce->RefToPosition()[2], weights );

    // Add distributed node forces
    if ( !IsFixedNode( m_vElements[ ptrPtForce->GetID() ].GetGlobalNodeNumberFromLocalNode( 0 ) ) ) {
        m_Forces[ m_vElements[ ptrPtForce->GetID() ].GetGlobalNodeNumberFromLocalNode( 0 ) ] += weights[0] * ptrPtForce->RefToForce();
    }
    if ( !IsFixedNode( m_vElements[ ptrPtForce->GetID() ].GetGlobalNodeNumberFromLocalNode( 1 ) ) ) {
        m_Forces[ m_vElements[ ptrPtForce->GetID() ].GetGlobalNodeNumberFromLocalNode( 1 ) ] += weights[1] * ptrPtForce->RefToForce();
    }
    if ( !IsFixedNode( m_vElements[ ptrPtForce->GetID() ].GetGlobalNodeNumberFromLocalNode( 2 ) ) ) {
        m_Forces[ m_vElements[ ptrPtForce->GetID() ].GetGlobalNodeNumberFromLocalNode( 2 ) ] += weights[2] * ptrPtForce->RefToForce();
    }
    if ( !IsFixedNode( m_vElements[ ptrPtForce->GetID() ].GetGlobalNodeNumberFromLocalNode( 3 ) ) ) {
        m_Forces[ m_vElements[ ptrPtForce->GetID() ].GetGlobalNodeNumberFromLocalNode( 3 ) ] += weights[3] * ptrPtForce->RefToForce();
    }
    if ( !IsFixedNode( m_vElements[ ptrPtForce->GetID() ].GetGlobalNodeNumberFromLocalNode( 4 ) ) ) {
        m_Forces[ m_vElements[ ptrPtForce->GetID() ].GetGlobalNodeNumberFromLocalNode( 4 ) ] += weights[4] * ptrPtForce->RefToForce();
    }
    if ( !IsFixedNode( m_vElements[ ptrPtForce->GetID() ].GetGlobalNodeNumberFromLocalNode( 5 ) ) ) {
        m_Forces[ m_vElements[ ptrPtForce->GetID() ].GetGlobalNodeNumberFromLocalNode( 5 ) ] += weights[5] * ptrPtForce->RefToForce();
    }
    if ( !IsFixedNode( m_vElements[ ptrPtForce->GetID() ].GetGlobalNodeNumberFromLocalNode( 6 ) ) ) {
        m_Forces[ m_vElements[ ptrPtForce->GetID() ].GetGlobalNodeNumberFromLocalNode( 6 ) ] += weights[6] * ptrPtForce->RefToForce();
    }
    if ( !IsFixedNode( m_vElements[ ptrPtForce->GetID() ].GetGlobalNodeNumberFromLocalNode( 7 ) ) ) {
        m_Forces[ m_vElements[ ptrPtForce->GetID() ].GetGlobalNodeNumberFromLocalNode( 7 ) ] += weights[7] * ptrPtForce->RefToForce();
    }

    /*
    // DEBUG
    std::cout << "weights: ";
    for ( int i = 0; i < 8; ++i ) {
        std::cout << "\t" << weights[i];
    }
    std::cout << "\n";
    //*/
}
//-----------------------------------------------------------------------------
#endif//TAPs_ADD_INTERACTIONS
//=============================================================================




//=============================================================================
// OpenGL
#if defined(__gl_h_) || defined(__GL_H__)
//-----------------------------------------------------------------------------
template <typename T>
void MeshHexahedra<T>::Draw ()
{
    if ( m_pTrx ) {
        glPushMatrix();
        m_pTrx->TransformByOpenGLForDrawing();
    }

    glPushAttrib( GL_ALL_ATTRIB_BITS );
    glDisable( GL_LIGHTING );
    std::vector< Hexahedron<T> >::iterator it = m_vElements.begin();
    while ( it != m_vElements.end() ) {
        it->Draw();
        ++it;
    }
    glPopAttrib();

    if ( m_pTrx ) {
        glPopMatrix();
    }
}

//-----------------------------------------------------------------------------
template <typename T>
void MeshHexahedra<T>::DrawUndeformedMesh ()
{
    if ( m_pTrx ) {
        glPushMatrix();
        m_pTrx->TransformByOpenGLForDrawing();
    }

    glPushAttrib( GL_ALL_ATTRIB_BITS );
    glDisable( GL_LIGHTING );
    std::vector< Hexahedron<T> >::iterator it = m_vElements.begin();
    while ( it != m_vElements.end() ) {
        it->DrawUndeformedMesh();
        ++it;
    }
    glPopAttrib();

    if ( m_pTrx ) {
        glPopMatrix();
    }
}

//-----------------------------------------------------------------------------
template <typename T>
void MeshHexahedra<T>::DrawDeformedMesh ()
{
    if ( m_pTrx ) {
        glPushMatrix();
        m_pTrx->TransformByOpenGLForDrawing();
    }

    glPushAttrib( GL_ALL_ATTRIB_BITS );
    glDisable( GL_LIGHTING );
    std::vector< Hexahedron<T> >::iterator it = m_vElements.begin();
    while ( it != m_vElements.end() ) {
        it->DrawDeformedMesh();
        ++it;
    }
    glPopAttrib();

    if ( m_pTrx ) {
        glPopMatrix();
    }
}

//-----------------------------------------------------------------------------
#ifdef  TAPs_ADD_INTERACTIONS
template <typename T>
void MeshHexahedra<T>::DrawInteractionForces ()
{
    if ( m_pTrx ) {
        glPushMatrix();
        m_pTrx->TransformByOpenGLForDrawing();
    }

    // Point forces
    ForceListsPtr<T>::TypeDef_PointForces::const_iterator pf = m_ForceLists.PointForces.begin();
    glPushAttrib( GL_ALL_ATTRIB_BITS );
    glLineWidth( 1 );
    glDisable( GL_LIGHTING );
    glBegin( GL_LINES );
    while ( pf != m_ForceLists.PointForces.end() ) {
        Vector3<T> pos = m_vElements[ (**pf).GetID() ].ComputeDeformedPositionBasedOnCoordinates( 
            (**pf).RefToPosition()[0], (**pf).RefToPosition()[1], (**pf).RefToPosition()[2] );

        //std::cout << *(m_vElements[ (**pf).GetID() ].DeformedNode( 6 )) << "\n";
        //std::cout << "pos: " << pos << "\n";
        //std::cout << "force: " << (**pf).RefToForce() << "\n";

        glColor3f( 0, 1, 0 );
        glVertex3fv( pos.GetDataFloat() );
        glColor3f( 1, 0, 0 );
        glVertex3fv( ((**pf).RefToForce()+pos).GetDataFloat() );
        ++pf;
    }
    glEnd();
    glPopAttrib();

    // CD forces
    pf = m_ForceLists.CDForces.begin();
    glPushAttrib( GL_ALL_ATTRIB_BITS );
    glLineWidth( 1 );
    glDisable( GL_LIGHTING );
    glBegin( GL_LINES );
    while ( pf != m_ForceLists.CDForces.end() ) {
        Vector3<T> pos = m_vElements[ (**pf).GetID() ].ComputeDeformedPositionBasedOnCoordinates( 
            (**pf).RefToPosition()[0], (**pf).RefToPosition()[1], (**pf).RefToPosition()[2] );

        //std::cout << *(m_vElements[ (**pf).GetID() ].DeformedNode( 6 )) << "\n";
        //std::cout << "pos: " << pos << "\n";
        //std::cout << "force: " << (**pf).RefToForce() << "\n";

        glColor3f( 0, 1, 0 );
        glVertex3fv( pos.GetDataFloat() );
        glColor3f( 1, 0, 0 );
        glVertex3fv( ((**pf).RefToForce()+pos).GetDataFloat() );
        ++pf;
    }
    glEnd();
    glPopAttrib();

    if ( m_pTrx ) {
        glPopMatrix();
    }
}
#endif//TAPs_ADD_INTERACTIONS
//-----------------------------------------------------------------------------
#endif // OpenGL
//=============================================================================




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