From 9db76bdba54c5e3535560e29c9b9112237aa0f1d Mon Sep 17 00:00:00 2001 From: Jiaming Guo Date: Wed, 15 Feb 2017 11:39:58 -0600 Subject: [PATCH] change flow object --- stim/biomodels/flow.h | 163 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------------------------------- 1 file changed, 92 insertions(+), 71 deletions(-) diff --git a/stim/biomodels/flow.h b/stim/biomodels/flow.h index 4990c02..bdb3c9d 100644 --- a/stim/biomodels/flow.h +++ b/stim/biomodels/flow.h @@ -92,82 +92,30 @@ namespace stim { } public: - std::vector > V; // velocity - std::vector > Q; // volume flow rate - std::vector > deltaP; // pressure drop - std::vector pressure; // pressure + T** C; // Conductance + std::vector > Q; // volume flow rate + std::vector QQ; // Q' vector + std::vector P; // initial pressure + std::vector pressure; // final pressure - flow() {} // default constructor - - // calculate the flow rate of 3D model(circle cross section) - void calculate_flow_rate(unsigned e, T r) { - stim::triple tmp_Q; - tmp_Q.first = V[e].first; // copy the vertices information - tmp_Q.second = V[e].second; - tmp_Q.third = V[e].third * stim::PI * pow(r, 2); // UNITS: uL/s - Q.push_back(tmp_Q); // push back the volume flow rate information for every edge - } - - // calculate the flow rate of 2D model(rectangular cross section) - void calculate_flow_rate(unsigned e, T r, T h) { - stim::triple tmp_Q; - tmp_Q.first = V[e].first; // copy the vertices information - tmp_Q.second = V[e].second; - tmp_Q.third = V[e].third * h * r; // UNITS: uL/s = mm^3/s - Q.push_back(tmp_Q); // push back the volume flow rate information for every edge - } - - // calculate the pressure drop of 3D model(circle cross section) - void calculate_deltaP(unsigned e, T u, T l, T r) { - stim::triple tmp_deltaP; - tmp_deltaP.first = V[e].first; // copy the vertices information - tmp_deltaP.second = V[e].second; - tmp_deltaP.third = (8 * u * l * Q[e].third) / (stim::PI * pow(r, 4)); // UNITS: g/mm/s^2 = Pa - deltaP.push_back(tmp_deltaP); // push back the volume flow rate information for every edge - } + //std::vector > V; // velocity + //std::vector > Q; // volume flow rate + //std::vector > deltaP; // pressure drop - // calculate the pressure drop of 2D model(rectangular cross section) - void calculate_deltaP(unsigned e, T u, T l, T r, T h) { - stim::triple tmp_deltaP; - tmp_deltaP.first = V[e].first; // copy the vertices information - tmp_deltaP.second = V[e].second; - tmp_deltaP.third = (12 * u * l * Q[e].third) / (stim::PI * h * pow(r, 3)); // UNITS: g/mm/s^2 = Pa - deltaP.push_back(tmp_deltaP); // push back the volume flow rate information for every edge - } + flow() {} // default constructor - // better way to do this??? - // find the maximum and minimum pressure positions - void find_max_min_pressure(size_t n_e, size_t n_v, unsigned& max, unsigned& min) { - std::vector P(n_v, FLT_MAX); - // set one to reference - P[Q[0].first] = 0.0; - unsigned first = 0; - unsigned second = 0; - // calculate all the relative pressure in brute force manner - for (unsigned e = 0; e < n_e; e++) { - // assuming the obj file stores in a straight order, in other words, like swc file - first = Q[e].first; - second = Q[e].second; - if (P[first] != FLT_MAX) // if pressure at start vertex is known - P[second] = P[first] - deltaP[e].third; - else if (P[second] != FLT_MAX) // if pressure at end vertex is known - P[first] = P[second] + deltaP[e].third; + void init(unsigned n_e, unsigned n_v) { + + C = new T*[n_v](); + for (unsigned i = 0; i < n_v; i++) { + C[i] = new T[n_v](); } - // find the maximum and minimum pressure position - auto m1 = std::max_element(P.begin(), P.end()); // temporarily max number - auto m2 = std::min_element(P.begin(), P.end()); // temporarily min number - - max = std::distance(P.begin(), m1); - min = std::distance(P.begin(), m2); - - T tmp_m = *m2; - // Now set the lowest pressure port to reference pressure(0.0 Pa) - for (unsigned i = 0; i < n_v; i++) - P[i] -= tmp_m; + QQ.resize(n_v); + P.resize(n_v); + pressure.resize(n_v); - for (unsigned i = 0; i < n_v; i++) - pressure.push_back(P[i]); + Q.resize(n_e); } void inversion(T** A, int order, T** C) { @@ -270,4 +218,77 @@ namespace stim { }; } -#endif \ No newline at end of file +#endif + + + +//// calculate the flow rate of 3D model(circle cross section) +//void calculate_flow_rate(unsigned e, T r) { +// stim::triple tmp_Q; +// tmp_Q.first = V[e].first; // copy the vertices information +// tmp_Q.second = V[e].second; +// tmp_Q.third = V[e].third * stim::PI * pow(r, 2); // UNITS: uL/s +// Q.push_back(tmp_Q); // push back the volume flow rate information for every edge +//} + +//// calculate the flow rate of 2D model(rectangular cross section) +//void calculate_flow_rate(unsigned e, T r, T h) { +// stim::triple tmp_Q; +// tmp_Q.first = V[e].first; // copy the vertices information +// tmp_Q.second = V[e].second; +// tmp_Q.third = V[e].third * h * r; // UNITS: uL/s = mm^3/s +// Q.push_back(tmp_Q); // push back the volume flow rate information for every edge +//} + +//// calculate the pressure drop of 3D model(circle cross section) +//void calculate_deltaP(unsigned e, T u, T l, T r) { +// stim::triple tmp_deltaP; +// tmp_deltaP.first = V[e].first; // copy the vertices information +// tmp_deltaP.second = V[e].second; +// tmp_deltaP.third = (8 * u * l * Q[e].third) / (stim::PI * pow(r, 4)); // UNITS: g/mm/s^2 = Pa +// deltaP.push_back(tmp_deltaP); // push back the volume flow rate information for every edge +//} + +//// calculate the pressure drop of 2D model(rectangular cross section) +//void calculate_deltaP(unsigned e, T u, T l, T r, T h) { +// stim::triple tmp_deltaP; +// tmp_deltaP.first = V[e].first; // copy the vertices information +// tmp_deltaP.second = V[e].second; +// tmp_deltaP.third = (12 * u * l * Q[e].third) / (h * pow(r, 3)); // UNITS: g/mm/s^2 = Pa +// deltaP.push_back(tmp_deltaP); // push back the volume flow rate information for every edge +//} + +//// better way to do this??? +//// find the maximum and minimum pressure positions +//void find_max_min_pressure(size_t n_e, size_t n_v, unsigned& max, unsigned& min) { +// std::vector P(n_v, FLT_MAX); +// // set one to reference +// P[Q[0].first] = 0.0; +// unsigned first = 0; +// unsigned second = 0; +// // calculate all the relative pressure in brute force manner +// for (unsigned e = 0; e < n_e; e++) { +// // assuming the obj file stores in a straight order, in other words, like swc file +// first = Q[e].first; +// second = Q[e].second; +// if (P[first] != FLT_MAX) // if pressure at start vertex is known +// P[second] = P[first] - deltaP[e].third; +// else if (P[second] != FLT_MAX) // if pressure at end vertex is known +// P[first] = P[second] + deltaP[e].third; +// } + +// // find the maximum and minimum pressure position +// auto m1 = std::max_element(P.begin(), P.end()); // temporarily max number +// auto m2 = std::min_element(P.begin(), P.end()); // temporarily min number + +// max = std::distance(P.begin(), m1); +// min = std::distance(P.begin(), m2); + +// T tmp_m = *m2; +// // Now set the lowest pressure port to reference pressure(0.0 Pa) +// for (unsigned i = 0; i < n_v; i++) +// P[i] -= tmp_m; + +// for (unsigned i = 0; i < n_v; i++) +// pressure.push_back(P[i]); +//} \ No newline at end of file -- libgit2 0.21.4