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


#include <cmath>
#include "mesh_const.h"
#include "mesh_globals.h"
#include "metric_field.h"
#include "size_function.h"
#include "integrators.h"

using namespace std;

/// an isotropic metric around a cylinder.
template <class V> class Cylindrical_Metric : public Metric_Field<V>
{

///
  double radius;
  double size;
  double thickness;

  V node1;
  V node2;
  Vector vect;


public :
///
  Cylindrical_Metric() :Metric_Field<V> ( 6 ),vect ( V::Dimension )
  {
    int i;
    radius=0;
    size=0;
    thickness=0;

    for ( i=0; i<V::Dimension; ++i )
    {
      node1[i]=0;
      node2[i]=0;
      vect[i]=0;
    }
  }

///
  Cylindrical_Metric ( double rad,double thk, double sz,V & n1,V & n2 ) :Metric_Field<V> ( 6 )
  {
    Set_Params ( rad,thk,sz,n1,n2 );
  }

///
  void Set_Params ( double rad, double thk,double sz, V & n1,V & n2 )
  {
    node1=n1;
    node2=n2;
    radius=rad;
    size=sz;
    thickness=thk;
    vect.Init ( n1,n2 );
  }

///
  virtual void operator () ( V& Ver,Metric_Tensor& MM )
  {
    Vector vect2 ( V::Dimension );
    Vector Vr ( V::Dimension );


    V Proj;
    double pscal,length,length2,r,l,sizel;
    vect2.Init ( node1,Ver );
    pscal=vect2*vect;
    length=vect.Norm();
    length2=vect2.Norm();

    if ( length2==0.0 )
    {
      Proj=node1;
      pscal=0;
      l=0;
    }
    else
    {
      pscal/=length;
      l=pscal;
      pscal/=length;

      for ( int i=0; i<V::Dimension; ++i )
      {
        Proj[i]=node1[i]+vect[i]*pscal;
      }
    }

    Vr.Init ( Proj,Ver );
    r=Vr.Norm();
    l=fabs ( l-length/2 );
    double delt=thickness;

    if ( size>delt )
    {
      delt=size;
    }

    if ( delt>radius )
    {
      delt=radius;
    }

    // now, r and l are the "cylindrical coordinates" in the cylinder.
    // 1- check if r and l are between radius-delt and radius+delt bzw length/2-size and length/2+size;

    if ( ( r<=radius+delt ) && ( r>=radius-delt ) && ( l<=length/2+delt ) )
    {
      sizel=size;
    }
    else

      // 2- inside the cylinder for both r and l
      if ( ( r<=radius-delt ) && ( l<=length/2+delt ) )
      {
        sizel=radius-delt-r+size;
      }
      else

        // 3- outside for r
        if ( ( r>=radius+delt ) && ( l<=length/2+delt ) )
        {
          sizel=r-radius-delt+size;
        }
    // 4- ouside for l
        else if ( ( r>=radius-delt ) && ( r<=radius+delt ) && ( l>=length/2+delt ) )
        {
          sizel=l-length/2-delt+size;
        }
        else

          // 5- outside both l and r
          if ( ( r>=radius+delt ) && ( l>=length/2+delt ) )
          {
            double rr=r- ( radius+delt );
            double ll=l- ( length/2+delt );
            sizel=sqrt ( rr*rr+ll*ll ) +size;
          }
          else

            // 6- outside for l but inside for r
            if ( ( r<=radius-delt ) && ( l>=length/2+delt ) )
            {
              double rr= ( radius-delt )-r;
              double ll=l- ( length/2+delt );
              sizel=sqrt ( rr*rr+ll*ll ) +size;
            }

    MM.Set_Isotropic_Size ( sizel );

  }


///
  virtual double Metric_Distance ( V& V1,V& V2 )
  {

    Size_Function<V> F ( V1,V2,*this );

    switch ( this->prec )
    {
      case 0 :
      {
        Super_Simple_Integrator<Size_Function<V> > S ( F );
        return S ( 0.,1.0 );
        break;
      }

      case 1 :
      {
        Simple_Integrator<Size_Function<V> > S ( F );
        return S ( 0.,1.0 );
        break;
      }

      default :
      {
        Simpson_Integrator<Size_Function<V> > S ( F );
        return S ( 0.,1.,this->prec );
        break;
      }
    }
  }
};


/// an anisotropic metric around a line.
template <class V> class Line_Metric : public Metric_Field<V>
{

///
  double sizet;
  double sizen;
  double thickness;
  V node1;
  V node2;
  Vector vect;

public :
///
  Line_Metric() :Metric_Field<V> ( 6 ),vect ( V::Dimension )
  {
    int i;
    sizet=0;
    sizen=0;

    for ( i=0; i<V::Dimension; ++i )
    {
      node1[i]=0;
      node2[i]=0;
      vect[i]=0;
      vect[i]=0;
    }
  }

///
  Line_Metric ( double szt,double szn,double thk,V & n1,V & n2 ) :Metric_Field<V> ( 6 ),vect ( V::Dimension )
  {
    Set_Params ( szt,szn,thk,n1,n2 );
  }

///
  void Set_Params ( double szt, double szn, double thk,V & n1,V & n2 )
  {
    node1=n1;
    node2=n2;
    sizet=szt;
    sizen=szn;
    thickness=thk;
    vect.Init ( n1,n2 );
  }

///
  virtual void operator () ( V& Ver,Metric_Tensor& MM )
  {
    Vector vect2 ( V::Dimension );
    Vector Vr ( V::Dimension );

    V Proj;
    double pscal,length,length2,r,l,size_n,size_tt,deltn,deltt;
    vect2.Init ( node1,Ver );
    pscal=vect2*vect;
    length=vect.Norm();
    length2=vect2.Norm();

    if ( length2==0.0 )
    {
      Proj=node1;
      pscal=0;
      l=0;
    }
    else
    {
      pscal/=length;
      l=pscal;
      pscal/=length;

      for ( int i=0; i<V::Dimension; ++i )
      {
        Proj[i]=node1[i]+vect[i]*pscal;
      }
    }

    Vr.Init ( Proj,Ver );
    r=Vr.Norm();
    l=fabs ( l-length/2 );

    deltn=thickness;

    if ( sizen>deltn )
    {
      deltn=sizen;
    }

    deltt=deltn;

    if ( sizet>deltt )
    {
      deltt=sizet;
    }

    // now, r and l are the "cylindrical coordinates" with the line.
    if ( ( r<=deltn ) && ( l<=length/2+deltt ) )
    {
      size_n=sizen;
      size_tt=sizet;
    }
    else

      // 2- over or below the line
      if ( ( r>deltn ) && ( l<=length/2+deltt ) )
      {
        size_n=r-deltn+sizen;
        size_tt=r-deltn+sizet;
      }
      else

        // 3- beside the line r<=deltn
        if ( ( r<=deltn ) && ( l>length/2+deltt ) )
        {
          size_n=l-length/2-deltt+sizen;
          size_tt=l-length/2-deltt+sizet;
        }
    // 3- beside the line r>deltn (4 quadrants)
        else if ( ( r>deltn ) && ( l>length/2+deltt ) )
        {
          double rr=r-deltn;
          double ll=l- ( length/2+deltt );
          size_n=sqrt ( rr*rr+ll*ll ) +sizen;
          size_tt=sqrt ( rr*rr+ll*ll ) +sizet;
        }

    Square_Matrix Mbasis ( V::Dimension );
    Square_Matrix EigenValues ( V::Dimension,0.0 );
    Square_Matrix Rot ( V::Dimension ),Rot2 ( V::Dimension ),Evec ( V::Dimension );
    Vector Eval ( V::Dimension );
    int i,j;

    EigenValues ( 0,0 ) =1/ ( size_tt*size_tt );

    for ( i=1; i<V::Dimension; i++ )
    {
      EigenValues ( i,i ) =1/ ( size_n*size_n );
    }

    for ( i=0; i<V::Dimension; i++ )
    {
      Mbasis ( 0,i ) =vect[i];
    }

    Mbasis.Complete_Base ( 1 );
    Rot.Mult ( EigenValues,Mbasis );
    Mbasis.Transpose();
    Rot2.Mult ( Mbasis,Rot );


    for ( i=0; i<V::Dimension; i++ )
      for ( j=i; j<V::Dimension; j++ )
      {
        MM ( i,j ) = ( Rot2 ( i,j ) +Rot2 ( j,i ) ) /2;
      }

    Evec.Eigen_Vectors ( MM,Eval );
  }


///
  virtual double Metric_Distance ( V& V1,V& V2 )
  {

    Size_Function<V> F ( V1,V2,*this );

    switch ( this->prec )
    {
      case 0 :
      {
        Super_Simple_Integrator<Size_Function<V> > S ( F );
        return S ( 0.,1.0 );
        break;
      }

      case 1 :
      {
        Simple_Integrator<Size_Function<V> > S ( F );
        return S ( 0.,1.0 );
        break;
      }

      default :
      {
        Simpson_Integrator<Size_Function<V> > S ( F );
        return S ( 0.,1.,this->prec );
        break;
      }
    }
  }
};


/// an anisotropic metric around a set of lines.
template <class V> class Multiline_Metric : public Metric_Field<V>
{

///
  double sizet;
  double sizen;
  double thickness;
  list<V> vertices;

public:

///
  Multiline_Metric() :Metric_Field<V> ( 6 )
  {
    sizet=0;
    sizen=0;
  }

///
  Multiline_Metric ( double szt,double szn,double thk,V tabV[],int num_pairs ) :Metric_Field<V> ( 6 )
  {
    Set_Params ( szt,szn,thk,tabV,num_pairs );
  }

///
  void Set_Params ( double szt, double szn,double thk, V tabV[],int num_pairs )
  {

    sizet=szt;
    sizen=szn;
    thickness=thk;
    vertices.clear();
    int k=0;

    for ( int i=0; i<num_pairs; i++ )
    {
      vertices.push_back ( tabV[k++] );
      vertices.push_back ( tabV[k++] );
    }
  }

/// returns the distance of the closest segment and its vertices.
  double FindClosest ( V& Ver,V& V1min,V &V2min )
  {
    V V1,V2;
    double d,dmin=-1,nrm,d2,d2min=-1;
    double param[2];

    typename list<V>::iterator it=vertices.begin();

    while ( it!=vertices.end() )
    {
      V1=* ( it++ );
      V2=* ( it++ );
      Vector Vec ( V::Dimension,V1,V2 );
      nrm=Vec.Norm();
      Is_On_Line<V> ( Ver,V1,V2,param );

      if ( ( param[0] ) >1.0 )
      {
        param[0]-=1.0;
      }
      else if ( param[0]<0.0 )
      {
        param[0]*=-1;
      }
      else
      {
        param[0]=0;
      }

      d=sqrt ( param[1]*param[1]+param[0]*param[0] ) *nrm; // norme "2"
      d2=fabs ( param[0] );
//          d=max(fabs(param[1]),fabs(param[0]))*nrm; // norme "1"
//          d=fabs(param[1])*nrm; // norme "infinie"

      if ( ! ( IsEqual ( d,dmin ) ) )
      {
        if ( ( d<dmin ) || ( dmin<0 ) )
        {
          dmin=d;
          d2min=d2;
          V1min=V1;
          V2min=V2;
        }
      }
      else
      {
        if ( ( d2<d2min ) || ( d2min<0 ) )
        {
          dmin=d;
          d2min=d2;
          V1min=V1;
          V2min=V2;
        }
      }
    }

    return max ( dmin,d2min );
  }

///
  virtual void operator () ( V& Ver,Metric_Tensor& MM )
  {
    V V1,V2;
    FindClosest ( Ver,V1,V2 );
    Line_Metric<V> LMetric ( sizet,sizen,thickness,V1,V2 );
    LMetric.Set_Precision ( this->Get_Precision() );
    LMetric ( Ver,MM );
  }


///
  virtual double Metric_Distance ( V& V1,V& V2 )
  {

    Size_Function<V> F ( V1,V2,*this );

    switch ( this->prec )
    {
      case 0 :
      {
        Super_Simple_Integrator<Size_Function<V> > S ( F );
        return S ( 0.,1.0 );
        break;
      }

      case 1 :
      {
        Simple_Integrator<Size_Function<V> > S ( F );
        return S ( 0.,1.0 );
        break;
      }

      default :
      {
        Simpson_Integrator<Size_Function<V> > S ( F );
        return S ( 0.,1.,this->prec );
        break;
      }
    }
  }
};

#endif // METRIC_LCM_TOOLS_H