/*
    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_H
#define SIMPLEX_H
//---------------------------------------------------------------------------

#include "mesh_const.h"

#include <list>

#include "vertex.h"
#include "element.h"
#include "simplex_base.h"
#include "topological_face.h"
#include "bipoint.h"

using namespace std;


/// simplex used in the mesh generator
template<int nbpts,int dim,class V=Vertex<dim> > class Simplex : public Simplex_Base<nbpts,dim,V> , public Element<dim,V>
{


public :
///
  typedef list<Simplex<nbpts,dim,V> > Container_Type;
///
  typedef typename Container_Type::iterator Iterator;
// there is no classification of simplexes
///
  typedef typename Simplex_Base<nbpts,dim,V>::Vertex_Iterator Vertex_Iterator;
  typedef V Vertex_Type;

private :
/// the zone to which the element belongs
  int Zone;


public :

/// init
  Simplex ( void )
  {
    for ( int j=0; j<nbpts; ++j )
    {
      this->vertex[j]= ( Vertex_Iterator ) NULL;
    }

    Zone=0;
  }

/// init with a used specified array of vertices
  Simplex ( Vertex_Iterator v[] )
  {
    Set_Vertices ( v );
    Zone=0;
  }

/// sets the zone to Z
  void Set_Zone ( int Z )
  {
    Zone=Z;
  }

/// gets the zone ID
  int Get_Zone ( void ) const
  {
    return ( Zone );
  }

/// sets the vertices composing the simplex
  void Set_Vertices ( Vertex_Iterator v[] )
  {
    Simplex_Base<nbpts,dim,Vertex_Type>::Set_Vertices ( v );
  }

/// show info
  void Display ( void )
  {
    cout << "El " << this->GetID() << " : ";

    for ( int i=0; i<nbpts; ++i )
    {
      cout << ( * ( this->vertex[i] ) ).GetID() << " " ;
    }

    cout << endl;
  }

  /** splits a simplex along one edge (tab) by a node (cutter)
   and get the resulting element in an user allocated array (tabfill) */
  void Split ( Vertex_Iterator tab[2],Vertex_Iterator cutter,Vertex_Iterator Tabfill[] )
  {
    int j,k=0;
    bool ok;

    for ( int i=0; i<2; ++i )
    {
      ok=false;

      for ( j=0; j<nbpts; ++j )
      {
        if ( tab[i]!=this->vertex[j] )
        {
          Tabfill[k++]=this->vertex[j];
        }
        else
        {
          ok=true;
        }
      }

      if ( ok )
      {
        Tabfill[k++]=cutter;
      }
      else
      {
        throw ( BAD_TOPOLOGY );
      }
    }
  }


  /* gets the generalized circumcircle paramerters
  center is the center, radius is the radius looked through the metric
  in the real space, the circle is an ellipsis, thus the name of the method*/
  void Ellipsoid ( Vertex_Type &Center,double &Radius, Metric_Tensor &M )
  {
    int i,j;
    const int dimension=Vertex_Type::Dimension;

    Square_Matrix SYS ( dimension ),SYSINV ( dimension );
    Vector X[nbpts-1], Dummy ( dimension );
    Vector XX ( dimension );
    Vector RHS ( dimension ),SOL ( dimension );

    for ( i=0; i<nbpts-1; ++i )
    {
      X[i].Resize ( dimension );
    }

    for ( i=0; i<nbpts-1; ++i )
    {

      for ( j=0; j<dimension; ++j )
      {
        XX[j]= ( * ( this->vertex[i] ) ) [j]- ( * ( this->vertex[nbpts-1] ) ) [j];
        ( X[i] ) [j]=XX[j]; // colonnes = vecteurs de base du simplexe
      }

      M.Mult ( XX,Dummy );

      for ( j=0; j<dimension; ++j )
      {
        SYS ( i,j ) =2*Dummy[j];
      }

      RHS[i]=Dummy*XX;
    }

// manque le code pour mettre les n vecteurs de base de l'espace complementaire en N dimensions...
// mais dans le cas 3D surfacique c'est facile, je prend le produit vectoriel...

    if ( ( dimension==3 ) && ( nbpts==3 ) )
    {

      SYS ( 2,0 ) = ( X[0] ) [1]* ( X[1] ) [2]- ( X[0] ) [2]* ( X[1] ) [1];
      SYS ( 2,1 ) = ( X[0] ) [2]* ( X[1] ) [0]- ( X[0] ) [0]* ( X[1] ) [2];
      SYS ( 2,2 ) = ( X[0] ) [0]* ( X[1] ) [1]- ( X[0] ) [1]* ( X[1] ) [0];
      RHS[2]=0;
    }
    else

// mais dans le cas 2D filaire c'est facile aussi

      if ( ( dimension==2 ) && ( nbpts==2 ) )
      {
        SYS ( 1,0 ) =   - ( X[0] ) [1];
        SYS ( 1,1 ) = ( X[0] ) [0];
        RHS[1]=0;
      }
      else

// mais dans le cas ND filaire c'est facile aussi mais moins simple ...

        if ( ( dimension>2 ) && ( nbpts==2 ) )
        {

          double piai=1;

          for ( i=0; i<dimension; ++i )
          {
            piai*= ( X[0] ) [i];
          }

          for ( i=1; i<dimension; ++i )
          {
            for ( j=0; j< ( i-1 ); ++j )
            {
              SYS ( i,j ) = 0.;
            }

            SYS ( i,i-1 ) = -piai/ ( X[0] ) [i-1];

            for ( j=i; j<dimension; ++j )
            {
              SYS ( i,j ) = 0.;
            }

            SYS ( i,dimension-1 ) =piai/ ( X[0] ) [dimension-1];
            RHS[i]=0;
          }
        }

        else if ( dimension!=nbpts-1 )
        {
          throw NOT_IMPLEMENTED;
        }

    SYS.Solve_Linear_System ( RHS,SOL );

    M.Mult ( SOL,Dummy );
    Radius=sqrt ( Dummy*SOL );

    for ( j=0; j<dimension; ++j )
    {
      SOL[j]+= ( * ( this->vertex[nbpts-1] ) ) [j];
      Center[j]=SOL[j];
    }
  }
};

#endif // SIMPLEX_H



/// some useful alias...Segment in 3D space
typedef Simplex<2,3,Vertex<3> > Segment;
/// Triangle in 3D space
typedef Simplex<3,3,Vertex<3> > Triangle;
/// Tetrahedron in 3D space
typedef Simplex<4,3,Vertex<3> > Tetrahedron;