codebase / stimlib

Commit 9db76bdba54c5e3535560e29c9b9112237aa0f1d

Authored by Jiaming Guo
1 parent 28b882ce

change flow object

Showing 1 changed file with 92 additions and 71 deletions   Show diff stats
Inline Side-by-side
stim/biomodels/flow.h
  Show/Hide comments   View file @9db76bd
... ... @@ -92,82 +92,30 @@ namespace stim {
92 92 }
93 93  
94 94 public:
95   - std::vector<typename stim::triple<unsigned, unsigned, T> > V; // velocity
96   - std::vector<typename stim::triple<unsigned, unsigned, T> > Q; // volume flow rate
97   - std::vector<typename stim::triple<unsigned, unsigned, T> > deltaP; // pressure drop
98   - std::vector<T> pressure; // pressure
  95 + T** C; // Conductance
  96 + std::vector<typename stim::triple<unsigned, unsigned, float> > Q; // volume flow rate
  97 + std::vector<T> QQ; // Q' vector
  98 + std::vector<T> P; // initial pressure
  99 + std::vector<T> pressure; // final pressure
99 100  
100   - flow() {} // default constructor
101   -
102   - // calculate the flow rate of 3D model(circle cross section)
103   - void calculate_flow_rate(unsigned e, T r) {
104   - stim::triple<unsigned, unsigned, T> tmp_Q;
105   - tmp_Q.first = V[e].first; // copy the vertices information
106   - tmp_Q.second = V[e].second;
107   - tmp_Q.third = V[e].third * stim::PI * pow(r, 2); // UNITS: uL/s
108   - Q.push_back(tmp_Q); // push back the volume flow rate information for every edge
109   - }
110   -
111   - // calculate the flow rate of 2D model(rectangular cross section)
112   - void calculate_flow_rate(unsigned e, T r, T h) {
113   - stim::triple<unsigned, unsigned, T> tmp_Q;
114   - tmp_Q.first = V[e].first; // copy the vertices information
115   - tmp_Q.second = V[e].second;
116   - tmp_Q.third = V[e].third * h * r; // UNITS: uL/s = mm^3/s
117   - Q.push_back(tmp_Q); // push back the volume flow rate information for every edge
118   - }
119   -
120   - // calculate the pressure drop of 3D model(circle cross section)
121   - void calculate_deltaP(unsigned e, T u, T l, T r) {
122   - stim::triple<unsigned, unsigned, T> tmp_deltaP;
123   - tmp_deltaP.first = V[e].first; // copy the vertices information
124   - tmp_deltaP.second = V[e].second;
125   - tmp_deltaP.third = (8 * u * l * Q[e].third) / (stim::PI * pow(r, 4)); // UNITS: g/mm/s^2 = Pa
126   - deltaP.push_back(tmp_deltaP); // push back the volume flow rate information for every edge
127   - }
  101 + //std::vector<typename stim::triple<unsigned, unsigned, T> > V; // velocity
  102 + //std::vector<typename stim::triple<unsigned, unsigned, T> > Q; // volume flow rate
  103 + //std::vector<typename stim::triple<unsigned, unsigned, T> > deltaP; // pressure drop
128 104  
129   - // calculate the pressure drop of 2D model(rectangular cross section)
130   - void calculate_deltaP(unsigned e, T u, T l, T r, T h) {
131   - stim::triple<unsigned, unsigned, T> tmp_deltaP;
132   - tmp_deltaP.first = V[e].first; // copy the vertices information
133   - tmp_deltaP.second = V[e].second;
134   - tmp_deltaP.third = (12 * u * l * Q[e].third) / (stim::PI * h * pow(r, 3)); // UNITS: g/mm/s^2 = Pa
135   - deltaP.push_back(tmp_deltaP); // push back the volume flow rate information for every edge
136   - }
  105 + flow() {} // default constructor
137 106  
138   - // better way to do this???
139   - // find the maximum and minimum pressure positions
140   - void find_max_min_pressure(size_t n_e, size_t n_v, unsigned& max, unsigned& min) {
141   - std::vector<T> P(n_v, FLT_MAX);
142   - // set one to reference
143   - P[Q[0].first] = 0.0;
144   - unsigned first = 0;
145   - unsigned second = 0;
146   - // calculate all the relative pressure in brute force manner
147   - for (unsigned e = 0; e < n_e; e++) {
148   - // assuming the obj file stores in a straight order, in other words, like swc file
149   - first = Q[e].first;
150   - second = Q[e].second;
151   - if (P[first] != FLT_MAX) // if pressure at start vertex is known
152   - P[second] = P[first] - deltaP[e].third;
153   - else if (P[second] != FLT_MAX) // if pressure at end vertex is known
154   - P[first] = P[second] + deltaP[e].third;
  107 + void init(unsigned n_e, unsigned n_v) {
  108 +
  109 + C = new T*[n_v]();
  110 + for (unsigned i = 0; i < n_v; i++) {
  111 + C[i] = new T[n_v]();
155 112 }
156 113  
157   - // find the maximum and minimum pressure position
158   - auto m1 = std::max_element(P.begin(), P.end()); // temporarily max number
159   - auto m2 = std::min_element(P.begin(), P.end()); // temporarily min number
160   -
161   - max = std::distance(P.begin(), m1);
162   - min = std::distance(P.begin(), m2);
163   -
164   - T tmp_m = *m2;
165   - // Now set the lowest pressure port to reference pressure(0.0 Pa)
166   - for (unsigned i = 0; i < n_v; i++)
167   - P[i] -= tmp_m;
  114 + QQ.resize(n_v);
  115 + P.resize(n_v);
  116 + pressure.resize(n_v);
168 117  
169   - for (unsigned i = 0; i < n_v; i++)
170   - pressure.push_back(P[i]);
  118 + Q.resize(n_e);
171 119 }
172 120  
173 121 void inversion(T** A, int order, T** C) {
... ... @@ -270,4 +218,77 @@ namespace stim {
270 218 };
271 219 }
272 220  
273   -#endif
274 221 \ No newline at end of file
  222 +#endif
  223 +
  224 +
  225 +
  226 +//// calculate the flow rate of 3D model(circle cross section)
  227 +//void calculate_flow_rate(unsigned e, T r) {
  228 +// stim::triple<unsigned, unsigned, T> tmp_Q;
  229 +// tmp_Q.first = V[e].first; // copy the vertices information
  230 +// tmp_Q.second = V[e].second;
  231 +// tmp_Q.third = V[e].third * stim::PI * pow(r, 2); // UNITS: uL/s
  232 +// Q.push_back(tmp_Q); // push back the volume flow rate information for every edge
  233 +//}
  234 +
  235 +//// calculate the flow rate of 2D model(rectangular cross section)
  236 +//void calculate_flow_rate(unsigned e, T r, T h) {
  237 +// stim::triple<unsigned, unsigned, T> tmp_Q;
  238 +// tmp_Q.first = V[e].first; // copy the vertices information
  239 +// tmp_Q.second = V[e].second;
  240 +// tmp_Q.third = V[e].third * h * r; // UNITS: uL/s = mm^3/s
  241 +// Q.push_back(tmp_Q); // push back the volume flow rate information for every edge
  242 +//}
  243 +
  244 +//// calculate the pressure drop of 3D model(circle cross section)
  245 +//void calculate_deltaP(unsigned e, T u, T l, T r) {
  246 +// stim::triple<unsigned, unsigned, T> tmp_deltaP;
  247 +// tmp_deltaP.first = V[e].first; // copy the vertices information
  248 +// tmp_deltaP.second = V[e].second;
  249 +// tmp_deltaP.third = (8 * u * l * Q[e].third) / (stim::PI * pow(r, 4)); // UNITS: g/mm/s^2 = Pa
  250 +// deltaP.push_back(tmp_deltaP); // push back the volume flow rate information for every edge
  251 +//}
  252 +
  253 +//// calculate the pressure drop of 2D model(rectangular cross section)
  254 +//void calculate_deltaP(unsigned e, T u, T l, T r, T h) {
  255 +// stim::triple<unsigned, unsigned, T> tmp_deltaP;
  256 +// tmp_deltaP.first = V[e].first; // copy the vertices information
  257 +// tmp_deltaP.second = V[e].second;
  258 +// tmp_deltaP.third = (12 * u * l * Q[e].third) / (h * pow(r, 3)); // UNITS: g/mm/s^2 = Pa
  259 +// deltaP.push_back(tmp_deltaP); // push back the volume flow rate information for every edge
  260 +//}
  261 +
  262 +//// better way to do this???
  263 +//// find the maximum and minimum pressure positions
  264 +//void find_max_min_pressure(size_t n_e, size_t n_v, unsigned& max, unsigned& min) {
  265 +// std::vector<T> P(n_v, FLT_MAX);
  266 +// // set one to reference
  267 +// P[Q[0].first] = 0.0;
  268 +// unsigned first = 0;
  269 +// unsigned second = 0;
  270 +// // calculate all the relative pressure in brute force manner
  271 +// for (unsigned e = 0; e < n_e; e++) {
  272 +// // assuming the obj file stores in a straight order, in other words, like swc file
  273 +// first = Q[e].first;
  274 +// second = Q[e].second;
  275 +// if (P[first] != FLT_MAX) // if pressure at start vertex is known
  276 +// P[second] = P[first] - deltaP[e].third;
  277 +// else if (P[second] != FLT_MAX) // if pressure at end vertex is known
  278 +// P[first] = P[second] + deltaP[e].third;
  279 +// }
  280 +
  281 +// // find the maximum and minimum pressure position
  282 +// auto m1 = std::max_element(P.begin(), P.end()); // temporarily max number
  283 +// auto m2 = std::min_element(P.begin(), P.end()); // temporarily min number
  284 +
  285 +// max = std::distance(P.begin(), m1);
  286 +// min = std::distance(P.begin(), m2);
  287 +
  288 +// T tmp_m = *m2;
  289 +// // Now set the lowest pressure port to reference pressure(0.0 Pa)
  290 +// for (unsigned i = 0; i < n_v; i++)
  291 +// P[i] -= tmp_m;
  292 +
  293 +// for (unsigned i = 0; i < n_v; i++)
  294 +// pressure.push_back(P[i]);
  295 +//}
275 296 \ No newline at end of file
... ...