fixed index marking bugs after network subdivided

Jiaming Guo
1 parent 4bb615eb
Showing 2 changed files with 110 additions and 24 deletions Show diff stats
flow.h
main.cu
@@ -67,7 +67,7 @@ namespace stim {
 #ifdef __CUDACC__
 	// for sphere
 	template <typename T>
-	__global__ void inside_sphere(const stim::sphere<T> *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<T> *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 <typename T>
-	__global__ void inside_cone(const stim::cone<T> *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<T> *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<T> P;														// initial pressure
-		std::vector<T> v;														// velocity
+		std::vector<T> v;														// average velocity along each edge
 		std::vector<typename stim::vec3<T> > main_feeder;						// inlet/outlet main feeder
 		std::vector<unsigned> pendant_vertex;
 		std::vector<typename stim::triple<unsigned, unsigned, T> > 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<stim::vec3<T> > back_vertex() {
+			std::vector<stim::vec3<T> > result;
+
+			for (unsigned i = 0; i < num_vertex; i++)
+				result.push_back(stim::vec3<T>(V[i][0], V[i][1], V[i][2]));
+
+			return result;
+		}
+
+		// get pendant vertices
 		std::vector<unsigned> get_pendant_vertex() {
 			std::vector<unsigned> 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<stim::vec3<T> > 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<T> & operator= (flow<T> rhs) {
+			E = rhs.E;
+			V = rhs.V;
+
+			return *this;
+		}
 	};
 }
@@ -40,6 +40,7 @@ unsigned num_vertex;									// number of vertex in the network
 std::vector<unsigned> pendant_vertex;					// list of pendant vertex index in GT
 std::vector<std::string> menu_option = { "simulation", "build inlet/outlet", "manufacture", "adjustment" };
 stim::flow<float> flow;									// flow object
+stim::flow<float> 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<stim::vec3<float> > 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<float> 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<std::string> 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);