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



#include <cmath>


#include "mesh_const.h"
#include "size_function.h"
#include "functionals.h"
#include "integrators.h"
#include "metric_field.h"

using namespace std;

/// An example of user metric
template <class V> class User_Metric : public Metric_Field<V>
{
///
  double K;
///
  double mu;
///
  double phi;
///
  double Ri;
///
  double Re;
///
  double Pinj;
///
  double dt;
///
  double sztheta;
///
  double desired_dr ( double r )
  {
    double Rm= ( Ri+Re ) /2;

// constant time step
    if ( r<=0 )
    {
      return -1;  // error
    }

    return ( K*Pinj*dt ) / ( mu*phi*r*log ( r/Ri ) );


// constant size
//      return (K*Pinj*dt)/(mu*phi*Rm*log(Rm/Ri));


// constant CFL condition
//      if (r<=0) return -1; // error
//      return (K*Pinj*dt)/(mu*phi*r*log(Rm/Ri));


// constant error
//      return (K*Pinj*dt*r)/(mu*phi*Rm*Rm*log(Rm/Ri));


  }



public :
///
  void Setf ( double RRi,double RRe,double KK,double mmu,double pphi,double PPinj,double ddt,double ssztheta )
  {
    K=KK;
    mu=mmu;
    phi=pphi,Ri=RRi;
    Re=RRe;
    Pinj=PPinj;
    dt=ddt;
    sztheta=ssztheta;
  }

///
  User_Metric ( double RRi,double RRe,double KK,double mmu,double pphi,double PPinj,double ddt,double ssztheta )
  {
    Setf ( RRi,RRe,KK,mmu,pphi,PPinj,ddt,ssztheta );
  }

///
  virtual void operator () ( V& Ver,Metric_Tensor& MM )
  {
    Anti_Symmetric_Matrix R ( V::Dimension );
    Symmetric_Matrix MT ( V::Dimension );

    double x=Ver[0];
    double y=Ver[1];
    double r=sqrt ( x*x+y*y );
    double cteta=x/r;
    double steta=y/r;
    double szr;
    double szt;

    szt=r*sztheta/Ri;
    Metric_Tensor M ( V::Dimension,1. );

    if ( r<=Ri )
    {
      szr= ( szt*3 );
    }
    else
    {
      szr=desired_dr ( r );

      if ( ( szr> ( szt*3 ) ) || ( szr<0 ) )
      {
        szr= ( szt*3 );
      }
    }

    if ( ( szt ) > ( 6*szr ) )
    {
      szt=6*szr;
    }

    M ( 0,0 ) = ( 1./ ( szr*szr ) );
    M ( 1,1 ) = ( 1./ ( szt*szt ) );

    if ( V::Dimension==3 )
    {
      M ( 2,2 ) = ( 1./ ( szr*szr ) );
    }

    // rotation
    R ( 0,0,cteta );
    R ( 1,1,cteta );
    R ( 0,1,steta );

    if ( V::Dimension==3 )
    {
      R ( 2,2,1 );
      R ( 1,2,0 );
      R ( 0,2,0 );
    }

    MT.tRSR ( M,R );

    for ( int i=0; i<V::Dimension; ++i )
      for ( int j=i; j<V::Dimension; ++j )
      {
        MM ( i,j ) =MT ( i,j );
      }
  }



///
  virtual double Metric_Distance ( V& V1,V& V2 )
  {
    Size_Function<V> F ( V1,V2,*this );
//      if (prec==0)
//      {
//          Super_Simple_Integrator<Size_Function<V> > S(F);
//          return S(0.,1.0);
//      }
//      else
//      {
    Simpson_Integrator<Size_Function<V> > S ( F );
    return S ( 0.,1.,10 );
//      }
  }

};


#endif // USER_METRIC_H