TAPsFEMMesh.cpp

Go to the documentation of this file.

/******************************************************************************
TAPsFEMMesh.cpp
******************************************************************************/
/******************************************************************************
SUKITTI PUNAK   (12/15/2009)
UPDATE          (10/27/2010)
******************************************************************************/
#include "TAPsFEMMesh.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.

#include "TAPsFEMMesh_withSparseMatrix.cpp"

BEGIN_NAMESPACE_TAPs__FEM
//=============================================================================
// Constructors
//-----------------------------------------------------------------------------
template <typename T>
Mesh<T>::Mesh ()
    : m_iIterationLimit( 20 ), m_tErrorThreshold( 1E-6 )
{}
//-----------------------------------------------------------------------------
//template <typename T>
//Mesh<T>::Mesh ( Mesh<T> const &orig )
//{}
//-----------------------------------------------------------------------------
template <typename T>
Mesh<T>::~Mesh ()
{}  
//-----------------------------------------------------------------------------
//template <typename T>
//std::string Mesh<T>::StrInfo () const
//{}
//-----------------------------------------------------------------------------
//=============================================================================
// Assignment Operator
//-----------------------------------------------------------------------------
//template <typename T>
//Mesh<T> & Mesh<T>::operator= ( Mesh<T> const &orig )
//{}

//=============================================================================
// Operations
//-----------------------------------------------------------------------------
template <typename T>
void Mesh<T>::SetSizeOfStiffnessMatrix ()
{
    Mstiff.ChangeNumOfRows( GetNumOfNodes() );
}
//-----------------------------------------------------------------------------
//=============================================================================




//=============================================================================
// Get/Set
//-----------------------------------------------------------------------------
template <typename T>
unsigned int Mesh<T>::GetNumOfNodes () const
{
    return m_vDeformedMeshNodes.size();
}

template <typename T>
Vector3<T> const & Mesh<T>::ForceAtNode ( unsigned int node ) const
{
    assert( node < GetNumOfNodes() );
    return m_Forces[node];
}

template <typename T>
Vector3<T> & Mesh<T>::ForceAtNode ( unsigned int node )
{
    assert( node < GetNumOfNodes() );
    return m_Forces[node];
}

template <typename T>
Vector3<T> const & Mesh<T>::FilterAtNode ( unsigned int node ) const
{
    assert( node < GetNumOfNodes() );
    return m_Filters[node];
}

template <typename T>
Vector3<T> & Mesh<T>::FilterAtNode ( unsigned int node )
{
    assert( node < GetNumOfNodes() );
    return m_Filters[node];
}

template <typename T>
Vector3<T> const & Mesh<T>::OrigPositionOfNode ( unsigned int node ) const
{
    assert( node < GetNumOfNodes() );
    return m_vUndeformedMeshNodes[node];
}

template <typename T>
Vector3<T> & Mesh<T>::OrigPositionOfNode ( unsigned int node )
{
    assert( node < GetNumOfNodes() );
    return m_vUndeformedMeshNodes[node];
}

template <typename T>
Vector3<T> const & Mesh<T>::CurrPositionOfNode ( unsigned int node ) const
{
    assert( node < GetNumOfNodes() );
    return m_vDeformedMeshNodes[node];
}

template <typename T>
Vector3<T> & Mesh<T>::CurrPositionOfNode ( unsigned int node )
{
    assert( node < GetNumOfNodes() );
    return m_vDeformedMeshNodes[node];
}

template <typename T>
void Mesh<T>::FixNode ( unsigned int node, bool b )
{
    assert( node < GetNumOfNodes() );
    m_FixedNodes[ node ] = b;
    //m_Filters[ node ].Clear();
}

template <typename T>
void Mesh<T>::FixNode ( unsigned int node )
{
    assert( node < GetNumOfNodes() );
    m_FixedNodes[ node ] = true;
    //m_Filters[ node ].Clear();
}

template <typename T>
void Mesh<T>::UnfixNode ( unsigned int node )
{
    assert( node < GetNumOfNodes() );
    m_FixedNodes[ node ] = false;
    //m_Filters[ node ].SetXYZ(1,1,1);
}

template <typename T>
void Mesh<T>::ToggleFixUnfixOfNode ( unsigned int node )
{
    assert( node < GetNumOfNodes() );
    m_FixedNodes[ node ] = !m_FixedNodes[ node ];
}

template <typename T>
bool Mesh<T>::IsFixedNode ( unsigned int node )
{
    assert( node < GetNumOfNodes() );
    return m_FixedNodes[ node ];
}

template <typename T>
void Mesh<T>::SetFixedDisplacementOfNode ( unsigned int node, Vector3<T> const & u )
{
    assert( node < GetNumOfNodes() );
    m_FixedDisplacements[ node ] = u;
}

template <typename T>
void Mesh<T>::SetFixedDisplacementOfNode ( unsigned int node, T x, T y, T z )
{
    assert( node < GetNumOfNodes() );
    m_FixedDisplacements[ node ].SetXYZ( x, y, z );
}


#ifdef  TAPs_SIM_DYNAMIC_SYSTEM
    template <typename T>
    T Mesh<T>::GetTimeStep () const
    {
        return m_tTimeStep;
    }
    template <typename T>
    void Mesh<T>::SetTimeStep ( T timestep )
    {
        m_tTimeStep = timestep;
    }

    template <typename T>
    Vector3<T> const & Mesh<T>::MassOfNode ( unsigned int node ) const
    {
        assert( node < GetNumOfNodes() );
        return m_vNodeMass[node];
    }
    template <typename T>
    Vector3<T> & Mesh<T>::MassOfNode ( unsigned int node )
    {
        assert( node < GetNumOfNodes() );
        return m_vNodeMass[node];
    }

    template <typename T>
    Vector3<T> const & Mesh<T>::CurrentVelocityOfNode ( unsigned int node ) const
    {
        assert( node < GetNumOfNodes() );
        return m_vNodeVelocities_curr[node];
    }
    template <typename T>
    Vector3<T> & Mesh<T>::CurrentVelocityOfNode ( unsigned int node )
    {
        assert( node < GetNumOfNodes() );
        return m_vNodeVelocities_curr[node];
    }

    template <typename T>
    Vector3<T> const & Mesh<T>::NextVelocityOfNode ( unsigned int node ) const
    {
        assert( node < GetNumOfNodes() );
        return m_vNodeVelocities_next[node];
    }
    template <typename T>
    Vector3<T> & Mesh<T>::NextVelocityOfNode ( unsigned int node )
    {
        assert( node < GetNumOfNodes() );
        return m_vNodeVelocities_next[node];
    }

    template <typename T>
    T Mesh<T>::GetDampingMass () const
    {
        return m_tNodeDamping_Alpha;
    }
    template <typename T>
    void Mesh<T>::SetDampingMass ( T d )
    {
        m_tNodeDamping_Alpha = d;
    }

    template <typename T>
    T Mesh<T>::GetDampingStiffness () const
    {
        return m_tNodeDamping_Beta;
    }
    template <typename T>
    void Mesh<T>::SetDampingStiffness ( T d )
    {
        m_tNodeDamping_Beta = d;
    }

    template <typename T>
    void Mesh<T>::ResetCurrentVelocities ()
    {
        for ( unsigned int i = 0; i < GetNumOfNodes(); ++i ) {
            m_vNodeVelocities[m_currVelIdx][i].Clear();
        }
    }

#endif//TAPs_SIM_DYNAMIC_SYSTEM
//-----------------------------------------------------------------------------
//=============================================================================


//-----------------------------------------------------------------------------
template <typename T>
void Mesh<T>::UpdateDisplacementVector ()
{
    for ( unsigned int i = 0; i < GetNumOfNodes(); ++i ) {
        m_u[i] = m_vDeformedMeshNodes[i] - m_vUndeformedMeshNodes[i];
    }
}

//-----------------------------------------------------------------------------
template <typename T>
void Mesh<T>::UpdateDisplacementVector_wFixedDisplacementNodes ()
{
    for ( unsigned int i = 0; i < GetNumOfNodes(); ++i ) {
        // Update with fixed displacements
        if ( m_FixedNodes[i] ) {
            m_u_out[i] = m_u[i] = m_FixedDisplacements[i];
        }
        // Update with displacements from simulations
        else {
            m_u[i] = m_vDeformedMeshNodes[i] - m_vUndeformedMeshNodes[i];
        }
    }
}

//-----------------------------------------------------------------------------
template <typename T>
void Mesh<T>::UpdateDeformedMeshNodes ()
{
    for ( unsigned int i = 0; i < GetNumOfNodes(); ++i ) {
        m_vDeformedMeshNodes[i] = m_vUndeformedMeshNodes[i] + m_u_out[i];
    }
}

//-------------------------------------------------------------------
template <typename T>
void Mesh<T>::AllocateDataForSimulation ( T mass )
{
    unsigned int size = GetNumOfNodes();

    // Set pointer to a transformation support object
    m_pTrx = NULL;

    // Temps for MPCGSolver
    v.resize( size );
    Pinv.resize( size );
    r.resize( size );
    s.resize( size );
    c.resize( size );
    q.resize( size );

    // For prescribed fixed displacements
    m_FixedNodes.resize( size );
    m_FixedDisplacements.resize( size );

    // The data are for SparseMatrix_Matrix3x3<T>::MPCGSolverSpecificForFEM function.
    m_Filters.clear();
    m_Filters.resize( size );
    m_Forces.clear();
    m_Forces.resize( size );
    m_u.clear();
    m_u.resize( size );
    m_u_out.clear();
    m_u_out.resize( size );
    m_Precond.clear();
    m_Precond.resize( size );

    // Initialize
    for ( unsigned int i = 0; i < size; ++i ) {
        m_FixedNodes[i] = false;
        m_Filters[i].SetXYZ(1,1,1);
        m_Precond[i].SetXYZ( 1,1,1 );
        
    }

#ifdef  TAPs_SIM_DYNAMIC_SYSTEM
    m_tTimeStep = T(0.001); // 1 ms
    m_vNodeMass.clear();
    m_vNodeMass.resize( size );
    m_currVelIdx = 0;
    m_nextVelIdx = 1;
    m_vNodeVelocities[m_currVelIdx].clear();
    m_vNodeVelocities[m_currVelIdx].resize( size );
    m_vNodeVelocities[m_nextVelIdx].clear();
    m_vNodeVelocities[m_nextVelIdx].resize( size );

    // Initialize
    T m = mass / size;
    for ( unsigned int i = 0; i < size; ++i ) {
        m_vNodeMass[i].SetXYZ(m,m,m);
    }

    m_tNodeDamping_Alpha = 0.1;
    m_tNodeDamping_Beta  = 0.1;

    // Copy stiffness matrix
    Mstiff_2 = Mstiff;
    m_A = Mstiff;

    m_b.clear();
    m_b.resize( size );
#endif//TAPs_SIM_DYNAMIC_SYSTEM

#ifdef  TAPs_ADD_INTERACTIONS
    m_ForceLists.ClearAllForces();
#endif//TAPs_ADD_INTERACTIONS

}

//-------------------------------------------------------------------
template <typename T>
void Mesh<T>::ClearAllNodeForces ()
{
    for ( unsigned int i = 0; i < GetNumOfNodes(); ++i ) {
        m_Forces[i].Clear();
    }
}

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




//=============================================================================
// For interactions
#ifdef  TAPs_ADD_INTERACTIONS

    //---------------------------------------------------------------
    template <typename T>
    void Mesh<T>::AddPointForce ( PointForce<T> * const ptrPtForce )
    {
        ptrPtForce->RefToForce() = m_pTrx->GetInverseMatrixTransform() * ptrPtForce->RefToForce();
        m_ForceLists.PointForces.push_back( ptrPtForce );
    }
    //---------------------------------------------------------------
    template <typename T>
    void Mesh<T>::AddPointForceInLocalCoordinates ( PointForce<T> * const ptrPtForce )
    {
        m_ForceLists.PointForces.push_back( ptrPtForce );
    }
    //---------------------------------------------------------------
    template <typename T>
    void Mesh<T>::RemovePointForce ( PointForce<T> * const ptrPtForce )
    {
        m_ForceLists.PointForces.remove( ptrPtForce );
    }

    //---------------------------------------------------------------
    template <typename T>
    void Mesh<T>::AddCDForce ( PointForce<T> * const ptrCDForce )
    {
        ptrCDForce->RefToForce() = m_pTrx->GetInverseMatrixTransform() * ptrCDForce->RefToForce();
        m_ForceLists.CDForces.push_back( ptrCDForce );
    }
    //---------------------------------------------------------------
    template <typename T>
    void Mesh<T>::AddCDForceInLocalCoordinates ( PointForce<T> * const ptrCDForce )
    {
        m_ForceLists.CDForces.push_back( ptrCDForce );
    }
    //---------------------------------------------------------------
    template <typename T>
    void Mesh<T>::RemoveCDForce ( PointForce<T> * const ptrCDForce )
    {
        m_ForceLists.CDForces.remove( ptrCDForce );
    }

    //---------------------------------------------------------------
    template <typename T>
    void Mesh<T>::AddConnectionForce ( PointForce<T> * const ptrConnectionForce )
    {
        ptrConnectionForce->RefToForce() = m_pTrx->GetInverseMatrixTransform() * ptrConnectionForce->RefToForce();
        m_ForceLists.ConnectionForces.push_back( ptrConnectionForce );
    }
    //---------------------------------------------------------------
    template <typename T>
    void Mesh<T>::AddConnectionForceInLocalCoordinates ( PointForce<T> * const ptrConnectionForce )
    {
        m_ForceLists.ConnectionForces.push_back( ptrConnectionForce );
    }
    //---------------------------------------------------------------
    template <typename T>
    void Mesh<T>::RemoveConnectionForce ( PointForce<T> * const ptrConnectionForce )
    {
        m_ForceLists.ConnectionForces.remove( ptrConnectionForce );
    }

    //---------------------------------------------------------------
    template <typename T>
    void Mesh<T>::ClearAllPointLoads ()
    {
        m_ForceLists.ClearAllPointForces();
    }
    //---------------------------------------------------------------
    template <typename T>
    void Mesh<T>::ClearAllCDLoads ()
    {
        m_ForceLists.ClearAllCDForces();
    }
    //---------------------------------------------------------------
    template <typename T>
    void Mesh<T>::ClearAllConnectionLoads ()
    {
        m_ForceLists.ClearAllConnectionForces();
    }
    //---------------------------------------------------------------
    template <typename T>
    void Mesh<T>::ClearAllAppliedLoads ()
    {
        m_ForceLists.ClearAllForces();
    }

    //---------------------------------------------------------------
    template <typename T>
    void Mesh<T>::ApplyAllLoads ()
    {
        for ( unsigned int i = 0; i < GetNumOfNodes(); ++i ) {
            m_Forces[i].Clear();
        }
        // Point forces
        ForceListsPtr<T>::TypeDef_PointForces::iterator af = m_ForceLists.PointForces.begin();
        while ( af != m_ForceLists.PointForces.end() ) {
            ApplyPointLoad( *af );
            ++af;
        }
        // CD forces
        af = m_ForceLists.CDForces.begin();
        while ( af != m_ForceLists.CDForces.end() ) {
            ApplyPointLoad( *af );
            ++af;
        }
        // Connection forces
        af = m_ForceLists.ConnectionForces.begin();
        while ( af != m_ForceLists.ConnectionForces.end() ) {
            ApplyPointLoad( *af );
            ++af;
        }
    }
//-----------------------------------------------------------------------------
#endif//TAPs_ADD_INTERACTIONS
//=============================================================================

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

        glPushAttrib( GL_ALL_ATTRIB_BITS );
        glDisable( GL_LIGHTING );
        glPointSize( 5 );
        glColor3f( 0, 1, 0 );
        glBegin( GL_POINTS );
        for ( unsigned int i = 0; i < GetNumOfNodes(); ++i ) {
            if ( m_FixedNodes[i] ) {
                glColor3f( 1, 0, 0 );
                glVertex3f( m_vDeformedMeshNodes[i][0], m_vDeformedMeshNodes[i][1], m_vDeformedMeshNodes[i][2] );
            }
            else if ( m_Filters[i][0] == 0 && m_Filters[i][1] == 0 && m_Filters[i][2] == 0 ) {
                //glPushMatrix();
                //glScalef( 5, 5, 5 );
                //glTranslatef( m_vUndeformedMeshNodes[i][0], m_vUndeformedMeshNodes[i][1], m_vUndeformedMeshNodes[i][2] );
                //OpenGL::OpenGLUsefulObj<T>::DrawWorldCoordinateAxes();
                //glPopMatrix();
                glColor3f( 0, 1, 0 );
                glVertex3f( m_vDeformedMeshNodes[i][0], m_vDeformedMeshNodes[i][1], m_vDeformedMeshNodes[i][2] );
            }
        }
        glEnd();
        glEnd();
        glPopAttrib();

        if ( m_pTrx ) {
            glPopMatrix();
        }
    }
//-----------------------------------------------------------------------------
#endif // OpenGL
//=============================================================================
END_NAMESPACE_TAPs__FEM
//34567890123456789012345678901234567890123456789012345678901234567890123456789
//--+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----