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



#include "mesh_const.h"
#include "metric_field.h"

/// metric composed of other metrics (metric intersection)
template <class V> class Multi_Metric_Field : public Metric_Field<V>
{

  list<Metric_Field<V>* > metrics;
  int num;
public :

///
  Multi_Metric_Field ( int i=5 ) :Metric_Field<V> ( i )
  {
    num=0;
  }

///
  void Add_Metric ( Metric_Field<V>& m )
  {
    list<Metric_Field<V>* >::iterator it;
    it=find ( metrics.begin(),metrics.end(),&m );

    if ( it==metrics.end() )
    {
      metrics.push_back ( &m );
      m.Set_Precision ( Get_Precision() );
      num++;
    }
  }

///
  void Erase_Metric ( Metric_Field<V>& m )
  {
    list<Metric_Field<V>* >::iterator it;
    it=find ( metrics.begin(),metrics.end(),&m );

    if ( it!=metrics.end() )
    {
      metrics.erase ( it );
      num--;
    }
  }

///
  void Clear_Metrics ( void )
  {
    metrics.clear();
    num=0;
  }

///
  virtual void Set_Precision ( int i )
  {
    Metric_Field<V>::Set_Precision ( i );
    list<Metric_Field<V>* >::iterator it;

    for ( it=metrics.begin(); it!=metrics.end(); it++ )
    {
      ( *it )->Set_Precision ( i );
    }
  }


/// metric intersections (works so-so).
  virtual void operator () ( V& Ver,Metric_Tensor& MM )
  {
    list<Metric_Field<V>* >::iterator it=metrics.begin();
    int i,j;

    if ( num==0 )
    {
      MM.Set ( 1 );
    }
    else if ( num==1 )
    {
      ( * ( *it ) ) ( Ver,MM );
    }
    else
    {
      Vector ValP1 ( V::Dimension );
      Vector ValP2 ( V::Dimension );
      Metric_Tensor M2 ( V::Dimension );
      Generic_Matrix VecP1 ( V::Dimension );
      Generic_Matrix VecP2 ( V::Dimension );
      Vector VecPi ( V::Dimension );
      Vector VecPi2 ( V::Dimension );
      Generic_Matrix P ( V::Dimension );
      Generic_Matrix VV ( V::Dimension ),VVV ( V::Dimension );
      double val;

      ( * ( *it ) ) ( Ver,MM );
      VecP1.Eigen_Vectors ( MM,ValP1 );

      while ( ++it!=metrics.end() )
      {
        ( * ( *it ) ) ( Ver,M2 );
        VecP2.Eigen_Vectors ( M2,ValP2 );

        for ( i=0; i<V::Dimension; i++ )
        {
          for ( j=0; j<V::Dimension; j++ )
          {
            VecPi[j]=VecP2 ( j,i );
          }

          M2.Mult ( VecPi,VecPi2 );
          val=VecPi*VecPi2;

          if ( val>ValP1[i] )
          {
            ValP1[i]=val;
          }

          for ( j=0; j<V::Dimension; j++ )
          {
            VecPi[j]=VecP1 ( j,i );
          }

          MM.Mult ( VecPi,VecPi2 );
          val=VecPi*VecPi2;

          if ( val>ValP2[i] )
          {
            ValP2[i]=val;
          }
        }

        VV.Set ( ValP1 );
        P.Mult ( VecP1,VV );
        VecP1.Transpose();  // it's orthogonal since metrics are symmetric and definite positive
        // so M^-1=M^t.
        VV.Mult ( P,VecP1 );
//              for (i=0;i<V::Dimension;i++)
//              {
//                  for (j=i;j<V::Dimension;j++)
//                  {
//                      MM(i,j)=VV(i,j)/2;
//                  }
//              }
        VVV.Set ( ValP2 );
        P.Mult ( VecP2,VVV );
        VecP2.Transpose();
        VVV.Mult ( P,VecP2 );

        if ( VV.Determinant() >VVV.Determinant() )
        {
          for ( i=0; i<V::Dimension; i++ )
          {
            for ( j=i; j<V::Dimension; j++ )
            {
              MM ( i,j ) +=VV ( i,j );
            }
          }

          VecP1.Transpose();
        }
        else
        {
          for ( i=0; i<V::Dimension; i++ )
          {
            for ( j=i; j<V::Dimension; j++ )
            {
              MM ( i,j ) +=VVV ( i,j );
            }
          }

          VecP2.Transpose();
          VecP1=VecP2;
          ValP1=ValP2;
        }

//              VecP1.Eigen_Vectors(MM,ValP1);
      }
    }
  };

///
  virtual double Metric_Distance ( V& V1,V& V2 )
  {
    double tmp;
    double max=0;
    list<Metric_Field<V>* >::iterator it;

    for ( it=metrics.begin(); it!=metrics.end(); it++ )
    {
      tmp= ( *it )->Metric_Distance ( V1,V2 );

      if ( tmp>max )
      {
        max=tmp;
      }
    }

    return max;
  };

///
  virtual ~Multi_Metric_Field()
  {
  }

};

#endif // MULTI_METRIC_FIELD_H