From 8334680aedb2270ba2dbd0924240d0e9fff3ad4b Mon Sep 17 00:00:00 2001 From: Jiaming Guo Date: Tue, 16 May 2017 15:34:56 -0500 Subject: [PATCH] add square wave connections --- flow.h | 657 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- main.cu | 93 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------- 2 files changed, 546 insertions(+), 204 deletions(-) diff --git a/flow.h b/flow.h index 6ba6abc..a9bb9ca 100644 --- a/flow.h +++ b/flow.h @@ -30,7 +30,7 @@ namespace stim { std::vector v; // vertices' indices std::vector > V; // vertices' coordinates T l; // length - T r; // radii + T r; // radius T deltaP; // pressure drop T Q; // volume flow rate }; @@ -38,15 +38,15 @@ namespace stim { template struct sphere { stim::vec3 c; // center of sphere - T r; // radii + T r; // radius }; template - struct cone { // radii changes gradually + struct cone { // radius changes gradually stim::vec3 c1; // center of geometry start hat stim::vec3 c2; // center of geometry end hat - T r1; // radii at start hat - T r2; // radii at end hat + T r1; // radius at start hat + T r2; // radius at end hat }; template @@ -105,7 +105,7 @@ namespace stim { float distance = FLT_MAX; float tmp_distance; - float rr; // radii at the surface where projection meets + float rr; // radius at the surface where projection meets for (unsigned i = 0; i < num; i++) { // find the nearest cylinder tmp_distance = ((world_pixel - E[i].c1).cross(world_pixel - E[i].c2)).len() / (E[i].c2 - E[i].c1).len(); @@ -241,6 +241,7 @@ namespace stim { public: + bool set = false; // flag indicates the pressure has been set std::vector P; // initial pressure std::vector v; // velocity std::vector > main_feeder; // inlet/outlet main feeder @@ -311,7 +312,7 @@ namespace stim { } } - // get the radii of vertex i + // get the radius of vertex i T get_radius(unsigned i) { unsigned tmp_e; // edge index @@ -596,7 +597,7 @@ namespace stim { // find two envelope caps for two spheres // @param cp1, cp2: list of points on the cap // @param center1, center2: center point of cap - // @param r1, r2: radii of cap + // @param r1, r2: radius of cap void find_envelope(std::vector > &cp1, std::vector > &cp2, stim::vec3 center1, stim::vec3 center2, float r1, float r2, GLint subdivision) { stim::vec3 tmp_d; @@ -700,7 +701,7 @@ namespace stim { // calculate the bigger circle d1 = t1 - center1; dot = d1.dot(tmp_d); - a = dot / (r1 * 1) * r1; // a = cos(alpha) * radii + a = dot / (r1 * 1) * r1; // a = cos(alpha) * radius new_c = center1 + a * tmp_d; new_r = sqrt(pow(r1, 2) - pow(a, 2)); new_u = t1 - new_c; @@ -803,7 +804,7 @@ namespace stim { std::vector > cp; T radius; - glColor3f(1.0f, 1.0f, 1.0f); + glColor3f(0.0f, 0.0f, 0.0f); for (unsigned i = 0; i < num_edge; i++) { // for every edge for (unsigned j = 0; j < E[i].size() - 1; j++) { // for every point on current edge tmp_d = E[i][j + 1] - E[i][j]; @@ -812,9 +813,9 @@ namespace stim { tmp_c.rotate(tmp_d); radius = (E[i].r(j + 1) + E[i].r(j)) / 2; if (v[i] > 0) // if flow flows from j to j+1 - head = center + tmp_d * sqrt(3) * radius; + head = center + tmp_d * 2 * sqrt(3) * radius; else - head = center - tmp_d * sqrt(3) * radius; + head = center - tmp_d * 2 * sqrt(3) * radius; stim::circle c(center, radius, tmp_d, tmp_c.U); cp = c.glpoints(subdivision); @@ -837,7 +838,7 @@ namespace stim { stim::vec3 U = bb.B; // get the top right corner width = U[2] - L[2] + 10.0f; - glColor3f(1.0f, 1.0f, 1.0f); + glColor3f(0.0f, 0.0f, 0.0f); for (unsigned i = 0; i < main_feeder.size(); i++) { // front face glBegin(GL_QUADS); @@ -896,6 +897,7 @@ namespace stim { if (!glIsList(dlist)) { dlist = glGenLists(1); glNewList(dlist, GL_COMPILE); + glColor3f(0.0f, 0.0f, 0.0f); for (unsigned i = 0; i < inlet.size(); i++) { glBegin(GL_LINE_STRIP); for (unsigned j = 0; j < inlet[i].V.size(); j++) @@ -915,7 +917,7 @@ namespace stim { } // draw the bridge as tubes - void tube_bridge(T subdivision, T radii = 5.0f) { + void tube_bridge(T subdivision, T radius = 5.0f) { if (!glIsList(dlist + 1)) { glNewList(dlist + 1, GL_COMPILE); @@ -924,13 +926,14 @@ namespace stim { stim::circle unit_c; // unit circle for finding the rotation start direction std::vector > cp1; std::vector > cp2; + glColor3f(0.0f, 0.0f, 0.0f); for (unsigned i = 0; i < inlet.size(); i++) { // render vertex as sphere for (unsigned j = 1; j < inlet[i].V.size() - 1; j++) { glPushMatrix(); glTranslatef(inlet[i].V[j][0], inlet[i].V[j][1], inlet[i].V[j][2]); - glutSolidSphere(radii, subdivision, subdivision); + glutSolidSphere(radius, subdivision, subdivision); glPopMatrix(); } // render edge as cylinder @@ -957,7 +960,7 @@ namespace stim { for (unsigned j = 1; j < outlet[i].V.size() - 1; j++) { glPushMatrix(); glTranslatef(outlet[i].V[j][0], outlet[i].V[j][1], outlet[i].V[j][2]); - glutSolidSphere(radii, subdivision, subdivision); + glutSolidSphere(radius, subdivision, subdivision); glPopMatrix(); } // render edge as cylinder @@ -1038,7 +1041,7 @@ namespace stim { // mark the vertex void mark_vertex() { - glColor3f(1.0f, 1.0f, 1.0f); + glColor3f(0.0f, 0.0f, 0.0f); for (unsigned i = 0; i < num_vertex; i++) { glRasterPos3f(V[i][0], V[i][1] + get_radius(i), V[i][2]); std::stringstream ss; @@ -1137,7 +1140,7 @@ namespace stim { // build connection between all inlets and outlets // connection will trail along one axis around the bounding box - void build_synthetic_connection(T viscosity, T radii = 5.0f) { + void build_synthetic_connection(T viscosity, T radius = 5.0f) { stim::vec3 L = bb.A; // get the bottom left corner stim::vec3 U = bb.B; // get the top right corner @@ -1162,7 +1165,7 @@ namespace stim { tmp_b.v.push_back(input[i].first); tmp_b.Q = input[i].third; tmp_b.l = (bus_v - mid_v).len() + (mid_v - tmp_v).len(); - tmp_b.r = radii; + tmp_b.r = radius; inlet.push_back(tmp_b); } @@ -1182,19 +1185,206 @@ namespace stim { tmp_b.v.push_back(output[i].first); tmp_b.Q = output[i].third; tmp_b.l = (bus_v - mid_v).len() + (mid_v - tmp_v).len(); - tmp_b.r = radii; + tmp_b.r = radius; outlet.push_back(tmp_b); } } + // find the number of U-shape or square-shape structure for extending length of connection + // @param t: width = t * radius + int find_number_square(T origin_l, T desire_l, T radius = 5.0f, int times = 4) { + + bool done = false; // flag indicates the current number of square shape structure is feasible + int n = origin_l / (times * 2 * radius); // number of square shape structure + T need_l = desire_l - origin_l; + T height; // height of the square shapce structure + + while (!done) { + height = need_l / (2 * n); // calculate the height + if (height > 2 * radius) { + done = true; + } + else { + n--; + } + } + + return n; + } + + // build square connections + void build_square_connection(int i, T width, T height, T origin_l, T desire_l, int n, int feeder, T threshold, bool z, bool left = true, bool up = true, int times = 4, T ratio = 0, T radius = 5.0f) { + + int coef_up = (up) ? 1 : -1; // y coefficient + int coef_left = (left) ? 1 : -1; // x coefficient + int coef_z = (z) ? 1 : -1; // z coefficient + int inverse = 1; // inverse flag + stim::vec3 cor_v; // corner vertex + + stim::vec3 tmp_v; + if (feeder == 1) + tmp_v = inlet[i].V[0]; + else if (feeder == 0) + tmp_v = outlet[i].V[0]; + + // check whether the height of connections is acceptable + if (height > threshold) { // acceptable + // re-calculate height + if (ratio > 0.0f && ratio <= 1.0f) { // use fragment if set + cor_v = tmp_v + stim::vec3(-coef_left * origin_l, 0, 0); // get the original corner vertex + desire_l = desire_l - origin_l * (1.0f - ratio / 1.0f); + origin_l = (T)origin_l * ratio / 1.0f; + n = find_number_square(origin_l, desire_l); + } + height = (desire_l - (1 + 2 * n) * origin_l) / std::pow(2 * n, 2); + // check whether the connections are good + while (height > threshold) { + n++; + width = (T)(origin_l) / (2 * n); + height = (desire_l - (1 + 2 * n) * origin_l) / std::pow(2 * n, 2); + // check whether it appears overlap + if (width < times * radius) { + n--; + width = (T)(origin_l) / (2 * n); + height = (desire_l - (1 + 2 * n) * origin_l) / std::pow(2 * n, 2); + break; + } + } + while (height < times * radius) { + n--; + width = (T)(origin_l) / (2 * n); + height = (desire_l - (1 + 2 * n) * origin_l) / std::pow(2 * n, 2); + } + + // cube-like structure construction + for (int j = 0; j < n; j++) { + // "up" + for (int k = 0; k < n; k++) { + // in + tmp_v = tmp_v + stim::vec3(0, 0, coef_z * height); + if (feeder == 1) + inlet[i].V.push_back(tmp_v); + else if (feeder == 0) + outlet[i].V.push_back(tmp_v); + // "up" + tmp_v = tmp_v + stim::vec3(0, inverse * coef_up * width, 0); + if (feeder == 1) + inlet[i].V.push_back(tmp_v); + else if (feeder == 0) + outlet[i].V.push_back(tmp_v); + // out + tmp_v = tmp_v + stim::vec3(0, 0, -coef_z * height); + if (feeder == 1) + inlet[i].V.push_back(tmp_v); + else if (feeder == 0) + outlet[i].V.push_back(tmp_v); + // "down" + tmp_v = tmp_v + stim::vec3(0, inverse * coef_up * width, 0); + if (feeder == 1) + inlet[i].V.push_back(tmp_v); + else if (feeder == 0) + outlet[i].V.push_back(tmp_v); + } + + // "left" + tmp_v = tmp_v + stim::vec3(-coef_left * width, 0, 0); + if (feeder == 1) + inlet[i].V.push_back(tmp_v); + else if (feeder == 0) + outlet[i].V.push_back(tmp_v); + + if (inverse == 1) // revert inverse + inverse = -1; + else + inverse = 1; + + // "down" + for (int k = 0; k < n; k++) { + + tmp_v = tmp_v + stim::vec3(0, 0, coef_z * height); + if (feeder == 1) + inlet[i].V.push_back(tmp_v); + else if (feeder == 0) + outlet[i].V.push_back(tmp_v); + + tmp_v = tmp_v + stim::vec3(0, inverse * coef_up * width, 0); + if (feeder == 1) + inlet[i].V.push_back(tmp_v); + else if (feeder == 0) + outlet[i].V.push_back(tmp_v); + + tmp_v = tmp_v + stim::vec3(0, 0, -coef_z * height); + if (feeder == 1) + inlet[i].V.push_back(tmp_v); + else if (feeder == 0) + outlet[i].V.push_back(tmp_v); + + tmp_v = tmp_v + stim::vec3(0, inverse * coef_up * width, 0); + if (feeder == 1) + inlet[i].V.push_back(tmp_v); + else if (feeder == 0) + outlet[i].V.push_back(tmp_v); + } + + // "left" + tmp_v = tmp_v + stim::vec3(-coef_left * width, 0, 0); + if (feeder == 1) + inlet[i].V.push_back(tmp_v); + else if (feeder == 0) + outlet[i].V.push_back(tmp_v); + + if (inverse == 1) // revert inverse + inverse = -1; + else + inverse = 1; + } + // if use fragment to do square wave connection, need to push_back the corner vertex + if (ratio > 0.0f && ratio <= 1.0f) { + if (feeder == 1) + inlet[i].V.push_back(cor_v); + else if (feeder == 0) + outlet[i].V.push_back(cor_v); + } + } + else { + for (int j = 0; j < n; j++) { + + // move in Z-shape + tmp_v = tmp_v + stim::vec3(0, coef_up * height, 0); + if (feeder == 1) + inlet[i].V.push_back(tmp_v); + else if (feeder == 0) + outlet[i].V.push_back(tmp_v); + + tmp_v = tmp_v + stim::vec3(-coef_left * width, 0, 0); + if (feeder == 1) + inlet[i].V.push_back(tmp_v); + else if (feeder == 0) + outlet[i].V.push_back(tmp_v); + + tmp_v = tmp_v + stim::vec3(0, -coef_up * height, 0); + if (feeder == 1) + inlet[i].V.push_back(tmp_v); + else if (feeder == 0) + outlet[i].V.push_back(tmp_v); + + tmp_v = tmp_v + stim::vec3(-coef_left * width, 0, 0); + if (feeder == 1) + inlet[i].V.push_back(tmp_v); + else if (feeder == 0) + outlet[i].V.push_back(tmp_v); + } + } + } + // automatically modify bridge to make it feasible - void modify_synthetic_connection(T viscosity, T rou, T radii = 5.0f) { + void modify_synthetic_connection(T viscosity, T rou, bool H, T threshold, T ratio = 0.0f, T radius = 5.0f) { glDeleteLists(dlist, 1); // delete display list for modify glDeleteLists(dlist + 1, 1); - // because of radii change at the port vertex, there will be a pressure drop at that port + // because of radius change at the port vertex, there will be a pressure drop at that port // it follows the bernoulli equation // p1 + 1/2*rou*v1^2 + rou*g*h1 = p2 + 1/2*rou*v2^2 + rou*g*h2 // Q1 = Q2 -> v1*r1^2 = v2*r2^2 @@ -1203,7 +1393,7 @@ namespace stim { for (unsigned i = 0; i < pendant_vertex.size(); i++) { idx = pendant_vertex[i]; T tmp_v = get_velocity(idx); // velocity at that pendant vertex - T ar = get_radius(idx) / radii; + T ar = get_radius(idx) / radius; new_pressure[idx] = pressure[idx] + 1.0f / 2.0f * rou * std::pow(tmp_v, 2) * (1.0f - std::pow(ar, 4)); } @@ -1213,217 +1403,302 @@ namespace stim { unsigned inlet_index; T tmp_p; for (unsigned i = 0; i < inlet.size(); i++) { - tmp_p = new_pressure[inlet[i].v[0]] + ((8 * viscosity * inlet[i].l * inlet[i].Q) / ((float)stim::PI * std::pow(radii, 4))); + tmp_p = new_pressure[inlet[i].v[0]] + ((8 * viscosity * inlet[i].l * inlet[i].Q) / ((float)stim::PI * std::pow(radius, 4))); if (tmp_p > source_pressure) { source_pressure = tmp_p; inlet_index = i; } } - // automatically modify inlet bridge to make it feasible - bool upper = false; // flag indicates the whether the port is upper than main feeder - bool invert = false; // there are two version of hilbert curve depends on starting position with respect to the cup - T new_l; - stim::vec3 bus_v; // the port point on the bus - stim::vec3 mid_v; // the original corner point - stim::vec3 tmp_v; // the pendant point - int order = 0; // order of hilbert curve (iteration) - for (unsigned i = 0; i < inlet.size(); i++) { - if (i != inlet_index) { - new_l = (source_pressure - new_pressure[inlet[i].v[0]]) * ((float)stim::PI * std::pow(radii, 4)) / (8 * viscosity * inlet[i].Q); + // find minimum pressure outlet port + T end_pressure = FLT_MAX; + unsigned outlet_index; + for (unsigned i = 0; i < outlet.size(); i++) { + tmp_p = new_pressure[outlet[i].v[0]] - ((8 * viscosity * outlet[i].l * outlet[i].Q) / ((float)stim::PI * std::pow(radius, 4))); + if (tmp_p < end_pressure) { + end_pressure = tmp_p; + outlet_index = i; + } + } - if (inlet[i].V[2][1] > main_feeder[0][1]) { // check out upper side of lower side - upper = true; - invert = false; - } - else { - upper = false; - invert = true; - } + // automatically modify inlet bridge using Hilbert curves + if (H) { + bool upper = false; // flag indicates the whether the port is upper than main feeder + bool invert = false; // there are two version of hilbert curve depends on starting position with respect to the cup + T new_l; + stim::vec3 bus_v; // the port point on the bus + stim::vec3 mid_v; // the original corner point + stim::vec3 tmp_v; // the pendant point + int order = 0; // order of hilbert curve (iteration) + for (unsigned i = 0; i < inlet.size(); i++) { + if (i != inlet_index) { + new_l = (source_pressure - new_pressure[inlet[i].v[0]]) * ((float)stim::PI * std::pow(radius, 4)) / (8 * viscosity * inlet[i].Q); - T origin_l = (inlet[i].V[1] - inlet[i].V[2]).len(); - T desire_l = new_l - (inlet[i].V[0] - inlet[i].V[1]).len(); - find_hilbert_order(origin_l, desire_l, order); - - bus_v = inlet[i].V[0]; - mid_v = inlet[i].V[1]; - tmp_v = inlet[i].V[2]; - inlet[i].V.clear(); - inlet[i].V.push_back(tmp_v); - inlet[i].l = new_l; - - 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); + if (inlet[i].V[2][1] > main_feeder[0][1]) { // check out upper side of lower side + upper = true; + invert = false; + } + else { + upper = false; + invert = true; + } + + T origin_l = (inlet[i].V[1] - inlet[i].V[2]).len(); + T desire_l = new_l - (inlet[i].V[0] - inlet[i].V[1]).len(); + find_hilbert_order(origin_l, desire_l, order); + + bus_v = inlet[i].V[0]; + mid_v = inlet[i].V[1]; + tmp_v = inlet[i].V[2]; + inlet[i].V.clear(); + inlet[i].V.push_back(tmp_v); + inlet[i].l = new_l; + + if (desire_l - origin_l < 2 * radius) { // 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 { // if the radii is not feasible - int count = 1; - while (dl <= 2 * radii) { - dl = origin_l / (std::pow(2, order - count) - 1); - count++; - } - count--; + 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 (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); + if (dl > 2 * radius) { // if the radius 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); - 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; + if (tmp_v[0] != mid_v[0]) + inlet[i].V.push_back(mid_v); + inlet[i].V.push_back(bus_v); + } + else { // if the radius is not feasible + int count = 1; + while (dl <= 2 * radius) { + 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 * radius) + std::cout << "infeasible connection between inlets!" << std::endl; + + if (upper) + hilbert_curve(i, &tmp_v[0], order, dl, 1, !invert, LEFT); + 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 + } + } - find_hilbert_order(origin_l, desire_l, order); + // automatically modify outlet bridge to make it feasible + for (unsigned i = 0; i < outlet.size(); i++) { + if (i != outlet_index) { + new_l = (new_pressure[outlet[i].v[0]] - end_pressure) * ((float)stim::PI * std::pow(radius, 4)) / (8 * viscosity * outlet[i].Q); - 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 (outlet[i].V[2][1] > main_feeder[1][1]) { + upper = true; + invert = true; + } + else { + upper = false; + invert = false; + } - if (upper) - hilbert_curve(i, &tmp_v[0], order, dl, 1, !invert, LEFT); - else - hilbert_curve(i, &tmp_v[0], order, dl, 1, !invert, RIGHT); + T origin_l = (outlet[i].V[1] - outlet[i].V[2]).len(); + T desire_l = new_l - (outlet[i].V[0] - outlet[i].V[1]).len(); + find_hilbert_order(origin_l, desire_l, order); + + bus_v = outlet[i].V[0]; + mid_v = outlet[i].V[1]; + tmp_v = outlet[i].V[2]; + outlet[i].V.clear(); + outlet[i].V.push_back(tmp_v); + outlet[i].l = new_l; + + if (desire_l - origin_l < 2 * radius) { // 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 (tmp_v[1] != bus_v[1]) - inlet[i].V.push_back(bus_v); + if (dl > 2 * radius) { // if the radius 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 radius is not feasible + int count = 1; + while (dl <= 2 * radius) { + 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 * radius) + 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, RIGHT); + + if (tmp_v[1] != bus_v[1]) + outlet[i].V.push_back(bus_v); + } } + std::reverse(outlet[i].V.begin(), outlet[i].V.end()); } - std::reverse(inlet[i].V.begin(), inlet[i].V.end()); // from bus to pendant vertex } } + // automatically modify inlet bridge using square shape constructions + else { + bool upper; // flag indicates the connection is upper than the bus + bool z; // flag indicates the connection direction along z-axis + T new_l; // new length + stim::vec3 bus_v; // the port point on the bus + stim::vec3 mid_v; // the original corner point + stim::vec3 tmp_v; // the pendant point + int n; + T width, height; // width and height of the square - // find minimum pressure outlet port - source_pressure = FLT_MAX; - unsigned outlet_index; - for (unsigned i = 0; i < outlet.size(); i++) { - tmp_p = new_pressure[outlet[i].v[0]] - ((8 * viscosity * outlet[i].l * outlet[i].Q) / ((float)stim::PI * std::pow(radii, 4))); - if (tmp_p < source_pressure) { - source_pressure = tmp_p; - outlet_index = i; - } - } + for (unsigned i = 0; i < inlet.size(); i++) { + if (i != inlet_index) { + new_l = (source_pressure - new_pressure[inlet[i].v[0]]) * ((float)stim::PI * std::pow(radius, 4)) / (8 * viscosity * inlet[i].Q); // calculate the new length of the connection - // automatically modify outlet bridge to make it feasible - 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); + bus_v = inlet[i].V[0]; + mid_v = inlet[i].V[1]; + tmp_v = inlet[i].V[2]; - if (outlet[i].V[2][1] > main_feeder[1][1]) { - upper = true; - invert = true; - } - else { - upper = false; - invert = false; - } + if (inlet[i].V[2][1] > main_feeder[0][1]) // check out upper side of lower side + upper = true; + else + upper = false; - T origin_l = (outlet[i].V[1] - outlet[i].V[2]).len(); - T desire_l = new_l - (outlet[i].V[0] - outlet[i].V[1]).len(); - find_hilbert_order(origin_l, desire_l, order); - - bus_v = outlet[i].V[0]; - mid_v = outlet[i].V[1]; - tmp_v = outlet[i].V[2]; - outlet[i].V.clear(); - outlet[i].V.push_back(tmp_v); - outlet[i].l = new_l; - - 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); + if (inlet[i].V[2][2] > main_feeder[0][2]) + z = true; + else + z = false; + + T origin_l = (inlet[i].V[1] - inlet[i].V[2]).len(); + T desire_l = new_l - (inlet[i].V[0] - inlet[i].V[1]).len(); + inlet[i].V.clear(); + inlet[i].V.push_back(tmp_v); + inlet[i].l = new_l; + + n = find_number_square(origin_l, desire_l); + + width = (T)origin_l / (2 * n); + height = (desire_l - origin_l) / (2 * n); + + build_square_connection(i, width, height, origin_l, desire_l, n, 1, threshold, z, true, upper, 10, ratio); + inlet[i].V.push_back(bus_v); + + std::reverse(inlet[i].V.begin(), inlet[i].V.end()); // from bus to pendant vertex } - 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--; + } + + for (unsigned i = 0; i < outlet.size(); i++) { + if (i != outlet_index) { + new_l = (new_pressure[outlet[i].v[0]] - end_pressure) * ((float)stim::PI * std::pow(radius, 4)) / (8 * viscosity * outlet[i].Q); // calculate the new length of the connection - 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); + bus_v = outlet[i].V[0]; + mid_v = outlet[i].V[1]; + tmp_v = outlet[i].V[2]; - 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; + if (outlet[i].V[2][1] > main_feeder[1][1]) // check out upper side of lower side + upper = true; + else + upper = false; - find_hilbert_order(origin_l, desire_l, order); + if (outlet[i].V[2][2] > main_feeder[1][2]) + z = true; + else + z = false; - 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; + T origin_l = (outlet[i].V[1] - outlet[i].V[2]).len(); + T desire_l = new_l - (outlet[i].V[0] - outlet[i].V[1]).len(); + outlet[i].V.clear(); + outlet[i].V.push_back(tmp_v); + outlet[i].l = new_l; - if (upper) - hilbert_curve(i, &tmp_v[0], order, dl, 0, !invert, LEFT); - else - hilbert_curve(i, &tmp_v[0], order, dl, 0, !invert, RIGHT); + n = find_number_square(origin_l, desire_l); - if (tmp_v[1] != bus_v[1]) - outlet[i].V.push_back(bus_v); - } + width = (T)origin_l / (2 * n); + height = (desire_l - origin_l) / (2 * n); + + build_square_connection(i, width, height, origin_l, desire_l, n, 0, threshold, z, false, upper, 10, ratio); + outlet[i].V.push_back(bus_v); + + std::reverse(outlet[i].V.begin(), outlet[i].V.end()); // from bus to pendant vertex } - std::reverse(outlet[i].V.begin(), outlet[i].V.end()); } } } // check current bridge to see feasibility - void check_synthetic_connection(T viscosity, T radii = 5.0f) { + void check_synthetic_connection(T viscosity, T radius = 5.0f) { T eps = 0.01f; - T source_pressure = pressure[inlet[0].v[0]] + (8 * viscosity * inlet[0].l * inlet[0].Q) / ((float)stim::PI * std::pow(radii, 4)); + T source_pressure = pressure[inlet[0].v[0]] + (8 * viscosity * inlet[0].l * inlet[0].Q) / ((float)stim::PI * std::pow(radius, 4)); T tmp_p; for (unsigned i = 1; i < inlet.size(); i++) { - tmp_p = pressure[inlet[i].v[0]] + (8 * viscosity * inlet[i].l * inlet[i].Q) / ((float)stim::PI * std::pow(radii, 4)); + tmp_p = pressure[inlet[i].v[0]] + (8 * viscosity * inlet[i].l * inlet[i].Q) / ((float)stim::PI * std::pow(radius, 4)); T delta = fabs(tmp_p - source_pressure); if (delta > eps) { std::cout << "Nonfeasible connection!" << std::endl; break; } } - source_pressure = pressure[outlet[0].v[0]] - (8 * viscosity * outlet[0].l * outlet[0].Q) / ((float)stim::PI * std::pow(radii, 4)); + source_pressure = pressure[outlet[0].v[0]] - (8 * viscosity * outlet[0].l * outlet[0].Q) / ((float)stim::PI * std::pow(radius, 4)); for (unsigned i = 1; i < outlet.size(); i++) { - tmp_p = pressure[outlet[i].v[0]] - (8 * viscosity * outlet[i].l * outlet[i].Q) / ((float)stim::PI * std::pow(radii, 4)); + tmp_p = pressure[outlet[i].v[0]] - (8 * viscosity * outlet[i].l * outlet[i].Q) / ((float)stim::PI * std::pow(radius, 4)); T delta = fabs(tmp_p - source_pressure); if (delta > eps) { std::cout << "Nonfeasible connection!" << std::endl; @@ -1436,7 +1711,7 @@ namespace stim { // prepare for image stack void preparation(T &Xl, T &Xr, T &Yt, T &Yb, T &Z, T length = 210.0f, T height = 10.0f) { - T max_radii = 0.0f; + T max_radius = 0.0f; T top = FLT_MIN; T bottom = FLT_MAX; @@ -1519,17 +1794,17 @@ namespace stim { top = A[i].c[1]; if (A[i].c[1] < bottom) bottom = A[i].c[1]; - if (A[i].r > max_radii) - max_radii = A[i].r; + if (A[i].r > max_radius) + max_radius = A[i].r; } Yt = top; // top bound y coordinate Yb = bottom; // bottom bound y coordinate - Z = (bb.B[2] - bb.A[2] + 2 * max_radii); // bounding box width(along z-axis) + Z = (bb.B[2] - bb.A[2] + 2 * max_radius); // bounding box width(along z-axis) } /// making image stack main function - void make_image_stack(T dx, T dy, T dz, std::string stackdir, T radii = 5.0f) { + void make_image_stack(T dx, T dy, T dz, std::string stackdir, T radius = 5.0f) { /// preparation for making image stack T X, Xl, Xr, Y, Yt, Yb, Z; diff --git a/main.cu b/main.cu index 5db4cda..a122423 100644 --- a/main.cu +++ b/main.cu @@ -54,10 +54,12 @@ float zoom_factor = 10.0f; // zooming factor float border_factor = 20.0f; // border float radii_factor = 1.0f; // radii changing factor GLint subdivision = 20; // slices and stacks -float default_radii = 5.0f; // default radii of network vertex +float default_radius = 5.0f; // default radii of network vertex float delta = 0.01f; // small discrepancy float eps = 20.0f; // epsilon threshold float max_pressure = 0.0f; // maximum pressure that the channel can bear +float height_threshold = 100.0f; // connection height constraint +float fragment_ratio = 0.0f; // fragment ratio // glut event parameters int mouse_x; // window x-coordinate @@ -73,6 +75,8 @@ unsigned pressure_index; // the index of vertex that is clicked // build inlet/outlet parameters bool build_inlet_outlet = false; // flag indicates building inlets and outlets bool modified_bridge = false; // flag indicates having modified inlet/outlet connection +bool hilbert_curve = false; // flag indicates enabling hilbert curves constructions +bool change_fragment = false; // flag indicates changing fragment for square wave connections // manufacture parameters bool manufacture = false; // flag indicates manufacture mode @@ -89,7 +93,7 @@ inline void get_background() { // set the initial radii flow.init(num_edge, num_vertex); // initialize flow object for (unsigned i = 0; i < num_edge; i++) - flow.set_r(i, default_radii); + flow.set_r(i, default_radius); } // convert from window coordinates to world coordinates @@ -116,6 +120,7 @@ inline void window_to_world(GLdouble &x, GLdouble &y, GLdouble &z) { // initialize flow object void flow_initialize() { + flow.set = true; stim::vec3 center = bb.center(); for (unsigned i = 0; i < pendant_vertex.size(); i++) { @@ -187,6 +192,7 @@ void glut_render() { glEnable(GL_DEPTH_TEST); glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); + glClearColor(1.0f, 1.0f, 1.0f, 1.0f); glut_projection(); glut_modelview(); @@ -238,7 +244,7 @@ void glut_render() { glFlush(); // pressure bar text - glColor3f(1.0f, 1.0f, 1.0f); + glColor3f(0.0f, 0.0f, 0.0f); glRasterPos2f(0.0f, vY - border_factor); std::stringstream ss_p; ss_p << "Pressure Bar"; @@ -251,7 +257,7 @@ void glut_render() { glRasterPos2f((border_factor * 1.5f), (border_factor + i * step)); std::stringstream ss_n; ss_n << (float)i * max_pressure / 10; - glutBitmapString(GLUT_BITMAP_TIMES_ROMAN_10, (const unsigned char*)(ss_n.str().c_str())); + glutBitmapString(GLUT_BITMAP_HELVETICA_18, (const unsigned char*)(ss_n.str().c_str())); } glPopMatrix(); glMatrixMode(GL_PROJECTION); @@ -286,7 +292,7 @@ void glut_render() { glFlush(); // pressure bar text - glColor3f(1.0f, 1.0f, 1.0f); + glColor3f(0.0f, 0.0f, 0.0f); glRasterPos2f(0.0f, vY - border_factor); std::stringstream ss_p; ss_p << "Velocity range"; @@ -299,13 +305,58 @@ void glut_render() { glRasterPos2f(border_factor * 1.5f, border_factor + i * step); std::stringstream ss_n; ss_n << min_v + i * (max_v - min_v) / 10; - glutBitmapString(GLUT_BITMAP_TIMES_ROMAN_10, (const unsigned char*)(ss_n.str().c_str())); + glutBitmapString(GLUT_BITMAP_HELVETICA_18, (const unsigned char*)(ss_n.str().c_str())); } glPopMatrix(); glMatrixMode(GL_PROJECTION); glPopMatrix(); } - + + // bring up a ratio bar on the left + if (change_fragment) { + + glMatrixMode(GL_PROJECTION); // set up the 2d viewport for mode text printing + glPushMatrix(); + glLoadIdentity(); + vX = glutGet(GLUT_WINDOW_WIDTH); // get the current window width + vY = glutGet(GLUT_WINDOW_HEIGHT); // get the current window height + glViewport(0, 0, vX, vY); // locate to left bottom corner + gluOrtho2D(0, vX, 0, vY); // define othogonal aspect + + glMatrixMode(GL_MODELVIEW); + glPushMatrix(); + glLoadIdentity(); + + glLineWidth(border_factor); + glBegin(GL_LINES); + glColor3f(0.0, 0.0, 1.0); // blue to red + glVertex2f(border_factor, border_factor); + glColor3f(1.0, 0.0, 0.0); + glVertex2f(border_factor, (vY - 2.0f * border_factor)); + glEnd(); + glFlush(); + + // pressure bar text + glColor3f(0.0f, 0.0f, 0.0f); + glRasterPos2f(0.0f, vY - border_factor); + std::stringstream ss_p; + ss_p << "Ratio bar"; + glutBitmapString(GLUT_BITMAP_HELVETICA_18, (const unsigned char*)(ss_p.str().c_str())); + + // pressure range text + float step = vY - 3.0f * border_factor; + step /= 10; + for (unsigned i = 0; i < 11; i++) { + glRasterPos2f((border_factor * 1.5f), (border_factor + i * step)); + std::stringstream ss_n; + ss_n << (float)i * 1.0f / 10; + glutBitmapString(GLUT_BITMAP_HELVETICA_18, (const unsigned char*)(ss_n.str().c_str())); + } + glPopMatrix(); + glMatrixMode(GL_PROJECTION); + glPopMatrix(); + } + glutSwapBuffers(); } @@ -318,7 +369,8 @@ void glut_menu(int value) { build_inlet_outlet = false; manufacture = false; modified_bridge = false; - flow_initialize(); + if (!flow.set) + flow_initialize(); flow_stable_state(); // main function of solving the linear system flow.print_flow(); @@ -331,7 +383,7 @@ void glut_menu(int value) { manufacture = false; if (!modified_bridge) { flow.set_main_feeder(); - flow.build_synthetic_connection(u); + flow.build_synthetic_connection(u, default_radius); } glut_set_menu(num, 3); @@ -341,7 +393,7 @@ void glut_menu(int value) { simulation = false; build_inlet_outlet = false; manufacture = true; - flow.check_synthetic_connection(u); + flow.check_synthetic_connection(u, default_radius); } glutPostRedisplay(); @@ -396,6 +448,13 @@ void glut_mouse(int button, int state, int x, int y) { flow.print_flow(); } } + else if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN && modified_bridge && change_fragment) { + if (y > 2 * border_factor || y < vY - border_factor) { // within the ratio bar range + fragment_ratio = (float)(vY - y - border_factor) / ((float)vY - border_factor) * 1.0f; + flow.modify_synthetic_connection(u, rou, hilbert_curve, height_threshold, fragment_ratio, default_radius); + change_fragment = false; + } + } } // define camera move based on mouse wheel move @@ -419,7 +478,7 @@ void glut_wheel(int wheel, int direction, int x, int y) { tmp_r = flow.get_radius(pressure_index); tmp_r -= radii_factor; if (tmp_r <= 0) - tmp_r = default_radii; + tmp_r = default_radius; } flow.set_radius(pressure_index, tmp_r); flow_stable_state(); @@ -464,12 +523,16 @@ void glut_keyboard(unsigned char key, int x, int y) { #endif } else if (build_inlet_outlet && !modified_bridge) { - flow.modify_synthetic_connection(u, rou); modified_bridge = true; + + if (hilbert_curve) + flow.modify_synthetic_connection(u, rou, hilbert_curve, height_threshold); + else + change_fragment = true; } break; } - + glutPostRedisplay(); } @@ -510,6 +573,7 @@ int main(int argc, char* argv[]) { args.add("maxpress", "maximum allowed pressure in g / units / s^2, default 2 is for blood when units = um", "2", "real value > 0"); args.add("viscosity", "set the viscosity of the fluid (in g / units / s), default .00001 is for blood when units = um", ".00001", "real value > 0"); args.add("rou", "set the desity of the fluid (in g / units^3), default 1.06*10^-12 is for blood when units = um", ".00000000000106", "real value > 0"); + args.add("hilbert", "enable hilbert curves connections", "0", "value 1 for enablement"); args.add("stackres", "spacing between pixel samples in each dimension(in units/pixel)", "1 1 1", "real value > 0"); args.add("stackdir", "set the directory of the output image stack", "", "any existing directory (ex. /home/name/network)"); @@ -540,6 +604,9 @@ int main(int argc, char* argv[]) { // normally the blood density in capillaries: 1060 kg/m^3 = 1.06*10^-12 g/um^3 rou = args["rou"].as_float(); + // check whether to enable hilbert curves or not + hilbert_curve = args["hilbert"].as_int(); + // get the vexel and image stack size dx = args["stackres"].as_float(0); dy = args["stackres"].as_float(1); -- libgit2 0.21.4