TAPsFEMTetrahedron.cpp

Go to the documentation of this file.

/******************************************************************************
TAPsFEMTetrahedron.cpp
******************************************************************************/
/******************************************************************************
SUKITTI PUNAK   (12/15/2009)
UPDATE          (05/01/2010)
******************************************************************************/
#include "TAPsFEMTetrahedron.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>
Tetrahedron<T>::Tetrahedron ( 
        std::vector< Vector3<T> > * deformedNodes,      
        std::vector< Vector3<T> > * undeformedNodes,    
        T YoungModulus,     
        T PoissionRatio,    
        unsigned int globalNode0,   
        unsigned int globalNode1,   
        unsigned int globalNode2,   
        unsigned int globalNode3    
) : Element( deformedNodes, undeformedNodes, YoungModulus, PoissionRatio )
{
    ids[0] = globalNode0;
    ids[1] = globalNode1;
    ids[2] = globalNode2;
    ids[3] = globalNode3;
    CalUndeformedVolume();  // the result volume is stored in Volume (the property (undeformed) for volume)
    UpdateElasticMatrixElements();  // update the elastic matrix due to Young's modulus, Poisson's ratio, and (undeformed) volume
    UpdateStrainDisplacementMatrixElements();   // update the elements of the strain-displacement matrix
    UpdateStiffnessMatrixElements();    // update all 16 sub stiffness matrices of this element
}
//-----------------------------------------------------------------------------
//template <typename T>
//Tetrahedron<T>::Tetrahedron ( Tetrahedron<T> const &orig )
//{}
//-----------------------------------------------------------------------------
template <typename T>
Tetrahedron<T>::~Tetrahedron ()
{}  
//-----------------------------------------------------------------------------
template <typename T>
std::string Tetrahedron<T>::StrInfo () const
{
    std::stringstream ss;
    ss  << "FEM::Tetrahedron<" << typeid(T).name() << ">:\n"
        << "\tDeformed - Undeformed Node 0: " << (*pX)[ids[0]] << " - " << (*pU)[ids[0]] << "\n"
        << "\tDeformed - Undeformed Node 1: " << (*pX)[ids[1]] << " - " << (*pU)[ids[1]] << "\n"
        << "\tDeformed - Undeformed Node 2: " << (*pX)[ids[2]] << " - " << (*pU)[ids[2]] << "\n"
        << "\tDeformed - Undeformed Node 3: " << (*pX)[ids[3]] << " - " << (*pU)[ids[3]] << "\n"
        << "\tDeformed - Undeformed Volume: " << CalDeformedVolume()  << " - " << Volume << "\n";
    return ss.str();
}
//-----------------------------------------------------------------------------
//=============================================================================
// Assignment Operator
//-----------------------------------------------------------------------------
//template <typename T>
//Tetrahedron<T> & Tetrahedron<T>::operator= ( Tetrahedron<T> const &orig )
//{}
//-----------------------------------------------------------------------------
//=============================================================================
// Get/Set Functions
//-----------------------------------------------------------------------------
template <typename T>
Matrix3x3<T> const & Tetrahedron<T>::GetStiffnessMatrix ( unsigned int n, unsigned int m ) const
{
    assert( 0 <= n && n <= 3 );
    assert( 0 <= m && m <= 3 );
    return Ke[n][m];
}   
//-----------------------------------------------------------------------------
//=============================================================================
// Operations
//-----------------------------------------------------------------------------
template <typename T>
void Tetrahedron<T>::SetNode ( unsigned int i, unsigned int globalNodeNumber )
{
    assert( 0 <= i && i <= 3 );
    ids[i] = globalNodeNumber;
}   
//-----------------------------------------------------------------------------
//template <typename T>
//void Tetrahedron<T>::ResetToDeformedState ()
//{
//} 
//-----------------------------------------------------------------------------
template <typename T>
void Tetrahedron<T>::ResetToUndeformedState ()
{
    (*pU)[ids[0]] = (*pX)[ids[0]];
    (*pU)[ids[1]] = (*pX)[ids[1]];
    (*pU)[ids[2]] = (*pX)[ids[2]];
    (*pU)[ids[3]] = (*pX)[ids[3]];
}   
//-----------------------------------------------------------------------------
template <typename T>
Vector3<T> * const Tetrahedron<T>::DeformedNode ( unsigned int i )
{
    assert( 0 <= i && i <= 3 );
    return &(*pU)[ids[i]];
}   
//-----------------------------------------------------------------------------
template <typename T>
Vector3<T> * const Tetrahedron<T>::UndeformedNode ( unsigned int i )
{
    assert( 0 <= i && i <= 3 );
    return &(*pX)[ids[i]];
}   
//-----------------------------------------------------------------------------
template <typename T>
Vector3<T> Tetrahedron<T>::ComputeDeformedPositionBasedOnCoordinates ( T coord0, T coord1, T coord2, T coord3 )
{
    Vector3<T> & T0  = (*pX)[ids[0]];
    Vector3<T> T10 = (*pX)[ids[1]] - T0;
    Vector3<T> T20 = (*pX)[ids[2]] - T0;
    Vector3<T> T30 = (*pX)[ids[3]] - T0;

    T x = T10[0]*coord1 + T20[0]*coord2 + T30[0]*coord3 + T0[0];
    T y = T10[1]*coord1 + T20[1]*coord2 + T30[1]*coord3 + T0[1];
    T z = T10[2]*coord1 + T20[2]*coord2 + T30[2]*coord3 + T0[2];

    return Vector3<T>( x, y, z );
}
//-----------------------------------------------------------------------------
template <typename T>
Vector3<T> Tetrahedron<T>::ComputeUndeformedPositionBasedOnCoordinates ( T coord0, T coord1, T coord2, T coord3 )
{
    Vector3<T> & T0  = (*pU)[ids[0]];
    Vector3<T> T10 = (*pU)[ids[1]] - T0;
    Vector3<T> T20 = (*pU)[ids[2]] - T0;
    Vector3<T> T30 = (*pU)[ids[3]] - T0;

    T x = T10[0]*coord1 + T20[0]*coord2 + T30[0]*coord3 + T0[0];
    T y = T10[1]*coord1 + T20[1]*coord2 + T30[1]*coord3 + T0[1];
    T z = T10[2]*coord1 + T20[2]*coord2 + T30[2]*coord3 + T0[2];

    return Vector3<T>( x, y, z );
}
//-----------------------------------------------------------------------------
template <typename T>
T   Tetrahedron<T>::CalDeformedVolume () const
{
    return ((((*pU)[ids[1]]-(*pU)[ids[0]]) ^ ((*pU)[ids[2]]-(*pU)[ids[0]])) * ((*pU)[ids[3]]-(*pU)[ids[0]])) / T(6);
}   
//-----------------------------------------------------------------------------
template <typename T>
T   Tetrahedron<T>::CalUndeformedVolume ()
{
    Volume = ((((*pX)[ids[1]]-(*pX)[ids[0]]) ^ ((*pX)[ids[2]]-(*pX)[ids[0]])) * ((*pX)[ids[3]]-(*pX)[ids[0]])) / T(6);
    return Volume;
}   
//-----------------------------------------------------------------------------
template <typename T>
void Tetrahedron<T>::UpdateStrainDisplacementMatrixElements ()
{
    // B = [ B0 B1 B2 B3 ]
    // where Bn (n = 0,1,2,3) are
    //      | bn  0  0 |
    //      |  0 cn  0 |
    // Bn = |  0  0 dn |
    //      | cn bn  0 |
    //      |  0 dn cn |
    //      | dn  0 bn |
    //
    // bn, cn, and dn are
    // | a0 b0 c0 d0 |                  |  1  1  1  1 |
    // | a1 b1 c1 d1 | = the inverse of | x0 x1 x2 x3 |
    // | a2 b2 c2 d2 |                  | y1 y2 y3 y4 |
    // | a3 b3 c3 d3 |                  | z1 z2 z3 z4 |
    //

    // bn, cn, and dn, where n = 1,2,3, can be computed more efficiently as
    // | b1 c1 d1 |                  | x1-x0 x2-x0 x3-x0 |
    // | b2 c2 d2 | = the inverse of | y1-y0 y2-y0 y3-y0 |
    // | b3 c3 d3 |                  | z1-z0 z2-z0 z3-z0 |
    //
    T e[9] = {
        ((*pX)[ids[1]])[0]-((*pX)[ids[0]])[0], ((*pX)[ids[2]])[0]-((*pX)[ids[0]])[0], ((*pX)[ids[3]])[0]-((*pX)[ids[0]])[0], 
        ((*pX)[ids[1]])[1]-((*pX)[ids[0]])[1], ((*pX)[ids[2]])[1]-((*pX)[ids[0]])[1], ((*pX)[ids[3]])[1]-((*pX)[ids[0]])[1], 
        ((*pX)[ids[1]])[2]-((*pX)[ids[0]])[2], ((*pX)[ids[2]])[2]-((*pX)[ids[0]])[2], ((*pX)[ids[3]])[2]-((*pX)[ids[0]])[2]
    };
    b[1] = e[4]*e[8] - e[5]*e[7];
    c[1] = e[2]*e[7] - e[1]*e[8];
    d[1] = e[1]*e[5] - e[2]*e[4];
    T detInv = T(1) / (e[0]*b[1] + e[3]*c[1] + e[6]*d[1]);
    b[1] *= detInv;
    c[1] *= detInv;
    d[1] *= detInv;
    b[2] = -( e[3]*e[8] - e[6]*e[5] ) * detInv;
    c[2] =  ( e[0]*e[8] - e[6]*e[2] ) * detInv;
    d[2] = -( e[0]*e[5] - e[3]*e[2] ) * detInv;
    b[3] =  ( e[3]*e[7] - e[6]*e[4] ) * detInv;
    c[3] = -( e[0]*e[7] - e[6]*e[1] ) * detInv;
    d[3] =  ( e[0]*e[4] - e[3]*e[1] ) * detInv;

    // Compute b0, c0, and d0 from
    // b0 = -b1-b2-b3
    // c0 = -c1-c2-c3
    // d0 = -d1-d2-d3
    b[0] = -(b[1]+b[2]+b[3]);
    c[0] = -(c[1]+c[2]+c[3]);
    d[0] = -(d[1]+d[2]+d[3]);

    /*
    // DEBUG
    std::cout << "B:\n" 
        << "\t" << b[0] << "\t" << c[0] << "\t" << d[0] << "\n"
        << "\t" << b[1] << "\t" << c[1] << "\t" << d[1] << "\n"
        << "\t" << b[2] << "\t" << c[2] << "\t" << d[2] << "\n"
        << "\t" << b[3] << "\t" << c[3] << "\t" << d[3] << "\n";

    //*
    Matrix3x3<T> test1( e );

    Matrix4x4<T> test2( 
        1, 1, 1, 1,
        ((*pX)[ids[0]])[0], ((*pX)[ids[1]])[0], ((*pX)[ids[2]])[0], ((*pX)[ids[3]])[0],
        ((*pX)[ids[0]])[1], ((*pX)[ids[1]])[1], ((*pX)[ids[2]])[1], ((*pX)[ids[3]])[1],
        ((*pX)[ids[0]])[2], ((*pX)[ids[1]])[2], ((*pX)[ids[2]])[2], ((*pX)[ids[3]])[2]
    );

    std::cout << test2.Inversed() << "\n";
    std::cout << test1.Inversed() << "\n";

    std::cout << "b[0] checked: " << test2[ 1] - b[0] << "\n";
    std::cout << "c[0] checked: " << test2[ 2] - c[0] << "\n";
    std::cout << "d[0] checked: " << test2[ 3] - d[0] << "\n";
    std::cout << "b[1] checked: " << test2[ 5] - b[1] << "\n";
    std::cout << "c[1] checked: " << test2[ 6] - c[1] << "\n";
    std::cout << "d[1] checked: " << test2[ 7] - d[1] << "\n";
    std::cout << "b[2] checked: " << test2[ 9] - b[2] << "\n";
    std::cout << "c[2] checked: " << test2[10] - c[2] << "\n";
    std::cout << "d[2] checked: " << test2[11] - d[2] << "\n";
    std::cout << "b[3] checked: " << test2[13] - b[3] << "\n";
    std::cout << "c[3] checked: " << test2[14] - c[3] << "\n";
    std::cout << "d[3] checked: " << test2[15] - d[3] << "\n";

    //std::cout << "\n4x4 - 3x3 RESULTS:\n";
    //std::cout << "b[1] checked: " << test2[ 5] - test1[0] << "\n";
    //std::cout << "c[1] checked: " << test2[ 6] - test1[1] << "\n";
    //std::cout << "d[1] checked: " << test2[ 7] - test1[2] << "\n";
    //std::cout << "b[2] checked: " << test2[ 9] - test1[3] << "\n";
    //std::cout << "c[2] checked: " << test2[10] - test1[4] << "\n";
    //std::cout << "d[2] checked: " << test2[11] - test1[5] << "\n";
    //std::cout << "b[3] checked: " << test2[13] - test1[6] << "\n";
    //std::cout << "c[3] checked: " << test2[14] - test1[7] << "\n";
    //std::cout << "d[3] checked: " << test2[15] - test1[8] << "\n";

    std::cout << "b[0] + b[1] + b[2] + b[3] = 0 = " << b[0]+b[1]+b[2]+b[3] << "\n";
    std::cout << "c[0] + c[1] + c[2] + c[3] = 0 = " << c[0]+c[1]+c[2]+c[3] << "\n";
    std::cout << "d[0] + d[1] + d[2] + d[3] = 0 = " << d[0]+d[1]+d[2]+d[3] << "\n";
    std::cout << "-----------------------------------------\n\n";
    //*/
}
//-----------------------------------------------------------------------------
template <typename T>
void Tetrahedron<T>::UpdateStiffnessMatrixElements ()
{
    // K^e_{nm} =
    //   | bn  0  0 | | e f f | | bm  0  0 |   | cn  0 dn | | g 0 0 | | cm bm  0 |
    //   |  0 cn  0 |*| f e f |*|  0 cm  0 | + | bn dn  0 |*| 0 g 0 |*|  0 dm cm |
    //   |  0  0 dn | | f f e | |  0  0 dm |   |  0 cn bn | | 0 0 g | | dm  0 bm |
    // =
    //   | bn*e*bm + cn*g*cm + dn*g*dm      bn*f*cm + cn*g*bm           bn*f*dm + dn*g*bm      |
    //   |      cn*f*bm + bn*g*cm      cn*e*cm + bn*g*bm + dn*g*dm      cn*f*dm + dn*g*cm      |
    //   |      dn*f*bm + bn*g*dm           dn*f*cm + cn*g*dm      dn*e*dm + cn*g*cm + bn*g*bm |

    T be[4] = { b[0]*Ee, b[1]*Ee, b[2]*Ee, b[3]*Ee };
    T bf[4] = { b[0]*Ef, b[1]*Ef, b[2]*Ef, b[3]*Ef };
    T bg[4] = { b[0]*Eg, b[1]*Eg, b[2]*Eg, b[3]*Eg };

    T ce[4] = { c[0]*Ee, c[1]*Ee, c[2]*Ee, c[3]*Ee };
    T cf[4] = { c[0]*Ef, c[1]*Ef, c[2]*Ef, c[3]*Ef };
    T cg[4] = { c[0]*Eg, c[1]*Eg, c[2]*Eg, c[3]*Eg };

    T de[4] = { d[0]*Ee, d[1]*Ee, d[2]*Ee, d[3]*Ee };
    T df[4] = { d[0]*Ef, d[1]*Ef, d[2]*Ef, d[3]*Ef };
    T dg[4] = { d[0]*Eg, d[1]*Eg, d[2]*Eg, d[3]*Eg };

    for ( int n = 0; n < 4; ++n ) {
        for ( int m = 0; m < 4; ++m ) {
            if ( n <= m ) {
                // 1st row
                Ke[n][m](0,0) = be[n]*b[m] + cg[n]*c[m] + dg[n]*d[m];
                Ke[n][m](0,1) = bf[n]*c[m] + cg[n]*b[m];
                Ke[n][m](0,2) = bf[n]*d[m] + dg[n]*b[m];
                // 2nd row
                Ke[n][m](1,0) = cf[n]*b[m] + bg[n]*c[m];
                Ke[n][m](1,1) = ce[n]*c[m] + bg[n]*b[m] + dg[n]*d[m];
                Ke[n][m](1,2) = cf[n]*d[m] + dg[n]*c[m];
                // 3rd row
                Ke[n][m](2,0) = df[n]*b[m] + bg[n]*d[m];
                Ke[n][m](2,1) = df[n]*c[m] + cg[n]*d[m];
                Ke[n][m](2,2) = de[n]*d[m] + cg[n]*c[m] + bg[n]*b[m];
            }
            // due to the symmetry of the linear tetrahedron's stiffness sub-matrices
            else {
                Ke[n][m] = Ke[m][n].GetTranspose();
            }
        }
    }
}
//-----------------------------------------------------------------------------
//=============================================================================




//=============================================================================
// OpenGL
#if defined(__gl_h_) || defined(__GL_H__)
//-----------------------------------------------------------------------------
template <typename T>
void Tetrahedron<T>::Draw ()
{
    glPushAttrib( GL_LINE_BIT );

    // Draw the undeformed element
    glLineWidth( 1 );
    glBegin( GL_LINE_LOOP );
        glColor3f( 1.0f, 0.0f, 0.0f );
        glVertex3fv( (*pX)[ids[0]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[1]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[3]].GetDataFloat() );
    glEnd();
    glBegin( GL_LINE_STRIP );
        glColor3f( 0.0f, 1.0f, 0.0f );
        glVertex3fv( (*pX)[ids[1]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[2]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[3]].GetDataFloat() );
    glEnd();
    glBegin( GL_LINES );
        glColor3f( 0.0f, 0.0f, 1.0f );
        glVertex3fv( (*pX)[ids[2]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[0]].GetDataFloat() );
        //glVertex3fv( (*pX)[ids[3]].GetDataFloat() );
    glEnd();

    // Draw the deformed element
    glLineWidth( 3 );
    glBegin( GL_LINE_LOOP );
        glColor3f( 1.0f, 0.0f, 0.0f );
        glVertex3fv( (*pU)[ids[0]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[1]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[3]].GetDataFloat() );
    glEnd();
    glBegin( GL_LINE_STRIP );
        glColor3f( 0.0f, 1.0f, 0.0f );
        glVertex3fv( (*pU)[ids[1]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[2]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[3]].GetDataFloat() );
    glEnd();
    glBegin( GL_LINES );
        glColor3f( 0.0f, 0.0f, 1.0f );
        glVertex3fv( (*pU)[ids[2]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[0]].GetDataFloat() );
        //glVertex3fv( (*pU)[ids[3]].GetDataFloat() );
    glEnd();

    glPopAttrib();
}


//-------------------------------------------------------------------
template <typename T>
void Tetrahedron<T>::DrawUndeformedMesh ()
{
    glPushAttrib( GL_LINE_BIT );
    // Draw the undeformed element
    glLineWidth( 1 );
    glBegin( GL_LINE_LOOP );
        glColor3f( 1.0f, 0.0f, 0.0f );
        glVertex3fv( (*pX)[ids[0]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[1]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[3]].GetDataFloat() );
    glEnd();
    glBegin( GL_LINE_STRIP );
        glColor3f( 0.0f, 1.0f, 0.0f );
        glVertex3fv( (*pX)[ids[1]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[2]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[3]].GetDataFloat() );
    glEnd();
    glBegin( GL_LINES );
        glColor3f( 0.0f, 0.0f, 1.0f );
        glVertex3fv( (*pX)[ids[2]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[0]].GetDataFloat() );
        //glVertex3fv( (*pX)[ids[3]].GetDataFloat() );
    glEnd();
    glPopAttrib();
}


//-------------------------------------------------------------------
template <typename T>
void Tetrahedron<T>::DrawDeformedMesh ()
{
    glPushAttrib( GL_LINE_BIT );
    // Draw the deformed element
    glLineWidth( 1 );
    glBegin( GL_LINE_LOOP );
        glColor3f( 1.0f, 0.0f, 0.0f );
        glVertex3fv( (*pU)[ids[0]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[1]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[3]].GetDataFloat() );
    glEnd();
    glBegin( GL_LINE_STRIP );
        glColor3f( 0.0f, 1.0f, 0.0f );
        glVertex3fv( (*pU)[ids[1]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[2]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[3]].GetDataFloat() );
    glEnd();
    glBegin( GL_LINES );
        glColor3f( 0.0f, 0.0f, 1.0f );
        glVertex3fv( (*pU)[ids[2]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[0]].GetDataFloat() );
        //glVertex3fv( (*pU)[ids[3]].GetDataFloat() );
    glEnd();
    glPopAttrib();
}


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