/*
    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.
*/
#include "mesh_globals.h"
#include "linear_algebra.h"
#include <cstdio>
#include <cstring>



//---------------------------------------------------------------------------


char* fortran_str ( double d, char buf[] )
{
  sprintf ( buf,"%0.16E",d );
  int sz=strlen ( buf );

  for ( int i=0; i<sz; i++ )
    if ( buf[i]=='E' )
    {
      buf[i]='D';
    }

  return buf;
}

bool Is_Inside_Triangle3D ( double n[3], double ele[3][3],double param[3] )
{

  int i,j;

  Vector V1 ( 3 ),V2 ( 3 ),SOL ( 3 );
  double Norme;
  Vector rhs ( 3,n );
  Square_Matrix M ( 3 );

  /*
  //ne pas utiliser sous peine de ne pas connaitre le vecteur param.
  //si le point ne se trouve pas dans la boite englobante

  //  1 - test de la boite englobante (rapidos)
      double max[3];
      double min[3];

      for (i=0;i<3;++i)
      {
          max[i]=(ele[0])[i];
          min[i]=max[i];
          for (j=1;j<3;++j)
          {
              double t=(ele[j])[i];
              if (max[i]<t) max[i]=t;
              if (min[i]>t) min[i]=t;
          }
      }

      for (i=0;i<3;++i)
      {
          if (IsLess<double>(n[i],min[i]))
              return(false);
          if (IsLess<double>(max[i],n[i]))
              return(false);
      }
  // 2- si on est effectivement dedans, tester si on est dans le patch (certain et lent)
  */

  for ( i=0; i<3; ++i )
    for ( j=0; j<3; ++j )
    {
      M ( i,j ) =ele[j][i];
    }

  for ( i=0; i<3; ++i )
  {
    rhs[i]-=M ( i,2 );

    for ( j=0; j<3; ++j )
    {
      M ( i,j )-=M ( i,2 );
    }
  }


  for ( i=0; i<3; ++i )
  {
    V1[i]=M ( i,0 );
    V2[i]=M ( i,1 );
  }

  SOL.Cross ( V1,V2 );

  Norme= ( V1.Norm() +V2.Norm() ) /2.0;


  for ( i=0; i<3; ++i )
  {
    M ( i,2 ) =SOL[i]/Norme;
  }

  LU_Matrix LU ( M );
  LU.Solve_Linear_System ( rhs,SOL );

  for ( i=0; i<3; ++i )
  {
    param[i]=SOL[i];
  }

  if (
    ( ( param[0]+param[1] ) <1.0001 ) &&
    ( param[0]>-0.0001 ) &&
    ( param[1]>-0.0001 ) &&
    ( fabs ( param[2] ) <0.0001 )
  )
  {
    return true;
  }
  else
  {
    return false;
  }
}