Commit 0a04ac9dffd53f742ed85362d6a9a03b2fd5eb33
Merge branch 'master' of https://git.stim.ee.uh.edu/codebase/stimlib
Showing
1 changed file
with
92 additions
and
71 deletions
Show diff stats
stim/biomodels/flow.h
@@ -92,82 +92,30 @@ namespace stim { | @@ -92,82 +92,30 @@ namespace stim { | ||
92 | } | 92 | } |
93 | 93 | ||
94 | public: | 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 | void inversion(T** A, int order, T** C) { | 121 | void inversion(T** A, int order, T** C) { |
@@ -270,4 +218,77 @@ namespace stim { | @@ -270,4 +218,77 @@ namespace stim { | ||
270 | }; | 218 | }; |
271 | } | 219 | } |
272 | 220 | ||
273 | -#endif | ||
274 | \ No newline at end of file | 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 | \ No newline at end of file | 296 | \ No newline at end of file |