Hermes2D  2.0
scalar_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/glew.h>
19 #include <GL/freeglut.h>
20 #include <algorithm>
21 #include <cmath>
22 #include <list>
23 #include "global.h"
24 #include "scalar_view.h"
25 
26 #define GL_BUFFER_OFFSET(i) ((char *)NULL + (i))
27 
28 /* constants */
29 #define MIN_CONT_STEP 1.0E-2
30 #define CONT_CHANGE 1.0E-2
31 #define D3DV_SCALE_STEP_COEF 1.1
32 
33 namespace Hermes
34 {
35  namespace Hermes2D
36  {
37  namespace Views
38  {
39  const int ScalarView::fovy = 50;
40  const double ScalarView::znear = 0.05;
41  const double ScalarView::zfar = 10;
42 
43  void ScalarView::init()
44  {
45  lin = new Linearizer;
46  pmode = mode3d = false;
47  normals = NULL;
48  panning = false;
49 
50  contours = false;
51  cont_orig = 0.0;
52  cont_step = 0.2;
53  cont_color[0] = 0.0f; cont_color[1] = 0.0f; cont_color[2] = 0.0f;
54 
55  show_edges = true;
56  show_aabb = false;
57  edges_color[0] = 0.5f; edges_color[1] = 0.4f; edges_color[2] = 0.4f;
58 
59  show_values = true;
60  lin_updated = false;
61  gl_coord_buffer = 0; gl_index_buffer = 0; gl_edge_inx_buffer = 0;
62 
63  do_zoom_to_fit = true;
64  is_constant = false;
65  }
66 
67 #ifndef _MSC_VER
68 
69  ScalarView::ScalarView(const char* title, WinGeom* wg) :
70  View(title, wg),
71  vertex_nodes(0),
72  pointed_vertex_node(NULL),
73  allow_node_selection(false),
74  pointed_node_widget(0),
75  selected_node_widget(0),
76  node_pixel_radius(10),
77  node_widget_vert_cnt(32),
78  element_id_widget(0),
79  show_element_info(false)
80  {
81  init();
82  }
83 #else
84  ScalarView::ScalarView(WinGeom* wg) :
85  View("ScalarView", wg),
86  vertex_nodes(0),
87  pointed_vertex_node(NULL),
88  allow_node_selection(false),
89  pointed_node_widget(0),
90  selected_node_widget(0),
91  node_pixel_radius(10),
92  node_widget_vert_cnt(32),
93  element_id_widget(0),
94  show_element_info(false)
95  {
96  init();
97  }
98 #endif
99 
100  ScalarView::ScalarView(char* title, WinGeom* wg) :
101  View(title, wg),
102  vertex_nodes(0),
103  pointed_vertex_node(NULL),
104  allow_node_selection(false),
105  pointed_node_widget(0),
106  selected_node_widget(0),
107  node_pixel_radius(10),
108  node_widget_vert_cnt(32),
109  element_id_widget(0),
110  show_element_info(false)
111  {
112  init();
113  }
114 
115  ScalarView::~ScalarView()
116  {
117  delete [] normals;
118  vertex_nodes.clear();
119  delete lin;
120  }
121 
122  void ScalarView::on_close()
123  {
124  //OpenGL clenaup
125  if(pointed_node_widget != 0)
126  {
127  glDeleteLists(pointed_node_widget, 1);
128  pointed_node_widget = 0;
129  }
130  if(selected_node_widget != 0)
131  {
132  glDeleteLists(selected_node_widget, 1);
133  selected_node_widget = 0;
134  }
135  if(element_id_widget != 0)
136  {
137  glDeleteLists(element_id_widget, 1);
138  element_id_widget = 0;
139  }
140  if(gl_coord_buffer != 0)
141  {
142  glDeleteBuffersARB(1, &gl_coord_buffer);
143  gl_coord_buffer = 0;
144  }
145  if(gl_index_buffer != 0)
146  {
147  glDeleteBuffersARB(1, &gl_index_buffer);
148  gl_index_buffer = 0;
149  }
150 
151  //call of parent implementation
152  View::on_close();
153  }
154 
155  void ScalarView::show(MeshFunction<double>* sln, double eps, int item,
156  MeshFunction<double>* xdisp, MeshFunction<double>* ydisp, double dmult)
157  {
158  // For preservation of the sln's active element. Will be set back after the visualization.
159  Element* active_element = sln->get_active_element();
160 
161  if(!range_auto)
162  lin->set_max_absolute_value(std::max(fabs(range_min), fabs(range_max)));
163 
164  lin->set_displacement(xdisp, ydisp, dmult);
165  lin->lock_data();
166 
167  lin->process_solution(sln, item, eps);
168 
169  update_mesh_info();
170 
171  // Initialize mesh nodes for displaying and selection.
172  init_vertex_nodes(sln->get_mesh());
173 
174  // Initialize element info.
175  init_element_info(sln->get_mesh());
176 
177  lin->unlock_data();
178 
179  create();
180  update_layout();
181  wait_for_draw();
182 
183  // FIXME: find out why this has to be called after wait_for_draw in order for the view to be reset initially.
184  reset_view(false); // setting true here makes the view always reset after calling 'show'; particularly in the adaptivity process,
185  // it would disallow the observation of the process from a manually set viewpoint.
186  refresh();
187 
188  // Now we reset the active element if it was set before the MeshFunction sln entered this method.
189  // Only for Solutions. This method may fail for filters, as they may not have RefMaps correctly set.
190  if(dynamic_cast<Solution<double>*>(sln) != NULL)
191  if(active_element != NULL)
192  // Also when there has not been a call to set_active_element since assignment to this MeshFunction,
193  // there is nothing to restore to.
194  if(active_element->active)
195  sln->set_active_element(active_element);
196  }
197 
198  void ScalarView::show_linearizer_data(double eps, int item)
199  {
200  update_mesh_info();
201 
202  create();
203  update_layout();
204  wait_for_draw();
205  // FIXME: find out why this has to be called after wait_for_draw in order for the view to be reset initially.
206  reset_view(false); // setting true here makes the view always reset after calling 'show'; particularly in the adaptivity process,
207  // it would disallow the observation of the process from a manually set viewpoint.
208  refresh();
209  }
210 
211  void ScalarView::update_mesh_info()
212  {
213  // Calculate normals if necessary.
214  if(mode3d)
215  calculate_normals(lin->get_vertices(), lin->get_num_vertices(), lin->get_triangles(), lin->get_num_triangles());
216  else
217  {
218  delete [] normals;
219  normals = NULL;
220  }
221 
222  // Get a range of vertex values (or use the range set by the user).
223  double vert_min = lin->get_min_value();
224  double vert_max = lin->get_max_value();
225  // Special case: constant function; offset the lower limit of range so that the domain is drawn under the
226  // function and also the scale is drawn correctly.
227  if((vert_max - vert_min) < 1e-8)
228  {
229  is_constant = true;
230  vert_min -= 0.5;
231  }
232 
233  if(range_auto)
234  {
235  range_min = vert_min;
236  range_max = vert_max;
237  }
238 
239  if(fabs(range_max - range_min) < 1e-8)
240  value_irange = 1.0;
241  else
242  value_irange = 1.0 / (range_max - range_min);
243 
244  // Calculate the axes-aligned bounding box in the xy-plane.
245  lin->calc_vertices_aabb(&vertices_min_x, &vertices_max_x, &vertices_min_y, &vertices_max_y);
246 
247  // Calculate average value.
248  value_range_avg = 0.0;
249  double3* verts = lin->get_vertices();
250  const int num_verts = lin->get_num_vertices();
251  for(int i = 0; i < num_verts; i++)
252  if(verts[i][2] > range_max)
253  value_range_avg += range_max;
254  else if(verts[i][2] < range_min)
255  value_range_avg += range_min;
256  else
257  value_range_avg += verts[i][2];
258  value_range_avg /= num_verts;
259 
260  lin_updated = true;
261  }
262 
263  bool ScalarView::compare_vertex_nodes_x(const VertexNodeInfo& a, const VertexNodeInfo& b)
264  {
265  return a.x < b.x;
266  }
267 
268  void ScalarView::init_vertex_nodes(const Mesh* mesh)
269  {
270  //clear all selections
271  pointed_vertex_node = NULL;
272  vertex_nodes.clear();
273 
274  //count a number of active nodes
275  int active_nodes = 0;
276  int max_node_id = mesh->get_max_node_id();
277  for(int i = 0; i < max_node_id; i++)
278  {
279  Node* node = mesh->get_node(i);
280  if(node->type == 0 && node->used) //inspect only vertex nodes
281  active_nodes++;
282  }
283 
284  //allocate
285  vertex_nodes.clear();
286  vertex_nodes.reserve(active_nodes);
287 
288  //copy
289  for(int i = 0; i < max_node_id; i++)
290  {
291  Node* node = mesh->get_node(i);
292  if(node->type == 0 && node->used)
293  vertex_nodes.push_back(VertexNodeInfo(node->id, (float)node->x, (float)node->y));
294  }
295 
296  //sort
297  std::sort(vertex_nodes.begin(), vertex_nodes.end(), compare_vertex_nodes_x);
298  }
299 
300  ScalarView::VertexNodeInfo* ScalarView::find_nearest_node_in_range(float x, float y, float radius)
301  {
302  VertexNodeInfo node_info(-1, x - radius, y); //right side of the widget
303  Hermes::vector<VertexNodeInfo>::iterator found_iter = std::lower_bound(vertex_nodes.begin(), vertex_nodes.end(), node_info, compare_vertex_nodes_x);
304  Hermes::vector<VertexNodeInfo>::iterator found_nearest = vertex_nodes.end();
305  float found_nearest_dist = -1;
306  while (found_iter != vertex_nodes.end() && abs(found_iter->x - x) <= radius)
307  {
308  if(abs(found_iter->y - y) <= radius)
309  {
310  float dist = std::min(found_iter->x - x, found_iter->y - y);
311  if(found_nearest_dist < 0 || dist < found_nearest_dist)
312  {
313  found_nearest_dist = dist;
314  found_nearest = found_iter;
315  }
316  }
317  ++found_iter;
318  }
319 
320  if(found_nearest != vertex_nodes.end())
321  return found_nearest.operator->();
322  else
323  return NULL;
324  }
325 
326  void ScalarView::draw_single_vertex_node(const VertexNodeInfo& node)
327  {
328  //prepare environment
329  glPushMatrix();
330  glTranslatef(node.x, node.y, 0.0f);
331  glScalef(1/(float)scale, 1/(float)scale, 1.0f);
332 
333  //prepare number
334  void* font = GLUT_BITMAP_HELVETICA_10;
335  unsigned char buffer[128];
336  sprintf((char*)buffer, "%d", node.id);
337  int width_px = glutBitmapLength(font, buffer);
338  int height_px = glutBitmapHeight(font);
339 
340  //draw target
341  if(width_px > (2*node_pixel_radius))
342  {
343  glPushMatrix();
344  float coef = width_px / (2.0f * node_pixel_radius);
345  glScalef(coef, coef, 1.0f);
346  glCallList(selected_node_widget);
347  glPopMatrix();
348  }
349  else
350  glCallList(selected_node_widget);
351 
352  //draw number
353  glTranslatef(-width_px/2.0f, -height_px/3.0f, 0.0f);
354  glColor3f(1.0f, 1.0f, 1.0f);
355  glRasterPos2f(0.0f, 0.0f);
356  glutBitmapString(font, buffer);
357 
358  //clear environment
359  glPopMatrix();
360  }
361 
362  void ScalarView::draw_vertex_nodes()
363  {
364  //create widgets
365  create_nodes_widgets();
366 
367  //prepare environment
368  glMatrixMode(GL_MODELVIEW);
369  glDisable(GL_TEXTURE_1D);
370  glDisable(GL_LIGHTING);
371 
372  //draw selected nodes
373  glDisable(GL_BLEND);
374  Hermes::vector<VertexNodeInfo>::const_iterator iter = vertex_nodes.begin();
375  while (iter != vertex_nodes.end())
376  {
377  if(iter->selected)
378  draw_single_vertex_node(*iter);
379  ++iter;
380  }
381 
382  //draw a node under cursor
383  if(pointed_vertex_node != NULL)
384  {
385  glEnable(GL_BLEND);
386  glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
387  glPushMatrix();
388  glTranslatef(pointed_vertex_node->x, pointed_vertex_node->y, 0.0f);
389  glScalef(1/(float)scale, 1/(float)scale, 1.0f);
390  glCallList(pointed_node_widget);
391  glPopMatrix();
392  glDisable(GL_BLEND);
393  }
394  }
395 
396  void ScalarView::create_nodes_widgets()
397  {
398  //pointed node
399  if(pointed_node_widget == 0)
400  {
401  pointed_node_widget = glGenLists(1);
402  glNewList(pointed_node_widget, GL_COMPILE);
403  {
404  //background
405  glColor4f(1.0f, 1.0f, 1.0f, 0.5f);
406  glBegin(GL_TRIANGLE_FAN);
407  glVertex2f(0.0f, 0.0f);
408  float radius = 1.3f * node_pixel_radius;
409  for(int i = 0; i < node_widget_vert_cnt; i++)
410  {
411  float angle = (float)((i / (float)(node_widget_vert_cnt-1)) * (2.0f * M_PI));
412  glVertex2f(radius * cos(angle), radius * sin(angle));
413  }
414  glEnd();
415 
416  //foreground
417  glColor4f(1.0f, 0.0f, 0.0f, 1.0f);
418  glBegin(GL_LINE_STRIP);
419  radius = (float)node_pixel_radius;
420  for(int i = 0; i < node_widget_vert_cnt; i++)
421  {
422  float angle = (float)((i / (float)(node_widget_vert_cnt-1)) * (2.0f * M_PI));
423  glVertex2f(radius * cos(angle), radius * sin(angle));
424  }
425  glEnd();
426  }
427  glEndList();
428  }
429 
430  //selected mesh node
431  if(selected_node_widget == 0)
432  {
433  selected_node_widget = glGenLists(1);
434  glNewList(selected_node_widget, GL_COMPILE);
435  {
436  //background
437  glColor4f(1.0f, 1.0f, 1.0f, 1.0f);
438  glBegin(GL_TRIANGLE_FAN);
439  glVertex2f(0.0f, 0.0f);
440  float radius = 1.3f * node_pixel_radius;
441  for(int i = 0; i < node_widget_vert_cnt; i++)
442  {
443  float angle = (float)((i / (float)(node_widget_vert_cnt-1)) * (2.0f * M_PI));
444  glVertex2f(radius * cos(angle), radius * sin(angle));
445  }
446  glEnd();
447 
448  //foreground
449  glColor4f(0.4f, 0.4f, 0.2f, 1.0f);
450  glBegin(GL_TRIANGLE_FAN);
451  glVertex2f(0.0f, 0.0f);
452  radius = (float)node_pixel_radius;
453  for(int i = 0; i < node_widget_vert_cnt; i++)
454  {
455  float angle = (float)((i / (float)(node_widget_vert_cnt-1)) * (2.0f * M_PI));
456  glVertex2f(radius * cos(angle), radius * sin(angle));
457  }
458  glEnd();
459  }
460  glEndList();
461  }
462  }
463 
464  void ScalarView::init_element_info(const Mesh* mesh)
465  {
466  //cleanup
467  element_infos.clear();
468 
469  //check how many active elements are neccessary and prepare space for them
470  element_infos.reserve(mesh->get_num_active_elements());
471 
472  //build element info
473  Element *element = NULL;
474  for_all_active_elements(element, mesh)
475  {
476  double sum_x = 0.0, sum_y = 0.0;
477  double max_x, max_y, min_x, min_y;
478  max_x = min_x = element->vn[0]->x;
479  max_y = min_y = element->vn[0]->y;
480  for(unsigned int i = 0; i < element->get_nvert(); i++)
481  {
482  sum_x += element->vn[i]->x;
483  sum_y += element->vn[i]->y;
484 
485  if(element->vn[i]->x > max_x)
486  max_x = element->vn[i]->x;
487  if(element->vn[i]->x < min_x)
488  min_x = element->vn[i]->x;
489  if(element->vn[i]->y > max_y)
490  max_y = element->vn[i]->y;
491  if(element->vn[i]->y < min_y)
492  min_y = element->vn[i]->y;
493  }
494  element_infos.push_back(ElementInfo(element->id,
495  (float)(sum_x / element->get_nvert()), (float)(sum_y / element->get_nvert()),
496  (float)(max_x - min_x), (float)(max_y - min_y)));
497  }
498  }
499 
500  Linearizer* ScalarView::get_linearizer()
501  {
502  return this->lin;
503  }
504 
505  void ScalarView::draw_element_infos_2d()
506  {
507  //create widgets
508  create_element_info_widgets();
509 
510  //prepare environment
511  glMatrixMode(GL_MODELVIEW);
512  glDisable(GL_TEXTURE_1D);
513  glDisable(GL_LIGHTING);
514  glDisable(GL_BLEND);
515 
516  //draw element IDs
517  Hermes::vector<ElementInfo>::const_iterator iter = element_infos.begin();
518  while (iter != element_infos.end())
519  {
520  //check element dimension in pixels
521  float width = (float)(iter->width * scale);
522  float height = (float)(iter->height * scale);
523 
524  //draw if AABB of element is large enough
525  if(width > 6*node_pixel_radius && height > 3*node_pixel_radius)
526  {
527  //prepare environment
528  glPushMatrix();
529  glTranslatef(iter->x, iter->y, 0.0f);
530  glScalef(1/(float)scale, 1/(float)scale, 1.0f);
531 
532  //prepare number
533  void* font = GLUT_BITMAP_HELVETICA_10;
534  unsigned char buffer[128];
535  sprintf((char*)buffer, "%d", iter->id);
536  int width_px = glutBitmapLength(font, buffer);
537  int height_px = glutBitmapHeight(font);
538 
539  //draw background
540  glCallList(element_id_widget);
541 
542  //draw number
543  glTranslatef(-width_px/2.0f, -height_px/3.0f, 0.0f);
544  glColor3f(1.0f, 1.0f, 1.0f);
545  glRasterPos2f(0.0f, 0.0f);
546  glutBitmapString(font, buffer);
547 
548  //clear environment
549  glPopMatrix();
550  }
551 
552  //advance to next
553  ++iter;
554  }
555  }
556 
557  void ScalarView::create_element_info_widgets()
558  {
559  if(element_id_widget == 0)
560  {
561  element_id_widget = glGenLists(1);
562  glNewList(element_id_widget, GL_COMPILE);
563  {
564  glBegin(GL_QUADS);
565 
566  //background
567  float radius = 2.0f * node_pixel_radius;
568  glColor4f(1.0f, 1.0f, 1.0f, 1.0f);
569  glVertex2f(-radius*1.1, radius*0.5);
570  glVertex2f( radius*1.1, radius*0.5);
571  glVertex2f( radius*1.1, -radius*0.5);
572  glVertex2f(-radius*1.1, -radius*0.5);
573 
574  //foreground
575  radius = 1.8f * node_pixel_radius;
576  glColor4f(0.2f, 0.2f, 0.4f, 1.0f);
577  glVertex2f(-radius*1.1, radius*0.5);
578  glVertex2f( radius*1.1, radius*0.5);
579  glVertex2f( radius*1.1, -radius*0.5);
580  glVertex2f(-radius*1.1, -radius*0.5);
581 
582  glEnd();
583  }
584  glEndList();
585  }
586  }
587 
588  void ScalarView::show_contours(double step, double orig)
589  {
590  if(step == 0.0)
591  throw Exceptions::ValueException("step", step, 0.0);
592  contours = true;
593  cont_orig = orig;
594  cont_step = step;
595  set_palette_filter(true);
596  refresh();
597  }
598 
599  static double my_ceil(double x)
600  {
601  double y = ceil(x);
602  if(y > x) return y;
603  return y + 1.0;
604  }
605 
606  void ScalarView::draw_tri_contours(double3* vert, int3* tri)
607  {
608  // sort the vertices by their value, keep track of the permutation sign
609  int i, idx[3], perm = 0;
610  memcpy(idx, tri, sizeof(idx));
611  for (i = 0; i < 2; i++)
612  {
613  if(vert[idx[0]][2] > vert[idx[1]][2]) { std::swap(idx[0], idx[1]); perm++; }
614  if(vert[idx[1]][2] > vert[idx[2]][2]) { std::swap(idx[1], idx[2]); perm++; }
615  }
616  if(fabs(vert[idx[0]][2] - vert[idx[2]][2]) < 1e-3 * fabs(cont_step)) return;
617 
618  // get the first (lowest) contour value
619  double val = vert[idx[0]][2];
620  val = my_ceil((val - cont_orig) / cont_step) * cont_step + cont_orig;
621 
622  int l1 = 0, l2 = 1;
623  int r1 = 0, r2 = 2;
624 
625  // Adjustment of the contour step.
626  /*
627  while(std::abs(val - vert[idx[r2]][2]) > 50 * cont_step)
628  cont_step *= 2;
629  while(std::abs(val - vert[idx[r2]][2]) < 2E-2 * cont_step)
630  cont_step /= 2;
631  */
632 
633  while (val < vert[idx[r2]][2])
634  {
635  double ld = vert[idx[l2]][2] - vert[idx[l1]][2];
636  double rd = vert[idx[r2]][2] - vert[idx[r1]][2];
637 
638  // draw a slice of the triangle
639  while (val < vert[idx[l2]][2])
640  {
641  double lt = (val - vert[idx[l1]][2]) / ld;
642  double rt = (val - vert[idx[r1]][2]) / rd;
643 
644  double x1 = (1.0 - lt) * vert[idx[l1]][0] + lt * vert[idx[l2]][0];
645  double y1 = (1.0 - lt) * vert[idx[l1]][1] + lt * vert[idx[l2]][1];
646  double x2 = (1.0 - rt) * vert[idx[r1]][0] + rt * vert[idx[r2]][0];
647  double y2 = (1.0 - rt) * vert[idx[r1]][1] + rt * vert[idx[r2]][1];
648 
649  if(perm & 1) { glVertex2d(x1, y1); glVertex2d(x2, y2); }
650  else { glVertex2d(x2, y2); glVertex2d(x1, y1); }
651 
652  val += cont_step;
653  }
654  l1 = 1;
655  l2 = 2;
656  }
657  }
658 
659  void ScalarView::prepare_gl_geometry()
660  {
661  if(lin_updated)
662  {
663  lin_updated = false;
664 
665  try
666  {
667  //get input data
668  int vert_cnt = lin->get_num_vertices();
669  double3* verts = lin->get_vertices();
670  int tri_cnt = lin->get_num_triangles();
671  int3* tris = lin->get_triangles();
672 
673  //check if extension is supported
674  if(!GLEW_ARB_vertex_buffer_object)
675  throw std::runtime_error("ARB_vertex_buffer_object not supported");
676 
677  //reallocate indices
678  if(gl_index_buffer == 0 || tri_cnt > max_gl_tris)
679  {
680  if(gl_index_buffer == 0)
681  glGenBuffersARB(1, &gl_index_buffer);
682  glBindBufferARB(GL_ELEMENT_ARRAY_BUFFER_ARB, gl_index_buffer);
683  glBufferDataARB(GL_ELEMENT_ARRAY_BUFFER_ARB, sizeof(GLint) * tri_cnt * 3, NULL, GL_DYNAMIC_DRAW_ARB);
684  GLenum err = glGetError();
685  if(err != GL_NO_ERROR)
686  throw std::runtime_error("unable to allocate vertex buffer: " + err);
687  max_gl_tris = tri_cnt;
688  }
689  else
690  {
691  glBindBufferARB(GL_ELEMENT_ARRAY_BUFFER_ARB, gl_index_buffer);
692  }
693 
694  //fill indices
695  GLuint* gl_triangle = (GLuint*)glMapBufferARB(GL_ELEMENT_ARRAY_BUFFER_ARB, GL_WRITE_ONLY_ARB);
696  if(gl_triangle == NULL)
697  throw std::runtime_error("unable to map index buffer: " + glGetError());
698  gl_tri_cnt = 0;
699  for(int i = 0; i < tri_cnt; i++)
700  {
701  const int3& triangle = tris[i];
702  const double3& vert_a = verts[triangle[0]];
703  const double3& vert_b = verts[triangle[1]];
704  const double3& vert_c = verts[triangle[2]];
705  if(finite(vert_a[2]) && finite(vert_b[2]) && finite(vert_c[2]))
706  {
707  gl_triangle[0] = (GLint)triangle[0];
708  gl_triangle[1] = (GLint)triangle[1];
709  gl_triangle[2] = (GLint)triangle[2];
710  gl_tri_cnt++;
711  gl_triangle += 3; //three indices per triangle
712  }
713  }
714  glUnmapBufferARB(GL_ELEMENT_ARRAY_BUFFER_ARB);
715 
716  //reallocate vertices
717  if(gl_coord_buffer == 0 || vert_cnt > max_gl_verts)
718  {
719  if(gl_coord_buffer == 0)
720  glGenBuffersARB(1, &gl_coord_buffer);
721  glBindBufferARB(GL_ARRAY_BUFFER_ARB, gl_coord_buffer);
722  glBufferDataARB(GL_ARRAY_BUFFER_ARB, sizeof(GLVertex2) * vert_cnt, NULL, GL_DYNAMIC_DRAW_ARB);
723  GLenum err = glGetError();
724  if(err != GL_NO_ERROR)
725  throw std::runtime_error("unable to allocate coord buffer: " + err);
726  max_gl_verts = vert_cnt;
727  }
728  else
729  {
730  glBindBufferARB(GL_ARRAY_BUFFER_ARB, gl_coord_buffer);
731  }
732 
733  //fill vertices
734  GLVertex2* gl_verts = (GLVertex2*)glMapBufferARB(GL_ARRAY_BUFFER_ARB, GL_WRITE_ONLY_ARB);
735  if(gl_verts == NULL)
736  throw std::runtime_error("unable to map coord buffer: " + glGetError());
737  for(int i = 0; i < vert_cnt; i++)
738  gl_verts[i] = GLVertex2((float)verts[i][0], (float)verts[i][1], (float)((verts[i][2] - range_min) * value_irange));
739  glUnmapBufferARB(GL_ARRAY_BUFFER_ARB);
740 
741  //allocate edge indices
742  if(gl_edge_inx_buffer == 0)
743  {
744  glGenBuffersARB(1, &gl_edge_inx_buffer);
745  glBindBufferARB(GL_ARRAY_BUFFER_ARB, gl_edge_inx_buffer);
746  glBufferDataARB(GL_ARRAY_BUFFER_ARB, sizeof(GLuint) * H2DV_GL_MAX_EDGE_BUFFER * 2, NULL, GL_DYNAMIC_DRAW_ARB);
747  GLenum err = glGetError();
748  if(err != GL_NO_ERROR) { //if it fails, no problem
749  glDeleteBuffersARB(1, &gl_edge_inx_buffer);
750  gl_edge_inx_buffer = 0;
751  }
752  }
753  glBindBufferARB(GL_ARRAY_BUFFER_ARB, 0);
754  }
755  catch(std::exception &e)
756  { //out-of-memory or any other failure
757  if(gl_coord_buffer) { glDeleteBuffersARB(1, &gl_coord_buffer); gl_coord_buffer = 0; }
758  if(gl_index_buffer) { glDeleteBuffersARB(1, &gl_index_buffer); gl_index_buffer = 0; }
759  if(gl_edge_inx_buffer) { glDeleteBuffersARB(1, &gl_edge_inx_buffer); gl_edge_inx_buffer = 0; }
760  }
761  }
762  }
763 
764  void ScalarView::draw_values_2d()
765  {
766  //set texture for coloring
767  glEnable(GL_TEXTURE_1D);
768  glBindTexture(GL_TEXTURE_1D, gl_pallete_tex_id);
769  glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_DECAL);
770  glColor3f(0, 1, 0);
771 
772  //set texture transformation matrix
773  glMatrixMode(GL_TEXTURE);
774  glLoadIdentity();
775  glTranslated(tex_shift, 0.0, 0.0);
776  glScaled(tex_scale, 0.0, 0.0);
777 
778  //render triangles
779  if(gl_coord_buffer == 0 || gl_index_buffer == 0)
780  { //render using the safe but slow methoold and slow method
781  //obtain data
782  int tri_cnt = lin->get_num_triangles();
783  const double3* vert = lin->get_vertices();
784  const int3* tris = lin->get_triangles();
785 
786  //render
787  glBegin(GL_TRIANGLES);
788  for (int i = 0; i < tri_cnt; i++)
789  {
790  const int3& triangle = tris[i];
791  const double3& vert_a = vert[triangle[0]];
792  const double3& vert_b = vert[triangle[1]];
793  const double3& vert_c = vert[triangle[2]];
794 
795  if(finite(vert_a[2]) && finite(vert_b[2]) && finite(vert_c[2]))
796  {
797  glTexCoord1d((vert_a[2] - range_min) * value_irange);
798  glVertex2d(vert_a[0], vert_a[1]);
799  glTexCoord1d((vert_b[2] - range_min) * value_irange);
800  glVertex2d(vert_b[0], vert_b[1]);
801  glTexCoord1d((vert_c[2] - range_min) * value_irange);
802  glVertex2d(vert_c[0], vert_c[1]);
803  }
804  }
805  glEnd();
806  }
807  else { //render using vertex buffer object
808  if(gl_tri_cnt > 0)
809  {
810  //bind vertices
811  glBindBufferARB(GL_ARRAY_BUFFER_ARB, gl_coord_buffer);
812  glVertexPointer(2, GL_FLOAT, sizeof(GLVertex2), GL_BUFFER_OFFSET(0));
813  glTexCoordPointer(1, GL_FLOAT, sizeof(GLVertex2), GL_BUFFER_OFFSET(GLVertex2::H2D_OFFSETOF_COORD));
814  glEnableClientState(GL_VERTEX_ARRAY);
815  glEnableClientState(GL_TEXTURE_COORD_ARRAY);
816 
817  //bind indices
818  glBindBufferARB(GL_ELEMENT_ARRAY_BUFFER_ARB, gl_index_buffer);
819 
820  //render
821  glDrawElements(GL_TRIANGLES, 3*gl_tri_cnt, GL_UNSIGNED_INT, GL_BUFFER_OFFSET(0));
822 
823  //GL cleanup
824  glDisableClientState(GL_VERTEX_ARRAY);
825  glDisableClientState(GL_TEXTURE_COORD_ARRAY);
826  }
827  }
828 
829  //GL clenaup
830  glMatrixMode(GL_TEXTURE); //switch-off texture transform
831  glLoadIdentity();
832  glMatrixMode(GL_MODELVIEW);
833  }
834 
835  void ScalarView::draw_edges_2d()
836  {
837  glColor3fv(edges_color);
838  bool displayed = false;
839  if(gl_edge_inx_buffer != 0) {//VBO
840  //prepage GL
841  glEnableClientState(GL_VERTEX_ARRAY);
842 
843  //map edge buffer
844  glBindBufferARB(GL_ARRAY_BUFFER_ARB, gl_edge_inx_buffer);
845  GLuint *gl_inx_buffer = (GLuint*)glMapBufferARB(GL_ARRAY_BUFFER_ARB, GL_WRITE_ONLY);
846  GLenum err = glGetError();
847  if(err == GL_NO_ERROR) {//render edges
848  unsigned int buffer_inx = 0;
849  int2* edges = lin->get_edges();
850  int* edge_markers = lin->get_edge_markers();
851  int edge_cnt = lin->get_num_edges();
852  for (int i = 0; i < edge_cnt; i++) //TODO: we could draw only left-right, top-bottom ones
853  {
854  const int2 &edge = edges[i];
855  const int &edge_marker = edge_markers[i];
856  if(show_edges || edge_marker != 0) { //select internal edges to draw
857  gl_inx_buffer[buffer_inx] = (GLuint)edge[0];
858  gl_inx_buffer[buffer_inx + 1] = (GLuint)edge[1];
859  buffer_inx += 2;
860  }
861 
862  //render buffer if it is full or if this is the last edge processed
863  if(buffer_inx == (2*H2DV_GL_MAX_EDGE_BUFFER) || (buffer_inx > 0 && i == (edge_cnt-1)))
864  {
865  glUnmapBufferARB(GL_ARRAY_BUFFER_ARB);
866 
867  //bind vertices and buffers
868  glBindBufferARB(GL_ARRAY_BUFFER_ARB, gl_coord_buffer);
869  glVertexPointer(2, GL_FLOAT, sizeof(GLVertex2), GL_BUFFER_OFFSET(0));
870  glBindBufferARB(GL_ELEMENT_ARRAY_BUFFER_ARB, gl_edge_inx_buffer);
871 
872  //render
873  glDrawElements(GL_LINES, buffer_inx, GL_UNSIGNED_INT, GL_BUFFER_OFFSET(0));
874 
875  //map edge buffer
876  glBindBufferARB(GL_ELEMENT_ARRAY_BUFFER_ARB, 0);
877  glBindBufferARB(GL_ARRAY_BUFFER_ARB, gl_edge_inx_buffer);
878  gl_inx_buffer = (GLuint*)glMapBufferARB(GL_ARRAY_BUFFER_ARB, GL_WRITE_ONLY);
879  buffer_inx = 0;
880  }
881  }
882  glUnmapBufferARB(GL_ARRAY_BUFFER_ARB);
883 
884  //GL cleanup
885  glDisableClientState(GL_VERTEX_ARRAY);
886 
887  displayed = true;
888  }
889  }
890 
891  //safe fallback
892  if(!displayed) { //VBO not suppored
893  glBegin(GL_LINES);
894  draw_edges(&draw_gl_edge, NULL, !show_edges);
895  glEnd();
896  }
897  }
898 
899  void ScalarView::draw_normals_3d()
900  {
901  double normal_xzscale = 1.0 / xzscale, normal_yscale = 1.0 / yscale;
902 
903  glPushAttrib(GL_ENABLE_BIT);
904  glPushMatrix();
905 
906  const int num_vert = lin->get_num_vertices();
907  double3* vert = lin->get_vertices();
908 
909  glDisable(GL_LIGHTING);
910  glDisable(GL_TEXTURE_1D);
911  glScaled(1, 1, -1);
912  glColor3f(0.8f, 0.5f, 0.5f);
913  glBegin(GL_LINES);
914  for(int i = 0; i < num_vert; i++)
915  {
916  double x = (vert[i][0] - xctr) * xzscale;
917  double y = (vert[i][2] - yctr) * yscale;
918  double z = (vert[i][1] - zctr) * xzscale;
919  glVertex3d(x, y, z);
920  glVertex3d(x + 0.1*normals[i][0]*normal_xzscale, y + 0.1*normals[i][2]*normal_yscale, z + 0.1*normals[i][1]*normal_xzscale);
921  }
922  glEnd();
923 
924  glPopMatrix();
925  glPopAttrib();
926  }
927 
928  void ScalarView::draw_gl_edge(int inx_vert_a, int inx_vert_b, ScalarView* viewer, void* param)
929  {
930  double3* verts = viewer->lin->get_vertices();
931  glVertex2d(verts[inx_vert_a][0], verts[inx_vert_a][1]);
932  glVertex2d(verts[inx_vert_b][0], verts[inx_vert_b][1]);
933  }
934 
935  void ScalarView::draw_edges(DrawSingleEdgeCallback draw_single_edge, void* param, bool boundary_only)
936  {
937  int2* edges = lin->get_edges();
938  int* edge_markers = lin->get_edge_markers();
939  int edge_cnt = lin->get_num_edges();
940  for (int i = 0; i < edge_cnt; i++) //TODO: we could draw only left-right, top-bottom ones
941  {
942  const int2 &edge = edges[i];
943  const int &edge_marker = edge_markers[i];
944  if(!boundary_only || edge_marker != 0) //select internal edges to draw
945  draw_single_edge(edge[0], edge[1], this, param);
946  }
947  }
948 
949 #define V0 vertices_min_x - xctr, range_min - yctr, -(vertices_min_y - zctr)
950 #define V1 vertices_max_x - xctr, range_min - yctr, -(vertices_min_y - zctr)
951 #define V2 vertices_max_x - xctr, range_min - yctr, -(vertices_max_y - zctr)
952 #define V3 vertices_min_x - xctr, range_min - yctr, -(vertices_max_y - zctr)
953 #define V4 vertices_min_x - xctr, range_max - yctr, -(vertices_min_y - zctr)
954 #define V5 vertices_max_x - xctr, range_max - yctr, -(vertices_min_y - zctr)
955 #define V6 vertices_max_x - xctr, range_max - yctr, -(vertices_max_y - zctr)
956 #define V7 vertices_min_x - xctr, range_max - yctr, -(vertices_max_y - zctr)
957 
958  void ScalarView::draw_aabb()
959  {
960  // Axis-aligned bounding box of the model.
961  GLdouble aabb[] =
962  {
963  V0, V1, V2, V3, // bottom
964  V0, V1, V5, V4, // front
965  V0, V3, V7, V4, // left
966  V1, V2, V6, V5, // right
967  V2, V3, V7, V6, // back
968  V4, V5, V6, V7 // top
969  };
970 
971  // Set the edge color.
972  glColor3fv(edges_color);
973 
974  // Make the cube wire frame.
975  glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
976 
977  // Scale the box appropriately.
978  glPushMatrix();
979  glScaled(xzscale, yscale, xzscale);
980 
981  // Activate and specify pointer to vertex array.
982  glEnableClientState(GL_VERTEX_ARRAY);
983  glVertexPointer(3, GL_DOUBLE, 0, aabb);
984 
985  // Draw the box (glDrawElements would be better, but do not work for me).
986  glDrawArrays(GL_QUADS, 0, 24);
987 
988  // Deactivate vertex arrays after drawing.
989  glDisableClientState(GL_VERTEX_ARRAY);
990 
991  // Recover the original model/view matrix
992  glPopMatrix();
993  }
994 
995  void ScalarView::on_display()
996  {
997  int i, j;
998 
999  // lock and get data
1000  lin->lock_data();
1001  int3* tris = lin->get_triangles();
1002  double3* vert = lin->get_vertices();
1003  int2* edges = lin->get_edges();
1004 
1005  glPolygonMode(GL_FRONT_AND_BACK, pmode ? GL_LINE : GL_FILL);
1006 
1007  //prepare vertices
1008  prepare_gl_geometry();
1009 
1010  if(!mode3d)
1011  {
1012  set_ortho_projection();
1013  glDisable(GL_LIGHTING);
1014  glDisable(GL_DEPTH_TEST);
1015  glDisable(GL_TEXTURE_2D);
1016  glDisable(GL_BLEND);
1017 
1018  // setup transformation (follows View::transform_x and View::transform_y)
1019  glMatrixMode(GL_MODELVIEW);
1020  glPushMatrix();
1021  glLoadIdentity();
1022  glTranslated(center_x, center_y, 0.0);
1023  glScaled(1.0, -1.0, 1.0);
1024  glTranslated(trans_x, trans_y, 0.0);
1025  glScaled(scale, scale, 1.0);
1026 
1027  // draw all triangles
1028  if(show_values)
1029  draw_values_2d();
1030 
1031  // draw contours
1032  glDisable(GL_TEXTURE_1D);
1033  if(contours)
1034  {
1035  tris = lin->get_contour_triangles();
1036  glColor3fv(cont_color);
1037  glBegin(GL_LINES);
1038  for (i = 0; i < lin->get_num_contour_triangles(); i++)
1039  {
1040  if(finite(vert[tris[i][0]][2]) && finite(vert[tris[i][1]][2]) && finite(vert[tris[i][2]][2]))
1041  {
1042  draw_tri_contours(vert, &tris[i]);
1043  }
1044  }
1045  tris = lin->get_triangles();
1046  glEnd();
1047  }
1048 
1049  // draw edges and boundary of mesh
1050  draw_edges_2d();
1051 
1052  //draw element IDS
1053  if(show_element_info)
1054  draw_element_infos_2d();
1055 
1056  //draw nodes and node info
1057  if(show_edges && allow_node_selection)
1058  draw_vertex_nodes();
1059 
1060  //cleanup
1061  glPopMatrix();
1062  }
1063  else
1064  {
1065  set_3d_projection(fovy, znear, zfar);
1066 
1067  glClear(GL_DEPTH_BUFFER_BIT);
1068  glEnable(GL_DEPTH_TEST);
1069 
1070  // Set camera transforamtion.
1071  glMatrixMode(GL_MODELVIEW);
1072  glLoadIdentity();
1073 
1074  // Initialize light and material.
1075  init_lighting();
1076 
1077  if(do_zoom_to_fit)
1078  {
1079  // If the user presses 'c' to center the view automatically, calculate how far to move
1080  // the visualized model away from the viewer so that it is completely visible in current window.
1081  ztrans = calculate_ztrans_to_fit_view();
1082  do_zoom_to_fit = false;
1083  }
1084 
1085  // Set model transformation.
1086  glTranslated(xtrans, ytrans, ztrans);
1087  glRotated(xrot, 1, 0, 0);
1088  glRotated(yrot, 0, 1, 0);
1089 
1090  // Draw the surface.
1091  glEnable(GL_LIGHTING);
1092  glEnable(GL_TEXTURE_1D);
1093  glBindTexture(GL_TEXTURE_1D, gl_pallete_tex_id);
1094  glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE);
1095  glEnable(GL_NORMALIZE);
1096  glEnable(GL_POLYGON_OFFSET_FILL);
1097  glPolygonOffset(1.0, 1.0);
1098  glBegin(GL_TRIANGLES);
1099  double normal_xzscale = 1.0 / xzscale, normal_yscale = 1.0 / yscale;
1100  for (i = 0; i < lin->get_num_triangles(); i++)
1101  {
1102  for (j = 0; j < 3; j++)
1103  {
1104  glNormal3d(normals[tris[i][j]][0] * normal_xzscale, normals[tris[i][j]][2] * normal_yscale, -normals[tris[i][j]][1] * normal_xzscale);
1105  glTexCoord2d((vert[tris[i][j]][2] - range_min) * value_irange * tex_scale + tex_shift, 0.0);
1106  glVertex3d((vert[tris[i][j]][0] - xctr) * xzscale,
1107  (vert[tris[i][j]][2] - yctr) * yscale,
1108  -(vert[tris[i][j]][1] - zctr) * xzscale);
1109  }
1110  }
1111  glEnd();
1112  glDisable(GL_POLYGON_OFFSET_FILL);
1113 
1114  // Draw edges.
1115  glDisable(GL_LIGHTING);
1116  glDisable(GL_TEXTURE_1D);
1117  if(show_edges)
1118  {
1119  glColor3fv(edges_color);
1120  glBegin(GL_LINES);
1121  for (i = 0; i < lin->get_num_edges(); i++)
1122  {
1123  glVertex3d((vert[edges[i][0]][0] - xctr) * xzscale,
1124  (vert[edges[i][0]][2] - yctr) * yscale,
1125  -(vert[edges[i][0]][1] - zctr) * xzscale);
1126  glVertex3d((vert[edges[i][1]][0] - xctr) * xzscale,
1127  (vert[edges[i][1]][2] - yctr) * yscale,
1128  -(vert[edges[i][1]][1] - zctr) * xzscale);
1129  }
1130  glEnd();
1131  }
1132 
1133  // Draw the whole bounding box or only the boundary edges.
1134  if(show_aabb)
1135  draw_aabb();
1136  else
1137  {
1138  glColor3fv(edges_color);
1139  glBegin(GL_LINES);
1140  for (i = 0; i < lin->get_num_edges(); i++)
1141  {
1142  // Outline of the domain boundary at the bottom of the plot or at the bottom user-defined limit
1143  double y_coord = (range_min - yctr) * yscale;
1144  int2& edge = edges[i];
1145  int& edge_marker = lin->get_edge_markers()[i];
1146  if(edge_marker)
1147  {
1148  glVertex3d((vert[edge[0]][0] - xctr) * xzscale, y_coord,
1149  -(vert[edge[0]][1] - zctr) * xzscale);
1150  glVertex3d((vert[edge[1]][0] - xctr) * xzscale, y_coord,
1151  -(vert[edge[1]][1] - zctr) * xzscale);
1152  }
1153  }
1154  glEnd();
1155  }
1156  }
1157 
1158  lin->unlock_data();
1159  }
1160 
1161  static inline void normalize(double& x, double& y, double& z)
1162  {
1163  double l = 1.0 / sqrt(sqr(x) + sqr(y) + sqr(z));
1164  x *= l; y *= l; z *= l;
1165  }
1166 
1167  void ScalarView::calculate_normals(double3* vert, int num_verts, int3* tris, int num_tris)
1168  {
1169  if(normals != NULL)
1170  delete [] normals;
1171  normals = new double3[num_verts];
1172  memset(normals, 0, sizeof(double3) * num_verts);
1173  for (int i = 0; i < num_tris; i++)
1174  {
1175  int3 &tri = tris[i];
1176  double ax = (vert[tri[1]][0] - vert[tri[0]][0]);
1177  double ay = (vert[tri[1]][1] - vert[tri[0]][1]);
1178  double az = (vert[tri[1]][2] - vert[tri[0]][2]);
1179 
1180  double bx = (vert[tri[2]][0] - vert[tri[0]][0]);
1181  double by = (vert[tri[2]][1] - vert[tri[0]][1]);
1182  double bz = (vert[tri[2]][2] - vert[tri[0]][2]);
1183 
1184  double nx = ay * bz - az * by;
1185  double ny = az * bx - ax * bz;
1186  double nz = ax * by - ay * bx;
1187  normalize(nx, ny, nz);
1188 
1189  for (int j = 0; j < 3; j++)
1190  {
1191  normals[tri[j]][0] += nx;
1192  normals[tri[j]][1] += ny;
1193  normals[tri[j]][2] += nz;
1194  }
1195  }
1196 
1197  for (int i = 0; i < num_verts; i++)
1198  normalize(normals[i][0], normals[i][1], normals[i][2]);
1199  }
1200 
1201  void ScalarView::update_layout()
1202  {
1203  View::update_layout();
1204  // (x, y, -z) coordinates (in the eye coordinate system) of the point that lies at the center of solution domain
1205  // and vertically at the position of the average value of all vertices. This point will be the origin for all
1206  // drawing and also the initial position of the camera.
1207  xctr = (vertices_max_x + vertices_min_x) / 2.0;
1208  yctr = value_range_avg;
1209  zctr = (vertices_max_y + vertices_min_y) / 2.0;
1210  }
1211 
1212  void ScalarView::reset_view(bool force_reset)
1213  {
1214  if(force_reset || view_not_reset) { // Reset 3d view.
1215  xrot = 40.0; yrot = 0.0;
1216  xtrans = ytrans = ztrans = 0.0;
1217  // Ensure that the model (before applying any transformations) may be translated along the view axis completely
1218  // into the viewing volume, and that it is not too flat (so that considerable amount of detail is visible).
1219  // Later on, we will determine the translation distance so that the model occupies as much of the view space as possible.
1220 
1221  // Set scaling into the (-1, 1) range on the x- and z- axes
1222  double max_radial = std::max(vertices_max_x - vertices_min_x, vertices_max_y - vertices_min_y);
1223  xzscale = 2.0 / max_radial;
1224 
1225  // Determine the largest y-axis distance of the function surface from the line of sight
1226  double max_axial = std::max(range_max - yctr, fabs(range_min - yctr));
1227 
1228  // We use the same scaling for the y-axis as for the x- and z- axes if after this scaling,
1229  // 1. the model could be translated so that it lies completely below the top clipping plane of the viewing frustum and
1230  // its farthest (away from the camera) corner is at least 1 unit before the far clipping plane (to allow some zooming out),
1231  // 2. the model's bounding box's part above (or below, whichever is bigger) the average function value (yctr) has dimensions
1232  // (-1, 1)x(0, height)x(-1, 1), where height > 0.1.
1233  // If this is not true, the model is either too tall or too flat and will be subject to different scaling
1234  // along the y-axis, such that it would fit to the viewing frustum at one third of its depth (i.e. at one third of
1235  // the total available zooming range).
1237  double tan_fovy_half = tan((double) fovy / 2.0 / 180.0 * M_PI);
1238  double max_allowed_height = (zfar-3)*tan_fovy_half;
1239  if(max_axial * xzscale > max_allowed_height || (max_axial * xzscale < 0.1 && !is_constant) )
1240  yscale = (znear + zfar) / 3.0 * tan_fovy_half * value_irange;
1241  else
1242  yscale = xzscale;
1243 
1244  do_zoom_to_fit = true;
1245  }
1246  View::reset_view(force_reset); // Reset 2d view.
1247  }
1248 
1249  double ScalarView::calculate_ztrans_to_fit_view()
1250  {
1251  // Axis-aligned bounding box of the model (homogeneous coordinates in the model space), divided into the bottom and top base.
1252  GLdouble aabb[2][16] =
1253  {
1254  {
1255  V0, 1,
1256  V1, 1,
1257  V2, 1,
1258  V3, 1
1259  },
1260  {
1261  V4, 1,
1262  V5, 1,
1263  V6, 1,
1264  V7, 1
1265  }
1266  };
1267  GLdouble aabb_base[16];
1268 
1269  // Tangents of the half of the vertical and horizontal field of view angles.
1270  double tan_fovy_half = tan((double) fovy / 2.0 / 180.0 * M_PI);
1271  double tan_fovx_half = tan_fovy_half * (double) output_width / output_height;
1272 
1273  // The viewpoint z-coordinate we are looking for.
1274  double optimal_viewpoint_pos = 0.0;
1275 
1276  // We will change the model transformation matrix, so save the current one.
1277  glPushMatrix();
1278 
1279  // Create the model transformation matrix with the default z-translation.
1280  glLoadIdentity();
1281  glTranslated(xtrans, ytrans, ztrans);
1282  glRotated(xrot, 1, 0, 0);
1283  glRotated(yrot, 0, 1, 0);
1284  glScaled(xzscale, yscale, xzscale);
1285 
1286  // As far as I know, OpenGL can only perform 4x4 matrix multiplication, so we find the 8 transformed bounding box corners by
1287  // first multiplying the transformation matrix with a matrix having as its 4 columns the coordinates of the bottom base corners
1288  // and then with a matrix created from the top base corners.
1289  for (int i = 0; i < 2; i++)
1290  {
1291  glPushMatrix();
1292  glMultMatrixd(aabb[i]);
1293  glGetDoublev( GL_MODELVIEW_MATRIX, aabb_base );
1294  glPopMatrix();
1295 
1296  // Go through the transformed corners and find the biggest distance of its perspective projection center to the origin.
1297  GLdouble *coord_ptr = &aabb_base[0];
1298  for (int j = 0; j < 4; j++)
1299  {
1300  double perspective_center_to_origin_dist = fabs(coord_ptr[0]) / tan_fovx_half + coord_ptr[2];
1301  if(perspective_center_to_origin_dist > optimal_viewpoint_pos)
1302  optimal_viewpoint_pos = perspective_center_to_origin_dist;
1303 
1304  perspective_center_to_origin_dist = fabs(coord_ptr[1]) / tan_fovy_half + coord_ptr[2];
1305  if(perspective_center_to_origin_dist > optimal_viewpoint_pos)
1306  optimal_viewpoint_pos = perspective_center_to_origin_dist;
1307 
1308  coord_ptr += 4;
1309  }
1310  }
1311 
1312  // Restore the original model transformation matrix.
1313  glPopMatrix();
1314 
1315  return -optimal_viewpoint_pos;
1316  }
1317 
1318  void ScalarView::set_vertical_scaling(double sc)
1319  {
1320  if(mode3d)
1321  yscale *= sc;
1322  else if(contours)
1323  cont_step *= sc;
1324  refresh();
1325  }
1326 
1327  void ScalarView::set_min_max_range(double min, double max)
1328  {
1330  if(fabs(max-min) < 1e-8)
1331  {
1332  this->warn("Range (%f, %f) is too narrow: adjusted to (%f, %f)", min, max, min-0.5, max);
1333  min -= 0.5;
1334  }
1335  View::set_min_max_range(min, max);
1336  }
1337 
1338  void ScalarView::init_lighting()
1339  {
1340  float light_specular[] = { 1.0f, 1.0f, 1.0f, 1.0f };
1341  float light_ambient[] = { 0.1f, 0.1f, 0.1f, 1.0f };
1342  float light_diffuse[] = { 1.0f, 1.0f, 1.0f, 1.0f };
1343 
1344  glEnable(GL_LIGHT0);
1345  glLightfv(GL_LIGHT0, GL_SPECULAR, light_specular);
1346  glLightfv(GL_LIGHT0, GL_AMBIENT, light_ambient);
1347  glLightfv(GL_LIGHT0, GL_DIFFUSE, light_diffuse);
1348  // glLightfv(GL_LIGHT0, GL_POSITION, light_position);
1349 
1350  float material_ambient[] = { 0.5f, 0.5f, 0.5f, 1.0f };
1351  float material_diffuse[] = { 0.8f, 0.8f, 0.8f, 1.0f };
1352  float material_specular[] = { 1.0f, 1.0f, 1.0f, 1.0f };
1353 
1354  glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT, material_ambient);
1355  glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, material_diffuse);
1356  glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, material_specular);
1357  glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 128);
1358  glDisable(GL_COLOR_MATERIAL);
1359 
1360  glShadeModel(GL_SMOOTH);
1361  glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, 1);
1362  if(GLEW_EXT_separate_specular_color)
1363  {
1364  glLightModeli(GL_LIGHT_MODEL_COLOR_CONTROL_EXT, GL_SEPARATE_SPECULAR_COLOR_EXT);
1365  }
1366  }
1367 
1368  void ScalarView::set_3d_mode(bool enable)
1369  {
1370  mode3d = enable; dragging = scaling = false;
1371  if(mode3d)
1372  {
1373  lin->lock_data();
1374  if(normals == NULL)
1375  calculate_normals(lin->get_vertices(), lin->get_num_vertices(), lin->get_triangles(), lin->get_num_triangles());
1376  lin->unlock_data();
1377  }
1378  refresh();
1379  }
1380 
1381  void ScalarView::on_key_down(unsigned char key, int x, int y)
1382  {
1383  switch (key)
1384  {
1385  case 'm':
1386  show_edges = !show_edges;
1387  refresh();
1388  break;
1389 
1390  case 'l':
1391  pmode = !pmode;
1392  refresh();
1393  break;
1394 
1395  case 'c':
1396  {
1397  reset_view(true);
1398  refresh();
1399  break;
1400  }
1401 
1402  case 'f':
1403  set_palette_filter(pal_filter != GL_LINEAR);
1404  break;
1405 
1406  case 'k':
1407  contours = !contours;
1408  refresh();
1409  break;
1410 
1411  case 'i':
1412  show_element_info = !show_element_info;
1413  if(!mode3d)
1414  refresh();
1415  break;
1416 
1417  case 'b':
1418  show_aabb = !show_aabb;
1419  if(mode3d)
1420  refresh();
1421  break;
1422 
1423  case '3':
1424  {
1425  mode3d = !mode3d;
1426  dragging = scaling = false;
1427  if(mode3d)
1428  {
1429  lin->lock_data();
1430  if(normals == NULL)
1431  calculate_normals(lin->get_vertices(), lin->get_num_vertices(), lin->get_triangles(), lin->get_num_triangles());
1432  lin->unlock_data();
1433  }
1434  refresh();
1435  break;
1436  }
1437 
1438  case '*':
1439  lin->lock_data();
1440  if(mode3d)
1441  yscale *= D3DV_SCALE_STEP_COEF;
1442  else if(contours)
1443  cont_step *= D3DV_SCALE_STEP_COEF;
1444  lin->unlock_data();
1445  refresh();
1446  break;
1447 
1448  case '/':
1449  lin->lock_data();
1450  if(mode3d)
1451  yscale /= D3DV_SCALE_STEP_COEF;
1452  else if(contours)
1453  cont_step /= D3DV_SCALE_STEP_COEF;
1454  lin->unlock_data();
1455  refresh();
1456  break;
1457 
1458  default:
1459  View::on_key_down(key, x, y);
1460  break;
1461  }
1462  }
1463 
1464 
1465  void ScalarView::on_mouse_move(int x, int y)
1466  {
1467  if(mode3d && (dragging || scaling || panning))
1468  {
1469  if(dragging)
1470  {
1471  yrot += 0.4 * (x - mouse_x);
1472  xrot += 0.4 * (y - mouse_y);
1473 
1474  if(xrot < -90) xrot = -90;
1475  else if(xrot > 90) xrot = 90;
1476  }
1477  else if(scaling)
1478  {
1479  ztrans += 0.01 * (mouse_y - y);
1480  if(ztrans > -0.25) ztrans = -0.25;
1481  else if(ztrans < -7) ztrans = -7;
1482  }
1483  else
1484  {
1485  xtrans += 0.002 * (x - mouse_x);
1486  ytrans += 0.002 * (mouse_y - y);
1487  }
1488 
1489  refresh();
1490  mouse_x = x;
1491  mouse_y = y;
1492  return;
1493  }
1494  else
1495  {
1496  if(!mode3d && show_edges && !dragging && !scaling && !panning)
1497  {
1498  if(allow_node_selection)
1499  {
1500  VertexNodeInfo* found_node = find_nearest_node_in_range((float)untransform_x(x), (float)untransform_y(y), (float)(node_pixel_radius / scale));
1501  if(found_node != pointed_vertex_node)
1502  pointed_vertex_node = found_node;
1503  }
1504  else
1505  pointed_vertex_node = NULL;
1506  refresh();
1507  }
1508  else
1509  {
1510  View::on_mouse_move(x, y);
1511  }
1512  }
1513  }
1514 
1515  void ScalarView::on_right_mouse_down(int x, int y)
1516  {
1517  //handle node selection
1518  if(allow_node_selection && pointed_vertex_node != NULL)
1519  {
1520  if(pointed_vertex_node->selected) //deselct node
1521  {
1522  pointed_vertex_node->selected = false;
1523  }
1524  else { //select node
1525  pointed_vertex_node->selected = true;
1526  }
1527  refresh();
1528  }
1529  else { //avoid mixing of functionality
1530  View::on_right_mouse_down(x, y);
1531  }
1532  }
1533 
1534  void ScalarView::on_middle_mouse_down(int x, int y)
1535  {
1536  if(!mode3d) return;
1537  dragging = scaling = false;
1538  panning = true;
1539  }
1540 
1541  void ScalarView::on_middle_mouse_up(int x, int y)
1542  {
1543  panning = false;
1544  }
1545 
1546 
1547  const char* ScalarView::get_help_text() const
1548  {
1549  return
1550  "ScalarView\n\n"
1551  "Controls:\n"
1552  " Left mouse - pan\n"
1553  " Right mouse - zoom\n"
1554  " 3 - toggle 3D mode\n"
1555  " C - center image\n"
1556  " F - toggle smooth palette\n"
1557  " H - render high-quality frame\n"
1558  " K - toggle contours\n"
1559  " M - toggle mesh\n"
1560  " B - toggle bounding box\n"
1561  " I - toggle element IDs (2d only)\n"
1562  " P - cycle palettes\n"
1563  " S - save screenshot\n"
1564  " * - increase contour density\n"
1565  " / - decrease contour density\n"
1566  " F1 - this help\n"
1567  " Esc, Q - quit\n\n"
1568  "3D mode:\n"
1569  " Left mouse - rotate\n"
1570  " Right mouse - zoom\n"
1571  " Middle mouse - pan\n"
1572  " * - increase Z scale\n"
1573  " / - decrease Z scale";
1574  }
1575  }
1576  }
1577 }
1578 #endif