 	}
 }
 
+void CheckCUDAError(const char *msg)
+{
+	cudaError_t err = cudaGetLastError();
+	if (cudaSuccess != err) {
+		fprintf(stderr, "Cuda error: %s: %s.\n", msg, cudaGetErrorString(err));
+		exit(EXIT_FAILURE);
+	}
+
+}
 __device__
 double MyAtomicAdd(double* address, double val)
 {

	{
		fout << i << ':' << host_test_polygon[i] << endl;
	}*/
}
	free(seeds_neighbors);
	free(polygon_buffer);
}

__global__
void compute_RVD_with_knn(double* seeds_pointer, int seeds_nb,
double* mesh_vertex, int mesh_vertex_nb,
int* mesh_facet, int mesh_facet_nb,
int* facet_center_neighbor_index, int* seeds_neighbor_index, int k, double* ret_seeds)
{
	int tid = blockIdx.x * blockDim.x + threadIdx.x;
	if (tid >= mesh_facet_nb) return;

	int* k_indexes = &facet_center_neighbor_index[tid * k];

	int f_idx1 = mesh_facet[tid * 3 + 0];
	int f_idx2 = mesh_facet[tid * 3 + 1];
	int f_idx3 = mesh_facet[tid * 3 + 2];

	int3 facet_index = { f_idx1, f_idx2, f_idx3 };

	double3 v1 = { mesh_vertex[facet_index.x * 3 + 0],
		mesh_vertex[facet_index.x * 3 + 1],
		mesh_vertex[facet_index.x * 3 + 2] };

	double3 v2 = { mesh_vertex[facet_index.y * 3 + 0],
		mesh_vertex[facet_index.y * 3 + 1],
		mesh_vertex[facet_index.y * 3 + 2] };

	double3 v3 = { mesh_vertex[facet_index.z * 3 + 0],
		mesh_vertex[facet_index.z * 3 + 1],
		mesh_vertex[facet_index.z * 3 + 2] };

	double* polygon = (double*)malloc(sizeof(double) * 10 * 4);
	int vertex_nb;

	/*
		get the nearest k points' index
		clip the facet
		*/
	for (int i = 0; i < k; ++i)
	{
		int current_seed = facet_center_neighbor_index[k * tid + i];
		vertex_nb = 3;

		//intialize the polygon with the 3 vertex of the current facet
		/*
		polygon pointer can be apart by several vertex
		a vertex is made up with x,y,z,w
		w : the weight of a vertex
		*/
		polygon[0] = v1.x; polygon[1] = v1.y; polygon[2] = v1.z; polygon[3] = 1.0;
		polygon[4] = v2.x; polygon[5] = v2.y; polygon[6] = v2.z; polygon[7] = 1.0;
		polygon[8] = v3.x; polygon[9] = v3.y; polygon[10] = v3.z; polygon[11] = 1.0;

		intersection_clip_facet_with_knn(polygon, vertex_nb, current_seed, seeds_pointer, seeds_nb, seeds_neighbor_index, k);

		//now we get the clipped polygon stored in "polygon"
		//take care of the sychonize
		//change the polygon data into "weight" and "position"
		if (vertex_nb == 0)
			break;

		double weight;
		double3 position;

		int _v1 = 0;
		int _v2, _v3;

		double3 pos1, pos2, pos3;
		double d1, d2, d3;
		int triangle_nb = vertex_nb - 2;

		double total_weight = 0.0;
		double3 centriodTimesWeight = { 0.0, 0.0, 0.0 };

		double current_weight = 0.0;
		double3 current_posTimesWeight = { 0.0, 0.0, 0.0 };

		for (int i = 1; i < vertex_nb - 1; ++i)
		{
			_v2 = i; _v3 = i + 1;

			pos1 = { polygon[_v1 * 4 + 0], polygon[_v1 * 4 + 1], polygon[_v1 * 4 + 2] };
			d1 = polygon[_v1 * 4 + 3];

			pos2 = { polygon[_v2 * 4 + 0], polygon[_v2 * 4 + 1], polygon[_v2 * 4 + 2] };
			d2 = polygon[_v2 * 4 + 3];

			pos3 = { polygon[_v3 * 4 + 0], polygon[_v3 * 4 + 1], polygon[_v3 * 4 + 2] };
			d3 = polygon[_v3 * 4 + 3];

			computeTriangleCentriod(pos1, pos2, pos3, d1, d2, d3, centriodTimesWeight, total_weight);

			current_weight += total_weight;
			current_posTimesWeight.x += centriodTimesWeight.x;
			current_posTimesWeight.y += centriodTimesWeight.y;
			current_posTimesWeight.z += centriodTimesWeight.z;

			total_weight = 0.0;
			centriodTimesWeight = { 0.0, 0.0, 0.0 };
		}

		//atomicAdd(&SeedsPolygon_nb[current_seed], 1);
		if (current_weight != 0 && triangle_nb > 0){
			current_posTimesWeight.x /= current_weight;
			current_posTimesWeight.y /= current_weight;
			current_posTimesWeight.z /= current_weight;

			current_weight /= triangle_nb;

			double3 temp_pos;

			temp_pos.x = current_posTimesWeight.x * current_weight;
			temp_pos.y = current_posTimesWeight.y * current_weight;
			temp_pos.z = current_posTimesWeight.z * current_weight;

			MyAtomicAdd(&SeedsInformation[current_seed * 4 + 0], temp_pos.x);
			MyAtomicAdd(&SeedsInformation[current_seed * 4 + 1], temp_pos.y);
			MyAtomicAdd(&SeedsInformation[current_seed * 4 + 2], temp_pos.z);
			MyAtomicAdd(&SeedsInformation[current_seed * 4 + 3], current_weight);
		}

		if (tid == 0)
		{
			for (int m = 0; m < seeds_nb * 4; ++m)
				ret_seeds[m] = m;

			ret_seeds[0] = v1.x;
			ret_seeds[1] = v1.y;
			ret_seeds[2] = v1.z;

			for (int m = 0; m < vertex_nb * 4; ++m)
			{
				ret_seeds[4 + m] = polygon[m];
			}
			ret_seeds[21] = vertex_nb;
		}
		if (i == 1)
		return;
	}
	__syncthreads();

	
	
	for (int i = 0; i < seeds_nb * 4; ++i)
	{
		ret_seeds[i] = SeedsInformation[i];
	}
	
	free(polygon);
	return;
}


/*
	a thread to compute the RVD in parallel way
