TAPsDeformMeshWithCut.hpp

Go to the documentation of this file.

/******************************************************************************
TAPsDeformMeshWithCut.hpp

DeformMeshWithCut class is a class for a deformable mesh modeled by a network of 
springs.

SUKITTI PUNAK   (02/01/2005)
UPDATE          (02/01/2005)
*******************************************************************************
+====================+
| DESCRIPTION / NOTE |
+====================+
The model is a mesh of m-by-n virtual masses, which each mass linked to 
its neighbors by massless springs of natural length greater than zeroes.
The linkage between neighbors comprised of three kinds of linking springs:
  -- Structural springs linking masses [m,n] and [m+1,n], and 
                                masses [m,n] and [m,n+1]
  -- Shear springs linking masses [m,n] and [m+1,n+1], and 
                           masses [m+1,n] and [m,n+1]
  -- Flexion springs linking masses [m,n] and [m+2,n], and 
                             masses [m,n] and [m,n+2]

  (0,0)---------(0,1)---------(0,2)
    |  \       /  |  \       /  |
    |    \   /    |    \   /    |
    |      +      |      +      |
    |    /   \    |    /   \    |
    |  /       \  |  /       \  |
  (1,0)---------(1,1)---------(1,2)
    |  \       /  |  \       /  |
    |    \   /    |    \   /    |
    |      +      |      +      |
    |    /   \    |    /   \    |
    |  /       \  |  /       \  |
  (2,0)---------(2,1)---------(2,2)

******************************************************************************/
#ifndef TAPs_DEFORM_MESH_WITH_CUT_HPP
#define TAPs_DEFORM_MESH_WITH_CUT_HPP

#include "../Physics/TAPsPhysics.hpp"   // for gravitational constant
#include "../Physics/TAPsSpring.hpp"    // SPT_Spring.h all ready include SPT_Particle.h
#include "../OpenGL/TAPsOpenGLSupport.hpp"

BEGIN_NAMESPACE_TAPs
//=============================================================================
template <typename T>
class DeformMeshWithCut : public /*virtual*/ OpenGL::OpenGLSupport {
//=============================================================================
    //-------------------------------------------------------------------------
    // put it through ostream
    friend std::ostream &operator<<( std::ostream &output, const DeformMeshWithCut<T> &Obj )
    {
        output  << "DeformMeshWithCut ==> "
                << "\n  number of rows:       " << Obj.m_iNoRows
                << "\n  number of columns:    " << Obj.m_iNoCols
                << "\n";
        return output;
    }

//-------------------------------------------------------------------------------------------------
private:
    // Data Members     ============================================================================
    int m_iNoRows;  // number of rows in y-direction    (y values)
    int m_iNoCols;  // number of columns in x-direction (x values)
    // Mass Points
    Particle<T> ***m_pcParticle;
    // object counters
    int m, n;   // m = row number; n = col number
//-------------------------------------------------------------------------------------------------
private:
    // Helper Functions     ========================================================================
    void CreateParticles( Vector3<T> const &topLeftPos, Vector3<T> const &bottomRightPos, T mass );
    void DeleteParticles();
//-------------------------------------------------------------------------------------------------
public:
    // Member Functions     ========================================================================
    // Constructor -------------------------------------------------------------
    DeformMeshWithCut(    
              int rows,         // number of rows (on world-y-axis)
              int cols,         // number of columns (on world-x-axis)
              const Vector3<T> &topLeftPos,         // start position (should be top left cornor)
              const Vector3<T> &bottomRightPos,     // end position (should be bottom right cornor)
              T mass,                                   // the whole mass
              BitmapImage * = NULL );       // Image Pointer
    //-------------------------------------------------------------------------
    ~DeformMeshWithCut();                           // destructor
    //-------------------------------------------------------------------------
    // CalParticleNormals( int i ) should runs faster than CalParticleNormals(),
    // since it calculates the normals more direct.
    void CalParticleNormals();          // calculate particle normals
    void CalParticleNormals( int i );   // calculate particle normals, int i is dummy
    void AddNormal( Vector3<T> & Dest, 
                    Vector3<T> const & A, Vector3<T> const & B, Vector3<T> const & C );


    // Calculate Forces --------------------------------------------------------
    /*
    void ClearForce();                      // clear force of particles to zeros
    void ForceDueToGravity( T g = SP_Physics::g );  // calculate (weight) force due to (acceleration of) gravity
    void ForceDueToViscousDamping();        // calculate (a loss) force due to viscous damping
    void ForceDueToViscousFluidMoving();    // calculate (wind)   force due to viscous fluid moving
    */
    // Calculate All Forces Acting on Particles
//  void CalForces( Vector3<T> vWind = Vector3<T>(0,0,0), T g = Physics_SI::K::g );

    // Display Routine(s) ------------------------------------------------------
    void DisplayGL( OpenGL::Enum::DrawMode DM = OpenGL::Enum::POLYGON );        // display cloth
    void DisplayGL_ParticleNormals();                       // display particle normals
    //-------------------------------------------------------------------------
    // Get/Set and Animation / Simulation Fns
    Vector3<T> GetPositionOfPointNumber( int row, int col ) const {
        assert( 0 <= row && row < m_iNoRows );
        assert( 0 <= col && col < m_iNoCols );
        return m_pcParticle[row][col]->GetPosition();
    }
    void SetPositionOfPointNumber( int row, int col, Vector3<T> & newPosition ) {
        assert( 0 <= row && row < m_iNoRows );
        assert( 0 <= col && col < m_iNoCols );
        //-----------------------------
        m_pcParticle[row][col]->SetPosition( newPosition );
    }
    void MovePositionOfPointNumber( int row, int col, Vector3<T> & addPosition ) {
        assert( 0 <= row && row < m_iNoRows );
        assert( 0 <= col && col < m_iNoCols );
        //-----------------------------
        m_pcParticle[row][col]->SetPosition( 
            m_pcParticle[row][col]->GetPosition() + addPosition );
    }
    //-------------------------------------------------------------------------
protected:
    //-------------------------------------------------------------------------
    // Pure Virtual Fns
    // Helper Fn
    // void DrawGL()
    virtual void DrawGL ( GLenum ) {};
    //-------------------------------------------------------------------------
public:
    // Get/Set Functions
    void GetNumberOfRowsAndColumns ( int & rows, int & cols ) const
    {   rows = m_iNoRows;   cols = m_iNoCols;   }
    bool GetParticleFixStatus( int r, int c ) const     { return m_pcParticle[r][c]->m_bFixed; }
    void SetParticleFixStatus( int r, int c, bool b )   { m_pcParticle[r][c]->m_bFixed = b; }

    // TEST
    void StepSimulation( T dt ) {
        T divider = 10.0;
        //CalForces( Vector3<T>( (rand()%1000)/divider, (rand()%10)/divider, (rand()%10)/divider ) );
        // Explicit Euler's Method for Integration for ODE Solver
        for ( m = 0; m < m_iNoRows; m++ ) {
            for ( n = 0; n < m_iNoCols; n++ ) {
                m_pcParticle[m][n]->m_vtAccel  = m_pcParticle[m][n]->m_vtForce / m_pcParticle[m][n]->m_tMass;
                m_pcParticle[m][n]->m_vtVeloc += m_pcParticle[m][n]->m_vtAccel * dt;
                m_pcParticle[m][n]->m_vtPosit += m_pcParticle[m][n]->m_vtVeloc * dt;
            }
        }
    }

    //===============================================================
    // DeformUpdate
    //          North
    //            |
    //     West---C---East
    //            |
    //          South
    void DeformUpdate() {
        //static T testValue = 1.6; // Good Number
        static T testValue = 4;
        static T damp  = testValue * testValue;
        static T freq1 = testValue * testValue;
        static T freq2 = testValue / 5;

        T x, y, z;
        T vx, vy, vz;
        T ax, ay, az;               // change in velocity -- i.e. acceleration * time unit
        Vector3<T> nWest, nEast, nNorth, nSouth;    // locations of neighbors
        Vector3<T> nNorthWest, nNorthEast, nSouthWest, nSouthEast;
        Vector3<T> neighborsAvg;
        //------------------------------------------------------------
        // Middle Vertices
        int iHalfNoCols = m_iNoCols / 2;
        for ( m = 1; m < m_iNoRows-1; ++m ) {
            for ( n = 1; n < m_iNoCols-1; ++n ) {       // Full plane
            //for ( n = 1; n < iHalfNoCols-1; ++n ) {       // Half plane
                if ( m_pcParticle[m][n]->GetFixStatus() )   continue;
                //-------------------------------
                // Get neighbors' positions
                nWest   = m_pcParticle[m][n-1]->GetPosition();
                nEast   = m_pcParticle[m][n+1]->GetPosition();
                nSouth  = m_pcParticle[m-1][n]->GetPosition();
                nNorth  = m_pcParticle[m+1][n]->GetPosition();
                //-------------------------------
                nNorthWest  = m_pcParticle[m+1][n-1]->GetPosition();
                nNorthEast  = m_pcParticle[m+1][n+1]->GetPosition();
                nSouthWest  = m_pcParticle[m-1][n-1]->GetPosition();
                nSouthEast  = m_pcParticle[m-1][n+1]->GetPosition();
                //-------------------------------
                //neighborsAvg = (nWest + nEast + nSouth + nNorth ) / 4.0;
                neighborsAvg = (nWest + nEast + nSouth + nNorth + nNorthWest + nNorthEast + nSouthWest + nSouthEast) / 8.0;
                m_pcParticle[m][n]->GetPosition().GetXYZ( x, y, z );
                m_pcParticle[m][n]->GetVelocity().GetXYZ( vx, vy, vz );

                //m_pcParticle[m][n]->SetPosition( m_pcParticle[m][n]->GetPosition() + Vector3<T>(0.01,0.01,0.01) );
                ax = 0.1 * ( 
                            freq1 * (neighborsAvg.GetX() - x) / m_pcParticle[m][n]->GetMass()
                            + freq2 * (neighborsAvg.GetX()*neighborsAvg.GetX() - x*x) 
                          //+ freq2 * (n * 1.0/*size/meshSize = 1.0*/ - x) // Not working for 3D
                          - damp * vx
                );
                ay = 0.1 * ( 
                            freq1 * (neighborsAvg.GetY() - y) / m_pcParticle[m][n]->GetMass()
                            + freq2 * (neighborsAvg.GetY()*neighborsAvg.GetY() - y*y) 
                          //+ freq2 * (m * 1.0/*size/meshSize = 1.0*/ - y) // Not working for 3D
                          - damp * vy
                );
                az = 0.1 * ( 
                            freq1 * (neighborsAvg.GetZ() - z) / m_pcParticle[m][n]->GetMass()
                            + freq2 * (neighborsAvg.GetZ()*neighborsAvg.GetZ() - z*z) 
                          //+ freq2 * (m * 1.0/*size/meshSize = 1.0*/ - z) // Not working for 3D
                          - damp * vz
                );
                vx += ax;
                vy += ay;
                vz += az;
                m_pcParticle[m][n]->SetVelocity( vx, vy, vz );
                x += vx;
                y += vy;
                z += vz;
                m_pcParticle[m][n]->SetPosition( x, y, z );
            }
            /*
            for ( n = iHalfNoCols+2; n < m_iNoRows-1; ++n ) {
                if ( m_pcParticle[m][n]->GetFixStatus() )   continue;
                //-------------------------------
                // Get neighbors' positions
                nWest   = m_pcParticle[m][n-1]->GetPosition();
                nEast   = m_pcParticle[m][n+1]->GetPosition();
                nSouth  = m_pcParticle[m-1][n]->GetPosition();
                nNorth  = m_pcParticle[m+1][n]->GetPosition();
                //-------------------------------
                nNorthWest  = m_pcParticle[m+1][n-1]->GetPosition();
                nNorthEast  = m_pcParticle[m+1][n+1]->GetPosition();
                nSouthWest  = m_pcParticle[m-1][n-1]->GetPosition();
                nSouthEast  = m_pcParticle[m-1][n+1]->GetPosition();
                //-------------------------------
                //neighborsAvg = (nWest + nEast + nSouth + nNorth ) / 4.0;
                neighborsAvg = (nWest + nEast + nSouth + nNorth + nNorthWest + nNorthEast + nSouthWest + nSouthEast) / 8.0;
                m_pcParticle[m][n]->GetPosition().GetXYZ( x, y, z );
                m_pcParticle[m][n]->GetVelocity().GetXYZ( vx, vy, vz );

                //m_pcParticle[m][n]->SetPosition( m_pcParticle[m][n]->GetPosition() + Vector3<T>(0.01,0.01,0.01) );
                ax = 0.1 * ( 
                            freq1 * (neighborsAvg.GetX() - x) / m_pcParticle[m][n]->GetMass()
                            + freq2 * (neighborsAvg.GetX()*neighborsAvg.GetX() - x*x) 
                          - damp * vx
                );
                ay = 0.1 * ( 
                            freq1 * (neighborsAvg.GetY() - y) / m_pcParticle[m][n]->GetMass()
                            + freq2 * (neighborsAvg.GetY()*neighborsAvg.GetY() - y*y) 
                          - damp * vy
                );
                az = 0.1 * ( 
                            freq1 * (neighborsAvg.GetZ() - z) / m_pcParticle[m][n]->GetMass()
                            + freq2 * (neighborsAvg.GetZ()*neighborsAvg.GetZ() - z*z) 
                          - damp * vz
                );
                vx += ax;
                vy += ay;
                vz += az;
                m_pcParticle[m][n]->SetVelocity( vx, vy, vz );
                x += vx;
                y += vy;
                z += vz;
                m_pcParticle[m][n]->SetPosition( x, y, z );
            }
            //*/
        }



        //{ m = 0;
        //  for ( n = 1; n < m_iNoCols-1; ++n ) {
        //      if ( m_pcParticle[m][n]->GetFixStatus() )   continue;
        //      //-------------------------------
        //      // Get neighbors' positions
        //      nWest   = m_pcParticle[m][n-1]->GetPosition();
        //      nEast   = m_pcParticle[m][n+1]->GetPosition();
        //      nNorth  = m_pcParticle[m+1][n]->GetPosition();
        //      //-------------------------------
        //      nNorthWest  = m_pcParticle[m+1][n-1]->GetPosition();
        //      nNorthEast  = m_pcParticle[m+1][n+1]->GetPosition();
        //      //-------------------------------
        //      //neighborsAvg = (nWest + nEast + nSouth + nNorth ) / 4.0;
        //      neighborsAvg = (nWest + nEast + nNorth + nNorthWest + nNorthEast) / 5.0;
        //      m_pcParticle[m][n]->GetPosition().GetXYZ( x, y, z );
        //      m_pcParticle[m][n]->GetVelocity().GetXYZ( vx, vy, vz );

        //      //m_pcParticle[m][n]->SetPosition( m_pcParticle[m][n]->GetPosition() + Vector3<T>(0.01,0.01,0.01) );
        //      ax = 0.1 * ( 
        //                  freq1 * (neighborsAvg.GetX() - x) / m_pcParticle[m][n]->GetMass()
        //                  + freq2 * (neighborsAvg.GetX()*neighborsAvg.GetX() - x*x) 
        //                //+ freq2 * (n * 1.0/*size/meshSize = 1.0*/ - x) // Not working for 3D
        //                - damp * vx
        //      );
        //      ay = 0.1 * ( 
        //                  freq1 * (neighborsAvg.GetY() - y) / m_pcParticle[m][n]->GetMass()
        //                  + freq2 * (neighborsAvg.GetY()*neighborsAvg.GetY() - y*y) 
        //                //+ freq2 * (m * 1.0/*size/meshSize = 1.0*/ - y) // Not working for 3D
        //                - damp * vy
        //      );
        //      az = 0.1 * ( 
        //                  freq1 * (neighborsAvg.GetZ() - z) / m_pcParticle[m][n]->GetMass()
        //                  + freq2 * (neighborsAvg.GetZ()*neighborsAvg.GetZ() - z*z) 
        //                //+ freq2 * (m * 1.0/*size/meshSize = 1.0*/ - z) // Not working for 3D
        //                - damp * vz
        //      );
        //      vx += ax;
        //      vy += ay;
        //      vz += az;
        //      m_pcParticle[m][n]->SetVelocity( vx, vy, vz );
        //      x += vx;
        //      y += vy;
        //      z += vz;
        //      m_pcParticle[m][n]->SetPosition( x, y, z );
        //  }
        //}
        //{ m = m_iNoRows-1;
        //  for ( n = 1; n < m_iNoCols-1; ++n ) {
        //      if ( m_pcParticle[m][n]->GetFixStatus() )   continue;
        //      //-------------------------------
        //      // Get neighbors' positions
        //      nWest   = m_pcParticle[m][n-1]->GetPosition();
        //      nEast   = m_pcParticle[m][n+1]->GetPosition();
        //      nSouth  = m_pcParticle[m-1][n]->GetPosition();
        //      //-------------------------------
        //      nSouthWest  = m_pcParticle[m-1][n-1]->GetPosition();
        //      nSouthEast  = m_pcParticle[m-1][n+1]->GetPosition();
        //      //-------------------------------
        //      //neighborsAvg = (nWest + nEast + nSouth + nNorth ) / 4.0;
        //      neighborsAvg = (nWest + nEast + nSouth + nSouthWest + nSouthEast) / 5.0;
        //      m_pcParticle[m][n]->GetPosition().GetXYZ( x, y, z );
        //      m_pcParticle[m][n]->GetVelocity().GetXYZ( vx, vy, vz );

        //      //m_pcParticle[m][n]->SetPosition( m_pcParticle[m][n]->GetPosition() + Vector3<T>(0.01,0.01,0.01) );
        //      ax = 0.1 * ( 
        //                  freq1 * (neighborsAvg.GetX() - x) / m_pcParticle[m][n]->GetMass()
        //                  + freq2 * (neighborsAvg.GetX()*neighborsAvg.GetX() - x*x) 
        //                //+ freq2 * (n * 1.0/*size/meshSize = 1.0*/ - x) // Not working for 3D
        //                - damp * vx
        //      );
        //      ay = 0.1 * ( 
        //                  freq1 * (neighborsAvg.GetY() - y) / m_pcParticle[m][n]->GetMass()
        //                  + freq2 * (neighborsAvg.GetY()*neighborsAvg.GetY() - y*y) 
        //                //+ freq2 * (m * 1.0/*size/meshSize = 1.0*/ - y) // Not working for 3D
        //                - damp * vy
        //      );
        //      az = 0.1 * ( 
        //                  freq1 * (neighborsAvg.GetZ() - z) / m_pcParticle[m][n]->GetMass()
        //                  + freq2 * (neighborsAvg.GetZ()*neighborsAvg.GetZ() - z*z) 
        //                //+ freq2 * (m * 1.0/*size/meshSize = 1.0*/ - z) // Not working for 3D
        //                - damp * vz
        //      );
        //      vx += ax;
        //      vy += ay;
        //      vz += az;
        //      m_pcParticle[m][n]->SetVelocity( vx, vy, vz );
        //      x += vx;
        //      y += vy;
        //      z += vz;
        //      m_pcParticle[m][n]->SetPosition( x, y, z );
        //  }
        //}
        //{ n = 0;
        //  for ( m = 1; m < m_iNoRows-1; ++m ) {
        //      if ( m_pcParticle[m][n]->GetFixStatus() )   continue;
        //      //-------------------------------
        //      // Get neighbors' positions
        //      nEast   = m_pcParticle[m][n+1]->GetPosition();
        //      nSouth  = m_pcParticle[m-1][n]->GetPosition();
        //      nNorth  = m_pcParticle[m+1][n]->GetPosition();
        //      //-------------------------------
        //      nNorthEast  = m_pcParticle[m+1][n+1]->GetPosition();
        //      nSouthEast  = m_pcParticle[m-1][n+1]->GetPosition();
        //      //-------------------------------
        //      //neighborsAvg = (nWest + nEast + nSouth + nNorth ) / 4.0;
        //      neighborsAvg = (nEast + nSouth + nNorth + nNorthEast + nSouthEast) / 5.0;
        //      m_pcParticle[m][n]->GetPosition().GetXYZ( x, y, z );
        //      m_pcParticle[m][n]->GetVelocity().GetXYZ( vx, vy, vz );

        //      //m_pcParticle[m][n]->SetPosition( m_pcParticle[m][n]->GetPosition() + Vector3<T>(0.01,0.01,0.01) );
        //      ax = 0.1 * ( 
        //                  freq1 * (neighborsAvg.GetX() - x) / m_pcParticle[m][n]->GetMass()
        //                  + freq2 * (neighborsAvg.GetX()*neighborsAvg.GetX() - x*x) 
        //                //+ freq2 * (n * 1.0/*size/meshSize = 1.0*/ - x) // Not working for 3D
        //                - damp * vx
        //      );
        //      ay = 0.1 * ( 
        //                  freq1 * (neighborsAvg.GetY() - y) / m_pcParticle[m][n]->GetMass()
        //                  + freq2 * (neighborsAvg.GetY()*neighborsAvg.GetY() - y*y) 
        //                //+ freq2 * (m * 1.0/*size/meshSize = 1.0*/ - y) // Not working for 3D
        //                - damp * vy
        //      );
        //      az = 0.1 * ( 
        //                  freq1 * (neighborsAvg.GetZ() - z) / m_pcParticle[m][n]->GetMass()
        //                  + freq2 * (neighborsAvg.GetZ()*neighborsAvg.GetZ() - z*z) 
        //                //+ freq2 * (m * 1.0/*size/meshSize = 1.0*/ - z) // Not working for 3D
        //                - damp * vz
        //      );
        //      vx += ax;
        //      vy += ay;
        //      vz += az;
        //      m_pcParticle[m][n]->SetVelocity( vx, vy, vz );
        //      x += vx;
        //      y += vy;
        //      z += vz;
        //      m_pcParticle[m][n]->SetPosition( x, y, z );
        //  }
        //}
        //{ n = m_iNoCols-1;
        //  for ( m = 1; m < m_iNoRows-1; ++m ) {
        //      if ( m_pcParticle[m][n]->GetFixStatus() )   continue;
        //      //-------------------------------
        //      // Get neighbors' positions
        //      nWest   = m_pcParticle[m][n-1]->GetPosition();
        //      nSouth  = m_pcParticle[m-1][n]->GetPosition();
        //      nNorth  = m_pcParticle[m+1][n]->GetPosition();
        //      //-------------------------------
        //      nNorthWest  = m_pcParticle[m+1][n-1]->GetPosition();
        //      nSouthWest  = m_pcParticle[m-1][n-1]->GetPosition();
        //      //-------------------------------
        //      //neighborsAvg = (nWest + nEast + nSouth + nNorth ) / 4.0;
        //      neighborsAvg = (nWest + nSouth + nNorth + nNorthWest + nSouthWest) / 5.0;
        //      m_pcParticle[m][n]->GetPosition().GetXYZ( x, y, z );
        //      m_pcParticle[m][n]->GetVelocity().GetXYZ( vx, vy, vz );

        //      //m_pcParticle[m][n]->SetPosition( m_pcParticle[m][n]->GetPosition() + Vector3<T>(0.01,0.01,0.01) );
        //      ax = 0.1 * ( 
        //                  freq1 * (neighborsAvg.GetX() - x) / m_pcParticle[m][n]->GetMass()
        //                  + freq2 * (neighborsAvg.GetX()*neighborsAvg.GetX() - x*x) 
        //                //+ freq2 * (n * 1.0/*size/meshSize = 1.0*/ - x) // Not working for 3D
        //                - damp * vx
        //      );
        //      ay = 0.1 * ( 
        //                  freq1 * (neighborsAvg.GetY() - y) / m_pcParticle[m][n]->GetMass()
        //                  + freq2 * (neighborsAvg.GetY()*neighborsAvg.GetY() - y*y) 
        //                //+ freq2 * (m * 1.0/*size/meshSize = 1.0*/ - y) // Not working for 3D
        //                - damp * vy
        //      );
        //      az = 0.1 * ( 
        //                  freq1 * (neighborsAvg.GetZ() - z) / m_pcParticle[m][n]->GetMass()
        //                  + freq2 * (neighborsAvg.GetZ()*neighborsAvg.GetZ() - z*z) 
        //                //+ freq2 * (m * 1.0/*size/meshSize = 1.0*/ - z) // Not working for 3D
        //                - damp * vz
        //      );
        //      vx += ax;
        //      vy += ay;
        //      vz += az;
        //      m_pcParticle[m][n]->SetVelocity( vx, vy, vz );
        //      x += vx;
        //      y += vy;
        //      z += vz;
        //      m_pcParticle[m][n]->SetPosition( x, y, z );
        //  }
        //}
        //{ m = 0;
        //  n = 0;
        //  if ( !m_pcParticle[m][n]->GetFixStatus() ) {
        //      //-------------------------------
        //      // Get neighbors' positions
        //      nEast   = m_pcParticle[m][n+1]->GetPosition();
        //      nNorth  = m_pcParticle[m+1][n]->GetPosition();
        //      //-------------------------------
        //      nNorthEast  = m_pcParticle[m+1][n+1]->GetPosition();
        //      //-------------------------------
        //      //neighborsAvg = (nWest + nEast + nSouth + nNorth ) / 4.0;
        //      neighborsAvg = (nEast + nNorth + nNorthEast) / 3.0;
        //      m_pcParticle[m][n]->GetPosition().GetXYZ( x, y, z );
        //      m_pcParticle[m][n]->GetVelocity().GetXYZ( vx, vy, vz );

        //      //m_pcParticle[m][n]->SetPosition( m_pcParticle[m][n]->GetPosition() + Vector3<T>(0.01,0.01,0.01) );
        //      ax = 0.1 * ( 
        //                  freq1 * (neighborsAvg.GetX() - x) / m_pcParticle[m][n]->GetMass()
        //                  + freq2 * (neighborsAvg.GetX()*neighborsAvg.GetX() - x*x) 
        //                //+ freq2 * (n * 1.0/*size/meshSize = 1.0*/ - x) // Not working for 3D
        //                - damp * vx
        //      );
        //      ay = 0.1 * ( 
        //                  freq1 * (neighborsAvg.GetY() - y) / m_pcParticle[m][n]->GetMass()
        //                  + freq2 * (neighborsAvg.GetY()*neighborsAvg.GetY() - y*y) 
        //                //+ freq2 * (m * 1.0/*size/meshSize = 1.0*/ - y) // Not working for 3D
        //                - damp * vy
        //      );
        //      az = 0.1 * ( 
        //                  freq1 * (neighborsAvg.GetZ() - z) / m_pcParticle[m][n]->GetMass()
        //                  + freq2 * (neighborsAvg.GetZ()*neighborsAvg.GetZ() - z*z) 
        //                //+ freq2 * (m * 1.0/*size/meshSize = 1.0*/ - z) // Not working for 3D
        //                - damp * vz
        //      );
        //      vx += ax;
        //      vy += ay;
        //      vz += az;
        //      m_pcParticle[m][n]->SetVelocity( vx, vy, vz );
        //      x += vx;
        //      y += vy;
        //      z += vz;
        //      m_pcParticle[m][n]->SetPosition( x, y, z );
        //  }
        //}
        //{ m = 0;
        //  n = m_iNoCols-1;
        //  if ( !m_pcParticle[m][n]->GetFixStatus() ) {
        //      //-------------------------------
        //      // Get neighbors' positions
        //      nWest   = m_pcParticle[m][n-1]->GetPosition();
        //      nNorth  = m_pcParticle[m+1][n]->GetPosition();
        //      //-------------------------------
        //      nNorthWest  = m_pcParticle[m+1][n-1]->GetPosition();
        //      //-------------------------------
        //      //neighborsAvg = (nWest + nEast + nSouth + nNorth ) / 4.0;
        //      neighborsAvg = (nWest + nNorth + nNorthWest) / 3.0;
        //      m_pcParticle[m][n]->GetPosition().GetXYZ( x, y, z );
        //      m_pcParticle[m][n]->GetVelocity().GetXYZ( vx, vy, vz );

        //      //m_pcParticle[m][n]->SetPosition( m_pcParticle[m][n]->GetPosition() + Vector3<T>(0.01,0.01,0.01) );
        //      ax = 0.1 * ( 
        //                  freq1 * (neighborsAvg.GetX() - x) / m_pcParticle[m][n]->GetMass()
        //                  + freq2 * (neighborsAvg.GetX()*neighborsAvg.GetX() - x*x) 
        //                //+ freq2 * (n * 1.0/*size/meshSize = 1.0*/ - x) // Not working for 3D
        //                - damp * vx
        //      );
        //      ay = 0.1 * ( 
        //                  freq1 * (neighborsAvg.GetY() - y) / m_pcParticle[m][n]->GetMass()
        //                  + freq2 * (neighborsAvg.GetY()*neighborsAvg.GetY() - y*y) 
        //                //+ freq2 * (m * 1.0/*size/meshSize = 1.0*/ - y) // Not working for 3D
        //                - damp * vy
        //      );
        //      az = 0.1 * ( 
        //                  freq1 * (neighborsAvg.GetZ() - z) / m_pcParticle[m][n]->GetMass()
        //                  + freq2 * (neighborsAvg.GetZ()*neighborsAvg.GetZ() - z*z) 
        //                //+ freq2 * (m * 1.0/*size/meshSize = 1.0*/ - z) // Not working for 3D
        //                - damp * vz
        //      );
        //      vx += ax;
        //      vy += ay;
        //      vz += az;
        //      m_pcParticle[m][n]->SetVelocity( vx, vy, vz );
        //      x += vx;
        //      y += vy;
        //      z += vz;
        //      m_pcParticle[m][n]->SetPosition( x, y, z );
        //  }
        //}
        //{ m = m_iNoRows-1;
        //  n = 0;
        //  if ( !m_pcParticle[m][n]->GetFixStatus() ) {
        //      //-------------------------------
        //      // Get neighbors' positions
        //      nEast   = m_pcParticle[m][n+1]->GetPosition();
        //      nSouth  = m_pcParticle[m-1][n]->GetPosition();
        //      //-------------------------------
        //      nSouthEast  = m_pcParticle[m-1][n+1]->GetPosition();
        //      //-------------------------------
        //      //neighborsAvg = (nWest + nEast + nSouth + nNorth ) / 4.0;
        //      neighborsAvg = (nEast + nSouth + nSouthEast) / 3.0;
        //      m_pcParticle[m][n]->GetPosition().GetXYZ( x, y, z );
        //      m_pcParticle[m][n]->GetVelocity().GetXYZ( vx, vy, vz );

        //      //m_pcParticle[m][n]->SetPosition( m_pcParticle[m][n]->GetPosition() + Vector3<T>(0.01,0.01,0.01) );
        //      ax = 0.1 * ( 
        //                  freq1 * (neighborsAvg.GetX() - x) / m_pcParticle[m][n]->GetMass()
        //                  + freq2 * (neighborsAvg.GetX()*neighborsAvg.GetX() - x*x) 
        //                //+ freq2 * (n * 1.0/*size/meshSize = 1.0*/ - x) // Not working for 3D
        //                - damp * vx
        //      );
        //      ay = 0.1 * ( 
        //                  freq1 * (neighborsAvg.GetY() - y) / m_pcParticle[m][n]->GetMass()
        //                  + freq2 * (neighborsAvg.GetY()*neighborsAvg.GetY() - y*y) 
        //                //+ freq2 * (m * 1.0/*size/meshSize = 1.0*/ - y) // Not working for 3D
        //                - damp * vy
        //      );
        //      az = 0.1 * ( 
        //                  freq1 * (neighborsAvg.GetZ() - z) / m_pcParticle[m][n]->GetMass()
        //                  + freq2 * (neighborsAvg.GetZ()*neighborsAvg.GetZ() - z*z) 
        //                //+ freq2 * (m * 1.0/*size/meshSize = 1.0*/ - z) // Not working for 3D
        //                - damp * vz
        //      );
        //      vx += ax;
        //      vy += ay;
        //      vz += az;
        //      m_pcParticle[m][n]->SetVelocity( vx, vy, vz );
        //      x += vx;
        //      y += vy;
        //      z += vz;
        //      m_pcParticle[m][n]->SetPosition( x, y, z );
        //  }
        //}
        //{ m = m_iNoRows-1;
        //  n = m_iNoCols-1;
        //  if ( !m_pcParticle[m][n]->GetFixStatus() ) {
        //      //-------------------------------
        //      // Get neighbors' positions
        //      nWest   = m_pcParticle[m][n-1]->GetPosition();
        //      nSouth  = m_pcParticle[m-1][n]->GetPosition();
        //      //-------------------------------
        //      nSouthWest  = m_pcParticle[m-1][n-1]->GetPosition();
        //      //-------------------------------
        //      //neighborsAvg = (nWest + nEast + nSouth + nNorth ) / 4.0;
        //      neighborsAvg = (nWest + nSouth + nSouthWest) / 3.0;
        //      m_pcParticle[m][n]->GetPosition().GetXYZ( x, y, z );
        //      m_pcParticle[m][n]->GetVelocity().GetXYZ( vx, vy, vz );

        //      //m_pcParticle[m][n]->SetPosition( m_pcParticle[m][n]->GetPosition() + Vector3<T>(0.01,0.01,0.01) );
        //      ax = 0.1 * ( 
        //                  freq1 * (neighborsAvg.GetX() - x) / m_pcParticle[m][n]->GetMass()
        //                  + freq2 * (neighborsAvg.GetX()*neighborsAvg.GetX() - x*x) 
        //                //+ freq2 * (n * 1.0/*size/meshSize = 1.0*/ - x) // Not working for 3D
        //                - damp * vx
        //      );
        //      ay = 0.1 * ( 
        //                  freq1 * (neighborsAvg.GetY() - y) / m_pcParticle[m][n]->GetMass()
        //                  + freq2 * (neighborsAvg.GetY()*neighborsAvg.GetY() - y*y) 
        //                //+ freq2 * (m * 1.0/*size/meshSize = 1.0*/ - y) // Not working for 3D
        //                - damp * vy
        //      );
        //      az = 0.1 * ( 
        //                  freq1 * (neighborsAvg.GetZ() - z) / m_pcParticle[m][n]->GetMass()
        //                  + freq2 * (neighborsAvg.GetZ()*neighborsAvg.GetZ() - z*z) 
        //                //+ freq2 * (m * 1.0/*size/meshSize = 1.0*/ - z) // Not working for 3D
        //                - damp * vz
        //      );
        //      vx += ax;
        //      vy += ay;
        //      vz += az;
        //      m_pcParticle[m][n]->SetVelocity( vx, vy, vz );
        //      x += vx;
        //      y += vy;
        //      z += vz;
        //      m_pcParticle[m][n]->SetPosition( x, y, z );
        //  }
        //}
        //------------------------------------------------------------
        CalParticleNormals();
    }
};
//CLASSCLASSCLASSCLASSCLASSCLASSCLASSCLASSCLASSCLASSCLASSCLASSCLASSCLASSCLASSCLASSCLASSCLASSCLASSCLA


//HELPERHELPERHELPERHELPERHELPERHELPERHELPERHELPERHELPERHELPERHELPERHELPERHELPERHELPERHELPERHELPERHE
//
//      DeformMeshWithCut   H E L P E R    F U N C T I O N    M E M B E R    D E F I N I T I O N S
//
//HELPERHELPERHELPERHELPERHELPERHELPERHELPERHELPERHELPERHELPERHELPERHELPERHELPERHELPERHELPERHELPERHE
//-------------------------------------------------------------------------------------------------
// Create Particles
template <typename T>
void DeformMeshWithCut<T>::CreateParticles( const Vector3<T> &topLeftPos, const Vector3<T> &bottomRightPos, T mass )
{
    // Local Variable
    Vector3<T> currentPosVector = topLeftPos;
    Vector3<T> diagonalVector = bottomRightPos - topLeftPos;
    T ptMass = mass / ( m_iNoCols * m_iNoRows );
    T colIncStep = diagonalVector[0] / (m_iNoCols-1);
    T rowIncStep = diagonalVector[1] / (m_iNoRows-1);
    //T zIncStep   = diagonalVector[2] / ((m_iNoCols + m_iNoRows)/2.0);
    T zIncStep   = diagonalVector[2] / (m_iNoRows-1);

    // Create Particle Pointers
    m_pcParticle = new Particle<T>**[m_iNoRows];
    for ( m = 0; m < m_iNoRows; m++ ) {
        m_pcParticle[m] = new Particle<T>*[m_iNoCols];
    }

    // Create Particles
    Vector3<T> zeroVector(0,0,0);
    for ( m = 0; m < m_iNoRows; m++ ) { // row loops
        for ( n = 0; n < m_iNoCols; n++ ) { // column loops
            m_pcParticle[m][n] = new Particle<T>( 
                                    ptMass,             // mass 
                                    currentPosVector,   // position
                                    zeroVector,         // velocity
                                    zeroVector,         // acceleration
                                    zeroVector );       // force
            currentPosVector[0] = currentPosVector[0] + colIncStep; // increment x position
        }
        currentPosVector[0] = topLeftPos[0];                    // reset x position
        currentPosVector[1] = currentPosVector[1] + rowIncStep; // increment y position
        currentPosVector[2] = currentPosVector[2] + zIncStep;   // increment z position
    }

    
    int iHalfNoCols = m_iNoCols / 2;
    T deepAdjustment[] = { -0.1, -0.4, -1.2, -0.4, -0.1 };
    for ( m = 5; m < m_iNoRows-5; m++ ) {   // row loops
        int iCount = 0;
        for ( n = iHalfNoCols-2; n <= iHalfNoCols+2; ++n ) {    // column loops
            m_pcParticle[m][n]->SetPosition( m_pcParticle[m][n]->GetPosition() + Vector3<T>( 0, deepAdjustment[iCount++], 0 ) );
        }
    }

    /*
    // Dump Particle Positions on Screen
    std::cout << "Particle Positions:\n";
    for ( m = 0; m < m_iNoRows; m++ ) { // row loops
        for ( n = 0; n < m_iNoCols; n++ ) { // column loops
            std::cout << "[" << m << "," << n << "] " << m_pcParticle[m][n]->m_vtPosit << " ";
        }
        std::cout << "\n";
    }
    //*/
}
//-------------------------------------------------------------------------------------------------
// Delete Particles
template <typename T>
void DeformMeshWithCut<T>::DeleteParticles()
{
    // Delete Particles
    for ( m = 0; m < m_iNoRows; m++ ) {
        for ( n = 0; n < m_iNoCols; n++ ) {
            delete m_pcParticle[m][n];
        }
        delete [] m_pcParticle[m];
    }
    delete [] m_pcParticle;
}
//HELPERHELPERHELPERHELPERHELPERHELPERHELPERHELPERHELPERHELPERHELPERHELPERHELPERHELPERHELPERHELPERHE


//MEMBERMEMBERMEMBERMEMBERMEMBERMEMBERMEMBERMEMBERMEMBERMEMBERMEMBERMEMBERMEMBERMEMBERMEMBERMEMBERME
//
//      DeformMeshWithCut   F U N C T I O N    M E M B E R    D E F I N I T I O N S
//
//MEMBERMEMBERMEMBERMEMBERMEMBERMEMBERMEMBERMEMBERMEMBERMEMBERMEMBERMEMBERMEMBERMEMBERMEMBERMEMBERME
//-------------------------------------------------------------------------------------------------
template <typename T>
DeformMeshWithCut<T>::DeformMeshWithCut( 
              int rows,             // number of rows (on world-y-axis)
              int cols,             // number of cols (on world-x-axis)
              const Vector3<T> &topLeftPos,         // start position (should be top left cornor)
              const Vector3<T> &bottomRightPos,     // end position (should be bottom right cornor)
              T mass,                                   // the whole mass
              BitmapImage *img )            // Image Pointer
    : OpenGLSupport( img ),     // parent class constructor
      m_iNoRows( rows ), m_iNoCols( cols )
{
    CreateParticles( topLeftPos, bottomRightPos, mass );
    CalParticleNormals();
}
//-------------------------------------------------------------------------------------------------
template <typename T>
DeformMeshWithCut<T>::~DeformMeshWithCut()
{
    DeleteParticles();
}
//-------------------------------------------------------------------------------------------------
// Calculate Particle Normals
template <typename T>
void DeformMeshWithCut<T>::CalParticleNormals()
{
/*
    -----------------
    | \  3  |  2  / |
    |   \   |   /   |
    | 4   \ | /  1  |   Normal of the middle (m,n) particle is the sum of the eight normals.
    |<------+------>|   E.g. normal 1 is the unit vector of the cross product of (m,n)-->(m,n+1) 
    | 5   / | \  8  |   vector with (m,n)-->(m-1,n+1) vector (rotated counter clockwise).
    |   /   |   \   |
    | /  6  |  7  \ |
    -----------------
*/
    // These particles have all 8 neighbors ------------------------------------
    for ( m = 1; m < m_iNoRows-1; m++ ) {
        for ( n = 1; n < m_iNoCols-1; n++ ) {
            // (1) (m,n)-->(m,n+1)   and (m,n)-->(m-1,n+1)
            m_pcParticle[m][n]->m_vtNormal 
                = ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m  ][n+1]->m_vtPosit) 
                  ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m-1][n+1]->m_vtPosit) ).Normalized();
            // (2) (m,n)-->(m-1,n+1) and (m,n)-->(m-1,n)
            m_pcParticle[m][n]->m_vtNormal 
                += ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m-1][n+1]->m_vtPosit) 
                   ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m-1][n  ]->m_vtPosit) ).Normalized();
            // (3) (m,n)-->(m-1,n)   and (m,n)-->(m-1,n-1)
            m_pcParticle[m][n]->m_vtNormal 
                += ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m-1][n  ]->m_vtPosit) 
                   ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m-1][n-1]->m_vtPosit) ).Normalized();
            // (4) (m,n)-->(m-1,n-1) and (m,n)-->(m,n-1)
            m_pcParticle[m][n]->m_vtNormal 
                += ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m-1][n-1]->m_vtPosit) 
                   ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m  ][n-1]->m_vtPosit) ).Normalized();
            // (5) (m,n)-->(m,n-1)   and (m,n)-->(m+1,n-1)
            m_pcParticle[m][n]->m_vtNormal 
                += ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m  ][n-1]->m_vtPosit) 
                   ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m+1][n-1]->m_vtPosit) ).Normalized();
            // (6) (m,n)-->(m+1,n-1) and (m,n)-->(m+1,n)
            m_pcParticle[m][n]->m_vtNormal 
                += ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m+1][n-1]->m_vtPosit) 
                   ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m+1][n  ]->m_vtPosit) ).Normalized();
            // (7) (m,n)-->(m+1,n)   and (m,n)-->(m+1,n+1)
            m_pcParticle[m][n]->m_vtNormal 
                += ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m+1][n  ]->m_vtPosit) 
                   ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m+1][n+1]->m_vtPosit) ).Normalized();
            // (8) (m,n)-->(m+1,n+1) and (m,n)-->(m,n+1)
            m_pcParticle[m][n]->m_vtNormal 
                += ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m+1][n+1]->m_vtPosit) 
                   ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m  ][n+1]->m_vtPosit) ).Normalized();
            // normalize the sum result vector
            m_pcParticle[m][n]->m_vtNormal.Normalized();
        }
    }

    //-------------------------------------------------------------------------
    // Normals of Particles on the Top Row (Except the Cornor Particle)
    m = 0;
    for ( n = 1; n < m_iNoCols-1; n++ ) {
        // (5) (m,n)-->(m,n-1)   and (m,n)-->(m+1,n-1)
        m_pcParticle[m][n]->m_vtNormal 
            = ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m  ][n-1]->m_vtPosit) 
              ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m+1][n-1]->m_vtPosit) ).Normalized();
        // (6) (m,n)-->(m+1,n-1) and (m,n)-->(m+1,n)
        m_pcParticle[m][n]->m_vtNormal 
            += ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m+1][n-1]->m_vtPosit) 
               ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m+1][n  ]->m_vtPosit) ).Normalized();
        // (7) (m,n)-->(m+1,n)   and (m,n)-->(m+1,n+1)
        m_pcParticle[m][n]->m_vtNormal 
            += ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m+1][n  ]->m_vtPosit) 
               ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m+1][n+1]->m_vtPosit) ).Normalized();
        // (8) (m,n)-->(m+1,n+1) and (m,n)-->(m,n+1)
        m_pcParticle[m][n]->m_vtNormal 
            += ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m+1][n+1]->m_vtPosit) 
               ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m  ][n+1]->m_vtPosit) ).Normalized();
        // normalize the sum result vector
        m_pcParticle[m][n]->m_vtNormal.Normalized();
    }
    //-------------------------------------------------------------------------
    // Normals of Particles on the Bottom Row (Except the Cornor Particle)
    m = m_iNoRows-1;
    for ( n = 1; n < m_iNoCols-1; n++ ) {
        // (1) (m,n)-->(m,n+1)   and (m,n)-->(m-1,n+1)
        m_pcParticle[m][n]->m_vtNormal 
            = ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m  ][n+1]->m_vtPosit) 
              ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m-1][n+1]->m_vtPosit) ).Normalized();
        // (2) (m,n)-->(m-1,n+1) and (m,n)-->(m-1,n)
        m_pcParticle[m][n]->m_vtNormal 
            += ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m-1][n+1]->m_vtPosit) 
               ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m-1][n  ]->m_vtPosit) ).Normalized();
        // (3) (m,n)-->(m-1,n)   and (m,n)-->(m-1,n-1)
        m_pcParticle[m][n]->m_vtNormal 
            += ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m-1][n  ]->m_vtPosit) 
               ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m-1][n-1]->m_vtPosit) ).Normalized();
        // (4) (m,n)-->(m-1,n-1) and (m,n)-->(m,n-1)
        m_pcParticle[m][n]->m_vtNormal 
            += ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m-1][n-1]->m_vtPosit) 
               ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m  ][n-1]->m_vtPosit) ).Normalized();
        // normalize the sum result vector
        m_pcParticle[m][n]->m_vtNormal.Normalized();
    }
    //-------------------------------------------------------------------------
    // Normals of Particles on the Left Column (Except the Cornor Particle)
    n = 0;
    for ( m = 1; m < m_iNoRows-1; m++ ) {
        // (7) (m,n)-->(m+1,n)   and (m,n)-->(m+1,n+1)
        m_pcParticle[m][n]->m_vtNormal 
            = ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m+1][n  ]->m_vtPosit) 
              ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m+1][n+1]->m_vtPosit) ).Normalized();
        // (8) (m,n)-->(m+1,n+1) and (m,n)-->(m,n+1)
        m_pcParticle[m][n]->m_vtNormal 
            += ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m+1][n+1]->m_vtPosit) 
               ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m  ][n+1]->m_vtPosit) ).Normalized();
        // (1) (m,n)-->(m,n+1)   and (m,n)-->(m-1,n+1)
        m_pcParticle[m][n]->m_vtNormal 
            += ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m  ][n+1]->m_vtPosit) 
               ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m-1][n+1]->m_vtPosit) ).Normalized();
        // (2) (m,n)-->(m-1,n+1) and (m,n)-->(m-1,n)
        m_pcParticle[m][n]->m_vtNormal 
            += ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m-1][n+1]->m_vtPosit) 
               ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m-1][n  ]->m_vtPosit) ).Normalized();
        // normalize the sum result vector
        m_pcParticle[m][n]->m_vtNormal.Normalized();
    }
    //-------------------------------------------------------------------------
    // Normals of Particles on the Right Row (Except the Cornor Particle)
    n = m_iNoCols-1;
    for ( m = 1; m < m_iNoRows-1; m++ ) {
        // (3) (m,n)-->(m-1,n)   and (m,n)-->(m-1,n-1)
        m_pcParticle[m][n]->m_vtNormal 
            = ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m-1][n  ]->m_vtPosit) 
              ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m-1][n-1]->m_vtPosit) ).Normalized();
        // (4) (m,n)-->(m-1,n-1) and (m,n)-->(m,n-1)
        m_pcParticle[m][n]->m_vtNormal 
            += ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m-1][n-1]->m_vtPosit) 
               ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m  ][n-1]->m_vtPosit) ).Normalized();
        // (5) (m,n)-->(m,n-1)   and (m,n)-->(m+1,n-1)
        m_pcParticle[m][n]->m_vtNormal 
            += ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m  ][n-1]->m_vtPosit) 
               ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m+1][n-1]->m_vtPosit) ).Normalized();
        // (6) (m,n)-->(m+1,n-1) and (m,n)-->(m+1,n)
        m_pcParticle[m][n]->m_vtNormal 
            += ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m+1][n-1]->m_vtPosit) 
               ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m+1][n  ]->m_vtPosit) ).Normalized();
        // normalize the sum result vector
        m_pcParticle[m][n]->m_vtNormal.Normalized();
    }
    
    // Top Left Cornor ---------------------------------------------------------
    m = 0; n = 0;
        // (7) (m,n)-->(m+1,n)   and (m,n)-->(m+1,n+1)
        m_pcParticle[m][n]->m_vtNormal 
            = ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m+1][n  ]->m_vtPosit) 
              ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m+1][n+1]->m_vtPosit) ).Normalized();
        // (8) (m,n)-->(m+1,n+1) and (m,n)-->(m,n+1)
        m_pcParticle[m][n]->m_vtNormal 
            += ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m+1][n+1]->m_vtPosit) 
               ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m  ][n+1]->m_vtPosit) ).Normalized();
        m_pcParticle[m][n]->m_vtNormal.Normalized();    // normalize the sum result vector
    // Bottom Left Cornor ------------------------------------------------------
    m = m_iNoRows-1; n = 0;
        // (1) (m,n)-->(m,n+1)   and (m,n)-->(m-1,n+1)
        m_pcParticle[m][n]->m_vtNormal 
            = ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m  ][n+1]->m_vtPosit) 
              ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m-1][n+1]->m_vtPosit) ).Normalized();
        // (2) (m,n)-->(m-1,n+1) and (m,n)-->(m-1,n)
        m_pcParticle[m][n]->m_vtNormal 
            += ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m-1][n+1]->m_vtPosit) 
               ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m-1][n  ]->m_vtPosit) ).Normalized();
        m_pcParticle[m][n]->m_vtNormal.Normalized();    // normalize the sum result vector
    // Bottom Right Cornor -----------------------------------------------------
    m = m_iNoRows-1; n = m_iNoCols-1;
        // (3) (m,n)-->(m-1,n)   and (m,n)-->(m-1,n-1)
        m_pcParticle[m][n]->m_vtNormal 
            = ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m-1][n  ]->m_vtPosit) 
              ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m-1][n-1]->m_vtPosit) ).Normalized();
        // (4) (m,n)-->(m-1,n-1) and (m,n)-->(m,n-1)
        m_pcParticle[m][n]->m_vtNormal 
            += ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m-1][n-1]->m_vtPosit) 
               ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m  ][n-1]->m_vtPosit) ).Normalized();
        m_pcParticle[m][n]->m_vtNormal.Normalized();    // normalize the sum result vector
    // Top Right Cornor --------------------------------------------------------
    m = 0; n = m_iNoCols-1;
        // (5) (m,n)-->(m,n-1)   and (m,n)-->(m+1,n-1)
        m_pcParticle[m][n]->m_vtNormal 
            = ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m  ][n-1]->m_vtPosit) 
              ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m+1][n-1]->m_vtPosit) ).Normalized();
        // (6) (m,n)-->(m+1,n-1) and (m,n)-->(m+1,n)
        m_pcParticle[m][n]->m_vtNormal 
            += ( (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m+1][n-1]->m_vtPosit) 
               ^ (m_pcParticle[m][n]->m_vtPosit - m_pcParticle[m+1][n  ]->m_vtPosit) ).Normalized();
        m_pcParticle[m][n]->m_vtNormal.Normalized();    // normalize the sum result vector
}
//-------------------------------------------------------------------------------------------------
// Calculate Particle Normals
template <typename T>
void DeformMeshWithCut<T>::CalParticleNormals( int i )
{
/*
    -----------------
    | \  3  |  2  / |
    |   \   |   /   |
    | 4   \ | /  1  |   Normal of the middle (m,n) particle is the sum of the eight normals.
    |<------+------>|   E.g. normal 1 is the unit vector of the cross product of (m,n)-->(m,n+1) 
    | 5   / | \  8  |   vector with (m,n)-->(m-1,n+1) vector (rotated counter clockwise).
    |   /   |   \   |
    | /  6  |  7  \ |
    -----------------
*/
    // These particles have all 8 neighbors ------------------------------------
    for ( m = 1; m < m_iNoRows-1; m++ ) {
        for ( n = 1; n < m_iNoCols-1; n++ ) {
            // (1) (m,n)-->(m,n+1)   and (m,n)-->(m-1,n+1)
            m_pcParticle[m  ][n  ]->m_vtNormal[0] = 0;
            m_pcParticle[m  ][n  ]->m_vtNormal[1] = 0;
            m_pcParticle[m  ][n  ]->m_vtNormal[2] = 0;
            AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                       m_pcParticle[m  ][n  ]->m_vtPosit,
                       m_pcParticle[m  ][n+1]->m_vtPosit,
                       m_pcParticle[m-1][n+1]->m_vtPosit );
            // (2) (m,n)-->(m-1,n+1) and (m,n)-->(m-1,n)
            AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                       m_pcParticle[m  ][n  ]->m_vtPosit,
                       m_pcParticle[m-1][n+1]->m_vtPosit,
                       m_pcParticle[m-1][n  ]->m_vtPosit );
            // (3) (m,n)-->(m-1,n)   and (m,n)-->(m-1,n-1)
            AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                       m_pcParticle[m  ][n  ]->m_vtPosit,
                       m_pcParticle[m-1][n  ]->m_vtPosit,
                       m_pcParticle[m-1][n-1]->m_vtPosit );
            // (4) (m,n)-->(m-1,n-1) and (m,n)-->(m,n-1)
            AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                       m_pcParticle[m  ][n  ]->m_vtPosit,
                       m_pcParticle[m-1][n-1]->m_vtPosit,
                       m_pcParticle[m  ][n-1]->m_vtPosit );
            // (5) (m,n)-->(m,n-1)   and (m,n)-->(m+1,n-1)
            AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                       m_pcParticle[m  ][n  ]->m_vtPosit,
                       m_pcParticle[m  ][n-1]->m_vtPosit,
                       m_pcParticle[m+1][n-1]->m_vtPosit );
            // (6) (m,n)-->(m+1,n-1) and (m,n)-->(m+1,n)
            AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                       m_pcParticle[m  ][n  ]->m_vtPosit,
                       m_pcParticle[m+1][n-1]->m_vtPosit,
                       m_pcParticle[m+1][n  ]->m_vtPosit );
            // (7) (m,n)-->(m+1,n)   and (m,n)-->(m+1,n+1)
            AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                       m_pcParticle[m  ][n  ]->m_vtPosit,
                       m_pcParticle[m+1][n  ]->m_vtPosit,
                       m_pcParticle[m+1][n+1]->m_vtPosit );
            // (8) (m,n)-->(m+1,n+1) and (m,n)-->(m,n+1)
            AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                       m_pcParticle[m  ][n  ]->m_vtPosit,
                       m_pcParticle[m+1][n+1]->m_vtPosit,
                       m_pcParticle[m  ][n+1]->m_vtPosit );
            // normalize the sum result vector
            m_pcParticle[m][n]->m_vtNormal.Normalized();
        }
    }

    //-------------------------------------------------------------------------
    // Normals of Particles on the Top Row (Except the Cornor Particle)
    m = 0;
    for ( n = 1; n < m_iNoCols-1; n++ ) {
        // (5) (m,n)-->(m,n-1)   and (m,n)-->(m+1,n-1)
        m_pcParticle[m][n]->m_vtNormal[0] = 0;
        m_pcParticle[m][n]->m_vtNormal[1] = 0;
        m_pcParticle[m][n]->m_vtNormal[2] = 0;
        AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                   m_pcParticle[m  ][n  ]->m_vtPosit,
                   m_pcParticle[m  ][n-1]->m_vtPosit,
                   m_pcParticle[m+1][n-1]->m_vtPosit );
        // (6) (m,n)-->(m+1,n-1) and (m,n)-->(m+1,n)
        AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                   m_pcParticle[m  ][n  ]->m_vtPosit,
                   m_pcParticle[m+1][n-1]->m_vtPosit,
                   m_pcParticle[m+1][n  ]->m_vtPosit );
        // (7) (m,n)-->(m+1,n)   and (m,n)-->(m+1,n+1)
        AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                   m_pcParticle[m  ][n  ]->m_vtPosit,
                   m_pcParticle[m+1][n  ]->m_vtPosit,
                   m_pcParticle[m+1][n+1]->m_vtPosit );
        // (8) (m,n)-->(m+1,n+1) and (m,n)-->(m,n+1)
        AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                   m_pcParticle[m  ][n  ]->m_vtPosit,
                   m_pcParticle[m+1][n+1]->m_vtPosit,
                   m_pcParticle[m  ][n+1]->m_vtPosit );
        // normalize the sum result vector
        m_pcParticle[m][n]->m_vtNormal.Normalized();
    }
    //-------------------------------------------------------------------------
    // Normals of Particles on the Bottom Row (Except the Cornor Particle)
    m = m_iNoRows-1;
    for ( n = 1; n < m_iNoCols-1; n++ ) {
        // (1) (m,n)-->(m,n+1)   and (m,n)-->(m-1,n+1)
        m_pcParticle[m][n]->m_vtNormal[0] = 0;
        m_pcParticle[m][n]->m_vtNormal[1] = 0;
        m_pcParticle[m][n]->m_vtNormal[2] = 0;
        AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                   m_pcParticle[m  ][n  ]->m_vtPosit,
                   m_pcParticle[m  ][n+1]->m_vtPosit,
                   m_pcParticle[m-1][n+1]->m_vtPosit );
        // (2) (m,n)-->(m-1,n+1) and (m,n)-->(m-1,n)
        AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                   m_pcParticle[m  ][n  ]->m_vtPosit,
                   m_pcParticle[m-1][n+1]->m_vtPosit,
                   m_pcParticle[m-1][n  ]->m_vtPosit );
        // (3) (m,n)-->(m-1,n)   and (m,n)-->(m-1,n-1)
        AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                   m_pcParticle[m  ][n  ]->m_vtPosit,
                   m_pcParticle[m-1][n  ]->m_vtPosit,
                   m_pcParticle[m-1][n-1]->m_vtPosit );
        // (4) (m,n)-->(m-1,n-1) and (m,n)-->(m,n-1)
        AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                   m_pcParticle[m  ][n  ]->m_vtPosit,
                   m_pcParticle[m-1][n-1]->m_vtPosit,
                   m_pcParticle[m  ][n-1]->m_vtPosit );
        // normalize the sum result vector
        m_pcParticle[m][n]->m_vtNormal.Normalized();
    }
    //-------------------------------------------------------------------------
    // Normals of Particles on the Left Column (Except the Cornor Particle)
    n = 0;
    for ( m = 1; m < m_iNoRows-1; m++ ) {
        // (7) (m,n)-->(m+1,n)   and (m,n)-->(m+1,n+1)
        m_pcParticle[m][n]->m_vtNormal[0] = 0;
        m_pcParticle[m][n]->m_vtNormal[1] = 0;
        m_pcParticle[m][n]->m_vtNormal[2] = 0;
        AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                   m_pcParticle[m  ][n  ]->m_vtPosit,
                   m_pcParticle[m+1][n  ]->m_vtPosit,
                   m_pcParticle[m+1][n+1]->m_vtPosit );
        // (8) (m,n)-->(m+1,n+1) and (m,n)-->(m,n+1)
        AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                   m_pcParticle[m  ][n  ]->m_vtPosit,
                   m_pcParticle[m+1][n+1]->m_vtPosit,
                   m_pcParticle[m  ][n+1]->m_vtPosit );
        // (1) (m,n)-->(m,n+1)   and (m,n)-->(m-1,n+1)
        AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                   m_pcParticle[m  ][n  ]->m_vtPosit,
                   m_pcParticle[m  ][n+1]->m_vtPosit,
                   m_pcParticle[m-1][n+1]->m_vtPosit );
        // (2) (m,n)-->(m-1,n+1) and (m,n)-->(m-1,n)
        AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                   m_pcParticle[m  ][n  ]->m_vtPosit,
                   m_pcParticle[m-1][n+1]->m_vtPosit,
                   m_pcParticle[m-1][n  ]->m_vtPosit );
        // normalize the sum result vector
        m_pcParticle[m][n]->m_vtNormal.Normalized();
    }
    //-------------------------------------------------------------------------
    // Normals of Particles on the Right Row (Except the Cornor Particle)
    n = m_iNoCols-1;
    for ( m = 1; m < m_iNoRows-1; m++ ) {
        // (3) (m,n)-->(m-1,n)   and (m,n)-->(m-1,n-1)
        m_pcParticle[m][n]->m_vtNormal[0] = 0;
        m_pcParticle[m][n]->m_vtNormal[1] = 0;
        m_pcParticle[m][n]->m_vtNormal[2] = 0;
        AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                   m_pcParticle[m  ][n  ]->m_vtPosit,
                   m_pcParticle[m-1][n  ]->m_vtPosit,
                   m_pcParticle[m-1][n-1]->m_vtPosit );
        // (4) (m,n)-->(m-1,n-1) and (m,n)-->(m,n-1)
        AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                   m_pcParticle[m  ][n  ]->m_vtPosit,
                   m_pcParticle[m-1][n-1]->m_vtPosit,
                   m_pcParticle[m  ][n-1]->m_vtPosit );
        // (5) (m,n)-->(m,n-1)   and (m,n)-->(m+1,n-1)
        AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                   m_pcParticle[m  ][n  ]->m_vtPosit,
                   m_pcParticle[m  ][n-1]->m_vtPosit,
                   m_pcParticle[m+1][n-1]->m_vtPosit );
        // (6) (m,n)-->(m+1,n-1) and (m,n)-->(m+1,n)
        AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                   m_pcParticle[m  ][n  ]->m_vtPosit,
                   m_pcParticle[m+1][n-1]->m_vtPosit,
                   m_pcParticle[m+1][n  ]->m_vtPosit );
        // normalize the sum result vector
        m_pcParticle[m][n]->m_vtNormal.Normalized();
    }
    
    // Top Left Cornor ---------------------------------------------------------
    m = 0; n = 0;
        // (7) (m,n)-->(m+1,n)   and (m,n)-->(m+1,n+1)
        m_pcParticle[m][n]->m_vtNormal[0] = 0;
        m_pcParticle[m][n]->m_vtNormal[1] = 0;
        m_pcParticle[m][n]->m_vtNormal[2] = 0;
        AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                   m_pcParticle[m  ][n  ]->m_vtPosit,
                   m_pcParticle[m+1][n  ]->m_vtPosit,
                   m_pcParticle[m+1][n+1]->m_vtPosit );
        // (8) (m,n)-->(m+1,n+1) and (m,n)-->(m,n+1)
        AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                   m_pcParticle[m  ][n  ]->m_vtPosit,
                   m_pcParticle[m+1][n+1]->m_vtPosit,
                   m_pcParticle[m  ][n+1]->m_vtPosit );
        m_pcParticle[m][n]->m_vtNormal.Normalized();    // normalize the sum result vector
    // Bottom Left Cornor ------------------------------------------------------
    m = m_iNoRows-1; n = 0;
        // (1) (m,n)-->(m,n+1)   and (m,n)-->(m-1,n+1)
        m_pcParticle[m][n]->m_vtNormal[0] = 0;
        m_pcParticle[m][n]->m_vtNormal[1] = 0;
        m_pcParticle[m][n]->m_vtNormal[2] = 0;
        AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                   m_pcParticle[m  ][n  ]->m_vtPosit,
                   m_pcParticle[m  ][n+1]->m_vtPosit,
                   m_pcParticle[m-1][n+1]->m_vtPosit );
        // (2) (m,n)-->(m-1,n+1) and (m,n)-->(m-1,n)
        AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                   m_pcParticle[m  ][n  ]->m_vtPosit,
                   m_pcParticle[m-1][n+1]->m_vtPosit,
                   m_pcParticle[m-1][n  ]->m_vtPosit );
        m_pcParticle[m][n]->m_vtNormal.Normalized();    // normalize the sum result vector
    // Bottom Right Cornor -----------------------------------------------------
    m = m_iNoRows-1; n = m_iNoCols-1;
        // (3) (m,n)-->(m-1,n)   and (m,n)-->(m-1,n-1)
        m_pcParticle[m][n]->m_vtNormal[0] = 0;
        m_pcParticle[m][n]->m_vtNormal[1] = 0;
        m_pcParticle[m][n]->m_vtNormal[2] = 0;
        AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                   m_pcParticle[m  ][n  ]->m_vtPosit,
                   m_pcParticle[m-1][n  ]->m_vtPosit,
                   m_pcParticle[m-1][n-1]->m_vtPosit );
        // (4) (m,n)-->(m-1,n-1) and (m,n)-->(m,n-1)
        AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                   m_pcParticle[m  ][n  ]->m_vtPosit,
                   m_pcParticle[m-1][n-1]->m_vtPosit,
                   m_pcParticle[m  ][n-1]->m_vtPosit );
        m_pcParticle[m][n]->m_vtNormal.Normalized();    // normalize the sum result vector
    // Top Right Cornor --------------------------------------------------------
    m = 0; n = m_iNoCols-1;
        // (5) (m,n)-->(m,n-1)   and (m,n)-->(m+1,n-1)
        m_pcParticle[m][n]->m_vtNormal[0] = 0;
        m_pcParticle[m][n]->m_vtNormal[1] = 0;
        m_pcParticle[m][n]->m_vtNormal[2] = 0;
        AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                   m_pcParticle[m  ][n  ]->m_vtPosit,
                   m_pcParticle[m  ][n-1]->m_vtPosit,
                   m_pcParticle[m+1][n-1]->m_vtPosit );
        // (6) (m,n)-->(m+1,n-1) and (m,n)-->(m+1,n)
        AddNormal( m_pcParticle[m  ][n  ]->m_vtNormal,
                   m_pcParticle[m  ][n  ]->m_vtPosit,
                   m_pcParticle[m+1][n-1]->m_vtPosit,
                   m_pcParticle[m+1][n  ]->m_vtPosit );
        m_pcParticle[m][n]->m_vtNormal.Normalized();    // normalize the sum result vector
}
/*
//-------------------------------------------------------------------------------------------------
// Clear Force of Particles to Zeros
template <typename T>
void DeformMeshWithCut<T>::ClearForce()
{
    T x = 0, y = 0, z = 0;
    for ( m = 0; m < m_iNoRows; m++ ) {
        for ( n = 0; n < m_iNoCols; n++ ) {
            m_pcParticle[m][n]->m_vtForce( x, y, z );
        }
    }
}
//-------------------------------------------------------------------------------------------------
// Calculate (Weight) Force due to (Acceleration of) Gravity
template <typename T>
void DeformMeshWithCut<T>::ForceDueToGravity( T g )
{
    
}
//-------------------------------------------------------------------------------------------------
// Calculate (a Loss) Force due to Viscous Damping
template <typename T>
void DeformMeshWithCut<T>::ForceDueToViscousDamping()
{
}
//-------------------------------------------------------------------------------------------------
// Calculate (Wind) Force due to Viscous Fluid Moving
template <typename T>
void DeformMeshWithCut<T>::ForceDueToViscousFluidMoving()
{
}
*/
//-------------------------------------------------------------------------------------------------
// Display Via OpenGL
template <typename T>
void DeformMeshWithCut<T>::DisplayGL( OpenGL::Enum::DrawMode DM )
{
    //int ptSize[1];
    glColor3f( 1.0f, 1.0f, 1.0f );

    DeformUpdate();
    /*
    for ( count = 0; count < 8; count++ ) {
        StepSimulation( 0.0005 );
    }
    */

    glPushAttrib( GL_ALL_ATTRIB_BITS );

    //SetMaterial( OpenGL::Enum::METAL_GOLD );
    //SetMaterial( OpenGL::Enum::METAL_SILVER );
    //SetMaterial( OpenGL::Enum::METAL_BRONZE );
    //SetMaterial( OpenGL::Enum::RED_01 );
    SetMaterial( OpenGL::Enum::WHITE_01 );
    ApplyMaterial( OpenGL::Enum::FRONT );
    SetMaterial( OpenGL::Enum::METAL_GOLD );
    material.ApplyMaterial( OpenGL::Enum::BACK );
    
    switch (DM) {
        // Draw Points for DeformMeshWithCut
        case OpenGL::Enum::POINT:
            // Set Point Size
            //glGetIntegerv( GL_POINT_SIZE, ptSize );
            //glPointSize( 2 );
            glBegin( GL_POINTS );
            for ( m = 0; m < m_iNoRows; m++ ) {
                for ( n = 0; n < m_iNoCols; n++ ) {
                    glVertex3f( m_pcParticle[m][n]->m_vtPosit[0], 
                                m_pcParticle[m][n]->m_vtPosit[1],
                                m_pcParticle[m][n]->m_vtPosit[2] );
                }
            }
            glEnd();

            //-------------------------------------------
            // Display Starting Points
            glDisable( GL_LIGHTING );
            glColor3f( 0.0f, 1.0f, 0.0f );
            glPointSize( 7 );
            m = 0; n = 0;
            glBegin( GL_POINTS );
                glVertex3f( m_pcParticle[m][n]->m_vtPosit[0], 
                            m_pcParticle[m][n]->m_vtPosit[1],
                            m_pcParticle[m][n]->m_vtPosit[2] );
            glEnd();
            glColor3f( 0.0f, 1.0f, 1.0f );
            glPointSize( 5 );
            m = 0; n = 1;
            glBegin( GL_POINTS );
                glVertex3f( m_pcParticle[m][n]->m_vtPosit[0], 
                            m_pcParticle[m][n]->m_vtPosit[1],
                            m_pcParticle[m][n]->m_vtPosit[2] );
            glEnd();
            glPointSize( 1 );
            glEnable( GL_LIGHTING );
            //-------------------------------------------

            // Restore Point Size
            //glPointSize( ptSize[0] );
            break;

        // Draw Wire Frame for DeformMeshWithCut
        case OpenGL::Enum::WIRE_FRAME:
            // Draw Row Lines
            for ( m = 0; m < m_iNoRows; m++ ) {
                glBegin( GL_LINE_STRIP );
                for ( n = 0; n < m_iNoCols; n++ ) {
                    glVertex3f( m_pcParticle[m][n]->m_vtPosit[0], 
                                m_pcParticle[m][n]->m_vtPosit[1],
                                m_pcParticle[m][n]->m_vtPosit[2] );
                }
                glEnd();
            }
            // Draw Column Lines
            for ( n = 0; n < m_iNoCols; n++ ) {
                glBegin( GL_LINE_STRIP );
                for ( m = 0; m < m_iNoRows; m++ ) {
                    glVertex3f( m_pcParticle[m][n]->m_vtPosit[0], 
                                m_pcParticle[m][n]->m_vtPosit[1],
                                m_pcParticle[m][n]->m_vtPosit[2] );
                }
                glEnd();
            }

            /*
            //-------------------------------------------
            // Display Starting Points
            glDisable( GL_LIGHTING );
            glColor3f( 0.0f, 1.0f, 0.0f );
            glPointSize( 7 );
            m = 0; n = 0;
            glBegin( GL_POINTS );
                glVertex3f( m_pcParticle[m][n]->m_vtPosit[0], 
                            m_pcParticle[m][n]->m_vtPosit[1],
                            m_pcParticle[m][n]->m_vtPosit[2] );
            glEnd();
            glColor3f( 0.0f, 1.0f, 1.0f );
            glPointSize( 5 );
            m = 0; n = 1;
            glBegin( GL_POINTS );
                glVertex3f( m_pcParticle[m][n]->m_vtPosit[0], 
                            m_pcParticle[m][n]->m_vtPosit[1],
                            m_pcParticle[m][n]->m_vtPosit[2] );
            glEnd();
            glPointSize( 1 );
            glEnable( GL_LIGHTING );
            //-------------------------------------------
            //*/

            /*
            // Draw Diagonal Lines
            glBegin( GL_LINES );
            for ( m = 0; m < m_iNoRows-1; m++ ) {
                for ( n = 0; n < m_iNoCols-1; n++ ) {
                    // Draw a line from [m][n] --- [m+1][n+1]
                    glVertex3f( m_pcParticle[m][n]->m_vtPosit[0], 
                                m_pcParticle[m][n]->m_vtPosit[1],
                                m_pcParticle[m][n]->m_vtPosit[2] );
                    glVertex3f( m_pcParticle[m+1][n+1]->m_vtPosit[0], 
                                m_pcParticle[m+1][n+1]->m_vtPosit[1],
                                m_pcParticle[m+1][n+1]->m_vtPosit[2] );
                    // Draw a line from [m+1][n] to [m][n+1]
                    glVertex3f( m_pcParticle[m+1][n]->m_vtPosit[0], 
                                m_pcParticle[m+1][n]->m_vtPosit[1],
                                m_pcParticle[m+1][n]->m_vtPosit[2] );
                    glVertex3f( m_pcParticle[m][n+1]->m_vtPosit[0], 
                                m_pcParticle[m][n+1]->m_vtPosit[1],
                                m_pcParticle[m][n+1]->m_vtPosit[2] );
                }
            }
            glEnd();
            //*/
            break;

        // Draw Polygons for DeformMeshWithCut
        case OpenGL::Enum::POLYGON:
            // With Texture
            if ( m_pcImg != NULL ) {

                //glPushAttrib( GL_ENABLE_BIT );
                //glEnable( GL_TEXTURE_2D );
                DrawTexture();

                float tr, tc;
                float trStep = 1.0 / (m_iNoRows-1);
                float tcStep = 1.0 / (m_iNoCols-1);
                for ( m = 0, tr = 1.0; m < m_iNoRows-1; m++, tr -= trStep ) {
                    glBegin( GL_QUAD_STRIP );
                    for ( n = 0, tc = 0.0; n < m_iNoCols; n++, tc += tcStep ) {
                        glNormal3f( m_pcParticle[m][n]->m_vtNormal[0], 
                                    m_pcParticle[m][n]->m_vtNormal[1],
                                    m_pcParticle[m][n]->m_vtNormal[2] );
                        glTexCoord2f(tc, tr);
                        glVertex3f( m_pcParticle[m][n]->m_vtPosit[0], 
                                    m_pcParticle[m][n]->m_vtPosit[1],
                                    m_pcParticle[m][n]->m_vtPosit[2] );

                        glNormal3f( m_pcParticle[m+1][n]->m_vtNormal[0], 
                                    m_pcParticle[m+1][n]->m_vtNormal[1],
                                    m_pcParticle[m+1][n]->m_vtNormal[2] );
                        glTexCoord2f(tc, tr-trStep);
                        glVertex3f( m_pcParticle[m+1][n]->m_vtPosit[0], 
                                    m_pcParticle[m+1][n]->m_vtPosit[1],
                                    m_pcParticle[m+1][n]->m_vtPosit[2] );
                    }
                    glEnd();
                }
                //glBindTexture( GL_TEXTURE_2D, 0 );
                //glPopAttrib();
                //glDisable( GL_TEXTURE_2D );
            }
            // Without Texture
            else {
                for ( m = 0; m < m_iNoRows-1; m++ ) {
                    glBegin( GL_QUAD_STRIP );
                    for ( n = 0; n < m_iNoCols; n++ ) {
                        glNormal3f( m_pcParticle[m][n]->m_vtNormal[0], 
                                    m_pcParticle[m][n]->m_vtNormal[1],
                                    m_pcParticle[m][n]->m_vtNormal[2] );
                        glVertex3f( m_pcParticle[m][n]->m_vtPosit[0], 
                                    m_pcParticle[m][n]->m_vtPosit[1],
                                    m_pcParticle[m][n]->m_vtPosit[2] );

                        glNormal3f( m_pcParticle[m+1][n]->m_vtNormal[0], 
                                    m_pcParticle[m+1][n]->m_vtNormal[1],
                                    m_pcParticle[m+1][n]->m_vtNormal[2] );
                        glVertex3f( m_pcParticle[m+1][n]->m_vtPosit[0], 
                                    m_pcParticle[m+1][n]->m_vtPosit[1],
                                    m_pcParticle[m+1][n]->m_vtPosit[2] );
                    }
                    glEnd();
                }
            }
            break;
    }

    // DEBUG
    //DisplayGL_ParticleNormals();

    /*
    //-------------------------------------------
    // Display Starting Points
    glDisable( GL_LIGHTING );
    glColor3f( 0.0f, 1.0f, 0.0f );
    glPointSize( 7 );
    m = 0; n = 0;
    glBegin( GL_POINTS );
        glVertex3f( m_pcParticle[m][n]->m_vtPosit[0], 
                    m_pcParticle[m][n]->m_vtPosit[1],
                    m_pcParticle[m][n]->m_vtPosit[2] );
    glEnd();
    glColor3f( 0.0f, 1.0f, 1.0f );
    glPointSize( 5 );
    m = 0; n = 1;
    glBegin( GL_POINTS );
        glVertex3f( m_pcParticle[m][n]->m_vtPosit[0], 
                    m_pcParticle[m][n]->m_vtPosit[1],
                    m_pcParticle[m][n]->m_vtPosit[2] );
    glEnd();
    glPointSize( 1 );
    glEnable( GL_LIGHTING );
    //-------------------------------------------
    //*/

    glPopAttrib();
    // END TEST
}
//-------------------------------------------------------------------------------------------------
// Display Particle Normals Via OpenGL
template <typename T>
void DeformMeshWithCut<T>::DisplayGL_ParticleNormals()
{
    T x, y, z;

    // Draw Normal of each particle
    glPushAttrib( GL_ENABLE_BIT | GL_CURRENT_BIT );
    glEnable( GL_TEXTURE_2D );
    //glDisable( GL_LIGHTING );
    glColor3f( 1.0f, 0.0f, 0.0f );
    glBegin( GL_LINES );
    for ( m = 0; m < m_iNoRows; m++ ) {
        for ( n = 0; n < m_iNoCols; n++ ) {
            m_pcParticle[m][n]->m_vtPosit.GetXYZ( x, y, z );
            glVertex3f( x, y, z );
            glVertex3f( x + m_pcParticle[m][n]->m_vtNormal[0], 
                        y + m_pcParticle[m][n]->m_vtNormal[1],
                        z + m_pcParticle[m][n]->m_vtNormal[2] );
        }
    }
    glEnd();
    glPopAttrib();
}
//================================================================================================
// FRIEND FUNCTIONS
//-------------------------------------------------------------------------------------------------
// put it through ostream
/*template <typename T>
ostream &operator<<( ostream &output, const DeformMeshWithCut<T> &Obj )
{
    output  << "DeformMeshWithCut ==> "
            << "\n  number of rows:       " << Obj.m_iNoRows
            << "\n  number of columns:    " << Obj.m_iNoCols
            << "\n  viscous damping:      " << Obj.m_tKvd
            << "\n  wind damping:         " << Obj.m_tKvfm
            << "\n  stiffness of structural, shear, and flexion springs: " 
                    << Obj.m_tKs[0] << ", " << Obj.m_tKs[1] << ", and " << Obj.m_tKs[2]
            << "\n  damping constant for structural, shear, and flexion springs: "
                    << Obj.m_tKd[0] << ", " << Obj.m_tKd[1] << ", and " << Obj.m_tKd[2]
            << "\n  length constrain for structural, shear, and flextion springs: "
                    << Obj.m_tLengthConstrain[0] << ", " 
                    << Obj.m_tLengthConstrain[1] << ", and " 
                    << Obj.m_tLengthConstrain[2]
            << "\n";

    return output;
}*/

template <typename T>
void DeformMeshWithCut<T>::AddNormal( Vector3<T> & Dest, 
                       Vector3<T> const & A, Vector3<T> const & B, Vector3<T> const & C )
{
    T U[3] = { A[0]-B[0], A[1]-B[1], A[2]-B[2] };
    T V[3] = { A[0]-C[0], A[1]-C[1], A[2]-C[2] };
    T W[3] = { U[1]*V[2] - U[2]*V[1], U[2]*V[0] - U[0]*V[2], U[0]*V[1] - U[1]*V[0] };
    T SqrtLength = Math<T>::Sqrt( W[0]*W[0] + W[1]*W[1] + W[2]*W[2] );
    Dest[0] += W[0]/SqrtLength;
    Dest[1] += W[1]/SqrtLength;
    Dest[2] += W[2]/SqrtLength;
}
//MEMBERMEMBERMEMBERMEMBERMEMBERMEMBERMEMBERMEMBERMEMBERMEMBERMEMBERMEMBERMEMBER

//=============================================================================
END_NAMESPACE_TAPs
//-----------------------------------------------------------------------------
#endif
//345678901234567890123456789012345678901234567890123456789012345678901234567890
//--+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8