Jiaming Guo / flow3

Commit f4105b89553b247a70fe04555ee2a41a792899a8

Authored by Jiaming Guo
1 parent 9e60c6a7

add new function: smooth colormap, manually change direct connection...etc

Showing 2 changed files with 573 additions and 46 deletions   Show diff stats
Inline Side-by-side
... ... @@ -238,6 +238,8 @@ namespace stim {
238 238 std::vector<typename stim::triple<unsigned, unsigned, float> > Q; // volume flow rate
239 239 std::vector<T> QQ; // Q' vector
240 240 std::vector<T> pressure; // final pressure
  241 + std::vector<typename std::vector<typename stim::vec3<T> > > in_backup; // inlet connection back up
  242 + std::vector<typename std::vector<typename stim::vec3<T> > > out_backup;
241 243  
242 244 public:
243 245  
... ... @@ -254,6 +256,10 @@ namespace stim {
254 256 std::vector<typename stim::cone<T> > B; // cone(cylinder) model for making image stack
255 257 std::vector<typename stim::cuboid<T> > CU; // cuboid model for making image stack
256 258 stim::gl_aaboundingbox<T> bb; // bounding box
  259 + std::vector<bool> inlet_feasibility; // list of flags indicate whether one inlet connection is feasible
  260 + std::vector<bool> outlet_feasibility;
  261 + std::vector<typename std::pair<stim::vec3<T>, stim::vec3<T> > > inbb; // inlet connection bounding box
  262 + std::vector<typename std::pair<stim::vec3<T>, stim::vec3<T> > > outbb; // outlet connection bounding box
257 263  
258 264 flow() {} // default constructor
259 265 ~flow() {
... ... @@ -804,8 +810,8 @@ namespace stim {
804 810 std::vector<typename stim::vec3<float> > cp2(subdivision + 1);
805 811 for (unsigned i = 0; i < num_edge; i++) { // for every edge
806 812 if (i == index) { // render in tranparency for direction indication
807   - glEnable(GL_BLEND); // enable color blend
808   - glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); // set blend function
  813 + glEnable(GL_BLEND); // enable color blend
  814 + glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); // set blend function
809 815 glDisable(GL_DEPTH_TEST);
810 816 glColor4f((float)color[i * 3 + 0] / 255, (float)color[i * 3 + 1] / 255, (float)color[i * 3 + 2] / 255, 0.5f);
811 817 }
... ... @@ -1052,19 +1058,33 @@ namespace stim {
1052 1058 }
1053 1059  
1054 1060 // draw the bridge as lines
1055   - void line_bridge() {
  1061 + void line_bridge(bool &redisplay) {
  1062 +
  1063 + if (redisplay)
  1064 + glDeleteLists(dlist, 1);
  1065 + redisplay = false;
1056 1066  
1057 1067 if (!glIsList(dlist)) {
1058 1068 dlist = glGenLists(1);
1059 1069 glNewList(dlist, GL_COMPILE);
1060   - glColor3f(0.0f, 0.0f, 0.0f);
  1070 +
  1071 + glLineWidth(50);
1061 1072 for (unsigned i = 0; i < inlet.size(); i++) {
  1073 + if (inlet_feasibility[i])
  1074 + glColor3f(0.0f, 0.0f, 0.0f); // feasible -> black
  1075 + else
  1076 + glColor3f(1.0f, 0.0f, 0.0f); // nonfeasible -> red
  1077 +
1062 1078 glBegin(GL_LINE_STRIP);
1063 1079 for (unsigned j = 0; j < inlet[i].V.size(); j++)
1064 1080 glVertex3f(inlet[i].V[j][0], inlet[i].V[j][1], inlet[i].V[j][2]);
1065 1081 glEnd();
1066 1082 }
1067 1083 for (unsigned i = 0; i < outlet.size(); i++) {
  1084 + if (outlet_feasibility[i])
  1085 + glColor3f(0.0f, 0.0f, 0.0f); // feasible -> black
  1086 + else
  1087 + glColor3f(1.0f, 0.0f, 0.0f); // nonfeasible -> red
1068 1088 glBegin(GL_LINE_STRIP);
1069 1089 for (unsigned j = 0; j < outlet[i].V.size(); j++)
1070 1090 glVertex3f(outlet[i].V[j][0], outlet[i].V[j][1], outlet[i].V[j][2]);
... ... @@ -1255,6 +1275,7 @@ namespace stim {
1255 1275 bool online = false; // flag indicates the point is on the line-segment
1256 1276 float a, b;
1257 1277  
  1278 + // inner network
1258 1279 for (unsigned i = 0; i < E.size(); i++) {
1259 1280 for (unsigned j = 0; j < E[i].size() - 1; j++) {
1260 1281 v1 = E[i][j + 1] - E[i][j]; // a -> b = ab
... ... @@ -1284,6 +1305,60 @@ namespace stim {
1284 1305 }
1285 1306 }
1286 1307  
  1308 + // inlet connection
  1309 + for (unsigned i = 0; i < inlet.size(); i++) {
  1310 + v1 = inlet[i].V[2] - inlet[i].V[1];
  1311 + v2 = v0 - inlet[i].V[1];
  1312 + v3 = v0 - inlet[i].V[2];
  1313 +
  1314 + tmp_d = v2.dot(v1); // av·ab
  1315 +
  1316 + // check the line relative position
  1317 + a = v2.dot(v1.norm());
  1318 + b = v3.dot(v1.norm());
  1319 + if (a < v1.len() && b < v1.len()) // if the length of projection fragment is longer than the line-segment
  1320 + online = true;
  1321 + else
  1322 + online = false;
  1323 +
  1324 + if (tmp_d <= 0.0 || tmp_d > std::pow(v1.len(), 2) && !online) // projection lies outside the line-segment
  1325 + continue;
  1326 + else {
  1327 + tmp_d = v1.cross(v2).len() / v1.len(); // perpendicular distance of point to segment: |v1 x v2| / |v1|
  1328 + if (tmp_d < d) {
  1329 + d = tmp_d;
  1330 + tmp_i = i;
  1331 + }
  1332 + }
  1333 + }
  1334 +
  1335 + // outlet connection
  1336 + for (unsigned i = 0; i < outlet.size(); i++) {
  1337 + v1 = outlet[i].V[2] - outlet[i].V[1];
  1338 + v2 = v0 - outlet[i].V[1];
  1339 + v3 = v0 - outlet[i].V[2];
  1340 +
  1341 + tmp_d = v2.dot(v1); // av·ab
  1342 +
  1343 + // check the line relative position
  1344 + a = v2.dot(v1.norm());
  1345 + b = v3.dot(v1.norm());
  1346 + if (a < v1.len() && b < v1.len()) // if the length of projection fragment is longer than the line-segment
  1347 + online = true;
  1348 + else
  1349 + online = false;
  1350 +
  1351 + if (tmp_d <= 0.0 || tmp_d > std::pow(v1.len(), 2) && !online) // projection lies outside the line-segment
  1352 + continue;
  1353 + else {
  1354 + tmp_d = v1.cross(v2).len() / v1.len(); // perpendicular distance of point to segment: |v1 x v2| / |v1|
  1355 + if (tmp_d < d) {
  1356 + d = tmp_d;
  1357 + tmp_i = i;
  1358 + }
  1359 + }
  1360 + }
  1361 +
1287 1362 eps += get_radius(tmp_i, tmp_j);
1288 1363  
1289 1364 if (d < eps) {
... ... @@ -1293,6 +1368,86 @@ namespace stim {
1293 1368  
1294 1369 return false;
1295 1370 }
  1371 + inline bool epsilon_edge(T x, T y, T z, T eps, unsigned &idx, unsigned &port) {
  1372 +
  1373 + T d = FLT_MAX;
  1374 + T tmp_d;
  1375 + unsigned tmp_i = 0;
  1376 + stim::vec3<T> v1;
  1377 + stim::vec3<T> v2;
  1378 + stim::vec3<T> v3;
  1379 + stim::vec3<T> v0 = stim::vec3<float>(x, y, z);
  1380 + bool online = false; // flag indicates the point is on the line-segment
  1381 + float a, b;
  1382 +
  1383 + // inlet connection
  1384 + for (unsigned i = 0; i < inlet.size(); i++) {
  1385 + for (unsigned j = 0; j < inlet[i].V.size() - 1; j++) {
  1386 + v1 = inlet[i].V[j + 1] - inlet[i].V[j];
  1387 + v2 = v0 - inlet[i].V[j];
  1388 + v3 = v0 - inlet[i].V[j + 1];
  1389 +
  1390 + tmp_d = v2.dot(v1); // av·ab
  1391 +
  1392 + // check the line relative position
  1393 + a = v2.dot(v1.norm());
  1394 + b = v3.dot(v1.norm());
  1395 + if (a < v1.len() && b < v1.len()) // if the length of projection fragment is longer than the line-segment
  1396 + online = true;
  1397 + else
  1398 + online = false;
  1399 +
  1400 + if (tmp_d <= 0.0 || tmp_d > std::pow(v1.len(), 2) && !online) // projection lies outside the line-segment
  1401 + continue;
  1402 + else {
  1403 + tmp_d = v1.cross(v2).len() / v1.len(); // perpendicular distance of point to segment: |v1 x v2| / |v1|
  1404 + if (tmp_d < d) {
  1405 + d = tmp_d;
  1406 + tmp_i = i;
  1407 + port = 0;
  1408 + }
  1409 + }
  1410 +
  1411 + }
  1412 + }
  1413 +
  1414 + // outlet connection
  1415 + for (unsigned i = 0; i < outlet.size(); i++) {
  1416 + for (unsigned j = 0; j < outlet[i].V.size() - 1; j++) {
  1417 + v1 = outlet[i].V[j + 1] - outlet[i].V[j];
  1418 + v2 = v0 - outlet[i].V[j];
  1419 + v3 = v0 - outlet[i].V[j + 1];
  1420 +
  1421 + tmp_d = v2.dot(v1); // av·ab
  1422 +
  1423 + // check the line relative position
  1424 + a = v2.dot(v1.norm());
  1425 + b = v3.dot(v1.norm());
  1426 + if (a < v1.len() && b < v1.len()) // if the length of projection fragment is longer than the line-segment
  1427 + online = true;
  1428 + else
  1429 + online = false;
  1430 +
  1431 + if (tmp_d <= 0.0 || tmp_d > std::pow(v1.len(), 2) && !online) // projection lies outside the line-segment
  1432 + continue;
  1433 + else {
  1434 + tmp_d = v1.cross(v2).len() / v1.len(); // perpendicular distance of point to segment: |v1 x v2| / |v1|
  1435 + if (tmp_d < d) {
  1436 + d = tmp_d;
  1437 + tmp_i = i;
  1438 + port = 1;
  1439 + }
  1440 + }
  1441 + }
  1442 + }
  1443 +
  1444 + if (d < eps) {
  1445 + idx = tmp_i;
  1446 + return true;
  1447 + }
  1448 +
  1449 + return false;
  1450 + }
1296 1451  
1297 1452 /// build main feeder connection
1298 1453 // set up main feeder and main port of both input and output
... ... @@ -1404,6 +1559,8 @@ namespace stim {
1404 1559  
1405 1560 outlet.push_back(tmp_b);
1406 1561 }
  1562 +
  1563 + backup();
1407 1564 }
1408 1565  
1409 1566 // find the number of U-shape or square-shape structure for extending length of connection
... ... @@ -1436,12 +1593,13 @@ namespace stim {
1436 1593 int coef_z = (z) ? 1 : -1; // z coefficient
1437 1594 int inverse = 1; // inverse flag
1438 1595 stim::vec3<T> cor_v; // corner vertex
1439   -
  1596 + std::pair<stim::vec3<T>, stim::vec3<T>> tmp_bb;
1440 1597 stim::vec3<T> tmp_v;
1441   - if (feeder == 1)
  1598 + if (feeder == 1)
1442 1599 tmp_v = inlet[i].V[0];
1443   - else if (feeder == 0)
  1600 + else if (feeder == 0)
1444 1601 tmp_v = outlet[i].V[0];
  1602 + tmp_bb.first = tmp_v;
1445 1603  
1446 1604 // check whether the height of connections is acceptable
1447 1605 if (height > threshold) { // acceptable
... ... @@ -1522,6 +1680,8 @@ namespace stim {
1522 1680 inlet[i].V.push_back(tmp_v);
1523 1681 else if (feeder == 0)
1524 1682 outlet[i].V.push_back(tmp_v);
  1683 + if (j == n - 1 && k == 0) // first time go "in"
  1684 + tmp_bb.second = tmp_v;
1525 1685  
1526 1686 tmp_v = tmp_v + stim::vec3<T>(0, inverse * coef_up * width, 0);
1527 1687 if (feeder == 1)
... ... @@ -1555,7 +1715,7 @@ namespace stim {
1555 1715 inverse = 1;
1556 1716 }
1557 1717 // if use fragment to do square wave connection, need to push_back the corner vertex
1558   - if (ratio > 0.0f && ratio <= 1.0f) {
  1718 + if (ratio > 0.0f && ratio < 1.0f) {
1559 1719 if (feeder == 1)
1560 1720 inlet[i].V.push_back(cor_v);
1561 1721 else if (feeder == 0)
... ... @@ -1565,25 +1725,30 @@ namespace stim {
1565 1725 else {
1566 1726 for (int j = 0; j < n; j++) {
1567 1727  
1568   - // move in Z-shape
  1728 + // up
1569 1729 tmp_v = tmp_v + stim::vec3<T>(0, coef_up * height, 0);
1570 1730 if (feeder == 1)
1571 1731 inlet[i].V.push_back(tmp_v);
1572 1732 else if (feeder == 0)
1573 1733 outlet[i].V.push_back(tmp_v);
1574 1734  
  1735 + // left
1575 1736 tmp_v = tmp_v + stim::vec3<T>(-coef_left * width, 0, 0);
1576 1737 if (feeder == 1)
1577 1738 inlet[i].V.push_back(tmp_v);
1578 1739 else if (feeder == 0)
1579 1740 outlet[i].V.push_back(tmp_v);
  1741 + if (j == n - 1)
  1742 + tmp_bb.second = tmp_v;
1580 1743  
  1744 + // down
1581 1745 tmp_v = tmp_v + stim::vec3<T>(0, -coef_up * height, 0);
1582 1746 if (feeder == 1)
1583 1747 inlet[i].V.push_back(tmp_v);
1584 1748 else if (feeder == 0)
1585 1749 outlet[i].V.push_back(tmp_v);
1586 1750  
  1751 + // left
1587 1752 tmp_v = tmp_v + stim::vec3<T>(-coef_left * width, 0, 0);
1588 1753 if (feeder == 1)
1589 1754 inlet[i].V.push_back(tmp_v);
... ... @@ -1591,6 +1756,10 @@ namespace stim {
1591 1756 outlet[i].V.push_back(tmp_v);
1592 1757 }
1593 1758 }
  1759 + if (feeder == 1)
  1760 + inbb[i] = tmp_bb;
  1761 + else if (feeder == 0)
  1762 + outbb[i] = tmp_bb;
1594 1763 }
1595 1764  
1596 1765 // automatically modify bridge to make it feasible
... ... @@ -1598,7 +1767,7 @@ namespace stim {
1598 1767  
1599 1768 glDeleteLists(dlist, 1); // delete display list for modify
1600 1769 glDeleteLists(dlist + 1, 1);
1601   -
  1770 +
1602 1771 // because of radius change at the port vertex, there will be a pressure drop at that port
1603 1772 // it follows the bernoulli equation
1604 1773 // p1 + 1/2*rou*v1^2 + rou*g*h1 = p2 + 1/2*rou*v2^2 + rou*g*h2
... ... @@ -1815,13 +1984,15 @@ namespace stim {
1815 1984 // automatically modify inlet bridge using square shape constructions
1816 1985 else {
1817 1986 bool upper; // flag indicates the connection is upper than the bus
1818   - bool z; // flag indicates the connection direction along z-axis
  1987 + bool z; // flag indicates the connection direction along z-axis
1819 1988 T new_l; // new length
1820 1989 stim::vec3<T> bus_v; // the port point on the bus
1821 1990 stim::vec3<T> mid_v; // the original corner point
1822 1991 stim::vec3<T> tmp_v; // the pendant point
1823 1992 int n;
1824 1993 T width, height; // width and height of the square
  1994 + inbb.resize(inlet.size()); // resize bounding box of inlets/outlets connections
  1995 + outbb.resize(outlet.size());
1825 1996  
1826 1997 for (unsigned i = 0; i < inlet.size(); i++) {
1827 1998 if (i != inlet_index) {
... ... @@ -1857,6 +2028,10 @@ namespace stim {
1857 2028  
1858 2029 std::reverse(inlet[i].V.begin(), inlet[i].V.end()); // from bus to pendant vertex
1859 2030 }
  2031 + else {
  2032 + inbb[i].first = inlet[i].V[2];
  2033 + inbb[i].second = inlet[i].V[1];
  2034 + }
1860 2035 }
1861 2036  
1862 2037 for (unsigned i = 0; i < outlet.size(); i++) {
... ... @@ -1893,35 +2068,159 @@ namespace stim {
1893 2068  
1894 2069 std::reverse(outlet[i].V.begin(), outlet[i].V.end()); // from bus to pendant vertex
1895 2070 }
  2071 + else {
  2072 + outbb[i].first = outlet[i].V[2];
  2073 + outbb[i].second = outlet[i].V[1];
  2074 + }
1896 2075 }
1897 2076 }
  2077 +
  2078 + check_special_connection(); // check special connections
1898 2079 }
1899 2080  
1900   - // check current bridge to see feasibility
1901   - void check_synthetic_connection(T viscosity, T radius = 5.0f) {
  2081 + /// check current connections to find overlapping
  2082 + // phase 1 check -> direct connection intersection
  2083 + void check_direct_connection() {
1902 2084  
1903   - T eps = 0.01f;
1904   - T source_pressure = pressure[inlet[0].v[0]] + (8 * viscosity * inlet[0].l * inlet[0].Q) / ((float)stim::PI * std::pow(radius, 4));
1905   - T tmp_p;
1906   - for (unsigned i = 1; i < inlet.size(); i++) {
1907   - tmp_p = pressure[inlet[i].v[0]] + (8 * viscosity * inlet[i].l * inlet[i].Q) / ((float)stim::PI * std::pow(radius, 4));
1908   - T delta = fabs(tmp_p - source_pressure);
1909   - if (delta > eps) {
1910   - std::cout << "Nonfeasible connection!" << std::endl;
1911   - break;
  2085 + unsigned num;
  2086 + // check inlet
  2087 + num = inlet.size();
  2088 + inlet_feasibility.resize(num, true);
  2089 + for (unsigned i = 0; i < num; i++) {
  2090 + for (unsigned j = 0; j < num; j++) {
  2091 + if (i != j) {
  2092 + if (inlet[i].V[0][2] == inlet[j].V[0][2]) {
  2093 + if ((inlet[i].V[1][0] <= inlet[j].V[1][0] && fabs(inlet[i].V[1][1]) <= fabs(inlet[j].V[1][1])) || (inlet[i].V[1][0] >= inlet[j].V[1][0] && fabs(inlet[i].V[1][1]) >= fabs(inlet[j].V[1][1]))) {
  2094 + inlet_feasibility[i] = false;
  2095 + break;
  2096 + }
  2097 + }
  2098 + else
  2099 + inlet_feasibility[i] = true;
  2100 + }
1912 2101 }
1913 2102 }
1914   - source_pressure = pressure[outlet[0].v[0]] - (8 * viscosity * outlet[0].l * outlet[0].Q) / ((float)stim::PI * std::pow(radius, 4));
1915   - for (unsigned i = 1; i < outlet.size(); i++) {
1916   - tmp_p = pressure[outlet[i].v[0]] - (8 * viscosity * outlet[i].l * outlet[i].Q) / ((float)stim::PI * std::pow(radius, 4));
1917   - T delta = fabs(tmp_p - source_pressure);
1918   - if (delta > eps) {
1919   - std::cout << "Nonfeasible connection!" << std::endl;
1920   - break;
  2103 +
  2104 + // check outlet
  2105 + num = outlet.size();
  2106 + outlet_feasibility.resize(num, true);
  2107 + for (unsigned i = 0; i < num; i++) {
  2108 + for (unsigned j = 0; j < num; j++) {
  2109 + if (i != j) {
  2110 + if (outlet[i].V[0][2] == outlet[j].V[0][2]) {
  2111 + if ((outlet[i].V[1][0] <= outlet[j].V[1][0] && fabs(outlet[i].V[1][1]) >= fabs(outlet[j].V[1][1])) || (outlet[i].V[1][0] >= outlet[j].V[1][0] && fabs(outlet[i].V[1][1]) <= fabs(outlet[j].V[1][1]))) {
  2112 + outlet_feasibility[i] = false;
  2113 + break;
  2114 + }
  2115 + }
  2116 + else
  2117 + outlet_feasibility[i] = true;
  2118 + }
  2119 + }
  2120 + }
  2121 + }
  2122 +
  2123 + // phase 2 check -> special connection intersection
  2124 + void check_special_connection(T radius = 5.0f) {
  2125 +
  2126 + // temp AABB centers and halfwidths
  2127 + stim::vec3<T> c1, c2;
  2128 + stim::vec3<T> r1, r2;
  2129 + // inlets' special connections checking
  2130 + for (unsigned i = 0; i < inbb.size(); i++) {
  2131 + for (unsigned j = 0; j < inbb.size(); j++) {
  2132 + if (j != i) {
  2133 + c1 = stim::vec3<T>((inbb[i].first + inbb[i].second) / 2);
  2134 + c2 = stim::vec3<T>((inbb[j].first + inbb[j].second) / 2);
  2135 + for (unsigned k = 0; k < 3; k++) {
  2136 + r1[k] = fabs(inbb[i].first[k] - inbb[i].second[k]) / 2;
  2137 + r2[k] = fabs(inbb[j].first[k] - inbb[j].second[k]) / 2;
  2138 + }
  2139 + // test AABBAABB
  2140 + if (fabs(c1[0] - c2[0]) > (r1[0] + r2[0] + 2 * radius) || fabs(c1[1] - c2[1]) > (r1[1] + r2[1] + 2 * radius) || fabs(c1[2] - c2[2]) > (r1[2] + r2[2] + 2 * radius))
  2141 + inlet_feasibility[i] = true;
  2142 + else
  2143 + inlet_feasibility[i] = false;
  2144 + }
  2145 + }
  2146 + }
  2147 +
  2148 + // outlets' special connections checking
  2149 + for (unsigned i = 0; i < outbb.size(); i++) {
  2150 + for (unsigned j = 0; j < outbb.size(); j++) {
  2151 + if (j != i) {
  2152 + c1 = stim::vec3<T>((outbb[i].first + outbb[i].second) / 2);
  2153 + c2 = stim::vec3<T>((outbb[j].first + outbb[j].second) / 2);
  2154 + for (unsigned k = 0; k < 3; k++) {
  2155 + r1[k] = fabs(outbb[i].first[k] - outbb[i].second[k]) / 2;
  2156 + r2[k] = fabs(outbb[j].first[k] - outbb[j].second[k]) / 2;
  2157 + }
  2158 + // test AABBAABB
  2159 + if (fabs(c1[0] - c2[0]) > (r1[0] + r2[0] + 2 * radius) || fabs(c1[1] - c2[1]) > (r1[1] + r2[1] + 2 * radius) || fabs(c1[2] - c2[2]) > (r1[2] + r2[2] + 2 * radius))
  2160 + outlet_feasibility[i] = true;
  2161 + else
  2162 + outlet_feasibility[i] = false;
  2163 + }
  2164 + }
  2165 + }
  2166 + }
  2167 +
  2168 + // clear synthetic connections
  2169 + void clear_synthetic_connection() {
  2170 +
  2171 + // restore direct synthetic connecions
  2172 + T l = 0.0f;
  2173 + for (unsigned i = 0; i < inlet.size(); i++) {
  2174 + inlet[i].V.clear();
  2175 + for (unsigned j = 0; j < in_backup[i].size(); j++) {
  2176 + inlet[i].V.push_back(in_backup[i][j]);
  2177 + if (j != in_backup[i].size() - 1)
  2178 + l += (in_backup[i][j + 1] - in_backup[i][j]).len();
  2179 + }
  2180 + inlet[i].l = l;
  2181 + l = 0.0f;
  2182 + }
  2183 + for (unsigned i = 0; i < outlet.size(); i++) {
  2184 + outlet[i].V.clear();
  2185 + for (unsigned j = 0; j < out_backup[i].size(); j++) {
  2186 + outlet[i].V.push_back(out_backup[i][j]);
  2187 + if (j != out_backup[i].size() - 1)
  2188 + l += (out_backup[i][j + 1] - out_backup[i][j]).len();
1921 2189 }
  2190 + outlet[i].l = l;
  2191 + l = 0.0f;
1922 2192 }
  2193 +
  2194 + // clear up inlets/outlets connection bounding box
  2195 + inbb.clear();
  2196 + outbb.clear();
1923 2197 }
1924 2198  
  2199 + // back up direct synthetic connection whenever modified
  2200 + void backup() {
  2201 +
  2202 + in_backup.clear();
  2203 + out_backup.clear();
  2204 +
  2205 + // back up direct synthetic connecions
  2206 + std::vector<typename stim::vec3<T> > V;
  2207 + for (unsigned i = 0; i < inlet.size(); i++) {
  2208 + for (unsigned j = 0; j < inlet[i].V.size(); j++) {
  2209 + V.push_back(inlet[i].V[j]);
  2210 + }
  2211 + in_backup.push_back(V);
  2212 + V.clear();
  2213 + }
  2214 + for (unsigned i = 0; i < outlet.size(); i++) {
  2215 + for (unsigned j = 0; j < outlet[i].V.size(); j++) {
  2216 + V.push_back(outlet[i].V[j]);
  2217 + }
  2218 + out_backup.push_back(V);
  2219 + V.clear();
  2220 + }
  2221 + }
  2222 +
  2223 +
1925 2224 /// make binary image stack
1926 2225 // prepare for image stack
1927 2226 void preparation(T &Xl, T &Xr, T &Yt, T &Yb, T &Z, T length = 210.0f, T height = 10.0f) {
... ...
... ... @@ -46,6 +46,7 @@ float min_v;
46 46 int mods; // special keyboard input
47 47 std::vector<unsigned char> color; // velocity color map
48 48 std::vector<int> velocity_bar; // velocity bar
  49 +float length = 210.0f; // cuboid length
49 50  
50 51 // hard-coded parameters
51 52 float camera_factor = 1.2f; // start point of the camera as a function of X and Y size
... ... @@ -64,7 +65,9 @@ float fragment_ratio = 0.0f; // fragment ratio
64 65 // glut event parameters
65 66 int mouse_x; // window x-coordinate
66 67 int mouse_y; // window y-coordinate
67   -bool LTbutton = false; // true means down while false means up
  68 +int picked_x; // picked window x-coordinate
  69 +int picked_y; // picked window y-coordinate
  70 +bool LTbutton = false; // true means down while false means up
68 71  
69 72 // simulation parameters
70 73 bool render_direction = false; // flag indicates rendering flow direction for one edge
... ... @@ -79,6 +82,14 @@ bool build_inlet_outlet = false; // flag indicates building inlets and outlets
79 82 bool modified_bridge = false; // flag indicates having modified inlet/outlet connection
80 83 bool hilbert_curve = false; // flag indicates enabling hilbert curves constructions
81 84 bool change_fragment = false; // flag indicates changing fragment for square wave connections
  85 +bool picked_connection = false; // flag indicates picked one connection
  86 +bool render_new_connection = false; // flag indicates rendering new line connection in trasparency
  87 +bool redisplay = false; // flag indicates redisplay rendering
  88 +bool connection_done = false; // flag indicates finishing connections
  89 +unsigned connection_index = -1; // the index of connection that is picked
  90 +unsigned port_index = 0; // inlet (0) or outlet (1)
  91 +stim::vec3<float> tmp_v1, tmp_v2; // temp vertex
  92 +int coef; // computational coefficient factor
82 93  
83 94 // manufacture parameters
84 95 bool manufacture = false; // flag indicates manufacture mode
... ... @@ -211,19 +222,41 @@ void glut_render() {
211 222 flow.mark_vertex();
212 223 //flow.glSolidCone(subdivision);
213 224 flow.glSolidCylinder(direction_index, color, subdivision);
214   - flow.glSolidCuboid();
  225 + flow.glSolidCuboid(length);
215 226 if (render_direction)
216 227 flow.glSolidCone(direction_index, subdivision);
217 228 }
218 229  
219 230 if (build_inlet_outlet) {
220   - flow.line_bridge();
  231 + flow.line_bridge(redisplay);
221 232 }
222 233  
223 234 if (manufacture) {
224 235 flow.glSolidCuboid();
225 236 flow.tube_bridge(subdivision);
226 237 }
  238 +
  239 + if (picked_connection && render_new_connection) {
  240 + glEnable(GL_BLEND);
  241 + glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
  242 + glColor4f(0.0f, 0.0f, 0.0f, 0.4f);
  243 + glBegin(GL_LINE_STRIP);
  244 + if (!port_index) {
  245 + glVertex3f(flow.inlet[connection_index].V[1][0], flow.inlet[connection_index].V[1][1], flow.inlet[connection_index].V[1][2]);
  246 + glVertex3f(tmp_v1[0], tmp_v1[1], tmp_v1[2]);
  247 + glVertex3f(tmp_v2[0], tmp_v2[1], tmp_v2[2]);
  248 + glVertex3f(flow.inlet[connection_index].V[2][0], flow.inlet[connection_index].V[2][1], flow.inlet[connection_index].V[2][2]);
  249 + }
  250 + else {
  251 + glVertex3f(flow.outlet[connection_index].V[1][0], flow.outlet[connection_index].V[1][1], flow.outlet[connection_index].V[1][2]);
  252 + glVertex3f(tmp_v1[0], tmp_v1[1], tmp_v1[2]);
  253 + glVertex3f(tmp_v2[0], tmp_v2[1], tmp_v2[2]);
  254 + glVertex3f(flow.outlet[connection_index].V[2][0], flow.outlet[connection_index].V[2][1], flow.outlet[connection_index].V[2][2]);
  255 + }
  256 + glEnd();
  257 + glFlush();
  258 + glDisable(GL_BLEND);
  259 + }
227 260  
228 261 // render bars
229 262 // bring up a pressure bar on left
... ... @@ -287,12 +320,13 @@ void glut_render() {
287 320 glLoadIdentity();
288 321  
289 322 float step = (vY - 3 * border_factor);
290   - step /= num_edge;
291   - for (unsigned i = 0; i < num_edge; i++) {
  323 + step /= BREWER_CTRL_PTS - 1;
  324 + for (unsigned i = 0; i < BREWER_CTRL_PTS - 1; i++) {
292 325 glLineWidth(border_factor);
293 326 glBegin(GL_LINES);
294   - glColor3f((float)color[velocity_bar[i] * 3 + 0] / 255, (float)color[velocity_bar[i] * 3 + 1] / 255, (float)color[velocity_bar[i] * 3 + 2] / 255);
  327 + glColor3f(BREWERCP[i * 4 + 0], BREWERCP[i * 4 + 1], BREWERCP[i * 4 + 2]);
295 328 glVertex2f(border_factor, border_factor + i * step);
  329 + glColor3f(BREWERCP[(i + 1) * 4 + 0], BREWERCP[(i + 1) * 4 + 1], BREWERCP[(i + 1) * 4 + 2]);
296 330 glVertex2f(border_factor, border_factor + (i + 1) * step);
297 331 glEnd();
298 332 }
... ... @@ -376,6 +410,7 @@ void glut_menu(int value) {
376 410 build_inlet_outlet = false;
377 411 manufacture = false;
378 412 modified_bridge = false;
  413 + connection_done = false;
379 414 if (!flow.set)
380 415 flow_initialize();
381 416 flow_stable_state(); // main function of solving the linear system
... ... @@ -388,9 +423,16 @@ void glut_menu(int value) {
388 423 simulation = false;
389 424 build_inlet_outlet = true;
390 425 manufacture = false;
391   - if (!modified_bridge) {
  426 + if (!modified_bridge && !connection_done) {
392 427 flow.set_main_feeder();
393 428 flow.build_synthetic_connection(u, default_radius);
  429 + flow.check_direct_connection(); // check whether direct connections intersect each other
  430 + connection_done = true;
  431 + }
  432 + else if (modified_bridge) {
  433 + modified_bridge = false;
  434 + redisplay = true;
  435 + flow.clear_synthetic_connection();
394 436 }
395 437  
396 438 glut_set_menu(num, 3);
... ... @@ -400,7 +442,6 @@ void glut_menu(int value) {
400 442 simulation = false;
401 443 build_inlet_outlet = false;
402 444 manufacture = true;
403   - flow.check_synthetic_connection(u, default_radius);
404 445 }
405 446  
406 447 glutPostRedisplay();
... ... @@ -426,40 +467,99 @@ void glut_motion(int x, int y) {
426 467 // defines passive mouse motion function
427 468 void glut_passive_motion(int x, int y) {
428 469  
429   - mouse_x = x;
430   - mouse_y = y;
  470 + mods = glutGetModifiers();
431 471  
432 472 // check whether the mouse point near to an edge
433 473 GLdouble posX, posY, posZ;
434 474 window_to_world(posX, posY, posZ); // get the world coordinates
435 475  
436   - if (simulation || build_inlet_outlet) {
  476 + if (simulation || build_inlet_outlet && !mods) {
437 477 bool flag = flow.epsilon_edge((float)posX, (float)posY, (float)posZ, eps, direction_index);
438 478 if (flag)
439 479 render_direction = true;
440 480 else {
441 481 if (render_direction) // if the direction is displaying currently, do a short delay
442   - Sleep(1000);
  482 + Sleep(300);
443 483 render_direction = false;
444 484 direction_index = -1;
445 485 }
446 486 }
447 487  
448   - glutPostRedisplay(); // re-draw the visualization
  488 + if (mods == GLUT_ACTIVE_SHIFT && picked_connection) {
  489 + render_new_connection = true;
  490 + unsigned i;
  491 + if (!port_index) {
  492 + tmp_v1 = stim::vec3<float>(flow.inlet[connection_index].V[1][0], flow.inlet[connection_index].V[1][1] + (float)(picked_y - y), flow.inlet[connection_index].V[1][2]);
  493 + tmp_v2 = stim::vec3<float>(flow.inlet[connection_index].V[2][0], flow.inlet[connection_index].V[2][1] + (float)(picked_y - y), flow.inlet[connection_index].V[2][2]);
  494 + i = flow.inlet[connection_index].V.size();
  495 + if (coef * tmp_v1[1] < coef * flow.inlet[connection_index].V[i - 1][1]) {
  496 + tmp_v1[1] = flow.inlet[connection_index].V[i - 1][1];
  497 + tmp_v2[1] = flow.inlet[connection_index].V[i - 1][1];
  498 + }
  499 + }
  500 + else {
  501 + tmp_v1 = stim::vec3<float>(flow.outlet[connection_index].V[1][0], flow.outlet[connection_index].V[1][1] + (float)(picked_y - y), flow.outlet[connection_index].V[1][2]);
  502 + tmp_v2 = stim::vec3<float>(flow.outlet[connection_index].V[2][0], flow.outlet[connection_index].V[2][1] + (float)(picked_y - y), flow.outlet[connection_index].V[2][2]);
  503 + i = flow.outlet[connection_index].V.size();
  504 + if (coef * tmp_v1[1] < coef * flow.outlet[connection_index].V[i - 1][1]) {
  505 + tmp_v1[1] = flow.outlet[connection_index].V[i - 1][1];
  506 + tmp_v2[1] = flow.outlet[connection_index].V[i - 1][1];
  507 + }
  508 + }
  509 + }
  510 + else if (mods == GLUT_ACTIVE_CTRL && picked_connection) {
  511 + render_new_connection = true;
  512 + if (!port_index) {
  513 + tmp_v1 = stim::vec3<float>(flow.inlet[connection_index].V[0][0] + (float)(x - picked_x), flow.inlet[connection_index].V[0][1], flow.inlet[connection_index].V[0][2]);
  514 + tmp_v2 = stim::vec3<float>(flow.inlet[connection_index].V[1][0] + (float)(x - picked_x), flow.inlet[connection_index].V[1][1], flow.inlet[connection_index].V[1][2]);
  515 + if (tmp_v1[0] < flow.main_feeder[port_index][0] - length / 2) {
  516 + tmp_v1[0] = flow.main_feeder[port_index][0] - length / 2;
  517 + tmp_v2[0] = flow.main_feeder[port_index][0] - length / 2;
  518 + }
  519 + else if (tmp_v1[0] > flow.main_feeder[port_index][0] + length / 2) {
  520 + tmp_v1[0] = flow.main_feeder[port_index][0] + length / 2;
  521 + tmp_v2[0] = flow.main_feeder[port_index][0] + length / 2;
  522 + }
  523 + }
  524 + else {
  525 + tmp_v1 = stim::vec3<float>(flow.outlet[connection_index].V[0][0] + (float)(x - picked_x), flow.outlet[connection_index].V[0][1], flow.outlet[connection_index].V[0][2]);
  526 + tmp_v2 = stim::vec3<float>(flow.outlet[connection_index].V[1][0] + (float)(x - picked_x), flow.outlet[connection_index].V[1][1], flow.outlet[connection_index].V[1][2]);
  527 + if (tmp_v1[0] > flow.main_feeder[port_index][0] + length / 2) {
  528 + tmp_v1[0] = flow.main_feeder[port_index][0] + length / 2;
  529 + tmp_v2[0] = flow.main_feeder[port_index][0] + length / 2;
  530 + }
  531 + else if (tmp_v1[0] < flow.main_feeder[port_index][0] - length / 2) {
  532 + tmp_v1[0] = flow.main_feeder[port_index][0] - length / 2;
  533 + tmp_v2[0] = flow.main_feeder[port_index][0] - length / 2;
  534 + }
  535 + }
  536 + }
  537 + else
  538 + render_new_connection = false;
  539 +
  540 + mouse_x = x;
  541 + mouse_y = y;
  542 +
  543 + glutPostRedisplay(); // re-draw the visualization
449 544 }
450 545  
451 546 // get click window coordinates
452 547 void glut_mouse(int button, int state, int x, int y) {
453 548  
  549 + mods = glutGetModifiers(); // get special keyboard input
  550 +
454 551 mouse_x = x;
455 552 mouse_y = y;
  553 + if (!mods) {
  554 + picked_connection = false;
  555 + render_new_connection = false;
  556 + }
456 557 if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN)
457 558 LTbutton = true;
458 559 else if (button == GLUT_LEFT_BUTTON && state == GLUT_UP)
459 560 LTbutton = false;
460 561  
461   - if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN && simulation && !to_select_pressure) {
462   -
  562 + if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN && !mods && simulation && !to_select_pressure) {
463 563 GLdouble posX, posY, posZ;
464 564 window_to_world(posX, posY, posZ); // get the world coordinates
465 565  
... ... @@ -470,7 +570,7 @@ void glut_mouse(int button, int state, int x, int y) {
470 570 to_select_pressure = true;
471 571 }
472 572 }
473   - else if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN && simulation && to_select_pressure) {
  573 + else if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN && !mods && simulation && to_select_pressure) {
474 574 if (y > 2 * border_factor || y < vY - border_factor) { // within the pressure bar range
475 575 to_select_pressure = false;
476 576 float tmp_pressure = (float)(vY - y - border_factor) / ((float)vY - border_factor) * max_pressure;
... ... @@ -480,13 +580,141 @@ void glut_mouse(int button, int state, int x, int y) {
480 580 flow.print_flow();
481 581 }
482 582 }
483   - else if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN && modified_bridge && change_fragment) {
  583 + else if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN && !mods && modified_bridge && change_fragment) {
484 584 if (y > 2 * border_factor || y < vY - border_factor) { // within the ratio bar range
485 585 fragment_ratio = (float)(vY - y - border_factor) / ((float)vY - border_factor) * 1.0f;
486 586 flow.modify_synthetic_connection(u, rou, hilbert_curve, height_threshold, fragment_ratio, default_radius);
487 587 change_fragment = false;
488 588 }
489 589 }
  590 + // move connections along y-axis
  591 + else if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN && mods == GLUT_ACTIVE_SHIFT && !modified_bridge && !picked_connection) {
  592 + GLdouble posX, posY, posZ;
  593 + window_to_world(posX, posY, posZ); // get the world coordinates
  594 +
  595 + bool flag = flow.epsilon_edge((float)posX, (float)posY, (float)posZ, eps, connection_index, port_index);
  596 + if (flag) {
  597 + picked_connection = true;
  598 + picked_x = x;
  599 + picked_y = y;
  600 + if (!port_index)
  601 + if (flow.inlet[connection_index].V[2][1] > flow.main_feeder[port_index][1])
  602 + coef = 1;
  603 + else
  604 + coef = -1;
  605 + else
  606 + if (flow.outlet[connection_index].V[2][1] > flow.main_feeder[port_index][1])
  607 + coef = 1;
  608 + else
  609 + coef = -1;
  610 + }
  611 + else
  612 + picked_connection = false;
  613 + }
  614 + else if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN && mods == GLUT_ACTIVE_SHIFT && !modified_bridge && render_new_connection) {
  615 + float l = 0.0f;
  616 + std::vector<typename stim::vec3<float> > V;
  617 + unsigned i;
  618 + if (!port_index) {
  619 + i = flow.inlet[connection_index].V.size();
  620 + if (tmp_v2[1] != flow.inlet[connection_index].V[i - 1][1]) {
  621 + V.resize(4);
  622 + V[0] = flow.inlet[connection_index].V[0];
  623 + V[1] = tmp_v1;
  624 + V[2] = tmp_v2;
  625 + V[3] = flow.inlet[connection_index].V[i - 1];
  626 + std::swap(flow.inlet[connection_index].V, V);
  627 + }
  628 + else {
  629 + V.resize(3);
  630 + V[0] = flow.inlet[connection_index].V[0];
  631 + V[1] = tmp_v1;
  632 + V[2] = tmp_v2;
  633 + std::swap(flow.inlet[connection_index].V, V);
  634 + }
  635 + // calculate new length
  636 + for (unsigned i = 0; i < flow.inlet[connection_index].V.size() - 1; i++) {
  637 + l += (flow.inlet[connection_index].V[i + 1] - flow.inlet[connection_index].V[i]).len();
  638 + }
  639 + flow.inlet[connection_index].l = l;
  640 + }
  641 + else {
  642 + i = flow.outlet[connection_index].V.size();
  643 + if (tmp_v2[1] != flow.outlet[connection_index].V[i - 1][1]) {
  644 + V.resize(4);
  645 + V[0] = flow.outlet[connection_index].V[0];
  646 + V[1] = tmp_v1;
  647 + V[2] = tmp_v2;
  648 + V[3] = flow.outlet[connection_index].V[i - 1];
  649 + std::swap(flow.outlet[connection_index].V, V);
  650 + }
  651 + else {
  652 + V.resize(3);
  653 + V[0] = flow.outlet[connection_index].V[0];
  654 + V[1] = tmp_v1;
  655 + V[2] = tmp_v2;
  656 + std::swap(flow.outlet[connection_index].V, V);
  657 + }
  658 + // calculate new length
  659 + for (unsigned i = 0; i < flow.outlet[connection_index].V.size() - 1; i++) {
  660 + l += (flow.outlet[connection_index].V[i + 1] - flow.outlet[connection_index].V[i]).len();
  661 + }
  662 + flow.outlet[connection_index].l = l;
  663 + }
  664 +
  665 + redisplay = true;
  666 + render_new_connection = false;
  667 + picked_connection = false;
  668 +
  669 + flow.check_direct_connection();
  670 + flow.backup(); // back up direct synthetic connections
  671 + }
  672 + // move connections along x-axis
  673 + else if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN && mods == GLUT_ACTIVE_CTRL && !modified_bridge && !picked_connection) {
  674 + GLdouble posX, posY, posZ;
  675 + window_to_world(posX, posY, posZ); // get the world coordinates
  676 +
  677 + bool flag = flow.epsilon_edge((float)posX, (float)posY, (float)posZ, eps, connection_index, port_index);
  678 + if (flag) {
  679 + picked_connection = true;
  680 + picked_x = x;
  681 + picked_y = y;
  682 + if (!port_index)
  683 + coef = 1;
  684 + else
  685 + coef = -1;
  686 + }
  687 + else
  688 + picked_connection = false;
  689 + }
  690 + else if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN && mods == GLUT_ACTIVE_CTRL && !modified_bridge && render_new_connection) {
  691 + float l = 0.0f;
  692 + if (!port_index) {
  693 + flow.inlet[connection_index].V[0] = tmp_v1;
  694 + flow.inlet[connection_index].V[1] = tmp_v2;
  695 + // calculate new length
  696 + for (unsigned i = 0; i < flow.inlet[connection_index].V.size() - 1; i++) {
  697 + l += (flow.inlet[connection_index].V[i + 1] - flow.inlet[connection_index].V[i]).len();
  698 + }
  699 + flow.inlet[connection_index].l = l;
  700 + }
  701 + else {
  702 + flow.outlet[connection_index].V[0] = tmp_v1;
  703 + flow.outlet[connection_index].V[1] = tmp_v2;
  704 + // calculate new length
  705 + for (unsigned i = 0; i < flow.outlet[connection_index].V.size() - 1; i++) {
  706 + l += (flow.outlet[connection_index].V[i + 1] - flow.outlet[connection_index].V[i]).len();
  707 + }
  708 + flow.outlet[connection_index].l = l;
  709 + }
  710 +
  711 + redisplay = true;
  712 + render_new_connection = false;
  713 + picked_connection = false;
  714 +
  715 + flow.check_direct_connection();
  716 + flow.backup();
  717 + }
490 718 }
491 719  
492 720 // define camera move based on mouse wheel move
... ...