17 #include <GL/freeglut.h>
19 #include "stream_view.h"
27 StreamView::StreamView(
const char* title, WinGeom* wg)
28 : View(title, wg), vec(NULL)
40 StreamView::StreamView(
char* title, WinGeom* wg)
41 : View(title, wg), vec(NULL)
58 throw Hermes::Exceptions::Exception(
"Identical solutions passed to the two-argument version of show(). This is most likely a mistake.");
59 show(xsln, ysln, marker, step, eps, H2D_FN_VAL_0, H2D_FN_VAL_0);
62 bool StreamView::is_in_triangle(
int idx,
double x,
double y, double3& bar)
64 double4* vert = vec->get_vertices();
65 int3* xtris = vec->get_triangles();
66 int3& tri = xtris[idx];
67 double x1 = vert[tri[0]][0], x2 = vert[tri[1]][0], x3 = vert[tri[2]][0];
68 double y1 = vert[tri[0]][1], y2 = vert[tri[1]][1], y3 = vert[tri[2]][1];
69 double jac = ((x1 - x3)*(y2 - y3) - (x2 - x3)*(y1 - y3));
70 double eps = jac * 1e-8;
71 double a = ((y2 - y3) * (x - x3) - (x2 - x3) * (y - y3));
72 double b = ((x1 - x3) * (y - y3) - (y1 - y3) * (x - x3));
73 double c = jac - a - b;
74 bar[0] = a / jac; bar[1] = b / jac; bar[2] = c / jac;
75 if((a >= - eps && a <= jac + eps) && (b >= - eps && b <= jac + eps) && (c >= - eps && c <= jac + eps))
81 void StreamView::add_element_to_tree(
Node* father,
int e_idx,
double x_min,
double x_max,
double y_min,
double y_max)
83 double4* vert = vec->get_vertices();
84 int3* xtris = vec->get_triangles();
85 if(father->leaf ==
true)
87 father->elements[father->num_elem++] = e_idx;
88 if(father->num_elem >= 100)
91 for (
int k = 0; k < 2; k++)
93 father->sons[k] =
new Node;
94 father->sons[k]->level = father->level + 1;
95 father->sons[k]->leaf =
true;
96 father->sons[k]->num_elem = 0;
98 for (
int i = 0; i < father->num_elem; i++)
99 add_element_to_tree(father, father->elements[i], x_min, x_max, y_min, y_max);
104 double x_mid = (x_min + x_max)/2;
105 double y_mid = (y_min + y_max)/2;
106 int3& tri = xtris[e_idx];
107 if(father->level % 2)
109 if(vert[tri[0]][1] <= y_mid || vert[tri[1]][1] <= y_mid || vert[tri[2]][1] <= y_mid)
110 add_element_to_tree(father->sons[0], e_idx, x_min, x_max, y_min, y_mid);
111 if(vert[tri[0]][1] >= y_mid || vert[tri[1]][1] >= y_mid || vert[tri[2]][1] >= y_mid)
112 add_element_to_tree(father->sons[1], e_idx, x_min, x_max, y_mid, y_max);
116 if(vert[tri[0]][0] <= x_mid || vert[tri[1]][0] <= x_mid || vert[tri[2]][0] <= x_mid)
117 add_element_to_tree(father->sons[0], e_idx, x_min, x_mid, y_min, y_max);
118 if(vert[tri[0]][0] >= x_mid || vert[tri[1]][0] >= x_mid || vert[tri[2]][0] >= x_mid)
119 add_element_to_tree(father->sons[1], e_idx, x_mid, x_max, y_min, y_max);
124 void StreamView::build_tree()
129 for (
int i = 0; i < vec->get_num_triangles(); i++)
131 add_element_to_tree(root, i, root_x_min, root_x_max, root_y_min, root_y_max);
135 int StreamView::find_triangle_in_tree(
double x,
double y,
Node* father,
double x_min,
double x_max,
double y_min,
double y_max, double3& bar)
137 if(father->leaf ==
true)
139 double4* vert = vec->get_vertices();
140 int3* xtris = vec->get_triangles();
141 for (
int idx = 0; idx < father->num_elem; idx++)
143 int i = father->elements[idx];
144 if(is_in_triangle(i, x, y, bar))
return i;
150 double x_mid = (x_max + x_min)/2;
151 double y_mid = (y_max + y_min)/2;
152 if(father->level % 2)
154 return find_triangle_in_tree(x, y, father->sons[0], x_min, x_max, y_min, y_mid, bar);
156 return find_triangle_in_tree(x, y, father->sons[1], x_min, x_max, y_mid, y_max, bar);
159 return find_triangle_in_tree(x, y, father->sons[0], x_min, x_mid, y_min, y_max, bar);
161 return find_triangle_in_tree(x, y, father->sons[1], x_mid, x_max, y_min, y_max, bar);
165 bool StreamView::get_solution_values(
double x,
double y,
double& xval,
double& yval)
167 double4* vert = vec->get_vertices();
168 int3* xtris = vec->get_triangles();
171 if((e_idx = find_triangle_in_tree(x, y, root, root_x_min, root_x_max, root_y_min, root_y_max, bar)) == -1)
return false;
172 int3& tri = xtris[e_idx];
173 xval = bar[0] * vert[tri[0]][2] + bar[1] * vert[tri[1]][2] + bar[2] * vert[tri[2]][2];
174 yval = bar[0] * vert[tri[0]][3] + bar[1] * vert[tri[1]][3] + bar[2] * vert[tri[2]][3];
178 void StreamView::delete_tree(
Node* father)
180 if(father->leaf ==
false)
182 delete_tree(father->sons[0]);
183 delete_tree(father->sons[1]);
188 int StreamView::create_streamline(
double x_start,
double y_start,
int idx)
190 double ODE_EPS = 1e-5;
191 double tau = initial_tau;
195 const int buffer_length = 5000;
196 double2*
buffer =
new double2[buffer_length];
199 bool almost_end_of_domain =
false;
201 double k1, k2, k3, k4, k5;
202 double l1, l2, l3, l4, l5;
203 double x1, x2, x3, x4, x5;
204 double y1, y2, y3, y4, y5;
208 if(get_solution_values(x, y, k1, l1) ==
false)
211 if(tau < min_tau)
break;
214 if(fabs(k1) / max_mag < 1e-5 && fabs(l1) / max_mag < 1e-5)
break;
220 if(k >= buffer_length)
break;
224 if(almost_end_of_domain)
226 almost_end_of_domain =
false;
228 if(tau < min_tau) { end =
true;
break; }
232 x1 = x + 1.0/3.0 * tau * k1; y1 = y + 1.0/3.0 * tau * l1;
233 if(get_solution_values(x1, y1, k2, l2) ==
false) { almost_end_of_domain =
true;
continue; }
234 x2 = x + 1.0/6.0 * tau * (k1 + k2); y2 = y + 1.0/6.0 * tau * (l1 + l2);
235 if(get_solution_values(x2, y2, k3, l3) ==
false) { almost_end_of_domain =
true;
continue; }
236 x3 = x + tau * (0.125 * k1 + 0.375 * k3); y3 = y + tau * (0.125 * l1 + 0.375 * l3);
237 if(get_solution_values(x3, y3, k4, l4) ==
false) { almost_end_of_domain =
true;
continue; }
238 x4 = x + tau * (0.5 * k1 - 1.5 * k3 + 2.0 * k4); y4 = y + tau * (0.5 * l1 - 1.5 * l3 + 2.0 * l4);
239 if(get_solution_values(x4, y4, k5, l5) ==
false) { almost_end_of_domain =
true;
continue; }
240 x5 = x + tau * 1.0/6.0 * (k1 + 4.0 * k4 + k5); y5 = y + tau * 1.0/6.0 * (l1 + 4.0 * l4 + l5);
243 double x_err = 1.0/5.0 * (x4 - x5) / (root_x_max - root_x_min);
244 double y_err = 1.0/5.0 * (y4 - y5) / (root_y_max - root_y_min);
245 double err = std::max(fabs(x_err), fabs(y_err));
248 tau_ok =
true; x = x5; y = y5;
252 tau_ok =
false; tau = tau/2;
253 if(tau < min_tau) tau = min_tau;
259 tau = 0.8 * tau * pow(ODE_EPS/err, 0.2);
260 if(tau > max_tau) tau = max_tau;
267 streamlines[idx] =
new double2[k];
268 memcpy(streamlines[idx], buffer, k*
sizeof(double2));
274 static double4* comp_vert;
275 static int compare(
const void* p1,
const void* p2)
277 const int3* e1 = ((
const int3*) p1);
278 const int3* e2 = ((
const int3*) p2);
279 double x1 = comp_vert[(*e1)[0]][0];
280 double y1 = comp_vert[(*e1)[0]][1];
281 double x2 = comp_vert[(*e2)[0]][0];
282 double y2 = comp_vert[(*e2)[0]][1];
283 if(x1 < x2)
return -1;
284 if(x1 == x2 && y1 < y2)
return -1;
285 if(x1 > x2)
return 1;
286 if(x1 == x2 && y1 > y2)
return 1;
287 if(x1 == x2 && y1 == y2)
return 0;
288 throw Hermes::Exceptions::Exception(
"internal error: reached end of non-void function");
292 int StreamView::find_initial_edge(
int num_edges, int3* edges)
295 double4* vert = vec->get_vertices();
296 for (i = 0; i < num_edges; i++)
300 for (j = 0; j < num_edges; j++)
301 if(vert[edges[j][1]][0] == vert[edges[i][0]][0] && vert[edges[j][1]][1] == vert[edges[i][0]][1])
303 if(j == num_edges)
return i;
309 static int find_next_edge(
int num_edges, int3* edges,
int b_idx)
313 int3* edge = (int3*) bsearch(&key, edges, num_edges,
sizeof(int3), compare);
320 void StreamView::find_initial_points(
int marker,
double step, double2*& initial_points)
323 int ne = vec->get_num_edges();
324 int2* edges = vec->get_edges();
325 int* edge_markers = vec->get_edge_markers();
326 double4* vert = vec->get_vertices();
327 int3* bnd_edges =
new int3[ne];
328 for (
int i = 0; i < ne; i++)
330 if(edge_markers[i] == marker)
332 bnd_edges[k][0] = edges[i][0];
333 bnd_edges[k][1] = edges[i][1];
334 bnd_edges[k++][2] = 0;
342 qsort(bnd_edges, num_edges,
sizeof(int3), compare);
345 int buffer_length = 1000;
346 initial_points =
new double2[buffer_length];
349 while ((idx = find_initial_edge(num_edges, bnd_edges)) != -1)
351 double tmp_step = step;
354 bnd_edges[idx][2] = 1;
355 double ax = vert[bnd_edges[idx][0]][0];
double bx = vert[bnd_edges[idx][1]][0];
356 double ay = vert[bnd_edges[idx][0]][1];
double by = vert[bnd_edges[idx][1]][1];
357 double len = sqrt(sqr(bx - ax) + sqr(by - ay));
358 double remaining_len = len;
double init_x = ax;
double init_y = ay;
359 while (tmp_step < remaining_len)
361 initial_points[k][0] = init_x + tmp_step * ((bx - ax) / len);
362 initial_points[k][1] = init_y + tmp_step * ((by - ay) / len);
363 remaining_len = remaining_len - tmp_step;
365 init_x = initial_points[k][0];
366 init_y = initial_points[k][1];
369 tmp_step = tmp_step - remaining_len;
371 while ((idx = find_next_edge(num_edges, bnd_edges, bnd_edges[idx][1])) != -1);
382 vec->process_solution(xsln, ysln, xitem, yitem, eps);
387 range_min = vec->get_min_value();
388 range_max = vec->get_max_value();
390 vec->calc_vertices_aabb(&vertices_min_x, &vertices_max_x, &vertices_min_y, &vertices_max_y);
393 double4* vert = vec->get_vertices();
394 for (
int i = 0; i < vec->get_num_vertices(); i++)
396 if(vert[i][0] < root_x_min) root_x_min = vert[i][0];
397 if(vert[i][0] > root_x_max) root_x_max = vert[i][0];
398 if(vert[i][1] < root_y_min) root_y_min = vert[i][1];
399 if(vert[i][1] > root_y_max) root_y_max = vert[i][1];
402 initial_tau = std::max(root_x_max - root_x_min, root_y_max - root_y_min) / 100;
403 max_tau = initial_tau * 10;
404 min_tau = initial_tau / 50;
405 max_mag = vec->get_max_value();
411 double2* initial_points;
412 find_initial_points(marker, step, initial_points);
414 streamlines = (double2**) malloc(
sizeof(double2*) * (num_stream));
415 streamlength = (
int*) malloc(
sizeof(
int) * (num_stream));
416 for (
int i = 0; i < num_stream; i++)
417 streamlength[i] = create_streamline(initial_points[i][0], initial_points[i][1], i);
419 delete [] initial_points;
430 void StreamView::add_streamline(
double x,
double y)
433 throw Hermes::Exceptions::Exception(
"Function add_streamline must be called after StreamView::show().");
435 streamlines = (double2**) realloc(streamlines,
sizeof(double2*) * (num_stream + 1));
436 streamlength = (
int*) realloc(streamlength,
sizeof(
int) * (num_stream + 1));
437 streamlength[num_stream] = create_streamline(x, y, num_stream);
442 static int n_vert(
int i) {
return (i + 1) % 3; }
443 static int p_vert(
int i) {
return (i + 2) % 3; }
445 void StreamView::on_display()
447 set_ortho_projection();
448 glDisable(GL_LIGHTING);
449 glDisable(GL_DEPTH_TEST);
450 glDisable(GL_TEXTURE_1D);
451 glPolygonMode(GL_FRONT_AND_BACK, pmode ? GL_LINE : GL_FILL);
456 int nv = vec->get_num_vertices();
457 double4* vert = vec->get_vertices();
458 double2* tvert =
new double2[nv];
460 for (i = 0; i < nv; i++)
462 tvert[i][0] = transform_x(vert[i][0]);
463 tvert[i][1] = transform_y(vert[i][1]);
467 double min = range_min, max = range_max;
468 if(range_auto) { min = vec->get_min_value(); max = vec->get_max_value(); }
469 double irange = 1.0 / (max - min);
471 if(fabs(min - max) < 1e-8) { irange = 1.0; min -= 0.5; }
474 int3* xtris = vec->get_triangles();
476 glEnable(GL_TEXTURE_1D);
477 glBindTexture(GL_TEXTURE_1D, gl_pallete_tex_id);
478 glBegin(GL_TRIANGLES);
479 glColor3f(0.95f, 0.95f, 0.95f);
480 for (i = 0; i < vec->get_num_triangles(); i++)
482 double mag = sqrt(sqr(vert[xtris[i][0]][2]) + sqr(vert[xtris[i][0]][3]));
483 glTexCoord2d((mag -min) * irange * tex_scale + tex_shift, 0.0);
484 glVertex2d(tvert[xtris[i][0]][0], tvert[xtris[i][0]][1]);
486 mag = sqrt(sqr(vert[xtris[i][1]][2]) + sqr(vert[xtris[i][1]][3]));
487 glTexCoord2d((mag -min) * irange * tex_scale + tex_shift, 0.0);
488 glVertex2d(tvert[xtris[i][1]][0], tvert[xtris[i][1]][1]);
490 mag = sqrt(sqr(vert[xtris[i][2]][2]) + sqr(vert[xtris[i][2]][3]));
491 glTexCoord2d((mag -min) * irange * tex_scale + tex_shift, 0.0);
492 glVertex2d(tvert[xtris[i][2]][0], tvert[xtris[i][2]][1]);
495 glDisable(GL_TEXTURE_1D);
498 glColor3f(0.5, 0.5, 0.5);
500 int2* edges = vec->get_edges();
501 for (i = 0; i < vec->get_num_edges(); i++)
503 if(lines || edges[i][2] != 0)
505 glVertex2d(tvert[edges[i][0]][0], tvert[edges[i][0]][1]);
506 glVertex2d(tvert[edges[i][1]][0], tvert[edges[i][1]][1]);
512 glColor3f(0.0, 0.0, 0.0);
513 for (i = 0; i < num_stream; i++)
515 glBegin(GL_LINE_STRIP);
517 for (
int j = 0; j < streamlength[i] - 1; j++)
519 glVertex2d(transform_x(streamlines[i][j][0]), transform_y(streamlines[i][j][1]));
520 glVertex2d(transform_x(streamlines[i][j + 1][0]), transform_y(streamlines[i][j + 1][1]));
529 void StreamView::on_mouse_move(
int x,
int y)
531 View::on_mouse_move(x, y);
534 void StreamView::on_key_down(
unsigned char key,
int x,
int y)
554 set_palette_filter(pal_filter != GL_LINEAR);
562 delete [] streamlines[num_stream];
568 View::on_key_down(key, x, y);
573 void StreamView::on_left_mouse_down(
int x,
int y)
575 View::on_left_mouse_down(x, y);
578 if(!scale_focused && glutGetModifiers() == GLUT_ACTIVE_CTRL)
581 streamlines = (double2**) realloc(streamlines,
sizeof(double2*) * (num_stream + 1));
582 streamlength = (
int*) realloc(streamlength,
sizeof(
int) * (num_stream + 1));
583 streamlength[num_stream] = create_streamline(untransform_x(x), untransform_y(y), num_stream);
589 const char* StreamView::get_help_text()
const
594 " Left mouse - pan\n"
595 " Right mouse - zoom\n"
596 " CTRL + Left mouse click - add steamline\n"
597 " CTRL z - delete last streamline\n"
598 " C - center image\n"
599 " F - toggle smooth palette\n"
600 " H - render high-quality frame\n"
602 " P - cycle palettes\n"
603 " S - save screenshot\n"
608 StreamView::~StreamView()
611 for (
int i = 0; i < num_stream; i++)
612 delete [] streamlines[i];
613 delete [] streamlines;
614 delete [] streamlength;