/*!\file:  InterpFromMeshToMesh2dxt.cpp
 * \brief  "thread" core code for interpolating values from a structured grid.
 */ 

#ifdef HAVE_CONFIG_H
	#include "config.h"
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include "./InterpFromMeshToMesh2dx.h"
#include "../shared/shared.h"

#undef __FUNCT__ 
#define __FUNCT__ "InterpFromMeshToMesh2dxt"
void* InterpFromMeshToMesh2dxt(void* vpthread_handle){

	/*gate variables :*/
	InterpFromMeshToMesh2dxThreadStruct* gate=NULL;
	pthread_handle* handle=NULL;
	int     my_thread;
	int     num_threads;
	int     i0;
	int     i1;
	int     debug;
	int     nels_data;
	double* x_data=NULL;
	double* y_data=NULL;
	double* index_data=NULL;
	int     nods_prime;
	double* x_prime=NULL;
	double* y_prime=NULL;
	double* data=NULL;
	int     interpolation_type;
	double  default_value;
	Vec     data_prime;
	double x_prime_min,x_prime_max;
	double y_prime_min,y_prime_max;

	/*local variables: */
	int    i,j;
	double x_tria_min,x_tria_max;
	double y_tria_min,y_tria_max;
	double area;
	double area_1,area_2,area_3;
	double data_value;
	int    step;

	/*recover handle and gate: */
	handle=(pthread_handle*)vpthread_handle;
	gate=(InterpFromMeshToMesh2dxThreadStruct*)handle->gate;
	my_thread=handle->id;
	num_threads=handle->num;

	/*recover parameters :*/
	debug=                gate->debug;
	nels_data=            gate->nels_data;
	x_data=               gate->x_data;
	y_data=               gate->y_data;
	index_data=           gate->index_data;
	nods_prime=           gate->nods_prime;
	x_prime=              gate->x_prime;
	y_prime=              gate->y_prime;
	data=                 gate->data;
	interpolation_type=   gate->interpolation_type;
	default_value=        gate->default_value;
	data_prime=           gate->data_prime;
	x_prime_min=          gate->x_prime_min;
	x_prime_max=          gate->x_prime_max;
	y_prime_min=          gate->y_prime_min;
	y_prime_max=          gate->y_prime_max;

	/*distribute elements across threads :*/
	step=(int)floor((double)nels_data/(double)num_threads);
	for(i=0;i<(my_thread+1);i++){
		i0=i*step;
		if(i==(num_threads-1))i1=nels_data;
		else i1=i0+step;
	}
	
	for (i=i0;i<i1;i++){

		/*display current iteration*/
		if (debug && (my_thread==0) && fmod((double)(i-i0),(double)100)==0) printf("\b\b\b\b\b\b\b%5.2lf %%",(double)(i-i0)/(i1-i0)*100);

		/*Get extrema coordinates of current elements*/
		x_tria_min=x_data[(int)index_data[3*i+0]-1]; x_tria_max=x_tria_min;
		y_tria_min=y_data[(int)index_data[3*i+0]-1]; y_tria_max=y_tria_min;
		
		for (j=1;j<3;j++){
			if(x_data[(int)index_data[3*i+j]-1]<x_tria_min) x_tria_min=x_data[(int)index_data[3*i+j]-1];
			if(x_data[(int)index_data[3*i+j]-1]>x_tria_max) x_tria_max=x_data[(int)index_data[3*i+j]-1];
			if(y_data[(int)index_data[3*i+j]-1]<y_tria_min) y_tria_min=y_data[(int)index_data[3*i+j]-1];
			if(y_data[(int)index_data[3*i+j]-1]>y_tria_max) y_tria_max=y_data[(int)index_data[3*i+j]-1];
		}

		/*if there is no point inside the domain, go to next iteration*/
		if ( x_prime_max < x_tria_min ) continue; 
		if ( x_prime_min > x_tria_max ) continue; 
		if ( y_prime_max < y_tria_min ) continue; 
		if ( y_prime_min > y_tria_max ) continue; 

		/*get area of the current element (Jacobian = 2 * area)*/
		//area =x2 * y3 - y2*x3 + x1 * y2 - y1 * x2 + x3 * y1 - y3 * x1;
		area=x_data[(int)index_data[3*i+1]-1]*y_data[(int)index_data[3*i+2]-1]-y_data[(int)index_data[3*i+1]-1]*x_data[(int)index_data[3*i+2]-1]
		  +  x_data[(int)index_data[3*i+0]-1]*y_data[(int)index_data[3*i+1]-1]-y_data[(int)index_data[3*i+0]-1]*x_data[(int)index_data[3*i+1]-1]
		  +  x_data[(int)index_data[3*i+2]-1]*y_data[(int)index_data[3*i+0]-1]-y_data[(int)index_data[3*i+2]-1]*x_data[(int)index_data[3*i+0]-1];

		/*loop over the prime nodes*/
		for (j=0;j<nods_prime;j++){

			/*if the current point is not in the triangle, continue*/
			if ( x_prime[j] < x_tria_min ) continue; 
			if ( x_prime[j] > x_tria_max ) continue; 
			if ( y_prime[j] < y_tria_min ) continue; 
			if ( y_prime[j] > y_tria_max ) continue; 

			/*Get first area coordinate = det(x-x3  x2-x3 ; y-y3   y2-y3)/area*/
			area_1=((x_prime[j]-x_data[(int)index_data[3*i+2]-1])*(y_data[(int)index_data[3*i+1]-1]-y_data[(int)index_data[3*i+2]-1]) 
						-  (y_prime[j]-y_data[(int)index_data[3*i+2]-1])*(x_data[(int)index_data[3*i+1]-1]-x_data[(int)index_data[3*i+2]-1]))/area;
			/*Get second area coordinate =det(x1-x3  x-x3 ; y1-y3   y-y3)/area*/
			area_2=((x_data[(int)index_data[3*i+0]-1]-x_data[(int)index_data[3*i+2]-1])*(y_prime[j]-y_data[(int)index_data[3*i+2]-1]) 
						- (y_data[(int)index_data[3*i+0]-1]-y_data[(int)index_data[3*i+2]-1])*(x_prime[j]-x_data[(int)index_data[3*i+2]-1]))/area;
			/*Get third area coordinate = 1-area1-area2*/
			area_3=1-area_1-area_2;

			/*is the current point in the current element?*/
			if (area_1>=0 && area_2>=0 && area_3>=0){

				/*Yes ! compute the value on the point*/
				if (interpolation_type==1){
					/*nodal interpolation*/
					data_value=area_1*data[(int)index_data[3*i+0]-1]+area_2*data[(int)index_data[3*i+1]-1]+area_3*data[(int)index_data[3*i+2]-1];
				}
				else{
					/*element interpolation*/
					data_value=data[i];
				}
				if (isnan(data_value)) data_value=default_value;

				/*insert value and go to the next point*/
				VecSetValue(data_prime,j,data_value,INSERT_VALUES);
			}
		}
	}
}