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

#include "gmshNode.h"
#include "gmshElement.h"
#include "FastMarching.h"
#include "SpeedFunctor.h"
#include "ExportFunctor.h"
#include "SPoint3.h"
#include <vector>
#include "GModel.h"
#include "GmshGlobal.h"
#include "SPointToDouble.h"

typedef FM::FastMarching_c<FM::gmshNode,FM::gmshElement,double> FastMarch_c;

class FMSpeedConstant : public FM::SpeedFunctor<FM::gmshElement>
{
public:
  FMSpeedConstant(double sp=1.0) : speed(sp) {}
  virtual double operator() ( const FM::gmshElement &ele ) const
  {
    return speed;
  }
  double speed;
};

class ExportMSH : public FM::ExportFunctor<FM::gmshNode,FastMarch_c::FieldPair>
{
public:
  ExportMSH (  GModel *m ) :f ( m ) {}
  virtual FastMarch_c::FieldPair& operator() ( const FM::gmshNode &N )
  {
    return ff[N.GetVertex()];
  }

  /**
   * Write a "_t.msh" file containing 1st element of FieldPairs
   * and a "_prop.msh" file containing 2nd element of FieldPairs
   * 
   * @param[in] name const char* prefix name of files
  */
  virtual void Write ( const char* name )
  {
    string n1(name),n2(name);
    n1+="_t.msh";
    n2+="_prop.msh";
    f->writeMSH(n1);
    f->writeMSH(n2);
    
    FILE *file1,*file2;
    file1=fopen(n1.c_str(),"a");
    file2=fopen(n2.c_str(),"a");
    fprintf(file1,"$NodeData\n1\n\"LS\"\n1\n0.0\n3\n0\n1\n%lu\n",ff.size());
    fprintf(file2,"$NodeData\n1\n\"LS\"\n1\n0.0\n3\n0\n1\n%lu\n",ff.size());
    map<FM::gmshNode,FastMarch_c::FieldPair>::iterator it,begin,end;
    begin=ff.begin();
    end=ff.end();

    for ( it=begin; it!=end; ++it )
    {
      MVertex *v= ( *it ).first.GetVertex();
      double val1 = ( *it ).second.first;
      double val2 = ( *it ).second.second;
      fprintf(file1,"%lu %f\n",v->getIndex(),val1);
      fprintf(file2,"%lu %f\n",v->getIndex(),val2);
    }
    fprintf(file1,"$EndNodeData\n");
    fprintf(file2,"$EndNodeData\n");
    fclose(file1);
    fclose(file2);
  }
  virtual void Write ( std::map<FM::gmshNode,FastMarch_c::FieldPair> &dest )
  {
    dest=ff;
  }
  
  virtual void Clear ( void )
  {
    ff.clear();
  }
private :
   GModel *f;
  map<FM::gmshNode,FastMarch_c::FieldPair> ff;

};


// init narrowband with distance to point
void InitPt ( GModel *m,SPoint3 pt,map<FM::gmshNode,FastMarch_c::FieldPair> &initvals,set<FM::gmshElement> &bndinit ); 
// init narrowband with distance to analytical iso-0 (and transport arbitrary function)
void Init ( GModel *m,const SPointToDouble& ls,const SPointToDouble& inittrans,map<FM::gmshNode,FastMarch_c::FieldPair> &initvals,set<FM::gmshElement> &bndinit );

class FMField
{
  /// Map liking nodes to its FieldPair containing data
  map<FM::gmshNode,FastMarch_c::FieldPair> lstime;

  FM::ExportFunctor<FM::gmshNode,FastMarch_c::FieldPair>* ef;

public:
  /**
   * Constructor
  */
  FMField(void) {};
  FMField(FastMarch_c &FM_) {Import(FM_);}
  FMField(const FMField &other) : lstime(other.lstime) {}
  virtual void SetExportFunctor ( FM::ExportFunctor<FM::gmshNode,FastMarch_c::FieldPair>& efp )
  {  ef=&efp;  }
  virtual void Export ( const char *name );
  /**
   * lstime of FM_ is copied into lstime of the field (Call Export (FastMarching.h) on FM_ with protected map lstime)
   * 
   * @param[in] FM_ FastMarch_c object whose data must be imported
  */
  virtual void Import (FastMarch_c &FM_ ) {FM_.Export(lstime);}
  /**
   * Copy into the present object the result from the substraction between other1 and other2.
   * IOW : compute other1 - other2 and save into present FMField
   * @param[in] other1 FMField 1st element of substraction
   * @param[in] other2 FMField substracted (2nd element)
   * 
  */
  virtual void Subtract(FMField &other1,FMField &other2);
  /**
   * Copy into the present object the result from the sum between other1 and other2.
   * IOW : compute other1 + other2 and save into present FMField
   * @param[in] other1 FMField 1st element
   * @param[in] other2 FMField 2nd element
   * 
  */
  virtual void Add(FMField &other1,FMField &other2);
  /**
   * Copy into the present object the curvature of other1
   * @param[in] other1 FMField whose curvature is computed
  */
  virtual void RCurvature(FMField &other1);
  virtual void Norm(FM::gmshElement e, SVector3 &ret);
//  virtual void double RCurvature(FM::gmshNode n);
  virtual void Gradient(FM::gmshElement e,FastMarch_c::FieldPair vec[3]);

  /**
   * Interpolate FieldPair at SPoint3 pt :
   * Find for SPoint3 the gmshElement concerned and interpolate. (Point and element unknown)
   * 
   * @param[in] pt SPoint3 to interpolate
   * @return FieldPair. If no element found, field (0,0).
  */
  virtual FastMarch_c::FieldPair operator()(SPoint3 pt);

  /**
   * Interpolate FieldPair at SPoint3 pt :
   * Find at SPoint3 pt in known gmshElement e. (element known, point unknown)
   * 
   * @param[in] e gmshElement where the point must be.
   * @param[in] pt SPoint3 to find.
   * @return FieldPair. If field (0,0), no element found. 
  */
  virtual FastMarch_c::FieldPair operator()(FM::gmshElement e,SPoint3 pt);

  /**
   * Interpolate FieldPair at SPoint3 pt :
   * Get FieldPair at a node. (element unknown, node known)
   * 
   * @param n gmshNode wich is the known point
   * @return FieldPair. If field (0,0), no element found.
  */
  virtual FastMarch_c::FieldPair operator()(FM::gmshNode n);

  /**
   * Look for element containing point (ugly and slow).
   * 
   * @param[in] pt SPoint3 point concerned.
   * @return gmshElement where the point is. Empty gmshElement if not found.
  */
  virtual FM::gmshElement GetElement(SPoint3 pt);
};


class FieldSet
{
public:
  /**
   * Add field to protected vector fields at the end
   * @param[in] f the field (FMField) to add
  */
  virtual void addField(const FMField &f) {fields.push_back(f);}
  /**
   * Clear protected vector fields
  */
  virtual void clearFields(void) {fields.clear();}
  /**
   * Get size of protected vector fields
   * @return int size of fields
  */
  virtual int getNumFields(void) {return fields.size();}
protected :
  /// Vector containing the fields computed
  std::vector<FMField> fields;
};

#endif //FM_FMFIELD_H