/*!\file InterpFromMeshToMesh2dx
 */

#include "./InterpFromMeshToMesh2dx.h"

#include "../../shared/shared.h"
#include "../../include/include.h"
#include "../../toolkits/toolkits.h"
#include "../../objects/objects.h"
#include "../modules.h"

using namespace bamg;
using namespace std;

int InterpFromMeshToMesh2dx(double** pdata_interp,double* index_data,double* x_data,double* y_data,int nods_data,int nels_data,
			double* data,int data_rows,int data_cols,double* x_interp,double* y_interp,int nods_interp,double* default_values,int num_default_values, DataSet* contours){
	
	/*Output*/
	double* data_interp=NULL;

	/*Intermediary*/
	R2     r;
	I2     I;
	int    i,j,k;
	int    it;
	int    i0,i1,i2;
	double areacoord[3];
	double aa,bb;
	double data_value;
	Icoor2 dete[3];
	int verbose=0;

	/*default values: */
	Vector*    vec_incontour=NULL;
	double*    incontour=NULL;
	bool   skip_bamg=false;

	/*Checks*/
	if (data_cols<=0){
		_error_("data provided has a negative number of columns");
	}
	if (data_rows!=nods_data && data_rows!=nels_data){
		_error_("data provided should have either %i or %i lines (not %i)",nods_data,nels_data,data_rows);
	}
	if((num_default_values) && (data_cols>1)){
		_error_("data provided can only have 1 column if a default value is provided");
	}
	
	/*If default values supplied, figure out which nodes are inside the contour, including the border of the contour: */
	if(num_default_values){
		ContourToNodesx( &vec_incontour,x_interp,y_interp,nods_interp,contours,1);
		incontour=vec_incontour->ToMPISerial();
	}

	/*Initialize output*/
	if (verbose) printf("Initializing output vector\n");
	data_interp=(double*)xmalloc(nods_interp*data_cols*sizeof(double));

	// read background mesh 
	if (verbose) printf("Reading mesh\n");
	Mesh Th(index_data,x_data,y_data,nods_data,nels_data); 
	Th.CreateSingleVertexToTriangleConnectivity();

	//Loop over output nodes
	if (verbose) printf("Loop over the nodes\n");
	for(i=0;i<nods_interp;i++){
		
		/*reset skip_bamg: */
		skip_bamg=false;

		/*figure out if we should skip bamg logic: */
		if(num_default_values){
			if(!incontour[i]){
				/*This node is not inside the contour. Skip Bamg logic and apply default value.: */
				skip_bamg=true;
			}
		}

		if(skip_bamg==false){

			//Get current point coordinates
			r.x=x_interp[i]; r.y=y_interp[i];
			I2 I=Th.R2ToI2(r);

			//Find triangle holding r/I
			Triangle &tb=*Th.TriangleFindFromCoord(I,dete);

			// internal point 
			if (tb.det>0){ 
				//Area coordinate
				areacoord[0]= (double) dete[0]/tb.det;
				areacoord[1]= (double) dete[1]/tb.det;
				areacoord[2]= (double) dete[2]/tb.det;
				//3 vertices of the triangle
				i0=Th.GetId(tb[0]);
				i1=Th.GetId(tb[1]);
				i2=Th.GetId(tb[2]);
				//triangle number
				it=Th.GetId(tb);
			}
			//external point
			else {
				//Get closest adjacent triangle (inside the mesh)
				AdjacentTriangle ta=CloseBoundaryEdge(I,&tb,aa,bb).Adj();
				int k=ta;
				Triangle &tc=*(Triangle*)ta;
				//Area coordinate
				areacoord[VerticesOfTriangularEdge[k][1]] = aa;
				areacoord[VerticesOfTriangularEdge[k][0]] = bb;
				areacoord[OppositeVertex[k]] = 1 - aa -bb;
				//3 vertices of the triangle
				i0=Th.GetId(tc[0]);
				i1=Th.GetId(tc[1]);
				i2=Th.GetId(tc[2]);
				//triangle number
				it=Th.GetId(tc);
			}
			
			if (data_rows==nods_data){
				for (j=0;j<data_cols;j++){
					data_interp[i*data_cols+j]=areacoord[0]*data[data_cols*i0+j]+areacoord[1]*data[data_cols*i1+j]+areacoord[2]*data[data_cols*i2+j];
				}
			}
			else{
				for (j=0;j<data_cols;j++){
					if (it<0 || it>=nels_data){
						_error_("Triangle number %i not in [0 %i], because not correctly implemented yet... interpolate on grid first",it,nels_data);
					}
					data_interp[i*data_cols+j]=data[data_cols*it+j];
				}
			}
		}
		else{
			if(num_default_values==1) data_interp[i]=default_values[0];
			else data_interp[i]=default_values[i];
		}
	}

	/*Assign output pointers:*/
	if (verbose) printf("Assigning output\n");
	*pdata_interp=data_interp;

	/*No error return*/
	return 1;
}
