diff --git a/flow.h b/flow.h index 2de6076..29dabc9 100644 --- a/flow.h +++ b/flow.h @@ -61,7 +61,7 @@ namespace stim { #ifdef __CUDACC__ // for sphere template - __global__ void inside_sphere(const stim::sphere *V, unsigned num, size_t *R, T *S, unsigned char *ptr, unsigned x, unsigned y, unsigned z) { + __global__ void inside_sphere(const stim::sphere *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 - __global__ void inside_cone(const stim::cone *E, unsigned num, size_t *R, T *S, unsigned char *ptr, unsigned x, unsigned y, unsigned z) { + __global__ void inside_cone(const stim::cone *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 - __global__ void inside_cuboid(const stim::cuboid *B, unsigned num, size_t *R, T *S, unsigned char *ptr, unsigned x, unsigned y, unsigned z) { + __global__ void inside_cuboid(const stim::cuboid *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 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 corner = stim::vec3(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(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 corner = stim::vec3(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(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)); -- libgit2 0.21.4