/*
    C++ Mesh Generation Library
    Copyright (c) 2000echet <eric at bechet dot ca>

    This file is part of the C++ Mesh Generation Library.

    See the NOTICE & LICENSE files for conditions.
*/
//---------------------------------------------------------------------------
#ifndef SIMPLEX_BASE_H
#define SIMPLEX_BASE_H
//---------------------------------------------------------------------------

#include "mesh_const.h"

#ifdef USE_HASH_TABLE
#include <hash_set>
#endif

#include <list>
#include <algorithm>
#include <vector>

#include "linear_algebra.h"
#include "vertex_base.h"
using namespace std;


/**class used to implement a simplex in any dimension, based on it's number of nodes,
space dimension, and vertex type*/
template<int nbpts,int dim,class V=Vertex_Base<dim> > class Simplex_Base
{
public:

///
  enum {Nb_Nodes=nbpts};
///
  enum {Nb_Faces=nbpts};
///
  typedef typename V::Iterator Vertex_Iterator;

protected :
///array of vertex iterators
  Vertex_Iterator vertex[nbpts];

public:
/// default constructor
  Simplex_Base ( void )
  {
  }

///Gets the normal if possible (depends on dim ans nbpts)
  double Normal ( Vector &norm );

///const Operator [] overloaded to return a vertex (non mutable)
  const Vertex_Iterator& operator [] ( int i ) const
  {
//      if ((i<nbpts)&&(i>=0)) return vertex[i]; else throw ERR_OUT_OF_BOUNDS;
    return vertex[i];
  }

///Operator [] overloaded to return a vertex (mutable)
  Vertex_Iterator& operator [] ( int i )
  {
//      if ((i<nbpts)&&(i>=0)) return vertex[i]; else throw ERR_OUT_OF_BOUNDS;
    return vertex[i];
  }

  /**gets the generalized volume of the simplex (det J in finite elements)
  it is implemented for space dim 1,2,3 and number of nodes 2,3,4
  otherwise returns 0*/
  double Volume ( void )
  {
    double vol=0;
    double tmp1,tmp2,tmp3,tmp4,tmp5;
    int i;

    switch ( nbpts )
    {
      case 1 :
        return 0.;
        break;

      case 2 :
      {
        if ( V::Dimension<1 )
        {
          return 0.;
        }

        //if # points = 2; compute the distance between the points
        for ( i=0; i<V::Dimension; ++i )
        {
          tmp1= ( *vertex[0] ) [i]- ( *vertex[1] ) [i];
          vol+=tmp1*tmp1;
        }

        if ( vol>=0. )
        {
          return sqrt ( vol );
        }
        else
        {
          return 0.;
        }

        break;
      }

      case 3 :
      {
        vol=0;
        //if the # of points = 3, get the triangle surface

        if ( V::Dimension==2 )
        {
          tmp1= ( *vertex[0] ) [0]- ( *vertex[1] ) [0];
          tmp2= ( *vertex[0] ) [1]- ( *vertex[2] ) [1];
          tmp3= ( *vertex[0] ) [1]- ( *vertex[1] ) [1];
          tmp4= ( *vertex[0] ) [0]- ( *vertex[2] ) [0];
          tmp5=tmp1*tmp2-tmp3*tmp4;
          vol=tmp5*tmp5;
          return sqrt ( vol ) /2.;
        }
        else if ( V::Dimension==3 )
          //This case is the same as above for a 3D representation
        {
          for ( i=0; i<3; ++i )
          {
            tmp1= ( *vertex[0] ) [i]- ( *vertex[1] ) [i];
            tmp2= ( *vertex[0] ) [ ( i+1 ) %3]- ( *vertex[2] ) [ ( i+1 ) %3];
            tmp3= ( *vertex[0] ) [ ( i+1 ) %3]- ( *vertex[1] ) [ ( i+1 ) %3];
            tmp4= ( *vertex[0] ) [i]- ( *vertex[2] ) [i];
            tmp5=tmp1*tmp2-tmp3*tmp4;
            vol+=tmp5*tmp5;
          }

          return sqrt ( vol ) /2.;
        }
        else
        {
          return 0.;
        }

        break;
      }

      case 4 :
      {
        if ( V::Dimension!=3 )
        {
          return 0.;
        }

        //this block estimates the volume of a tetrahedric element

        for ( i=0; i<3; ++i )
        {
          tmp1= ( *vertex[0] ) [i]- ( *vertex[1] ) [i];
          tmp2= ( *vertex[0] ) [ ( i+1 ) %3]- ( *vertex[2] ) [ ( i+1 ) %3];
          tmp3= ( *vertex[0] ) [ ( i+1 ) %3]- ( *vertex[3] ) [ ( i+1 ) %3];
          tmp4= ( *vertex[0] ) [ ( i+2 ) %3]- ( *vertex[2] ) [ ( i+2 ) %3];
          tmp5= ( *vertex[0] ) [ ( i+2 ) %3]- ( *vertex[3] ) [ ( i+2 ) %3];
          vol+=tmp1* ( tmp2*tmp5-tmp4*tmp3 );
        }

        return fabs ( vol ) /6;
        break;
      }

      default :
        return 0.;

    }
  }

  /**Simplex_Base overloaded constructor (passes a pointer to Vertex_Iterator
  and calls Set_vertices function describes later)*/
  Simplex_Base ( Vertex_Iterator v[] )
  {
    Set_Vertices ( v );
  }


///Gets vertices in an user allocated array
  void Get_Vertices ( Vertex_Iterator v[] ) const
  {
    for ( int i=0; i<nbpts; ++i )
    {
      v[i]=vertex[i];
    }
  }

///Returns # of vertices
  const int Get_Number_Vertices ( void ) const
  {
    return nbpts;
  }

///sets vertices from an user allocated array
  void Set_Vertices ( Vertex_Iterator v[] )
  {
    for ( int j=0; j<nbpts; ++j )
    {
      vertex[j]=v[j];
    }
  }

//le code suivant est debile tant que les noeuds ne sont pas classes dans le simplexe
//a coder differemnt si besoin est (en cas de classification des simplexes dans une structure de donnees type "set")

  /*  bool operator == (const Simplex_Base<nbpts,V>& S) const
      {
          bool tmp=true;
          for (int i=0;i<nbpts;++i)
              if (vertex[i]!=S.vertex[i])
              {
                  tmp=false;
                  break;
              }
          return tmp;
      }
  */
  /*  bool operator < (const Simplex_Base<nbpts,V>& S) const
      {
          for (int i=0;i<nbpts;++i)
          {
              if (!(&(*(vertex[i]))<&*(S.vertex[i]))) return false;
              if (!(vertex[i]==S.vertex[i])) return true;
          }
          return false;
      }
  */

};


template<int nbpts,int dim,class V>
double Simplex_Base<nbpts,dim,V>::Normal ( Vector &norm )
{
  int delta=V::Dimension- ( nbpts-1 );
  int i;

  switch ( delta )
  {
    case 1 :
    {
      switch ( nbpts )
      {
        case 2:
        {
          Vector V1 ( dim,*vertex[0],*vertex[1] );
          norm[0]=-V1[1];
          norm[1]=V1[0];
          norm.Normalize();
          break;
        }

        case 3:
        {
          Vector V1 ( dim,*vertex[0],*vertex[1] );
          Vector V2 ( dim,*vertex[0],*vertex[2] );
          norm.Cross ( V1,V2 );
          norm.Normalize();
          break;
        }

        default:
          for ( i=0; i<V::Dimension; ++i )
          {
            norm[i]=0;
          }
      }

      break;
    }

    default :
      for ( i=0; i<V::Dimension; ++i )
      {
        norm[i]=0;
      }
  }

  return 0.0;
}


#endif // SIMPLEX_BASE_H