TAPsFEMMeshHexahedraGen.cpp

Go to the documentation of this file.

/******************************************************************************
TAPsFEMMeshHexahedraGen.cpp
******************************************************************************/
/******************************************************************************
SUKITTI PUNAK   (02/18/2010)
UPDATE          (04/16/2010)
******************************************************************************/
#include "TAPsFEMMeshHexahedraGen.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
//=============================================================================
//template <typename T> std::vector< Hexahedron<T> > * MeshHexahedraGen<T>::m_pListOfElements = NULL;
//template <typename T> std::vector< Vector3<T> > * MeshHexahedraGen<T>::m_pListOfDeformedMeshNodes = NULL;
//template <typename T> std::vector< Vector3<T> > * MeshHexahedraGen<T>::m_pListOfUndeformedMeshNodes = NULL;
//-----------------------------------------------------------------------------

//template <typename T>
//void MeshHexahedraGen<T>::SetAccessors ( MeshHexahedra<T> & rHexFemMesh )
//{
//  m_pListOfElements               = &( rHexFemMesh.RetListOfElements() );
//  m_pListOfDeformedMeshNodes      = &( rHexFemMesh.RetListOfDeformedMeshNodes() );
//  m_pListOfUndeformedMeshNodes    = &( rHexFemMesh.RetListOfUndeformedMeshNodes() );
//}
//template <typename T>
//void MeshHexahedraGen<T>::ClearAccessors ()
//{
//  m_pListOfElements               = NULL;
//  m_pListOfDeformedMeshNodes      = NULL;
//  m_pListOfUndeformedMeshNodes    = NULL;
//}

//-----------------------------------------------------------------------------
template <typename T, typename DATA>
void MeshHexahedraGen<T,DATA>::OneElementMesh ( MeshHexahedra<T> & rHexFemMesh, T mass, T YoungsModulus, T PoissonsRatio )
{
    //    3---------2      ^ y
    //   /|        /|      |
    //  / |       / |      |
    // 7---------6  |      |
    // |  0- - - | -1      +------>x
    // | /       | /      /
    // |/        |/      /
    // 4---------5      v z

    std::vector< Hexahedron<T> > & E = rHexFemMesh.RetListOfElements();
    std::vector< Vector3<T> > & Nd = rHexFemMesh.RetListOfDeformedMeshNodes();
    std::vector< Vector3<T> > & Nu = rHexFemMesh.RetListOfUndeformedMeshNodes();

    // Clear the mesh
    E.clear();
    Nd.clear();
    Nu.clear();

    /*
    // 8 deformed nodes
    Nd.push_back( Vector3<T>(-1,-1,-1) );
    Nd.push_back( Vector3<T>( 1,-1,-1) );
    Nd.push_back( Vector3<T>( 1, 1,-1) );
    Nd.push_back( Vector3<T>(-1, 1,-1) );
    Nd.push_back( Vector3<T>(-1,-1, 1) );
    Nd.push_back( Vector3<T>( 1,-1, 1) );
    Nd.push_back( Vector3<T>( 1, 1, 1) );
    Nd.push_back( Vector3<T>(-1, 1, 1) );
    // 8 undeformed nodes
    Nu.push_back( Vector3<T>(-1,-1,-1) );
    Nu.push_back( Vector3<T>( 1,-1,-1) );
    Nu.push_back( Vector3<T>( 1, 1,-1) );
    Nu.push_back( Vector3<T>(-1, 1,-1) );
    Nu.push_back( Vector3<T>(-1,-1, 1) );
    Nu.push_back( Vector3<T>( 1,-1, 1) );
    Nu.push_back( Vector3<T>( 1, 1, 1) );
    Nu.push_back( Vector3<T>(-1, 1, 1) );
    //*/

    //*
    //-----------------------------------------------------
    // TEST
    // This data gives
    //       |  0    0.5   0     |
    //   J = | -0.5        0     |
    //       |  0    0     0.375 |
    // 8 deformed nodes
    Nd.push_back( Vector3<T>(1,0,0) );
    Nd.push_back( Vector3<T>(1,1,0) );
    Nd.push_back( Vector3<T>(0,1,0) );
    Nd.push_back( Vector3<T>(0,0,0) );
    Nd.push_back( Vector3<T>(1,0,0.75) );
    Nd.push_back( Vector3<T>(1,1,0.75) );
    Nd.push_back( Vector3<T>(0,1,0.75) );
    Nd.push_back( Vector3<T>(0,0,0.75) );
    // 8 undeformed nodes
    Nu.push_back( Vector3<T>(1,0,0) );
    Nu.push_back( Vector3<T>(1,1,0) );
    Nu.push_back( Vector3<T>(0,1,0) );
    Nu.push_back( Vector3<T>(0,0,0) );
    Nu.push_back( Vector3<T>(1,0,0.75) );
    Nu.push_back( Vector3<T>(1,1,0.75) );
    Nu.push_back( Vector3<T>(0,1,0.75) );
    Nu.push_back( Vector3<T>(0,0,0.75) );
    //-----------------------------------------------------
    //*/

    E.push_back( Hexahedron<T>( 
        &Nd,            // the 8 deformed nodes
        &Nu,            // the 8 undeformed nodes
        YoungsModulus,  // Young's modulus
        PoissionsRatio, // Poisson's ratio
        0, 1, 2, 3, 4, 5, 6, 7 )
    );

    rHexFemMesh.SetSizeOfStiffnessMatrix();
    rHexFemMesh.AssembleStiffnessMatrix();

    // Allocate data for simulation
    rHexFemMesh.AllocateDataForSimulation( mass );
}

//-----------------------------------------------------------------------------
template <typename T, typename DATA>
void MeshHexahedraGen<T,DATA>::TwoElementMeshJoinByNode ( MeshHexahedra<T> & rHexFemMesh, T mass, T YoungsModulus, T PoissonsRatio )
{
    // JOIN AT PT 0
    //                 3---------2      ^ y
    //         11-----/| -10    /|      |
    //         /|    / | /|    / |      |
    //        / |   7---------6  |      |
    //      14------| -0- - - | -1      +------>x
    //       |  8- -| /| -9   | /      /
    //       | /    |/ | /    |/      /
    //       |/     4---------5      v z
    //      12---------13

    std::vector< Hexahedron<T> > & E = rHexFemMesh.RetListOfElements();
    std::vector< Vector3<T> > & Nd = rHexFemMesh.RetListOfDeformedMeshNodes();
    std::vector< Vector3<T> > & Nu = rHexFemMesh.RetListOfUndeformedMeshNodes();

    // Clear the mesh
    E.clear();
    Nd.clear();
    Nu.clear();

    // 15 deformed nodes
    // 1st element                          // Pt#
    Nd.push_back( Vector3<T>( 0, 0, 0) );   // 0
    Nd.push_back( Vector3<T>( 1, 0, 0) );   // 1
    Nd.push_back( Vector3<T>( 1, 1, 0) );   // 2
    Nd.push_back( Vector3<T>( 0, 1, 0) );   // 3
    Nd.push_back( Vector3<T>( 0, 0, 1) );   // 4
    Nd.push_back( Vector3<T>( 1, 0, 1) );   // 5
    Nd.push_back( Vector3<T>( 1, 1, 1) );   // 6
    Nd.push_back( Vector3<T>( 0, 1, 1) );   // 7
    // 2nd element
    Nd.push_back( Vector3<T>(-1,-1,-1) );   // 8
    Nd.push_back( Vector3<T>( 0,-1,-1) );   // 9
    Nd.push_back( Vector3<T>( 0, 0,-1) );   // 10
    Nd.push_back( Vector3<T>(-1, 0,-1) );   // 11
    Nd.push_back( Vector3<T>(-1,-1, 0) );   // 12
    Nd.push_back( Vector3<T>( 0,-1, 0) );   // 13
    //Nd.push_back( Vector3<T>( 0, 0, 0) ); // == 0
    Nd.push_back( Vector3<T>(-1, 0, 0) );   // 14

    // 15 undeformed nodes
    // 1st element                          // Pt#
    Nu.push_back( Vector3<T>( 0, 0, 0) );   // 0
    Nu.push_back( Vector3<T>( 1, 0, 0) );   // 1
    Nu.push_back( Vector3<T>( 1, 1, 0) );   // 2
    Nu.push_back( Vector3<T>( 0, 1, 0) );   // 3
    Nu.push_back( Vector3<T>( 0, 0, 1) );   // 4
    Nu.push_back( Vector3<T>( 1, 0, 1) );   // 5
    Nu.push_back( Vector3<T>( 1, 1, 1) );   // 6
    Nu.push_back( Vector3<T>( 0, 1, 1) );   // 7
    // 2nd element
    Nu.push_back( Vector3<T>(-1,-1,-1) );   // 8
    Nu.push_back( Vector3<T>( 0,-1,-1) );   // 9
    Nu.push_back( Vector3<T>( 0, 0,-1) );   // 10
    Nu.push_back( Vector3<T>(-1, 0,-1) );   // 11
    Nu.push_back( Vector3<T>(-1,-1, 0) );   // 12
    Nu.push_back( Vector3<T>( 0,-1, 0) );   // 13
    //Nu.push_back( Vector3<T>( 0, 0, 0) ); // == 0
    Nu.push_back( Vector3<T>(-1, 0, 0) );   // 14

    // 1st element
    E.push_back( Hexahedron<T>( 
        &Nd,            // the 8 deformed nodes
        &Nu,            // the 8 undeformed nodes
        YoungsModulus,  // Young's modulus
        PoissionsRatio, // Poisson's ratio
        0, 1, 2, 3, 4, 5, 6, 7 )
    );
    // 2nd element
    E.push_back( Hexahedron<T>( 
        &Nd,            // the 8 deformed nodes
        &Nu,            // the 8 undeformed nodes
        YoungsModulus,  // Young's modulus
        PoissionsRatio, // Poisson's ratio
        8, 9, 10, 11, 12, 13, 0, 14 )
    );

    rHexFemMesh.SetSizeOfStiffnessMatrix();
    rHexFemMesh.AssembleStiffnessMatrix();

    // Allocate data for simulation
    rHexFemMesh.AllocateDataForSimulation( mass );
}

//-----------------------------------------------------------------------------
template <typename T, typename DATA>
void MeshHexahedraGen<T,DATA>::TwoElementMeshJoinByEdge ( MeshHexahedra<T> & rHexFemMesh, T mass, T YoungsModulus, T PoissonsRatio )
{
    // JOIN AT EDGE 0-4
    //                 3---------2      ^ y
    //                /|        /|      |
    //               / |       / |      |
    //              7---------6  |      |
    //      10------| -0- - - | -1      +------>x
    //      /|      | /|      | /      /
    //     / |      |/ |      |/      /
    //   13---------4---------5      v z
    //    |  8------|--9
    //    | /       | /
    //    |/        |/
    //   11---------12

    std::vector< Hexahedron<T> > & E = rHexFemMesh.RetListOfElements();
    std::vector< Vector3<T> > & Nd = rHexFemMesh.RetListOfDeformedMeshNodes();
    std::vector< Vector3<T> > & Nu = rHexFemMesh.RetListOfUndeformedMeshNodes();

    // Clear the mesh
    E.clear();
    Nd.clear();
    Nu.clear();

    // 15 deformed nodes
    // 1st element                          // Pt#
    Nd.push_back( Vector3<T>( 0, 0, 0) );   // 0
    Nd.push_back( Vector3<T>( 1, 0, 0) );   // 1
    Nd.push_back( Vector3<T>( 1, 1, 0) );   // 2
    Nd.push_back( Vector3<T>( 0, 1, 0) );   // 3
    Nd.push_back( Vector3<T>( 0, 0, 1) );   // 4
    Nd.push_back( Vector3<T>( 1, 0, 1) );   // 5
    Nd.push_back( Vector3<T>( 1, 1, 1) );   // 6
    Nd.push_back( Vector3<T>( 0, 1, 1) );   // 7
    // 2nd element
    Nd.push_back( Vector3<T>(-1,-1, 0) );   // 8
    Nd.push_back( Vector3<T>( 0,-1, 0) );   // 9
    //Nd.push_back( Vector3<T>( 0, 0, 0) ); // == 0
    Nd.push_back( Vector3<T>(-1, 0, 0) );   // 10
    Nd.push_back( Vector3<T>(-1,-1, 1) );   // 11
    Nd.push_back( Vector3<T>( 0,-1, 1) );   // 12
    //Nd.push_back( Vector3<T>( 0, 0, 1) ); // == 4
    Nd.push_back( Vector3<T>(-1, 0, 1) );   // 13

    // 15 undeformed nodes
    // 1st element                          // Pt#
    Nu.push_back( Vector3<T>( 0, 0, 0) );   // 0
    Nu.push_back( Vector3<T>( 1, 0, 0) );   // 1
    Nu.push_back( Vector3<T>( 1, 1, 0) );   // 2
    Nu.push_back( Vector3<T>( 0, 1, 0) );   // 3
    Nu.push_back( Vector3<T>( 0, 0, 1) );   // 4
    Nu.push_back( Vector3<T>( 1, 0, 1) );   // 5
    Nu.push_back( Vector3<T>( 1, 1, 1) );   // 6
    Nu.push_back( Vector3<T>( 0, 1, 1) );   // 7
    // 2nd element
    Nu.push_back( Vector3<T>(-1,-1, 0) );   // 8
    Nu.push_back( Vector3<T>( 0,-1, 0) );   // 9
    //Nu.push_back( Vector3<T>( 0, 0, 0) ); // == 0
    Nu.push_back( Vector3<T>(-1, 0, 0) );   // 10
    Nu.push_back( Vector3<T>(-1,-1, 1) );   // 11
    Nu.push_back( Vector3<T>( 0,-1, 1) );   // 12
    //Nu.push_back( Vector3<T>( 0, 0, 1) ); // == 4
    Nu.push_back( Vector3<T>(-1, 0, 1) );   // 13

    // 1st element
    E.push_back( Hexahedron<T>( 
        &Nd,            // the 8 deformed nodes
        &Nu,            // the 8 undeformed nodes
        YoungsModulus,  // Young's modulus
        PoissionsRatio, // Poisson's ratio
        0, 1, 2, 3, 4, 5, 6, 7 )
    );
    // 2nd element
    E.push_back( Hexahedron<T>( 
        &Nd,            // the 8 deformed nodes
        &Nu,            // the 8 undeformed nodes
        YoungsModulus,  // Young's modulus
        PoissionsRatio, // Poisson's ratio
        8, 9, 0, 10, 11, 12, 4, 13 )
    );

    rHexFemMesh.SetSizeOfStiffnessMatrix();
    rHexFemMesh.AssembleStiffnessMatrix();

    // Allocate data for simulation
    rHexFemMesh.AllocateDataForSimulation( mass );
}

//-----------------------------------------------------------------------------
template <typename T, typename DATA>
void MeshHexahedraGen<T,DATA>::TwoElementMeshJoinByFace ( MeshHexahedra<T> & rHexFemMesh, T mass, T YoungsModulus, T PoissonsRatio )
{
    // JOIN AT FACE 4-5-6-7
    //                 3---------2      ^ y
    //                /|        /|      |
    //               / |       / |      |
    //              7---------6  |      |
    //             /|  0- - -/| -1      +------>x
    //            / | /     / | /      /
    //          11---------10 |/      /
    //           |  4------|--5      v z
    //           | /       | /
    //           |/        |/
    //           8---------9

    std::vector< Hexahedron<T> > & E = rHexFemMesh.RetListOfElements();
    std::vector< Vector3<T> > & Nd = rHexFemMesh.RetListOfDeformedMeshNodes();
    std::vector< Vector3<T> > & Nu = rHexFemMesh.RetListOfUndeformedMeshNodes();

    // Clear the mesh
    E.clear();
    Nd.clear();
    Nu.clear();

    // 15 deformed nodes
    // 1st element                          // Pt#
    Nd.push_back( Vector3<T>( 0, 0, 0) );   // 0
    Nd.push_back( Vector3<T>( 1, 0, 0) );   // 1
    Nd.push_back( Vector3<T>( 1, 1, 0) );   // 2
    Nd.push_back( Vector3<T>( 0, 1, 0) );   // 3
    Nd.push_back( Vector3<T>( 0, 0, 1) );   // 4
    Nd.push_back( Vector3<T>( 1, 0, 1) );   // 5
    Nd.push_back( Vector3<T>( 1, 1, 1) );   // 6
    Nd.push_back( Vector3<T>( 0, 1, 1) );   // 7
    // 2nd element
    Nd.push_back( Vector3<T>( 0, 0, 2) );   // 8
    Nd.push_back( Vector3<T>( 1, 0, 2) );   // 9
    Nd.push_back( Vector3<T>( 1, 1, 2) );   // 10
    Nd.push_back( Vector3<T>( 0, 1, 2) );   // 11

    // 15 undeformed nodes
    // 1st element                          // Pt#
    Nu.push_back( Vector3<T>( 0, 0, 0) );   // 0
    Nu.push_back( Vector3<T>( 1, 0, 0) );   // 1
    Nu.push_back( Vector3<T>( 1, 1, 0) );   // 2
    Nu.push_back( Vector3<T>( 0, 1, 0) );   // 3
    Nu.push_back( Vector3<T>( 0, 0, 1) );   // 4
    Nu.push_back( Vector3<T>( 1, 0, 1) );   // 5
    Nu.push_back( Vector3<T>( 1, 1, 1) );   // 6
    Nu.push_back( Vector3<T>( 0, 1, 1) );   // 7
    // 2nd element
    Nu.push_back( Vector3<T>( 0, 0, 2) );   // 8
    Nu.push_back( Vector3<T>( 1, 0, 2) );   // 9
    Nu.push_back( Vector3<T>( 1, 1, 2) );   // 10
    Nu.push_back( Vector3<T>( 0, 1, 2) );   // 11

    // 1st element
    E.push_back( Hexahedron<T>( 
        &Nd,            // the 8 deformed nodes
        &Nu,            // the 8 undeformed nodes
        YoungsModulus,  // Young's modulus
        PoissionsRatio, // Poisson's ratio
        0, 1, 2, 3, 4, 5, 6, 7 )
    );
    // 2nd element
    E.push_back( Hexahedron<T>( 
        &Nd,            // the 8 deformed nodes
        &Nu,            // the 8 undeformed nodes
        YoungsModulus,  // Young's modulus
        PoissionsRatio, // Poisson's ratio
        4, 5, 6, 7, 8, 9, 10, 11 )
    );

    rHexFemMesh.SetSizeOfStiffnessMatrix();
    rHexFemMesh.AssembleStiffnessMatrix();

    // Allocate data for simulation
    rHexFemMesh.AllocateDataForSimulation( mass );
}



//-----------------------------------------------------------------------------
template <typename T, typename DATA>
void MeshHexahedraGen<T,DATA>::Bar (
    MeshHexahedra<T> & rHexFemMesh,
    unsigned int numOfElements,
    T width,
    T height,
    T depth,
    T mass,
    T YoungsModulus,
    T PoissonsRatio
)
{
    // JOIN AT FACE 4-5-6-7
    //                 3---------2      ^ y
    //                /|        /|      |
    //               / |       / |      |
    //              7---------6  |      |
    //             /|  0- - -/| -1      +------>x
    //            / | /     / | /      /
    //          11---------10 |/      /
    //           |  4------|--5      v z
    //           | /       | /
    //           |/        |/
    //           8---------9

    std::vector< Hexahedron<T> > & E = rHexFemMesh.RetListOfElements();
    std::vector< Vector3<T> > & Nd = rHexFemMesh.RetListOfDeformedMeshNodes();
    std::vector< Vector3<T> > & Nu = rHexFemMesh.RetListOfUndeformedMeshNodes();

    // Clear the mesh
    E.clear();
    Nd.clear();
    Nu.clear();

    T depthStep = depth / numOfElements;
    T w = width, h = height, d = T(0);

    // 1st element (deformed nodes)         // Pt#
    Nd.push_back( Vector3<T>(-w,-h, d) );   // 0
    Nd.push_back( Vector3<T>( w,-h, d) );   // 1
    Nd.push_back( Vector3<T>( w, h, d) );   // 2
    Nd.push_back( Vector3<T>(-w, h, d) );   // 3

    Nu.push_back( Vector3<T>(-w,-h, d) );   // 0
    Nu.push_back( Vector3<T>( w,-h, d) );   // 1
    Nu.push_back( Vector3<T>( w, h, d) );   // 2
    Nu.push_back( Vector3<T>(-w, h, d) );   // 3

    // Next elements
    for ( unsigned int i = 0, n=0; i < numOfElements; ++i, n+=4 ) {
        d += depthStep;

        Nd.push_back( Vector3<T>(-w,-h, d) );   // 4
        Nd.push_back( Vector3<T>( w,-h, d) );   // 5
        Nd.push_back( Vector3<T>( w, h, d) );   // 6
        Nd.push_back( Vector3<T>(-w, h, d) );   // 7

        Nu.push_back( Vector3<T>(-w,-h, d) );   // 4
        Nu.push_back( Vector3<T>( w,-h, d) );   // 5
        Nu.push_back( Vector3<T>( w, h, d) );   // 6
        Nu.push_back( Vector3<T>(-w, h, d) );   // 7

        E.push_back( Hexahedron<T>( 
            &Nd,            // the 8 deformed nodes
            &Nu,            // the 8 undeformed nodes
            YoungsModulus,  // Young's modulus
            PoissonsRatio,  // Poisson's ratio
            n, n+1, n+2, n+3, n+4, n+5, n+6, n+7 )
        );
    }

    rHexFemMesh.SetSizeOfStiffnessMatrix();
    rHexFemMesh.AssembleStiffnessMatrix();

    // Allocate data for simulation
    rHexFemMesh.AllocateDataForSimulation( mass );
}



//-----------------------------------------------------------------------------
template <typename T, typename DATA>
bool MeshHexahedraGen<T,DATA>::MeshFromGridGenerator (
        MeshHexahedra<T> & rHexFemMesh,     
        GridGenerator<T,DATA> & rGridGen,   
        T   mass,           
        T   YoungsModulus,  
        T   PoissonsRatio   
)
{

    if ( !rGridGen.IsGridGenerated() ) {
        return false;
    }

    // Grid layout (right-handed coordinates)
    //      ^ y
    //      |
    //      |
    //      0-----> x       x = column
    //     /                y = row
    //    /                 z = depth
    //   v z
    //
    //  Where X := #columns-1, Y := #rows-1, and Z := #slices-1.
    //
    //                   O-------O-------O---------------O-------O
    //                  /|      /|      /|      ...     /|      /|
    //                 / |(0,Y,0)|     / |      ...    / |(X,Y,0)|
    //                /  |    /  |    /  |      ...   /  |    /  |
    //               /   O---/---O---/---O-----------/---O---/---O
    //              /    |  /    |  /    |      ... /    |  /   /|
    //             /     | /     | /     |      .../     | /   / |
    //            /      |/      |/      |      ../      |/   /  |
    //           /       /-------/-------O-------/-------/---/---O
    //          /       /|      /               HIGHEST / POINT /|
    //         O-------O-------O---------------O-------@   /   / |
    //         |       |       |      ...      | .     |  /   /  |
    //         |(0,Y,Z)|       |      ...      |(X,Y,Z)| /   /   |
    //         |       |       |      ...      | .     |/   /    |
    //         O-------O-------O---------------O-------O   /     |
    //         |       |       |      ...      | .     |  /      |
    //         |       |       |      ...      | .     | /       |
    //         |       |       |      ...      | .     |/        |
    //         O-------O-------O---------------O-------O         |
    //         |         |             .         .     |         |
    //         |         O-------O-----.-O-------------|-O-------0
    //         |        /|       |     . |      ...    | |      /|
    //         |       / |(0,0,0)|     . |      ...    | |(X,0,0)|
    //         |      /  |       |     . |      ...    | |    /  |
    //         |     /   @-------O-----.-O-------------|-O---/---O
    //         | LOWEST / POINT /      ./              |/   /   /
    //         |   /   /       /       .               |   /   /
    //         |  /   /       /       /.              /|  /   /
    //         | /   /       /       / .             / | /   /
    //         |/   /       /       /  .            /  |/   /
    //         O-------O-------O---------------O-------O   /
    //         |  /    |  /    |  /   ...      |  /    |  /
    //         |(0,0,Z)| /     | /    ...      |(X,0,Z)| /
    //         |/      |/      |/     ...      |/      |/
    //         O-------O-------O---------------O-------O
    //
    //  where D := Z;  R := Y;  C := X;
    //
    //  The n-th Z-slice
    //
    //  (nRC+(R-1)+C+0)--(nRC+(R-1)+C+1)--(nRC+(R-1)+C+2)-- ... --(nRC+(R+1)C+C-1)
    //          |                |                |                       |
    //          |                |                |                       |
    //          |                |                |                       |
    //         ...              ...              ...        ...          ...
    //          |                |                |                       |
    //          |                |                |                       |
    //          |                |                |                       |
    //     (nRC+2C+0)-------(nRC+2C+1)-------(nRC+2C+2)---- ... -----(nRC+2C+C-1)
    //          |                |                |                       |
    //          |                |                |                       |
    //          |                |                |                       |
    //      (nRC+C+0)--------(nRC+C+1)--------(nRC+C+2)---- ... ------(nRC+C+C-1)
    //          |                |                |                       |
    //          |                |                |                       |
    //          |                |                |                       |
    //       (nRC+0)----------(nRC+1)----------(nRC+2)----- ... -------(nRC+C-1)
    //
    //
    //
    //                 @---------@      ^ y
    //                /|        /|      |
    //               / | SLICE /A|      |
    //              @---------@  |      |
    //             /|  @- - -/| -@      +------>x
    //            / | SLICE /B| /      /
    //           @---------@  |/      /
    //           |  @------|--@      v z
    //           | /SLICE C| /
    //           |/        |/
    //           @---------@


    std::vector< Hexahedron<T> > & E = rHexFemMesh.RetListOfElements();
    std::vector< Vector3<T> > & Nd = rHexFemMesh.RetListOfDeformedMeshNodes();
    std::vector< Vector3<T> > & Nu = rHexFemMesh.RetListOfUndeformedMeshNodes();

    // Clear the mesh
    E.clear();
    Nd.clear();
    Nu.clear();

    // Z goes first follows by Y then X.
    // So the slice is in Z-plane with Y as rows and X as columns
    T **** const GridVertices = rGridGen.ReturnPtrToVertexPositionData();
    unsigned int z,y,x;
    unsigned int C = rGridGen.GetGridSizeX();
    unsigned int R = rGridGen.GetGridSizeY();
    unsigned int D = rGridGen.GetGridSizeZ();

    // Add mesh vertices from the grid data
    for ( z = 0; z < D; ++z ) {
        for ( y = 0; y < R; ++y ) {
            for ( x = 0; x < C; ++x ) {
                Vector3<T> V( GridVertices[x][y][z][0], GridVertices[x][y][z][1], GridVertices[x][y][z][2] );
                Nd.push_back( V );
                Nu.push_back( V );
            }
        }
    }

    // Add hexahedron element from the mesh vertices
    unsigned int idx[4];
    unsigned int RC  = R*C;
    unsigned int SLICE = 0;
    unsigned int ROW;
    for ( z = 0; z < D-1; ++z ) {
        ROW = 0;
        for ( y = 0; y < R-1; ++y ) {
            for ( x = 0; x < C-1; ++x ) {
                // Pt 1 to 8
                //
                //    idx1 3---------2      ^ y
                //        /|        /|      |
                //       / |       / |      |
                // idx3 7---------6  |      |
                // idx0 |  0- - - | -1      +------>x
                //      | /       | /      /
                //      |/        |/      /
                // idx2 4---------5      v z
                idx[0] = SLICE + ROW + x;
                idx[1] = idx[0] + C;
                idx[2] = idx[0] + RC;
                idx[3] = idx[2] + C;

                E.push_back( Hexahedron<T>( 
                    &Nd,            // the 8 deformed nodes
                    &Nu,            // the 8 undeformed nodes
                    YoungsModulus,  // Young's modulus
                    PoissonsRatio,  // Poisson's ratio
                    idx[0], idx[0]+1, idx[1]+1, idx[1], idx[2], idx[2]+1, idx[3]+1, idx[3] )
                );
            }
            ROW += C;
        }
        SLICE += RC;
    }

    //std::cout << "set size of stiffness matrix" << std::endl;
    rHexFemMesh.SetSizeOfStiffnessMatrix();
    //std::cout << "assemble stiffness matrix" << std::endl;
    rHexFemMesh.AssembleStiffnessMatrix();

    // Allocate data for simulation
    //std::cout << "Allocate data for simulation" << std::endl;
    rHexFemMesh.AllocateDataForSimulation( mass );

    // DEBUG
    std::cout << "Generated Hexahedral FEM mesh has " << rHexFemMesh.GetNumOfNodes() 
        << " nodes and " << rHexFemMesh.GetNumOfElements() << " elements.\n";

    return true;
}



//-----------------------------------------------------------------------------
template <typename T, typename DATA>
bool MeshHexahedraGen<T,DATA>::MeshFromGridGenerator (
        MeshHexahedra<T> & rHexFemMesh,     
        std::vector< HexBarycentricCoordsForVertexRef<T,DATA> > & rBaryCoords, 
        //std::vector< unsigned int > & rElementID, //!< identify the element that the barycentric coordinates are embedded in
        GridGenerator<T,DATA> & rGridGen,   
        T   mass,           
        T   YoungsModulus,  
        T   PoissonsRatio,  
        bool bGenElementsInsideVolume,  
        int *** pElementsGridLocations  
)
{

    if ( !rGridGen.IsGridGenerated() ) {
        return false;
    }

    // Grid layout (right-handed coordinates)
    //      ^ y
    //      |
    //      |
    //      0-----> x       x = column
    //     /                y = row
    //    /                 z = depth
    //   v z
    //
    //  Where X := #columns-1, Y := #rows-1, and Z := #slices-1.
    //
    //                   O-------O-------O---------------O-------O
    //                  /|      /|      /|      ...     /|      /|
    //                 / |(0,Y,0)|     / |      ...    / |(X,Y,0)|
    //                /  |    /  |    /  |      ...   /  |    /  |
    //               /   O---/---O---/---O-----------/---O---/---O
    //              /    |  /    |  /    |      ... /    |  /   /|
    //             /     | /     | /     |      .../     | /   / |
    //            /      |/      |/      |      ../      |/   /  |
    //           /       /-------/-------O-------/-------/---/---O
    //          /       /|      /               HIGHEST / POINT /|
    //         O-------O-------O---------------O-------@   /   / |
    //         |       |       |      ...      | .     |  /   /  |
    //         |(0,Y,Z)|       |      ...      |(X,Y,Z)| /   /   |
    //         |       |       |      ...      | .     |/   /    |
    //         O-------O-------O---------------O-------O   /     |
    //         |       |       |      ...      | .     |  /      |
    //         |       |       |      ...      | .     | /       |
    //         |       |       |      ...      | .     |/        |
    //         O-------O-------O---------------O-------O         |
    //         |         |             .         .     |         |
    //         |         O-------O-----.-O-------------|-O-------0
    //         |        /|       |     . |      ...    | |      /|
    //         |       / |(0,0,0)|     . |      ...    | |(X,0,0)|
    //         |      /  |       |     . |      ...    | |    /  |
    //         |     /   @-------O-----.-O-------------|-O---/---O
    //         | LOWEST / POINT /      ./              |/   /   /
    //         |   /   /       /       .               |   /   /
    //         |  /   /       /       /.              /|  /   /
    //         | /   /       /       / .             / | /   /
    //         |/   /       /       /  .            /  |/   /
    //         O-------O-------O---------------O-------O   /
    //         |  /    |  /    |  /   ...      |  /    |  /
    //         |(0,0,Z)| /     | /    ...      |(X,0,Z)| /
    //         |/      |/      |/     ...      |/      |/
    //         O-------O-------O---------------O-------O
    //
    //  where D := Z;  R := Y;  C := X;
    //
    //  The n-th Z-slice
    //
    //  (nRC+(R-1)+C+0)--(nRC+(R-1)+C+1)--(nRC+(R-1)+C+2)-- ... --(nRC+(R+1)C+C-1)
    //          |                |                |                       |
    //          |                |                |                       |
    //          |                |                |                       |
    //         ...              ...              ...        ...          ...
    //          |                |                |                       |
    //          |                |                |                       |
    //          |                |                |                       |
    //     (nRC+2C+0)-------(nRC+2C+1)-------(nRC+2C+2)---- ... -----(nRC+2C+C-1)
    //          |                |                |                       |
    //          |                |                |                       |
    //          |                |                |                       |
    //      (nRC+C+0)--------(nRC+C+1)--------(nRC+C+2)---- ... ------(nRC+C+C-1)
    //          |                |                |                       |
    //          |                |                |                       |
    //          |                |                |                       |
    //       (nRC+0)----------(nRC+1)----------(nRC+2)----- ... -------(nRC+C-1)
    //
    //
    //
    //                 @---------@      ^ y
    //                /|        /|      |
    //               / | SLICE /A|      |
    //              @---------@  |      |
    //             /|  @- - -/| -@      +------>x
    //            / | SLICE /B| /      /
    //           @---------@  |/      /
    //           |  @------|--@      v z
    //           | /SLICE C| /
    //           |/        |/
    //           @---------@

    std::vector< Hexahedron<T> > & E = rHexFemMesh.RetListOfElements();
    std::vector< Vector3<T> > & Nd = rHexFemMesh.RetListOfDeformedMeshNodes();
    std::vector< Vector3<T> > & Nu = rHexFemMesh.RetListOfUndeformedMeshNodes();

    // Clear the mesh
    E.clear();
    Nd.clear();
    Nu.clear();

    // Z goes first follows by Y then X.
    // So the slice is in Z-plane with Y as rows and X as columns
    T **** const GridVertices = rGridGen.ReturnPtrToVertexPositionData();
    unsigned int z,y,x;
    unsigned int C = rGridGen.GetGridSizeX();
    unsigned int R = rGridGen.GetGridSizeY();
    unsigned int D = rGridGen.GetGridSizeZ();
    unsigned int RC  = R*C;
    unsigned int C_1 = C-1;
    unsigned int R_1 = R-1;
    unsigned int D_1 = D-1;
    unsigned int RC_1 = R_1*C_1;

    unsigned int count = 0;

    // convert std::vector< GenBarycentricCoordsForVertexRef<T,DATA> > 
    // to std::vector< std::vector< GenBarycentricCoordsForVertexRef<T,DATA> > >
    std::vector< std::vector<unsigned int> > TransferCoords( D_1*RC_1 );
    std::vector< HexBarycentricCoordsForVertexRef<T,DATA> >::iterator Coords = rGridGen.GetHexBarycentricCoords().begin();
    while ( Coords != rGridGen.GetHexBarycentricCoords().end() ) {
        Coords->GetGridLocations( x, y, z );
        unsigned int location = z*RC_1 + y*C_1 + x;
        TransferCoords[ location ].push_back( count++ );
        ++Coords;
    }

    //====================================================================
    // Add ONLY mesh vertices THAT EMBEDED IN FEM MESH from the grid data
    //--------------------------------------------------------------------

    // Copy to rBaryCoords (barycentric of polygonal embedded in FEM mesh)
    rBaryCoords = rGridGen.GetHexBarycentricCoords();

    // Allocate temp data to store the valid FEM mesh vertices and elements
    struct TDAT {
        TDAT() : validvertex(false), validelement(false) {}
        bool validvertex;
        unsigned int vertexid;
        bool validelement;
    };
    TDAT *** LOC = new TDAT**[D];
    unsigned int DRC = D*R*C;
    for ( z = 0; z < D; ++z ) {
        LOC[z] = new TDAT*[R];
        for ( y = 0; y < R; ++y ) {
            LOC[z][y] = new TDAT[C];
            // Initialize vertexid to DRC --> same as mark it as the vertex is not created yet
            for ( unsigned int x = 0; x < C; ++x ) {
                LOC[z][y][x].vertexid = DRC;
            }
        }
    }

    // (SURFACE) MESH DATA
    // Mark the valid vertices (location) in temp data based on the barycentric coordinates stored in rGridGen
    Coords = rGridGen.GetHexBarycentricCoords().begin();
    while ( Coords != rGridGen.GetHexBarycentricCoords().end() ) {
        Coords->GetGridLocations( x, y, z );
        // One valid element
        LOC[z][y][x].validelement = true;
        // Eight valid vertices (0 as valid; ULONG_MAX as notvalid)
        LOC[z  ][y  ][x  ].validvertex = true;
        LOC[z  ][y  ][x+1].validvertex = true;
        LOC[z  ][y+1][x+1].validvertex = true;
        LOC[z  ][y+1][x  ].validvertex = true;
        LOC[z+1][y  ][x  ].validvertex = true;
        LOC[z+1][y  ][x+1].validvertex = true;
        LOC[z+1][y+1][x+1].validvertex = true;
        LOC[z+1][y+1][x  ].validvertex = true;
        ++Coords;
    }

    // (VOLUME) GRID DATA
    if ( bGenElementsInsideVolume ) {
        // Mark the valid vertices (location) in temp data based on the grid vertex data stored in rGridGen
        GridGenerator<T,DATA>::VertexFlag *** const vertexFlag = rGridGen.ReturnPtrToVertexFlagData();
        // Mark the valid vertices (location) in temp data
        for ( x = 0; x < C; ++x ) {
            for ( y = 0; y < R; ++y ) {
                for ( z = 0; z < D; ++z ) {
                    if (    vertexFlag[x][y][z] == GridGenerator<T,DATA>::RIGHT_OUTSIDE_BOUNDARY
                        ||  vertexFlag[x][y][z] == GridGenerator<T,DATA>::ON_BOUNDARY
                        ||  vertexFlag[x][y][z] == GridGenerator<T,DATA>::RIGHT_INSIDE_BOUNDARY
                        ||  vertexFlag[x][y][z] == GridGenerator<T,DATA>::INSIDE_MODEL )
                    {
                        LOC[z][y][x].validvertex = true;
                    }
                }
            }
        }
        // Mark the valid elements (location) in temp data
        for ( x = 0; x < C_1; ++x ) {
            for ( y = 0; y < R_1; ++y ) {
                for ( z = 0; z < D_1; ++z ) {
                    if (    LOC[z  ][y  ][x  ].validvertex
                        &&  LOC[z  ][y  ][x+1].validvertex
                        &&  LOC[z  ][y+1][x+1].validvertex
                        &&  LOC[z  ][y+1][x  ].validvertex
                        &&  LOC[z+1][y  ][x  ].validvertex
                        &&  LOC[z+1][y  ][x+1].validvertex
                        &&  LOC[z+1][y+1][x+1].validvertex
                        &&  LOC[z+1][y+1][x  ].validvertex )
                    {
                        LOC[z][y][x].validelement = true;
                    }
                }
            }
        }
    }

    // Add the valid vertices to the FEM mesh data for vertices
    count = 0;
    for ( z = 0; z < D_1; ++z ) {
        for ( y = 0; y < R_1; ++y ) {
            for ( x = 0; x < C_1; ++x ) {
                if ( LOC[z][y][x].validelement ) {
                    // v0
                    if ( LOC[z][y][x].validvertex && LOC[z][y][x].vertexid == DRC ) {
                        Vector3<T> V( GridVertices[x][y][z][0], GridVertices[x][y][z][1], GridVertices[x][y][z][2] );
                        Nd.push_back( V );
                        Nu.push_back( V );
                        LOC[z][y][x].vertexid = count++;
                    }
                    // v1
                    if ( LOC[z][y][x+1].validvertex && LOC[z][y][x+1].vertexid == DRC ) {
                        Vector3<T> V( GridVertices[x+1][y][z][0], GridVertices[x+1][y][z][1], GridVertices[x+1][y][z][2] );
                        Nd.push_back( V );
                        Nu.push_back( V );
                        LOC[z][y][x+1].vertexid = count++;
                    }
                    // v2
                    if ( LOC[z][y+1][x+1].validvertex && LOC[z][y+1][x+1].vertexid == DRC ) {
                        Vector3<T> V( GridVertices[x+1][y+1][z][0], GridVertices[x+1][y+1][z][1], GridVertices[x+1][y+1][z][2] );
                        Nd.push_back( V );
                        Nu.push_back( V );
                        LOC[z][y+1][x+1].vertexid = count++;
                    }
                    // v3
                    if ( LOC[z][y+1][x].validvertex && LOC[z][y+1][x].vertexid == DRC ) {
                        Vector3<T> V( GridVertices[x][y+1][z][0], GridVertices[x][y+1][z][1], GridVertices[x][y+1][z][2] );
                        Nd.push_back( V );
                        Nu.push_back( V );
                        LOC[z][y+1][x].vertexid = count++;
                    }
                    // v4
                    if ( LOC[z+1][y][x].validvertex && LOC[z+1][y][x].vertexid == DRC ) {
                        Vector3<T> V( GridVertices[x][y][z+1][0], GridVertices[x][y][z+1][1], GridVertices[x][y][z+1][2] );
                        Nd.push_back( V );
                        Nu.push_back( V );
                        LOC[z+1][y][x].vertexid = count++;
                    }
                    // v5
                    if ( LOC[z+1][y][x+1].validvertex && LOC[z+1][y][x+1].vertexid == DRC  ) {
                        Vector3<T> V( GridVertices[x+1][y][z+1][0], GridVertices[x+1][y][z+1][1], GridVertices[x+1][y][z+1][2] );
                        Nd.push_back( V );
                        Nu.push_back( V );
                        LOC[z+1][y][x+1].vertexid = count++;
                    }
                    // v6
                    if ( LOC[z+1][y+1][x+1].validvertex && LOC[z+1][y+1][x+1].vertexid == DRC  ) {
                        Vector3<T> V( GridVertices[x+1][y+1][z+1][0], GridVertices[x+1][y+1][z+1][1], GridVertices[x+1][y+1][z+1][2] );
                        Nd.push_back( V );
                        Nu.push_back( V );
                        LOC[z+1][y+1][x+1].vertexid = count++;
                    }
                    // v7
                    if ( LOC[z+1][y+1][x].validvertex && LOC[z+1][y+1][x].vertexid == DRC  ) {
                        Vector3<T> V( GridVertices[x][y+1][z+1][0], GridVertices[x][y+1][z+1][1], GridVertices[x][y+1][z+1][2] );
                        Nd.push_back( V );
                        Nu.push_back( V );
                        LOC[z+1][y+1][x].vertexid = count++;
                    }
                }
            }
        }
    }

    // Add hexahedron element from the mesh vertices
    count = 0;
    //rElementID.resize( rBaryCoords.size() );

    for ( z = 0; z < D_1; ++z ) {
        for ( y = 0; y < R_1; ++y ) {
            for ( x = 0; x < C_1; ++x ) {
                if ( LOC[z][y][x].validelement ) {
                    // Pt 0 to 7
                    //
                    //         3---------2      ^ y
                    //        /|        /|      |
                    //       / |       / |      |
                    //      7---------6  |      |
                    //      |  0- - - | -1      +------>x
                    //      | /       | /      /
                    //      |/        |/      /
                    //      4---------5      v z
                    E.push_back( Hexahedron<T>( 
                        &Nd,            // the 8 deformed nodes
                        &Nu,            // the 8 undeformed nodes
                        YoungsModulus,  // Young's modulus
                        PoissonsRatio,  // Poisson's ratio
                        LOC[z  ][y  ][x  ].vertexid, 
                        LOC[z  ][y  ][x+1].vertexid, 
                        LOC[z  ][y+1][x+1].vertexid, 
                        LOC[z  ][y+1][x  ].vertexid, 
                        LOC[z+1][y  ][x  ].vertexid, 
                        LOC[z+1][y  ][x+1].vertexid, 
                        LOC[z+1][y+1][x+1].vertexid, 
                        LOC[z+1][y+1][x  ].vertexid )
                    );
                    // Record the element id for identifying barycentric coordinate locations
                    unsigned int loc = z*RC_1 + y*C_1 + x;
                    for ( unsigned int i = 0; i < TransferCoords[loc].size(); ++i ) {
                        //rElementID[ TransferCoords[loc][i] ] = count;
                        rBaryCoords[ TransferCoords[loc][i] ].SetElementID( count );
                    }

                    // Store the element ID in the grid location xyz -- width,height,depth
                    if ( pElementsGridLocations ) {
                        pElementsGridLocations[x][y][z] = count;
                    }

                    ++count;

                }
            }
        }
    }

    // Set pointers to hexahedra
    // MUST BE DONE AFTER ALL HEXAHEDRAL ELEMENTS ARE GENERATED!
    // OTHERWISE THE POINTERS WILL NOT BE VALID DUE TO
    // THE AUTOMATICAL RESIZING OF STD::VECTOR!
    for ( unsigned int i = 0; i < rBaryCoords.size(); ++i ) {
        rBaryCoords[i].SetPtrToHexahedron( &E[ rBaryCoords[i].GetElementID() ] );
    }

    // Set the barycentric pointer for each vertex
    // MUST BE DONE AFTER ALL BARY COORDS ARE GENERATED!
    // OTHERWISE THE POINTERS WILL NOT BE VALID DUE TO
    // THE AUTOMATICAL RESIZING OF STD::VECTOR!
    count = 0;
    TAPs::OpenGL::HalfEdgeModel<T> * pHEModel = dynamic_cast< TAPs::OpenGL::HalfEdgeModel<T> * >( rGridGen.GetModelUseToGenGrid() );
    HEVertex<T> * pVertex = pHEModel->GetVertexList()->Head();
    while ( pVertex != NULL ) {
        pVertex->SetBarycentricPtr( &rBaryCoords[count] );
        //std::cout << count << "\tTet\t" << pVertex << "\t" << pVertex->GetBarycentricPtr() << "\t" << *(pVertex->GetBarycentricPtr()) << "\n";
        pVertex = pVertex->Next();
        ++count;
    }

    // Delete temp data
    for ( z = 0; z < D; ++z ) {
        for ( y = 0; y < R; ++y ) {
            delete [] LOC[z][y];
        }
        delete [] LOC[z];
    }
    delete [] LOC;

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

    rHexFemMesh.SetSizeOfStiffnessMatrix();
    rHexFemMesh.AssembleStiffnessMatrix();

    // Allocate data for simulation
    rHexFemMesh.AllocateDataForSimulation( mass );

    // DEBUG
    std::cout << "Generated Hexahedral FEM mesh has " << rHexFemMesh.GetNumOfNodes() 
        << " nodes and " << rHexFemMesh.GetNumOfElements() << " elements.\n";
    std::cout << std::flush;

    return true;
}


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