Blame view

flow.h 99.1 KB
9191c39e   Jiaming Guo   first version of ...
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
  #ifndef FLOW3_H
  #define FLOW3_H
  
  #include <algorithm>
  
  //STIM include
  #include <stim/parser/arguments.h>
  #include <stim/visualization/gl_network.h>
  #include <stim/visualization/colormap.h>
  #include <stim/math/matrix.h>
  #include <stim/visualization/gl_aaboundingbox.h>
  #include <stim/ui/progressbar.h>
  #include <stim/grids/image_stack.h>
  
  #ifdef __CUDACC__
  #include <cublas_v2.h>
  #include <stim/cuda/cudatools/error.h>
  #endif
  
  namespace stim {
  	template <typename A, typename B, typename C>
  	struct triple {
  		A first;
  		B second;
  		C third;
  	};
  
  	template <typename T>
  	struct bridge {
  		std::vector<unsigned> v;				// vertices' indices
  		std::vector<typename stim::vec3<T> > V;	// vertices' coordinates
  		T l;		// length
8334680a   Jiaming Guo   add square wave c...
33
  		T r;		// radius
9191c39e   Jiaming Guo   first version of ...
34
35
36
37
38
39
40
  		T deltaP;	// pressure drop
  		T Q;		// volume flow rate
  	};
  
  	template <typename T>
  	struct sphere {
  		stim::vec3<T> c;		// center of sphere
8334680a   Jiaming Guo   add square wave c...
41
  		T r;					// radius
9191c39e   Jiaming Guo   first version of ...
42
43
44
  	};
  
  	template <typename T>
8334680a   Jiaming Guo   add square wave c...
45
  	struct cone {				// radius changes gradually
9191c39e   Jiaming Guo   first version of ...
46
47
  		stim::vec3<T> c1;		// center of geometry start hat
  		stim::vec3<T> c2;		// center of geometry end hat
8334680a   Jiaming Guo   add square wave c...
48
49
  		T r1;					// radius at start hat
  		T r2;					// radius at end hat
9191c39e   Jiaming Guo   first version of ...
50
51
52
53
54
55
56
57
58
59
  	};
  
  	template <typename T>
  	struct cuboid {
  		stim::vec3<T> c;
  		T l;					// length
  		T w;					// width
  		T h;					// height
  	};
  
b61addd8   Jiaming Guo   add adjustment fe...
60
61
62
63
64
65
  	template <typename T>
  	struct circuit {
  		std::vector<typename std::pair<unsigned, unsigned> > v;		// end vertex index
  		std::vector<T> r;											// branch resistence
  	};
  
9191c39e   Jiaming Guo   first version of ...
66
67
68
69
  	/// indicator function
  #ifdef __CUDACC__
  	// for sphere
  	template <typename T>
8f6c2057   Jiaming Guo   fixed index marki...
70
  	__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) {
9191c39e   Jiaming Guo   first version of ...
71
72
73
74
75
76
77
78
79
80
  
  		unsigned ix = blockDim.x * blockIdx.x + threadIdx.x;
  		unsigned iy = blockDim.y * blockIdx.y + threadIdx.y;
  
  		if (ix >= R[1] || iy >= R[2]) return;		// avoid seg-fault
  
  		// find world_pixel coordinates
  		stim::vec3<T> world_pixel;
  		world_pixel[0] = (T)ix * S[1] - x;			// translate origin to center of the network
  		world_pixel[1] = (T)iy * S[2] - y;
4bb615eb   Jiaming Guo   fixed generating ...
81
  		world_pixel[2] = ((T)z - Z / 2.0f) * S[3];	// ???center of box minus half width
9191c39e   Jiaming Guo   first version of ...
82
83
84
85
86
87
88
89
90
91
92
93
94
95
  
  		float distance = FLT_MAX;
  		float tmp_distance;
  		unsigned idx;
  
  		for (unsigned i = 0; i < num; i++) {
  			tmp_distance = (V[i].c - world_pixel).len();
  			if (tmp_distance <= distance) {
  				distance = tmp_distance;
  				idx = i;
  			}
  		}
  		if (distance <= V[idx].r)
  			ptr[(R[2] - 1 - iy) * R[0] * R[1] + ix * R[0]] = 255;
8f6c2057   Jiaming Guo   fixed index marki...
96
97
98
99
100
  		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;
  		}
9191c39e   Jiaming Guo   first version of ...
101
102
103
104
  	}
  
  	// for cone
  	template <typename T>
8f6c2057   Jiaming Guo   fixed index marki...
105
  	__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) {
9191c39e   Jiaming Guo   first version of ...
106
107
108
109
110
111
112
113
114
  
  		unsigned ix = blockDim.x * blockIdx.x + threadIdx.x;
  		unsigned iy = blockDim.y * blockIdx.y + threadIdx.y;
  
  		if (ix >= R[1] || iy >= R[2]) return;			// avoid segfault
  
  		stim::vec3<T> world_pixel;
  		world_pixel[0] = (T)ix * S[1] - x;
  		world_pixel[1] = (T)iy * S[2] - y;
4bb615eb   Jiaming Guo   fixed generating ...
115
  		world_pixel[2] = ((T)z - Z / 2.0f) * S[3];
9191c39e   Jiaming Guo   first version of ...
116
117
118
  
  		float distance = FLT_MAX;
  		float tmp_distance;
8334680a   Jiaming Guo   add square wave c...
119
  		float rr;										// radius at the surface where projection meets
9191c39e   Jiaming Guo   first version of ...
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
  
  		for (unsigned i = 0; i < num; i++) {			// find the nearest cylinder
  			tmp_distance = ((world_pixel - E[i].c1).cross(world_pixel - E[i].c2)).len() / (E[i].c2 - E[i].c1).len();
  			if (tmp_distance <= distance) {
  				// we only focus on point to line segment
  				// check to see whether projection is lying outside the line segment
  				float a = (world_pixel - E[i].c1).dot((E[i].c2 - E[i].c1).norm());
  				float b = (world_pixel - E[i].c2).dot((E[i].c1 - E[i].c2).norm());
  				float length = (E[i].c1 - E[i].c2).len();
  				if (a <= length && b <= length) {		// projection lying inside the line segment
  					distance = tmp_distance;
  					rr = E[i].r1 + (E[i].r2 - E[i].r1) * a / (length);		// linear change
  				}
  			}
  		}
  		if (distance <= rr)
  			ptr[(R[2] - 1 - iy) * R[0] * R[1] + ix * R[0]] = 255;
8f6c2057   Jiaming Guo   fixed index marki...
137
138
139
140
141
  		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;
  		}
9191c39e   Jiaming Guo   first version of ...
142
143
144
145
  	}
  
  	// for source bus
  	template <typename T>
4bb615eb   Jiaming Guo   fixed generating ...
146
  	__global__ void inside_cuboid(const stim::cuboid<T> *B, size_t num, T Z, size_t *R, T *S, unsigned char *ptr, int x, int y, int z) {
9191c39e   Jiaming Guo   first version of ...
147
148
149
150
151
152
153
154
155
  
  		unsigned ix = blockDim.x * blockIdx.x + threadIdx.x;
  		unsigned iy = blockDim.y * blockIdx.y + threadIdx.y;
  
  		if (ix >= R[1] || iy >= R[2]) return;			// avoid segfault
  
  		stim::vec3<T> world_pixel;
  		world_pixel[0] = (T)ix * S[1] - x;
  		world_pixel[1] = (T)iy * S[2] - y;
4bb615eb   Jiaming Guo   fixed generating ...
156
  		world_pixel[2] = ((T)z - Z / 2.0f) * S[3];
9191c39e   Jiaming Guo   first version of ...
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
  
  		for (unsigned i = 0; i < num; i++) {
  			bool left_outside = false;					// flag indicates point is outside the left bound
  			bool right_outside = false;
  
  			stim::vec3<T> tmp = B[i].c;
  			stim::vec3<T> L = stim::vec3<T>(tmp[0] - B[i].l / 2.0f, tmp[1] - B[i].h / 2.0f, tmp[2] - B[i].w / 2.0f);
  			stim::vec3<T> U = stim::vec3<T>(tmp[0] + B[i].l / 2.0f, tmp[1] + B[i].h / 2.0f, tmp[2] + B[i].w / 2.0f);
  
  			for (unsigned d = 0; d < 3; d++) {
  				if (world_pixel[d] < L[d])				// if the point is less than the minimum bound
  					left_outside = true;
  				if (world_pixel[d] > U[d])				// if the point is greater than the maximum bound
  					right_outside = true;
  			}
  			if (!left_outside && !right_outside)
  				ptr[(R[2] - 1 - iy) * R[0] * R[1] + ix * R[0]] = 255;
  		}
  	}
  #endif
  
  	template <typename T>
9b69bb4c   Jiaming Guo   add diplay list t...
179
  	class flow : public stim::gl_network<T> {
9191c39e   Jiaming Guo   first version of ...
180
181
  
  	private:
9b69bb4c   Jiaming Guo   add diplay list t...
182
  
9191c39e   Jiaming Guo   first version of ...
183
184
  		unsigned num_edge;
  		unsigned num_vertex;
9b69bb4c   Jiaming Guo   add diplay list t...
185
186
187
  		GLuint dlist;					// display list for inlets/outlets connections
  
  		enum direction { UP, LEFT, DOWN, RIGHT };
9191c39e   Jiaming Guo   first version of ...
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
  
  		// calculate the cofactor of elemen[row][col]
  		void get_minor(T** src, T** dest, int row, int col, int order) {
  
  			// index of element to be copied
  			int rowCount = 0;
  			int colCount = 0;
  
  			for (int i = 0; i < order; i++) {
  				if (i != row) {
  					colCount = 0;
  					for (int j = 0; j < order; j++) {
  						// when j is not the element
  						if (j != col) {
  							dest[rowCount][colCount] = src[i][j];
  							colCount++;
  						}
  					}
  					rowCount++;
  				}
  			}
  		}
  
  		// calculate the det()
  		T determinant(T** mat, int order) {
  
  			// degenate case when n = 1
  			if (order == 1)
  				return mat[0][0];
  
  			T det = 0.0;		// determinant value
  
  								// allocate the cofactor matrix
  			T** minor = (T**)malloc((order - 1) * sizeof(T*));
  			for (int i = 0; i < order - 1; i++)
  				minor[i] = (T*)malloc((order - 1) * sizeof(T));
  
  
  			for (int i = 0; i < order; i++) {
  
  				// get minor of element(0, i)
  				get_minor(mat, minor, 0, i, order);
  
  				// recursion
  				det += (i % 2 == 1 ? -1.0 : 1.0) * mat[0][i] * determinant(minor, order - 1);
  			}
  
  			// release memory
  			for (int i = 0; i < order - 1; i++)
  				free(minor[i]);
  			free(minor);
  
  			return det;
  		}
  
  	protected:
  
  		using stim::network<T>::E;
  		using stim::network<T>::V;
  		using stim::network<T>::get_start_vertex;
  		using stim::network<T>::get_end_vertex;
  		using stim::network<T>::get_r;
  		using stim::network<T>::get_average_r;
  		using stim::network<T>::get_l;
9b69bb4c   Jiaming Guo   add diplay list t...
252
  
9191c39e   Jiaming Guo   first version of ...
253
254
255
256
  		T** C;																	// Conductance
  		std::vector<typename stim::triple<unsigned, unsigned, float> > Q;		// volume flow rate
  		std::vector<T> QQ;														// Q' vector
  		std::vector<T> pressure;												// final pressure
f4105b89   Jiaming Guo   add new function:...
257
258
  		std::vector<typename std::vector<typename stim::vec3<T> > > in_backup;	// inlet connection back up
  		std::vector<typename std::vector<typename stim::vec3<T> > > out_backup;
b61addd8   Jiaming Guo   add adjustment fe...
259
  		std::string units;														// length units
9191c39e   Jiaming Guo   first version of ...
260
261
262
  
  	public:
  
8334680a   Jiaming Guo   add square wave c...
263
  		bool set = false;														// flag indicates the pressure has been set
9191c39e   Jiaming Guo   first version of ...
264
  		std::vector<T> P;														// initial pressure
8f6c2057   Jiaming Guo   fixed index marki...
265
  		std::vector<T> v;														// average velocity along each edge
9191c39e   Jiaming Guo   first version of ...
266
267
268
269
270
271
272
273
274
275
  		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
  		std::vector<typename stim::triple<unsigned, unsigned, T> > output;
  		std::vector<typename stim::bridge<T> > inlet;							// input bridge
  		std::vector<typename stim::bridge<T> > outlet;							// output bridge
  		std::vector<typename stim::sphere<T> > A;			// sphere model for making image stack
  		std::vector<typename stim::cone<T> > B;				// cone(cylinder) model for making image stack
  		std::vector<typename stim::cuboid<T> > CU;			// cuboid model for making image stack
  		stim::gl_aaboundingbox<T> bb;						// bounding box
f4105b89   Jiaming Guo   add new function:...
276
277
278
279
  		std::vector<bool> inlet_feasibility;				// list of flags indicate whether one inlet connection is feasible
  		std::vector<bool> outlet_feasibility;
  		std::vector<typename std::pair<stim::vec3<T>, stim::vec3<T> > > inbb;	// inlet connection bounding box
  		std::vector<typename std::pair<stim::vec3<T>, stim::vec3<T> > > outbb;	// outlet connection bounding box
b61addd8   Jiaming Guo   add adjustment fe...
280
281
  		T Ps;												// source and end pressure
  		T Pe;
9191c39e   Jiaming Guo   first version of ...
282
283
284
285
286
287
  
  		flow() {}				// default constructor
  		~flow() {
  			for (unsigned i = 0; i < num_vertex; i++)
  				delete[] C[i];
  			delete[] C;
b61addd8   Jiaming Guo   add adjustment fe...
288
  		}		// default destructor
9191c39e   Jiaming Guo   first version of ...
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
  
  		void init(unsigned n_e, unsigned n_v) {
  
  			num_edge = n_e;
  			num_vertex = n_v;
  
  			C = new T*[n_v]();
  			for (unsigned i = 0; i < n_v; i++) {
  				C[i] = new T[n_v]();
  			}
  
  			QQ.resize(n_v);
  			P.resize(n_v);
  			pressure.resize(n_v);
  
  			Q.resize(n_e);
  			v.resize(n_e);
  		}
  
  		void clear() {
  
  			for (unsigned i = 0; i < num_vertex; i++) {
  				QQ[i] = 0;
  				pressure[i] = 0;
  				for (unsigned j = 0; j < num_vertex; j++) {
  					C[i][j] = 0;
  				}
  			}
  			main_feeder.clear();
  			input.clear();
  			output.clear();
  			inlet.clear();
  			outlet.clear();
9b69bb4c   Jiaming Guo   add diplay list t...
322
323
324
325
326
  
  			if (glIsList(dlist)) {
  				glDeleteLists(dlist, 1);					// delete display list for modify
  				glDeleteLists(dlist + 1, 1);
  			}
9191c39e   Jiaming Guo   first version of ...
327
328
  		}
  
b61addd8   Jiaming Guo   add adjustment fe...
329
330
331
332
  		void set_units(std::string u) {
  			units = u;
  		}
  
9191c39e   Jiaming Guo   first version of ...
333
334
  		// copy radius from cylinder to flow
  		void set_radius(unsigned i, T radius) {
9b69bb4c   Jiaming Guo   add diplay list t...
335
  
9191c39e   Jiaming Guo   first version of ...
336
337
338
339
340
341
342
343
  			for (unsigned j = 0; j < num_edge; j++) {
  				if (E[j].v[0] == i)
  					E[j].cylinder<T>::set_r(0, radius);
  				else if (E[j].v[1] == i)
  					E[j].cylinder<T>::set_r(E[j].size() - 1, radius);
  			}
  		}
  
8334680a   Jiaming Guo   add square wave c...
344
  		// get the radius of vertex i
9191c39e   Jiaming Guo   first version of ...
345
  		T get_radius(unsigned i) {
9b69bb4c   Jiaming Guo   add diplay list t...
346
  
9191c39e   Jiaming Guo   first version of ...
347
348
349
350
351
352
353
354
355
  			unsigned tmp_e;				// edge index
  			unsigned tmp_v;				// vertex index in that edge
  			for (unsigned j = 0; j < num_edge; j++) {
  				if (E[j].v[0] == i) {
  					tmp_e = j;
  					tmp_v = 0;
  				}
  				else if (E[j].v[1] == i) {
  					tmp_e = j;
6332d1e7   Jiaming Guo   fixed warnings
356
  					tmp_v = (unsigned)(E[j].size() - 1);
9191c39e   Jiaming Guo   first version of ...
357
358
359
360
361
362
  				}
  			}
  
  			return E[tmp_e].r(tmp_v);
  		}
  
8f6c2057   Jiaming Guo   fixed index marki...
363
  		// get the radius of point j on edge i
9e60c6a7   Jiaming Guo   fixed minor bug i...
364
365
366
367
  		T get_radius(unsigned i, unsigned j) {
  			return E[i].r(j);
  		}
  
8f6c2057   Jiaming Guo   fixed index marki...
368
369
370
371
372
373
374
375
376
377
378
  		// 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
ce1a6f5e   Jiaming Guo   add functions for...
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
  		std::vector<unsigned> get_pendant_vertex() {
  			std::vector<unsigned> result;
  			int count = 0;
  
  			for (unsigned i = 0; i < V.size(); i++) {			// for every vertex
  				for (unsigned j = 0; j < E.size(); j++) {		// for every edge
  					if (i == E[j].v[0] || i == E[j].v[1])		// check whether current vertex terminates one edge
  						count++;
  				}
  				if (count == 1) 								// is pendant vertex
  					result.push_back(i);
  				count = 0;										// reset count
  			}
  
  			return result;
  		}
  
9191c39e   Jiaming Guo   first version of ...
396
397
  		// get the velocity of pendant vertex i
  		T get_velocity(unsigned i) {
9b69bb4c   Jiaming Guo   add diplay list t...
398
  
9191c39e   Jiaming Guo   first version of ...
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
  			unsigned tmp_e;				// edge index
  			for (unsigned j = 0; j < num_edge; j++) {
  				if (E[j].v[0] == i) {
  					tmp_e = j;
  					break;
  				}
  				else if (E[j].v[1] == i) {
  					tmp_e = j;
  					break;
  				}
  			}
  
  			return v[tmp_e];
  		}
  
  		// set pressure at specifi vertex
  		void set_pressure(unsigned i, T value) {
  			P[i] = value;
  		}
  
ce1a6f5e   Jiaming Guo   add functions for...
419
420
421
422
423
424
  		// extrct the largest connected component
  		void extract_lcc() {
  			
  
  		}
  
9191c39e   Jiaming Guo   first version of ...
425
426
427
428
429
430
431
  		// solve the linear system to get stable flow state
  		void solve_flow(T viscosity) {
  
  			// clear up last time simulation
  			clear();
  
  			// get the pendant vertex indices
ce1a6f5e   Jiaming Guo   add functions for...
432
  			pendant_vertex = get_pendant_vertex();
9191c39e   Jiaming Guo   first version of ...
433
434
435
436
437
438
439
440
  
  			// get bounding box
  			bb = (*this).boundingbox();
  
  			// set the conductance matrix of flow object
  			unsigned start_vertex = 0;
  			unsigned end_vertex = 0;
  			for (unsigned i = 0; i < num_edge; i++) {
6332d1e7   Jiaming Guo   fixed warnings
441
442
  				start_vertex = (unsigned)get_start_vertex(i);		// get the start vertex index of current edge
  				end_vertex = (unsigned)get_end_vertex(i);			// get the end vertex index of current edge
9191c39e   Jiaming Guo   first version of ...
443
  
b61addd8   Jiaming Guo   add adjustment fe...
444
  				C[start_vertex][end_vertex] = -((T)stim::PI * std::pow(get_average_r(i), 4)) / (8 * u * get_l(i));
9191c39e   Jiaming Guo   first version of ...
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
  
  				C[end_vertex][start_vertex] = C[start_vertex][end_vertex];
  			}
  			// set the diagonal to the negative sum of row element
  			float sum = 0.0;
  			for (unsigned i = 0; i < num_vertex; i++) {
  				for (unsigned j = 0; j < num_vertex; j++) {
  					sum += C[i][j];
  				}
  				C[i][i] = -sum;
  				sum = 0.0;
  			}
  
  			// get the Q' vector QQ
  			// matrix manipulation to zero out the conductance matrix as defined by the boundary values that were enterd
  			for (unsigned i = 0; i < num_vertex; i++) {
  				if (P[i] != 0) {			// for every dangle vertex
  					for (unsigned j = 0; j < num_vertex; j++) {
  						if (j == i) {
  							QQ[i] = C[i][i] * P[i];
  						}
  						else {
  							C[i][j] = 0;
  							QQ[j] = QQ[j] - C[j][i] * P[i];
  							C[j][i] = 0;
  						}
  					}
  				}
  			}
  
  			// get the inverse of conductance matrix
  			stim::matrix<float> _C(num_vertex, num_vertex);
  			inversion(C, num_vertex, _C.data());
9b69bb4c   Jiaming Guo   add diplay list t...
478
  
9191c39e   Jiaming Guo   first version of ...
479
480
481
482
483
484
485
486
  			// get the pressure in the network
  			for (unsigned i = 0; i < num_vertex; i++) {
  				for (unsigned j = 0; j < num_vertex; j++) {
  					pressure[i] += _C(i, j) * QQ[j];
  				}
  			}
  
  			// get the flow state from known pressure
b61addd8   Jiaming Guo   add adjustment fe...
487
488
489
  			T start_pressure = 0.0;
  			T end_pressure = 0.0;
  			T deltaP = 0.0;
9191c39e   Jiaming Guo   first version of ...
490
  			for (unsigned i = 0; i < num_edge; i++) {
6332d1e7   Jiaming Guo   fixed warnings
491
492
  				start_vertex = (unsigned)get_start_vertex(i);
  				end_vertex = (unsigned)get_end_vertex(i);
9191c39e   Jiaming Guo   first version of ...
493
494
495
496
497
498
  				start_pressure = pressure[start_vertex];		// get the start vertex pressure of current edge
  				end_pressure = pressure[end_vertex];			// get the end vertex pressure of current edge
  				deltaP = start_pressure - end_pressure;				// deltaP = Pa - Pb
  
  				Q[i].first = start_vertex;
  				Q[i].second = end_vertex;
9b69bb4c   Jiaming Guo   add diplay list t...
499
  
b61addd8   Jiaming Guo   add adjustment fe...
500
501
  				Q[i].third = ((T)stim::PI * std::pow(get_average_r(i), 4) * deltaP) / (8 * u * get_l(i));
  				v[i] = Q[i].third / ((T)stim::PI * std::pow(get_average_r(i), 2));
9191c39e   Jiaming Guo   first version of ...
502
503
504
505
506
  			}
  		}
  
  		// get the brewer color map based on velocity
  		void get_color_map(T& max_v, T& min_v, std::vector<unsigned char>& color, std::vector<unsigned> pendant_vertex) {
9b69bb4c   Jiaming Guo   add diplay list t...
507
  
6332d1e7   Jiaming Guo   fixed warnings
508
509
  			size_t num_edge = Q.size();
  			size_t num_vertex = QQ.size();
9191c39e   Jiaming Guo   first version of ...
510
511
512
513
514
515
516
517
518
  
  			// find the absolute maximum velocity and minimum velocity
  			std::vector<float> abs_V(num_edge);
  			for (unsigned i = 0; i < num_edge; i++) {
  				abs_V[i] = std::fabsf(v[i]);
  			}
  
  			max_v = *std::max_element(abs_V.begin(), abs_V.end());
  			min_v = *std::min_element(abs_V.begin(), abs_V.end());
9b69bb4c   Jiaming Guo   add diplay list t...
519
  
9191c39e   Jiaming Guo   first version of ...
520
521
522
  			// get the color map based on velocity range along the network
  			color.clear();
  			if (pendant_vertex.size() == 2 && num_edge - num_vertex + 1 <= 0) 		// only one inlet and one outlet
6ea69c93   Jiaming Guo   change bus size
523
  				color.resize(num_edge * 3, (unsigned char)128);
9191c39e   Jiaming Guo   first version of ...
524
525
526
527
528
529
530
531
532
533
  			else {
  				color.resize(num_edge * 3);
  				stim::cpu2cpu<float>(&abs_V[0], &color[0], num_edge, min_v, max_v, stim::cmBrewer);
  			}
  		}
  
  		// print flow
  		void print_flow() {
  
  			// show the pressure information in console box
b61addd8   Jiaming Guo   add adjustment fe...
534
  			std::cout << "PRESSURE(g/" << units << "/s^2):" << std::endl;
9191c39e   Jiaming Guo   first version of ...
535
536
537
538
  			for (unsigned i = 0; i < num_vertex; i++) {
  				std::cout << "[" << i << "] " << pressure[i] << std::endl;
  			}
  			// show the flow rate information in console box
b61addd8   Jiaming Guo   add adjustment fe...
539
  			std::cout << "VOLUME FLOW RATE(" << units << "^3/s):" << std::endl;
9191c39e   Jiaming Guo   first version of ...
540
541
542
543
544
545
546
547
548
549
  			for (unsigned i = 0; i < num_edge; i++) {
  				std::cout << "(" << Q[i].first << "," << Q[i].second << ")" << Q[i].third << std::endl;
  			}
  		}
  
  		/// helper function
  		// find hilbert curve order
  		// @param: current direct length between two vertices
  		// @param: desire length
  		void find_hilbert_order(T l, T d, int &order) {
9b69bb4c   Jiaming Guo   add diplay list t...
550
  
9191c39e   Jiaming Guo   first version of ...
551
552
553
554
555
556
  			bool flag = false;
  			int o = 1;
  			T tmp;					// temp of length
  			while (!flag) {
  				// convert from cartesian length to hilbert length
  				// l -> l * (4 ^ order - 1)/(2 ^ order - 1)
6332d1e7   Jiaming Guo   fixed warnings
557
  				tmp = (T)(l * (std::pow(4, o) - 1) / (std::pow(2, o) - 1));
9191c39e   Jiaming Guo   first version of ...
558
559
560
561
562
563
564
565
  				if (tmp >= d)
  					flag = true;
  				else
  					o++;
  			}
  			order = o;
  		}
  
9e60c6a7   Jiaming Guo   fixed minor bug i...
566
  		// move hilbert curves
9191c39e   Jiaming Guo   first version of ...
567
  		void move(unsigned i, T *c, direction dir, T dl, int feeder, bool invert) {
9b69bb4c   Jiaming Guo   add diplay list t...
568
  
9191c39e   Jiaming Guo   first version of ...
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
  			int cof = (invert) ? -1 : 1;
  
  			switch (dir) {
  			case UP:
  				c[1] += dl;
  				break;
  			case LEFT:
  				c[0] -= cof * dl;
  				break;
  			case DOWN:
  				c[1] -= dl;
  				break;
  			case RIGHT:
  				c[0] += cof * dl;
  				break;
  			}
  
  			stim::vec3<T> tmp;
  			for (unsigned i = 0; i < 3; i++)
  				tmp[i] = c[i];
9b69bb4c   Jiaming Guo   add diplay list t...
589
  
9191c39e   Jiaming Guo   first version of ...
590
591
592
593
594
595
  			if (feeder == 1)					// inlet main feeder
  				inlet[i].V.push_back(tmp);
  			else if (feeder == 0)				// outlet main feeder
  				outlet[i].V.push_back(tmp);
  		}
  
9e60c6a7   Jiaming Guo   fixed minor bug i...
596
  		// form hilbert curves
9191c39e   Jiaming Guo   first version of ...
597
  		void hilbert_curve(unsigned i, T *c, int order, T dl, int feeder, bool invert, direction dir = DOWN) {
9b69bb4c   Jiaming Guo   add diplay list t...
598
599
  
  			if (order == 1) {
9191c39e   Jiaming Guo   first version of ...
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
  				switch (dir) {
  				case UP:
  					move(i, c, DOWN, dl, feeder, invert);
  					move(i, c, RIGHT, dl, feeder, invert);
  					move(i, c, UP, dl, feeder, invert);
  					break;
  				case LEFT:
  					move(i, c, RIGHT, dl, feeder, invert);
  					move(i, c, DOWN, dl, feeder, invert);
  					move(i, c, LEFT, dl, feeder, invert);
  					break;
  				case DOWN:
  					move(i, c, UP, dl, feeder, invert);
  					move(i, c, LEFT, dl, feeder, invert);
  					move(i, c, DOWN, dl, feeder, invert);
  					break;
  				case RIGHT:
  					move(i, c, LEFT, dl, feeder, invert);
  					move(i, c, UP, dl, feeder, invert);
  					move(i, c, RIGHT, dl, feeder, invert);
  					break;
  				}
9b69bb4c   Jiaming Guo   add diplay list t...
622
  
9191c39e   Jiaming Guo   first version of ...
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
  			}
  			else if (order > 1) {
  				switch (dir) {
  				case UP:
  					hilbert_curve(i, c, order - 1, dl, feeder, invert, LEFT);
  					move(i, c, DOWN, dl, feeder, invert);
  					hilbert_curve(i, c, order - 1, dl, feeder, invert, UP);
  					move(i, c, RIGHT, dl, feeder, invert);
  					hilbert_curve(i, c, order - 1, dl, feeder, invert, UP);
  					move(i, c, UP, dl, feeder, invert);
  					hilbert_curve(i, c, order - 1, dl, feeder, invert, RIGHT);
  					break;
  				case LEFT:
  					hilbert_curve(i, c, order - 1, dl, feeder, invert, UP);
  					move(i, c, RIGHT, dl, feeder, invert);
  					hilbert_curve(i, c, order - 1, dl, feeder, invert, LEFT);
  					move(i, c, DOWN, dl, feeder, invert);
  					hilbert_curve(i, c, order - 1, dl, feeder, invert, LEFT);
  					move(i, c, LEFT, dl, feeder, invert);
  					hilbert_curve(i, c, order - 1, dl, feeder, invert, DOWN);
  					break;
  				case DOWN:
  					hilbert_curve(i, c, order - 1, dl, feeder, invert, RIGHT);
  					move(i, c, UP, dl, feeder, invert);
  					hilbert_curve(i, c, order - 1, dl, feeder, invert, DOWN);
  					move(i, c, LEFT, dl, feeder, invert);
  					hilbert_curve(i, c, order - 1, dl, feeder, invert, DOWN);
  					move(i, c, DOWN, dl, feeder, invert);
  					hilbert_curve(i, c, order - 1, dl, feeder, invert, LEFT);
  					break;
  				case RIGHT:
  					hilbert_curve(i, c, order - 1, dl, feeder, invert, DOWN);
  					move(i, c, LEFT, dl, feeder, invert);
  					hilbert_curve(i, c, order - 1, dl, feeder, invert, RIGHT);
  					move(i, c, UP, dl, feeder, invert);
  					hilbert_curve(i, c, order - 1, dl, feeder, invert, RIGHT);
  					move(i, c, RIGHT, dl, feeder, invert);
  					hilbert_curve(i, c, order - 1, dl, feeder, invert, UP);
  					break;
  				}
  			}
  		}
  
9191c39e   Jiaming Guo   first version of ...
666
667
668
669
  		/// render function
  		// find two envelope caps for two spheres
  		// @param cp1, cp2: list of points on the cap
  		// @param center1, center2: center point of cap
8334680a   Jiaming Guo   add square wave c...
670
  		// @param r1, r2: radius of cap
9191c39e   Jiaming Guo   first version of ...
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
  		void find_envelope(std::vector<typename stim::vec3<float> > &cp1, std::vector<typename stim::vec3<float> > &cp2, stim::vec3<float> center1, stim::vec3<float> center2, float r1, float r2, GLint subdivision) {
  
  			stim::vec3<float> tmp_d;
  			if (r1 == r2) {						// two vertices have the same radius
  				tmp_d = center2 - center1;		// calculate the direction vector
  				tmp_d = tmp_d.norm();
  				stim::circle<float> tmp_c;		// in order to get zero direction vector
  				tmp_c.rotate(tmp_d);
  
  				stim::circle<float> c1(center1, r1, tmp_d, tmp_c.U);
  				stim::circle<float> c2(center2, r2, tmp_d, tmp_c.U);
  				cp1 = c1.glpoints(subdivision);
  				cp2 = c2.glpoints(subdivision);
  			}
  			else {
  				if (r1 < r2) {					// switch index, we always want r1 to be larger than r2
  					stim::vec3<float> tmp_c = center2;
  					center2 = center1;
  					center1 = tmp_c;
  					float tmp_r = r2;
  					r2 = r1;
  					r1 = tmp_r;
  				}
  				tmp_d = center2 - center1;		// bigger one points to smaller one
  				tmp_d = tmp_d.norm();
  
  				float D = (center1 - center2).len();
  				stim::vec3<float> exp;
  				exp[0] = (center2[0] * r1 - center1[0] * r2) / (r1 - r2);
  				exp[1] = (center2[1] * r1 - center1[1] * r2) / (r1 - r2);
  
  				stim::vec3<float> t1, t2, t3, t4;
  				t1[2] = t2[2] = center1[2];		// decide the specific plane to work on
  				t3[2] = t4[2] = center2[2];
  
  				// first two
  				t1[0] = pow(r1, 2)*(exp[0] - center1[0]);
  				t1[0] += r1*(exp[1] - center1[1])*sqrt(pow((exp[0] - center1[0]), 2) + pow((exp[1] - center1[1]), 2) - pow(r1, 2));
  				t1[0] /= (pow((exp[0] - center1[0]), 2) + pow((exp[1] - center1[1]), 2));
  				t1[0] += center1[0];
  
  				t2[0] = pow(r1, 2)*(exp[0] - center1[0]);
  				t2[0] -= r1*(exp[1] - center1[1])*sqrt(pow((exp[0] - center1[0]), 2) + pow((exp[1] - center1[1]), 2) - pow(r1, 2));
  				t2[0] /= (pow((exp[0] - center1[0]), 2) + pow((exp[1] - center1[1]), 2));
  				t2[0] += center1[0];
  
  				t1[1] = pow(r1, 2)*(exp[1] - center1[1]);
  				t1[1] -= r1*(exp[0] - center1[0])*sqrt(pow((exp[0] - center1[0]), 2) + pow((exp[1] - center1[1]), 2) - pow(r1, 2));
  				t1[1] /= (pow((exp[0] - center1[0]), 2) + pow((exp[1] - center1[1]), 2));
  				t1[1] += center1[1];
  
  				t2[1] = pow(r1, 2)*(exp[1] - center1[1]);
  				t2[1] += r1*(exp[0] - center1[0])*sqrt(pow((exp[0] - center1[0]), 2) + pow((exp[1] - center1[1]), 2) - pow(r1, 2));
  				t2[1] /= (pow((exp[0] - center1[0]), 2) + pow((exp[1] - center1[1]), 2));
  				t2[1] += center1[1];
  
  				// check the correctness of the points
  				//float s = (center1[1] - t1[1])*(exp[1] - t1[1]) / ((t1[0] - center1[0])*(t1[0] - exp[0]));
  				//if (s != 1) {			// swap t1[1] and t2[1]
  				//	float tmp_t = t2[1];
  				//	t2[1] = t1[1];
  				//	t1[1] = tmp_t;
  				//}
  
  				// second two
  				t3[0] = pow(r2, 2)*(exp[0] - center2[0]);
  				t3[0] += r2*(exp[1] - center2[1])*sqrt(pow((exp[0] - center2[0]), 2) + pow((exp[1] - center2[1]), 2) - pow(r2, 2));
  				t3[0] /= (pow((exp[0] - center2[0]), 2) + pow((exp[1] - center2[1]), 2));
  				t3[0] += center2[0];
  
  				t4[0] = pow(r2, 2)*(exp[0] - center2[0]);
  				t4[0] -= r2*(exp[1] - center2[1])*sqrt(pow((exp[0] - center2[0]), 2) + pow((exp[1] - center2[1]), 2) - pow(r2, 2));
  				t4[0] /= (pow((exp[0] - center2[0]), 2) + pow((exp[1] - center2[1]), 2));
  				t4[0] += center2[0];
  
  				t3[1] = pow(r2, 2)*(exp[1] - center2[1]);
  				t3[1] -= r2*(exp[0] - center2[0])*sqrt(pow((exp[0] - center2[0]), 2) + pow((exp[1] - center2[1]), 2) - pow(r2, 2));
  				t3[1] /= (pow((exp[0] - center2[0]), 2) + pow((exp[1] - center2[1]), 2));
  				t3[1] += center2[1];
  
  				t4[1] = pow(r2, 2)*(exp[1] - center2[1]);
  				t4[1] += r2*(exp[0] - center2[0])*sqrt(pow((exp[0] - center2[0]), 2) + pow((exp[1] - center2[1]), 2) - pow(r2, 2));
  				t4[1] /= (pow((exp[0] - center2[0]), 2) + pow((exp[1] - center2[1]), 2));
  				t4[1] += center2[1];
  
  				// check the correctness of the points
  				//s = (center2[1] - t3[1])*(exp[1] - t3[1]) / ((t3[0] - center2[0])*(t3[0] - exp[0]));
  				//if (s != 1) {			// swap t1[1] and t2[1]
  				//	float tmp_t = t4[1];
  				//	t4[1] = t3[1];
  				//	t3[1] = tmp_t;
  				//}
  
  				stim::vec3<float> d1;
  				float dot;
  				float a;
  				float new_r;
  				stim::vec3<float> new_u;
  				stim::vec3<float> new_c;
  
  				// calculate the bigger circle
  				d1 = t1 - center1;
  				dot = d1.dot(tmp_d);
8334680a   Jiaming Guo   add square wave c...
774
  				a = dot / (r1 * 1) * r1;			// a = cos(alpha) * radius
9191c39e   Jiaming Guo   first version of ...
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
  				new_c = center1 + a * tmp_d;
  				new_r = sqrt(pow(r1, 2) - pow(a, 2));
  				new_u = t1 - new_c;
  
  				stim::circle<float> c1(new_c, new_r, tmp_d, new_u);
  				cp1 = c1.glpoints(subdivision);
  
  				// calculate the smaller circle
  				d1 = t3 - center2;
  				dot = d1.dot(tmp_d);
  				a = dot / (r2 * 1) * r2;
  				new_c = center2 + a * tmp_d;
  				new_r = sqrt(pow(r2, 2) - pow(a, 2));
  				new_u = t3 - new_c;
  
  				stim::circle<float> c2(new_c, new_r, tmp_d, new_u);
  				cp2 = c2.glpoints(subdivision);
  			}
  		}
  
  		// draw solid sphere at every vertex
ce1a6f5e   Jiaming Guo   add functions for...
796
  		void glSolidSphere(T max_pressure, GLint subdivision, T scale = 1.0f) {
9b69bb4c   Jiaming Guo   add diplay list t...
797
  
9191c39e   Jiaming Guo   first version of ...
798
799
  			// waste?
  			for (unsigned i = 0; i < num_edge; i++) {
9e60c6a7   Jiaming Guo   fixed minor bug i...
800
801
802
803
804
805
806
807
808
809
  				// draw the starting vertex
  				if (P[E[i].v[0]] != 0) {
  					stim::vec3<float> new_color;
  					new_color[0] = (P[E[i].v[0]] / max_pressure) > 0.5f ? 1.0f : 2.0f * P[E[i].v[0]] / max_pressure;						// red
  					new_color[1] = 0.0f;																									// green
  					new_color[2] = (P[E[i].v[0]] / max_pressure) > 0.5f ? 1.0f - 2.0f * (P[E[i].v[0]] / max_pressure - 0.5f) : 1.0f;		// blue
  					glColor3f(new_color[0], new_color[1], new_color[2]);
  
  					glPushMatrix();
  					glTranslatef(E[i][0][0], E[i][0][1], E[i][0][2]);
6ea69c93   Jiaming Guo   change bus size
810
  					glutSolidSphere(get_r(i, 0) * scale, subdivision, subdivision);
9e60c6a7   Jiaming Guo   fixed minor bug i...
811
812
  					glPopMatrix();
  				}
6ea69c93   Jiaming Guo   change bus size
813
814
815
816
817
818
819
820
821
822
823
824
  				else {
  					glEnable(GL_BLEND);											// enable color blend
  					glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);			// set blend function
  					glDisable(GL_DEPTH_TEST);
  					glColor4f(0.7f, 0.7f, 0.7f, 0.7f);							// gray color
  					glPushMatrix();
  					glTranslatef(E[i][0][0], E[i][0][1], E[i][0][2]);
  					glutSolidSphere(get_r(i, 0) * scale, subdivision, subdivision);
  					glPopMatrix();
  					glDisable(GL_BLEND);
  					glEnable(GL_DEPTH_TEST);
  				}
9e60c6a7   Jiaming Guo   fixed minor bug i...
825
826
827
828
829
830
831
832
  
  				// draw the ending vertex
  				if (P[E[i].v[1]] != 0) {
  					stim::vec3<float> new_color;
  					new_color[0] = (P[E[i].v[1]] / max_pressure) > 0.5f ? 1.0f : 2.0f * P[E[i].v[1]] / max_pressure;						// red
  					new_color[1] = 0.0f;																									// green
  					new_color[2] = (P[E[i].v[1]] / max_pressure) > 0.5f ? 1.0f - 2.0f * (P[E[i].v[1]] / max_pressure - 0.5f) : 1.0f;		// blue
  					glColor3f(new_color[0], new_color[1], new_color[2]);
9191c39e   Jiaming Guo   first version of ...
833
834
  
  					glPushMatrix();
9e60c6a7   Jiaming Guo   fixed minor bug i...
835
  					glTranslatef(E[i][E[i].size() - 1][0], E[i][E[i].size() - 1][1], E[i][E[i].size() - 1][2]);
6332d1e7   Jiaming Guo   fixed warnings
836
  					glutSolidSphere(get_r(i, (unsigned)(E[i].size() - 1)) * scale, subdivision, subdivision);
9191c39e   Jiaming Guo   first version of ...
837
838
  					glPopMatrix();
  				}
6ea69c93   Jiaming Guo   change bus size
839
840
841
842
843
844
845
  				else {
  					glEnable(GL_BLEND);											// enable color blend
  					glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);			// set blend function
  					glDisable(GL_DEPTH_TEST);
  					glColor4f(0.7f, 0.7f, 0.7f, 0.7f);							// gray color
  					glPushMatrix();
  					glTranslatef(E[i][E[i].size() - 1][0], E[i][E[i].size() - 1][1], E[i][E[i].size() - 1][2]);
6332d1e7   Jiaming Guo   fixed warnings
846
  					glutSolidSphere(get_r(i, (unsigned)(E[i].size() - 1)) * scale, subdivision, subdivision);
6ea69c93   Jiaming Guo   change bus size
847
848
849
850
  					glPopMatrix();
  					glDisable(GL_BLEND);
  					glEnable(GL_DEPTH_TEST);
  				}
9191c39e   Jiaming Guo   first version of ...
851
852
853
854
  			}
  		}
  
  		// draw edges as series of cylinders
ce1a6f5e   Jiaming Guo   add functions for...
855
  		void glSolidCylinder(unsigned index, std::vector<unsigned char> color, GLint subdivision, T scale = 1.0f) {
9191c39e   Jiaming Guo   first version of ...
856
857
858
859
  
  			stim::vec3<float> tmp_d;
  			stim::vec3<float> center1;
  			stim::vec3<float> center2;
9e60c6a7   Jiaming Guo   fixed minor bug i...
860
  			stim::circle<float> tmp_c;
9191c39e   Jiaming Guo   first version of ...
861
862
863
864
865
  			float r1;
  			float r2;
  			std::vector<typename stim::vec3<float> > cp1(subdivision + 1);
  			std::vector<typename stim::vec3<float> > cp2(subdivision + 1);
  			for (unsigned i = 0; i < num_edge; i++) {							// for every edge
9e60c6a7   Jiaming Guo   fixed minor bug i...
866
  				if (i == index) {												// render in tranparency for direction indication
f4105b89   Jiaming Guo   add new function:...
867
868
  					glEnable(GL_BLEND);											// enable color blend
  					glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);			// set blend function
9e60c6a7   Jiaming Guo   fixed minor bug i...
869
  					glDisable(GL_DEPTH_TEST);
6ea69c93   Jiaming Guo   change bus size
870
  					glColor4f((float)color[i * 3 + 0] / 255, (float)color[i * 3 + 1] / 255, (float)color[i * 3 + 2] / 255, 0.3f);
9e60c6a7   Jiaming Guo   fixed minor bug i...
871
872
873
  				}
  				else 
  					glColor3f((float)color[i * 3 + 0] / 255, (float)color[i * 3 + 1] / 255, (float)color[i * 3 + 2] / 255);
8f6c2057   Jiaming Guo   fixed index marki...
874
  
9191c39e   Jiaming Guo   first version of ...
875
876
877
878
  				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];
  
6ea69c93   Jiaming Guo   change bus size
879
880
  					r1 = get_r(i, j) * scale;
  					r2 = get_r(i, j + 1) * scale;
9b69bb4c   Jiaming Guo   add diplay list t...
881
  
9e60c6a7   Jiaming Guo   fixed minor bug i...
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
  					//// calculate the envelope caps
  					//find_envelope(cp1, cp2, center1, center2, r1, r2, subdivision);
  					if (j == 0) {
  						if (E[i].size() == 2)
  							find_envelope(cp1, cp2, center1, center2, r1, r2, subdivision);
  						else {
  							tmp_d = center2 - center1;
  							tmp_d = tmp_d.norm();
  							tmp_c.rotate(tmp_d);
  							stim::circle<float> c1(center1, r1, tmp_d, tmp_c.U);
  							cp1 = c1.glpoints(subdivision);
  							tmp_d = (E[i][j + 2] - center2) + (center2 - center1);
  							tmp_d = tmp_d.norm();
  							tmp_c.rotate(tmp_d);
  							stim::circle<float> c2(center2, r2, tmp_d, tmp_c.U);
  							cp2 = c2.glpoints(subdivision);
  						}
  					}
  					else if (j == E[i].size() - 2) {
  						tmp_d = (center2 - center1) + (center1 - E[i][j - 1]);
  						tmp_d = tmp_d.norm();
  						tmp_c.rotate(tmp_d);
  						stim::circle<float> c1(center1, r1, tmp_d, tmp_c.U);
  						cp1 = c1.glpoints(subdivision);
  						tmp_d = center2 - center1;
  						tmp_d = tmp_d.norm();
  						tmp_c.rotate(tmp_d);
  						stim::circle<float> c2(center2, r2, tmp_d, tmp_c.U);
  						cp2 = c2.glpoints(subdivision);
  					} 
  					else {
  						tmp_d = (center2 - center1) + (center1 - E[i][j - 1]);
  						tmp_d = tmp_d.norm();
  						tmp_c.rotate(tmp_d);
  						stim::circle<float> c1(center1, r1, tmp_d, tmp_c.U);
  						cp1 = c1.glpoints(subdivision);
  						tmp_d = (E[i][j + 2] - center2) + (center2 - center1);
  						tmp_d = tmp_d.norm();
  						tmp_c.rotate(tmp_d);
  						stim::circle<float> c2(center2, r2, tmp_d, tmp_c.U);
  						cp2 = c2.glpoints(subdivision);
  					}
9191c39e   Jiaming Guo   first version of ...
924
925
926
927
928
929
930
931
  
  					glBegin(GL_QUAD_STRIP);
  					for (unsigned j = 0; j < cp1.size(); j++) {
  						glVertex3f(cp1[j][0], cp1[j][1], cp1[j][2]);
  						glVertex3f(cp2[j][0], cp2[j][1], cp2[j][2]);
  					}
  					glEnd();
  				}
9e60c6a7   Jiaming Guo   fixed minor bug i...
932
933
934
935
936
  				if (i == index) {
  					glDisable(GL_BLEND);
  					glEnable(GL_DEPTH_TEST);
  				}
  					
9191c39e   Jiaming Guo   first version of ...
937
938
  			}
  			glFlush();
9191c39e   Jiaming Guo   first version of ...
939
940
  		}
  
6ea69c93   Jiaming Guo   change bus size
941
  		// draw the flow direction as cone, the size of the cone depends on the length of that edge
ce1a6f5e   Jiaming Guo   add functions for...
942
  		void glSolidCone(unsigned i, GLint subdivision, T scale = 1.0f, T threshold = 0.01f) {
9e60c6a7   Jiaming Guo   fixed minor bug i...
943
944
945
946
947
  
  			stim::vec3<T> tmp_d;									// direction
  			stim::vec3<T> center;									// cone hat center
  			stim::vec3<T> head;										// cone hat top
  			stim::circle<T> tmp_c;
6ea69c93   Jiaming Guo   change bus size
948
  			T h;													// height base of the cone
9e60c6a7   Jiaming Guo   fixed minor bug i...
949
950
951
  			std::vector<typename stim::vec3<T> > cp;
  			T radius;
  
6332d1e7   Jiaming Guo   fixed warnings
952
  			glColor3f(0.0f, 0.0f, 0.0f);							// lime color
9e60c6a7   Jiaming Guo   fixed minor bug i...
953
  
6332d1e7   Jiaming Guo   fixed warnings
954
  			size_t index = E[i].size() / 2 - 1;
9e60c6a7   Jiaming Guo   fixed minor bug i...
955
  			tmp_d = E[i][index + 1] - E[i][index];
6332d1e7   Jiaming Guo   fixed warnings
956
  			h = tmp_d.len() / 3.0f;									// get the height base by factor 3
9e60c6a7   Jiaming Guo   fixed minor bug i...
957
958
959
  			tmp_d = tmp_d.norm();
  			center = (E[i][index + 1] + E[i][index]) / 2;
  			tmp_c.rotate(tmp_d);
6332d1e7   Jiaming Guo   fixed warnings
960
961
  			radius = (E[i].r((unsigned)(index + 1)) + E[i].r((unsigned)index)) / 2 * scale;
  			radius = (T)((h / sqrt(3) < radius) ? h / sqrt(3) : radius);	// update radius
ce1a6f5e   Jiaming Guo   add functions for...
962
  			if (v[i] > threshold)
6ea69c93   Jiaming Guo   change bus size
963
  				head = center + tmp_d * h;
ce1a6f5e   Jiaming Guo   add functions for...
964
  			else if (v[i] < -threshold)
6ea69c93   Jiaming Guo   change bus size
965
  				head = center - tmp_d * h;
9e60c6a7   Jiaming Guo   fixed minor bug i...
966
967
968
969
  			
  			stim::circle<float> c(center, radius, tmp_d, tmp_c.U);
  			cp = c.glpoints(subdivision);
  
ce1a6f5e   Jiaming Guo   add functions for...
970
971
972
973
974
975
976
977
  			if (v[i] > threshold || v[i] < -threshold) {
  				glBegin(GL_TRIANGLE_FAN);
  				glVertex3f(head[0], head[1], head[2]);
  				for (unsigned k = 0; k < cp.size(); k++)
  					glVertex3f(cp[k][0], cp[k][1], cp[k][2]);
  				glEnd();
  				glFlush();
  			}
9e60c6a7   Jiaming Guo   fixed minor bug i...
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
  			// draw a cone for every edge to indicate 
  			//for (unsigned j = 0; j < E[i].size() - 1; j++) {	// for every point on current edge
  			//	tmp_d = E[i][j + 1] - E[i][j];
  			//	tmp_d = tmp_d.norm();
  			//	center = (E[i][j + 1] + E[i][j]) / 2;
  			//	tmp_c.rotate(tmp_d);
  			//	radius = (E[i].r(j + 1) + E[i].r(j)) / 2;
  			//	if (v[i] > 0)									// if flow flows from j to j+1
  			//		head = center + tmp_d * 2 * sqrt(3) * radius;
  			//	else
  			//		head = center - tmp_d * 2 * sqrt(3) * radius;
  
  			//	stim::circle<float> c(center, radius, tmp_d, tmp_c.U);
  			//	cp = c.glpoints(subdivision);
  
  			//	glBegin(GL_TRIANGLE_FAN);
  			//	glVertex3f(head[0], head[1], head[2]);
  			//	for (unsigned k = 0; k < cp.size(); k++)
  			//		glVertex3f(cp[k][0], cp[k][1], cp[k][2]);
  			//	glEnd();
  			//}
  			//glFlush();
  		}
ce1a6f5e   Jiaming Guo   add functions for...
1001
  		void glSolidCone(GLint subdivision, T scale = 1.0f, T threhold = 0.01f) {
9b69bb4c   Jiaming Guo   add diplay list t...
1002
  
9191c39e   Jiaming Guo   first version of ...
1003
1004
1005
1006
1007
  			stim::vec3<T> tmp_d;									// direction
  			stim::vec3<T> center;									// cone hat center
  			stim::vec3<T> head;										// cone hat top
  			stim::circle<T> tmp_c;
  			std::vector<typename stim::vec3<T> > cp;
6ea69c93   Jiaming Guo   change bus size
1008
  			T h;
9191c39e   Jiaming Guo   first version of ...
1009
1010
  			T radius;
  
6ea69c93   Jiaming Guo   change bus size
1011
  			glColor3f(0.600f, 0.847f, 0.788f);
9e60c6a7   Jiaming Guo   fixed minor bug i...
1012
  			// draw a cone for every edge to indicate 
9191c39e   Jiaming Guo   first version of ...
1013
  			for (unsigned i = 0; i < num_edge; i++) {				// for every edge
9e60c6a7   Jiaming Guo   fixed minor bug i...
1014
1015
1016
  				unsigned k1 = E[i].size() / 2 - 1;					// start and end index
  				unsigned k2 = E[i].size() / 2;
  				tmp_d = E[i][k2] - E[i][k1];
6ea69c93   Jiaming Guo   change bus size
1017
  				h = tmp_d.len() / 3.0f;								// get the height base by factor 3
9e60c6a7   Jiaming Guo   fixed minor bug i...
1018
1019
1020
  				tmp_d = tmp_d.norm();
  				center = (E[i][k2] + E[i][k1]) / 2;
  				tmp_c.rotate(tmp_d);
ce1a6f5e   Jiaming Guo   add functions for...
1021
  				radius = (E[i].r(k2) + E[i].r(k1)) / 2 * scale;
6ea69c93   Jiaming Guo   change bus size
1022
  				radius = (h / sqrt(3) < radius) ? h / sqrt(3) : radius;	// update radius by height base if necessary
ce1a6f5e   Jiaming Guo   add functions for...
1023
  				if (v[i] > threhold)										// if flow flows from k1 to k2
6ea69c93   Jiaming Guo   change bus size
1024
  					head = center + tmp_d * h;
ce1a6f5e   Jiaming Guo   add functions for...
1025
  				else if(v[i] < -threhold)
6ea69c93   Jiaming Guo   change bus size
1026
  					head = center - tmp_d * h;
9e60c6a7   Jiaming Guo   fixed minor bug i...
1027
1028
1029
  				stim::circle<float> c(center, radius, tmp_d, tmp_c.U);
  				cp = c.glpoints(subdivision);
  
ce1a6f5e   Jiaming Guo   add functions for...
1030
1031
1032
1033
1034
1035
1036
  				if (v[i] > threhold || v[i] < -threhold) {
  					glBegin(GL_TRIANGLE_FAN);
  					glVertex3f(head[0], head[1], head[2]);
  					for (unsigned k = 0; k < cp.size(); k++)
  						glVertex3f(cp[k][0], cp[k][1], cp[k][2]);
  					glEnd();
  				}
9e60c6a7   Jiaming Guo   fixed minor bug i...
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
  				//for (unsigned j = 0; j < E[i].size() - 1; j++) {	// for every point on current edge
  				//	tmp_d = E[i][j + 1] - E[i][j];
  				//	tmp_d = tmp_d.norm();
  				//	center = (E[i][j + 1] + E[i][j]) / 2;
  				//	tmp_c.rotate(tmp_d);
  				//	radius = (E[i].r(j + 1) + E[i].r(j)) / 2;
  				//	if (v[i] > 0)									// if flow flows from j to j+1
  				//		head = center + tmp_d * 2 * sqrt(3) * radius;
  				//	else
  				//		head = center - tmp_d * 2 * sqrt(3) * radius;
  
  				//	stim::circle<float> c(center, radius, tmp_d, tmp_c.U);
  				//	cp = c.glpoints(subdivision);
  
  				//	glBegin(GL_TRIANGLE_FAN);
  				//	glVertex3f(head[0], head[1], head[2]);
  				//	for (unsigned k = 0; k < cp.size(); k++)
  				//		glVertex3f(cp[k][0], cp[k][1], cp[k][2]);
  				//	glEnd();
  				//}
9191c39e   Jiaming Guo   first version of ...
1057
1058
1059
1060
1061
  			}
  			glFlush();
  		}
  
  		// draw main feeder as solid cube
ce1a6f5e   Jiaming Guo   add functions for...
1062
  		void glSolidCuboid(GLint subdivision, bool manufacture = false, T length = 40.0f, T height = 10.0f) {
9b69bb4c   Jiaming Guo   add diplay list t...
1063
  
9191c39e   Jiaming Guo   first version of ...
1064
  			T width;
4bb615eb   Jiaming Guo   fixed generating ...
1065
1066
1067
  			stim::gl_aaboundingbox<T> BB = (*this).boundingbox();
  			stim::vec3<T> L = BB.A;						// get the bottom left corner
  			stim::vec3<T> U = BB.B;						// get the top right corner
9191c39e   Jiaming Guo   first version of ...
1068
1069
  			width = U[2] - L[2] + 10.0f;
  
b61addd8   Jiaming Guo   add adjustment fe...
1070
1071
1072
1073
  			if (manufacture)
  				glColor3f(0.0f, 0.0f, 0.0f);			// black color
  			else
  				glColor3f(0.5f, 0.5f, 0.5f);			// gray color
9191c39e   Jiaming Guo   first version of ...
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
  			for (unsigned i = 0; i < main_feeder.size(); i++) {
  				// front face
  				glBegin(GL_QUADS);
  				glVertex3f(main_feeder[i][0] - length / 2, main_feeder[i][1] - height / 2, main_feeder[i][2] - width / 2);
  				glVertex3f(main_feeder[i][0] + length / 2, main_feeder[i][1] - height / 2, main_feeder[i][2] - width / 2);
  				glVertex3f(main_feeder[i][0] + length / 2, main_feeder[i][1] + height / 2, main_feeder[i][2] - width / 2);
  				glVertex3f(main_feeder[i][0] - length / 2, main_feeder[i][1] + height / 2, main_feeder[i][2] - width / 2);
  				glEnd();
  
  				// back face
  				glBegin(GL_QUADS);
  				glVertex3f(main_feeder[i][0] - length / 2, main_feeder[i][1] - height / 2, main_feeder[i][2] + width / 2);
  				glVertex3f(main_feeder[i][0] + length / 2, main_feeder[i][1] - height / 2, main_feeder[i][2] + width / 2);
  				glVertex3f(main_feeder[i][0] + length / 2, main_feeder[i][1] + height / 2, main_feeder[i][2] + width / 2);
  				glVertex3f(main_feeder[i][0] - length / 2, main_feeder[i][1] + height / 2, main_feeder[i][2] + width / 2);
  				glEnd();
  
  				// top face
  				glBegin(GL_QUADS);
  				glVertex3f(main_feeder[i][0] - length / 2, main_feeder[i][1] + height / 2, main_feeder[i][2] - width / 2);
  				glVertex3f(main_feeder[i][0] + length / 2, main_feeder[i][1] + height / 2, main_feeder[i][2] - width / 2);
  				glVertex3f(main_feeder[i][0] + length / 2, main_feeder[i][1] + height / 2, main_feeder[i][2] + width / 2);
  				glVertex3f(main_feeder[i][0] - length / 2, main_feeder[i][1] + height / 2, main_feeder[i][2] + width / 2);
  				glEnd();
  
  				// bottom face
  				glBegin(GL_QUADS);
  				glVertex3f(main_feeder[i][0] - length / 2, main_feeder[i][1] - height / 2, main_feeder[i][2] - width / 2);
  				glVertex3f(main_feeder[i][0] + length / 2, main_feeder[i][1] - height / 2, main_feeder[i][2] - width / 2);
  				glVertex3f(main_feeder[i][0] + length / 2, main_feeder[i][1] - height / 2, main_feeder[i][2] + width / 2);
  				glVertex3f(main_feeder[i][0] - length / 2, main_feeder[i][1] - height / 2, main_feeder[i][2] + width / 2);
  				glEnd();
  
  				// left face
  				glBegin(GL_QUADS);
  				glVertex3f(main_feeder[i][0] - length / 2, main_feeder[i][1] - height / 2, main_feeder[i][2] - width / 2);
  				glVertex3f(main_feeder[i][0] - length / 2, main_feeder[i][1] - height / 2, main_feeder[i][2] + width / 2);
  				glVertex3f(main_feeder[i][0] - length / 2, main_feeder[i][1] + height / 2, main_feeder[i][2] + width / 2);
  				glVertex3f(main_feeder[i][0] - length / 2, main_feeder[i][1] + height / 2, main_feeder[i][2] - width / 2);
  				glEnd();
  
  				// right face
  				glBegin(GL_QUADS);
  				glVertex3f(main_feeder[i][0] + length / 2, main_feeder[i][1] - height / 2, main_feeder[i][2] - width / 2);
  				glVertex3f(main_feeder[i][0] + length / 2, main_feeder[i][1] + height / 2, main_feeder[i][2] - width / 2);
  				glVertex3f(main_feeder[i][0] + length / 2, main_feeder[i][1] + height / 2, main_feeder[i][2] + width / 2);
  				glVertex3f(main_feeder[i][0] + length / 2, main_feeder[i][1] - height / 2, main_feeder[i][2] + width / 2);
  				glEnd();
  			}
9191c39e   Jiaming Guo   first version of ...
1123
1124
  		}
  
057ea93b   Jiaming Guo   add flow velocity...
1125
  		// draw flow velocity field, glyph
ce1a6f5e   Jiaming Guo   add functions for...
1126
  		void glyph(std::vector<unsigned char> color, GLint subdivision, T scale = 1.0f, bool frame = false, T r = 4.0f, T threshold = 0.01f) {
057ea93b   Jiaming Guo   add flow velocity...
1127
1128
1129
1130
1131
1132
1133
1134
1135
  			
  			// v1----v2-->v3
  			T k = 4.0f;							// quartering
  			stim::vec3<T> v1, v2, v3;			// three point
  			stim::vec3<T> d;					// direction vector
  			stim::circle<float> tmp_c;
  			std::vector<typename stim::vec3<float> > cp1(subdivision + 1);
  			std::vector<typename stim::vec3<float> > cp2(subdivision + 1);
  
ce1a6f5e   Jiaming Guo   add functions for...
1136
  			// rendering the arrows
057ea93b   Jiaming Guo   add flow velocity...
1137
1138
1139
1140
1141
  			for (unsigned i = 0; i < num_edge; i++) {				// for every edge
  				glColor3f((float)color[i * 3 + 0] / 255.0f, (float)color[i * 3 + 1] / 255.0f, (float)color[i * 3 + 2] / 255.0f);
  				for (unsigned j = 0; j < E[i].size() - 1; j++) {	// for every point on that edge
  
  					// consider the velocity valuence
ce1a6f5e   Jiaming Guo   add functions for...
1142
  					if (v[i] > threshold) {			// positive, from start point to end point
057ea93b   Jiaming Guo   add flow velocity...
1143
1144
1145
  						v1 = E[i][j];
  						v3 = E[i][j + 1];
  					}
ce1a6f5e   Jiaming Guo   add functions for...
1146
  					else if (v[i] < -threshold) {		// negative, from end point to start point
057ea93b   Jiaming Guo   add flow velocity...
1147
1148
1149
  						v1 = E[i][j + 1];
  						v3 = E[i][j];
  					}
ce1a6f5e   Jiaming Guo   add functions for...
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
  
  					if (v[i] > threshold || v[i] < -threshold) {
  						d = v3 - v1;
  						// place the arrow in the middel of one edge
  						v2 = v1 + (1.0f / k * 2.0f) * d;			// looks like =->=
  						v1 = v1 + (1.0f / k) * d;
  						v3 = v3 - (1.0f / k) * d;
  						d = d.norm();
  						tmp_c.rotate(d);
  
  						// render the cylinder part
  						stim::circle<T> c1(v1, r / 2 * scale, d, tmp_c.U);
  						cp1 = c1.glpoints(subdivision);
  						stim::circle<T> c2(v2, r / 2 * scale, d, tmp_c.U);
  						cp2 = c2.glpoints(subdivision);
  
  						glBegin(GL_QUAD_STRIP);
  						for (unsigned k = 0; k < cp1.size(); k++) {
  							glVertex3f(cp1[k][0], cp1[k][1], cp1[k][2]);
  							glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]);
  						}
  						glEnd();
  
  						// render the cone part
  						stim::circle<T> c3(v2, r * scale, d, tmp_c.U);
  						cp2 = c3.glpoints(subdivision);
  						glBegin(GL_TRIANGLE_FAN);
  						glVertex3f(v3[0], v3[1], v3[2]);
  						for (unsigned k = 0; k < cp2.size(); k++)
  							glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]);
  						glEnd();
057ea93b   Jiaming Guo   add flow velocity...
1181
  					}
ce1a6f5e   Jiaming Guo   add functions for...
1182
1183
  				}
  			}
057ea93b   Jiaming Guo   add flow velocity...
1184
  
ce1a6f5e   Jiaming Guo   add functions for...
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
  			// rendering frames
  			if (frame) {
  				frame = false;
  				stim::vec3<float> center1;
  				stim::vec3<float> center2;
  				stim::vec3<float> tmp_d;			// flow direction
  				float r1, r2;
  
  				for (unsigned i = 0; i < num_edge; i++) {
  					for (unsigned j = 0; j < E[i].size() - 1; j++) {
  						center1 = E[i][j];
  						center2 = E[i][j + 1];
  
  						r1 = get_r(i, j) * scale;
  						r2 = get_r(i, j + 1) * scale;
  
  						subdivision = 5;									// rough frames
  
  						if (j == 0) {
  							if (E[i].size() == 2)
  								find_envelope(cp1, cp2, center1, center2, r1, r2, subdivision);
  							else {
  								tmp_d = center2 - center1;
  								tmp_d = tmp_d.norm();
  								tmp_c.rotate(tmp_d);
  								stim::circle<float> c1(center1, r1, tmp_d, tmp_c.U);
  								cp1 = c1.glpoints(subdivision);
  								tmp_d = (E[i][j + 2] - center2) + (center2 - center1);
  								tmp_d = tmp_d.norm();
  								tmp_c.rotate(tmp_d);
  								stim::circle<float> c2(center2, r2, tmp_d, tmp_c.U);
  								cp2 = c2.glpoints(subdivision);
  							}
  						}
  						else if (j == E[i].size() - 2) {
  							tmp_d = (center2 - center1) + (center1 - E[i][j - 1]);
  							tmp_d = tmp_d.norm();
  							tmp_c.rotate(tmp_d);
  							stim::circle<float> c1(center1, r1, tmp_d, tmp_c.U);
  							cp1 = c1.glpoints(subdivision);
  							tmp_d = center2 - center1;
  							tmp_d = tmp_d.norm();
  							tmp_c.rotate(tmp_d);
  							stim::circle<float> c2(center2, r2, tmp_d, tmp_c.U);
  							cp2 = c2.glpoints(subdivision);
  						}
  						else {
  							tmp_d = (center2 - center1) + (center1 - E[i][j - 1]);
  							tmp_d = tmp_d.norm();
  							tmp_c.rotate(tmp_d);
  							stim::circle<float> c1(center1, r1, tmp_d, tmp_c.U);
  							cp1 = c1.glpoints(subdivision);
  							tmp_d = (E[i][j + 2] - center2) + (center2 - center1);
  							tmp_d = tmp_d.norm();
  							tmp_c.rotate(tmp_d);
  							stim::circle<float> c2(center2, r2, tmp_d, tmp_c.U);
  							cp2 = c2.glpoints(subdivision);
  						}
  
  						glColor3f(140/255.0f, 81/255.0f, 10/255.0f);
  						glBegin(GL_LINES);
  						for (unsigned k = 0; k < cp1.size(); k++) {
  							glVertex3f(cp1[k][0], cp1[k][1], cp1[k][2]);
  							glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]);
  						}
  						glEnd();
  					}
057ea93b   Jiaming Guo   add flow velocity...
1252
1253
  				}
  			}
057ea93b   Jiaming Guo   add flow velocity...
1254
1255
  		}
  
b61addd8   Jiaming Guo   add adjustment fe...
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
  		// display the total volume flow rate
  		void display_flow_rate(T in, T out) {
  			
  			glMatrixMode(GL_PROJECTION);									// set up the 2d viewport for mode text printing
  			glPushMatrix();
  			glLoadIdentity();
  			int X = glutGet(GLUT_WINDOW_WIDTH);								// get the current window width
  			int Y = glutGet(GLUT_WINDOW_HEIGHT);							// get the current window height
  			glViewport(0, 0, X, Y);											// locate to left bottom corner
  			gluOrtho2D(0, X, 0, Y);											// define othogonal aspect
  			glColor3f(0.8f, 0.0f, 0.0f);									// using red to show mode
  			glMatrixMode(GL_MODELVIEW);
  			glPushMatrix();
  			glLoadIdentity();
  
6332d1e7   Jiaming Guo   fixed warnings
1271
  			glRasterPos2f((GLfloat)(X / 2), (GLfloat)5.0f);					// hard coded position!!!!!
b61addd8   Jiaming Guo   add adjustment fe...
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
  			std::stringstream ss_p;
  			ss_p << "Q = ";				// Q = * um^3/s
  			ss_p << in;					
  			ss_p << " ";
  			ss_p << units;
  			ss_p << "^3/s";
  			glutBitmapString(GLUT_BITMAP_TIMES_ROMAN_24, (const unsigned char*)(ss_p.str().c_str()));
  
  			glPopMatrix();
  			glMatrixMode(GL_PROJECTION);
  			glPopMatrix();
  		}
  
0224d2ef   Jiaming Guo   fixed square wave...
1285
  		// draw the bridge as lines or arrows
ce1a6f5e   Jiaming Guo   add functions for...
1286
  		void line_bridge(bool &redisplay, T r = 4.0f) {
f4105b89   Jiaming Guo   add new function:...
1287
  
ce1a6f5e   Jiaming Guo   add functions for...
1288
  			if (redisplay) {							// check to see whether the display list needs to be updated
f4105b89   Jiaming Guo   add new function:...
1289
  				glDeleteLists(dlist, 1);
ce1a6f5e   Jiaming Guo   add functions for...
1290
1291
  				redisplay = false;
  			}
9b69bb4c   Jiaming Guo   add diplay list t...
1292
1293
1294
1295
  
  			if (!glIsList(dlist)) {
  				dlist = glGenLists(1);
  				glNewList(dlist, GL_COMPILE);
f4105b89   Jiaming Guo   add new function:...
1296
  
ce1a6f5e   Jiaming Guo   add functions for...
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
  				//// render flow direction arrows
  				//if (arrow) {
  				//	// v1----v2-->v3
  				//	T k = 4.0f;							// quartering
  				//	stim::vec3<T> v1, v2, v3;			// three point
  				//	stim::vec3<T> d;					// direction vector
  				//	stim::circle<float> tmp_c;
  				//	std::vector<typename stim::vec3<float> > cp1(subdivision + 1);
  				//	std::vector<typename stim::vec3<float> > cp2(subdivision + 1);
  
  				//	// inlet, right-going
  				//	for (unsigned i = 0; i < inlet.size(); i++) {
  				//		if (inlet_feasibility[i])
  				//			glColor3f(0.0f, 0.0f, 0.0f);			// feasible -> black
  				//		else
  				//			glColor3f(1.0f, 0.0f, 0.0f);			// nonfeasible -> red
  				//		for (unsigned j = 0; j < inlet[i].V.size() - 1; j++) {
  				//			v1 = inlet[i].V[j];
  				//			v3 = inlet[i].V[j + 1];
  				//			d = v3 - v1;
  				//			// place the arrow in the middel of one edge
  				//			v2 = v1 + (1.0f / k * 2.0f) * d;			// looks like =->=
  				//			v1 = v1 + (1.0f / k) * d;
  				//			v3 = v3 - (1.0f / k) * d;
  				//			d = d.norm();
  				//			tmp_c.rotate(d);
  
  				//			// render the cylinder part
  				//			stim::circle<T> c1(v1, r / 2, d, tmp_c.U);
  				//			cp1 = c1.glpoints(subdivision);
  				//			stim::circle<T> c2(v2, r / 2, d, tmp_c.U);
  				//			cp2 = c2.glpoints(subdivision);
  
  				//			glBegin(GL_QUAD_STRIP);
  				//			for (unsigned k = 0; k < cp1.size(); k++) {
  				//				glVertex3f(cp1[k][0], cp1[k][1], cp1[k][2]);
  				//				glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]);
  				//			}
  				//			glEnd();
  
  				//			// render the cone part
  				//			stim::circle<T> c3(v2, r, d, tmp_c.U);
  				//			cp2 = c3.glpoints(subdivision);
  				//			glBegin(GL_TRIANGLE_FAN);
  				//			glVertex3f(v3[0], v3[1], v3[2]);
  				//			for (unsigned k = 0; k < cp2.size(); k++)
  				//				glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]);
  				//			glEnd();
  				//		}
  				//	}
  
  				//	// outlet, right-going
  				//	for (unsigned i = 0; i < outlet.size(); i++) {
  				//		if (outlet_feasibility[i])
  				//			glColor3f(0.0f, 0.0f, 0.0f);			// feasible -> black
  				//		else
  				//			glColor3f(1.0f, 0.0f, 0.0f);			// nonfeasible -> red
  				//		for (unsigned j = 0; j < outlet[i].V.size() - 1; j++) {
  				//			v1 = outlet[i].V[j + 1];
  				//			v3 = outlet[i].V[j];
  				//			d = v3 - v1;
  				//			// place the arrow in the middel of one edge
  				//			v2 = v1 + (1.0f / k * 2.0f) * d;			// looks like =->=
  				//			v1 = v1 + (1.0f / k) * d;
  				//			v3 = v3 - (1.0f / k) * d;
  				//			d = d.norm();
  				//			tmp_c.rotate(d);
  
  				//			// render the cylinder part
  				//			stim::circle<T> c1(v1, r / 2, d, tmp_c.U);
  				//			cp1 = c1.glpoints(subdivision);
  				//			stim::circle<T> c2(v2, r / 2, d, tmp_c.U);
  				//			cp2 = c2.glpoints(subdivision);
  
  				//			glBegin(GL_QUAD_STRIP);
  				//			for (unsigned k = 0; k < cp1.size(); k++) {
  				//				glVertex3f(cp1[k][0], cp1[k][1], cp1[k][2]);
  				//				glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]);
  				//			}
  				//			glEnd();
  
  				//			// render the cone part
  				//			stim::circle<T> c3(v2, r, d, tmp_c.U);
  				//			cp2 = c3.glpoints(subdivision);
  				//			glBegin(GL_TRIANGLE_FAN);
  				//			glVertex3f(v3[0], v3[1], v3[2]);
  				//			for (unsigned k = 0; k < cp2.size(); k++)
  				//				glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]);
  				//			glEnd();
  				//		}
  				//	}
  
  				//	// render transparent lines as indexing
  				//	glEnable(GL_BLEND);											// enable color blend
  				//	glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);			// set blend function
  				//	glDisable(GL_DEPTH_TEST);
0224d2ef   Jiaming Guo   fixed square wave...
1393
1394
1395
  					glLineWidth(5);
  					for (unsigned i = 0; i < inlet.size(); i++) {
  						if (inlet_feasibility[i])
ce1a6f5e   Jiaming Guo   add functions for...
1396
1397
  							glColor3f(0.0f, 0.0f, 0.0f);
  				//			glColor4f(0.0f, 0.0f, 0.0f, 0.2f);		// feasible -> black
0224d2ef   Jiaming Guo   fixed square wave...
1398
  						else
ce1a6f5e   Jiaming Guo   add functions for...
1399
1400
  							glColor3f(1.0f, 0.0f, 0.0f);
  				//			glColor4f(1.0f, 0.0f, 0.0f, 0.2f);		// nonfeasible -> red
0224d2ef   Jiaming Guo   fixed square wave...
1401
1402
1403
1404
1405
1406
1407
1408
  
  						glBegin(GL_LINE_STRIP);
  						for (unsigned j = 0; j < inlet[i].V.size(); j++)
  							glVertex3f(inlet[i].V[j][0], inlet[i].V[j][1], inlet[i].V[j][2]);
  						glEnd();
  					}
  					for (unsigned i = 0; i < outlet.size(); i++) {
  						if (outlet_feasibility[i])
ce1a6f5e   Jiaming Guo   add functions for...
1409
1410
  							glColor3f(0.0f, 0.0f, 0.0f);
  				//			glColor4f(0.0f, 0.0f, 0.0f, 0.2f);		// feasible -> black
0224d2ef   Jiaming Guo   fixed square wave...
1411
  						else
ce1a6f5e   Jiaming Guo   add functions for...
1412
1413
  							glColor3f(1.0f, 0.0f, 0.0f);
  				//			glColor4f(1.0f, 0.0f, 0.0f, 0.2f);		// nonfeasible -> red
0224d2ef   Jiaming Guo   fixed square wave...
1414
1415
1416
1417
1418
  						glBegin(GL_LINE_STRIP);
  						for (unsigned j = 0; j < outlet[i].V.size(); j++)
  							glVertex3f(outlet[i].V[j][0], outlet[i].V[j][1], outlet[i].V[j][2]);
  						glEnd();
  					}
ce1a6f5e   Jiaming Guo   add functions for...
1419
1420
1421
  				//	glDisable(GL_BLEND);
  				//	glEnable(GL_DEPTH_TEST);
  				//}
9b69bb4c   Jiaming Guo   add diplay list t...
1422
  				glEndList();
9191c39e   Jiaming Guo   first version of ...
1423
  			}
9b69bb4c   Jiaming Guo   add diplay list t...
1424
  			glCallList(dlist);
9191c39e   Jiaming Guo   first version of ...
1425
1426
1427
  		}
  
  		// draw the bridge as tubes
6332d1e7   Jiaming Guo   fixed warnings
1428
  		void tube_bridge(bool &redisplay, GLint subdivision, T scale = 1.0f, T radius = 5.0f) {
ce1a6f5e   Jiaming Guo   add functions for...
1429
1430
1431
1432
1433
  
  			if (redisplay) {
  				glDeleteLists(dlist + 1, 1);
  				redisplay = false;
  			}
9191c39e   Jiaming Guo   first version of ...
1434
  
9b69bb4c   Jiaming Guo   add diplay list t...
1435
1436
1437
1438
1439
1440
1441
  			if (!glIsList(dlist + 1)) {
  				glNewList(dlist + 1, GL_COMPILE);
  
  				stim::vec3<T> dir;							// direction vector
  				stim::circle<T> unit_c;						// unit circle for finding the rotation start direction
  				std::vector<typename stim::vec3<T> > cp1;
  				std::vector<typename stim::vec3<T> > cp2;
8334680a   Jiaming Guo   add square wave c...
1442
  				glColor3f(0.0f, 0.0f, 0.0f);
9b69bb4c   Jiaming Guo   add diplay list t...
1443
1444
1445
1446
1447
1448
  
  				for (unsigned i = 0; i < inlet.size(); i++) {
  					// render vertex as sphere
  					for (unsigned j = 1; j < inlet[i].V.size() - 1; j++) {
  						glPushMatrix();
  						glTranslatef(inlet[i].V[j][0], inlet[i].V[j][1], inlet[i].V[j][2]);
ce1a6f5e   Jiaming Guo   add functions for...
1449
  						glutSolidSphere(radius * scale, subdivision, subdivision);
9b69bb4c   Jiaming Guo   add diplay list t...
1450
1451
1452
1453
1454
1455
1456
  						glPopMatrix();
  					}
  					// render edge as cylinder
  					for (unsigned j = 0; j < inlet[i].V.size() - 1; j++) {
  						dir = inlet[i].V[j] - inlet[i].V[j + 1];
  						dir = dir.norm();
  						unit_c.rotate(dir);
ce1a6f5e   Jiaming Guo   add functions for...
1457
1458
  						stim::circle<T> c1(inlet[i].V[j], inlet[i].r * scale, dir, unit_c.U);
  						stim::circle<T> c2(inlet[i].V[j + 1], inlet[i].r * scale, dir, unit_c.U);
9b69bb4c   Jiaming Guo   add diplay list t...
1459
1460
1461
1462
1463
1464
1465
1466
1467
  						cp1 = c1.glpoints(subdivision);
  						cp2 = c2.glpoints(subdivision);
  
  						glBegin(GL_QUAD_STRIP);
  						for (unsigned k = 0; k < cp1.size(); k++) {
  							glVertex3f(cp1[k][0], cp1[k][1], cp1[k][2]);
  							glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]);
  						}
  						glEnd();
9191c39e   Jiaming Guo   first version of ...
1468
  					}
9191c39e   Jiaming Guo   first version of ...
1469
  				}
9191c39e   Jiaming Guo   first version of ...
1470
  
9b69bb4c   Jiaming Guo   add diplay list t...
1471
1472
1473
1474
1475
  				for (unsigned i = 0; i < outlet.size(); i++) {
  					// render vertex as sphere
  					for (unsigned j = 1; j < outlet[i].V.size() - 1; j++) {
  						glPushMatrix();
  						glTranslatef(outlet[i].V[j][0], outlet[i].V[j][1], outlet[i].V[j][2]);
ce1a6f5e   Jiaming Guo   add functions for...
1476
  						glutSolidSphere(radius * scale, subdivision, subdivision);
9b69bb4c   Jiaming Guo   add diplay list t...
1477
1478
1479
1480
1481
1482
1483
  						glPopMatrix();
  					}
  					// render edge as cylinder
  					for (unsigned j = 0; j < outlet[i].V.size() - 1; j++) {
  						dir = outlet[i].V[j] - outlet[i].V[j + 1];
  						dir = dir.norm();
  						unit_c.rotate(dir);
ce1a6f5e   Jiaming Guo   add functions for...
1484
1485
  						stim::circle<T> c1(outlet[i].V[j], outlet[i].r * scale, dir, unit_c.U);
  						stim::circle<T> c2(outlet[i].V[j + 1], outlet[i].r * scale, dir, unit_c.U);
9b69bb4c   Jiaming Guo   add diplay list t...
1486
1487
1488
1489
1490
1491
1492
1493
1494
  						cp1 = c1.glpoints(subdivision);
  						cp2 = c2.glpoints(subdivision);
  
  						glBegin(GL_QUAD_STRIP);
  						for (unsigned k = 0; k < cp1.size(); k++) {
  							glVertex3f(cp1[k][0], cp1[k][1], cp1[k][2]);
  							glVertex3f(cp2[k][0], cp2[k][1], cp2[k][2]);
  						}
  						glEnd();
9191c39e   Jiaming Guo   first version of ...
1495
  					}
9191c39e   Jiaming Guo   first version of ...
1496
  				}
9b69bb4c   Jiaming Guo   add diplay list t...
1497
  				glEndList();
9191c39e   Jiaming Guo   first version of ...
1498
  			}
9b69bb4c   Jiaming Guo   add diplay list t...
1499
1500
  			glCallList(dlist + 1);
  		}	
9191c39e   Jiaming Guo   first version of ...
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
  
  		// draw gradient color bounding box outside the object
  		void bounding_box() {
  
  			stim::vec3<T> L = bb.A;						// get the bottom left corner
  			stim::vec3<T> U = bb.B;						// get the top right corner
  			
  			glLineWidth(1);
  			// front face of the box (in L[2])
  			glBegin(GL_LINE_LOOP);
  			glColor3f(0.0f, 0.0f, 0.0f);
  			glVertex3f(L[0], L[1], L[2]);
  			glColor3f(0.0f, 1.0f, 0.0f);
  			glVertex3f(L[0], U[1], L[2]);
  			glColor3f(1.0f, 1.0f, 0.0f);
  			glVertex3f(U[0], U[1], L[2]);
  			glColor3f(1.0f, 0.0f, 0.0f);
  			glVertex3f(U[0], L[1], L[2]);
  			glEnd();
  
  			// back face of the box (in U[2])
  			glBegin(GL_LINE_LOOP);
  			glColor3f(1.0f, 1.0f, 1.0f);
  			glVertex3f(U[0], U[1], U[2]);
  			glColor3f(0.0f, 1.0f, 1.0f);
  			glVertex3f(L[0], U[1], U[2]);
  			glColor3f(0.0f, 0.0f, 1.0f);
  			glVertex3f(L[0], L[1], U[2]);
  			glColor3f(1.0f, 0.0f, 1.0f);
  			glVertex3f(U[0], L[1], U[2]);
  			glEnd();
  
  			// fill out the rest of the lines to connect the two faces
  			glBegin(GL_LINES);
  			glColor3f(0.0f, 1.0f, 0.0f);
  			glVertex3f(L[0], U[1], L[2]);
  			glColor3f(0.0f, 1.0f, 1.0f);
  			glVertex3f(L[0], U[1], U[2]);
  			glColor3f(1.0f, 1.0f, 1.0f);
  			glVertex3f(U[0], U[1], U[2]);
  			glColor3f(1.0f, 1.0f, 0.0f);
  			glVertex3f(U[0], U[1], L[2]);
  			glColor3f(1.0f, 0.0f, 0.0f);
  			glVertex3f(U[0], L[1], L[2]);
  			glColor3f(1.0f, 0.0f, 1.0f);
  			glVertex3f(U[0], L[1], U[2]);
  			glColor3f(0.0f, 0.0f, 1.0f);
  			glVertex3f(L[0], L[1], U[2]);
  			glColor3f(0.0f, 0.0f, 0.0f);
  			glVertex3f(L[0], L[1], L[2]);
  			glEnd();
  		}
  
  		// mark the vertex
8f6c2057   Jiaming Guo   fixed index marki...
1555
  		void mark_vertex(std::vector<stim::vec3<T> > vertex, T scale = 1.0f) {
9191c39e   Jiaming Guo   first version of ...
1556
  			
8334680a   Jiaming Guo   add square wave c...
1557
  			glColor3f(0.0f, 0.0f, 0.0f);
8f6c2057   Jiaming Guo   fixed index marki...
1558
1559
  			for (unsigned i = 0; i < vertex.size(); i++) {
  				glRasterPos3f(vertex[i][0], vertex[i][1] + get_radius(i) * scale, vertex[i][2]);		// place position
9191c39e   Jiaming Guo   first version of ...
1560
  				std::stringstream ss;
8f6c2057   Jiaming Guo   fixed index marki...
1561
  				ss << i;																				// store index value
9191c39e   Jiaming Guo   first version of ...
1562
1563
1564
1565
  				glutBitmapString(GLUT_BITMAP_HELVETICA_18, (const unsigned char*)(ss.str().c_str()));
  			}
  		}
  
ce1a6f5e   Jiaming Guo   add functions for...
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
  		// mark the edge
  		void mark_edge() {
  
  			glColor3f(0.0f, 1.0f, 0.0f);
  			for (unsigned i = 0; i < num_edge; i++) {
  				glRasterPos3f((V[E[i].v[0]] + V[E[i].v[1]])[0]/2, (V[E[i].v[0]] + V[E[i].v[1]])[1] / 2, (V[E[i].v[0]] + V[E[i].v[1]])[2] / 2);
  				std::stringstream ss;
  				ss << i;
  				glutBitmapString(GLUT_BITMAP_HELVETICA_18, (const unsigned char*)(ss.str().c_str()));
  			}
  		}
  
9191c39e   Jiaming Guo   first version of ...
1578
1579
  		// find the nearest vertex of current click position
  		// return true and a value if found
6ea69c93   Jiaming Guo   change bus size
1580
  		inline bool epsilon_vertex(T x, T y, T z, T eps, T scale, unsigned& v) {
9191c39e   Jiaming Guo   first version of ...
1581
  
9e60c6a7   Jiaming Guo   fixed minor bug i...
1582
1583
  			T d = FLT_MAX;										// minimum distance between 2 vertices
  			T tmp_d = 0.0f;										// temporary stores distance for loop
9191c39e   Jiaming Guo   first version of ...
1584
  			unsigned tmp_i = 0;									// temporary stores connection index for loop
9e60c6a7   Jiaming Guo   fixed minor bug i...
1585
  			stim::vec3<T> tmp_v;								// temporary stores current loop point
9191c39e   Jiaming Guo   first version of ...
1586
1587
1588
  			d = FLT_MAX;										// set to max of float number
  
  			for (unsigned i = 0; i < V.size(); i++) {
9e60c6a7   Jiaming Guo   fixed minor bug i...
1589
  				tmp_v = stim::vec3<T>(x, y, z);
9191c39e   Jiaming Guo   first version of ...
1590
1591
1592
1593
1594
1595
1596
1597
  	
  				tmp_v = tmp_v - V[i];							// calculate a vector between two vertices
  				tmp_d = tmp_v.len();							// calculate length of that vector
  				if (tmp_d < d) {
  					d = tmp_d;									// if found a nearer vertex 
  					tmp_i = i;									// get the index of that vertex
  				}
  			}
6ea69c93   Jiaming Guo   change bus size
1598
  			eps += get_radius(tmp_i) * scale;					// increase epsilon accordingly
9191c39e   Jiaming Guo   first version of ...
1599
1600
1601
1602
1603
1604
1605
1606
  			if (d < eps) {										// if current click is close to any vertex
  				v = tmp_i;										// copy the extant vertex's index to v
  				return true;
  			}
  
  			return false;
  		}
  
9e60c6a7   Jiaming Guo   fixed minor bug i...
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
  		// find the nearest inlet/outlet connection line of current click position
  		// ab -> line segment, v -> point
  		// return true and a value if found
  		inline bool epsilon_edge(T x, T y, T z, T eps, unsigned &idx) {
  
  			T d = FLT_MAX;
  			T tmp_d;
  			unsigned tmp_i = 0;
  			unsigned tmp_j = 0;
  			stim::vec3<T> v1;
  			stim::vec3<T> v2;
  			stim::vec3<T> v3;
  			stim::vec3<T> v0 = stim::vec3<float>(x, y, z);
  			bool online = false;					// flag indicates the point is on the line-segment
  			float a, b;
  
f4105b89   Jiaming Guo   add new function:...
1623
  			// inner network
9e60c6a7   Jiaming Guo   fixed minor bug i...
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
  			for (unsigned i = 0; i < E.size(); i++) {
  				for (unsigned j = 0; j < E[i].size() - 1; j++) {
  					v1 = E[i][j + 1] - E[i][j];		// a -> b = ab
  					v2 = v0 - E[i][j];				// a -> v = av
  					v3 = v0 - E[i][j + 1];			// b -> v = bv
  
  					tmp_d = v2.dot(v1);				// av·ab
  
  					// check the line relative position
  					a = v2.dot(v1.norm());
  					b = v3.dot(v1.norm());
  					if (a < v1.len() && b < v1.len())		// if the length of projection fragment is longer than the line-segment
  						online = true;
  					else
  						online = false;
  
6ea69c93   Jiaming Guo   change bus size
1640
  					if (tmp_d <= 0.0 || tmp_d >= std::pow(v1.len(), 2) && !online)	// projection lies outside the line-segment
9e60c6a7   Jiaming Guo   fixed minor bug i...
1641
1642
1643
1644
1645
1646
1647
1648
1649
1650
1651
1652
  						continue;
  					else {
  						tmp_d = v1.cross(v2).len() / v1.len();						// perpendicular distance of point to segment: |v1 x v2| / |v1|
  						if (tmp_d < d) {
  							d = tmp_d;
  							tmp_i = i;
  							tmp_j = j;
  						}
  					} 
  				}
  			}
  
9e60c6a7   Jiaming Guo   fixed minor bug i...
1653
1654
1655
1656
1657
1658
1659
1660
1661
  			eps += get_radius(tmp_i, tmp_j);
  
  			if (d < eps) {
  				idx = tmp_i;
  				return true;
  			}
  
  			return false;
  		}
f4105b89   Jiaming Guo   add new function:...
1662
1663
1664
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694
1695
1696
1697
1698
1699
1700
1701
1702
1703
1704
1705
1706
1707
1708
1709
1710
1711
1712
1713
1714
1715
1716
1717
1718
1719
1720
1721
1722
1723
1724
1725
1726
1727
1728
1729
1730
1731
1732
1733
1734
1735
1736
1737
1738
1739
1740
1741
  		inline bool epsilon_edge(T x, T y, T z, T eps, unsigned &idx, unsigned &port) {
  
  			T d = FLT_MAX;
  			T tmp_d;
  			unsigned tmp_i = 0;
  			stim::vec3<T> v1;
  			stim::vec3<T> v2;
  			stim::vec3<T> v3;
  			stim::vec3<T> v0 = stim::vec3<float>(x, y, z);
  			bool online = false;					// flag indicates the point is on the line-segment
  			float a, b;
  
  			// inlet connection
  			for (unsigned i = 0; i < inlet.size(); i++) {
  				for (unsigned j = 0; j < inlet[i].V.size() - 1; j++) {
  					v1 = inlet[i].V[j + 1] - inlet[i].V[j];
  					v2 = v0 - inlet[i].V[j];
  					v3 = v0 - inlet[i].V[j + 1];
  
  					tmp_d = v2.dot(v1);				// av·ab
  
  					// check the line relative position
  					a = v2.dot(v1.norm());
  					b = v3.dot(v1.norm());
  					if (a < v1.len() && b < v1.len())		// if the length of projection fragment is longer than the line-segment
  						online = true;
  					else
  						online = false;
  
  					if (tmp_d <= 0.0 || tmp_d > std::pow(v1.len(), 2) && !online)	// projection lies outside the line-segment
  						continue;
  					else {
  						tmp_d = v1.cross(v2).len() / v1.len();						// perpendicular distance of point to segment: |v1 x v2| / |v1|
  						if (tmp_d < d) {
  							d = tmp_d;
  							tmp_i = i;
  							port = 0;
  						}
  					}
  				
  				}
  			}
  
  			// outlet connection
  			for (unsigned i = 0; i < outlet.size(); i++) {
  				for (unsigned j = 0; j < outlet[i].V.size() - 1; j++) {
  					v1 = outlet[i].V[j + 1] - outlet[i].V[j];
  					v2 = v0 - outlet[i].V[j];
  					v3 = v0 - outlet[i].V[j + 1];
  
  					tmp_d = v2.dot(v1);				// av·ab
  
  					// check the line relative position
  					a = v2.dot(v1.norm());
  					b = v3.dot(v1.norm());
  					if (a < v1.len() && b < v1.len())		// if the length of projection fragment is longer than the line-segment
  						online = true;
  					else
  						online = false;
  
  					if (tmp_d <= 0.0 || tmp_d > std::pow(v1.len(), 2) && !online)	// projection lies outside the line-segment
  						continue;
  					else {
  						tmp_d = v1.cross(v2).len() / v1.len();						// perpendicular distance of point to segment: |v1 x v2| / |v1|
  						if (tmp_d < d) {
  							d = tmp_d;
  							tmp_i = i;
  							port = 1;
  						}
  					}
  				}
  			}
  
  			if (d < eps) {
  				idx = tmp_i;
  				return true;
  			}
  
  			return false;
  		}
9e60c6a7   Jiaming Guo   fixed minor bug i...
1742
  
9191c39e   Jiaming Guo   first version of ...
1743
1744
  		/// build main feeder connection
  		// set up main feeder and main port of both input and output
15fc312f   Jiaming Guo   fixed modes conflict
1745
  		void set_main_feeder(T border = 120.0f) {
9191c39e   Jiaming Guo   first version of ...
1746
1747
1748
1749
1750
1751
1752
1753
  			
  			// 0 means outgoing while 1 means incoming
  			stim::vec3<T> inlet_main_feeder;
  			stim::vec3<T> outlet_main_feeder;
  
  			inlet_main_feeder = stim::vec3<T>(bb.A[0] - border, bb.center()[1], bb.center()[2]);
  			outlet_main_feeder = stim::vec3<T>(bb.B[0] + border, bb.center()[1], bb.center()[2]);
  			
b61addd8   Jiaming Guo   add adjustment fe...
1754
  			main_feeder.push_back(inlet_main_feeder);		// 0->inlet, 1->outlet
9191c39e   Jiaming Guo   first version of ...
1755
1756
1757
1758
  			main_feeder.push_back(outlet_main_feeder);
  
  			// find both input and output vertex
  			stim::triple<unsigned, unsigned, float> tmp;
6332d1e7   Jiaming Guo   fixed warnings
1759
  			size_t N = pendant_vertex.size();				// get the number of dangle vertex
9191c39e   Jiaming Guo   first version of ...
1760
1761
1762
1763
1764
1765
1766
1767
1768
1769
1770
1771
1772
1773
1774
1775
1776
1777
1778
1779
1780
1781
1782
1783
1784
1785
1786
1787
1788
1789
1790
1791
1792
1793
1794
1795
1796
1797
1798
1799
1800
1801
1802
1803
  			unsigned idx = 0;
  			for (unsigned i = 0; i < N; i++) {				// for every boundary vertex
  				idx = pendant_vertex[i];
  				for (unsigned j = 0; j < num_edge; j++) {	// for every edge
  					if (Q[j].first == idx) {			// starting vertex
  						if (Q[j].third > 0) {			// flow comes in
  							tmp.first = idx;
  							tmp.second = j;
  							tmp.third = Q[j].third;
  							input.push_back(tmp);
  							break;
  						}
  						// their might be a degenerate case that it equals to 0?
  						else if (Q[j].third < 0) {		// flow comes out
  							tmp.first = idx;
  							tmp.second = j;
  							tmp.third = -Q[j].third;
  							output.push_back(tmp);
  							break;
  						}
  					}
  					else if (Q[j].second == idx) {		// ending vertex
  						if (Q[j].third > 0) {			// flow comes in
  							tmp.first = idx;
  							tmp.second = j;
  							tmp.third = Q[j].third;
  							output.push_back(tmp);
  							break;
  						}
  						// their might be a degenerate case that it equals to 0?
  						else if (Q[j].third < 0) {		// flow comes out
  							tmp.first = idx;
  							tmp.second = j;
  							tmp.third = -Q[j].third;
  							input.push_back(tmp);
  							break;
  						}
  					}
  				}
  			}
  		}
  
  		// build connection between all inlets and outlets
  		// connection will trail along one axis around the bounding box
8334680a   Jiaming Guo   add square wave c...
1804
  		void build_synthetic_connection(T viscosity, T radius = 5.0f) {
9191c39e   Jiaming Guo   first version of ...
1805
1806
1807
1808
1809
1810
1811
1812
1813
  			
  			stim::vec3<T> L = bb.A;						// get the bottom left corner
  			stim::vec3<T> U = bb.B;						// get the top right corner
  			T box_length = U[0] - L[0];
  			T x0, dx;
  
  			stim::vec3<T> tmp_v;						// start vertex
  			stim::vec3<T> mid_v;						// middle point of the bridge
  			stim::vec3<T> bus_v;						// point on the bus
6ea69c93   Jiaming Guo   change bus size
1814
  			x0 = main_feeder[0][0] + 15.0f;				// assume bus length is 40.0f
9191c39e   Jiaming Guo   first version of ...
1815
1816
1817
  			for (unsigned i = 0; i < input.size(); i++) {
  				
  				tmp_v = V[input[i].first];
6ea69c93   Jiaming Guo   change bus size
1818
  				dx = 30.0f * ((tmp_v[0] - L[0]) / box_length);		// the socket position depends on proximity
9191c39e   Jiaming Guo   first version of ...
1819
1820
1821
1822
1823
1824
1825
1826
1827
1828
  				bus_v = stim::vec3<T>(x0 - dx, main_feeder[0][1], tmp_v[2]);
  				mid_v = stim::vec3<T>(x0 - dx, tmp_v[1], tmp_v[2]);
  
  				stim::bridge<T> tmp_b;
  				tmp_b.V.push_back(bus_v);
  				tmp_b.V.push_back(mid_v);
  				tmp_b.V.push_back(tmp_v);
  				tmp_b.v.push_back(input[i].first);
  				tmp_b.Q = input[i].third;
  				tmp_b.l = (bus_v - mid_v).len() + (mid_v - tmp_v).len();
8334680a   Jiaming Guo   add square wave c...
1829
  				tmp_b.r = radius;
9191c39e   Jiaming Guo   first version of ...
1830
1831
1832
1833
  
  				inlet.push_back(tmp_b);
  			}
  
6ea69c93   Jiaming Guo   change bus size
1834
  			x0 = main_feeder[1][0] - 15.0f;
9191c39e   Jiaming Guo   first version of ...
1835
1836
1837
  			for (unsigned i = 0; i < output.size(); i++) {
  
  				tmp_v = V[output[i].first];
6ea69c93   Jiaming Guo   change bus size
1838
  				dx = 30.0f * ((U[0] - tmp_v[0]) / box_length);		// the socket position depends on proximity
9191c39e   Jiaming Guo   first version of ...
1839
1840
1841
1842
1843
1844
1845
1846
1847
1848
  				bus_v = stim::vec3<T>(x0 + dx, main_feeder[1][1], tmp_v[2]);
  				mid_v = stim::vec3<T>(x0 + dx, tmp_v[1], tmp_v[2]);
  
  				stim::bridge<T> tmp_b;
  				tmp_b.V.push_back(bus_v);
  				tmp_b.V.push_back(mid_v);
  				tmp_b.V.push_back(tmp_v);
  				tmp_b.v.push_back(output[i].first);
  				tmp_b.Q = output[i].third;
  				tmp_b.l = (bus_v - mid_v).len() + (mid_v - tmp_v).len();
8334680a   Jiaming Guo   add square wave c...
1849
  				tmp_b.r = radius;
9191c39e   Jiaming Guo   first version of ...
1850
1851
1852
  
  				outlet.push_back(tmp_b);
  			}
f4105b89   Jiaming Guo   add new function:...
1853
1854
  
  			backup();
9191c39e   Jiaming Guo   first version of ...
1855
1856
  		}
  
8334680a   Jiaming Guo   add square wave c...
1857
1858
  		// find the number of U-shape or square-shape structure for extending length of connection
  		// @param t: width = t * radius
6332d1e7   Jiaming Guo   fixed warnings
1859
  		int find_number_square(T origin_l, T desire_l, int times = 10, T radius = 5.0f) {
8334680a   Jiaming Guo   add square wave c...
1860
1861
  			
  			bool done = false;						// flag indicates the current number of square shape structure is feasible
6332d1e7   Jiaming Guo   fixed warnings
1862
  			int n = (int)(origin_l / (times * 4 * radius));	// number of square shape structure
8334680a   Jiaming Guo   add square wave c...
1863
1864
1865
1866
1867
1868
1869
1870
1871
1872
1873
1874
1875
1876
1877
1878
1879
  			T need_l = desire_l - origin_l;
  			T height;								// height of the square shapce structure
  
  			while (!done) {
  				height = need_l / (2 * n);			// calculate the height
  				if (height > 2 * radius) {
  					done = true;
  				}
  				else {
  					n--;
  				}
  			}
  			
  			return n;
  		}
  
  		// build square connections
0224d2ef   Jiaming Guo   fixed square wave...
1880
  		void build_square_connection(int i, T width, T height, T origin_l, T desire_l, int n, int feeder, T threshold, bool z, bool left = true, bool up = true, int times = 10, T ratio = 0, T radius = 5.0f) {
8334680a   Jiaming Guo   add square wave c...
1881
1882
1883
1884
1885
1886
  		
  			int coef_up = (up) ? 1 : -1;				// y coefficient
  			int coef_left = (left) ? 1 : -1;			// x coefficient
  			int coef_z = (z) ? 1 : -1;					// z coefficient
  			int inverse = 1;							// inverse flag
  			stim::vec3<T> cor_v;						// corner vertex
f4105b89   Jiaming Guo   add new function:...
1887
  			std::pair<stim::vec3<T>, stim::vec3<T>> tmp_bb;
8334680a   Jiaming Guo   add square wave c...
1888
  			stim::vec3<T> tmp_v;
f4105b89   Jiaming Guo   add new function:...
1889
  			if (feeder == 1) 
70c0b942   Jiaming Guo   fixed bug when he...
1890
  				tmp_v = inlet[i].V[inlet[i].V.size() - 1];
f4105b89   Jiaming Guo   add new function:...
1891
  			else if (feeder == 0) 
70c0b942   Jiaming Guo   fixed bug when he...
1892
  				tmp_v = outlet[i].V[outlet[i].V.size() - 1];
f4105b89   Jiaming Guo   add new function:...
1893
  			tmp_bb.first = tmp_v;
8334680a   Jiaming Guo   add square wave c...
1894
  
b61addd8   Jiaming Guo   add adjustment fe...
1895
1896
1897
1898
1899
1900
1901
1902
1903
  			// pre-set fragments
  			if (ratio) {
  				T tmp_d, tmp_l;														// back ups
  				tmp_d = desire_l;	
  				tmp_l = origin_l;
  
  				cor_v = tmp_v + stim::vec3<T>(-coef_left * origin_l, 0, 0);			// get the original corner vertex
  				desire_l = desire_l - origin_l * (1.0f - ratio / 1.0f);
  				origin_l = (T)origin_l * ratio / 1.0f;
b710cea4   Jiaming Guo   fixed saving imag...
1904
  				n = find_number_square(origin_l, desire_l, times);
b61addd8   Jiaming Guo   add adjustment fe...
1905
1906
1907
1908
  				
  				width = (T)origin_l / (2 * n);										// updates
  				height = (desire_l - origin_l) / (2 * n);
  				
057ea93b   Jiaming Guo   add flow velocity...
1909
1910
  				// there are cases that the fragment can not satisfy the requirement for width
  				if (width < times * radius || n == 0) {								// check feasibility
4bb615eb   Jiaming Guo   fixed generating ...
1911
1912
  					ratio = 1.0f;													// load original lengths
  					desire_l = tmp_d;												// loading back-ups										
b61addd8   Jiaming Guo   add adjustment fe...
1913
1914
1915
  					origin_l = tmp_l;
  
  					std::cout << "Warning: current ratio is not feasible, use full original line." << std::endl;
b710cea4   Jiaming Guo   fixed saving imag...
1916
  					n = find_number_square(origin_l, desire_l, times);
b61addd8   Jiaming Guo   add adjustment fe...
1917
1918
1919
  
  					width = (T)origin_l / (2 * n);									// updates
  					height = (desire_l - origin_l) / (2 * n);
8334680a   Jiaming Guo   add square wave c...
1920
  				}
b61addd8   Jiaming Guo   add adjustment fe...
1921
1922
1923
1924
1925
  			}
  
  			// check whether it needs 3D square-wave-like connections
  			if (height > threshold) {					// enbale 3D connections
  				
6332d1e7   Jiaming Guo   fixed warnings
1926
  				height = (T)((desire_l - (1 + 2 * n) * origin_l) / std::pow(2 * n, 2));	// compute new height in 3D structure
b61addd8   Jiaming Guo   add adjustment fe...
1927
  				while (height > threshold) {			// increase order to decrease height
8334680a   Jiaming Guo   add square wave c...
1928
1929
  					n++;
  					width = (T)(origin_l) / (2 * n);
6332d1e7   Jiaming Guo   fixed warnings
1930
  					height = (T)((desire_l - (1 + 2 * n) * origin_l) / std::pow(2 * n, 2));
b61addd8   Jiaming Guo   add adjustment fe...
1931
  					// check whether it appears overlap, if it appears we choose last time height even if it is larger than threshold
8334680a   Jiaming Guo   add square wave c...
1932
1933
1934
  					if (width < times * radius) {
  						n--;
  						width = (T)(origin_l) / (2 * n);
6332d1e7   Jiaming Guo   fixed warnings
1935
  						height = (T)((desire_l - (1 + 2 * n) * origin_l) / std::pow(2 * n, 2));
8334680a   Jiaming Guo   add square wave c...
1936
1937
1938
  						break;
  					}
  				}
b61addd8   Jiaming Guo   add adjustment fe...
1939
1940
  
  				// check overlap in terms of height, has potential risk when both height and width are less than times * radius.
8334680a   Jiaming Guo   add square wave c...
1941
1942
1943
  				while (height < times * radius) {
  					n--;
  					width = (T)(origin_l) / (2 * n);
6332d1e7   Jiaming Guo   fixed warnings
1944
  					height = (T)((desire_l - (1 + 2 * n) * origin_l) / std::pow(2 * n, 2));
8334680a   Jiaming Guo   add square wave c...
1945
1946
  				}
  
4bb615eb   Jiaming Guo   fixed generating ...
1947
  				/// degenerated case, use original extending method
b61addd8   Jiaming Guo   add adjustment fe...
1948
  				if (n == 0) {
4bb615eb   Jiaming Guo   fixed generating ...
1949
1950
1951
1952
1953
1954
1955
1956
1957
1958
1959
1960
1961
  					height = (T)((desire_l - origin_l) / 2);
  					// "up"
  					tmp_v = tmp_v + stim::vec3<T>(0, coef_up * height, 0);
  					if (feeder == 1)
  						inlet[i].V.push_back(tmp_v);
  					else if (feeder == 0)
  						outlet[i].V.push_back(tmp_v);
  					// "left"
  					tmp_v = tmp_v + stim::vec3<T>(-coef_left * origin_l, 0, 0);
  					if (feeder == 1)
  						inlet[i].V.push_back(tmp_v);
  					else if (feeder == 0)
  						outlet[i].V.push_back(tmp_v);
b61addd8   Jiaming Guo   add adjustment fe...
1962
1963
  				}
  			
8334680a   Jiaming Guo   add square wave c...
1964
1965
1966
1967
1968
1969
1970
1971
1972
1973
1974
1975
1976
1977
1978
1979
1980
1981
1982
1983
1984
1985
  				// cube-like structure construction
  				for (int j = 0; j < n; j++) {
  					// "up"
  					for (int k = 0; k < n; k++) {
  						// in
  						tmp_v = tmp_v + stim::vec3<T>(0, 0, coef_z * height);
  						if (feeder == 1)
  							inlet[i].V.push_back(tmp_v);
  						else if (feeder == 0)
  							outlet[i].V.push_back(tmp_v);
  						// "up"
  						tmp_v = tmp_v + stim::vec3<T>(0, inverse * coef_up * width, 0);
  						if (feeder == 1)
  							inlet[i].V.push_back(tmp_v);
  						else if (feeder == 0)
  							outlet[i].V.push_back(tmp_v);
  						// out
  						tmp_v = tmp_v + stim::vec3<T>(0, 0, -coef_z * height);
  						if (feeder == 1)
  							inlet[i].V.push_back(tmp_v);
  						else if (feeder == 0)
  							outlet[i].V.push_back(tmp_v);
b61addd8   Jiaming Guo   add adjustment fe...
1986
  						// "up"
8334680a   Jiaming Guo   add square wave c...
1987
1988
1989
1990
1991
1992
1993
1994
1995
1996
1997
1998
1999
2000
2001
2002
2003
2004
2005
2006
2007
  						tmp_v = tmp_v + stim::vec3<T>(0, inverse * coef_up * width, 0);
  						if (feeder == 1)
  							inlet[i].V.push_back(tmp_v);
  						else if (feeder == 0)
  							outlet[i].V.push_back(tmp_v);
  					}
  
  					// "left"
  					tmp_v = tmp_v + stim::vec3<T>(-coef_left * width, 0, 0);
  					if (feeder == 1)
  						inlet[i].V.push_back(tmp_v);
  					else if (feeder == 0)
  						outlet[i].V.push_back(tmp_v);
  
  					if (inverse == 1)					// revert inverse
  						inverse = -1;
  					else
  						inverse = 1;
  
  					// "down"
  					for (int k = 0; k < n; k++) {
b61addd8   Jiaming Guo   add adjustment fe...
2008
  						// in
8334680a   Jiaming Guo   add square wave c...
2009
2010
2011
2012
2013
  						tmp_v = tmp_v + stim::vec3<T>(0, 0, coef_z * height);
  						if (feeder == 1)
  							inlet[i].V.push_back(tmp_v);
  						else if (feeder == 0)
  							outlet[i].V.push_back(tmp_v);
b61addd8   Jiaming Guo   add adjustment fe...
2014
  						// get the bounding box edge
f4105b89   Jiaming Guo   add new function:...
2015
2016
  						if (j == n - 1 && k == 0)		// first time go "in"
  							tmp_bb.second = tmp_v;
b61addd8   Jiaming Guo   add adjustment fe...
2017
  						// "down"
8334680a   Jiaming Guo   add square wave c...
2018
2019
2020
2021
2022
  						tmp_v = tmp_v + stim::vec3<T>(0, inverse * coef_up * width, 0);
  						if (feeder == 1)
  							inlet[i].V.push_back(tmp_v);
  						else if (feeder == 0)
  							outlet[i].V.push_back(tmp_v);
b61addd8   Jiaming Guo   add adjustment fe...
2023
  						// out
8334680a   Jiaming Guo   add square wave c...
2024
2025
2026
2027
2028
  						tmp_v = tmp_v + stim::vec3<T>(0, 0, -coef_z * height);
  						if (feeder == 1)
  							inlet[i].V.push_back(tmp_v);
  						else if (feeder == 0)
  							outlet[i].V.push_back(tmp_v);
b61addd8   Jiaming Guo   add adjustment fe...
2029
  						// "down"
8334680a   Jiaming Guo   add square wave c...
2030
2031
2032
2033
2034
2035
2036
2037
2038
2039
2040
2041
2042
2043
2044
2045
2046
2047
2048
2049
  						tmp_v = tmp_v + stim::vec3<T>(0, inverse * coef_up * width, 0);
  						if (feeder == 1)
  							inlet[i].V.push_back(tmp_v);
  						else if (feeder == 0)
  							outlet[i].V.push_back(tmp_v);
  					}
  
  					// "left"
  					tmp_v = tmp_v + stim::vec3<T>(-coef_left * width, 0, 0);
  					if (feeder == 1)
  						inlet[i].V.push_back(tmp_v);
  					else if (feeder == 0)
  						outlet[i].V.push_back(tmp_v);
  
  					if (inverse == 1)					// revert inverse
  						inverse = -1;
  					else
  						inverse = 1;
  				}
  				// if use fragment to do square wave connection, need to push_back the corner vertex
f4105b89   Jiaming Guo   add new function:...
2050
  				if (ratio > 0.0f && ratio < 1.0f) {
8334680a   Jiaming Guo   add square wave c...
2051
2052
2053
2054
2055
2056
  					if (feeder == 1)
  						inlet[i].V.push_back(cor_v);
  					else if (feeder == 0)
  						outlet[i].V.push_back(cor_v);
  				}
  			}
b61addd8   Jiaming Guo   add adjustment fe...
2057
  			// use 2D square-wave-like connections
8334680a   Jiaming Guo   add square wave c...
2058
  			else {
70c0b942   Jiaming Guo   fixed bug when he...
2059
2060
2061
  				if (height < times * radius) {			// if height is too small, decrease n and re-calculate height and width
  					height = times * radius;
  					T need_l = desire_l - origin_l;
6332d1e7   Jiaming Guo   fixed warnings
2062
  					n = (int)(need_l / (2 * height));
b61addd8   Jiaming Guo   add adjustment fe...
2063
  					if (n == 0)							// degenerated case
70c0b942   Jiaming Guo   fixed bug when he...
2064
2065
2066
2067
  						n = 1;
  					height = need_l / (2 * n);
  					width = origin_l / (2 * n);
  				}
8334680a   Jiaming Guo   add square wave c...
2068
2069
  				for (int j = 0; j < n; j++) {
  
f4105b89   Jiaming Guo   add new function:...
2070
  					// up
8334680a   Jiaming Guo   add square wave c...
2071
2072
2073
2074
2075
2076
  					tmp_v = tmp_v + stim::vec3<T>(0, coef_up * height, 0);
  					if (feeder == 1)
  						inlet[i].V.push_back(tmp_v);
  					else if (feeder == 0)
  						outlet[i].V.push_back(tmp_v);
  
f4105b89   Jiaming Guo   add new function:...
2077
  					// left
8334680a   Jiaming Guo   add square wave c...
2078
2079
2080
2081
2082
  					tmp_v = tmp_v + stim::vec3<T>(-coef_left * width, 0, 0);
  					if (feeder == 1)
  						inlet[i].V.push_back(tmp_v);
  					else if (feeder == 0)
  						outlet[i].V.push_back(tmp_v);
f4105b89   Jiaming Guo   add new function:...
2083
2084
  					if (j == n - 1)
  						tmp_bb.second = tmp_v;
8334680a   Jiaming Guo   add square wave c...
2085
  
f4105b89   Jiaming Guo   add new function:...
2086
  					// down
8334680a   Jiaming Guo   add square wave c...
2087
2088
2089
2090
2091
2092
  					tmp_v = tmp_v + stim::vec3<T>(0, -coef_up * height, 0);
  					if (feeder == 1)
  						inlet[i].V.push_back(tmp_v);
  					else if (feeder == 0)
  						outlet[i].V.push_back(tmp_v);
  
f4105b89   Jiaming Guo   add new function:...
2093
  					// left
8334680a   Jiaming Guo   add square wave c...
2094
2095
2096
2097
2098
2099
  					tmp_v = tmp_v + stim::vec3<T>(-coef_left * width, 0, 0);
  					if (feeder == 1)
  						inlet[i].V.push_back(tmp_v);
  					else if (feeder == 0)
  						outlet[i].V.push_back(tmp_v);
  				}
b61addd8   Jiaming Guo   add adjustment fe...
2100
2101
2102
2103
2104
2105
2106
  				// if use fragment to do square wave connection, need to push_back the corner vertex
  				if (ratio > 0.0f && ratio < 1.0f) {
  					if (feeder == 1)
  						inlet[i].V.push_back(cor_v);
  					else if (feeder == 0)
  						outlet[i].V.push_back(cor_v);
  				}
8334680a   Jiaming Guo   add square wave c...
2107
  			}
f4105b89   Jiaming Guo   add new function:...
2108
2109
2110
2111
  			if (feeder == 1)
  				inbb[i] = tmp_bb;
  			else if (feeder == 0)
  				outbb[i] = tmp_bb;
8334680a   Jiaming Guo   add square wave c...
2112
2113
  		}
  
9191c39e   Jiaming Guo   first version of ...
2114
  		// automatically modify bridge to make it feasible
6332d1e7   Jiaming Guo   fixed warnings
2115
  		void modify_synthetic_connection(T viscosity, T rou, bool H, T threshold, T &in, T &out, int times = 10, T ratio = 0.0f, T radius = 5.0f) {
9b69bb4c   Jiaming Guo   add diplay list t...
2116
2117
2118
  
  			glDeleteLists(dlist, 1);					// delete display list for modify
  			glDeleteLists(dlist + 1, 1);
f4105b89   Jiaming Guo   add new function:...
2119
  			
8334680a   Jiaming Guo   add square wave c...
2120
  			// because of radius change at the port vertex, there will be a pressure drop at that port
9191c39e   Jiaming Guo   first version of ...
2121
2122
2123
2124
2125
2126
2127
2128
  			// it follows the bernoulli equation
  			// p1 + 1/2*rou*v1^2 + rou*g*h1 = p2 + 1/2*rou*v2^2 + rou*g*h2
  			// Q1 = Q2 -> v1*r1^2 = v2*r2^2
  			std::vector<T> new_pressure = pressure;
  			unsigned idx;
  			for (unsigned i = 0; i < pendant_vertex.size(); i++) {
  				idx = pendant_vertex[i];
  				T tmp_v = get_velocity(idx);			// velocity at that pendant vertex
8334680a   Jiaming Guo   add square wave c...
2129
  				T ar = get_radius(idx) / radius;
9191c39e   Jiaming Guo   first version of ...
2130
2131
2132
2133
2134
2135
2136
2137
2138
  				new_pressure[idx] = pressure[idx] + 1.0f / 2.0f * rou * std::pow(tmp_v, 2) * (1.0f - std::pow(ar, 4));
  			}
  
  			// increase r -> increase Q -> decrease l
  			// find maximum pressure inlet port
  			T source_pressure = FLT_MIN;	// source pressure
  			unsigned inlet_index;
  			T tmp_p;
  			for (unsigned i = 0; i < inlet.size(); i++) {
8334680a   Jiaming Guo   add square wave c...
2139
  				tmp_p = new_pressure[inlet[i].v[0]] + ((8 * viscosity * inlet[i].l * inlet[i].Q) / ((float)stim::PI * std::pow(radius, 4)));
9191c39e   Jiaming Guo   first version of ...
2140
2141
2142
2143
2144
  				if (tmp_p > source_pressure) {
  					source_pressure = tmp_p;
  					inlet_index = i;
  				}
  			}
b61addd8   Jiaming Guo   add adjustment fe...
2145
  			Ps = source_pressure;
9191c39e   Jiaming Guo   first version of ...
2146
  
8334680a   Jiaming Guo   add square wave c...
2147
2148
2149
2150
2151
2152
2153
2154
2155
2156
  			// find minimum pressure outlet port
  			T end_pressure = FLT_MAX;
  			unsigned outlet_index;
  			for (unsigned i = 0; i < outlet.size(); i++) {
  				tmp_p = new_pressure[outlet[i].v[0]] - ((8 * viscosity * outlet[i].l * outlet[i].Q) / ((float)stim::PI * std::pow(radius, 4)));
  				if (tmp_p < end_pressure) {
  					end_pressure = tmp_p;
  					outlet_index = i;
  				}
  			}
b61addd8   Jiaming Guo   add adjustment fe...
2157
  			Pe = end_pressure;
9191c39e   Jiaming Guo   first version of ...
2158
  
8334680a   Jiaming Guo   add square wave c...
2159
2160
2161
2162
2163
2164
2165
2166
2167
2168
2169
2170
  			// automatically modify inlet bridge using Hilbert curves
  			if (H) {
  				bool upper = false;						// flag indicates the whether the port is upper than main feeder
  				bool invert = false;					// there are two version of hilbert curve depends on starting position with respect to the cup
  				T new_l;
  				stim::vec3<T> bus_v;					// the port point on the bus
  				stim::vec3<T> mid_v;					// the original corner point
  				stim::vec3<T> tmp_v;					// the pendant point
  				int order = 0;							// order of hilbert curve (iteration)
  				for (unsigned i = 0; i < inlet.size(); i++) {
  					if (i != inlet_index) {
  						new_l = (source_pressure - new_pressure[inlet[i].v[0]]) * ((float)stim::PI * std::pow(radius, 4)) / (8 * viscosity * inlet[i].Q);
9191c39e   Jiaming Guo   first version of ...
2171
  
8334680a   Jiaming Guo   add square wave c...
2172
2173
2174
2175
2176
2177
2178
2179
2180
2181
2182
2183
2184
2185
2186
2187
2188
2189
2190
2191
2192
2193
2194
2195
2196
2197
  						if (inlet[i].V[2][1] > main_feeder[0][1]) {		// check out upper side of lower side
  							upper = true;
  							invert = false;
  						}
  						else {
  							upper = false;
  							invert = true;
  						}
  
  						T origin_l = (inlet[i].V[1] - inlet[i].V[2]).len();
  						T desire_l = new_l - (inlet[i].V[0] - inlet[i].V[1]).len();
  						find_hilbert_order(origin_l, desire_l, order);
  
  						bus_v = inlet[i].V[0];
  						mid_v = inlet[i].V[1];
  						tmp_v = inlet[i].V[2];
  						inlet[i].V.clear();
  						inlet[i].V.push_back(tmp_v);
  						inlet[i].l = new_l;
  
  						if (desire_l - origin_l < 2 * radius) {	// do not need to use hilbert curve, just increase the length by draging out
  							T d = new_l - inlet[i].l;
  							stim::vec3<T> corner = stim::vec3<T>(tmp_v[0], tmp_v[1] + d / 2.0f * (tmp_v[1] > main_feeder[0][1] ? 1 : -1), tmp_v[2]);
  							inlet[i].V.push_back(corner);
  							corner = stim::vec3<T>(mid_v[0], mid_v[1] + d / 2.0f * (tmp_v[1] > main_feeder[0][1] ? 1 : -1), mid_v[2]);
  							inlet[i].V.push_back(corner);
6765b32b   Jiaming Guo   add hilbert curve
2198
2199
  							inlet[i].V.push_back(bus_v);
  						}
8334680a   Jiaming Guo   add square wave c...
2200
  						else {
6332d1e7   Jiaming Guo   fixed warnings
2201
2202
  							T fragment = (T)((desire_l - origin_l) / ((std::pow(4, order) - 1) / (std::pow(2, order) - 1) - 1));	// the length of the opening of cup 		
  							T dl = (T)(fragment / (std::pow(2, order) - 1));											// unit cup length
6765b32b   Jiaming Guo   add hilbert curve
2203
  
8334680a   Jiaming Guo   add square wave c...
2204
2205
2206
2207
2208
  							if (dl > 2 * radius) {				// if the radius is feasible
  								if (upper)
  									hilbert_curve(i, &tmp_v[0], order, dl, 1, invert, DOWN);
  								else
  									hilbert_curve(i, &tmp_v[0], order, dl, 1, invert, UP);
6765b32b   Jiaming Guo   add hilbert curve
2209
  
8334680a   Jiaming Guo   add square wave c...
2210
2211
2212
2213
2214
2215
2216
  								if (tmp_v[0] != mid_v[0])
  									inlet[i].V.push_back(mid_v);
  								inlet[i].V.push_back(bus_v);
  							}
  							else {								// if the radius is not feasible
  								int count = 1;
  								while (dl <= 2 * radius) {
6332d1e7   Jiaming Guo   fixed warnings
2217
  									dl = (T)(origin_l / (std::pow(2, order - count) - 1));
8334680a   Jiaming Guo   add square wave c...
2218
2219
2220
2221
2222
2223
2224
2225
2226
  									count++;
  								}
  								count--;
  
  								if (upper)
  									hilbert_curve(i, &tmp_v[0], order - count, dl, 1, invert, DOWN);
  								else
  									hilbert_curve(i, &tmp_v[0], order - count, dl, 1, invert, UP);
  
6332d1e7   Jiaming Guo   fixed warnings
2227
  								desire_l -= (T)(origin_l * ((std::pow(4, order - count) - 1) / (std::pow(2, order - count) - 1)));
8334680a   Jiaming Guo   add square wave c...
2228
2229
2230
2231
2232
  								origin_l = (bus_v - mid_v).len();
  								desire_l += origin_l;
  
  								find_hilbert_order(origin_l, desire_l, order);
  
6332d1e7   Jiaming Guo   fixed warnings
2233
2234
  								fragment = (T)((desire_l - origin_l) / ((std::pow(4, order) - 1) / (std::pow(2, order) - 1) - 1));
  								dl = (T)(fragment / (std::pow(2, order) - 1));
8334680a   Jiaming Guo   add square wave c...
2235
2236
2237
2238
2239
2240
2241
2242
2243
2244
2245
2246
2247
2248
2249
  								if (dl < 2 * radius)
  									std::cout << "infeasible connection between inlets!" << std::endl;
  
  								if (upper)
  									hilbert_curve(i, &tmp_v[0], order, dl, 1, !invert, LEFT);
  								else
  									hilbert_curve(i, &tmp_v[0], order, dl, 1, !invert, RIGHT);
  
  								if (tmp_v[1] != bus_v[1])
  									inlet[i].V.push_back(bus_v);
  							}
  						}
  						std::reverse(inlet[i].V.begin(), inlet[i].V.end());			// from bus to pendant vertex
  					}
  				}
6765b32b   Jiaming Guo   add hilbert curve
2250
  
8334680a   Jiaming Guo   add square wave c...
2251
2252
2253
2254
  				// automatically modify outlet bridge to make it feasible
  				for (unsigned i = 0; i < outlet.size(); i++) {
  					if (i != outlet_index) {
  						new_l = (new_pressure[outlet[i].v[0]] - end_pressure) * ((float)stim::PI * std::pow(radius, 4)) / (8 * viscosity * outlet[i].Q);
6765b32b   Jiaming Guo   add hilbert curve
2255
  
8334680a   Jiaming Guo   add square wave c...
2256
2257
2258
2259
2260
2261
2262
2263
  						if (outlet[i].V[2][1] > main_feeder[1][1]) {
  							upper = true;
  							invert = true;
  						}
  						else {
  							upper = false;
  							invert = false;
  						}
6765b32b   Jiaming Guo   add hilbert curve
2264
  
8334680a   Jiaming Guo   add square wave c...
2265
2266
2267
2268
2269
2270
2271
2272
2273
2274
2275
2276
2277
2278
2279
2280
2281
2282
2283
2284
  						T origin_l = (outlet[i].V[1] - outlet[i].V[2]).len();
  						T desire_l = new_l - (outlet[i].V[0] - outlet[i].V[1]).len();
  						find_hilbert_order(origin_l, desire_l, order);
  
  						bus_v = outlet[i].V[0];
  						mid_v = outlet[i].V[1];
  						tmp_v = outlet[i].V[2];
  						outlet[i].V.clear();
  						outlet[i].V.push_back(tmp_v);
  						outlet[i].l = new_l;
  
  						if (desire_l - origin_l < 2 * radius) {	// do not need to use hilbert curve, just increase the length by draging out
  							T d = new_l - outlet[i].l;
  							stim::vec3<T> corner = stim::vec3<T>(tmp_v[0], tmp_v[1] + d / 2.0f * (tmp_v[1] > main_feeder[0][1] ? 1 : -1), tmp_v[2]);
  							outlet[i].V.push_back(corner);
  							corner = stim::vec3<T>(mid_v[0], mid_v[1] + d / 2.0f * (tmp_v[1] > main_feeder[0][1] ? 1 : -1), mid_v[2]);
  							outlet[i].V.push_back(corner);
  							outlet[i].V.push_back(bus_v);
  						}
  						else {
6332d1e7   Jiaming Guo   fixed warnings
2285
2286
  							T fragment = (T)((desire_l - origin_l) / ((std::pow(4, order) - 1) / (std::pow(2, order) - 1) - 1));	// the length of the opening of cup 		
  							T dl = (T)(fragment / (std::pow(2, order) - 1));											// unit cup length
6765b32b   Jiaming Guo   add hilbert curve
2287
  
8334680a   Jiaming Guo   add square wave c...
2288
2289
2290
2291
2292
2293
2294
2295
2296
2297
2298
2299
2300
  							if (dl > 2 * radius) {				// if the radius is feasible
  								if (upper)
  									hilbert_curve(i, &tmp_v[0], order, dl, 0, invert, DOWN);
  								else
  									hilbert_curve(i, &tmp_v[0], order, dl, 0, invert, UP);
  
  								if (tmp_v[0] != mid_v[0])
  									outlet[i].V.push_back(mid_v);
  								outlet[i].V.push_back(bus_v);
  							}
  							else {								// if the radius is not feasible
  								int count = 1;
  								while (dl <= 2 * radius) {
6332d1e7   Jiaming Guo   fixed warnings
2301
  									dl = (T)(origin_l / (std::pow(2, order - count) - 1));
8334680a   Jiaming Guo   add square wave c...
2302
2303
2304
2305
2306
2307
2308
2309
2310
  									count++;
  								}
  								count--;
  
  								if (upper)
  									hilbert_curve(i, &tmp_v[0], order - count, dl, 0, invert, DOWN);
  								else
  									hilbert_curve(i, &tmp_v[0], order - count, dl, 0, invert, UP);
  
6332d1e7   Jiaming Guo   fixed warnings
2311
  								desire_l -= (T)(origin_l * ((std::pow(4, order - count) - 1) / (std::pow(2, order - count) - 1)));
8334680a   Jiaming Guo   add square wave c...
2312
2313
2314
2315
2316
  								origin_l = (bus_v - mid_v).len();
  								desire_l += origin_l;
  
  								find_hilbert_order(origin_l, desire_l, order);
  
6332d1e7   Jiaming Guo   fixed warnings
2317
2318
  								fragment = (T)((desire_l - origin_l) / ((std::pow(4, order) - 1) / (std::pow(2, order) - 1) - 1));
  								dl = (T)(fragment / (std::pow(2, order) - 1));
8334680a   Jiaming Guo   add square wave c...
2319
2320
2321
2322
2323
2324
2325
2326
2327
2328
2329
  								if (dl < 2 * radius)
  									std::cout << "infeasible connection between outlets!" << std::endl;
  
  								if (upper)
  									hilbert_curve(i, &tmp_v[0], order, dl, 0, !invert, LEFT);
  								else
  									hilbert_curve(i, &tmp_v[0], order, dl, 0, !invert, RIGHT);
  
  								if (tmp_v[1] != bus_v[1])
  									outlet[i].V.push_back(bus_v);
  							}
6765b32b   Jiaming Guo   add hilbert curve
2330
  						}
8334680a   Jiaming Guo   add square wave c...
2331
  						std::reverse(outlet[i].V.begin(), outlet[i].V.end());
6765b32b   Jiaming Guo   add hilbert curve
2332
  					}
9191c39e   Jiaming Guo   first version of ...
2333
2334
  				}
  			}
8334680a   Jiaming Guo   add square wave c...
2335
2336
2337
  			// automatically modify inlet bridge using square shape constructions
  			else {
  				bool upper;								// flag indicates the connection is upper than the bus
f4105b89   Jiaming Guo   add new function:...
2338
  				bool z;									// flag indicates the connection direction along z-axis
8334680a   Jiaming Guo   add square wave c...
2339
2340
2341
2342
2343
2344
  				T new_l;								// new length
  				stim::vec3<T> bus_v;					// the port point on the bus
  				stim::vec3<T> mid_v;					// the original corner point
  				stim::vec3<T> tmp_v;					// the pendant point
  				int n;
  				T width, height;						// width and height of the square
f4105b89   Jiaming Guo   add new function:...
2345
2346
  				inbb.resize(inlet.size());				// resize bounding box of inlets/outlets connections
  				outbb.resize(outlet.size());
9191c39e   Jiaming Guo   first version of ...
2347
  
8334680a   Jiaming Guo   add square wave c...
2348
2349
2350
  				for (unsigned i = 0; i < inlet.size(); i++) {
  					if (i != inlet_index) {
  						new_l = (source_pressure - new_pressure[inlet[i].v[0]]) * ((float)stim::PI * std::pow(radius, 4)) / (8 * viscosity * inlet[i].Q);	// calculate the new length of the connection 
9191c39e   Jiaming Guo   first version of ...
2351
  
8334680a   Jiaming Guo   add square wave c...
2352
2353
  						bus_v = inlet[i].V[0];
  						mid_v = inlet[i].V[1];
70c0b942   Jiaming Guo   fixed bug when he...
2354
  						tmp_v = inlet[i].V[2];										// not always pendant vertex
6765b32b   Jiaming Guo   add hilbert curve
2355
  
8334680a   Jiaming Guo   add square wave c...
2356
2357
2358
2359
  						if (inlet[i].V[2][1] > main_feeder[0][1]) 					// check out upper side of lower side
  							upper = true;
  						else
  							upper = false;
9191c39e   Jiaming Guo   first version of ...
2360
  
8334680a   Jiaming Guo   add square wave c...
2361
2362
2363
2364
2365
2366
2367
  						if (inlet[i].V[2][2] > main_feeder[0][2])
  							z = true;
  						else
  							z = false;
  
  						T origin_l = (inlet[i].V[1] - inlet[i].V[2]).len();
  						T desire_l = new_l - (inlet[i].V[0] - inlet[i].V[1]).len();
70c0b942   Jiaming Guo   fixed bug when he...
2368
2369
2370
2371
2372
2373
2374
2375
2376
2377
2378
  						if (inlet[i].V.size() != 3) {
  							desire_l = new_l - (inlet[i].V[0] - inlet[i].V[1]).len() - (inlet[i].V[2] - inlet[i].V[3]).len();
  							stim::vec3<T> tmp = inlet[i].V[3];
  							inlet[i].V.clear();
  							inlet[i].V.push_back(tmp);
  							inlet[i].V.push_back(tmp_v);
  						}
  						else {
  							inlet[i].V.clear();
  							inlet[i].V.push_back(tmp_v);
  						}
8334680a   Jiaming Guo   add square wave c...
2379
2380
  						inlet[i].l = new_l;
  
6332d1e7   Jiaming Guo   fixed warnings
2381
  						n = find_number_square(origin_l, desire_l, times);
8334680a   Jiaming Guo   add square wave c...
2382
2383
2384
2385
  
  						width = (T)origin_l / (2 * n);
  						height = (desire_l - origin_l) / (2 * n);
  
0224d2ef   Jiaming Guo   fixed square wave...
2386
  						build_square_connection(i, width, height, origin_l, desire_l, n, 1, threshold, z, true, upper, 2, ratio);
8334680a   Jiaming Guo   add square wave c...
2387
2388
2389
  						inlet[i].V.push_back(bus_v);
  
  						std::reverse(inlet[i].V.begin(), inlet[i].V.end());			// from bus to pendant vertex
6765b32b   Jiaming Guo   add hilbert curve
2390
  					}
f4105b89   Jiaming Guo   add new function:...
2391
2392
2393
2394
  					else {
  						inbb[i].first = inlet[i].V[2];
  						inbb[i].second = inlet[i].V[1];
  					}
8334680a   Jiaming Guo   add square wave c...
2395
2396
2397
2398
2399
  				}
  
  				for (unsigned i = 0; i < outlet.size(); i++) {
  					if (i != outlet_index) {
  						new_l = (new_pressure[outlet[i].v[0]] - end_pressure) * ((float)stim::PI * std::pow(radius, 4)) / (8 * viscosity * outlet[i].Q);	// calculate the new length of the connection 
6765b32b   Jiaming Guo   add hilbert curve
2400
  
8334680a   Jiaming Guo   add square wave c...
2401
2402
2403
  						bus_v = outlet[i].V[0];
  						mid_v = outlet[i].V[1];
  						tmp_v = outlet[i].V[2];
6765b32b   Jiaming Guo   add hilbert curve
2404
  
8334680a