/*!\file:  InterpFromMesh2dxt.cpp
 * \brief  thread core for InterpFromMesh2dxt code
 */ 

#include "./InterpFromMesh2dx.h"
#include "../../shared/shared.h"

void* InterpFromMesh2dxt(void* vpthread_handle){

	/*gate variables :*/
	InterpFromMesh2dxThreadStruct* gate=NULL;
	pthread_handle* handle=NULL;
	int     my_thread;
	int     num_threads;
	
	int interpolation_type;
	bool debug;
	int  nels_data;
	double* index_data=NULL;
	double* x_data=NULL;
	double* y_data=NULL;
	double* data=NULL;
	double xmin,xmax;
	double ymin,ymax;
	int    nods_prime;
	Vec    data_prime=NULL;
	double* x_prime=NULL;
	double* y_prime=NULL;
	double* default_values=NULL;
	int     num_default_values;
	double*    incontour=NULL;

	/*intermediary: */
	int     i0;
	int     i1;
	int     i,j;
	double  area;
	double  area_1,area_2,area_3;
	double  data_value;

	/*recover handle and gate: */
	handle=(pthread_handle*)vpthread_handle;
	gate=(InterpFromMesh2dxThreadStruct*)handle->gate;
	my_thread=handle->id;
	num_threads=handle->num;
	
	/*recover parameters :*/
	interpolation_type=gate->interpolation_type;
	debug=gate->debug;
	nels_data=gate->nels_data;
	index_data=gate->index_data;
	x_data=gate->x_data;
	y_data=gate->y_data;
	data=gate->data;
	xmin=gate->xmin;
	xmax=gate->xmax;
	ymin=gate->ymin;
	ymax=gate->ymax;
	nods_prime=gate->nods_prime;
	data_prime=gate->data_prime;
	x_prime=gate->x_prime;
	y_prime=gate->y_prime;
	default_values=gate->default_values;
	num_default_values=gate->num_default_values;
	incontour=gate->incontour;

	/*partition loop across threads: */
	PartitionRange(&i0,&i1,nels_data,num_threads,my_thread);

	/*Loop over the elements*/
	if (debug && my_thread==0) printf("      interpolation progress:   %5.2lf %%",0.0);

	for (i=i0;i<i1;i++){

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

		/*if there is no point inside the domain, go to next iteration*/
		if ( (x_data[(int)index_data[3*i+0]-1]<xmin) && (x_data[(int)index_data[3*i+1]-1]<xmin) && (x_data[(int)index_data[3*i+2]-1]<xmin)) continue;
		if ( (x_data[(int)index_data[3*i+0]-1]>xmax) && (x_data[(int)index_data[3*i+1]-1]>xmax) && (x_data[(int)index_data[3*i+2]-1]>xmax)) continue;
		if ( (y_data[(int)index_data[3*i+0]-1]<ymin) && (y_data[(int)index_data[3*i+1]-1]<ymin) && (y_data[(int)index_data[3*i+2]-1]<ymin)) continue;
		if ( (y_data[(int)index_data[3*i+0]-1]>ymax) && (y_data[(int)index_data[3*i+1]-1]>ymax) && (y_data[(int)index_data[3*i+2]-1]>ymax)) 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(incontour[j]){

				/*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)){
						if(num_default_values==1) data_value=default_values[0];
						else data_value=default_values[j];
					}

					/*insert value and go to the next point*/
					VecSetValue(data_prime,j,data_value,INSERT_VALUES);
				}
			}
		}
	}
	if (debug && my_thread==0) printf("\b\b\b\b\b\b\b%5.2lf %%\n",100.0);
	
	return NULL;

}
