 	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

	{
		fout << i << ':' << host_test_polygon[i] << endl;
	}*/
}
	}
}

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)
{
