Jiaming Guo / flow3

Commit 8334680aedb2270ba2dbd0924240d0e9fff3ad4b

Authored by Jiaming Guo
1 parent 9b69bb4c

add square wave connections

Showing 2 changed files with 546 additions and 204 deletions   Show diff stats
Inline Side-by-side
@@ -30,7 +30,7 @@ namespace stim { @@ -30,7 +30,7 @@ namespace stim {
30 std::vector<unsigned> v; // vertices' indices 30 std::vector<unsigned> v; // vertices' indices
31 std::vector<typename stim::vec3<T> > V; // vertices' coordinates 31 std::vector<typename stim::vec3<T> > V; // vertices' coordinates
32 T l; // length 32 T l; // length
33 - T r; // radii 33 + T r; // radius
34 T deltaP; // pressure drop 34 T deltaP; // pressure drop
35 T Q; // volume flow rate 35 T Q; // volume flow rate
36 }; 36 };
@@ -38,15 +38,15 @@ namespace stim { @@ -38,15 +38,15 @@ namespace stim {
38 template <typename T> 38 template <typename T>
39 struct sphere { 39 struct sphere {
40 stim::vec3<T> c; // center of sphere 40 stim::vec3<T> c; // center of sphere
41 - T r; // radii 41 + T r; // radius
42 }; 42 };
43 43
44 template <typename T> 44 template <typename T>
45 - struct cone { // radii changes gradually 45 + struct cone { // radius changes gradually
46 stim::vec3<T> c1; // center of geometry start hat 46 stim::vec3<T> c1; // center of geometry start hat
47 stim::vec3<T> c2; // center of geometry end hat 47 stim::vec3<T> c2; // center of geometry end hat
48 - T r1; // radii at start hat  
49 - T r2; // radii at end hat 48 + T r1; // radius at start hat
  49 + T r2; // radius at end hat
50 }; 50 };
51 51
52 template <typename T> 52 template <typename T>
@@ -105,7 +105,7 @@ namespace stim { @@ -105,7 +105,7 @@ namespace stim {
105 105
106 float distance = FLT_MAX; 106 float distance = FLT_MAX;
107 float tmp_distance; 107 float tmp_distance;
108 - float rr; // radii at the surface where projection meets 108 + float rr; // radius at the surface where projection meets
109 109
110 for (unsigned i = 0; i < num; i++) { // find the nearest cylinder 110 for (unsigned i = 0; i < num; i++) { // find the nearest cylinder
111 tmp_distance = ((world_pixel - E[i].c1).cross(world_pixel - E[i].c2)).len() / (E[i].c2 - E[i].c1).len(); 111 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 { @@ -241,6 +241,7 @@ namespace stim {
241 241
242 public: 242 public:
243 243
  244 + bool set = false; // flag indicates the pressure has been set
244 std::vector<T> P; // initial pressure 245 std::vector<T> P; // initial pressure
245 std::vector<T> v; // velocity 246 std::vector<T> v; // velocity
246 std::vector<typename stim::vec3<T> > main_feeder; // inlet/outlet main feeder 247 std::vector<typename stim::vec3<T> > main_feeder; // inlet/outlet main feeder
@@ -311,7 +312,7 @@ namespace stim { @@ -311,7 +312,7 @@ namespace stim {
311 } 312 }
312 } 313 }
313 314
314 - // get the radii of vertex i 315 + // get the radius of vertex i
315 T get_radius(unsigned i) { 316 T get_radius(unsigned i) {
316 317
317 unsigned tmp_e; // edge index 318 unsigned tmp_e; // edge index
@@ -596,7 +597,7 @@ namespace stim { @@ -596,7 +597,7 @@ namespace stim {
596 // find two envelope caps for two spheres 597 // find two envelope caps for two spheres
597 // @param cp1, cp2: list of points on the cap 598 // @param cp1, cp2: list of points on the cap
598 // @param center1, center2: center point of cap 599 // @param center1, center2: center point of cap
599 - // @param r1, r2: radii of cap 600 + // @param r1, r2: radius of cap
600 void find_envelope(std::vector<typename stim::vec3<float> > &cp1, std::vector<typename stim::vec3<float> > &cp2, stim::vec3<float> center1, stim::vec3<float> center2, float r1, float r2, GLint subdivision) { 601 void find_envelope(std::vector<typename stim::vec3<float> > &cp1, std::vector<typename stim::vec3<float> > &cp2, stim::vec3<float> center1, stim::vec3<float> center2, float r1, float r2, GLint subdivision) {
601 602
602 stim::vec3<float> tmp_d; 603 stim::vec3<float> tmp_d;
@@ -700,7 +701,7 @@ namespace stim { @@ -700,7 +701,7 @@ namespace stim {
700 // calculate the bigger circle 701 // calculate the bigger circle
701 d1 = t1 - center1; 702 d1 = t1 - center1;
702 dot = d1.dot(tmp_d); 703 dot = d1.dot(tmp_d);
703 - a = dot / (r1 * 1) * r1; // a = cos(alpha) * radii 704 + a = dot / (r1 * 1) * r1; // a = cos(alpha) * radius
704 new_c = center1 + a * tmp_d; 705 new_c = center1 + a * tmp_d;
705 new_r = sqrt(pow(r1, 2) - pow(a, 2)); 706 new_r = sqrt(pow(r1, 2) - pow(a, 2));
706 new_u = t1 - new_c; 707 new_u = t1 - new_c;
@@ -803,7 +804,7 @@ namespace stim { @@ -803,7 +804,7 @@ namespace stim {
803 std::vector<typename stim::vec3<T> > cp; 804 std::vector<typename stim::vec3<T> > cp;
804 T radius; 805 T radius;
805 806
806 - glColor3f(1.0f, 1.0f, 1.0f); 807 + glColor3f(0.0f, 0.0f, 0.0f);
807 for (unsigned i = 0; i < num_edge; i++) { // for every edge 808 for (unsigned i = 0; i < num_edge; i++) { // for every edge
808 for (unsigned j = 0; j < E[i].size() - 1; j++) { // for every point on current edge 809 for (unsigned j = 0; j < E[i].size() - 1; j++) { // for every point on current edge
809 tmp_d = E[i][j + 1] - E[i][j]; 810 tmp_d = E[i][j + 1] - E[i][j];
@@ -812,9 +813,9 @@ namespace stim { @@ -812,9 +813,9 @@ namespace stim {
812 tmp_c.rotate(tmp_d); 813 tmp_c.rotate(tmp_d);
813 radius = (E[i].r(j + 1) + E[i].r(j)) / 2; 814 radius = (E[i].r(j + 1) + E[i].r(j)) / 2;
814 if (v[i] > 0) // if flow flows from j to j+1 815 if (v[i] > 0) // if flow flows from j to j+1
815 - head = center + tmp_d * sqrt(3) * radius; 816 + head = center + tmp_d * 2 * sqrt(3) * radius;
816 else 817 else
817 - head = center - tmp_d * sqrt(3) * radius; 818 + head = center - tmp_d * 2 * sqrt(3) * radius;
818 819
819 stim::circle<float> c(center, radius, tmp_d, tmp_c.U); 820 stim::circle<float> c(center, radius, tmp_d, tmp_c.U);
820 cp = c.glpoints(subdivision); 821 cp = c.glpoints(subdivision);
@@ -837,7 +838,7 @@ namespace stim { @@ -837,7 +838,7 @@ namespace stim {
837 stim::vec3<T> U = bb.B; // get the top right corner 838 stim::vec3<T> U = bb.B; // get the top right corner
838 width = U[2] - L[2] + 10.0f; 839 width = U[2] - L[2] + 10.0f;
839 840
840 - glColor3f(1.0f, 1.0f, 1.0f); 841 + glColor3f(0.0f, 0.0f, 0.0f);
841 for (unsigned i = 0; i < main_feeder.size(); i++) { 842 for (unsigned i = 0; i < main_feeder.size(); i++) {
842 // front face 843 // front face
843 glBegin(GL_QUADS); 844 glBegin(GL_QUADS);
@@ -896,6 +897,7 @@ namespace stim { @@ -896,6 +897,7 @@ namespace stim {
896 if (!glIsList(dlist)) { 897 if (!glIsList(dlist)) {
897 dlist = glGenLists(1); 898 dlist = glGenLists(1);
898 glNewList(dlist, GL_COMPILE); 899 glNewList(dlist, GL_COMPILE);
  900 + glColor3f(0.0f, 0.0f, 0.0f);
899 for (unsigned i = 0; i < inlet.size(); i++) { 901 for (unsigned i = 0; i < inlet.size(); i++) {
900 glBegin(GL_LINE_STRIP); 902 glBegin(GL_LINE_STRIP);
901 for (unsigned j = 0; j < inlet[i].V.size(); j++) 903 for (unsigned j = 0; j < inlet[i].V.size(); j++)
@@ -915,7 +917,7 @@ namespace stim { @@ -915,7 +917,7 @@ namespace stim {
915 } 917 }
916 918
917 // draw the bridge as tubes 919 // draw the bridge as tubes
918 - void tube_bridge(T subdivision, T radii = 5.0f) { 920 + void tube_bridge(T subdivision, T radius = 5.0f) {
919 921
920 if (!glIsList(dlist + 1)) { 922 if (!glIsList(dlist + 1)) {
921 glNewList(dlist + 1, GL_COMPILE); 923 glNewList(dlist + 1, GL_COMPILE);
@@ -924,13 +926,14 @@ namespace stim { @@ -924,13 +926,14 @@ namespace stim {
924 stim::circle<T> unit_c; // unit circle for finding the rotation start direction 926 stim::circle<T> unit_c; // unit circle for finding the rotation start direction
925 std::vector<typename stim::vec3<T> > cp1; 927 std::vector<typename stim::vec3<T> > cp1;
926 std::vector<typename stim::vec3<T> > cp2; 928 std::vector<typename stim::vec3<T> > cp2;
  929 + glColor3f(0.0f, 0.0f, 0.0f);
927 930
928 for (unsigned i = 0; i < inlet.size(); i++) { 931 for (unsigned i = 0; i < inlet.size(); i++) {
929 // render vertex as sphere 932 // render vertex as sphere
930 for (unsigned j = 1; j < inlet[i].V.size() - 1; j++) { 933 for (unsigned j = 1; j < inlet[i].V.size() - 1; j++) {
931 glPushMatrix(); 934 glPushMatrix();
932 glTranslatef(inlet[i].V[j][0], inlet[i].V[j][1], inlet[i].V[j][2]); 935 glTranslatef(inlet[i].V[j][0], inlet[i].V[j][1], inlet[i].V[j][2]);
933 - glutSolidSphere(radii, subdivision, subdivision); 936 + glutSolidSphere(radius, subdivision, subdivision);
934 glPopMatrix(); 937 glPopMatrix();
935 } 938 }
936 // render edge as cylinder 939 // render edge as cylinder
@@ -957,7 +960,7 @@ namespace stim { @@ -957,7 +960,7 @@ namespace stim {
957 for (unsigned j = 1; j < outlet[i].V.size() - 1; j++) { 960 for (unsigned j = 1; j < outlet[i].V.size() - 1; j++) {
958 glPushMatrix(); 961 glPushMatrix();
959 glTranslatef(outlet[i].V[j][0], outlet[i].V[j][1], outlet[i].V[j][2]); 962 glTranslatef(outlet[i].V[j][0], outlet[i].V[j][1], outlet[i].V[j][2]);
960 - glutSolidSphere(radii, subdivision, subdivision); 963 + glutSolidSphere(radius, subdivision, subdivision);
961 glPopMatrix(); 964 glPopMatrix();
962 } 965 }
963 // render edge as cylinder 966 // render edge as cylinder
@@ -1038,7 +1041,7 @@ namespace stim { @@ -1038,7 +1041,7 @@ namespace stim {
1038 // mark the vertex 1041 // mark the vertex
1039 void mark_vertex() { 1042 void mark_vertex() {
1040 1043
1041 - glColor3f(1.0f, 1.0f, 1.0f); 1044 + glColor3f(0.0f, 0.0f, 0.0f);
1042 for (unsigned i = 0; i < num_vertex; i++) { 1045 for (unsigned i = 0; i < num_vertex; i++) {
1043 glRasterPos3f(V[i][0], V[i][1] + get_radius(i), V[i][2]); 1046 glRasterPos3f(V[i][0], V[i][1] + get_radius(i), V[i][2]);
1044 std::stringstream ss; 1047 std::stringstream ss;
@@ -1137,7 +1140,7 @@ namespace stim { @@ -1137,7 +1140,7 @@ namespace stim {
1137 1140
1138 // build connection between all inlets and outlets 1141 // build connection between all inlets and outlets
1139 // connection will trail along one axis around the bounding box 1142 // connection will trail along one axis around the bounding box
1140 - void build_synthetic_connection(T viscosity, T radii = 5.0f) { 1143 + void build_synthetic_connection(T viscosity, T radius = 5.0f) {
1141 1144
1142 stim::vec3<T> L = bb.A; // get the bottom left corner 1145 stim::vec3<T> L = bb.A; // get the bottom left corner
1143 stim::vec3<T> U = bb.B; // get the top right corner 1146 stim::vec3<T> U = bb.B; // get the top right corner
@@ -1162,7 +1165,7 @@ namespace stim { @@ -1162,7 +1165,7 @@ namespace stim {
1162 tmp_b.v.push_back(input[i].first); 1165 tmp_b.v.push_back(input[i].first);
1163 tmp_b.Q = input[i].third; 1166 tmp_b.Q = input[i].third;
1164 tmp_b.l = (bus_v - mid_v).len() + (mid_v - tmp_v).len(); 1167 tmp_b.l = (bus_v - mid_v).len() + (mid_v - tmp_v).len();
1165 - tmp_b.r = radii; 1168 + tmp_b.r = radius;
1166 1169
1167 inlet.push_back(tmp_b); 1170 inlet.push_back(tmp_b);
1168 } 1171 }
@@ -1182,19 +1185,206 @@ namespace stim { @@ -1182,19 +1185,206 @@ namespace stim {
1182 tmp_b.v.push_back(output[i].first); 1185 tmp_b.v.push_back(output[i].first);
1183 tmp_b.Q = output[i].third; 1186 tmp_b.Q = output[i].third;
1184 tmp_b.l = (bus_v - mid_v).len() + (mid_v - tmp_v).len(); 1187 tmp_b.l = (bus_v - mid_v).len() + (mid_v - tmp_v).len();
1185 - tmp_b.r = radii; 1188 + tmp_b.r = radius;
1186 1189
1187 outlet.push_back(tmp_b); 1190 outlet.push_back(tmp_b);
1188 } 1191 }
1189 } 1192 }
1190 1193
  1194 + // find the number of U-shape or square-shape structure for extending length of connection
  1195 + // @param t: width = t * radius
  1196 + int find_number_square(T origin_l, T desire_l, T radius = 5.0f, int times = 4) {
  1197 +
  1198 + bool done = false; // flag indicates the current number of square shape structure is feasible
  1199 + int n = origin_l / (times * 2 * radius); // number of square shape structure
  1200 + T need_l = desire_l - origin_l;
  1201 + T height; // height of the square shapce structure
  1202 +
  1203 + while (!done) {
  1204 + height = need_l / (2 * n); // calculate the height
  1205 + if (height > 2 * radius) {
  1206 + done = true;
  1207 + }
  1208 + else {
  1209 + n--;
  1210 + }
  1211 + }
  1212 +
  1213 + return n;
  1214 + }
  1215 +
  1216 + // build square connections
  1217 + 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) {
  1218 +
  1219 + int coef_up = (up) ? 1 : -1; // y coefficient
  1220 + int coef_left = (left) ? 1 : -1; // x coefficient
  1221 + int coef_z = (z) ? 1 : -1; // z coefficient
  1222 + int inverse = 1; // inverse flag
  1223 + stim::vec3<T> cor_v; // corner vertex
  1224 +
  1225 + stim::vec3<T> tmp_v;
  1226 + if (feeder == 1)
  1227 + tmp_v = inlet[i].V[0];
  1228 + else if (feeder == 0)
  1229 + tmp_v = outlet[i].V[0];
  1230 +
  1231 + // check whether the height of connections is acceptable
  1232 + if (height > threshold) { // acceptable
  1233 + // re-calculate height
  1234 + if (ratio > 0.0f && ratio <= 1.0f) { // use fragment if set
  1235 + cor_v = tmp_v + stim::vec3<T>(-coef_left * origin_l, 0, 0); // get the original corner vertex
  1236 + desire_l = desire_l - origin_l * (1.0f - ratio / 1.0f);
  1237 + origin_l = (T)origin_l * ratio / 1.0f;
  1238 + n = find_number_square(origin_l, desire_l);
  1239 + }
  1240 + height = (desire_l - (1 + 2 * n) * origin_l) / std::pow(2 * n, 2);
  1241 + // check whether the connections are good
  1242 + while (height > threshold) {
  1243 + n++;
  1244 + width = (T)(origin_l) / (2 * n);
  1245 + height = (desire_l - (1 + 2 * n) * origin_l) / std::pow(2 * n, 2);
  1246 + // check whether it appears overlap
  1247 + if (width < times * radius) {
  1248 + n--;
  1249 + width = (T)(origin_l) / (2 * n);
  1250 + height = (desire_l - (1 + 2 * n) * origin_l) / std::pow(2 * n, 2);
  1251 + break;
  1252 + }
  1253 + }
  1254 + while (height < times * radius) {
  1255 + n--;
  1256 + width = (T)(origin_l) / (2 * n);
  1257 + height = (desire_l - (1 + 2 * n) * origin_l) / std::pow(2 * n, 2);
  1258 + }
  1259 +
  1260 + // cube-like structure construction
  1261 + for (int j = 0; j < n; j++) {
  1262 + // "up"
  1263 + for (int k = 0; k < n; k++) {
  1264 + // in
  1265 + tmp_v = tmp_v + stim::vec3<T>(0, 0, coef_z * height);
  1266 + if (feeder == 1)
  1267 + inlet[i].V.push_back(tmp_v);
  1268 + else if (feeder == 0)
  1269 + outlet[i].V.push_back(tmp_v);
  1270 + // "up"
  1271 + tmp_v = tmp_v + stim::vec3<T>(0, inverse * coef_up * width, 0);
  1272 + if (feeder == 1)
  1273 + inlet[i].V.push_back(tmp_v);
  1274 + else if (feeder == 0)
  1275 + outlet[i].V.push_back(tmp_v);
  1276 + // out
  1277 + tmp_v = tmp_v + stim::vec3<T>(0, 0, -coef_z * height);
  1278 + if (feeder == 1)
  1279 + inlet[i].V.push_back(tmp_v);
  1280 + else if (feeder == 0)
  1281 + outlet[i].V.push_back(tmp_v);
  1282 + // "down"
  1283 + tmp_v = tmp_v + stim::vec3<T>(0, inverse * coef_up * width, 0);
  1284 + if (feeder == 1)
  1285 + inlet[i].V.push_back(tmp_v);
  1286 + else if (feeder == 0)
  1287 + outlet[i].V.push_back(tmp_v);
  1288 + }
  1289 +
  1290 + // "left"
  1291 + tmp_v = tmp_v + stim::vec3<T>(-coef_left * width, 0, 0);
  1292 + if (feeder == 1)
  1293 + inlet[i].V.push_back(tmp_v);
  1294 + else if (feeder == 0)
  1295 + outlet[i].V.push_back(tmp_v);
  1296 +
  1297 + if (inverse == 1) // revert inverse
  1298 + inverse = -1;
  1299 + else
  1300 + inverse = 1;
  1301 +
  1302 + // "down"
  1303 + for (int k = 0; k < n; k++) {
  1304 +
  1305 + tmp_v = tmp_v + stim::vec3<T>(0, 0, coef_z * height);
  1306 + if (feeder == 1)
  1307 + inlet[i].V.push_back(tmp_v);
  1308 + else if (feeder == 0)
  1309 + outlet[i].V.push_back(tmp_v);
  1310 +
  1311 + tmp_v = tmp_v + stim::vec3<T>(0, inverse * coef_up * width, 0);
  1312 + if (feeder == 1)
  1313 + inlet[i].V.push_back(tmp_v);
  1314 + else if (feeder == 0)
  1315 + outlet[i].V.push_back(tmp_v);
  1316 +
  1317 + tmp_v = tmp_v + stim::vec3<T>(0, 0, -coef_z * height);
  1318 + if (feeder == 1)
  1319 + inlet[i].V.push_back(tmp_v);
  1320 + else if (feeder == 0)
  1321 + outlet[i].V.push_back(tmp_v);
  1322 +
  1323 + tmp_v = tmp_v + stim::vec3<T>(0, inverse * coef_up * width, 0);
  1324 + if (feeder == 1)
  1325 + inlet[i].V.push_back(tmp_v);
  1326 + else if (feeder == 0)
  1327 + outlet[i].V.push_back(tmp_v);
  1328 + }
  1329 +
  1330 + // "left"
  1331 + tmp_v = tmp_v + stim::vec3<T>(-coef_left * width, 0, 0);
  1332 + if (feeder == 1)
  1333 + inlet[i].V.push_back(tmp_v);
  1334 + else if (feeder == 0)
  1335 + outlet[i].V.push_back(tmp_v);
  1336 +
  1337 + if (inverse == 1) // revert inverse
  1338 + inverse = -1;
  1339 + else
  1340 + inverse = 1;
  1341 + }
  1342 + // if use fragment to do square wave connection, need to push_back the corner vertex
  1343 + if (ratio > 0.0f && ratio <= 1.0f) {
  1344 + if (feeder == 1)
  1345 + inlet[i].V.push_back(cor_v);
  1346 + else if (feeder == 0)
  1347 + outlet[i].V.push_back(cor_v);
  1348 + }
  1349 + }
  1350 + else {
  1351 + for (int j = 0; j < n; j++) {
  1352 +
  1353 + // move in Z-shape
  1354 + tmp_v = tmp_v + stim::vec3<T>(0, coef_up * height, 0);
  1355 + if (feeder == 1)
  1356 + inlet[i].V.push_back(tmp_v);
  1357 + else if (feeder == 0)
  1358 + outlet[i].V.push_back(tmp_v);
  1359 +
  1360 + tmp_v = tmp_v + stim::vec3<T>(-coef_left * width, 0, 0);
  1361 + if (feeder == 1)
  1362 + inlet[i].V.push_back(tmp_v);
  1363 + else if (feeder == 0)
  1364 + outlet[i].V.push_back(tmp_v);
  1365 +
  1366 + tmp_v = tmp_v + stim::vec3<T>(0, -coef_up * height, 0);
  1367 + if (feeder == 1)
  1368 + inlet[i].V.push_back(tmp_v);
  1369 + else if (feeder == 0)
  1370 + outlet[i].V.push_back(tmp_v);
  1371 +
  1372 + tmp_v = tmp_v + stim::vec3<T>(-coef_left * width, 0, 0);
  1373 + if (feeder == 1)
  1374 + inlet[i].V.push_back(tmp_v);
  1375 + else if (feeder == 0)
  1376 + outlet[i].V.push_back(tmp_v);
  1377 + }
  1378 + }
  1379 + }
  1380 +
1191 // automatically modify bridge to make it feasible 1381 // automatically modify bridge to make it feasible
1192 - void modify_synthetic_connection(T viscosity, T rou, T radii = 5.0f) { 1382 + void modify_synthetic_connection(T viscosity, T rou, bool H, T threshold, T ratio = 0.0f, T radius = 5.0f) {
1193 1383
1194 glDeleteLists(dlist, 1); // delete display list for modify 1384 glDeleteLists(dlist, 1); // delete display list for modify
1195 glDeleteLists(dlist + 1, 1); 1385 glDeleteLists(dlist + 1, 1);
1196 1386
1197 - // because of radii change at the port vertex, there will be a pressure drop at that port 1387 + // because of radius change at the port vertex, there will be a pressure drop at that port
1198 // it follows the bernoulli equation 1388 // it follows the bernoulli equation
1199 // p1 + 1/2*rou*v1^2 + rou*g*h1 = p2 + 1/2*rou*v2^2 + rou*g*h2 1389 // p1 + 1/2*rou*v1^2 + rou*g*h1 = p2 + 1/2*rou*v2^2 + rou*g*h2
1200 // Q1 = Q2 -> v1*r1^2 = v2*r2^2 1390 // Q1 = Q2 -> v1*r1^2 = v2*r2^2
@@ -1203,7 +1393,7 @@ namespace stim { @@ -1203,7 +1393,7 @@ namespace stim {
1203 for (unsigned i = 0; i < pendant_vertex.size(); i++) { 1393 for (unsigned i = 0; i < pendant_vertex.size(); i++) {
1204 idx = pendant_vertex[i]; 1394 idx = pendant_vertex[i];
1205 T tmp_v = get_velocity(idx); // velocity at that pendant vertex 1395 T tmp_v = get_velocity(idx); // velocity at that pendant vertex
1206 - T ar = get_radius(idx) / radii; 1396 + T ar = get_radius(idx) / radius;
1207 new_pressure[idx] = pressure[idx] + 1.0f / 2.0f * rou * std::pow(tmp_v, 2) * (1.0f - std::pow(ar, 4)); 1397 new_pressure[idx] = pressure[idx] + 1.0f / 2.0f * rou * std::pow(tmp_v, 2) * (1.0f - std::pow(ar, 4));
1208 } 1398 }
1209 1399
@@ -1213,217 +1403,302 @@ namespace stim { @@ -1213,217 +1403,302 @@ namespace stim {
1213 unsigned inlet_index; 1403 unsigned inlet_index;
1214 T tmp_p; 1404 T tmp_p;
1215 for (unsigned i = 0; i < inlet.size(); i++) { 1405 for (unsigned i = 0; i < inlet.size(); i++) {
1216 - tmp_p = new_pressure[inlet[i].v[0]] + ((8 * viscosity * inlet[i].l * inlet[i].Q) / ((float)stim::PI * std::pow(radii, 4))); 1406 + tmp_p = new_pressure[inlet[i].v[0]] + ((8 * viscosity * inlet[i].l * inlet[i].Q) / ((float)stim::PI * std::pow(radius, 4)));
1217 if (tmp_p > source_pressure) { 1407 if (tmp_p > source_pressure) {
1218 source_pressure = tmp_p; 1408 source_pressure = tmp_p;
1219 inlet_index = i; 1409 inlet_index = i;
1220 } 1410 }
1221 } 1411 }
1222 1412
1223 - // automatically modify inlet bridge to make it feasible  
1224 - bool upper = false; // flag indicates the whether the port is upper than main feeder  
1225 - bool invert = false; // there are two version of hilbert curve depends on starting position with respect to the cup  
1226 - T new_l;  
1227 - stim::vec3<T> bus_v; // the port point on the bus  
1228 - stim::vec3<T> mid_v; // the original corner point  
1229 - stim::vec3<T> tmp_v; // the pendant point  
1230 - int order = 0; // order of hilbert curve (iteration)  
1231 - for (unsigned i = 0; i < inlet.size(); i++) {  
1232 - if (i != inlet_index) {  
1233 - new_l = (source_pressure - new_pressure[inlet[i].v[0]]) * ((float)stim::PI * std::pow(radii, 4)) / (8 * viscosity * inlet[i].Q); 1413 + // find minimum pressure outlet port
  1414 + T end_pressure = FLT_MAX;
  1415 + unsigned outlet_index;
  1416 + for (unsigned i = 0; i < outlet.size(); i++) {
  1417 + tmp_p = new_pressure[outlet[i].v[0]] - ((8 * viscosity * outlet[i].l * outlet[i].Q) / ((float)stim::PI * std::pow(radius, 4)));
  1418 + if (tmp_p < end_pressure) {
  1419 + end_pressure = tmp_p;
  1420 + outlet_index = i;
  1421 + }
  1422 + }
1234 1423
1235 - if (inlet[i].V[2][1] > main_feeder[0][1]) { // check out upper side of lower side  
1236 - upper = true;  
1237 - invert = false;  
1238 - }  
1239 - else {  
1240 - upper = false;  
1241 - invert = true;  
1242 - } 1424 + // automatically modify inlet bridge using Hilbert curves
  1425 + if (H) {
  1426 + bool upper = false; // flag indicates the whether the port is upper than main feeder
  1427 + bool invert = false; // there are two version of hilbert curve depends on starting position with respect to the cup
  1428 + T new_l;
  1429 + stim::vec3<T> bus_v; // the port point on the bus
  1430 + stim::vec3<T> mid_v; // the original corner point
  1431 + stim::vec3<T> tmp_v; // the pendant point
  1432 + int order = 0; // order of hilbert curve (iteration)
  1433 + for (unsigned i = 0; i < inlet.size(); i++) {
  1434 + if (i != inlet_index) {
  1435 + new_l = (source_pressure - new_pressure[inlet[i].v[0]]) * ((float)stim::PI * std::pow(radius, 4)) / (8 * viscosity * inlet[i].Q);
1243 1436
1244 - T origin_l = (inlet[i].V[1] - inlet[i].V[2]).len();  
1245 - T desire_l = new_l - (inlet[i].V[0] - inlet[i].V[1]).len();  
1246 - find_hilbert_order(origin_l, desire_l, order);  
1247 -  
1248 - bus_v = inlet[i].V[0];  
1249 - mid_v = inlet[i].V[1];  
1250 - tmp_v = inlet[i].V[2];  
1251 - inlet[i].V.clear();  
1252 - inlet[i].V.push_back(tmp_v);  
1253 - inlet[i].l = new_l;  
1254 -  
1255 - if (desire_l - origin_l < 2 * radii) { // do not need to use hilbert curve, just increase the length by draging out  
1256 - T d = new_l - inlet[i].l;  
1257 - 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]);  
1258 - inlet[i].V.push_back(corner);  
1259 - 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]);  
1260 - inlet[i].V.push_back(corner);  
1261 - inlet[i].V.push_back(bus_v);  
1262 - }  
1263 - else {  
1264 - T fragment = (desire_l - origin_l) / ((std::pow(4, order) - 1) / (std::pow(2, order) - 1) - 1); // the length of the opening of cup  
1265 - T dl = fragment / (std::pow(2, order) - 1); // unit cup length  
1266 -  
1267 - if (dl > 2 * radii) { // if the radii is feasible  
1268 - if (upper)  
1269 - hilbert_curve(i, &tmp_v[0], order, dl, 1, invert, DOWN);  
1270 - else  
1271 - hilbert_curve(i, &tmp_v[0], order, dl, 1, invert, UP);  
1272 -  
1273 - if (tmp_v[0] != mid_v[0])  
1274 - inlet[i].V.push_back(mid_v); 1437 + if (inlet[i].V[2][1] > main_feeder[0][1]) { // check out upper side of lower side
  1438 + upper = true;
  1439 + invert = false;
  1440 + }
  1441 + else {
  1442 + upper = false;
  1443 + invert = true;
  1444 + }
  1445 +
  1446 + T origin_l = (inlet[i].V[1] - inlet[i].V[2]).len();
  1447 + T desire_l = new_l - (inlet[i].V[0] - inlet[i].V[1]).len();
  1448 + find_hilbert_order(origin_l, desire_l, order);
  1449 +
  1450 + bus_v = inlet[i].V[0];
  1451 + mid_v = inlet[i].V[1];
  1452 + tmp_v = inlet[i].V[2];
  1453 + inlet[i].V.clear();
  1454 + inlet[i].V.push_back(tmp_v);
  1455 + inlet[i].l = new_l;
  1456 +
  1457 + if (desire_l - origin_l < 2 * radius) { // do not need to use hilbert curve, just increase the length by draging out
  1458 + T d = new_l - inlet[i].l;
  1459 + 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]);
  1460 + inlet[i].V.push_back(corner);
  1461 + 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]);
  1462 + inlet[i].V.push_back(corner);
1275 inlet[i].V.push_back(bus_v); 1463 inlet[i].V.push_back(bus_v);
1276 } 1464 }
1277 - else { // if the radii is not feasible  
1278 - int count = 1;  
1279 - while (dl <= 2 * radii) {  
1280 - dl = origin_l / (std::pow(2, order - count) - 1);  
1281 - count++;  
1282 - }  
1283 - count--; 1465 + else {
  1466 + T fragment = (desire_l - origin_l) / ((std::pow(4, order) - 1) / (std::pow(2, order) - 1) - 1); // the length of the opening of cup
  1467 + T dl = fragment / (std::pow(2, order) - 1); // unit cup length
1284 1468
1285 - if (upper)  
1286 - hilbert_curve(i, &tmp_v[0], order - count, dl, 1, invert, DOWN);  
1287 - else  
1288 - hilbert_curve(i, &tmp_v[0], order - count, dl, 1, invert, UP); 1469 + if (dl > 2 * radius) { // if the radius is feasible
  1470 + if (upper)
  1471 + hilbert_curve(i, &tmp_v[0], order, dl, 1, invert, DOWN);
  1472 + else
  1473 + hilbert_curve(i, &tmp_v[0], order, dl, 1, invert, UP);
1289 1474
1290 - desire_l -= origin_l * ((std::pow(4, order - count) - 1) / (std::pow(2, order - count) - 1));  
1291 - origin_l = (bus_v - mid_v).len();  
1292 - desire_l += origin_l; 1475 + if (tmp_v[0] != mid_v[0])
  1476 + inlet[i].V.push_back(mid_v);
  1477 + inlet[i].V.push_back(bus_v);
  1478 + }
  1479 + else { // if the radius is not feasible
  1480 + int count = 1;
  1481 + while (dl <= 2 * radius) {
  1482 + dl = origin_l / (std::pow(2, order - count) - 1);
  1483 + count++;
  1484 + }
  1485 + count--;
  1486 +
  1487 + if (upper)
  1488 + hilbert_curve(i, &tmp_v[0], order - count, dl, 1, invert, DOWN);
  1489 + else
  1490 + hilbert_curve(i, &tmp_v[0], order - count, dl, 1, invert, UP);
  1491 +
  1492 + desire_l -= origin_l * ((std::pow(4, order - count) - 1) / (std::pow(2, order - count) - 1));
  1493 + origin_l = (bus_v - mid_v).len();
  1494 + desire_l += origin_l;
  1495 +
  1496 + find_hilbert_order(origin_l, desire_l, order);
  1497 +
  1498 + fragment = (desire_l - origin_l) / ((std::pow(4, order) - 1) / (std::pow(2, order) - 1) - 1);
  1499 + dl = fragment / (std::pow(2, order) - 1);
  1500 + if (dl < 2 * radius)
  1501 + std::cout << "infeasible connection between inlets!" << std::endl;
  1502 +
  1503 + if (upper)
  1504 + hilbert_curve(i, &tmp_v[0], order, dl, 1, !invert, LEFT);
  1505 + else
  1506 + hilbert_curve(i, &tmp_v[0], order, dl, 1, !invert, RIGHT);
  1507 +
  1508 + if (tmp_v[1] != bus_v[1])
  1509 + inlet[i].V.push_back(bus_v);
  1510 + }
  1511 + }
  1512 + std::reverse(inlet[i].V.begin(), inlet[i].V.end()); // from bus to pendant vertex
  1513 + }
  1514 + }
1293 1515
1294 - find_hilbert_order(origin_l, desire_l, order); 1516 + // automatically modify outlet bridge to make it feasible
  1517 + for (unsigned i = 0; i < outlet.size(); i++) {
  1518 + if (i != outlet_index) {
  1519 + new_l = (new_pressure[outlet[i].v[0]] - end_pressure) * ((float)stim::PI * std::pow(radius, 4)) / (8 * viscosity * outlet[i].Q);
1295 1520
1296 - fragment = (desire_l - origin_l) / ((std::pow(4, order) - 1) / (std::pow(2, order) - 1) - 1);  
1297 - dl = fragment / (std::pow(2, order) - 1);  
1298 - if (dl < 2 * radii)  
1299 - std::cout << "infeasible connection between inlets!" << std::endl; 1521 + if (outlet[i].V[2][1] > main_feeder[1][1]) {
  1522 + upper = true;
  1523 + invert = true;
  1524 + }
  1525 + else {
  1526 + upper = false;
  1527 + invert = false;
  1528 + }
1300 1529
1301 - if (upper)  
1302 - hilbert_curve(i, &tmp_v[0], order, dl, 1, !invert, LEFT);  
1303 - else  
1304 - hilbert_curve(i, &tmp_v[0], order, dl, 1, !invert, RIGHT); 1530 + T origin_l = (outlet[i].V[1] - outlet[i].V[2]).len();
  1531 + T desire_l = new_l - (outlet[i].V[0] - outlet[i].V[1]).len();
  1532 + find_hilbert_order(origin_l, desire_l, order);
  1533 +
  1534 + bus_v = outlet[i].V[0];
  1535 + mid_v = outlet[i].V[1];
  1536 + tmp_v = outlet[i].V[2];
  1537 + outlet[i].V.clear();
  1538 + outlet[i].V.push_back(tmp_v);
  1539 + outlet[i].l = new_l;
  1540 +
  1541 + if (desire_l - origin_l < 2 * radius) { // do not need to use hilbert curve, just increase the length by draging out
  1542 + T d = new_l - outlet[i].l;
  1543 + 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]);
  1544 + outlet[i].V.push_back(corner);
  1545 + 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]);
  1546 + outlet[i].V.push_back(corner);
  1547 + outlet[i].V.push_back(bus_v);
  1548 + }
  1549 + else {
  1550 + T fragment = (desire_l - origin_l) / ((std::pow(4, order) - 1) / (std::pow(2, order) - 1) - 1); // the length of the opening of cup
  1551 + T dl = fragment / (std::pow(2, order) - 1); // unit cup length
1305 1552
1306 - if (tmp_v[1] != bus_v[1])  
1307 - inlet[i].V.push_back(bus_v); 1553 + if (dl > 2 * radius) { // if the radius is feasible
  1554 + if (upper)
  1555 + hilbert_curve(i, &tmp_v[0], order, dl, 0, invert, DOWN);
  1556 + else
  1557 + hilbert_curve(i, &tmp_v[0], order, dl, 0, invert, UP);
  1558 +
  1559 + if (tmp_v[0] != mid_v[0])
  1560 + outlet[i].V.push_back(mid_v);
  1561 + outlet[i].V.push_back(bus_v);
  1562 + }
  1563 + else { // if the radius is not feasible
  1564 + int count = 1;
  1565 + while (dl <= 2 * radius) {
  1566 + dl = origin_l / (std::pow(2, order - count) - 1);
  1567 + count++;
  1568 + }
  1569 + count--;
  1570 +
  1571 + if (upper)
  1572 + hilbert_curve(i, &tmp_v[0], order - count, dl, 0, invert, DOWN);
  1573 + else
  1574 + hilbert_curve(i, &tmp_v[0], order - count, dl, 0, invert, UP);
  1575 +
  1576 + desire_l -= origin_l * ((std::pow(4, order - count) - 1) / (std::pow(2, order - count) - 1));
  1577 + origin_l = (bus_v - mid_v).len();
  1578 + desire_l += origin_l;
  1579 +
  1580 + find_hilbert_order(origin_l, desire_l, order);
  1581 +
  1582 + fragment = (desire_l - origin_l) / ((std::pow(4, order) - 1) / (std::pow(2, order) - 1) - 1);
  1583 + dl = fragment / (std::pow(2, order) - 1);
  1584 + if (dl < 2 * radius)
  1585 + std::cout << "infeasible connection between outlets!" << std::endl;
  1586 +
  1587 + if (upper)
  1588 + hilbert_curve(i, &tmp_v[0], order, dl, 0, !invert, LEFT);
  1589 + else
  1590 + hilbert_curve(i, &tmp_v[0], order, dl, 0, !invert, RIGHT);
  1591 +
  1592 + if (tmp_v[1] != bus_v[1])
  1593 + outlet[i].V.push_back(bus_v);
  1594 + }
1308 } 1595 }
  1596 + std::reverse(outlet[i].V.begin(), outlet[i].V.end());
1309 } 1597 }
1310 - std::reverse(inlet[i].V.begin(), inlet[i].V.end()); // from bus to pendant vertex  
1311 } 1598 }
1312 } 1599 }
  1600 + // automatically modify inlet bridge using square shape constructions
  1601 + else {
  1602 + bool upper; // flag indicates the connection is upper than the bus
  1603 + bool z; // flag indicates the connection direction along z-axis
  1604 + T new_l; // new length
  1605 + stim::vec3<T> bus_v; // the port point on the bus
  1606 + stim::vec3<T> mid_v; // the original corner point
  1607 + stim::vec3<T> tmp_v; // the pendant point
  1608 + int n;
  1609 + T width, height; // width and height of the square
1313 1610
1314 - // find minimum pressure outlet port  
1315 - source_pressure = FLT_MAX;  
1316 - unsigned outlet_index;  
1317 - for (unsigned i = 0; i < outlet.size(); i++) {  
1318 - tmp_p = new_pressure[outlet[i].v[0]] - ((8 * viscosity * outlet[i].l * outlet[i].Q) / ((float)stim::PI * std::pow(radii, 4)));  
1319 - if (tmp_p < source_pressure) {  
1320 - source_pressure = tmp_p;  
1321 - outlet_index = i;  
1322 - }  
1323 - } 1611 + for (unsigned i = 0; i < inlet.size(); i++) {
  1612 + if (i != inlet_index) {
  1613 + 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
1324 1614
1325 - // automatically modify outlet bridge to make it feasible  
1326 - for (unsigned i = 0; i < outlet.size(); i++) {  
1327 - if (i != outlet_index) {  
1328 - new_l = (new_pressure[outlet[i].v[0]] - source_pressure) * ((float)stim::PI * std::pow(radii, 4)) / (8 * viscosity * outlet[i].Q); 1615 + bus_v = inlet[i].V[0];
  1616 + mid_v = inlet[i].V[1];
  1617 + tmp_v = inlet[i].V[2];
1329 1618
1330 - if (outlet[i].V[2][1] > main_feeder[1][1]) {  
1331 - upper = true;  
1332 - invert = true;  
1333 - }  
1334 - else {  
1335 - upper = false;  
1336 - invert = false;  
1337 - } 1619 + if (inlet[i].V[2][1] > main_feeder[0][1]) // check out upper side of lower side
  1620 + upper = true;
  1621 + else
  1622 + upper = false;
1338 1623
1339 - T origin_l = (outlet[i].V[1] - outlet[i].V[2]).len();  
1340 - T desire_l = new_l - (outlet[i].V[0] - outlet[i].V[1]).len();  
1341 - find_hilbert_order(origin_l, desire_l, order);  
1342 -  
1343 - bus_v = outlet[i].V[0];  
1344 - mid_v = outlet[i].V[1];  
1345 - tmp_v = outlet[i].V[2];  
1346 - outlet[i].V.clear();  
1347 - outlet[i].V.push_back(tmp_v);  
1348 - outlet[i].l = new_l;  
1349 -  
1350 - if (desire_l - origin_l < 2 * radii) { // do not need to use hilbert curve, just increase the length by draging out  
1351 - T d = new_l - outlet[i].l;  
1352 - 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]);  
1353 - outlet[i].V.push_back(corner);  
1354 - 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]);  
1355 - outlet[i].V.push_back(corner);  
1356 - outlet[i].V.push_back(bus_v); 1624 + if (inlet[i].V[2][2] > main_feeder[0][2])
  1625 + z = true;
  1626 + else
  1627 + z = false;
  1628 +
  1629 + T origin_l = (inlet[i].V[1] - inlet[i].V[2]).len();
  1630 + T desire_l = new_l - (inlet[i].V[0] - inlet[i].V[1]).len();
  1631 + inlet[i].V.clear();
  1632 + inlet[i].V.push_back(tmp_v);
  1633 + inlet[i].l = new_l;
  1634 +
  1635 + n = find_number_square(origin_l, desire_l);
  1636 +
  1637 + width = (T)origin_l / (2 * n);
  1638 + height = (desire_l - origin_l) / (2 * n);
  1639 +
  1640 + build_square_connection(i, width, height, origin_l, desire_l, n, 1, threshold, z, true, upper, 10, ratio);
  1641 + inlet[i].V.push_back(bus_v);
  1642 +
  1643 + std::reverse(inlet[i].V.begin(), inlet[i].V.end()); // from bus to pendant vertex
1357 } 1644 }
1358 - else {  
1359 - T fragment = (desire_l - origin_l) / ((std::pow(4, order) - 1) / (std::pow(2, order) - 1) - 1); // the length of the opening of cup  
1360 - T dl = fragment / (std::pow(2, order) - 1); // unit cup length  
1361 -  
1362 - if (dl > 2 * radii) { // if the radii is feasible  
1363 - if (upper)  
1364 - hilbert_curve(i, &tmp_v[0], order, dl, 0, invert, DOWN);  
1365 - else  
1366 - hilbert_curve(i, &tmp_v[0], order, dl, 0, invert, UP);  
1367 -  
1368 - if (tmp_v[0] != mid_v[0])  
1369 - outlet[i].V.push_back(mid_v);  
1370 - outlet[i].V.push_back(bus_v);  
1371 - }  
1372 - else { // if the radii is not feasible  
1373 - int count = 1;  
1374 - while (dl <= 2 * radii) {  
1375 - dl = origin_l / (std::pow(2, order - count) - 1);  
1376 - count++;  
1377 - }  
1378 - count--; 1645 + }
  1646 +
  1647 + for (unsigned i = 0; i < outlet.size(); i++) {
  1648 + if (i != outlet_index) {
  1649 + 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
1379 1650
1380 - if (upper)  
1381 - hilbert_curve(i, &tmp_v[0], order - count, dl, 0, invert, DOWN);  
1382 - else  
1383 - hilbert_curve(i, &tmp_v[0], order - count, dl, 0, invert, UP); 1651 + bus_v = outlet[i].V[0];
  1652 + mid_v = outlet[i].V[1];
  1653 + tmp_v = outlet[i].V[2];
1384 1654
1385 - desire_l -= origin_l * ((std::pow(4, order - count) - 1) / (std::pow(2, order - count) - 1));  
1386 - origin_l = (bus_v - mid_v).len();  
1387 - desire_l += origin_l; 1655 + if (outlet[i].V[2][1] > main_feeder[1][1]) // check out upper side of lower side
  1656 + upper = true;
  1657 + else
  1658 + upper = false;
1388 1659
1389 - find_hilbert_order(origin_l, desire_l, order); 1660 + if (outlet[i].V[2][2] > main_feeder[1][2])
  1661 + z = true;
  1662 + else
  1663 + z = false;
1390 1664
1391 - fragment = (desire_l - origin_l) / ((std::pow(4, order) - 1) / (std::pow(2, order) - 1) - 1);  
1392 - dl = fragment / (std::pow(2, order) - 1);  
1393 - if (dl < 2 * radii)  
1394 - std::cout << "infeasible connection between outlets!" << std::endl; 1665 + T origin_l = (outlet[i].V[1] - outlet[i].V[2]).len();
  1666 + T desire_l = new_l - (outlet[i].V[0] - outlet[i].V[1]).len();
  1667 + outlet[i].V.clear();
  1668 + outlet[i].V.push_back(tmp_v);
  1669 + outlet[i].l = new_l;
1395 1670
1396 - if (upper)  
1397 - hilbert_curve(i, &tmp_v[0], order, dl, 0, !invert, LEFT);  
1398 - else  
1399 - hilbert_curve(i, &tmp_v[0], order, dl, 0, !invert, RIGHT); 1671 + n = find_number_square(origin_l, desire_l);
1400 1672
1401 - if (tmp_v[1] != bus_v[1])  
1402 - outlet[i].V.push_back(bus_v);  
1403 - } 1673 + width = (T)origin_l / (2 * n);
  1674 + height = (desire_l - origin_l) / (2 * n);
  1675 +
  1676 + build_square_connection(i, width, height, origin_l, desire_l, n, 0, threshold, z, false, upper, 10, ratio);
  1677 + outlet[i].V.push_back(bus_v);
  1678 +
  1679 + std::reverse(outlet[i].V.begin(), outlet[i].V.end()); // from bus to pendant vertex
1404 } 1680 }
1405 - std::reverse(outlet[i].V.begin(), outlet[i].V.end());  
1406 } 1681 }
1407 } 1682 }
1408 } 1683 }
1409 1684
1410 // check current bridge to see feasibility 1685 // check current bridge to see feasibility
1411 - void check_synthetic_connection(T viscosity, T radii = 5.0f) { 1686 + void check_synthetic_connection(T viscosity, T radius = 5.0f) {
1412 1687
1413 T eps = 0.01f; 1688 T eps = 0.01f;
1414 - T source_pressure = pressure[inlet[0].v[0]] + (8 * viscosity * inlet[0].l * inlet[0].Q) / ((float)stim::PI * std::pow(radii, 4)); 1689 + T source_pressure = pressure[inlet[0].v[0]] + (8 * viscosity * inlet[0].l * inlet[0].Q) / ((float)stim::PI * std::pow(radius, 4));
1415 T tmp_p; 1690 T tmp_p;
1416 for (unsigned i = 1; i < inlet.size(); i++) { 1691 for (unsigned i = 1; i < inlet.size(); i++) {
1417 - tmp_p = pressure[inlet[i].v[0]] + (8 * viscosity * inlet[i].l * inlet[i].Q) / ((float)stim::PI * std::pow(radii, 4)); 1692 + tmp_p = pressure[inlet[i].v[0]] + (8 * viscosity * inlet[i].l * inlet[i].Q) / ((float)stim::PI * std::pow(radius, 4));
1418 T delta = fabs(tmp_p - source_pressure); 1693 T delta = fabs(tmp_p - source_pressure);
1419 if (delta > eps) { 1694 if (delta > eps) {
1420 std::cout << "Nonfeasible connection!" << std::endl; 1695 std::cout << "Nonfeasible connection!" << std::endl;
1421 break; 1696 break;
1422 } 1697 }
1423 } 1698 }
1424 - source_pressure = pressure[outlet[0].v[0]] - (8 * viscosity * outlet[0].l * outlet[0].Q) / ((float)stim::PI * std::pow(radii, 4)); 1699 + source_pressure = pressure[outlet[0].v[0]] - (8 * viscosity * outlet[0].l * outlet[0].Q) / ((float)stim::PI * std::pow(radius, 4));
1425 for (unsigned i = 1; i < outlet.size(); i++) { 1700 for (unsigned i = 1; i < outlet.size(); i++) {
1426 - tmp_p = pressure[outlet[i].v[0]] - (8 * viscosity * outlet[i].l * outlet[i].Q) / ((float)stim::PI * std::pow(radii, 4)); 1701 + tmp_p = pressure[outlet[i].v[0]] - (8 * viscosity * outlet[i].l * outlet[i].Q) / ((float)stim::PI * std::pow(radius, 4));
1427 T delta = fabs(tmp_p - source_pressure); 1702 T delta = fabs(tmp_p - source_pressure);
1428 if (delta > eps) { 1703 if (delta > eps) {
1429 std::cout << "Nonfeasible connection!" << std::endl; 1704 std::cout << "Nonfeasible connection!" << std::endl;
@@ -1436,7 +1711,7 @@ namespace stim { @@ -1436,7 +1711,7 @@ namespace stim {
1436 // prepare for image stack 1711 // prepare for image stack
1437 void preparation(T &Xl, T &Xr, T &Yt, T &Yb, T &Z, T length = 210.0f, T height = 10.0f) { 1712 void preparation(T &Xl, T &Xr, T &Yt, T &Yb, T &Z, T length = 210.0f, T height = 10.0f) {
1438 1713
1439 - T max_radii = 0.0f; 1714 + T max_radius = 0.0f;
1440 T top = FLT_MIN; 1715 T top = FLT_MIN;
1441 T bottom = FLT_MAX; 1716 T bottom = FLT_MAX;
1442 1717
@@ -1519,17 +1794,17 @@ namespace stim { @@ -1519,17 +1794,17 @@ namespace stim {
1519 top = A[i].c[1]; 1794 top = A[i].c[1];
1520 if (A[i].c[1] < bottom) 1795 if (A[i].c[1] < bottom)
1521 bottom = A[i].c[1]; 1796 bottom = A[i].c[1];
1522 - if (A[i].r > max_radii)  
1523 - max_radii = A[i].r; 1797 + if (A[i].r > max_radius)
  1798 + max_radius = A[i].r;
1524 } 1799 }
1525 1800
1526 Yt = top; // top bound y coordinate 1801 Yt = top; // top bound y coordinate
1527 Yb = bottom; // bottom bound y coordinate 1802 Yb = bottom; // bottom bound y coordinate
1528 - Z = (bb.B[2] - bb.A[2] + 2 * max_radii); // bounding box width(along z-axis) 1803 + Z = (bb.B[2] - bb.A[2] + 2 * max_radius); // bounding box width(along z-axis)
1529 } 1804 }
1530 1805
1531 /// making image stack main function 1806 /// making image stack main function
1532 - void make_image_stack(T dx, T dy, T dz, std::string stackdir, T radii = 5.0f) { 1807 + void make_image_stack(T dx, T dy, T dz, std::string stackdir, T radius = 5.0f) {
1533 1808
1534 /// preparation for making image stack 1809 /// preparation for making image stack
1535 T X, Xl, Xr, Y, Yt, Yb, Z; 1810 T X, Xl, Xr, Y, Yt, Yb, Z;
@@ -54,10 +54,12 @@ float zoom_factor = 10.0f; // zooming factor @@ -54,10 +54,12 @@ float zoom_factor = 10.0f; // zooming factor
54 float border_factor = 20.0f; // border 54 float border_factor = 20.0f; // border
55 float radii_factor = 1.0f; // radii changing factor 55 float radii_factor = 1.0f; // radii changing factor
56 GLint subdivision = 20; // slices and stacks 56 GLint subdivision = 20; // slices and stacks
57 -float default_radii = 5.0f; // default radii of network vertex 57 +float default_radius = 5.0f; // default radii of network vertex
58 float delta = 0.01f; // small discrepancy 58 float delta = 0.01f; // small discrepancy
59 float eps = 20.0f; // epsilon threshold 59 float eps = 20.0f; // epsilon threshold
60 float max_pressure = 0.0f; // maximum pressure that the channel can bear 60 float max_pressure = 0.0f; // maximum pressure that the channel can bear
  61 +float height_threshold = 100.0f; // connection height constraint
  62 +float fragment_ratio = 0.0f; // fragment ratio
61 63
62 // glut event parameters 64 // glut event parameters
63 int mouse_x; // window x-coordinate 65 int mouse_x; // window x-coordinate
@@ -73,6 +75,8 @@ unsigned pressure_index; // the index of vertex that is clicked @@ -73,6 +75,8 @@ unsigned pressure_index; // the index of vertex that is clicked
73 // build inlet/outlet parameters 75 // build inlet/outlet parameters
74 bool build_inlet_outlet = false; // flag indicates building inlets and outlets 76 bool build_inlet_outlet = false; // flag indicates building inlets and outlets
75 bool modified_bridge = false; // flag indicates having modified inlet/outlet connection 77 bool modified_bridge = false; // flag indicates having modified inlet/outlet connection
  78 +bool hilbert_curve = false; // flag indicates enabling hilbert curves constructions
  79 +bool change_fragment = false; // flag indicates changing fragment for square wave connections
76 80
77 // manufacture parameters 81 // manufacture parameters
78 bool manufacture = false; // flag indicates manufacture mode 82 bool manufacture = false; // flag indicates manufacture mode
@@ -89,7 +93,7 @@ inline void get_background() { @@ -89,7 +93,7 @@ inline void get_background() {
89 // set the initial radii 93 // set the initial radii
90 flow.init(num_edge, num_vertex); // initialize flow object 94 flow.init(num_edge, num_vertex); // initialize flow object
91 for (unsigned i = 0; i < num_edge; i++) 95 for (unsigned i = 0; i < num_edge; i++)
92 - flow.set_r(i, default_radii); 96 + flow.set_r(i, default_radius);
93 } 97 }
94 98
95 // convert from window coordinates to world coordinates 99 // convert from window coordinates to world coordinates
@@ -116,6 +120,7 @@ inline void window_to_world(GLdouble &amp;x, GLdouble &amp;y, GLdouble &amp;z) { @@ -116,6 +120,7 @@ inline void window_to_world(GLdouble &amp;x, GLdouble &amp;y, GLdouble &amp;z) {
116 // initialize flow object 120 // initialize flow object
117 void flow_initialize() { 121 void flow_initialize() {
118 122
  123 + flow.set = true;
119 stim::vec3<float> center = bb.center(); 124 stim::vec3<float> center = bb.center();
120 125
121 for (unsigned i = 0; i < pendant_vertex.size(); i++) { 126 for (unsigned i = 0; i < pendant_vertex.size(); i++) {
@@ -187,6 +192,7 @@ void glut_render() { @@ -187,6 +192,7 @@ void glut_render() {
187 192
188 glEnable(GL_DEPTH_TEST); 193 glEnable(GL_DEPTH_TEST);
189 glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); 194 glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  195 + glClearColor(1.0f, 1.0f, 1.0f, 1.0f);
190 glut_projection(); 196 glut_projection();
191 glut_modelview(); 197 glut_modelview();
192 198
@@ -238,7 +244,7 @@ void glut_render() { @@ -238,7 +244,7 @@ void glut_render() {
238 glFlush(); 244 glFlush();
239 245
240 // pressure bar text 246 // pressure bar text
241 - glColor3f(1.0f, 1.0f, 1.0f); 247 + glColor3f(0.0f, 0.0f, 0.0f);
242 glRasterPos2f(0.0f, vY - border_factor); 248 glRasterPos2f(0.0f, vY - border_factor);
243 std::stringstream ss_p; 249 std::stringstream ss_p;
244 ss_p << "Pressure Bar"; 250 ss_p << "Pressure Bar";
@@ -251,7 +257,7 @@ void glut_render() { @@ -251,7 +257,7 @@ void glut_render() {
251 glRasterPos2f((border_factor * 1.5f), (border_factor + i * step)); 257 glRasterPos2f((border_factor * 1.5f), (border_factor + i * step));
252 std::stringstream ss_n; 258 std::stringstream ss_n;
253 ss_n << (float)i * max_pressure / 10; 259 ss_n << (float)i * max_pressure / 10;
254 - glutBitmapString(GLUT_BITMAP_TIMES_ROMAN_10, (const unsigned char*)(ss_n.str().c_str())); 260 + glutBitmapString(GLUT_BITMAP_HELVETICA_18, (const unsigned char*)(ss_n.str().c_str()));
255 } 261 }
256 glPopMatrix(); 262 glPopMatrix();
257 glMatrixMode(GL_PROJECTION); 263 glMatrixMode(GL_PROJECTION);
@@ -286,7 +292,7 @@ void glut_render() { @@ -286,7 +292,7 @@ void glut_render() {
286 glFlush(); 292 glFlush();
287 293
288 // pressure bar text 294 // pressure bar text
289 - glColor3f(1.0f, 1.0f, 1.0f); 295 + glColor3f(0.0f, 0.0f, 0.0f);
290 glRasterPos2f(0.0f, vY - border_factor); 296 glRasterPos2f(0.0f, vY - border_factor);
291 std::stringstream ss_p; 297 std::stringstream ss_p;
292 ss_p << "Velocity range"; 298 ss_p << "Velocity range";
@@ -299,13 +305,58 @@ void glut_render() { @@ -299,13 +305,58 @@ void glut_render() {
299 glRasterPos2f(border_factor * 1.5f, border_factor + i * step); 305 glRasterPos2f(border_factor * 1.5f, border_factor + i * step);
300 std::stringstream ss_n; 306 std::stringstream ss_n;
301 ss_n << min_v + i * (max_v - min_v) / 10; 307 ss_n << min_v + i * (max_v - min_v) / 10;
302 - glutBitmapString(GLUT_BITMAP_TIMES_ROMAN_10, (const unsigned char*)(ss_n.str().c_str())); 308 + glutBitmapString(GLUT_BITMAP_HELVETICA_18, (const unsigned char*)(ss_n.str().c_str()));
303 } 309 }
304 glPopMatrix(); 310 glPopMatrix();
305 glMatrixMode(GL_PROJECTION); 311 glMatrixMode(GL_PROJECTION);
306 glPopMatrix(); 312 glPopMatrix();
307 } 313 }
308 - 314 +
  315 + // bring up a ratio bar on the left
  316 + if (change_fragment) {
  317 +
  318 + glMatrixMode(GL_PROJECTION); // set up the 2d viewport for mode text printing
  319 + glPushMatrix();
  320 + glLoadIdentity();
  321 + vX = glutGet(GLUT_WINDOW_WIDTH); // get the current window width
  322 + vY = glutGet(GLUT_WINDOW_HEIGHT); // get the current window height
  323 + glViewport(0, 0, vX, vY); // locate to left bottom corner
  324 + gluOrtho2D(0, vX, 0, vY); // define othogonal aspect
  325 +
  326 + glMatrixMode(GL_MODELVIEW);
  327 + glPushMatrix();
  328 + glLoadIdentity();
  329 +
  330 + glLineWidth(border_factor);
  331 + glBegin(GL_LINES);
  332 + glColor3f(0.0, 0.0, 1.0); // blue to red
  333 + glVertex2f(border_factor, border_factor);
  334 + glColor3f(1.0, 0.0, 0.0);
  335 + glVertex2f(border_factor, (vY - 2.0f * border_factor));
  336 + glEnd();
  337 + glFlush();
  338 +
  339 + // pressure bar text
  340 + glColor3f(0.0f, 0.0f, 0.0f);
  341 + glRasterPos2f(0.0f, vY - border_factor);
  342 + std::stringstream ss_p;
  343 + ss_p << "Ratio bar";
  344 + glutBitmapString(GLUT_BITMAP_HELVETICA_18, (const unsigned char*)(ss_p.str().c_str()));
  345 +
  346 + // pressure range text
  347 + float step = vY - 3.0f * border_factor;
  348 + step /= 10;
  349 + for (unsigned i = 0; i < 11; i++) {
  350 + glRasterPos2f((border_factor * 1.5f), (border_factor + i * step));
  351 + std::stringstream ss_n;
  352 + ss_n << (float)i * 1.0f / 10;
  353 + glutBitmapString(GLUT_BITMAP_HELVETICA_18, (const unsigned char*)(ss_n.str().c_str()));
  354 + }
  355 + glPopMatrix();
  356 + glMatrixMode(GL_PROJECTION);
  357 + glPopMatrix();
  358 + }
  359 +
309 glutSwapBuffers(); 360 glutSwapBuffers();
310 } 361 }
311 362
@@ -318,7 +369,8 @@ void glut_menu(int value) { @@ -318,7 +369,8 @@ void glut_menu(int value) {
318 build_inlet_outlet = false; 369 build_inlet_outlet = false;
319 manufacture = false; 370 manufacture = false;
320 modified_bridge = false; 371 modified_bridge = false;
321 - flow_initialize(); 372 + if (!flow.set)
  373 + flow_initialize();
322 flow_stable_state(); // main function of solving the linear system 374 flow_stable_state(); // main function of solving the linear system
323 flow.print_flow(); 375 flow.print_flow();
324 376
@@ -331,7 +383,7 @@ void glut_menu(int value) { @@ -331,7 +383,7 @@ void glut_menu(int value) {
331 manufacture = false; 383 manufacture = false;
332 if (!modified_bridge) { 384 if (!modified_bridge) {
333 flow.set_main_feeder(); 385 flow.set_main_feeder();
334 - flow.build_synthetic_connection(u); 386 + flow.build_synthetic_connection(u, default_radius);
335 } 387 }
336 388
337 glut_set_menu(num, 3); 389 glut_set_menu(num, 3);
@@ -341,7 +393,7 @@ void glut_menu(int value) { @@ -341,7 +393,7 @@ void glut_menu(int value) {
341 simulation = false; 393 simulation = false;
342 build_inlet_outlet = false; 394 build_inlet_outlet = false;
343 manufacture = true; 395 manufacture = true;
344 - flow.check_synthetic_connection(u); 396 + flow.check_synthetic_connection(u, default_radius);
345 } 397 }
346 398
347 glutPostRedisplay(); 399 glutPostRedisplay();
@@ -396,6 +448,13 @@ void glut_mouse(int button, int state, int x, int y) { @@ -396,6 +448,13 @@ void glut_mouse(int button, int state, int x, int y) {
396 flow.print_flow(); 448 flow.print_flow();
397 } 449 }
398 } 450 }
  451 + else if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN && modified_bridge && change_fragment) {
  452 + if (y > 2 * border_factor || y < vY - border_factor) { // within the ratio bar range
  453 + fragment_ratio = (float)(vY - y - border_factor) / ((float)vY - border_factor) * 1.0f;
  454 + flow.modify_synthetic_connection(u, rou, hilbert_curve, height_threshold, fragment_ratio, default_radius);
  455 + change_fragment = false;
  456 + }
  457 + }
399 } 458 }
400 459
401 // define camera move based on mouse wheel move 460 // define camera move based on mouse wheel move
@@ -419,7 +478,7 @@ void glut_wheel(int wheel, int direction, int x, int y) { @@ -419,7 +478,7 @@ void glut_wheel(int wheel, int direction, int x, int y) {
419 tmp_r = flow.get_radius(pressure_index); 478 tmp_r = flow.get_radius(pressure_index);
420 tmp_r -= radii_factor; 479 tmp_r -= radii_factor;
421 if (tmp_r <= 0) 480 if (tmp_r <= 0)
422 - tmp_r = default_radii; 481 + tmp_r = default_radius;
423 } 482 }
424 flow.set_radius(pressure_index, tmp_r); 483 flow.set_radius(pressure_index, tmp_r);
425 flow_stable_state(); 484 flow_stable_state();
@@ -464,12 +523,16 @@ void glut_keyboard(unsigned char key, int x, int y) { @@ -464,12 +523,16 @@ void glut_keyboard(unsigned char key, int x, int y) {
464 #endif 523 #endif
465 } 524 }
466 else if (build_inlet_outlet && !modified_bridge) { 525 else if (build_inlet_outlet && !modified_bridge) {
467 - flow.modify_synthetic_connection(u, rou);  
468 modified_bridge = true; 526 modified_bridge = true;
  527 +
  528 + if (hilbert_curve)
  529 + flow.modify_synthetic_connection(u, rou, hilbert_curve, height_threshold);
  530 + else
  531 + change_fragment = true;
469 } 532 }
470 break; 533 break;
471 } 534 }
472 - 535 +
473 glutPostRedisplay(); 536 glutPostRedisplay();
474 } 537 }
475 538
@@ -510,6 +573,7 @@ int main(int argc, char* argv[]) { @@ -510,6 +573,7 @@ int main(int argc, char* argv[]) {
510 args.add("maxpress", "maximum allowed pressure in g / units / s^2, default 2 is for blood when units = um", "2", "real value > 0"); 573 args.add("maxpress", "maximum allowed pressure in g / units / s^2, default 2 is for blood when units = um", "2", "real value > 0");
511 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"); 574 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");
512 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"); 575 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");
  576 + args.add("hilbert", "enable hilbert curves connections", "0", "value 1 for enablement");
513 args.add("stackres", "spacing between pixel samples in each dimension(in units/pixel)", "1 1 1", "real value > 0"); 577 args.add("stackres", "spacing between pixel samples in each dimension(in units/pixel)", "1 1 1", "real value > 0");
514 args.add("stackdir", "set the directory of the output image stack", "", "any existing directory (ex. /home/name/network)"); 578 args.add("stackdir", "set the directory of the output image stack", "", "any existing directory (ex. /home/name/network)");
515 579
@@ -540,6 +604,9 @@ int main(int argc, char* argv[]) { @@ -540,6 +604,9 @@ int main(int argc, char* argv[]) {
540 // normally the blood density in capillaries: 1060 kg/m^3 = 1.06*10^-12 g/um^3 604 // normally the blood density in capillaries: 1060 kg/m^3 = 1.06*10^-12 g/um^3
541 rou = args["rou"].as_float(); 605 rou = args["rou"].as_float();
542 606
  607 + // check whether to enable hilbert curves or not
  608 + hilbert_curve = args["hilbert"].as_int();
  609 +
543 // get the vexel and image stack size 610 // get the vexel and image stack size
544 dx = args["stackres"].as_float(0); 611 dx = args["stackres"].as_float(0);
545 dy = args["stackres"].as_float(1); 612 dy = args["stackres"].as_float(1);