add hilbert curve

Jiaming Guo
1 parent 9191c39e
Showing 1 changed file with 123 additions and 46 deletions Show diff stats
flow.h
@@ -61,7 +61,7 @@ namespace stim {
 #ifdef __CUDACC__
 	// for sphere
 	template <typename T>
-	__global__ void inside_sphere(const stim::sphere<T> *V, unsigned num, size_t *R, T *S, unsigned char *ptr, unsigned x, unsigned y, unsigned z) {
+	__global__ void inside_sphere(const stim::sphere<T> *V, unsigned num, size_t *R, T *S, unsigned char *ptr, int x, int y, int z) {
 		unsigned ix = blockDim.x * blockIdx.x + threadIdx.x;
 		unsigned iy = blockDim.y * blockIdx.y + threadIdx.y;
@@ -91,7 +91,7 @@ namespace stim {
 	// for cone
 	template <typename T>
-	__global__ void inside_cone(const stim::cone<T> *E, unsigned num, size_t *R, T *S, unsigned char *ptr, unsigned x, unsigned y, unsigned z) {
+	__global__ void inside_cone(const stim::cone<T> *E, unsigned num, size_t *R, T *S, unsigned char *ptr, int x, int y, int z) {
 		unsigned ix = blockDim.x * blockIdx.x + threadIdx.x;
 		unsigned iy = blockDim.y * blockIdx.y + threadIdx.y;
@@ -127,7 +127,7 @@ namespace stim {
 	// for source bus
 	template <typename T>
-	__global__ void inside_cuboid(const stim::cuboid<T> *B, unsigned num, size_t *R, T *S, unsigned char *ptr, unsigned x, unsigned y, unsigned z) {
+	__global__ void inside_cuboid(const stim::cuboid<T> *B, unsigned num, size_t *R, T *S, unsigned char *ptr, int x, int y, int z) {
 		unsigned ix = blockDim.x * blockIdx.x + threadIdx.x;
 		unsigned iy = blockDim.y * blockIdx.y + threadIdx.y;
@@ -586,17 +586,6 @@ namespace stim {
 			}
 		}
-		void draw_hilbert_curve(unsigned i, T *c, int order, T l, T d, bool upper, int feeder, bool invert) {
-			
-			T fragment = (d - l) / ((std::pow(4, order) - 1) / (std::pow(2, order) - 1) - 1);	// the length of the opening of cup 		
-			T dl = fragment / (std::pow(2, order) - 1);											// unit cup length
-			
-			if (upper)
-				hilbert_curve(i, c, order, dl, feeder, invert, DOWN);
-			else
-				hilbert_curve(i, c, order, dl, feeder, invert, UP);
-		}
-
 		/// render function
 		// find two envelope caps for two spheres
 		// @param cp1, cp2: list of points on the cap
@@ -1071,7 +1060,7 @@ namespace stim {
 		/// build main feeder connection
 		// set up main feeder and main port of both input and output
-		void set_main_feeder(T border = 200.0f) {
+		void set_main_feeder(T border = 400.0f) {
 			// 0 means outgoing while 1 means incoming
 			stim::vec3<T> inlet_main_feeder;
@@ -1242,11 +1231,62 @@ namespace stim {
 					inlet[i].V.push_back(tmp_v);
 					inlet[i].l = new_l;
-					draw_hilbert_curve(i, &tmp_v[0], order, origin_l, desire_l, upper, 1, invert);
-					if (tmp_v[0] != mid_v[0])
-						inlet[i].V.push_back(mid_v);
-					inlet[i].V.push_back(bus_v);
-					std::reverse(inlet[i].V.begin(), inlet[i].V.end());
+					if (desire_l - origin_l < 2 * radii) {	// do not need to use hilbert curve, just increase the length by draging out
+						T d = new_l - inlet[i].l;
+						stim::vec3<T> corner = stim::vec3<T>(tmp_v[0], tmp_v[1] + d / 2.0f * (tmp_v[1] > main_feeder[0][1] ? 1 : -1), tmp_v[2]);
+						inlet[i].V.push_back(corner);
+						corner = stim::vec3<T>(mid_v[0], mid_v[1] + d / 2.0f * (tmp_v[1] > main_feeder[0][1] ? 1 : -1), mid_v[2]);
+						inlet[i].V.push_back(corner);
+						inlet[i].V.push_back(bus_v);
+					}
+					else {
+						T fragment = (desire_l - origin_l) / ((std::pow(4, order) - 1) / (std::pow(2, order) - 1) - 1);	// the length of the opening of cup 		
+						T dl = fragment / (std::pow(2, order) - 1);											// unit cup length
+
+						if (dl > 2 * radii) {				// if the radii is feasible
+							if (upper)
+								hilbert_curve(i, &tmp_v[0], order, dl, 1, invert, DOWN);
+							else
+								hilbert_curve(i, &tmp_v[0], order, dl, 1, invert, UP);
+
+							if (tmp_v[0] != mid_v[0])
+								inlet[i].V.push_back(mid_v);
+							inlet[i].V.push_back(bus_v);
+						}
+						else {								// if the radii is not feasible
+							int count = 1;
+							while (dl <= 2 * radii) {
+								dl = origin_l / (std::pow(2, order - count) - 1);
+								count++;
+							}
+							count--;
+
+							if (upper)
+								hilbert_curve(i, &tmp_v[0], order - count, dl, 1, invert, DOWN);
+							else
+								hilbert_curve(i, &tmp_v[0], order - count, dl, 1, invert, UP);
+
+							desire_l -= origin_l * ((std::pow(4, order - count) - 1) / (std::pow(2, order - count) - 1));
+							origin_l = (bus_v - mid_v).len();
+							desire_l += origin_l;
+
+							find_hilbert_order(origin_l, desire_l, order);
+
+							fragment = (desire_l - origin_l) / ((std::pow(4, order) - 1) / (std::pow(2, order) - 1) - 1);
+							dl = fragment / (std::pow(2, order) - 1);
+							if (dl < 2 * radii)
+								std::cout << "infeasible connection between inlets!" << std::endl;
+
+							if (upper)
+								hilbert_curve(i, &tmp_v[0], order, dl, 1, !invert, RIGHT);
+							else
+								hilbert_curve(i, &tmp_v[0], order, dl, 1, !invert, RIGHT);
+
+							if (tmp_v[1] != bus_v[1])
+								inlet[i].V.push_back(bus_v);
+						}
+					}
+					std::reverse(inlet[i].V.begin(), inlet[i].V.end());			// from bus to pendant vertex
 				}
 			}
@@ -1265,7 +1305,7 @@ namespace stim {
 			for (unsigned i = 0; i < outlet.size(); i++) {
 				if (i != outlet_index) {
 					new_l = (new_pressure[outlet[i].v[0]] - source_pressure) * ((float)stim::PI * std::pow(radii, 4)) / (8 * viscosity * outlet[i].Q);
-					
+
 					if (outlet[i].V[2][1] > main_feeder[1][1]) {
 						upper = true;
 						invert = true;
@@ -1286,10 +1326,61 @@ namespace stim {
 					outlet[i].V.push_back(tmp_v);
 					outlet[i].l = new_l;
-					draw_hilbert_curve(i, &tmp_v[0], order, origin_l, desire_l, upper, 0, invert);
-					if (tmp_v[0] != mid_v[0])
-						outlet[i].V.push_back(mid_v);
-					outlet[i].V.push_back(bus_v);
+					if (desire_l - origin_l < 2 * radii) {	// do not need to use hilbert curve, just increase the length by draging out
+						T d = new_l - outlet[i].l;
+						stim::vec3<T> corner = stim::vec3<T>(tmp_v[0], tmp_v[1] + d / 2.0f * (tmp_v[1] > main_feeder[0][1] ? 1 : -1), tmp_v[2]);
+						outlet[i].V.push_back(corner);
+						corner = stim::vec3<T>(mid_v[0], mid_v[1] + d / 2.0f * (tmp_v[1] > main_feeder[0][1] ? 1 : -1), mid_v[2]);
+						outlet[i].V.push_back(corner);
+						outlet[i].V.push_back(bus_v);
+					}
+					else {
+						T fragment = (desire_l - origin_l) / ((std::pow(4, order) - 1) / (std::pow(2, order) - 1) - 1);	// the length of the opening of cup 		
+						T dl = fragment / (std::pow(2, order) - 1);											// unit cup length
+
+						if (dl > 2 * radii) {				// if the radii is feasible
+							if (upper)
+								hilbert_curve(i, &tmp_v[0], order, dl, 0, invert, DOWN);
+							else
+								hilbert_curve(i, &tmp_v[0], order, dl, 0, invert, UP);
+
+							if (tmp_v[0] != mid_v[0])
+								outlet[i].V.push_back(mid_v);
+							outlet[i].V.push_back(bus_v);
+						}
+						else {								// if the radii is not feasible
+							int count = 1;
+							while (dl <= 2 * radii) {
+								dl = origin_l / (std::pow(2, order - count) - 1);
+								count++;
+							}
+							count--;
+
+							if (upper)
+								hilbert_curve(i, &tmp_v[0], order - count, dl, 0, invert, DOWN);
+							else
+								hilbert_curve(i, &tmp_v[0], order - count, dl, 0, invert, UP);
+
+							desire_l -= origin_l * ((std::pow(4, order - count) - 1) / (std::pow(2, order - count) - 1));
+							origin_l = (bus_v - mid_v).len();
+							desire_l += origin_l;
+
+							find_hilbert_order(origin_l, desire_l, order);
+
+							fragment = (desire_l - origin_l) / ((std::pow(4, order) - 1) / (std::pow(2, order) - 1) - 1);
+							dl = fragment / (std::pow(2, order) - 1);
+							if (dl < 2 * radii)
+								std::cout << "infeasible connection between outlets!" << std::endl;
+
+							if (upper)
+								hilbert_curve(i, &tmp_v[0], order, dl, 0, !invert, LEFT);
+							else
+								hilbert_curve(i, &tmp_v[0], order, dl, 0, !invert, LEFT);
+
+							if (tmp_v[1] != bus_v[1])
+								outlet[i].V.push_back(bus_v);
+						}
+					}
 					std::reverse(outlet[i].V.begin(), outlet[i].V.end());
 				}
 			}
@@ -1365,31 +1456,17 @@ namespace stim {
 			// secondly push back outside connection
 			for (unsigned i = 0; i < inlet.size(); i++) {
-				if (inlet[i].V.size() == 3) {
-					new_sphere.c = inlet[i].V[1];
+				for (unsigned j = 1; j < inlet[i].V.size() - 1; j++) {
+					new_sphere.c = inlet[i].V[j];
 					new_sphere.r = inlet[i].r;
 					A.push_back(new_sphere);
 				}
-				else {
-					new_sphere.c = inlet[i].V[1];
-					new_sphere.r = inlet[i].r;
-					A.push_back(new_sphere);
-					new_sphere.c = inlet[i].V[2];
-					A.push_back(new_sphere);
-				}
 			}
 			for (unsigned i = 0; i < outlet.size(); i++) {
-				if (outlet[i].V.size() == 3) {
-					new_sphere.c = outlet[i].V[1];
-					new_sphere.r = outlet[i].r;
-					A.push_back(new_sphere);
-				}
-				else {
-					new_sphere.c = outlet[i].V[1];
+				for (unsigned j = 1; j < outlet[i].V.size() - 1; j++) {
+					new_sphere.c = outlet[i].V[j];
 					new_sphere.r = outlet[i].r;
 					A.push_back(new_sphere);
-					new_sphere.c = outlet[i].V[2];
-					A.push_back(new_sphere);
 				}
 			}
@@ -1491,9 +1568,9 @@ namespace stim {
 			unsigned p = 0;																// percentage of progress
 			for (unsigned i = 0; i < size_z; i++) {
-				unsigned x = 0 - (unsigned)Xl;					// translate whole network(including inlet/outlet) to origin
-				unsigned y = 0 - (unsigned)Yb;
-				unsigned z = i + (unsigned)center[2];				// box symmetric along z-axis
+				int x = 0 - (int)Xl;					// translate whole network(including inlet/outlet) to origin
+				int y = 0 - (int)Yb;
+				int z = i + (int)center[2];				// box symmetric along z-axis
 				// allocate image slice memory
 				unsigned char* d_ptr;
 				unsigned char* ptr = (unsigned char*)malloc(num * sizeof(unsigned char));