/*!\file Bamgx
 * \brief: use Bamg capabilities.
 */
#include "./Bamgx.h"
#include "../../classes/bamg/bamgobjects.h"
#include "../../shared/shared.h"
#include "../../include/include.h"
#include "../../toolkits/toolkits.h"

using namespace bamg;
using namespace std;

int Bamgx(BamgMesh* bamgmesh_out,BamgGeom* bamggeom_out,BamgMesh* bamgmesh_in,BamgGeom* bamggeom_in,BamgOpts* bamgopts){

	/*Bamg options*/
	int    maxnbv;
	double coef;
	int    verbosity;
	int    nbsmooth;

	/*intermediary*/
	int i;
	int noerr=1;
	double costheta=2;
	double hminaniso=1e-100; 
	Mesh* Thr=NULL;
	Mesh* Thb=NULL;

	/*Bamg options*/
	nbsmooth =bamgopts->nbsmooth;
	coef     =bamgopts->coeff;
	maxnbv   =bamgopts->maxnbv;
	verbosity=bamgopts->verbose;

	// no metric -> no smoothing 
	if (bamgopts->metric==NULL) nbsmooth=0; 

	/*If no mesh in input, generate one*/
	if(bamgmesh_in->TrianglesSize[0]==0){
		/*Mesh generation {{{*/

		//Step1: generate geometry Gh
		if (verbosity>0) _printLine_("Construction of a mesh from a given geometry");
		if (verbosity>1) _printLine_("   Processing geometry...");
		Geometry Gh(bamggeom_in,bamgopts);

		//get hmin and hmax from geometry to generate the metric
		bamgopts->hmin = Max(bamgopts->hmin,Gh.MinimalHmin());
		bamgopts->hmax = Min(bamgopts->hmax,Gh.MaximalHmax());

		//build metric using geometry
		if (verbosity>1) _printLine_("   Generating Metric...");
		for(i=0;i<Gh.nbv;i++){
			Metric M=Gh[i];
			EigenMetric Vp(M/coef);
			Vp.Maxh(bamgopts->hmax);
			Vp.Minh(bamgopts->hmin);
			Gh.vertices[i].m = Vp;
		}

		//generate mesh
		if (verbosity>1) _printLine_("   Generating Mesh...");
		Mesh Th(maxnbv,Gh,bamgopts);

		//Split corners if requested
		if(bamgopts->splitcorners) Th.SplitInternalEdgeWithBorderVertices();

		//Renumbering
		Th.TrianglesRenumberBySubDomain();

		//Crack mesh if requested
		if(bamgopts->Crack) Th.CrackMesh(bamgopts);

		//Build output
		if (verbosity>1) _printLine_("   Write Mesh...");
		Th.WriteMesh(bamgmesh_out,bamgopts);
		if (verbosity>1) _printLine_("   Write Geometry...");
		Gh.WriteGeometry(bamggeom_out,bamgopts);

		//clean up
	//	delete &Th;
	//	delete &Gh;
		/*}}}*/
	}
	else{
		/*Anisotropic mesh adaptation {{{*/

		// read background mesh 
		if (verbosity>0) _printLine_("Anisotropic mesh adaptation");
		if (verbosity>1) _printLine_("   Processing initial mesh and geometry...");
		Mesh BTh(bamggeom_in,bamgmesh_in,bamgopts); 

		//Make Quadtree from background mesh
		BTh.MakeBamgQuadtree();

		//Bound hmin and hmax
		bamgopts->hmin=Max(bamgopts->hmin,BTh.MinimalHmin());
		bamgopts->hmax=Min(bamgopts->hmax,BTh.MaximalHmax());

		//Generate initial metric
		if (bamgopts->metric){
			if (verbosity>1) _printLine_("   Processing Metric...");
			BTh.ReadMetric(bamgopts);
		}
		else { 
			if (verbosity>1) _printLine_("   Generating initial Metric...");
			Metric Mhmax(bamgopts->hmax);
			for (int iv=0;iv<BTh.nbv;iv++) BTh[iv].m = Mhmax;
		}

		//use present fields to generate metric if present
		if (bamgopts->field){
			if (verbosity>1) _printLine_("   Merge metric with field provided...");
			BTh.AddMetric(bamgopts);
		}

		// change using hVertices if provided
		if(bamgopts->hVertices && bamgopts->hVerticesSize[0]==BTh.nbv){
			if (verbosity>1) _printLine_("   Merging Metric with hVertices...");
			for (i=0;i<BTh.nbv;i++){
				if (!xIsNan<IssmPDouble>(bamgopts->hVertices[i])){
					BTh[i].m=Metric((float)bamgopts->hVertices[i]);
				}
			}
		}

		// change using hminVertices if provided
		if (bamgopts->hminVertices){
			if (verbosity>1) _printLine_("   Merging Metric with hminVertices...");
			for (i=0;i<BTh.nbv;i++){
				if (!xIsNan<IssmPDouble>(bamgopts->hminVertices[i])){
					Metric M=BTh.vertices[i].m;
					EigenMetric Vp(M/coef);
					Vp.Minh(bamgopts->hminVertices[i]);
					BTh.vertices[i].m=Vp;
				}
			}
		}

		// change using hmaxVertices if provided
		if (bamgopts->hmaxVertices){
			if (verbosity>1) _printLine_("   Merging Metric with hmaxVertices...");
			for (i=0;i<BTh.nbv;i++){
				if (!xIsNan<IssmPDouble>(bamgopts->hmaxVertices[i])){
					Metric M=BTh.vertices[i].m;
					EigenMetric Vp(M/coef);
					Vp.Maxh(bamgopts->hmaxVertices[i]);
					BTh.vertices[i].m=Vp;
				}
			}
		}

		//Add geometry metric if provided
		if(bamgopts->geometricalmetric) BTh.AddGeometryMetric(bamgopts);

		//Smoothe metric
		BTh.SmoothMetric(bamgopts->gradation);

		//Control element subdivision
		BTh.MaxSubDivision(bamgopts->maxsubdiv);

		//Bound anisotropy
		BTh.BoundAnisotropy(bamgopts->anisomax,hminaniso);

		//Build new mesh
		if (verbosity>1) _printLine_("   Generating Mesh...");
		Thr=&BTh,Thb=0;
		Mesh & Th( *(0 ?  new Mesh(*Thr,&Thr->Gh,Thb,maxnbv) :  new Mesh(maxnbv,BTh,bamgopts,bamgopts->KeepVertices)));
		if (Thr != &BTh) delete Thr;

		//Make quadrangle if requested
		if(costheta<=1.0) Th.MakeQuadrangles(costheta);
		//if () Th.SplitElement(2);

		//Split corners if requested
		if(bamgopts->splitcorners) Th.SplitInternalEdgeWithBorderVertices();

		//Renumber by subdomain
		Th.TrianglesRenumberBySubDomain();

		//Smooth vertices
		if(nbsmooth>0) Th.SmoothingVertex(nbsmooth,bamgopts->omega);

		//display info
		if(verbosity>0) {
			if (Th.nbt-Th.nbtout-Th.nbq*2){
				_printLine_("   new number of triangles = " << (Th.nbt-Th.nbtout-Th.nbq*2));
			}
			if (Th.nbq ){
				_printLine_("   new number of quads = " << Th.nbq);
			}
		}

		//Build output
		if (verbosity>1) _printLine_("   Write Mesh...");
		Th.WriteMesh(bamgmesh_out,bamgopts);
		if (verbosity>1) _printLine_("   Write Geometry...");
		Th.Gh.WriteGeometry(bamggeom_out,bamgopts);
		if (verbosity>1) _printLine_("   Write Metric...");
		BTh.WriteMetric(bamgopts);

		/*clean up*/
		delete &Th;
		//delete &BTh;
		/*}}}*/
	}

	/*No error return*/
	if (verbosity>1) _printLine_("   Exiting Bamg.");
	return noerr;

}
