Hermes2D  2.0
vector_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 
18 #include <GL/freeglut.h>
19 #include "global.h"
20 #include "vector_view.h"
21 
22 namespace Hermes
23 {
24  namespace Hermes2D
25  {
26  namespace Views
27  {
28  VectorView::VectorView(const char* title, WinGeom* wg)
29  : View(title, wg)
30  {
31  vec = new Vectorizer;
32  gx = gy = 0.0;
33  gs = 20.0;
34  hexa = true;
35  mode = 0;
36  lines = false;
37  pmode = false;
38  length_coef = 1.0;
39  }
40 
41  VectorView::VectorView(char* title, WinGeom* wg)
42  : View(title, wg)
43  {
44  vec = new Vectorizer;
45  gx = gy = 0.0;
46  gs = 20.0;
47  hexa = true;
48  mode = 0;
49  lines = false;
50  pmode = false;
51  length_coef = 1.0;
52  }
53 
54  VectorView::~VectorView()
55  {
56  delete vec;
57  }
58 
59  void VectorView::show(MeshFunction<double>* vsln, double eps)
60  {
61  if(vec == NULL)
62  vec = new Vectorizer;
63  if(vsln->get_num_components() < 2)
64  throw Hermes::Exceptions::Exception("The single-argument version of show() is only for vector-valued solutions.");
65  show(vsln, vsln, eps, H2D_FN_VAL_0, H2D_FN_VAL_1);
66  }
67 
68  void VectorView::show(MeshFunction<double>* xsln, MeshFunction<double>* ysln, double eps)
69  {
70  if(vec == NULL)
71  vec = new Vectorizer;
72  if(xsln == ysln)
73  this->warn("Identical solutions passed to the two-argument version of show(). Most likely this is a mistake.");
74  show(xsln, ysln, eps, H2D_FN_VAL_0, H2D_FN_VAL_0);
75  }
76 
77  void VectorView::show(MeshFunction<double>* xsln, MeshFunction<double>* ysln, double eps, int xitem, int yitem)
78  {
79  if(vec == NULL)
80  vec = new Vectorizer;
81 
82  vec->lock_data();
83  vec->process_solution(xsln, ysln, xitem, yitem, eps);
84  if(range_auto)
85  {
86  range_min = vec->get_min_value();
87  range_max = vec->get_max_value();
88  }
89  vec->calc_vertices_aabb(&vertices_min_x, &vertices_max_x, &vertices_min_y, &vertices_max_y);
90  vec->unlock_data();
91 
92  create();
93  update_layout();
94  reset_view(false);
95 
96  refresh();
97 
98  wait_for_draw();
99  }
100 
101  static int n_vert(int i) { return (i + 1) % 3; }
102  static int p_vert(int i) { return (i + 2) % 3; }
103 
104  void VectorView::set_mode(int mode)
105  {
106  this->mode = mode % 3;
107  refresh();
108  }
109 
110  Vectorizer* VectorView::get_vectorizer()
111  {
112  return this->vec;
113  }
114 
115  void VectorView::plot_arrow(double x, double y, double xval, double yval, double max, double min, double gs)
116  {
117  if(mode == 1)
118  glColor3f(0.0f, 0.0f, 0.0f);
119  else
120  glColor3f(0.5f, 0.5f, 0.5f);
121 
122  // magnitude
123  double Real_mag = sqrt(sqr(xval) + sqr(yval));
124  double mag = Real_mag;
125  if(Real_mag > max) mag = max;
126  double length = mag/max * gs * length_coef;
127  double width = 0.1 * gs;
128  if(mode == 1) width *= 1.2;
129  double xnew = x + gs * xval * mag / (max * Real_mag) * length_coef;
130  double ynew = y - gs * yval * mag / (max * Real_mag) * length_coef;
131 
132  if((mag)/(max - min) < 1e-5)
133  {
134  glTranslated(x, y, 0.0);
135 
136  glBegin(GL_QUADS);
137  glVertex2d( width, width);
138  glVertex2d( width, -width);
139  glVertex2d(-width, -width);
140  glVertex2d(-width, width);
141  glEnd();
142  }
143  else
144  {
145  glBegin(GL_LINES);
146  glVertex2d(x, y);
147  glVertex2d(xnew, ynew);
148  glEnd();
149 
150  glTranslated(x, y, 0.0);
151  glRotated(atan2(-yval, xval) * 180.0/M_PI, 0.0, 0.0, 1.0);
152 
153  glBegin(GL_TRIANGLES);
154  glVertex2d(length + 3 * width, 0.0);
155  glVertex2d(length - 2 * width, width);
156  glVertex2d(length - 2 * width, -width);
157  glEnd();
158  }
159 
160  glLoadIdentity();
161 
162  if(mode == 1)
163  {
164  float color[3];
165  get_palette_color((mag - min)/(max - min), color); // 0.0 -- 1.0
166  glColor3f(color[0], color[1], color[2]);
167 
168  if(mag/(max - min) < 1e-5)
169  {
170  glBegin(GL_QUADS);
171  glVertex2d( width, width);
172  glVertex2d( width, -width);
173  glVertex2d(-width, -width);
174  glVertex2d(-width, width);
175  glEnd();
176  }
177  else
178  {
179  glBegin(GL_LINES);
180  glVertex2d(x, y);
181  glVertex2d(xnew, ynew);
182  glEnd();
183 
184  glTranslated(x - 1, y, 0.0);
185  glRotated(atan2(-yval, xval) * 180.0/M_PI, 0.0, 0.0, 1.0);
186 
187  glBegin(GL_TRIANGLES);
188  glVertex2d(length + 3 * width, 0.0);
189  glVertex2d(length - 2 * width, width);
190  glVertex2d(length - 2 * width, -width);
191  glEnd();
192 
193  glLoadIdentity();
194  }
195  }
196  }
197 
198  void VectorView::on_display()
199  {
200  set_ortho_projection();
201  glDisable(GL_LIGHTING);
202  glDisable(GL_DEPTH_TEST);
203  glDisable(GL_TEXTURE_1D);
204  glPolygonMode(GL_FRONT_AND_BACK, pmode ? GL_LINE : GL_FILL);
205 
206  // initial grid point and grid step
207  double gt = gs;
208  if(hexa) gt *= sqrt(3.0)/2.0;
209 
210  double max_length = 0.0;
211 
212  // transform all vertices
213  vec->lock_data();
214  int i;
215  int nv = vec->get_num_vertices();
216  double4* vert = vec->get_vertices();
217  double2* tvert = new double2[nv];
218 
219  for (i = 0; i < nv; i++)
220  {
221  tvert[i][0] = transform_x(vert[i][0]);
222  tvert[i][1] = transform_y(vert[i][1]);
223 
224  // find max length of vectors
225  double length = sqr(vert[i][2]) + sqr(vert[i][3]);
226  if(length > max_length) max_length = length;
227  }
228  max_length = sqrt(max_length);
229 
230  // value range
231  double min = range_min, max = range_max;
232  if(range_auto) { min = vec->get_min_value(); max = vec->get_max_value(); }
233  double irange = 1.0 / (max - min);
234  // special case: constant solution
235  if(fabs(min - max) < 1e-8) { irange = 1.0; min -= 0.5; }
236 
237  // draw all triangles
238  int3* xtris = vec->get_triangles();
239 
240  if(mode != 1) glEnable(GL_TEXTURE_1D);
241  glBindTexture(GL_TEXTURE_1D, gl_pallete_tex_id);
242  glBegin(GL_TRIANGLES);
243  glColor3f(0.95f, 0.95f, 0.95f);
244  for (i = 0; i < vec->get_num_triangles(); i++)
245  {
246  double mag = sqrt(sqr(vert[xtris[i][0]][2]) + sqr(vert[xtris[i][0]][3]));
247  glTexCoord2d((mag -min) * irange * tex_scale + tex_shift, 0.0);
248  glVertex2d(tvert[xtris[i][0]][0], tvert[xtris[i][0]][1]);
249 
250  mag = sqrt(sqr(vert[xtris[i][1]][2]) + sqr(vert[xtris[i][1]][3]));
251  glTexCoord2d((mag -min) * irange * tex_scale + tex_shift, 0.0);
252  glVertex2d(tvert[xtris[i][1]][0], tvert[xtris[i][1]][1]);
253 
254  mag = sqrt(sqr(vert[xtris[i][2]][2]) + sqr(vert[xtris[i][2]][3]));
255  glTexCoord2d((mag -min) * irange * tex_scale + tex_shift, 0.0);
256  glVertex2d(tvert[xtris[i][2]][0], tvert[xtris[i][2]][1]);
257  }
258  glEnd();
259  glDisable(GL_TEXTURE_1D);
260 
261  // draw all edges
262  /*if(mode == 0) glColor3f(0.3, 0.3, 0.3);
263  else*/ glColor3f(0.5, 0.5, 0.5);
264  glBegin(GL_LINES);
265  int2* edges = vec->get_edges();
266  for (i = 0; i < vec->get_num_edges(); i++)
267  {
268  if(lines || edges[i][2] != 0)
269  {
270  glVertex2d(tvert[edges[i][0]][0], tvert[edges[i][0]][1]);
271  glVertex2d(tvert[edges[i][1]][0], tvert[edges[i][1]][1]);
272  }
273  }
274  glEnd();
275 
276  // draw dashed edges
277  if(lines)
278  {
279  glEnable(GL_LINE_STIPPLE);
280  glLineStipple(1, 0xCCCC);
281  glBegin(GL_LINES);
282  int2* dashes = vec->get_dashes();
283  for (i = 0; i < vec->get_num_dashes(); i++)
284  {
285  glVertex2d(tvert[dashes[i][0]][0], tvert[dashes[i][0]][1]);
286  glVertex2d(tvert[dashes[i][1]][0], tvert[dashes[i][1]][1]);
287  }
288  glEnd();
289  glDisable(GL_LINE_STIPPLE);
290  }
291 
292  // draw arrows
293  if(mode != 2)
294  {
295  for (i = 0; i < vec->get_num_triangles(); i++)
296  {
297  double miny = 1e100;
298  int idx = -1, k, l1, l2, r2, r1, s;
299  double lry, x;
300  double mr, ml, lx, rx, xval, yval;
301 
302  double wh = output_height + gt, ww = output_width + gs;
303  if((tvert[xtris[i][0]][0] < -gs) && (tvert[xtris[i][1]][0] < -gs) && (tvert[xtris[i][2]][0] < -gs)) continue;
304  if((tvert[xtris[i][0]][0] > ww) && (tvert[xtris[i][1]][0] > ww) && (tvert[xtris[i][2]][0] > ww)) continue;
305  if((tvert[xtris[i][0]][1] < -gt) && (tvert[xtris[i][1]][1] < -gt) && (tvert[xtris[i][2]][1] < -gt)) continue;
306  if((tvert[xtris[i][0]][1] > wh) && (tvert[xtris[i][1]][1] > wh) && (tvert[xtris[i][2]][1] > wh)) continue;
307 
308  // find vertex with min y-coordinate
309  for (k = 0; k < 3; k++)
310  if(tvert[xtris[i][k]][1] < miny)
311  miny = tvert[xtris[i][idx = k]][1];
312  l1 = r1 = xtris[i][idx];
313  l2 = xtris[i][n_vert(idx)];
314  r2 = xtris[i][p_vert(idx)];
315 
316  // plane of x and y values on triangle
317  double a[2], b[2], c[2], d[2];
318  for (int n = 0; n < 2; n++)
319  {
320  a[n] = (tvert[l1][1] - tvert[l2][1])*(vert[r1][2 +n] - vert[r2][2 + n]) - (vert[l1][2 + n] - vert[l2][2 + n])*(tvert[r1][1] - tvert[r2][1]);
321  b[n] = (vert[l1][2 + n] - vert[l2][2 + n])*(tvert[r1][0] - tvert[r2][0]) - (tvert[l1][0] - tvert[l2][0])*(vert[r1][2 + n] - vert[r2][2 + n]);
322  c[n] = (tvert[l1][0] - tvert[l2][0])*(tvert[r1][1] - tvert[r2][1]) - (tvert[l1][1] - tvert[l2][1])*(tvert[r1][0] - tvert[r2][0]);
323  d[n] = -a[n] * tvert[l1][0] - b[n] * tvert[l1][1] - c[n] * vert[l1][2 + n];
324  a[n] /= c[n]; b[n] /= c[n]; d[n] /= c[n];
325  }
326 
327  s = (int) ceil((tvert[l1][1] - gy)/gt); // first step
328  lry = gy + s*gt;
329  bool shift = hexa && (s & 1);
330 
331  // if there are two points with min y-coordinate, switch to the next segment
332  if((tvert[l1][1] == tvert[l2][1]) || (tvert[r1][1] == tvert[r2][1]))
333  {
334  if(tvert[l1][1] == tvert[l2][1])
335  {l1 = l2; l2 = r2;}
336  else if(tvert[r1][1] == tvert[r2][1])
337  {r1 = r2; r2 = l2;}
338  }
339 
340  // slope of the left and right segment
341  ml = (tvert[l1][0] - tvert[l2][0])/(tvert[l1][1] - tvert[l2][1]);
342  mr = (tvert[r1][0] - tvert[r2][0])/(tvert[r1][1] - tvert[r2][1]);
343  // x-coordinates of the endpoints of the first line
344  lx = tvert[l1][0] + ml * (lry - (tvert[l1][1]));
345  rx = tvert[r1][0] + mr * (lry - (tvert[r1][1]));
346 
347  if(lry < -gt)
348  {
349  k = (int) floor(-lry/gt);
350  lry += gt * k;
351  lx += k * ml * gt;
352  rx += k * mr * gt;
353  }
354 
355  // while we are in triangle
356  while (((lry < tvert[l2][1]) || (lry < tvert[r2][1])) && (lry < wh))
357  {
358  // while we are in the segment
359  while (((lry <= tvert[l2][1]) && (lry <= tvert[r2][1])) && (lry < wh))
360  {
361  double gz = gx;
362  if(shift) gz -= 0.5*gs;
363  s = (int) ceil((lx - gz)/gs);
364  x = gz + s*gs;
365  if(hexa) shift = !shift;
366 
367  if(x < -gs)
368  {
369  k = (int) floor(-x/gs);
370  x += gs * k;
371  }
372  // go along the line
373  while ((x < rx) && (x < ww))
374  {
375  // plot the arrow
376  xval = -a[0]*x - b[0]*lry - d[0];
377  yval = -a[1]*x - b[1]*lry - d[1];
378  plot_arrow(x, lry, xval, yval, max, min, gs);
379  x += gs;
380  }
381  // move to the next line
382  lx += ml*gt;
383  rx += mr*gt;
384  lry += gt;
385  }
386  // change segment
387  if(lry >= tvert[l2][1])
388  {
389  l1 = l2; l2 = r2;
390  ml = (tvert[l1][0] - tvert[l2][0])/(tvert[l1][1] - tvert[l2][1]);
391  lx = tvert[l1][0] + ml * (lry - (tvert[l1][1]));
392  }
393  else
394  {
395  r1 = r2; r2 = l2;
396  mr = (tvert[r1][0] - tvert[r2][0])/(tvert[r1][1] - tvert[r2][1]);
397  rx = tvert[r1][0] + mr * (lry - (tvert[r1][1]));
398  }
399  }
400  }
401  }
402 
403  delete [] tvert;
404  vec->unlock_data();
405  }
406 
407  void VectorView::on_mouse_move(int x, int y)
408  {
409  if(dragging)
410  {
411  gx += (x - mouse_x);
412  gy += (y - mouse_y);
413  }
414  View::on_mouse_move(x, y);
415  }
416 
417  void VectorView::on_key_down(unsigned char key, int x, int y)
418  {
419  switch (key)
420  {
421  case 'm':
422  lines = !lines;
423  refresh();
424  break;
425 
426  case 'l':
427  pmode = !pmode;
428  refresh();
429  break;
430 
431  case 'c':
432  reset_view(true);
433  refresh();
434  break;
435 
436  case 'f':
437  set_palette_filter(pal_filter != GL_LINEAR);
438  break;
439 
440  case 'x':
441  set_grid_type(!hexa);
442  break;
443 
444  case 'b':
445  mode++;
446  if(mode > 2) mode = 0;
447  refresh();
448  break;
449 
450  case '*':
451  case '/':
452  if(key == '*') length_coef *= 1.1; else length_coef /= 1.1;
453  refresh();
454  break;
455 
456  default:
457  View::on_key_down(key, x, y);
458  break;
459  }
460  }
461 
462  const char* VectorView::get_help_text() const
463  {
464  return
465  "VectorView\n\n"
466  "Controls:\n"
467  " Left mouse - pan\n"
468  " Right mouse - zoom\n"
469  " B - toggle view mode (type of arrows x no arrows)\n"
470  " * - extend arrows\n"
471  " / - shorten arrows\n"
472  " C - center image\n"
473  " F - toggle smooth palette\n"
474  " X - toggle hexagonal grid\n"
475  " H - render high-quality frame\n"
476  " M - toggle mesh\n"
477  " P - cycle palettes\n"
478  " S - save screenshot\n"
479  " F1 - this help\n"
480  " Esc, Q - quit";
481  }
482  }
483  }
484 }
485 #endif