TAPsFEMMeshTetrahedra.cpp

Go to the documentation of this file.

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

//-----------------------------------------------------------------------------
template <typename T>
void MeshTetrahedra<T>::Init ()
{
    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,1) );
    m_vUndeformedMeshNodes.push_back( Vector3<T>(0,1,0) );
    m_vUndeformedMeshNodes.push_back( Vector3<T>(1,0,0) );
    m_vUndeformedMeshNodes.push_back( Vector3<T>(0,0,0) );

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

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


//=============================================================================
// Operations
//-----------------------------------------------------------------------------
//template <typename T>
//void MeshTetrahedra<T>::ResetToDeformedState ()
//{
//  std::vector< Tetrahedron<T> >::iterator it = m_vElements.begin();
//  while ( it != m_vElements.end() ) {
//      it->ResetToDeformedState();
//      ++it;
//  }
//} 
//-----------------------------------------------------------------------------
template <typename T>
void MeshTetrahedra<T>::ResetToUndeformedState ()
{
    std::vector< Tetrahedron<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 MeshTetrahedra<T>::AssembleStiffnessMatrix ()
{
    #ifdef  TAPs_USE_SYMMETRIC_MATRIX
        Mstiff.SetDirectAccessors();
        std::vector< Tetrahedron<T> >::iterator it = m_vElements.begin();
        while( it != m_vElements.end() ) {
            for ( int n = 0; n < 4; ++n ) {
                for ( int m = 0; m < 4; ++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< Tetrahedron<T> >::iterator it = m_vElements.begin();
        while( it != m_vElements.end() ) {
            for ( int n = 0; n < 4; ++n ) {
                for ( int m = 0; m < 4; ++m ) {
                    Mstiff.DirectAccess( it->ids[n] )->UpdateAdd( it->ids[m], it->GetStiffnessMatrix( n, m ) );
                }
            }
            ++it;
        }
    #endif//TAPs_USE_SYMMETRIC_MATRIX
    //*/

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




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

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

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




//=============================================================================
// Simulation
//-----------------------------------------------------------------------------
template <typename T>
void MeshTetrahedra<T>::AdvanceSimulationStatic ()
{
    /*
    // TEST DATA -- CHECK PASSED
    // For one tetrahedron
    // these forces
    m_Forces[0].SetXYZ(  90.1,  99.2,  40.0 );
    m_Forces[1].SetXYZ( -2.9,   19.4,  17.2 );
    m_Forces[2].SetXYZ( -29.4, -50.4, -12.0 );
    m_Forces[3].SetXYZ( -57.8, -68.2, -45.2 );
    // results in the displacements of node 0 to 3 as follows:
    //  Vector<double,3>( 0.354976, 0.0602062, 0.100089 )'
    //  Vector<double,3>( -0.0882922, -0.000343599, 0.0135531 )'
    //  Vector<double,3>( -0.0649275, 0.0043777, -0.0704632 )'
    //  Vector<double,3>( -0.201756, -0.0642403, -0.0431789 )'
    //*/

    /*
    // TEST DATA -- CHECK PASSED
    // For one tetrahedron, these forces
    m_Forces[0].SetXYZ(   6,  104,  30 );
    m_Forces[1].SetXYZ( -18,   44,  18 );
    m_Forces[2].SetXYZ(  12,  -72, -12 );
    m_Forces[3].SetXYZ(   0,  -76, -36 );
    // with the filters set to
    m_Filter[0].SetXYZ(0,0,0);
    m_Filter[1].SetXYZ(0,1,0);
    m_Filter[2].SetXYZ(0,0,0);
    m_Filter[3].SetXYZ(0,0,0);
    // results in the displacements of node 0 to 3 as follows:
    //  Vector<double,3>( 0, 0, 0 );
    //  Vector<double,3>( 0, 1, 0 );
    //  Vector<double,3>( 0, 0, 0 );
    //  Vector<double,3>( 0, 0, 0 );
    //*/

    /*
    // TEST -- Whatever the value of f_1 is the displacement will be very high (the system is blown up)
    T f_1 = 5.0f;
    m_Forces[0].SetXYZ( f_1, f_1, f_1 );
    m_Forces[1].SetXYZ( f_1, f_1, f_1 );
    m_Forces[2].SetXYZ( f_1, f_1, f_1 );
    m_Forces[3].SetXYZ( f_1, f_1, f_1 );
    //*/

    Mesh<T>::AdvanceSimulationStatic();

    /*
    // TEST
    static int limit = 1;
    for ( int i = 0; i < limit; ++i ) {
        std::vector< Vector3<T> > & defNodes = RetListOfDeformedMeshNodes();
        T o = T(-0.00001);
        defNodes[0] += Vector3<T>( o, o, o );
    }
    //*/
}

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


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

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


//=============================================================================
// Checking
//-----------------------------------------------------------------------------
template <typename T>
std::string MeshTetrahedra<T>::CheckVolumeOfEachElement ()
{
    unsigned int count = 0;
    std::stringstream ss;
    std::vector< Tetrahedron<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 MeshTetrahedra<T>::ApplyPointLoad ( PointForce<T> const * const ptrPtForce )
{
    // The force stored in or referred by the PointFoce is 
    // in FEM local coordinates, not in world coordinates.
    // Hence, do need for coordinate conversion.

    // Add distributed node forces
    if ( !IsFixedNode( m_vElements[ ptrPtForce->GetID() ].GetGlobalNodeNumberFromLocalNode( 0 ) ) ) {
        // RefToPosition is barycentric coordinates
        T weight = T(1) - ptrPtForce->RefToPosition()[0] - ptrPtForce->RefToPosition()[1] - ptrPtForce->RefToPosition()[2];
        m_Forces[ m_vElements[ ptrPtForce->GetID() ].GetGlobalNodeNumberFromLocalNode( 0 ) ] += weight * ptrPtForce->RefToForce();
    }
    if ( !IsFixedNode( m_vElements[ ptrPtForce->GetID() ].GetGlobalNodeNumberFromLocalNode( 1 ) ) ) {
        m_Forces[ m_vElements[ ptrPtForce->GetID() ].GetGlobalNodeNumberFromLocalNode( 1 ) ] += ptrPtForce->RefToPosition()[0] * ptrPtForce->RefToForce();
    }
    if ( !IsFixedNode( m_vElements[ ptrPtForce->GetID() ].GetGlobalNodeNumberFromLocalNode( 2 ) ) ) {
        m_Forces[ m_vElements[ ptrPtForce->GetID() ].GetGlobalNodeNumberFromLocalNode( 2 ) ] += ptrPtForce->RefToPosition()[1] * ptrPtForce->RefToForce();
    }
    if ( !IsFixedNode( m_vElements[ ptrPtForce->GetID() ].GetGlobalNodeNumberFromLocalNode( 3 ) ) ) {
        m_Forces[ m_vElements[ ptrPtForce->GetID() ].GetGlobalNodeNumberFromLocalNode( 3 ) ] += ptrPtForce->RefToPosition()[2] * ptrPtForce->RefToForce();
    }
}
//-----------------------------------------------------------------------------
#endif//TAPs_ADD_INTERACTIONS
//=============================================================================




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

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

    if ( m_pTrx ) {
        glPopMatrix();
    }
}

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

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

    if ( m_pTrx ) {
        glPopMatrix();
    }
}

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

    glPushAttrib( GL_ALL_ATTRIB_BITS );
    glDisable( GL_LIGHTING );
    std::vector< Tetrahedron<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 MeshTetrahedra<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() ].ComputeUndeformedPositionBasedOnCoordinates( 
            0, // dummy value
            (**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 );
    //glTranslatef( 0, 0, 3 );
    glBegin( GL_LINES );
    while ( pf != m_ForceLists.CDForces.end() ) {
        Vector3<T> pos = m_vElements[ (**pf).GetID() ].ComputeUndeformedPositionBasedOnCoordinates( 
            0, // dummy value
            (**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 << (**pf).RefToPosition() << "\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----+----