From 8f6c2057eeb0eb428641e9cd53a7f437139dbdfb Mon Sep 17 00:00:00 2001 From: Jiaming Guo Date: Thu, 14 Sep 2017 14:33:43 -0500 Subject: [PATCH] fixed index marking bugs after network subdivided --- flow.h | 62 +++++++++++++++++++++++++++++++++++++++++++++++++++----------- main.cu | 72 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------- 2 files changed, 110 insertions(+), 24 deletions(-) diff --git a/flow.h b/flow.h index f505aa5..0c43ffb 100644 --- a/flow.h +++ b/flow.h @@ -67,7 +67,7 @@ namespace stim { #ifdef __CUDACC__ // for sphere template - __global__ void inside_sphere(const stim::sphere *V, size_t num, T Z,size_t *R, T *S, unsigned char *ptr, int x, int y, int z) { + __global__ void inside_sphere(const stim::sphere *V, size_t num, T Z,size_t *R, T *S, unsigned char *ptr, int x, int y, int z, T std = 1.0f) { unsigned ix = blockDim.x * blockIdx.x + threadIdx.x; unsigned iy = blockDim.y * blockIdx.y + threadIdx.y; @@ -93,11 +93,16 @@ namespace stim { } if (distance <= V[idx].r) ptr[(R[2] - 1 - iy) * R[0] * R[1] + ix * R[0]] = 255; + else { + T g = gaussianFunction(std::pow(distance - V[idx].r, 2), std); + if(g > 0.1f) + ptr[(R[2] - 1 - iy) * R[0] * R[1] + ix * R[0]] = 255; + } } // for cone template - __global__ void inside_cone(const stim::cone *E, size_t num, T Z, size_t *R, T *S, unsigned char *ptr, int x, int y, int z) { + __global__ void inside_cone(const stim::cone *E, size_t num, T Z, size_t *R, T *S, unsigned char *ptr, int x, int y, int z, T std = 1.0f) { unsigned ix = blockDim.x * blockIdx.x + threadIdx.x; unsigned iy = blockDim.y * blockIdx.y + threadIdx.y; @@ -129,6 +134,11 @@ namespace stim { } if (distance <= rr) ptr[(R[2] - 1 - iy) * R[0] * R[1] + ix * R[0]] = 255; + else { + T g = gaussianFunction(std::pow(distance - rr, 2), std); + if (g > 0.1f) + ptr[(R[2] - 1 - iy) * R[0] * R[1] + ix * R[0]] = 255; + } } // for source bus @@ -252,7 +262,7 @@ namespace stim { bool set = false; // flag indicates the pressure has been set std::vector P; // initial pressure - std::vector v; // velocity + std::vector v; // average velocity along each edge std::vector > main_feeder; // inlet/outlet main feeder std::vector pendant_vertex; std::vector > input; // first one store which vertex, second one stores which edge, third one stores in/out volume flow rate of that vertex @@ -350,12 +360,22 @@ namespace stim { return E[tmp_e].r(tmp_v); } - // get the radius of index j of edge i + // get the radius of point j on edge i T get_radius(unsigned i, unsigned j) { return E[i].r(j); } - // get the pendant vertices + // back up vertices + std::vector > back_vertex() { + std::vector > result; + + for (unsigned i = 0; i < num_vertex; i++) + result.push_back(stim::vec3(V[i][0], V[i][1], V[i][2])); + + return result; + } + + // get pendant vertices std::vector get_pendant_vertex() { std::vector result; int count = 0; @@ -851,7 +871,7 @@ namespace stim { } else glColor3f((float)color[i * 3 + 0] / 255, (float)color[i * 3 + 1] / 255, (float)color[i * 3 + 2] / 255); - + for (unsigned j = 0; j < E[i].size() - 1; j++) { // for every point on the edge center1 = E[i][j]; center2 = E[i][j + 1]; @@ -1532,13 +1552,13 @@ namespace stim { } // mark the vertex - void mark_vertex(T scale = 1.0f) { + void mark_vertex(std::vector > vertex, T scale = 1.0f) { glColor3f(0.0f, 0.0f, 0.0f); - for (unsigned i = 0; i < num_vertex; i++) { - glRasterPos3f(V[i][0], V[i][1] + get_radius(i) * scale, V[i][2]); + for (unsigned i = 0; i < vertex.size(); i++) { + glRasterPos3f(vertex[i][0], vertex[i][1] + get_radius(i) * scale, vertex[i][2]); // place position std::stringstream ss; - ss << i; + ss << i; // store index value glutBitmapString(GLUT_BITMAP_HELVETICA_18, (const unsigned char*)(ss.str().c_str())); } } @@ -2877,7 +2897,18 @@ namespace stim { f_file.close(); } - /// Calculate the inverse of A and store the result in C + /// check whether current network needs to be subdivided + bool check_subdivision() { + + for (size_t i = 0; i < num_edge; i++) { + if (E[i].size() > 2) { + return true; + } + } + return false; + } + + /// calculate the inverse of A and store the result in C void inversion(T** A, int order, T* C) { #ifdef __CUDACC__ @@ -2970,6 +3001,15 @@ namespace stim { free(minor); #endif } + + /// arithmetic + // assignment + flow & operator= (flow rhs) { + E = rhs.E; + V = rhs.V; + + return *this; + } }; } diff --git a/main.cu b/main.cu index e0af866..e5305f7 100644 --- a/main.cu +++ b/main.cu @@ -40,6 +40,7 @@ unsigned num_vertex; // number of vertex in the network std::vector pendant_vertex; // list of pendant vertex index in GT std::vector menu_option = { "simulation", "build inlet/outlet", "manufacture", "adjustment" }; stim::flow flow; // flow object +stim::flow backup; // flow backup float move_pace; // camera moving parameter float u; // viscosity float rou; // density @@ -79,7 +80,7 @@ int mouse_x; // window x-coordinate int mouse_y; // window y-coordinate int picked_x; // picked window x-coordinate int picked_y; // picked window y-coordinate -bool LTbutton = false; // true means down while false means up +bool LTbutton = false; // true means down while false means up // simulation parameters bool render_direction = false; // flag indicates rendering flow direction for one edge @@ -89,8 +90,10 @@ bool to_select_pressure = false; // flag indicates having selected a vertex to m bool mark_index = true; // flag indicates marking the index near the vertex bool glyph_mode = false; // flag indicates rendering glyph for flow velocity field bool frame_mode = false; // flag indicates rendering filament framing structrues +bool subdivided = false; // flag indicates subdivision status unsigned pressure_index; // the index of vertex that is clicked unsigned direction_index = UINT_MAX;// the index of edge that is pointed at +std::vector > back_vertex; // vertex back up for marking indices // build inlet/outlet parameters bool build_inlet_outlet = false; // flag indicates building inlets and outlets @@ -169,6 +172,8 @@ void flow_initialize() { flow.set = true; stim::vec3 center = bb.center(); + flow.P.clear(); + flow.P.resize(num_vertex, 0); // clear up initialized pressure for (unsigned i = 0; i < pendant_vertex.size(); i++) { if (flow.get_vertex(pendant_vertex[i])[0] <= center[0]) @@ -271,7 +276,7 @@ void glut_render() { if (!glyph_mode) { flow.glSolidSphere(max_pressure, subdivision, scale); if (mark_index) - flow.mark_vertex(scale); + flow.mark_vertex(back_vertex, scale); //flow.glSolidCone(subdivision); flow.glSolidCylinder(direction_index, color, subdivision, scale); } @@ -288,7 +293,7 @@ void glut_render() { if (!glyph_mode) { flow.glSolidSphere(max_pressure, subdivision, scale); if (mark_index) { - flow.mark_vertex(scale); + flow.mark_vertex(back_vertex, scale); //flow.mark_edge(); } //flow.glSolidCone(subdivision); @@ -299,7 +304,7 @@ void glut_render() { flow.glyph(color, subdivision, scale, frame_mode); } flow.glSolidCuboid(subdivision, manufacture, length); // render bus source - if (render_direction) // render the flow direction of the vertex pointed + if (render_direction && !glyph_mode) // render the flow direction of the vertex pointed flow.glSolidCone(direction_index, subdivision, scale); } @@ -495,6 +500,7 @@ void glut_menu(int value) { // first time if (!flow.set) { // only first time simulation called "simulation", ^_^ get_background(); // get the graph information + back_vertex = flow.back_vertex(); // vertex back up for marking indices flow_initialize(); // initialize flow condition menu_option[0] = "resimulation"; } @@ -584,6 +590,7 @@ void glut_passive_motion(int x, int y) { render_direction = false; direction_index = -1; } + undo = true; } if (mods == GLUT_ACTIVE_SHIFT && picked_connection) { @@ -888,7 +895,7 @@ void glut_keyboard(unsigned char key, int x, int y) { // register different keyboard operation switch (key) { - + // saving network flow profile case 's': flow.save_network(); @@ -903,7 +910,39 @@ void glut_keyboard(unsigned char key, int x, int y) { flow.saveNwt(filename); break; } - + + // subdivide current network for more detailed calculation + case 'd': { + // subdivide current network due to the limitation of current computation if needed + if (!subdivided && simulation && !glyph_mode) { + subdivided = true; + + // check whether current network can be subdivided + if (flow.check_subdivision()) { + flow.subdivision(); + get_background(); + } + + flow_initialize(); // initialize flow condition + // resimulation + flow_stable_state(); // main function of solving the linear system + flow.print_flow(); + undo = true; + } + else if (subdivided && simulation && !glyph_mode) { + subdivided = false; + flow = backup; // load back up + + get_background(); + flow_initialize(); + // resimulation + flow_stable_state(); // main function of solving the linear system + flow.print_flow(); + undo = true; + } + break; + } + // flow vector field visualization, Glyphs case 'f': if (glyph_mode && !manufacture && (simulation || build_inlet_outlet)) { @@ -923,7 +962,7 @@ void glut_keyboard(unsigned char key, int x, int y) { } undo = true; break; - + // filaments around arrows case 'g': if (glyph_mode) { @@ -932,8 +971,9 @@ void glut_keyboard(unsigned char key, int x, int y) { else frame_mode = true; } + undo = true; break; - + // open/close index marks case 'e': if (mark_index) @@ -1047,12 +1087,18 @@ int main(int argc, char* argv[]) { } else { // load network from user std::vector tmp = stim::parser::split(args.arg(0), '.'); - if ("obj" == tmp[1]) + if ("obj" == tmp[1]) { flow.load_obj(args.arg(0)); - else if ("nwt" == tmp[1]) // stim network binary format - flow.loadNwt(args.arg(0)); - else if ("swc" == tmp[1]) - flow.load_swc(args.arg(0)); + backup.load_obj(args.arg(0)); + } + else if ("nwt" == tmp[1]) { // stim network binary format + flow.loadNwt(args.arg(0)); + backup.loadNwt(args.arg(0)); + } + else if ("swc" == tmp[1]) { + flow.load_swc(args.arg(0)); + backup.load_swc(args.arg(0)); + } else { std::cout << "Invalid file type" << std::endl; std::exit(1); -- libgit2 0.21.4