Jiaming Guo / flow3

Commit b61addd86534d3a979f98d8c777a6077a6359066

Authored by Jiaming Guo
1 parent b94dd1e0

add adjustment feature. fix utilzing aabb for testing overlapping bugs. add save function.

Showing 2 changed files with 356 additions and 98 deletions   Show diff stats
Inline Side-by-side
... ... @@ -57,6 +57,12 @@ namespace stim {
57 57 T h; // height
58 58 };
59 59  
  60 + template <typename T>
  61 + struct circuit {
  62 + std::vector<typename std::pair<unsigned, unsigned> > v; // end vertex index
  63 + std::vector<T> r; // branch resistence
  64 + };
  65 +
60 66 /// indicator function
61 67 #ifdef __CUDACC__
62 68 // for sphere
... ... @@ -240,6 +246,7 @@ namespace stim {
240 246 std::vector<T> pressure; // final pressure
241 247 std::vector<typename std::vector<typename stim::vec3<T> > > in_backup; // inlet connection back up
242 248 std::vector<typename std::vector<typename stim::vec3<T> > > out_backup;
  249 + std::string units; // length units
243 250  
244 251 public:
245 252  
... ... @@ -260,13 +267,15 @@ namespace stim {
260 267 std::vector<bool> outlet_feasibility;
261 268 std::vector<typename std::pair<stim::vec3<T>, stim::vec3<T> > > inbb; // inlet connection bounding box
262 269 std::vector<typename std::pair<stim::vec3<T>, stim::vec3<T> > > outbb; // outlet connection bounding box
  270 + T Ps; // source and end pressure
  271 + T Pe;
263 272  
264 273 flow() {} // default constructor
265 274 ~flow() {
266 275 for (unsigned i = 0; i < num_vertex; i++)
267 276 delete[] C[i];
268 277 delete[] C;
269   - }
  278 + } // default destructor
270 279  
271 280 void init(unsigned n_e, unsigned n_v) {
272 281  
... ... @@ -307,6 +316,10 @@ namespace stim {
307 316 }
308 317 }
309 318  
  319 + void set_units(std::string u) {
  320 + units = u;
  321 + }
  322 +
310 323 // copy radius from cylinder to flow
311 324 void set_radius(unsigned i, T radius) {
312 325  
... ... @@ -384,7 +397,7 @@ namespace stim {
384 397 start_vertex = get_start_vertex(i); // get the start vertex index of current edge
385 398 end_vertex = get_end_vertex(i); // get the end vertex index of current edge
386 399  
387   - C[start_vertex][end_vertex] = -((float)stim::PI * std::pow(get_average_r(i), 4)) / (8 * u * get_l(i));
  400 + C[start_vertex][end_vertex] = -((T)stim::PI * std::pow(get_average_r(i), 4)) / (8 * u * get_l(i));
388 401  
389 402 C[end_vertex][start_vertex] = C[start_vertex][end_vertex];
390 403 }
... ... @@ -427,9 +440,9 @@ namespace stim {
427 440 }
428 441  
429 442 // get the flow state from known pressure
430   - float start_pressure = 0.0;
431   - float end_pressure = 0.0;
432   - float deltaP = 0.0;
  443 + T start_pressure = 0.0;
  444 + T end_pressure = 0.0;
  445 + T deltaP = 0.0;
433 446 for (unsigned i = 0; i < num_edge; i++) {
434 447 start_vertex = get_start_vertex(i);
435 448 end_vertex = get_end_vertex(i);
... ... @@ -440,8 +453,8 @@ namespace stim {
440 453 Q[i].first = start_vertex;
441 454 Q[i].second = end_vertex;
442 455  
443   - Q[i].third = ((float)stim::PI * std::pow(get_average_r(i), 4) * deltaP) / (8 * u * get_l(i));
444   - v[i] = Q[i].third / ((float)stim::PI * std::pow(get_average_r(i), 2));
  456 + Q[i].third = ((T)stim::PI * std::pow(get_average_r(i), 4) * deltaP) / (8 * u * get_l(i));
  457 + v[i] = Q[i].third / ((T)stim::PI * std::pow(get_average_r(i), 2));
445 458 }
446 459 }
447 460  
... ... @@ -474,12 +487,12 @@ namespace stim {
474 487 void print_flow() {
475 488  
476 489 // show the pressure information in console box
477   - std::cout << "PRESSURE(g/um/s^2):" << std::endl;
  490 + std::cout << "PRESSURE(g/" << units << "/s^2):" << std::endl;
478 491 for (unsigned i = 0; i < num_vertex; i++) {
479 492 std::cout << "[" << i << "] " << pressure[i] << std::endl;
480 493 }
481 494 // show the flow rate information in console box
482   - std::cout << "VOLUME FLOW RATE(um^3/s):" << std::endl;
  495 + std::cout << "VOLUME FLOW RATE(" << units << "^3/s):" << std::endl;
483 496 for (unsigned i = 0; i < num_edge; i++) {
484 497 std::cout << "(" << Q[i].first << "," << Q[i].second << ")" << Q[i].third << std::endl;
485 498 }
... ... @@ -1000,14 +1013,17 @@ namespace stim {
1000 1013 }
1001 1014  
1002 1015 // draw main feeder as solid cube
1003   - void glSolidCuboid(T length = 40.0f, T height = 10.0f) {
  1016 + void glSolidCuboid(bool manufacture = false, T length = 40.0f, T height = 10.0f) {
1004 1017  
1005 1018 T width;
1006 1019 stim::vec3<T> L = bb.A; // get the bottom left corner
1007 1020 stim::vec3<T> U = bb.B; // get the top right corner
1008 1021 width = U[2] - L[2] + 10.0f;
1009 1022  
1010   - glColor3f(0.0f, 0.0f, 0.0f);
  1023 + if (manufacture)
  1024 + glColor3f(0.0f, 0.0f, 0.0f); // black color
  1025 + else
  1026 + glColor3f(0.5f, 0.5f, 0.5f); // gray color
1011 1027 for (unsigned i = 0; i < main_feeder.size(); i++) {
1012 1028 // front face
1013 1029 glBegin(GL_QUADS);
... ... @@ -1060,6 +1076,35 @@ namespace stim {
1060 1076 glFlush();
1061 1077 }
1062 1078  
  1079 + // display the total volume flow rate
  1080 + void display_flow_rate(T in, T out) {
  1081 +
  1082 + glMatrixMode(GL_PROJECTION); // set up the 2d viewport for mode text printing
  1083 + glPushMatrix();
  1084 + glLoadIdentity();
  1085 + int X = glutGet(GLUT_WINDOW_WIDTH); // get the current window width
  1086 + int Y = glutGet(GLUT_WINDOW_HEIGHT); // get the current window height
  1087 + glViewport(0, 0, X, Y); // locate to left bottom corner
  1088 + gluOrtho2D(0, X, 0, Y); // define othogonal aspect
  1089 + glColor3f(0.8f, 0.0f, 0.0f); // using red to show mode
  1090 + glMatrixMode(GL_MODELVIEW);
  1091 + glPushMatrix();
  1092 + glLoadIdentity();
  1093 +
  1094 + glRasterPos2f(X / 2, 5); // hard coded position!!!!!
  1095 + std::stringstream ss_p;
  1096 + ss_p << "Q = "; // Q = * um^3/s
  1097 + ss_p << in;
  1098 + ss_p << " ";
  1099 + ss_p << units;
  1100 + ss_p << "^3/s";
  1101 + glutBitmapString(GLUT_BITMAP_TIMES_ROMAN_24, (const unsigned char*)(ss_p.str().c_str()));
  1102 +
  1103 + glPopMatrix();
  1104 + glMatrixMode(GL_PROJECTION);
  1105 + glPopMatrix();
  1106 + }
  1107 +
1063 1108 // draw the bridge as lines
1064 1109 void line_bridge(bool &redisplay) {
1065 1110  
... ... @@ -1409,7 +1454,7 @@ namespace stim {
1409 1454 inlet_main_feeder = stim::vec3<T>(bb.A[0] - border, bb.center()[1], bb.center()[2]);
1410 1455 outlet_main_feeder = stim::vec3<T>(bb.B[0] + border, bb.center()[1], bb.center()[2]);
1411 1456  
1412   - main_feeder.push_back(inlet_main_feeder);
  1457 + main_feeder.push_back(inlet_main_feeder); // 0->inlet, 1->outlet
1413 1458 main_feeder.push_back(outlet_main_feeder);
1414 1459  
1415 1460 // find both input and output vertex
... ... @@ -1550,22 +1595,42 @@ namespace stim {
1550 1595 tmp_v = outlet[i].V[outlet[i].V.size() - 1];
1551 1596 tmp_bb.first = tmp_v;
1552 1597  
1553   - // check whether the height of connections is acceptable
1554   - if (height > threshold) { // acceptable
1555   - // re-calculate height
1556   - if (ratio > 0.0f && ratio <= 1.0f) { // use fragment if set
1557   - cor_v = tmp_v + stim::vec3<T>(-coef_left * origin_l, 0, 0); // get the original corner vertex
1558   - desire_l = desire_l - origin_l * (1.0f - ratio / 1.0f);
1559   - origin_l = (T)origin_l * ratio / 1.0f;
  1598 + // pre-set fragments
  1599 + if (ratio) {
  1600 + T tmp_d, tmp_l; // back ups
  1601 + tmp_d = desire_l;
  1602 + tmp_l = origin_l;
  1603 +
  1604 + cor_v = tmp_v + stim::vec3<T>(-coef_left * origin_l, 0, 0); // get the original corner vertex
  1605 + desire_l = desire_l - origin_l * (1.0f - ratio / 1.0f);
  1606 + origin_l = (T)origin_l * ratio / 1.0f;
  1607 + n = find_number_square(origin_l, desire_l);
  1608 +
  1609 + width = (T)origin_l / (2 * n); // updates
  1610 + height = (desire_l - origin_l) / (2 * n);
  1611 +
  1612 + if (width < times * radius) { // check feasibility
  1613 + ratio = 0.0f; // load
  1614 + desire_l = tmp_d;
  1615 + origin_l = tmp_l;
  1616 +
  1617 + std::cout << "Warning: current ratio is not feasible, use full original line." << std::endl;
1560 1618 n = find_number_square(origin_l, desire_l);
  1619 +
  1620 + width = (T)origin_l / (2 * n); // updates
  1621 + height = (desire_l - origin_l) / (2 * n);
1561 1622 }
1562   - height = (desire_l - (1 + 2 * n) * origin_l) / std::pow(2 * n, 2);
1563   - // check whether the connections are good
1564   - while (height > threshold) {
  1623 + }
  1624 +
  1625 + // check whether it needs 3D square-wave-like connections
  1626 + if (height > threshold) { // enbale 3D connections
  1627 +
  1628 + height = (desire_l - (1 + 2 * n) * origin_l) / std::pow(2 * n, 2); // compute new height in 3D structure
  1629 + while (height > threshold) { // increase order to decrease height
1565 1630 n++;
1566 1631 width = (T)(origin_l) / (2 * n);
1567 1632 height = (desire_l - (1 + 2 * n) * origin_l) / std::pow(2 * n, 2);
1568   - // check whether it appears overlap
  1633 + // check whether it appears overlap, if it appears we choose last time height even if it is larger than threshold
1569 1634 if (width < times * radius) {
1570 1635 n--;
1571 1636 width = (T)(origin_l) / (2 * n);
... ... @@ -1573,12 +1638,21 @@ namespace stim {
1573 1638 break;
1574 1639 }
1575 1640 }
  1641 +
  1642 + // check overlap in terms of height, has potential risk when both height and width are less than times * radius.
1576 1643 while (height < times * radius) {
1577 1644 n--;
1578 1645 width = (T)(origin_l) / (2 * n);
1579 1646 height = (desire_l - (1 + 2 * n) * origin_l) / std::pow(2 * n, 2);
1580 1647 }
1581 1648  
  1649 + // degenerated case
  1650 + if (n == 0) {
  1651 + n = 1;
  1652 + width = (T)(origin_l) / (2 * n);
  1653 + height = (desire_l - (1 + 2 * n) * origin_l) / std::pow(2 * n, 2);
  1654 + }
  1655 +
1582 1656 // cube-like structure construction
1583 1657 for (int j = 0; j < n; j++) {
1584 1658 // "up"
... ... @@ -1601,7 +1675,7 @@ namespace stim {
1601 1675 inlet[i].V.push_back(tmp_v);
1602 1676 else if (feeder == 0)
1603 1677 outlet[i].V.push_back(tmp_v);
1604   - // "down"
  1678 + // "up"
1605 1679 tmp_v = tmp_v + stim::vec3<T>(0, inverse * coef_up * width, 0);
1606 1680 if (feeder == 1)
1607 1681 inlet[i].V.push_back(tmp_v);
... ... @@ -1623,27 +1697,28 @@ namespace stim {
1623 1697  
1624 1698 // "down"
1625 1699 for (int k = 0; k < n; k++) {
1626   -
  1700 + // in
1627 1701 tmp_v = tmp_v + stim::vec3<T>(0, 0, coef_z * height);
1628 1702 if (feeder == 1)
1629 1703 inlet[i].V.push_back(tmp_v);
1630 1704 else if (feeder == 0)
1631 1705 outlet[i].V.push_back(tmp_v);
  1706 + // get the bounding box edge
1632 1707 if (j == n - 1 && k == 0) // first time go "in"
1633 1708 tmp_bb.second = tmp_v;
1634   -
  1709 + // "down"
1635 1710 tmp_v = tmp_v + stim::vec3<T>(0, inverse * coef_up * width, 0);
1636 1711 if (feeder == 1)
1637 1712 inlet[i].V.push_back(tmp_v);
1638 1713 else if (feeder == 0)
1639 1714 outlet[i].V.push_back(tmp_v);
1640   -
  1715 + // out
1641 1716 tmp_v = tmp_v + stim::vec3<T>(0, 0, -coef_z * height);
1642 1717 if (feeder == 1)
1643 1718 inlet[i].V.push_back(tmp_v);
1644 1719 else if (feeder == 0)
1645 1720 outlet[i].V.push_back(tmp_v);
1646   -
  1721 + // "down"
1647 1722 tmp_v = tmp_v + stim::vec3<T>(0, inverse * coef_up * width, 0);
1648 1723 if (feeder == 1)
1649 1724 inlet[i].V.push_back(tmp_v);
... ... @@ -1671,12 +1746,13 @@ namespace stim {
1671 1746 outlet[i].V.push_back(cor_v);
1672 1747 }
1673 1748 }
  1749 + // use 2D square-wave-like connections
1674 1750 else {
1675 1751 if (height < times * radius) { // if height is too small, decrease n and re-calculate height and width
1676 1752 height = times * radius;
1677 1753 T need_l = desire_l - origin_l;
1678 1754 n = need_l / (2 * height);
1679   - if (n == 0)
  1755 + if (n == 0) // degenerated case
1680 1756 n = 1;
1681 1757 height = need_l / (2 * n);
1682 1758 width = origin_l / (2 * n);
... ... @@ -1713,6 +1789,13 @@ namespace stim {
1713 1789 else if (feeder == 0)
1714 1790 outlet[i].V.push_back(tmp_v);
1715 1791 }
  1792 + // if use fragment to do square wave connection, need to push_back the corner vertex
  1793 + if (ratio > 0.0f && ratio < 1.0f) {
  1794 + if (feeder == 1)
  1795 + inlet[i].V.push_back(cor_v);
  1796 + else if (feeder == 0)
  1797 + outlet[i].V.push_back(cor_v);
  1798 + }
1716 1799 }
1717 1800 if (feeder == 1)
1718 1801 inbb[i] = tmp_bb;
... ... @@ -1721,7 +1804,7 @@ namespace stim {
1721 1804 }
1722 1805  
1723 1806 // automatically modify bridge to make it feasible
1724   - void modify_synthetic_connection(T viscosity, T rou, bool H, T threshold, T ratio = 0.0f, T radius = 5.0f) {
  1807 + void modify_synthetic_connection(T viscosity, T rou, bool H, T threshold, T &in, T &out, T ratio = 0.0f, T radius = 5.0f) {
1725 1808  
1726 1809 glDeleteLists(dlist, 1); // delete display list for modify
1727 1810 glDeleteLists(dlist + 1, 1);
... ... @@ -1751,6 +1834,7 @@ namespace stim {
1751 1834 inlet_index = i;
1752 1835 }
1753 1836 }
  1837 + Ps = source_pressure;
1754 1838  
1755 1839 // find minimum pressure outlet port
1756 1840 T end_pressure = FLT_MAX;
... ... @@ -1762,6 +1846,7 @@ namespace stim {
1762 1846 outlet_index = i;
1763 1847 }
1764 1848 }
  1849 + Pe = end_pressure;
1765 1850  
1766 1851 // automatically modify inlet bridge using Hilbert curves
1767 1852 if (H) {
... ... @@ -1990,7 +2075,7 @@ namespace stim {
1990 2075 width = (T)origin_l / (2 * n);
1991 2076 height = (desire_l - origin_l) / (2 * n);
1992 2077  
1993   - build_square_connection(i, width, height, origin_l, desire_l, n, 1, threshold, z, true, upper, 10, ratio);
  2078 + build_square_connection(i, width, height, origin_l, desire_l, n, 1, threshold, z, true, upper, 5, ratio);
1994 2079 inlet[i].V.push_back(bus_v);
1995 2080  
1996 2081 std::reverse(inlet[i].V.begin(), inlet[i].V.end()); // from bus to pendant vertex
... ... @@ -2039,7 +2124,7 @@ namespace stim {
2039 2124 width = (T)origin_l / (2 * n);
2040 2125 height = (desire_l - origin_l) / (2 * n);
2041 2126  
2042   - build_square_connection(i, width, height, origin_l, desire_l, n, 0, threshold, z, false, upper, 10, ratio);
  2127 + build_square_connection(i, width, height, origin_l, desire_l, n, 0, threshold, z, false, upper, 5, ratio);
2043 2128 outlet[i].V.push_back(bus_v);
2044 2129  
2045 2130 std::reverse(outlet[i].V.begin(), outlet[i].V.end()); // from bus to pendant vertex
... ... @@ -2051,6 +2136,13 @@ namespace stim {
2051 2136 }
2052 2137 }
2053 2138  
  2139 + // save in-/out- volume flow rate
  2140 + in = out = 0.0f;
  2141 + for (unsigned i = 0; i < inlet.size(); i++)
  2142 + in += inlet[i].Q;
  2143 + for (unsigned i = 0; i < outlet.size(); i++)
  2144 + out += outlet[i].Q;
  2145 +
2054 2146 check_special_connection(); // check special connections
2055 2147 }
2056 2148  
... ... @@ -2060,23 +2152,22 @@ namespace stim {
2060 2152  
2061 2153 unsigned num;
2062 2154 // check inlet
2063   - num = inlet.size();
2064   - inlet_feasibility.resize(num, true);
  2155 + num = inlet.size(); // get the number of inlets
  2156 + inlet_feasibility.resize(num, true); // initialization
2065 2157 for (unsigned i = 0; i < num; i++) {
2066 2158 for (unsigned j = 0; j < num; j++) {
2067 2159 if (i != j) {
2068   - if (inlet[i].V[0][2] == inlet[j].V[0][2]) {
2069   - 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]))) {
  2160 + if (inlet[i].V[0][1] == inlet[j].V[0][1]) {
  2161 + 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][1] - main_feeder[0][1]) * (inlet[j].V[1][1] - main_feeder[0][1])) > 0 ? 1 : 0)) {
2070 2162 inlet_feasibility[i] = false;
2071 2163 break;
2072 2164 }
  2165 + else
  2166 + inlet_feasibility[i] = true;
2073 2167 }
2074   - else
2075   - inlet_feasibility[i] = true;
2076 2168 }
2077 2169 }
2078 2170 }
2079   -
2080 2171 // check outlet
2081 2172 num = outlet.size();
2082 2173 outlet_feasibility.resize(num, true);
... ... @@ -2084,7 +2175,7 @@ namespace stim {
2084 2175 for (unsigned j = 0; j < num; j++) {
2085 2176 if (i != j) {
2086 2177 if (outlet[i].V[0][2] == outlet[j].V[0][2]) {
2087   - 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]))) {
  2178 + 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][1] - main_feeder[1][1]) * (outlet[j].V[1][1] - main_feeder[1][1])) > 0 ? 1 : 0)) {
2088 2179 outlet_feasibility[i] = false;
2089 2180 break;
2090 2181 }
... ... @@ -2196,10 +2287,71 @@ namespace stim {
2196 2287 }
2197 2288 }
2198 2289  
  2290 + /// adjustment in order to match microfluidics experiments
  2291 + void adjust(T in, T out, T &Rt, T nQ, T viscosity, T radius = 5.0f) {
  2292 +
  2293 + // compute total resistance
  2294 + Rt = (Ps - Pe) / in;
  2295 + Pe = 0.0f;
  2296 +
  2297 + Ps = Rt * nQ;
2199 2298  
2200   - /// make binary image stack
  2299 + // adjust synthetic connections velocity flow rate. (linear scale)
  2300 + T k = nQ / in; // linear factor
  2301 + for (unsigned i = 0; i < inlet.size(); i++) {
  2302 + inlet[i].Q *= k;
  2303 + input[i].third *= k;
  2304 + }
  2305 + for (unsigned i = 0; i < outlet.size(); i++) {
  2306 + outlet[i].Q *= k;
  2307 + output[i].third *= k;
  2308 + }
  2309 +
  2310 + /// simulate inner network flow
  2311 + // clear up initialized pressure
  2312 + P.resize(num_vertex);
  2313 + for (unsigned i = 0; i < pendant_vertex.size(); i++) {
  2314 + unsigned index = UINT_MAX;
  2315 + for (unsigned j = 0; j < inlet.size(); j++) {
  2316 + if (inlet[j].v[0] == pendant_vertex[i]) {
  2317 + index = j;
  2318 + break;
  2319 + }
  2320 + }
  2321 + if (index != UINT_MAX) {
  2322 + P[inlet[index].v[0]] = Ps - ((T)8 * viscosity * inlet[index].l * inlet[index].Q / (stim::PI * std::pow(radius, 4)));
  2323 + }
  2324 + }
  2325 +
  2326 + for (unsigned i = 0; i < pendant_vertex.size(); i++) {
  2327 + unsigned index = UINT_MAX;
  2328 + for (unsigned j = 0; j < outlet.size(); j++) {
  2329 + if (outlet[j].v[0] == pendant_vertex[i]) {
  2330 + index = j;
  2331 + break;
  2332 + }
  2333 + }
  2334 + if (index != UINT_MAX) {
  2335 + P[outlet[index].v[0]] = Pe + ((T)8 * viscosity * outlet[index].l * outlet[index].Q / (stim::PI * std::pow(radius, 4)));
  2336 + }
  2337 + }
  2338 +
  2339 + // clear up previous simulation except synthetic connection parts
  2340 + for (unsigned i = 0; i < num_vertex; i++) {
  2341 + QQ[i] = 0;
  2342 + pressure[i] = 0;
  2343 + for (unsigned j = 0; j < num_vertex; j++) {
  2344 + C[i][j] = 0;
  2345 + }
  2346 + }
  2347 +
  2348 + // re-simulation
  2349 + solve_flow(viscosity);
  2350 + }
  2351 +
  2352 + /// make full-synthetic binary image stack
2201 2353 // prepare for image stack
2202   - void preparation(T &Xl, T &Xr, T &Yt, T &Yb, T &Z, T length = 40.0f, T height = 10.0f) {
  2354 + void preparation(T &Xl, T &Xr, T &Yt, T &Yb, T &Z, bool prototype = false, T length = 40.0f, T height = 10.0f, T radius = 5.0f) {
2203 2355  
2204 2356 T max_radius = 0.0f;
2205 2357 T top = FLT_MIN;
... ... @@ -2225,21 +2377,23 @@ namespace stim {
2225 2377 CU.push_back(new_cuboid);
2226 2378  
2227 2379 // take every point as sphere, every line as cone
2228   - for (unsigned i = 0; i < num_edge; i++) {
2229   - for (unsigned j = 0; j < E[i].size(); j++) {
2230   - new_sphere.c = E[i][j];
2231   - new_sphere.r = E[i].r(j);
2232   - A.push_back(new_sphere);
2233   - if (j != E[i].size() - 1) {
2234   - new_cone.c1 = E[i][j];
2235   - new_cone.c2 = E[i][j + 1];
2236   - new_cone.r1 = E[i].r(j);
2237   - new_cone.r2 = E[i].r(j + 1);
2238   - B.push_back(new_cone);
  2380 + if (!prototype) {
  2381 + for (unsigned i = 0; i < num_edge; i++) {
  2382 + for (unsigned j = 0; j < E[i].size(); j++) {
  2383 + new_sphere.c = E[i][j];
  2384 + new_sphere.r = E[i].r(j);
  2385 + A.push_back(new_sphere);
  2386 + if (j != E[i].size() - 1) {
  2387 + new_cone.c1 = E[i][j];
  2388 + new_cone.c2 = E[i][j + 1];
  2389 + new_cone.r1 = E[i].r(j);
  2390 + new_cone.r2 = E[i].r(j + 1);
  2391 + B.push_back(new_cone);
  2392 + }
2239 2393 }
2240 2394 }
2241 2395 }
2242   -
  2396 +
2243 2397 // secondly push back outside connection
2244 2398 for (unsigned i = 0; i < inlet.size(); i++) {
2245 2399 for (unsigned j = 1; j < inlet[i].V.size() - 1; j++) {
... ... @@ -2288,34 +2442,38 @@ namespace stim {
2288 2442 max_radius = A[i].r;
2289 2443 }
2290 2444  
2291   - Yt = top; // top bound y coordinate
2292   - Yb = bottom; // bottom bound y coordinate
  2445 + Yt = top + 2 * radius; // top bound y coordinate
  2446 + Yb = bottom + 2 * radius; // bottom bound y coordinate
2293 2447 Z = (bb.B[2] - bb.A[2] + 2 * max_radius); // bounding box width(along z-axis)
2294 2448 }
2295 2449  
2296 2450 /// making image stack main function
2297   - void make_image_stack(T dx, T dy, T dz, std::string stackdir, T radius = 5.0f) {
  2451 + void make_image_stack(stim::image_stack<unsigned char, T> I, T dx, T dy, T dz, std::string stackdir, bool prototype = false, T radius = 5.0f) {
2298 2452  
2299 2453 /// preparation for making image stack
2300 2454 T X, Xl, Xr, Y, Yt, Yb, Z;
2301   - preparation(Xl, Xr, Yt, Yb, Z);
  2455 + preparation(Xl, Xr, Yt, Yb, Z, prototype);
2302 2456 X = Xr - Xl; // bounding box length(along x-axis)
2303 2457 Y = Yt - Yb; // bounding box height(along y-axis)
2304   -
2305   - /// make
2306   - stim::image_stack<unsigned char, T> I;
2307   - T size_x, size_y, size_z;
2308   -
2309 2458 stim::vec3<T> center = bb.center(); // get the center of bounding box
  2459 + T size_x, size_y, size_z;
2310 2460  
2311   - size_x = X / dx + 1; // set the size of image
2312   - size_y = Y / dy + 1;
2313   - size_z = Z / dz + 1;
2314   -
2315   - /// initialize image stack object
2316   - I.init(1, size_x, size_y, size_z);
2317   - I.set_dim(dx, dy, dz);
2318   -
  2461 + if (!prototype) {
  2462 + /// make
  2463 + stim::image_stack<unsigned char, T> I;
  2464 + size_x = X / dx + 1; // set the size of image
  2465 + size_y = Y / dy + 1;
  2466 + size_z = Z / dz + 1;
  2467 + /// initialize image stack object
  2468 + I.init(1, size_x, size_y, size_z);
  2469 + I.set_dim(dx, dy, dz);
  2470 + }
  2471 + else {
  2472 + size_x = I.nx();
  2473 + size_y = I.ny();
  2474 + size_z = I.nz();
  2475 + }
  2476 +
2319 2477 // because of lack of memory, we have to computer one slice of stack per time
2320 2478 // allocate vertex, edge and bus
2321 2479 stim::sphere<T> *d_V;
... ... @@ -2363,6 +2521,8 @@ namespace stim {
2363 2521 memset(ptr, 0, num * sizeof(unsigned char));
2364 2522  
2365 2523 HANDLE_ERROR(cudaMalloc((void**)&d_ptr, num * sizeof(unsigned char)));
  2524 + if (prototype) // load prototype image stack if provided
  2525 + HANDLE_ERROR(cudaMemcpy(d_ptr, &I.data()[i * num], num * sizeof(unsigned char), cudaMemcpyHostToDevice));
2366 2526  
2367 2527 cudaDeviceProp prop;
2368 2528 cudaGetDeviceProperties(&prop, 0); // get cuda device properties structure
... ... @@ -2403,6 +2563,28 @@ namespace stim {
2403 2563 I.save_images(stackdir + "/image????.bmp");
2404 2564 }
2405 2565  
  2566 + /// save network flow profile
  2567 + void save_network() {
  2568 +
  2569 + // save the pressure information to CSV file
  2570 + std::string p_filename = "profile/pressure.csv";
  2571 + std::ofstream p_file;
  2572 + p_file.open(p_filename.c_str());
  2573 + p_file << "Vertex, Pressure(g/" << units << "/s^2)" << std::endl;
  2574 + for (unsigned i = 0; i < num_vertex; i++)
  2575 + p_file << i << "," << pressure[i] << std::endl;
  2576 + p_file.close();
  2577 +
  2578 + // save the flow information to CSV file
  2579 + std::string f_filename = "profile/flow_rate.csv";
  2580 + std::ofstream f_file;
  2581 + f_file.open(f_filename.c_str());
  2582 + f_file << "Edge, Volume flow rate(" << units << "^3/s)" << std::endl;
  2583 + for (unsigned i = 0; i < num_edge; i++)
  2584 + f_file << Q[i].first << "->" << Q[i].second << "," << Q[i].third << std::endl;
  2585 + f_file.close();
  2586 + }
  2587 +
2406 2588 /// Calculate the inverse of A and store the result in C
2407 2589 void inversion(T** A, int order, T* C) {
2408 2590  
... ...
... ... @@ -23,10 +23,12 @@
23 23 #include <stim/visualization/camera.h>
24 24 #include <stim/visualization/colormap.h>
25 25 #include <stim/cuda/cudatools/error.h>
  26 +#include <stim/grids/image_stack.h>
26 27  
27 28  
28 29 //********************parameter setting********************
29 30 // overall parameters
  31 +std::string units; // units used in this program
30 32 int vX, vY;
31 33 float dx, dy, dz; // x, y and z image scaling(units/pixel)
32 34 std::string stackdir = ""; // directory where image stacks will be stored
... ... @@ -36,7 +38,7 @@ stim::camera cam; // camera object
36 38 unsigned num_edge; // number of edges in the network
37 39 unsigned num_vertex; // number of vertex in the network
38 40 std::vector<unsigned> pendant_vertex; // list of pendant vertex index in GT
39   -std::vector<std::string> menu_option = { "simulation", "build inlet/outlet", "manufacture" };
  41 +std::vector<std::string> menu_option = { "simulation", "build inlet/outlet", "manufacture", "adjustment" };
40 42 stim::flow<float> flow; // flow object
41 43 float move_pace; // camera moving parameter
42 44 float u; // viscosity
... ... @@ -48,6 +50,13 @@ std::vector&lt;unsigned char&gt; color; // velocity color map
48 50 std::vector<int> velocity_bar; // velocity bar
49 51 float length = 40.0f; // cuboid length
50 52 float scale = 1.0f; // render scale factor
  53 +bool image_stack = false; // flag indicates an image stack been loaded
  54 +stim::image_stack<unsigned char, float> S; // image stack
  55 +float binary_threshold = 128; // threshold for binary transformation
  56 +float in = 0.0f; // total input volume flow rate
  57 +float out = 0.0f;
  58 +float Rt = 0.0f; // total resistance
  59 +float Qn = 0.0f; // required input total volume flow rate
51 60  
52 61 // hard-coded parameters
53 62 float camera_factor = 1.2f; // start point of the camera as a function of X and Y size
... ... @@ -134,6 +143,21 @@ inline void window_to_world(GLdouble &amp;x, GLdouble &amp;y, GLdouble &amp;z) {
134 143 gluUnProject(winX, winY, winZ, modelview, projection, viewport, &x, &y, &z);
135 144 }
136 145  
  146 +// convert current image stack into a binary mask
  147 +#ifdef __CUDACC__
  148 +template <typename T, typename F>
  149 +__global__ void binary_transform(unsigned N, T* ptr, F threshold) {
  150 +
  151 + unsigned ix = blockDim.x * blockIdx.x + threadIdx.x;
  152 + if (ix >= N) return; // avoid seg-fault
  153 +
  154 + if (ptr[ix] >= threshold) // binary transformation
  155 + ptr[ix] = 255;
  156 + else
  157 + ptr[ix] = 0;
  158 +}
  159 +#endif
  160 +
137 161  
138 162 //********************simulation function**********************
139 163 // initialize flow object
... ... @@ -163,6 +187,16 @@ void flow_stable_state() {
163 187 std::sort(velocity_bar.begin(), velocity_bar.end(), [&](int x, int y) {return abs(flow.v[x]) < abs(flow.v[y]); });
164 188 }
165 189  
  190 +// adjustment on input volume flow rate and corresponding flow simulation
  191 +void adjustment() {
  192 +
  193 + system("CLS"); // clear up console box
  194 + std::cout << "Please enter the input total volume flow rate: " << std::endl;
  195 + std::cin >> Qn;
  196 +
  197 + flow.adjust(in, out, Rt, Qn, u);
  198 +}
  199 +
166 200  
167 201 //********************glut function********************
168 202 // dynamically set menu
... ... @@ -226,7 +260,7 @@ void glut_render() {
226 260 flow.mark_vertex(scale);
227 261 //flow.glSolidCone(subdivision);
228 262 flow.glSolidCylinder(direction_index, color, scale, subdivision);
229   - flow.glSolidCuboid(length);
  263 + flow.glSolidCuboid(manufacture, length);
230 264 if (render_direction)
231 265 flow.glSolidCone(direction_index, scale, subdivision);
232 266 }
... ... @@ -236,7 +270,7 @@ void glut_render() {
236 270 }
237 271  
238 272 if (manufacture) {
239   - flow.glSolidCuboid();
  273 + flow.glSolidCuboid(manufacture, length);
240 274 flow.tube_bridge(subdivision);
241 275 }
242 276  
... ... @@ -402,6 +436,9 @@ void glut_render() {
402 436 glPopMatrix();
403 437 }
404 438  
  439 + if (build_inlet_outlet)
  440 + flow.display_flow_rate(in, out);
  441 +
405 442 glutSwapBuffers();
406 443 }
407 444  
... ... @@ -416,7 +453,7 @@ void glut_menu(int value) {
416 453 modified_bridge = false;
417 454 connection_done = false;
418 455 // first time
419   - if (!flow.set) {
  456 + if (!flow.set) { // only first time simulation called "simulation", ^_^
420 457 flow_initialize();
421 458 menu_option[0] = "resimulation";
422 459 }
... ... @@ -444,7 +481,7 @@ void glut_menu(int value) {
444 481 flow.clear_synthetic_connection();
445 482 }
446 483  
447   - glut_set_menu(num, 3);
  484 + glut_set_menu(num, 4);
448 485 }
449 486  
450 487 if (value == 3) {
... ... @@ -453,6 +490,16 @@ void glut_menu(int value) {
453 490 manufacture = true;
454 491 }
455 492  
  493 + if (value == 4) {
  494 + simulation = true;
  495 + build_inlet_outlet = false;
  496 + manufacture = false;
  497 +
  498 + adjustment(); // adjust network flow accordingly
  499 +
  500 + glut_set_menu(num, 1);
  501 + }
  502 +
456 503 glutPostRedisplay();
457 504 }
458 505  
... ... @@ -580,21 +627,24 @@ void glut_mouse(int button, int state, int x, int y) {
580 627 }
581 628 }
582 629 else if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN && !mods && simulation && to_select_pressure) {
583   - if (y > 2 * border_factor || y < vY - border_factor) { // within the pressure bar range
  630 + if (y >= 2 * border_factor && y <= vY - border_factor) { // within the pressure bar range
584 631 to_select_pressure = false;
585   - float tmp_pressure = (float)(vY - y - border_factor) / ((float)vY - border_factor) * max_pressure;
  632 + float tmp_pressure = (float)(vY - y - border_factor) / ((float)vY - 3.0f * border_factor) * max_pressure;
586 633 flow.set_pressure(pressure_index, tmp_pressure);
587   -
588 634 //flow_stable_state(); // main function of solving the linear system
589 635 //flow.print_flow();
590 636 }
591 637 }
592 638 else if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN && !mods && modified_bridge && change_fragment) {
593   - if (y > 2 * border_factor || y < vY - border_factor) { // within the ratio bar range
594   - fragment_ratio = (float)(vY - y - border_factor) / ((float)vY - border_factor) * 1.0f;
595   - flow.modify_synthetic_connection(u, rou, hilbert_curve, height_threshold, fragment_ratio, default_radius);
596   - change_fragment = false;
597   - }
  639 + if (y >= 2 * border_factor && y <= vY - border_factor) // within the ratio bar range
  640 + fragment_ratio = (float)(vY - y - border_factor) / ((float)vY - 3.0f * border_factor) * 1.0f;
  641 + else if (y < 2 * border_factor)
  642 + fragment_ratio = 1.0f;
  643 + else if (y > vY - border_factor)
  644 + fragment_ratio = 0.0f;
  645 +
  646 + change_fragment = false;
  647 + flow.modify_synthetic_connection(u, rou, hilbert_curve, height_threshold, in, out, fragment_ratio, default_radius);
598 648 }
599 649 // move connections along y-axis
600 650 else if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN && mods == GLUT_ACTIVE_SHIFT && !modified_bridge && !picked_connection) {
... ... @@ -789,17 +839,12 @@ void glut_keyboard(unsigned char key, int x, int y) {
789 839  
790 840 // register different keyboard operation
791 841 switch (key) {
792   -
793   - // zooming
794   - case 'w': // if keyboard 'w' is pressed, then move closer
795   - move_pace = zoom_factor;
796   - cam.Push(move_pace);
797   - break;
798   - case 's': // if keyboard 's' is pressed, then leave farther
799   - move_pace = -zoom_factor;
800   - cam.Push(move_pace);
  842 +
  843 + // saving network flow profile
  844 + case 's':
  845 + flow.save_network();
801 846 break;
802   -
  847 +
803 848 // open/close index marks
804 849 case 'e':
805 850 if (mark_index)
... ... @@ -812,7 +857,8 @@ void glut_keyboard(unsigned char key, int x, int y) {
812 857 case 'm':
813 858 if (manufacture) {
814 859 #ifdef __CUDACC__
815   - flow.make_image_stack(dx, dy, dz, stackdir);
  860 + flow.make_image_stack(S, dx, dy, dz, stackdir, image_stack);
  861 +
816 862 #else
817 863 std::cout << "You need to have a gpu to make image stack, sorry." << std::endl;
818 864 #endif
... ... @@ -821,7 +867,7 @@ void glut_keyboard(unsigned char key, int x, int y) {
821 867 modified_bridge = true;
822 868  
823 869 if (hilbert_curve)
824   - flow.modify_synthetic_connection(u, rou, hilbert_curve, height_threshold);
  870 + flow.modify_synthetic_connection(u, rou, hilbert_curve, height_threshold, in, out, fragment_ratio);
825 871 else
826 872 change_fragment = true;
827 873 }
... ... @@ -865,7 +911,7 @@ void glut_initialize() {
865 911 void advertise() {
866 912 std::cout << std::endl << std::endl;
867 913 std::cout << " =======================================================================================" << std::endl;
868   - std::cout << "|Thank you for using the synthetic microvascular model generator for microfluidics tool!|" << std::endl;
  914 + std::cout << "|Thank you for using the AFAN tool! |" << std::endl;
869 915 std::cout << "|Scalable Tissue Imaging and Modeling (STIM) Lab, University of Houston |" << std::endl;
870 916 std::cout << "|Developers: Jiaming Guo, David Mayerich |" << std::endl;
871 917 std::cout << "|Source: https://git.stim.ee.uh.edu/Jack/flow3.git |" << std::endl;
... ... @@ -879,11 +925,12 @@ int main(int argc, char* argv[]) {
879 925  
880 926 // add arguments
881 927 args.add("help", "prints the help");
882   - //args.add("network", "load network from .obj or .swc file");
  928 + args.add("units", "string indicating units of length for output measurements (ex. velocity)", "um", "text string");
883 929 args.add("maxpress", "maximum allowed pressure in g / units / s^2, default 2 is for blood when units = um", "2", "real value > 0");
884 930 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");
885 931 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");
886 932 args.add("hilbert", "activate hilbert curves connections", "0", "value 1 for enablement");
  933 + args.add("stack", "load the image stack");
887 934 args.add("stackres", "spacing between pixel samples in each dimension(in units/pixel)", "1 1 1", "real value > 0");
888 935 args.add("stackdir", "set the directory of the output image stack", "", "any existing directory (ex. /home/name/network)");
889 936  
... ... @@ -912,6 +959,10 @@ int main(int argc, char* argv[]) {
912 959 }
913 960 get_background();
914 961  
  962 + // get the units to work on
  963 + units = args["units"].as_string();
  964 + flow.set_units(units);
  965 +
915 966 // blood pressure in capillaries range from 15 - 35 torr
916 967 // 1 torr = 133.3 Pa
917 968 max_pressure = args["maxpress"].as_float();
... ... @@ -926,6 +977,31 @@ int main(int argc, char* argv[]) {
926 977 // check whether to enable hilbert curves or not
927 978 hilbert_curve = args["hilbert"].as_int();
928 979  
  980 + // load image stack if provided
  981 + if (args["stack"].is_set()) {
  982 + image_stack = true;
  983 + S.load_images(args["stack"].as_string());
  984 + // binary transformation
  985 +#ifdef __CUDACC__
  986 + unsigned N = S.samples(); // number of pixels loaded
  987 + unsigned char* d_S; // image stack stored in device
  988 + unsigned char* h_S = (unsigned char*)malloc(N * sizeof(unsigned char)); // image stack stored in host
  989 + cudaMalloc((void**)&d_S, N * sizeof(unsigned char));
  990 + cudaMemcpy(d_S, S.data(), N * sizeof(unsigned char), cudaMemcpyHostToDevice);
  991 +
  992 +
  993 + unsigned thread = 1024;
  994 + unsigned block = N / thread + 1;
  995 + binary_transform <<<block, thread>>> (N, d_S, binary_threshold);
  996 +
  997 + cudaMemcpy(h_S, d_S, N * sizeof(unsigned char), cudaMemcpyDeviceToHost);
  998 +
  999 + S.copy(h_S);
  1000 +
  1001 + S.save_images("image????.bmp");
  1002 +#endif
  1003 + }
  1004 +
929 1005 // get the vexel and image stack size
930 1006 dx = args["stackres"].as_float(0);
931 1007 dy = args["stackres"].as_float(1);
... ...