// ORIG-DATE:     Fev 2010
// -*- Mode : c++ -*-
//
// SUMMARY  : liaison medit freefem++ : adaptmesh in 3d 
// USAGE    : LGPL      
// ORG      : LJLL Universite Pierre et Marie Curie, Paris,  FRANCE 
// AUTHOR   : Jacques Morice
// Modif    : F. hecht  : Frederic.hecht@upmc.fr
// E-MAIL   : jacques.morice@ann.jussieu.fr

//
//   for automatic  compilation with ff-c++
//ff-c++-LIBRARY-dep:   mmg3d-v4
//ff-c++-cpp-dep: 
//  

/* 
 This file is part of Freefem++
 
 Freefem++ is free software; you can redistribute it and/or modify
 it under the terms of the GNU Lesser General Public License as published by
 the Free Software Foundation; either version 2.1 of the License, or
 (at your option) any later version.
 
 Freefem++  is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU Lesser General Public License for more details.
 
 You should have received a copy of the GNU Lesser General Public License
 along with Freefem++; if not, write to the Free Software
 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

 Thank to the ARN ()  FF2A3 grant
 ref:ANR-07-CIS7-002-01 
 */


//  ./ff-c++ mmg3dv4.cpp  -I../download/include/mmg3d/  -lmmg3d4 

#include "ff++.hpp" 
#include "msh3.hpp"
//#define ADAPTLIBRARY"
#include "dataff.h"

using namespace  Fem2D;
//using namespace  mmg3d;

inline void add2(int *k,int n, int a)
{
    for (int i=0;i<n;++i)
      k[i] += a; 
}
void set_mesh(void *dataff,int *data,int ldata)
{
    DataFF *dff=(DataFF *) dataff;
    
    int nnv= data[ff_id_vertex];
    int nnt= data[ff_id_tet];
    int nnbe= data[ff_id_tria];
    Vertex3 *vv = new Vertex3[nnv];
    Tet *tt= new Tet[nnt]; 
    Triangle3 *bb=new Triangle3[nnbe];
    
    Mesh3 * pTh= new Mesh3();
    pTh->nv = nnv;
    pTh->nt = nnt;
    pTh->nbe =nnbe;
    
    pTh->vertices = vv;
    pTh->elements = tt;
    pTh->borderelements = bb;
    
    pTh->mes=0.;
    pTh->mesb=0.;
    
    dff->mesh= pTh; 
    if(verbosity>5) cout << " Set_mesh nv=" << nnv << " nTet " << nnt<< " NTria " <<nnbe << endl;
}
void end_mesh(void *dataff)
{
    DataFF *dff=(DataFF *) dataff;

    Mesh3 &  Th = *(Mesh3 *) dff->mesh;
    
    Th.mes=0.;
    Th.mesb=0.;
    
    for (int i=0;i<Th.nbe;i++)  
	Th.mesb += Th.be(i).mesure();  

    //  Add FH to be consitant we all constructor ...  July 09
    Th.BuildBound();
    //Th.Save("TTTh.mesh");
    if(verbosity>5)
      cout << "end_mesh:  Th.mes = " << Th.mes << " Th.mesb = " << Th.mesb << endl;
    if(Th.nt > 0){ 
	Th.BuildAdj();
	Th.Buildbnormalv();  
	Th.BuildjElementConteningVertex();  
    }
    //  end add       
    
    if(verbosity>1)
	cout << "  -- End of Construct  mesh3: mesure = " << Th.mes << " border mesure " << Th.mesb <<  endl;  
    ffassert(Th.mes>=0); // add F. Hecht sep 2009.
        
}
void set_v(void *dataff,int i,double *xyz,int lab)
{
  i--; 
  DataFF *dff=(DataFF *) dataff;
  ffassert(dff->mesh);
  Mesh3 &  Th= * (Mesh3 *) dff->mesh;
  Th.vertices[i].x = xyz[0];
  Th.vertices[i].y = xyz[1];
  Th.vertices[i].z = xyz[2];
  Th.vertices[i].lab = lab;
  if(verbosity>10) 
    cout << " set_v3 " <<i << " " <<  xyz[0] << " "  << xyz[1] << " "  << xyz[02] << " "  << lab  << endl;
  
}
void set_elmt(void *dataff,int id,int i,int *k,int lab)
{
  i--;
  int n=0;
    DataFF *dff=(DataFF *) dataff;
    ffassert(dff->mesh);
    Mesh3 &  Th= *(Mesh3 *) dff->mesh;
      if( id == 2) 
      {
	n=3;
	Mesh3::BorderElement & K(Th.be(i));
	add2(k,3,-1);		  
	K.set(Th.vertices,k,lab);
      }
	
    else if ( id== 3) 
      {
	n=4; 
	Mesh3::Element & K(Th.t(i));
	add2(k,4,-1);
	K.set(Th.vertices,k,lab);
      }
    else
      {
	cout << " unknows id = " << id << " not 2 or 3 " << endl;
	ffassert(0); 
      }

  if(verbosity>10) 
    {
      cout << " set_ele"<< n << " "  <<i << " "  ;
      for(int j=0;j<n;j++)
	cout << k[j] << " ";
      cout << lab << endl;
    }

}
void get_mesh(void *dataff,int *data,int ldata)
{
    DataFF *dff=(DataFF *) dataff;
    assert(ldata>5);
    for(int i=0;i<ldata;++i)
	data[i]=0;
    ffassert(dff->mesh);
    Mesh3 &  Th= *(Mesh3*) dff->mesh;
    data[ff_id_vertex]=  Th.nv;
    data[ff_id_tria]=  Th.nbe;
    data[ff_id_tet]=  Th.nt;
    if(verbosity>9) 
      cout << " get_mesh " << Th.nv << " "<< Th.nbe << " "<< Th.nt << endl; 
    
}
void get_v3(void *dataff,int i,double *xyz,int *lab)
{
  i--;
  DataFF *dff=(DataFF *) dataff;
  ffassert(dff->mesh);
  Mesh3 &  Th= *(Mesh3*) dff->mesh;
  xyz[0] = Th.vertices[i].x ;
  xyz[1] = Th.vertices[i].y ;
  xyz[2] = Th.vertices[i].z ;
  *lab = Th.vertices[i].lab ;
  if(verbosity>10) 
    cout << " get_v3 " <<i << " " <<  xyz[0] << " "  << xyz[1] << " "  << xyz[2] << " "  << *lab  << endl;
    
  
}
void get_elmt(void *dataff,int id,int i,int *k,int *lab)
{
  i--;
  DataFF *dff=(DataFF *) dataff;
  ffassert(dff->mesh);
  Mesh3 &  Th= *(Mesh3*)dff->mesh;
  int n =0;
  if( id == 2) 
    {
      n=3; 
      Mesh3::BorderElement & K(Th.be(i));
      for (int j=0; j<n;++j)
	k[j]=Th(K[j]); 
      
      *lab = K.lab;
    }
  
  else if ( id== 3) 
    {
      n=4;
	Mesh3::Element & K(Th.t(i));
        for (int j=0; j<n;++j)
	  k[j]=Th(K[j]); 	
	*lab=K.lab;
    } 
  else {
    cout << "  id != 2, 3 , id = = "<< id  << endl;  
    ffassert(0); 
  }
  add2(k,n,+1);
  
  if(verbosity>10) 
    {
      cout << " get_ele"<< n << " "  <<i << " "  ;
      for(int j=0;j<n;j++)
	cout << k[j] << " ";
      cout << *lab << endl;
    }
  
}



class mmg3d_Op: public E_F0mps 
{
public:
  Expression eTh,xx,yy,zz;
  static const int n_name_param = 5; // 
  static basicAC_F0::name_and_type name_param[] ;
  Expression nargs[n_name_param];
  KN_<long>   karg(int i,Stack stack) const
  { return nargs[i] ? GetAny<KN_<long> >( (*nargs[i])(stack) ): KN_<long>((long*) 0,0L) ;}

  double  arg(int i,Stack stack,double a ) const{ return nargs[i] ? GetAny< double >( (*nargs[i])(stack) ): a;}
  string  arg(int i,Stack stack,const char * a ) const{ return nargs[i] ? *GetAny< string * >( (*nargs[i])(stack) ): a;}
  long  arg(int i,Stack stack, long a ) const{ return nargs[i] ? GetAny< long >( (*nargs[i])(stack) ): a;}
  
  
public:
  mmg3d_Op(const basicAC_F0 &  args ,Expression tth) 
    : eTh(tth),xx(0),yy(0),zz(0) 
  {
    if(verbosity >1) 
      cout << "mmg3d v4 "<< endl;
    args.SetNameParam(n_name_param,name_param,nargs);
    
    const E_Array * a1=0 ;
    if(nargs[1])  a1  = dynamic_cast<const E_Array *>(nargs[1]);
  
    if(a1) {
      if(a1->size() !=3) 
	CompileError("mmg3d(Th,displacement=[X,Y,Z],) ");
      xx=to<double>( (*a1)[0]); 
      yy=to<double>( (*a1)[1]);
      zz=to<double>( (*a1)[2]);
    }    
    else if ( nargs[1] ) 
      CompileError("mmg3d(Th,displacement=[X,Y,Z], .... ) ");
  } 
  
  AnyType operator()(Stack stack)  const ;
};


basicAC_F0::name_and_type  mmg3d_Op::name_param[]= {
  {  "metric", &typeid(KN<double> *)},   // 0
  {  "displacement", &typeid(E_Array)},  // 1
  {  "displVect", &typeid(KN_<double>)},   // 2
  {  "opt", &typeid(string*)}, // 3
  {  "Mb",&typeid(long)} // 4
};

class mmg3d_ff : public OneOperator { public:  
     mmg3d_ff() : OneOperator( atype<pmesh3>(), atype<pmesh3>() ) {}
  
  E_F0 * code(const basicAC_F0 & args) const 
  { 
	return  new  mmg3d_Op( args,t[0]->CastTo(args[0]) ); 
  }
};

AnyType mmg3d_Op::operator()(Stack stack)  const 
{
  // initialisation
  MeshPoint *mp(MeshPointStack(stack)) , mps=*mp;
  Mesh3 * pTh= GetAny<Mesh3 *>((*eTh)(stack));
  ffassert( pTh );
  Mesh3 &Th3=*pTh;
    string sarg= arg(3,stack,"");
    DataFF dff;
    dff.memory= arg(4,stack,128L); // 128 Mb .. ????  
    ffassert( dff.memory < 2048 );// 2 GiGa bytes   limite of integer .. 
    dff.typesol=0;
    dff.np=pTh->nv;
    dff.mesh=pTh;
    dff.meshname="Th";
    dff.imprim=verbosity;
    dff.sol=0;
    dff.mov=0;
    dff.set_mesh =set_mesh;
    dff.end_mesh =end_mesh;
    dff.set_v = set_v;
    dff.set_elmt = set_elmt;
    dff.get_mesh = get_mesh;
    dff.get_v3 = get_v3;
    dff.get_elmt= get_elmt;

  KN<double> *pmetric=0;
  
  if( nargs[0] ){ 
    pmetric = GetAny< KN<double> *>( (*nargs[0])(stack) );
    ffassert(pmetric) ;
  }

  KN<double>  cmetric; 
  if(pmetric)
    {
      int m=pmetric->N();
      if( m == Th3.nv*6) // 
	{
	  cmetric = (*pmetric);
	  dff.typesol=6;
	  dff.np = Th3.nv;
	  for(int i=0;i <m; i += 6)
	    std::swap(cmetric[i+2],cmetric[i+3]); 
	  dff.sol= &cmetric[0];
	  dff.solname="metrix-aniso";
	}
      else if ( m == Th3.nv) 
	{
	  dff.typesol = 1;
	  dff.np = Th3.nv;
	  dff.sol = &(*pmetric)[0];
	  dff.solname="metrix-iso";
	}
      else
	ExecError(" mmg3d v4: incompatibility  metric array  mesh ");
    }


 
 
  bool BoolMoving=0;
  KN<double> Moving(0);
  
  if( nargs[1] || nargs[2] ){
    BoolMoving=1;
    if( nargs[2] ){
      Moving = GetAny<double>( (*nargs[2])(stack) );
      assert( Moving.N() == 3*Th3.nv );
      dff.movename="move";
      dff.mov = Moving;
      if( Moving.N() != 3*Th3.nv ){ cerr << " Displacement vector is of size 3*Th.nv" << endl;
	ExecError(" mmg3d v4");} 
    }
    else{ 
      MeshPoint *mp3(MeshPointStack(stack));
      Moving.resize(3*Th3.nv);
      for( int i=0; i<Th3.nv; ++i){
	mp3->set( Th3.vertices[i].x, Th3.vertices[i].y, Th3.vertices[i].z );
	if(xx) Moving[3*i]   = GetAny<double>((*xx)(stack)); 
	if(yy) Moving[3*i+1] = GetAny<double>((*yy)(stack));
	if(zz) Moving[3*i+2] = GetAny<double>((*zz)(stack));  
      }
      dff.mov = Moving;
      dff.movename="move";
    }
    //if(verbosity > 2) 
    if(verbosity >2) cout << "displacement vector is loading" << endl;
  }

    int argc=1;
    char * argv[1000];
    char ff[10]="ff++";
    KN<char> args(sarg.size()+1);
    argv[0] = ff; 
    argv[1] = & args[0]; 
    strcpy(args, sarg.c_str());
    //cout << sarg << " == " << &args[0] << endl;
    char cc='\0';
    for(int i=0;i<args.N();cc=args[i++])
      {
	if (isspace(args[i]) &&  cc!='\\' ) cc=args[i]='\0';		
	else if ( !cc ) argv[argc++]=& args[i];	
	ffassert(argc<1000);
      }
    
 //   for(int i=0;i< argc; ++i)
//	cout << " *** arg " << i << " " << argv[i] << endl;
	    		
  int res= mainmmg3d(argc,argv,&dff);
  Mesh3 * pTh3 = (Mesh3 *) dff.mesh;
  if( res > 0){
      dff.mesh=0;
    cout << " problem of remeshing with mmg3d :: error" <<  res << endl; 
  }

  if(! pTh3 ){
	  cout << " problem of remeshing with mmg3d v 4 (no mesh)  :: error" <<  res << endl; 
      ExecError(" Error mmg3d" );} 
  else {
      //  end build of TH3...
      if(verbosity > 10) cout << "buildGtree" << endl;
      pTh3->BuildGTree();
  }
  
  *mp=mps;
  Add2StackOfPtr2FreeRC(stack,pTh3);
  return pTh3;
}


/*  class Init1 { public:
  Init1();
};

$1 */

static void Load_Init(){  // le constructeur qui ajoute la fonction "splitmesh3"  a freefem++ 
  
  //if (verbosity)
  if(verbosity) cout << " load: mmg3d  " << endl;
  
  Global.Add("mmg3d","(",new mmg3d_ff);
  
}


#define  WITH_NO_INIT
#include "msh3.hpp" 
LOADFUNC(Load_Init)