Hermes2D  2.0
stream_view.cpp
1 // This file is part of Hermes2D.
2 //
3 // Hermes2D is free software: you can redistribute it and/or modify
4 // it under the terms of the GNU General Public License as published by
5 // the Free Software Foundation, either version 2 of the License, or
6 // (at your option) any later version.
7 //
8 // Hermes2D is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 // GNU General Public License for more details.
12 //
13 // You should have received a copy of the GNU General Public License
14 // along with Hermes2D. If not, see <http://www.gnu.org/licenses/>.
15 
16 #ifndef NOGLUT
17 #include <GL/freeglut.h>
18 #include "global.h"
19 #include "stream_view.h"
20 
21 namespace Hermes
22 {
23  namespace Hermes2D
24  {
25  namespace Views
26  {
27  StreamView::StreamView(const char* title, WinGeom* wg)
28  : View(title, wg), vec(NULL)
29  {
30  lines = false;
31  pmode = false;
32  num_stream = 0;
33  root_x_min = 1e100;
34  root_y_min = 1e100;
35  root_x_max = -1e100;
36  root_y_max = -1e100;
37  root = NULL;
38  }
39 
40  StreamView::StreamView(char* title, WinGeom* wg)
41  : View(title, wg), vec(NULL)
42  {
43  lines = false;
44  pmode = false;
45  num_stream = 0;
46  root_x_min = 1e100;
47  root_y_min = 1e100;
48  root_x_max = -1e100;
49  root_y_max = -1e100;
50  root = NULL;
51  }
52 
53  void StreamView::show(MeshFunction<double>* xsln, MeshFunction<double>* ysln, int marker, double step, double eps)
54  {
55  if(this->vec == NULL)
56  this->vec = new Vectorizer;
57  if(xsln == ysln)
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);
60  }
61 
62  bool StreamView::is_in_triangle(int idx, double x, double y, double3& bar)
63  {
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))
76  return true;
77  else
78  return false;
79  }
80 
81  void StreamView::add_element_to_tree(Node* father, int e_idx, double x_min, double x_max, double y_min, double y_max)
82  {
83  double4* vert = vec->get_vertices();
84  int3* xtris = vec->get_triangles();
85  if(father->leaf == true)
86  {
87  father->elements[father->num_elem++] = e_idx;
88  if(father->num_elem >= 100) // too many elements
89  {
90  father->leaf = false;
91  for (int k = 0; k < 2; k++)
92  {
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;
97  }
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);
100  }
101  }
102  else
103  {
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) // level = 1, 3, 5, ...
108  {
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);
113  }
114  else // level 0, 2, 4, ...
115  {
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);
120  }
121  }
122  }
123 
124  void StreamView::build_tree()
125  {
126  root->leaf = true;
127  root->level = 0;
128  root->num_elem = 0;
129  for (int i = 0; i < vec->get_num_triangles(); i++)
130  {
131  add_element_to_tree(root, i, root_x_min, root_x_max, root_y_min, root_y_max);
132  }
133  }
134 
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)
136  {
137  if(father->leaf == true)
138  {
139  double4* vert = vec->get_vertices();
140  int3* xtris = vec->get_triangles();
141  for (int idx = 0; idx < father->num_elem; idx++)
142  {
143  int i = father->elements[idx];
144  if(is_in_triangle(i, x, y, bar)) return i;
145  }
146  return -1;
147  }
148  else
149  {
150  double x_mid = (x_max + x_min)/2;
151  double y_mid = (y_max + y_min)/2;
152  if(father->level % 2) // level = 1, 3, 5, ...
153  if(y <= y_mid)
154  return find_triangle_in_tree(x, y, father->sons[0], x_min, x_max, y_min, y_mid, bar);
155  else
156  return find_triangle_in_tree(x, y, father->sons[1], x_min, x_max, y_mid, y_max, bar);
157  else // level 0, 2, 4, ...
158  if(x <= x_mid)
159  return find_triangle_in_tree(x, y, father->sons[0], x_min, x_mid, y_min, y_max, bar);
160  else
161  return find_triangle_in_tree(x, y, father->sons[1], x_mid, x_max, y_min, y_max, bar);
162  }
163  }
164 
165  bool StreamView::get_solution_values(double x, double y, double& xval, double& yval)
166  {
167  double4* vert = vec->get_vertices();
168  int3* xtris = vec->get_triangles();
169  double3 bar;
170  int e_idx;
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];
175  return true;
176  }
177 
178  void StreamView::delete_tree(Node* father)
179  {
180  if(father->leaf == false)
181  {
182  delete_tree(father->sons[0]);
183  delete_tree(father->sons[1]);
184  }
185  delete [] father;
186  }
187 
188  int StreamView::create_streamline(double x_start, double y_start, int idx)
189  {
190  double ODE_EPS = 1e-5;
191  double tau = initial_tau;
192  double x = x_start;
193  double y = y_start;
194  int k = 0;
195  const int buffer_length = 5000;
196  double2* buffer = new double2[buffer_length];
197  bool tau_ok = false;
198  bool end = false;
199  bool almost_end_of_domain = false;
200 
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;
205 
206  while(1)
207  {
208  if(get_solution_values(x, y, k1, l1) == false) // point (x, y) does not lie in the domain
209  {
210  tau = tau/2; // draw streamline to the end of the domain
211  if(tau < min_tau) break;
212  continue;
213  }
214  if(fabs(k1) / max_mag < 1e-5 && fabs(l1) / max_mag < 1e-5) break; // stop streamline when zero solution
215 
216  // add new point to steamline
217  buffer[k][0] = x;
218  buffer[k][1] = y;
219  k++;
220  if(k >= buffer_length) break;
221 
222  do
223  {
224  if(almost_end_of_domain) // draw streamline to the end of the domain
225  {
226  almost_end_of_domain = false;
227  tau = tau/2;
228  if(tau < min_tau) { end = true; break; }
229  }
230 
231  // Merson's adaptive Runge-Kutta method
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);
241 
242  // error according to Merson
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));
246  if(err < ODE_EPS)
247  {
248  tau_ok = true; x = x5; y = y5;
249  }
250  else
251  {
252  tau_ok = false; tau = tau/2;
253  if(tau < min_tau) tau = min_tau;
254  }
255 
256  if(err < ODE_EPS/32)
257  {
258  // new tau according to Merson
259  tau = 0.8 * tau * pow(ODE_EPS/err, 0.2);
260  if(tau > max_tau) tau = max_tau;
261  }
262  }
263  while (!tau_ok);
264  if(end) break; // get out from both while cycles
265  }
266 
267  streamlines[idx] = new double2[k];
268  memcpy(streamlines[idx], buffer, k*sizeof(double2));
269  delete [] buffer;
270 
271  return k;
272  }
273 
274  static double4* comp_vert;
275  static int compare(const void* p1, const void* p2)
276  {
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");
289  return 0;
290  }
291 
292  int StreamView::find_initial_edge(int num_edges, int3* edges)
293  {
294  int i, j;
295  double4* vert = vec->get_vertices();
296  for (i = 0; i < num_edges; i++)
297  {
298  if(edges[i][2] == 0) // not visited yet
299  {
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])
302  break;
303  if(j == num_edges) return i;
304  }
305  }
306  return -1;
307  }
308 
309  static int find_next_edge(int num_edges, int3* edges, int b_idx)
310  {
311  int3 key;
312  key[0] = b_idx;
313  int3* edge = (int3*) bsearch(&key, edges, num_edges, sizeof(int3), compare);
314  if(edge == NULL)
315  return -1; // not found
316  else
317  return edge - edges;
318  }
319 
320  void StreamView::find_initial_points(int marker, double step, double2*& initial_points)
321  {
322  int k = 0;
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++)
329  {
330  if(edge_markers[i] == marker)
331  {
332  bnd_edges[k][0] = edges[i][0];
333  bnd_edges[k][1] = edges[i][1];
334  bnd_edges[k++][2] = 0; // not visited
335  }
336  }
337  int num_edges = k;
338 
339  // sort edges by their first vertex coordinates
340  // sort first by x, then by y
341  comp_vert = vert;
342  qsort(bnd_edges, num_edges, sizeof(int3), compare);
343 
344  // create list of initial points
345  int buffer_length = 1000;
346  initial_points = new double2[buffer_length];
347  k = 0;
348  int idx;
349  while ((idx = find_initial_edge(num_edges, bnd_edges)) != -1)
350  {
351  double tmp_step = step;
352  do
353  {
354  bnd_edges[idx][2] = 1; // visited
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)
360  {
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;
364  tmp_step = step;
365  init_x = initial_points[k][0];
366  init_y = initial_points[k][1];
367  k++;
368  }
369  tmp_step = tmp_step - remaining_len;
370  }
371  while ((idx = find_next_edge(num_edges, bnd_edges, bnd_edges[idx][1])) != -1);
372  }
373  num_stream = k;
374 
375  delete [] bnd_edges;
376  }
377 
378  void StreamView::show(MeshFunction<double>* xsln, MeshFunction<double>* ysln, int marker, double step, double eps, int xitem, int yitem)
379  {
380  if(vec == NULL)
381  vec = new Vectorizer;
382  vec->process_solution(xsln, ysln, xitem, yitem, eps);
383 
384  vec->lock_data();
385  if(range_auto)
386  {
387  range_min = vec->get_min_value();
388  range_max = vec->get_max_value();
389  }
390  vec->calc_vertices_aabb(&vertices_min_x, &vertices_max_x, &vertices_min_y, &vertices_max_y);
391 
392  // create streamlines
393  double4* vert = vec->get_vertices();
394  for (int i = 0; i < vec->get_num_vertices(); i++)
395  {
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];
400  }
401 
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();
406 
407  this->tick();
408  root = new Node;
409  build_tree();
410 
411  double2* initial_points;
412  find_initial_points(marker, step, initial_points);
413 
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);
418 
419  delete [] initial_points;
420 
421  vec->unlock_data();
422 
423  create();
424  update_layout();
425  reset_view(false);
426  refresh();
427  wait_for_draw();
428  }
429 
430  void StreamView::add_streamline(double x, double y)
431  {
432  if(root == NULL)
433  throw Hermes::Exceptions::Exception("Function add_streamline must be called after StreamView::show().");
434  this->tick();
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);
438  num_stream++;
439  refresh();
440  }
441 
442  static int n_vert(int i) { return (i + 1) % 3; }
443  static int p_vert(int i) { return (i + 2) % 3; }
444 
445  void StreamView::on_display()
446  {
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);
452 
453  // transform all vertices
454  vec->lock_data();
455  int i;
456  int nv = vec->get_num_vertices();
457  double4* vert = vec->get_vertices();
458  double2* tvert = new double2[nv];
459 
460  for (i = 0; i < nv; i++)
461  {
462  tvert[i][0] = transform_x(vert[i][0]);
463  tvert[i][1] = transform_y(vert[i][1]);
464  }
465 
466  // value range
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);
470  // special case: constant solution
471  if(fabs(min - max) < 1e-8) { irange = 1.0; min -= 0.5; }
472 
473  // draw all triangles
474  int3* xtris = vec->get_triangles();
475 
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++)
481  {
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]);
485 
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]);
489 
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]);
493  }
494  glEnd();
495  glDisable(GL_TEXTURE_1D);
496 
497  // draw all edges
498  glColor3f(0.5, 0.5, 0.5);
499  glBegin(GL_LINES);
500  int2* edges = vec->get_edges();
501  for (i = 0; i < vec->get_num_edges(); i++)
502  {
503  if(lines || edges[i][2] != 0)
504  {
505  glVertex2d(tvert[edges[i][0]][0], tvert[edges[i][0]][1]);
506  glVertex2d(tvert[edges[i][1]][0], tvert[edges[i][1]][1]);
507  }
508  }
509  glEnd();
510 
511  // draw streamlines
512  glColor3f(0.0, 0.0, 0.0);
513  for (i = 0; i < num_stream; i++)
514  {
515  glBegin(GL_LINE_STRIP);
516  int k = 0;
517  for (int j = 0; j < streamlength[i] - 1; j++)
518  {
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]));
521  }
522  glEnd();
523  }
524 
525  delete [] tvert;
526  vec->unlock_data();
527  }
528 
529  void StreamView::on_mouse_move(int x, int y)
530  {
531  View::on_mouse_move(x, y);
532  }
533 
534  void StreamView::on_key_down(unsigned char key, int x, int y)
535  {
536  switch (key)
537  {
538  case 'm':
539  lines = !lines;
540  refresh();
541  break;
542 
543  case 'l':
544  pmode = !pmode;
545  refresh();
546  break;
547 
548  case 'c':
549  reset_view(true);
550  refresh();
551  break;
552 
553  case 'f':
554  set_palette_filter(pal_filter != GL_LINEAR);
555  break;
556 
557  // delete last streamline
558  case 26: // ctrl z
559  if(num_stream > 0)
560  {
561  num_stream--;
562  delete [] streamlines[num_stream];
563  refresh();
564  }
565  break;
566 
567  default:
568  View::on_key_down(key, x, y);
569  break;
570  }
571  }
572 
573  void StreamView::on_left_mouse_down(int x, int y)
574  {
575  View::on_left_mouse_down(x, y);
576 
577  // adding streamline (initial point set at (x, y))
578  if(!scale_focused && glutGetModifiers() == GLUT_ACTIVE_CTRL)
579  {
580  this->tick();
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);
584  num_stream++;
585  refresh();
586  }
587  }
588 
589  const char* StreamView::get_help_text() const
590  {
591  return
592  "StreamView\n\n"
593  "Controls:\n"
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"
601  " M - toggle mesh\n"
602  " P - cycle palettes\n"
603  " S - save screenshot\n"
604  " F1 - this help\n"
605  " Esc, Q - quit";
606  }
607 
608  StreamView::~StreamView()
609  {
610  delete_tree(root);
611  for (int i = 0; i < num_stream; i++)
612  delete [] streamlines[i];
613  delete [] streamlines;
614  delete [] streamlength;
615  delete vec;
616  }
617  }
618  }
619 }
620 
621 #endif