TAPsFEMHexahedron.cpp

Go to the documentation of this file.

/******************************************************************************
TAPsFEMHexahedron.cpp
******************************************************************************/
/******************************************************************************
SUKITTI PUNAK   (12/30/2009)
UPDATE          (07/01/2010)
******************************************************************************/
#include "TAPsFEMHexahedron.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
//=============================================================================
// Static Data Members
//-----------------------------------------------------------------------------

template <typename T> bool Hexahedron<T>::g_Is1PtGaussPointInit = false;
template <typename T> bool Hexahedron<T>::g_Is2PtGaussPointInit = false;
template <typename T> T Hexahedron<T>::g_1PtGauss_Der[8*3];
template <typename T> T Hexahedron<T>::g_2PtGauss_Der[8*3*8];

    // For finding an approximated parametric coordinates from a global coordinates
template <typename T> const int Hexahedron<T>::MAX_ITERATION = 10;
template <typename T> const T   Hexahedron<T>::ERROR_THRESHOLD = 1E-3;
template <typename T> const T   Hexahedron<T>::NEG_BOUND = -1 - ERROR_THRESHOLD;
template <typename T> const T   Hexahedron<T>::POS_BOUND =  1 + ERROR_THRESHOLD;
template <typename T> const T   Hexahedron<T>::DIVERGENCE_THRESHOLD = 1E6;


//-----------------------------------------------------------------------------
// Static Member Functions
//-----------------------------------------------------------------------------

template <typename T>
void Hexahedron<T>::CalAndStoreShapeFnsPartialDerivatives ( T store[], int ptNumber, T xi, T eta, T zeta )
{
    unsigned int idx = ptNumber * 24;
    InterpolateShapeFunctionDerivatives( xi, eta, zeta, &(store[idx]) );
}


template <typename T>
void Hexahedron<T>::UpdateShapeFnsPartialDerivatives_1PtGuassRule ()
{
    if ( g_Is1PtGaussPointInit )    return;
    g_Is1PtGaussPointInit = true;

    // 1-point Gauss quadrature rule has 2.0 weight at 0 point
    CalAndStoreShapeFnsPartialDerivatives( g_1PtGauss_Der, 0, 0, 0, 0 );
}


template <typename T>
void Hexahedron<T>::UpdateShapeFnsPartialDerivatives_2PtGuassRule ()
{
    if ( g_Is2PtGaussPointInit )    return;
    g_Is2PtGaussPointInit = true;

    // 2-point Gauss quadrature rule has 1.0 wight for both -sqrt(1/3) and +sqrt(1/3) points
    T pt2 = sqrt(T(1)/T(3));
    T pt1 = -pt2;

    CalAndStoreShapeFnsPartialDerivatives( g_2PtGauss_Der, 0, pt1, pt1, pt1 );
    CalAndStoreShapeFnsPartialDerivatives( g_2PtGauss_Der, 1, pt1, pt1, pt2 );
    CalAndStoreShapeFnsPartialDerivatives( g_2PtGauss_Der, 2, pt1, pt2, pt1 );
    CalAndStoreShapeFnsPartialDerivatives( g_2PtGauss_Der, 3, pt1, pt2, pt2 );
    CalAndStoreShapeFnsPartialDerivatives( g_2PtGauss_Der, 4, pt2, pt1, pt1 );
    CalAndStoreShapeFnsPartialDerivatives( g_2PtGauss_Der, 5, pt2, pt1, pt2 );
    CalAndStoreShapeFnsPartialDerivatives( g_2PtGauss_Der, 6, pt2, pt2, pt1 );
    CalAndStoreShapeFnsPartialDerivatives( g_2PtGauss_Der, 7, pt2, pt2, pt2 );

    // The order of points does not effect the computation!
    //CalAndStoreShapeFnsPartialDerivatives( g_2PtGauss_Der, 0, pt1, pt1, pt1 );
    //CalAndStoreShapeFnsPartialDerivatives( g_2PtGauss_Der, 1, pt2, pt1, pt1 );
    //CalAndStoreShapeFnsPartialDerivatives( g_2PtGauss_Der, 2, pt2, pt2, pt1 );
    //CalAndStoreShapeFnsPartialDerivatives( g_2PtGauss_Der, 3, pt1, pt2, pt1 );
    //CalAndStoreShapeFnsPartialDerivatives( g_2PtGauss_Der, 4, pt1, pt1, pt2 );
    //CalAndStoreShapeFnsPartialDerivatives( g_2PtGauss_Der, 5, pt2, pt1, pt2 );
    //CalAndStoreShapeFnsPartialDerivatives( g_2PtGauss_Der, 6, pt2, pt2, pt2 );
    //CalAndStoreShapeFnsPartialDerivatives( g_2PtGauss_Der, 7, pt1, pt2, pt2 );
}


template <typename T>
Matrix3x3<T> Hexahedron<T>::CalJacobian ( T GaussStore[], int GaussPtNumber )
{
    int i = GaussPtNumber*24;
    int j = i + 8;
    int k = j + 8;

    return Matrix3x3<T>(
        // x_i * \frac{\partial{N_i}\partial{xi}}
        (*pX)[ids[0]][0]*GaussStore[i+0] + 
        (*pX)[ids[1]][0]*GaussStore[i+1] + 
        (*pX)[ids[2]][0]*GaussStore[i+2] + 
        (*pX)[ids[3]][0]*GaussStore[i+3] + 
        (*pX)[ids[4]][0]*GaussStore[i+4] + 
        (*pX)[ids[5]][0]*GaussStore[i+5] + 
        (*pX)[ids[6]][0]*GaussStore[i+6] + 
        (*pX)[ids[7]][0]*GaussStore[i+7], 
        // y_i * \frac{\partial{N_i}\partial{xi}}
        (*pX)[ids[0]][1]*GaussStore[i+0] + 
        (*pX)[ids[1]][1]*GaussStore[i+1] + 
        (*pX)[ids[2]][1]*GaussStore[i+2] + 
        (*pX)[ids[3]][1]*GaussStore[i+3] + 
        (*pX)[ids[4]][1]*GaussStore[i+4] + 
        (*pX)[ids[5]][1]*GaussStore[i+5] + 
        (*pX)[ids[6]][1]*GaussStore[i+6] + 
        (*pX)[ids[7]][1]*GaussStore[i+7], 
        // z_i * \frac{\partial{N_i}\partial{xi}}
        (*pX)[ids[0]][2]*GaussStore[i+0] + 
        (*pX)[ids[1]][2]*GaussStore[i+1] + 
        (*pX)[ids[2]][2]*GaussStore[i+2] + 
        (*pX)[ids[3]][2]*GaussStore[i+3] + 
        (*pX)[ids[4]][2]*GaussStore[i+4] + 
        (*pX)[ids[5]][2]*GaussStore[i+5] + 
        (*pX)[ids[6]][2]*GaussStore[i+6] + 
        (*pX)[ids[7]][2]*GaussStore[i+7], 

        // x_i * \frac{\partial{N_i}\partial{eta}}
        (*pX)[ids[0]][0]*GaussStore[j+0] + 
        (*pX)[ids[1]][0]*GaussStore[j+1] + 
        (*pX)[ids[2]][0]*GaussStore[j+2] + 
        (*pX)[ids[3]][0]*GaussStore[j+3] + 
        (*pX)[ids[4]][0]*GaussStore[j+4] + 
        (*pX)[ids[5]][0]*GaussStore[j+5] + 
        (*pX)[ids[6]][0]*GaussStore[j+6] + 
        (*pX)[ids[7]][0]*GaussStore[j+7], 
        // y_i * \frac{\partial{N_i}\partial{eta}}
        (*pX)[ids[0]][1]*GaussStore[j+0] + 
        (*pX)[ids[1]][1]*GaussStore[j+1] + 
        (*pX)[ids[2]][1]*GaussStore[j+2] + 
        (*pX)[ids[3]][1]*GaussStore[j+3] + 
        (*pX)[ids[4]][1]*GaussStore[j+4] + 
        (*pX)[ids[5]][1]*GaussStore[j+5] + 
        (*pX)[ids[6]][1]*GaussStore[j+6] + 
        (*pX)[ids[7]][1]*GaussStore[j+7], 
        // z_i * \frac{\partial{N_i}\partial{eta}}
        (*pX)[ids[0]][2]*GaussStore[j+0] + 
        (*pX)[ids[1]][2]*GaussStore[j+1] + 
        (*pX)[ids[2]][2]*GaussStore[j+2] + 
        (*pX)[ids[3]][2]*GaussStore[j+3] + 
        (*pX)[ids[4]][2]*GaussStore[j+4] + 
        (*pX)[ids[5]][2]*GaussStore[j+5] + 
        (*pX)[ids[6]][2]*GaussStore[j+6] + 
        (*pX)[ids[7]][2]*GaussStore[j+7], 

        // x_i * \frac{\partial{N_i}\partial{zeta}}
        (*pX)[ids[0]][0]*GaussStore[k+0] + 
        (*pX)[ids[1]][0]*GaussStore[k+1] + 
        (*pX)[ids[2]][0]*GaussStore[k+2] + 
        (*pX)[ids[3]][0]*GaussStore[k+3] + 
        (*pX)[ids[4]][0]*GaussStore[k+4] + 
        (*pX)[ids[5]][0]*GaussStore[k+5] + 
        (*pX)[ids[6]][0]*GaussStore[k+6] + 
        (*pX)[ids[7]][0]*GaussStore[k+7], 
        // y_i * \frac{\partial{N_i}\partial{zeta}}
        (*pX)[ids[0]][1]*GaussStore[k+0] + 
        (*pX)[ids[1]][1]*GaussStore[k+1] + 
        (*pX)[ids[2]][1]*GaussStore[k+2] + 
        (*pX)[ids[3]][1]*GaussStore[k+3] + 
        (*pX)[ids[4]][1]*GaussStore[k+4] + 
        (*pX)[ids[5]][1]*GaussStore[k+5] + 
        (*pX)[ids[6]][1]*GaussStore[k+6] + 
        (*pX)[ids[7]][1]*GaussStore[k+7], 
        // z_i * \frac{\partial{N_i}\partial{zeta}}
        (*pX)[ids[0]][2]*GaussStore[k+0] + 
        (*pX)[ids[1]][2]*GaussStore[k+1] + 
        (*pX)[ids[2]][2]*GaussStore[k+2] + 
        (*pX)[ids[3]][2]*GaussStore[k+3] + 
        (*pX)[ids[4]][2]*GaussStore[k+4] + 
        (*pX)[ids[5]][2]*GaussStore[k+5] + 
        (*pX)[ids[6]][2]*GaussStore[k+6] + 
        (*pX)[ids[7]][2]*GaussStore[k+7] 
    );
}

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


//=============================================================================
// Constructors
//-----------------------------------------------------------------------------
template <typename T>
Hexahedron<T>::Hexahedron ( 
        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,   
        unsigned int globalNode4,   
        unsigned int globalNode5,   
        unsigned int globalNode6,   
        unsigned int globalNode7    
) : Element( deformedNodes, undeformedNodes, YoungModulus, PoissionRatio )
{
    ids[0] = globalNode0;
    ids[1] = globalNode1;
    ids[2] = globalNode2;
    ids[3] = globalNode3;
    ids[4] = globalNode4;
    ids[5] = globalNode5;
    ids[6] = globalNode6;
    ids[7] = globalNode7;
    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

    // Strain Displacement Matrices are computed and never stored
    // Therefore, they are computed inside UpdateStiffnessMatrixElements function
    //UpdateStrainDisplacementMatrixElements(); // update the elements of the strain-displacement matrix

    UpdateStiffnessMatrixElements();    // update all 16 sub stiffness matrices of this element
}
//-----------------------------------------------------------------------------
//template <typename T>
//Hexahedron<T>::Hexahedron ( Hexahedron<T> const &orig )
//{}
//-----------------------------------------------------------------------------
template <typename T>
Hexahedron<T>::~Hexahedron ()
{}  
//-----------------------------------------------------------------------------
template <typename T>
std::string Hexahedron<T>::StrInfo () const
{
    std::stringstream ss;
    ss  << "FEM::Hexahedron<" << typeid(T).name() << ">:\n"
        << "\tDeformed - Undeformed Node 0: " << (*pU)[ids[0]] << " - " << (*pX)[ids[0]] << "\n"
        << "\tDeformed - Undeformed Node 1: " << (*pU)[ids[1]] << " - " << (*pX)[ids[1]] << "\n"
        << "\tDeformed - Undeformed Node 2: " << (*pU)[ids[2]] << " - " << (*pX)[ids[2]] << "\n"
        << "\tDeformed - Undeformed Node 3: " << (*pU)[ids[3]] << " - " << (*pX)[ids[3]] << "\n"
        << "\tDeformed - Undeformed Node 4: " << (*pU)[ids[4]] << " - " << (*pX)[ids[4]] << "\n"
        << "\tDeformed - Undeformed Node 5: " << (*pU)[ids[5]] << " - " << (*pX)[ids[5]] << "\n"
        << "\tDeformed - Undeformed Node 6: " << (*pU)[ids[6]] << " - " << (*pX)[ids[6]] << "\n"
        << "\tDeformed - Undeformed Node 7: " << (*pU)[ids[7]] << " - " << (*pX)[ids[7]] << "\n"
        << "\tDeformed - Undeformed Volume: " << CalDeformedVolume()  << " - " << Volume << "\n";
    return ss.str();
}
//-----------------------------------------------------------------------------
//=============================================================================
// Assignment Operator
//-----------------------------------------------------------------------------
//template <typename T>
//Hexahedron<T> & Hexahedron<T>::operator= ( Hexahedron<T> const &orig )
//{}
//-----------------------------------------------------------------------------
//=============================================================================
// Get/Set Functions
//-----------------------------------------------------------------------------
template <typename T>
Matrix3x3<T> const & Hexahedron<T>::GetStiffnessMatrix ( unsigned int n, unsigned int m ) const
{
    assert( 0 <= n && n <= 7 );
    assert( 0 <= m && m <= 7 );
    return Ke[n][m];
}   
//-----------------------------------------------------------------------------
//=============================================================================
// Operations
//-----------------------------------------------------------------------------
template <typename T>
void Hexahedron<T>::SetNode ( unsigned int i, unsigned int globalNodeNumber )
{
    assert( 0 <= i && i <= 7 );
    ids[i] = globalNodeNumber;
}   
//-----------------------------------------------------------------------------
//template <typename T>
//void Hexahedron<T>::ResetToDeformedState ()
//{
//} 
//-----------------------------------------------------------------------------
template <typename T>
void Hexahedron<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]];
    (*pU)[ids[4]] = (*pX)[ids[4]];
    (*pU)[ids[5]] = (*pX)[ids[5]];
    (*pU)[ids[6]] = (*pX)[ids[6]];
    (*pU)[ids[7]] = (*pX)[ids[7]];
}   
//-----------------------------------------------------------------------------
template <typename T>
Vector3<T> * const Hexahedron<T>::DeformedNode ( unsigned int i )
{
    assert( 0 <= i && i <= 7 );
    return &(*pU)[ids[i]];
}   
//-----------------------------------------------------------------------------
template <typename T>
Vector3<T> * const Hexahedron<T>::UndeformedNode ( unsigned int i )
{
    assert( 0 <= i && i <= 7 );
    return &(*pX)[ids[i]];
}
//-----------------------------------------------------------------------------
template <typename T>
Vector3<T> Hexahedron<T>::ComputeDeformedPositionBasedOnCoordinates ( T xi, T eta, T zeta )
{
    T s0 = 1 - xi;
    T s1 = 1 + xi;
    T t0 = 1 - eta;
    T t1 = 1 + eta;
    T r0 = 1 - zeta;
    T r1 = 1 + zeta;

    T N0 = s0 * t0 * r0;
    T N1 = s1 * t0 * r0;
    T N2 = s1 * t1 * r0;
    T N3 = s0 * t1 * r0;
    T N4 = s0 * t0 * r1;
    T N5 = s1 * t0 * r1;
    T N6 = s1 * t1 * r1;
    T N7 = s0 * t1 * r1;

    T x = 0.125 * ( (*pU)[ids[0]][0]*N0 + (*pU)[ids[1]][0]*N1 + (*pU)[ids[2]][0]*N2 + (*pU)[ids[3]][0]*N3 
                  + (*pU)[ids[4]][0]*N4 + (*pU)[ids[5]][0]*N5 + (*pU)[ids[6]][0]*N6 + (*pU)[ids[7]][0]*N7 );
    T y = 0.125 * ( (*pU)[ids[0]][1]*N0 + (*pU)[ids[1]][1]*N1 + (*pU)[ids[2]][1]*N2 + (*pU)[ids[3]][1]*N3 
                  + (*pU)[ids[4]][1]*N4 + (*pU)[ids[5]][1]*N5 + (*pU)[ids[6]][1]*N6 + (*pU)[ids[7]][1]*N7 );
    T z = 0.125 * ( (*pU)[ids[0]][2]*N0 + (*pU)[ids[1]][2]*N1 + (*pU)[ids[2]][2]*N2 + (*pU)[ids[3]][2]*N3 
                  + (*pU)[ids[4]][2]*N4 + (*pU)[ids[5]][2]*N5 + (*pU)[ids[6]][2]*N6 + (*pU)[ids[7]][2]*N7 );

    return Vector3<T>( x, y, z );
}
//-----------------------------------------------------------------------------
template <typename T>
Vector3<T> Hexahedron<T>::ComputeUndeformedPositionBasedOnCoordinates ( T xi, T eta, T zeta )
{
    T s0 = 1 - xi;
    T s1 = 1 + xi;
    T t0 = 1 - eta;
    T t1 = 1 + eta;
    T r0 = 1 - zeta;
    T r1 = 1 + zeta;

    T N0 = s0 * t0 * r0;
    T N1 = s1 * t0 * r0;
    T N2 = s1 * t1 * r0;
    T N3 = s0 * t1 * r0;
    T N4 = s0 * t0 * r1;
    T N5 = s1 * t0 * r1;
    T N6 = s1 * t1 * r1;
    T N7 = s0 * t1 * r1;

    T x = 0.125 * ( (*pX)[ids[0]][0]*N0 + (*pX)[ids[1]][0]*N1 + (*pX)[ids[2]][0]*N2 + (*pX)[ids[3]][0]*N3 
                  + (*pX)[ids[4]][0]*N4 + (*pX)[ids[5]][0]*N5 + (*pX)[ids[6]][0]*N6 + (*pX)[ids[7]][0]*N7 );
    T y = 0.125 * ( (*pX)[ids[0]][1]*N0 + (*pX)[ids[1]][1]*N1 + (*pX)[ids[2]][1]*N2 + (*pX)[ids[3]][1]*N3 
                  + (*pX)[ids[4]][1]*N4 + (*pX)[ids[5]][1]*N5 + (*pX)[ids[6]][1]*N6 + (*pX)[ids[7]][1]*N7 );
    T z = 0.125 * ( (*pX)[ids[0]][2]*N0 + (*pX)[ids[1]][2]*N1 + (*pX)[ids[2]][2]*N2 + (*pX)[ids[3]][2]*N3 
                  + (*pX)[ids[4]][2]*N4 + (*pX)[ids[5]][2]*N5 + (*pX)[ids[6]][2]*N6 + (*pX)[ids[7]][2]*N7 );

    return Vector3<T>( x, y, z );
}
//-----------------------------------------------------------------------------
template <typename T>
int Hexahedron<T>::EstimateParametricCoords (
    T x, T y, T z,          
    T &xi, T &eta, T &zeta  
)
{
    // Using Newton's method

    // Initialize parametric coordinates
    xi = eta = zeta = 0.5;
    T px = xi;
    T py = eta;
    T pz = zeta;
    Vector3<T> C[4];
    //Vector3<T> Point;
    T weights[8];
    T derivs[24];
    bool bFoundAGoodOne = false;

    // Iteration loop
    for ( int iteration = 0; !bFoundAGoodOne && iteration < MAX_ITERATION; ++iteration ) {
        //Point = ComputeDeformedPositionBasedOnCoordinates( xi, eta, zeta );
        InterpolateShapeFunctions( xi, eta, zeta, weights );
        InterpolateShapeFunctionDerivatives( xi, eta, zeta, derivs );

        // Newton function
        //-------------------------------------------------
        for ( int i = 0; i < 4; ++i ) {
            C[i].Clear();   // clear data for matrix
        }
        for ( int p = 0; p < 8; ++p ) {
            C[0] += weights[p]    * (*pU)[ids[p]];
            C[1] += derivs[p    ] * (*pU)[ids[p]];
            C[2] += derivs[p + 8] * (*pU)[ids[p]];
            C[3] += derivs[p +16] * (*pU)[ids[p]];
        }
        C[0][0] -= x;
        C[0][1] -= y;
        C[0][2] -= z;

        // Compute determinants and find better estimation
        T det = Matrix3x3<T>(   C[1][0], C[2][0], C[3][0],
                                C[1][1], C[2][1], C[3][1],
                                C[1][2], C[2][2], C[3][2] ).GetDeterminant();

        if ( fabs( det ) < 1E-20 ) {
            return -1;  // couldn't find a good estimation
        }

        xi   = px - Matrix3x3<T>(   C[0][0], C[2][0], C[3][0],
                                    C[0][1], C[2][1], C[3][1],
                                    C[0][2], C[2][2], C[3][2] ).GetDeterminant() / det;
        eta  = py - Matrix3x3<T>(   C[1][0], C[0][0], C[3][0],
                                    C[1][1], C[0][1], C[3][1],
                                    C[1][2], C[0][2], C[3][2] ).GetDeterminant() / det;
        zeta = pz - Matrix3x3<T>(   C[1][0], C[2][0], C[0][0],
                                    C[1][1], C[2][1], C[0][1],
                                    C[1][2], C[2][2], C[0][2] ).GetDeterminant() / det;
        //-------------------------------------------------

        // Check for divergence
        if (    fabs(xi) > DIVERGENCE_THRESHOLD
            ||  fabs(eta) > DIVERGENCE_THRESHOLD
            ||  fabs(zeta) > DIVERGENCE_THRESHOLD )
        {
            return -2;  // failed to find a good estimation due to divergence
        }

        // If less than the error threshold, stop.
        if (    fabs(xi-px) < ERROR_THRESHOLD
            &&  fabs(eta-py) < ERROR_THRESHOLD
            &&  fabs(zeta-pz) < ERROR_THRESHOLD )
        {
            bFoundAGoodOne = true;
        }
        // Else, repeat the Newton function.
        else {
            px = xi;
            py = eta;
            pz = zeta;
        }
    }

    // If the max iteration is reached before a good estimation is found.
    if ( !bFoundAGoodOne ) {
        return -3;  // failed to find a good estimation before the max iteration is reached
    }

    // Check whether the estimated point is inside or outside the element
    //InterpolateShapeFunctions( xi, eta, zeta, weights );
    if (    NEG_BOUND <= xi   && xi   <= POS_BOUND
        &&  NEG_BOUND <= eta  && eta  <= POS_BOUND
        &&  NEG_BOUND <= zeta && zeta <= POS_BOUND )
    {
        return 1;   // the global point is inside the hexahedron.
    }
    else {
        return 0;   // the global point is outside the hexahedron.
    }
}
//-----------------------------------------------------------------------------
template <typename T>
T   Hexahedron<T>::CalDeformedVolume () const
{
    std::cout << "Hexahedron<T>::CalDeformedVolume isn't implemented yet!!!\n";
    return 1;
}   
//-----------------------------------------------------------------------------
template <typename T>
T   Hexahedron<T>::CalUndeformedVolume ()
{
    // Assume the undeformed hexahedron is a parallelepiped, then
    // Volume = AxB.C
    Vector3<T> crossProd = (((*pX)[ids[1]]-(*pX)[ids[0]])^((*pX)[ids[3]]-(*pX)[ids[0]]));
    Volume = (((*pX)[ids[1]]-(*pX)[ids[0]])^((*pX)[ids[3]]-(*pX)[ids[0]])) * ((*pX)[ids[4]]-(*pX)[ids[0]]);
    return Volume;
}   
//-----------------------------------------------------------------------------
template <typename T>
void Hexahedron<T>::UpdateStrainDisplacementMatrixElements ()
{
    // Too many variables to be kept for updating
    // So the update is changed to instant computation and added into 
    // UpdateStiffnessMatrixElements.

    // compute the elements of B
    // which are Nn_dx, Nn_dy, and Nn_dz for n = 0...7

    // B = [ B0 B1 B2 B3 B4 B5 B6 B7 ]
    // where Bn (n = 0,1,2,3,4,5,6,7) are
    //      | bn  0  0 |
    //      |  0 cn  0 |
    // Bn = |  0  0 dn |
    //      | cn bn  0 |
    //      |  0 dn cn |
    //      | dn  0 bn |
    // where 
    // | bn |            | Nn_dxi   |
    // | cn | = J_n^{-1} | Nn_deta  |
    // | dn |            | Nn_dzeta |

    // PICK ONE OF THESE
    //UpdateStrainDisplacementMatrixElements_1PtGauss();
    //UpdateStrainDisplacementMatrixElements_2PtGauss();
}
//-----------------------------------------------------------------------------
template <typename T>
void Hexahedron<T>::UpdateStiffnessMatrixElements ()
{
    // PICK ONE OF THESE
    //UpdateStiffnessMatrixElements_1PtGauss();
    UpdateStiffnessMatrixElements_2PtGauss();
}
//-----------------------------------------------------------------------------
template <typename T>
void Hexahedron<T>::UpdateStiffnessMatrixElements_1PtGauss ()
{
    T b[8], c[8], d[8];
    Matrix3x3<T> J;

    // Calculate J and its inverse
    J = CalJacobian( g_1PtGauss_Der, 0 );
    T detJ = J.GetDeterminant();
    Matrix3x3<T> Jinv = J.GetInverse();

    // Calculate B[n]
    for ( int n=0; n<8; ++n ) {
        // | bn |            | Nn_dxi   |
        // | cn | = J^{-1} | Nn_deta  |
        // | dn |            | Nn_dzeta |
        b[n] = Jinv(0,0)*g_1PtGauss_Der[n] + Jinv(0,1)*g_1PtGauss_Der[n+8] + Jinv(0,2)*g_1PtGauss_Der[n+16];
        c[n] = Jinv(1,0)*g_1PtGauss_Der[n] + Jinv(1,1)*g_1PtGauss_Der[n+8] + Jinv(1,2)*g_1PtGauss_Der[n+16];
        d[n] = Jinv(2,0)*g_1PtGauss_Der[n] + Jinv(2,1)*g_1PtGauss_Der[n+8] + Jinv(2,2)*g_1PtGauss_Der[n+16];
    }

    // K^e_{nm} = Wi * Wj * Wk
    //   *
    //         | bn  0  0 cn  0 dn |
    //  Bn^T = |  0 cn  0 bn dn  0 |
    //         |  0  0 dn  0 cn bn |
    //   *
    //      | e f f 0 0 0 |
    //      | f e f 0 0 0 |
    //  E = | f f e 0 0 0 |
    //      | 0 0 0 g 0 0 |
    //      | 0 0 0 0 g 0 |
    //      | 0 0 0 0 0 g |
    //   *
    //       | bn  0  0 |
    //       |  0 cn  0 |
    //  Bn = |  0  0 dn |
    //       | cn bn  0 |
    //       |  0 dn cn |
    //       | dn  0 bn |
    //   *
    //  det(J)
    //
    // which is equivalent to
    //
    // K^e_{nm} = Wi * Wj * Wk *
    //    | 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 |
    //   * det(J)
    //  = Wi * Wj * Wk *
    //    | 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 |
    //   * det(J)
    T be[8] = { b[0]*Ee, b[1]*Ee, b[2]*Ee, b[3]*Ee, b[4]*Ee, b[5]*Ee, b[6]*Ee, b[7]*Ee };
    T bf[8] = { b[0]*Ef, b[1]*Ef, b[2]*Ef, b[3]*Ef, b[4]*Ef, b[5]*Ef, b[6]*Ef, b[7]*Ef };
    T bg[8] = { b[0]*Eg, b[1]*Eg, b[2]*Eg, b[3]*Eg, b[4]*Eg, b[5]*Eg, b[6]*Eg, b[7]*Eg };

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

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

    // Calculate K^e_{nm}
    for ( int n = 0; n < 8; ++n ) {
        for ( int m = 0; m < 8; ++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];

                Ke[n][m] *= 8 * detJ;   // multiply by weights and detJ
            }
            // due to the symmetry of the linear tetrahedron's stiffness sub-matrices
            else {
                Ke[n][m] = Ke[m][n].GetTranspose();
            }
        }
    }
}
//-----------------------------------------------------------------------------
template <typename T>
void Hexahedron<T>::UpdateStiffnessMatrixElements_2PtGauss ()
{
    // Clear K^e_{nm}
    for ( int n = 0; n < 8; ++n ) {
        for ( int m = 0; m < 8; ++m ) {
            Ke[n][m].MakeZero();
        }
    }

    // 2-pt Gauss quadrature for 3D becomes 8 points
    T b[8], c[8], d[8];

    for ( int pt = 0, idx = 0; pt < 8; ++pt, idx+=24 ) {

        // Calculate J[n] and B[n]
        Matrix3x3<T> J = CalJacobian( g_2PtGauss_Der, pt );
        T detJ = J.GetDeterminant();
        Matrix3x3<T> Jinv = J.GetInverse();

        for ( int n=0, idx2=idx; n<8; ++n, ++idx2 ) {
            // | bn |            | Nn_dxi   |
            // | cn | = J^{-1} | Nn_deta  |
            // | dn |            | Nn_dzeta |
            b[n] = Jinv(0,0)*g_2PtGauss_Der[idx2] + Jinv(0,1)*g_2PtGauss_Der[idx2+8] + Jinv(0,2)*g_2PtGauss_Der[idx2+16];
            c[n] = Jinv(1,0)*g_2PtGauss_Der[idx2] + Jinv(1,1)*g_2PtGauss_Der[idx2+8] + Jinv(1,2)*g_2PtGauss_Der[idx2+16];
            d[n] = Jinv(2,0)*g_2PtGauss_Der[idx2] + Jinv(2,1)*g_2PtGauss_Der[idx2+8] + Jinv(2,2)*g_2PtGauss_Der[idx2+16];
        }

        // K^e_{nm} = Wi * Wj * Wk
        //   *
        //         | bn  0  0 cn  0 dn |
        //  Bn^T = |  0 cn  0 bn dn  0 |
        //         |  0  0 dn  0 cn bn |
        //   *
        //      | e f f 0 0 0 |
        //      | f e f 0 0 0 |
        //  E = | f f e 0 0 0 |
        //      | 0 0 0 g 0 0 |
        //      | 0 0 0 0 g 0 |
        //      | 0 0 0 0 0 g |
        //   *
        //       | bn  0  0 |
        //       |  0 cn  0 |
        //  Bn = |  0  0 dn |
        //       | cn bn  0 |
        //       |  0 dn cn |
        //       | dn  0 bn |
        //   *
        //  det(J)
        //
        // which is equivalent to
        //
        // K^e_{nm} = Wi * Wj * Wk *
        //    | 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 |
        //   * det(J)
        //  = Wi * Wj * Wk *
        //    | 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 |
        //   * det(J)
        T be[8] = { b[0]*Ee, b[1]*Ee, b[2]*Ee, b[3]*Ee, b[4]*Ee, b[5]*Ee, b[6]*Ee, b[7]*Ee };
        T bf[8] = { b[0]*Ef, b[1]*Ef, b[2]*Ef, b[3]*Ef, b[4]*Ef, b[5]*Ef, b[6]*Ef, b[7]*Ef };
        T bg[8] = { b[0]*Eg, b[1]*Eg, b[2]*Eg, b[3]*Eg, b[4]*Eg, b[5]*Eg, b[6]*Eg, b[7]*Eg };

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

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

        // Calculate K^e_{nm}
        for ( int n = 0; n < 8; ++n ) {
            for ( int m = 0; m < 8; ++m ) {
                if ( n <= m ) {
                    // 1st row
                    Ke[n][m](0,0) += (be[n]*b[m] + cg[n]*c[m] + dg[n]*d[m]) * detJ;
                    Ke[n][m](0,1) += (bf[n]*c[m] + cg[n]*b[m]) * detJ;
                    Ke[n][m](0,2) += (bf[n]*d[m] + dg[n]*b[m]) * detJ;
                    // 2nd row
                    Ke[n][m](1,0) += (cf[n]*b[m] + bg[n]*c[m]) * detJ;
                    Ke[n][m](1,1) += (ce[n]*c[m] + bg[n]*b[m] + dg[n]*d[m]) * detJ;
                    Ke[n][m](1,2) += (cf[n]*d[m] + dg[n]*c[m]) * detJ;
                    // 3rd row
                    Ke[n][m](2,0) += (df[n]*b[m] + bg[n]*d[m]) * detJ;
                    Ke[n][m](2,1) += (df[n]*c[m] + cg[n]*d[m]) * detJ;
                    Ke[n][m](2,2) += (de[n]*d[m] + cg[n]*c[m] + bg[n]*b[m]) * detJ;
                }
                // due to the symmetry of the linear tetrahedron's stiffness sub-matrices
                // will be done after this loop
            }
        }
    }

    // Calculate K^e_{nm} from symmetry
    for ( int n = 0; n < 8; ++n ) {
        for ( int m = 0; m < 8; ++m ) {
            if ( n <= m ) {
                //i.e. Ke[n][m] *= 1*1*1;   // multiply by weights
            }
            // due to the symmetry of the linear tetrahedron's stiffness sub-matrices
            else {
                Ke[n][m] = Ke[m][n].GetTranspose();
            }
        }
    }
}

//-----------------------------------------------------------------------------
template <typename T>
void Hexahedron<T>::InterpolateShapeFunctions (
    T xi,           
    T eta,          
    T zeta,         
    T weights[8]    
)
{
    T S = T(0.125);
    T a_0 = 1 - xi;  T b_0 = 1 - eta;  T c_0 = 1 - zeta;
    T a_1 = 1 + xi;  T b_1 = 1 + eta;  T c_1 = 1 + zeta;

    weights[0] = S * a_0 * b_0 * c_0;
    weights[1] = S * a_1 * b_0 * c_0;
    weights[2] = S * a_1 * b_1 * c_0;
    weights[3] = S * a_0 * b_1 * c_0;
    weights[4] = S * a_0 * b_0 * c_1;
    weights[5] = S * a_1 * b_0 * c_1;
    weights[6] = S * a_1 * b_1 * c_1;
    weights[7] = S * a_0 * b_1 * c_1;
}

//-----------------------------------------------------------------------------
template <typename T>
void Hexahedron<T>::InterpolateShapeFunctionDerivatives (
    T xi,           
    T eta,          
    T zeta,         
    T derivs[24]    
)
{
    T S = T(0.125);
    T a_0 = 1 - xi;  T b_0 = 1 - eta;  T c_0 = 1 - zeta;
    T a_1 = 1 + xi;  T b_1 = 1 + eta;  T c_1 = 1 + zeta;

    // N_0 to N_7 partial diff by xi
    derivs[ 0] = -S * b_0 * c_0;
    derivs[ 1] =  S * b_0 * c_0;
    derivs[ 2] =  S * b_1 * c_0;
    derivs[ 3] = -S * b_1 * c_0;
    derivs[ 4] = -S * b_0 * c_1;
    derivs[ 5] =  S * b_0 * c_1;
    derivs[ 6] =  S * b_1 * c_1;
    derivs[ 7] = -S * b_1 * c_1;

    // N_0 to N_7 partial diff by eta
    derivs[ 8] = -S * a_0 * c_0;
    derivs[ 9] = -S * a_1 * c_0;
    derivs[10] =  S * a_1 * c_0;
    derivs[11] =  S * a_0 * c_0;
    derivs[12] = -S * a_0 * c_1;
    derivs[13] = -S * a_1 * c_1;
    derivs[14] =  S * a_1 * c_1;
    derivs[15] =  S * a_0 * c_1;
    
    // N_0 to N_7 partial diff by zeta
    derivs[16] = -S * a_0 * b_0;
    derivs[17] = -S * a_1 * b_0;
    derivs[18] = -S * a_1 * b_1;
    derivs[19] = -S * a_0 * b_1;
    derivs[20] =  S * a_0 * b_0;
    derivs[21] =  S * a_1 * b_0;
    derivs[22] =  S * a_1 * b_1;
    derivs[23] =  S * a_0 * b_1;
}
//-----------------------------------------------------------------------------
//=============================================================================




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

    // Draw the undeformed element
    glLineWidth( 1 );
    //glPushMatrix();
    //glTranslatef( 0.5, 0.5, 0.5 );
    glBegin( GL_LINE_LOOP );
        glColor3f( 1.0f, 0.0f, 0.0f );
        glVertex3fv( (*pX)[ids[0]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[1]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[2]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[3]].GetDataFloat() );
    glEnd();
    glBegin( GL_LINES );
        glColor3f( 0.0f, 1.0f, 0.0f );
        glVertex3fv( (*pX)[ids[0]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[4]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[1]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[5]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[2]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[6]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[3]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[7]].GetDataFloat() );
    glEnd();
    glBegin( GL_LINE_LOOP );
        glColor3f( 0.0f, 0.0f, 1.0f );
        glVertex3fv( (*pX)[ids[4]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[5]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[6]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[7]].GetDataFloat() );
    glEnd();
    //glPopMatrix();

    // 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[2]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[3]].GetDataFloat() );
    glEnd();
    glBegin( GL_LINES );
        glColor3f( 0.0f, 1.0f, 0.0f );
        glVertex3fv( (*pU)[ids[0]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[4]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[1]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[5]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[2]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[6]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[3]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[7]].GetDataFloat() );
    glEnd();
    glBegin( GL_LINE_LOOP );
        glColor3f( 0.0f, 0.0f, 1.0f );
        glVertex3fv( (*pU)[ids[4]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[5]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[6]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[7]].GetDataFloat() );
    glEnd();

    glPopAttrib();
}


//-------------------------------------------------------------------
template <typename T>
void Hexahedron<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[2]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[3]].GetDataFloat() );
    glEnd();
    glBegin( GL_LINES );
        glColor3f( 0.0f, 1.0f, 0.0f );
        glVertex3fv( (*pX)[ids[0]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[4]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[1]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[5]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[2]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[6]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[3]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[7]].GetDataFloat() );
    glEnd();
    glBegin( GL_LINE_LOOP );
        glColor3f( 0.0f, 0.0f, 1.0f );
        glVertex3fv( (*pX)[ids[4]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[5]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[6]].GetDataFloat() );
        glVertex3fv( (*pX)[ids[7]].GetDataFloat() );
    glEnd();
    glPopAttrib();
}


//-------------------------------------------------------------------
template <typename T>
void Hexahedron<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[2]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[3]].GetDataFloat() );
    glEnd();
    glBegin( GL_LINES );
        glColor3f( 0.0f, 1.0f, 0.0f );
        glVertex3fv( (*pU)[ids[0]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[4]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[1]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[5]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[2]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[6]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[3]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[7]].GetDataFloat() );
    glEnd();
    glBegin( GL_LINE_LOOP );
        glColor3f( 0.0f, 0.0f, 1.0f );
        glVertex3fv( (*pU)[ids[4]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[5]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[6]].GetDataFloat() );
        glVertex3fv( (*pU)[ids[7]].GetDataFloat() );
    glEnd();
    glPopAttrib();
}
//-----------------------------------------------------------------------------
#endif // OpenGL
//=============================================================================
//-----------------------------------------------------------------------------
//=============================================================================
END_NAMESPACE_TAPs__FEM
//34567890123456789012345678901234567890123456789012345678901234567890123456789
//--+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----